英文原版地址:http://www.cs.rochester.edu/~bh/cs400/using_lapack.html
以下是翻译部分:
很多人都因为不知道怎么使用CLAPACK而不得不转向MATLAB。其实CLPACK只是对LAPACK库的一个比较好的(pretty)封装。在这里,我向大家展示(使用CLAPACK库)的技巧。当我说CLAPACK库的时候,大部分时候我指的是LAPACK库。前者(CLAPACK)是后者(LAPACK)的一个f2c版本的(f2c: fortran to C, LAPACK库是用fortran写的)。
CLAPACK函数的命名规则
CLAPACK的函数都采用XYYZZZ的命名格式,X表示函数期望接收的数据类型。下表是X可以采用的类型及其意义的映射关系。
S | 实数,对应C中的float |
D | 双精度,对应C中的double |
C | 复数 |
Z | 双精度复数 |
字母YY则告诉函数将要处理的矩阵类型。GE和TR,分别表示一般的矩阵和三角矩阵。ZZZ表示真正进行的运算,例如:SV 函数式解Ax = b类型系统的的简单驱动函数(这里驱动函数(drivers),意思是指底层进行运算的函数),LS 是指用线性平方驱动函数来解(aX = b类型的系统)。
解Ax = b类型的线性系统的函数
我们将会使用简单驱动函数或者线性平方驱动函数来解双精度类型的的一般矩阵,所以函数的名字是(或者以这个开头)DGESV或者DGELS。DGESV函数的C原型是:
void dgesv_(const int *N, const int *nrhs, double *A, const int *lda, int *ipiv, double *b, const int *ldb, int *info);
注意,所有的参数都是以引用传递的方式进行传递的。每个参数对应的意义如下:
N | 矩阵A的列数 |
nrhs | 矩阵b的列数,通常是1 |
lda | 矩阵A的行数,(Leading Dimension of A) |
IPIV | 主元的索引(pivot indices) |
ldb | 矩阵B的行数 |
更详细的信息,请参看使用手册,DGELS的原型是:
void dgels_(const char *trans, const int *M, const int *N, const int *nrhs, double *A, const int *lda, double *b, const int *ldb, double *work, const int *lwork, int *info);
你可以猜一下这些参数的意义,然后查看手册。
进行SVD(奇异值分解)的函数
你可以能在某些时候需要进行SVD(奇异值分解)。计算SVD的函数是XGESVD。dgesvd_函数的C原型是:
dgesvd_(const char* jobu, const char* jobvt, const int *M, const int* N, double *A, const int* lda, double *U, const int *ldu, double *VT, const int *ldvt, double *work, const int *lwork, const int *info);
示例代码
#include < stdio.h> void dgesv_(const int *N, const int *nrhs, double *A, const int *lda, int *ipiv, double *b, const int *ldb, int *info); void dgels_(const char *trans, const int *M, const int *N, const int *nrhs, double *A, const int *lda, double *b, const int *ldb, double *work, const int * lwork, int *info); int main(void) { /* 3x3 matrix A * 76 25 11 * 27 89 51 * 18 60 32 */ double A[9] = {76, 27, 18, 25, 89, 60, 11, 51, 32}; double b[3] = {10, 7, 43}; int N = 3; int nrhs = 1; int lda = 3; int ipiv[3]; int ldb = 3; int info; dgesv_(&N, &nrhs, A, &lda, ipiv, b, &ldb, &info); if(info == 0) /* succeed */ printf("The solution is %lf %lf %lf\n", b[0], b[1], b[2]); else fprintf(stderr, "dgesv_ fails %d\n", info); return info; } |
使用如下命令进行编译:
cc lapack tut.c -L/usr/grads/lib -llapack -lblas -lF77 -lI77
注意:矩阵A是列主序矩阵。(这是)Fortan另外一个让人感到意外的地方。