LAPACK函数介绍不完全版

为方便线性代数运算，现将LAPACK中的函数介绍如下：

1.函数的命名规则：

LAPACK里的每个函数名已经说明了该函数的使用规则。所有函数都是以XYYZZZ的形式命名，对于某些函数，没有第六个字符，只是XYYZZ的形式。

第一个字母X代表以下的数据类型：

S REAL，单精度实数

D DOUBL

E PRECISION，双精度实数

C COMPLEX，单精度复数

Z COMPLEX*16 或DOUBLE COMPLEX

注：

在新版LAPACK中含有使用重复迭代法的函数DSGESV和ZCDESV。

头2个字母表示使用的精度：

DS 输入数据是double双精度，算法使用单精度

ZC 输入数据是complex*16，算法使用complex单精度复数

接下面两个字母YY代表数组的类型。

BD bidiagonal，双对角矩阵

DI diagonal，对角矩阵

GB general band，一般带状矩阵

GE general (i.e., unsymmetric, in some cases rectangular)，一般情形（即非对称，在有些情形下为矩形）

GG general matrices, generalized problem (i.e., a pair of general matrices)，一般矩阵，广义问题（即一对一般矩阵）

GT general tridiagonal，一般三对角矩阵

HB (complex) Hermitian band，（复数）厄尔米特带状阵HE (complex) Hermitian，（复数）厄尔米特矩阵

HG upper Hessenberg matrix, generalized problem (i.e a Hessenberg and a triangular matrix)，上海森伯格矩阵，广义问题（即一个海森伯格矩阵和一个三角矩阵）

HP (complex) Hermitian, packed storage，（复数）压缩储存的厄尔米特矩阵HS upper Hessenberg，上海森博格矩阵

OP (real) orthogonal, packed storage，（实数）压缩储存的正交阵

OR (real) orthogonal，（实数）正交阵

PB symmetric or Hermitian positive definite band，对称或厄尔米特正定带状矩阵

PO symmetric or Hermitian positive definite，对称或厄尔米特正定矩阵

PP symmetric or Hermitian positive definite, packed storage，压缩储存的对称或厄尔米特正定矩阵

PT symmetric or Hermitian positive definite tridiagonal，对称或厄尔米特正定三对角阵

SB (real) symmetric band，（实数）对称带状阵

SP symmetric, packed storage，压缩储存的对称阵

ST (real) symmetric tridiagonal，（实数）对称三对角阵

SY symmetric，对称阵

TB triangular band，三角形带状矩阵

TG triangular matrices, generalized problem (i.e., a pair of triangular matrices)，三角形矩阵，广义问题（即一对三角形阵）

TP triangular, packed storage，压缩储存的三角形阵

TR triangular (or in some cases quasi-triangular)，三角形阵（在某些情形下为类三角形阵）

TZ trapezoidal，梯形阵

UN (complex) unitary，（复数）酉矩阵

UP (complex) unitary, packed storage，（复数）压缩储存的酉矩阵

最后三个字母ZZZ代表计算方法。比如，SGEBRD是一个单精度函数，用于把一个实数一般阵压缩为双对角阵（a bidiagonal reduction，即BRD）。

2.函数讲解

2.1 dgesv_（）函数用来求解对称矩阵问题，否则真的要出错！

int dgesv_(

integer *n,

integer *nrhs,

doublereal *a,

integer *lda,

integer *ipiv,

doublereal *b,

integer *ldb,

integer *info

)；

从名字的意义上可以看出是用来解决双精度一般型的线性方程（组）的问题。

DGESV是用来求解实数的线性方程组AX=B的。A是N×N型矩阵，X和B是N×NRHS 型矩阵。

参数介绍：

N (input) INTEGER

. 线性方程组的个数，例如A矩阵的行数。N >= 0.

NRHS (input) INTEGER

右边矩阵的尺寸, 例如：B矩阵的列数. NRHS >= 0.

A (input/output) 双精度数组, 尺寸为LDA×N。

输入时为矩阵A的系数。

输出时, L和U是来自A = P*L*U的因式分解;L对角线的元素不被存储。.

LDA (input) INTEGER

数组A的主尺寸LDA >= max(1,N).

IPIV (output) INTEGER数组, 尺寸是N维。

The pivot indices that define the permutation matrix P;

row i of the matrix was interchanged with row IPIV(i).

B (input/output) 双精度数组, 尺寸为N×NRHS 。

输入时, 输入的是N×NRHS的B数组。

输出时, 如果INFO = 0, 输出N×NRHS为数组，即为方程组的解。

LDB (input) INTEGER

B数组的主尺寸。LDB >= max(1,N).

INFO (output) INTEGER

= 0: 成功退出。

< 0: 如果INFO = -i, 第i各自变量是一个不可接受的数值。

> 0: 如果INFO = i, U(i,i) 为0。因式分解已经完成，但是因式U是一个单

数，所以是不能够给出答案的。

与之对应的单精度方式为：sgesv_（）方程，只是把双精度的化为单精度就可以,另外一个就是zgesv _（）方程，它解的是复数形式。

2.2 dgeev_（）

对于非对称矩阵特征向量密集复杂的问题应用此函数，此函数的工能是求特征值。

d表示double。

ge表示general，说明是普通矩阵，按照列主序存储。

ev表示eigenvector（猜测），表达的是函数的功能。

int dgeev_(

char *jobvl,

char *jobvr,

integer *n,

doublereal *a,

integer *lda,

doublereal *wr,

doublereal *wi,

doublereal *vl,

integer *ldvl,

doublereal *vr,

integer *ldvr,

doublereal *work,

integer *lwork,

integer *info

)；

ZGEEV是用来求解N×N的非对称矩阵A，特征向量，或左面的特征向量或右边的特征向量是不对称的。

矩阵A的有特征向量满足：

A*V（j）= λ(j)*V(j); λ(j)是A的特征向量。

矩阵A的左特征向量满足：

u(j)**H * A = λ(j) * u(j)**H ；u(j)**H 表示的意思是u(j)的共轭转置。

参数介绍:

JOBVL (input) CHARACTER*1

= 'N': A的左特征向量没有被计算;

= 'V':A的左特征向量被计算了。

JOBVR (input) CHARACTER*1

= 'N': A的右特征向量没有被计算;

= 'V':A的右特征向量被计算了。

N (input) INTEGER

A矩阵的维数. N >= 0.

A (input/output) DOUBLE PRECISION数组, 尺寸是LDA×N。

输入时，A是N×N

输出时，A被覆盖重写。

LDA (input) INTEGER

矩阵A的主尺寸。LDA >= max(1,N).

WR (output) DOUBLE PRECISION数组，是N维的。

可以这样分配空间：

doublereal* wr = (doublereal*)malloc( sizeof(doublereal) * n) ;

WI (output) DOUBLE PRECISION 数组，是N维的。

可以这样分配空间：

doublereal* wi = (doublereal*)malloc( sizeof(doublereal) * n) ;

WR和WI 包含各自特征向量的实部和虚部分

VL (output) DOUBLE PRECISION 数组，尺寸为LDVL×N

doublereal* vl = (doublereal*)malloc( sizeof(doublereal) * n * ldvl) ;

如果JOBVL = 'V', 左特征向量u(j) 被存储one

after another in the columns of VL, in the same order

as their eigenvalues.

如果JOBVL = 'N', VL是不被引用的。

如果第j特征向量是real,则u(j) = VL(:,j),

the j-th column of VL.

If the j-th and (j+1)-st eigenvalues form a complex

conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and

u(j+1) = VL(:,j) - i*VL(:,j+1).

LDVL (input) INTEGER

矩阵VL的主尺寸。LDVL >= 1;如果JOBVL = 'V', LDVL >= N.

VR (output) DOUBLE PRECISION 数组, 尺寸为LDVR×N

doublereal* vr = (doublereal*)malloc( sizeof(doublereal) * n * ldvr) ;

如果JOBVR = 'V', 有特征向量v(j) are stored one

after another in the columns of VR, in the same order

as their eigenvalues.

If JOBVR = 'N', VR is not referenced.

If the j-th eigenvalue is real, then v(j) = VR(:,j),

the j-th column of VR.

If the j-th and (j+1)-st eigenvalues form a complex

conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and

v(j+1) = VR(:,j) - i*VR(:,j+1).

LDVR (input) INTEGER

矩阵VR的主尺寸. LDVR >= 1;

如果JOBVR = 'V', LDVR >= N.

WORK (workspace/output) DOUBLE PRECISION array,尺寸(MAX(1,LWORK)) 输出时，如果INFO = 0, WORK(1)返回LWORK.

doublereal *work = (doublereal*)malloc( sizeof(doublereal) * lwork) ;

LWORK (input) INTEGER

WORK的尺寸. LWORK >= max(1,3*N), 如果

JOBVL = 'V' 或者JOBVR = 'V', LWORK >= 4*N. 为了性能，LWORK 通常会大一点。

如果LWORK = -1, then a workspace query is assumed; the routine

only calculates the optimal size of the WORK array, returns

this value as the first entry of the WORK array, and no error

message related to LWORK is issued by XERBLA.

INFO (output) INTEGER

= 0: 成功退出

< 0: 如果INFO = -i, 第i各自变量是一个不可接受的数值。

> 0: 如果INFO = i, QR算法没能够得出所有的特征向量，没有被计算出的特征向量元素i+1:N of WR 和WI包含有融合的特征值。.

一个例子：

采用Matlab来验算：

可见计算结果是真确的！

网上查找到的一个例子：

CLAPACK的dgeev_求特征值

从这一个函数也可以认识LAPACK的一些习惯问题，

进而可以使用其他的函数。

》》特征值和特征向量

我们平时使用的一般都是右特征向量，定义为：

A*alpha = lamda * alpha

另外还有一种定义的左特征向量：

beta*A = lamda * beta

而该函数被设计为可以选择求出：

只有特征值，左边特征向量，右边特征向量

》》关于参数的解释

从头文件里面直接摘录的原型如下：

int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *

a, integer *lda, doublereal *wr,doublereal *wi, doublereal *vl,

integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work,

integer *lwork, integer *info);

首先，最明显的是，所有参数都按照指针传入。

然后这套函数库有个共同的习惯，

即要求调用者来处理空间，

包括提供返回值的空间，

用来计算的临时空间。

函数命名：

d表示double。

ge表示general，说明是普通矩阵，按照列主序存储。

ev表示eigenvector吧（疑问语气），表达的是函数的功能。

用到的类型解释：

char*，是字符串类型，但LAPACK的函数只关注该字符串的第一个字符integer,就是C里面的int，一般用来指定维度

doublereal，就是C里面的double，参数类型和函数前缀统一（这里是d）各个参数：

char *jobvl, 第一个字符为'N'，表示不求左特征向量，为'V'表示要求

char *jobvr, 同上，是对右特征向量的选项

integer *n, 矩阵的列数

doublereal *

a, 存储A矩阵的空间（列主序！！）

integer *lda, A矩阵的行数

doublereal *wr, 返回的特征值，实部

doublereal *wi, 返回的特征值，虚部

doublereal *vl, 左特征向量的存储空间，大小为ldvl*n

integer *ldvl, 左特征向量存储空间的行数

doublereal *vr, 右特征向量的存储空间，大小为ldvr*n

integer *ldvr, 右特征向量存储空间的行数

doublereal *work, 工作空间，一般至少要4*n

integer *lwork, 工作空间的大小

integer *info，返回信息，0表示成功，-i表示第i个参数有问题，+i表示执行错误

具体的看一下代码就明白了，

实在觉得恼火的，照搬着用也行。。

示例代码：

#include < stdio.h>

#pragma comment(lib , "blas.lib")

#pragma comment(lib , "clapack.lib")

void *malloc(size_t n) ;

#include "f2c.h"

#include "clapack.h"

int main(void)

{

/* 3x3 matrix A

* 76 25 11

* 27 89 51

* 18 60 32

*/

doublereal A[9] = {76, 27, 18, 25, 89, 60, 11, 51, 32};

integer info ;

int i , j ;

char jobvl = 'V' ;

char jobvr = 'V' ;

integer n = 3 ;

doublereal *a = A ;

integer lda = 3 ;

doublereal* wr = (doublereal*)malloc( sizeof(doublereal) * n) ;

doublereal* wi = (doublereal*)malloc( sizeof(doublereal) * n) ;

integer ldvr = 3 ;

doublereal* vr = (doublereal*)malloc( sizeof(doublereal) * n * ldvr) ;

integer ldvl = 3 ;

doublereal* vl = (doublereal*)malloc( sizeof(doublereal) * n * ldvl) ;

integer lwork = n * 4 ;

doublereal *work = (doublereal*)malloc( sizeof(doublereal) * lwork) ;

dgeev_(&jobvl, &jobvr, &n, a, &lda, wr, wi, vl , &ldvl , vr, &ldvr, work, &lwork, &info) ;

printf("info:%d\n" , info) ;

printf("D = \n") ;

for ( i = 0 ; i < n ; i ++ ){

for ( j = 0 ; j < n ; j ++ ){

if ( i == j ) printf("%10.5f" , wr[i]) ;

else printf("%10.5f" , 0.0) ;

}

printf("\n") ;

}

printf("Vl = \n") ;

for ( i = 0 ; i < n ; i ++ ){

for ( j = 0 ; j < n ; j ++ ){

printf("%10.5f" , vl[n * j + i]) ;

}

printf("\n") ;

}

printf("Vr = \n") ;

for ( i = 0 ; i < n ; i ++ ){

for ( j = 0 ; j < n ; j ++ ){

printf("%10.5f" , vr[n * j + i]) ;

}

printf("\n") ;

}

return info;

}

2.3 int sgetrf_() 通过行交换的方式将M×N的矩阵A进行LU分解。int sgetrf_(

integer *m,

integer *n,

real *a,

integer *lda,

integer *ipiv,

integer *info

);

矩阵变换有如下形式：

A = P * L * U

P矩阵为置换矩阵，L时下三角形式，U是上三角形式。

参数介绍：

M (input) INTEGER

矩阵A的行数。M >= 0.

N (input) INTEGER

矩阵A的列数。N >= 0.

A (input/output) REAL array, A（LDA×N)

输入时，是M×N的将被分解的矩阵。

输出时，A = P*L*U中的L,U;该单位对角元素的L不存储

LDA (input) INTEGER

矩阵A的主尺寸。LDA >= max(1,M).

IPIV (output) INTEGER array, dimension (min(M,N))

指点指数，for 1 <= i <= min(M,N), 矩阵的第i行与IPIV(i)行被交换。

integer* IPIV= (integer*)malloc(sizeof(integer)*min（M,N）);这样不会被报错。

INFO (output) INTEGER

= 0: 成功退出。

上述只是个简单的入门教程。下面我们将详细的介绍LAPACK+Users’+Guide的部分内容（实用的内容）。

3.函数例程：

在这个部分里面介绍LAPACK中的函数例程。更深一步的内容在下一部分。

3.1线性方程（组）

两种类型的韩函数解决线性方程（组）。

·一个简单的函数例程（以-SV结尾），用来解AX=B，通过转化A和B并将结果放入B中。·一个扩展的函数例程（以-SVX结尾），用来将解决下面的方程组：

-----解决B

A H＝，除非A是对称的或厄密共轭

X

A T＝或B

X

下面给个表格：