Windows系统下在C语言程序中调用CLAPACK进行矩阵运算

1、CLAPACK安装

Clapack安装步骤参考：如何在VC中调用CLAPACK(适用于VS2010)

简要流程如下：

首先从主页上下载CLAPACK包。http://www.netlib.org/clapack/ 上有很多版本。选择http://www.netlib.org/clapack/CLAPACK-3.1.1-VisualStudio.zip

下载解压后，我们可以看到如下目录结构:
/SRC             CLAPACK的源代码
/BLAS            BLAS的源代码
/F2CLIBS       F2C的源代码

/LIB                 编译后的库文件.lib
/INCLUDE       头文件
/TESTING       一些使用范例程序
/INSTALL        里面有UNIX和其他平台下安装的文件和方法。lawn81.pdf文件是CLAPACK的使用说明，大家使用的时候可以看看。

这里我要提醒大家。虽然软件包里已经有编译好的.lib文件，就在/LIB中。但是我建议大家不要直接使用他们，还是自己编译一下再用才保险。很多人就是因为直接使用他们而出错的。这点非常重要！如果你就是直接调用他们失败的，那不妨先自己编译一下再试试。

另外，CLAPACK需要F2CLIBS库，并且使用了blas和tmglib这两个外部的库。这几个库都已经包含在了CLAPACK的安装包中。所以，大家无需另外下载。当然，使用前我们还是要重新编译一下的，原因上面已经说过了。

接下来，编译他们。
首先双击 clapack.vcproj打开工程项目文件。工程中已经包含了所有的子项目。
我们根据个人需要，编译成debug模式，或者release模式。为了能在自己的程序调用中方便的进行debug，我这里就以debug模式为例说明。我的系统是win32。如果你是64位系统也是支持的，操作方法类似。

首先编译F2CLIBS，用于将fortran转换为c语言。
选择libf2c子项目。直接生成之。编译过程中可能会有一些warning，不要理会他们。编译成功后，输出文件是libf2cd.lib。这里的d就是debug模式，如果是release模式就是libf2c.lib。输出文件默认路径是/LIB文件夹。注意，/LIB/Win32下已经有一些lib了。大家最好把他们都先删除了，以免新旧文件混淆。

接着编译tmglib。这是一个矩阵库。
这个库在TESTING/MATGEN里面。选择他生成就好了。
输出文件还是在/LIB里面。文件名是tmglibd.lib。

然后是编译blas，选择项目blas, 编译之。输出文件BLASd.lib。

最后是编译CLAPACK，生成clapackd.lib.

2、CLAPACK调用

项目→属性→

C/C++→常规→附加包含目录：添加/INCLUDE目录。

连接器→附加库目录：添加/LIB/Win32目录。

输入→附加依赖项：添加libf2cd.lib BLASd.lib clapackd.lib tmglibd.lib等

3、CLAPACK测试

矩阵SVD分解

#include "stdio.h"    
#include "f2c.h"  
#include "clapack.h"  

#define SIZE 4  

int main()
{
    char JOBU;
    char JOBVT;

    int i;

    integer M = SIZE;
    integer N = SIZE;

    integer LDA = M;
    integer LDU = M;
    integer LDVT = N;
    integer LWORK;
    integer INFO;

    integer mn = min(M, N);
    integer MN = max(M, N);

    double a[SIZE*SIZE] = { 16.0, 5.0, 9.0, 4.0, 2.0, 11.0, 7.0, 14.0, 3.0, 10.0, 6.0, 15.0, 13.0, 8.0, 12.0, 1.0 };
    double s[SIZE];
    double wk[201];
    double uu[SIZE*SIZE];
    double vt[SIZE*SIZE];

    JOBU = 'A';
    JOBVT = 'A';
    LWORK = 201;

    /* Subroutine int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
    doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
    ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
    integer *info)
    */

    dgesvd_(&JOBU, &JOBVT, &M, &N, a, &LDA, s, uu,
        &LDU, vt, &LDVT, wk, &LWORK, &INFO);

    if (INFO > 0) {
        printf("The algorithm computing SVD failed to converge\n");
        exit(1);
    }

    for (i = 0; i< SIZE; i++) {
        printf("\n s[ %d ] = %f", i, s[i]);
    }
    getchar();
    return 0;
}

矩阵相乘：

#include "stdio.h"    
#include "f2c.h"  
#include "clapack.h"  
#include "blaswrap.h" 

int main() 
{
    char transa = 'T', transb = 'T';
    integer M = 2, N = 2, K = 2, LDA = K, LDB = N, LDC = M;
    double alpha = 1.0, A[4] = { 1, 2, 3, 4}, B[4] = { 5, 6, 7, 8 }, beta = 0.0, C[4]; 
    //下面的函数的意思是C = 1.0 * T(A) * T(B) + 0 * C，其中T()表示将某个矩阵转置
    //注意此时得到的C是按列存放的
    dgemm_(&transa, &transb, &M, &N, &K, &alpha, A, &LDA, B, &LDB, &beta, C, &LDC);
    
    return 0;
}