logo资料库

LAPACK函数介绍.pdf

第1页 / 共11页
第2页 / 共11页
第3页 / 共11页
第4页 / 共11页
第5页 / 共11页
第6页 / 共11页
第7页 / 共11页
第8页 / 共11页
资料共11页,剩余部分请下载后查看
为方便线性代数运算,现将LAPACK 中的函数介绍如下: 1.函数的命名规则: LAPACK 里的每个函数名已经说明了该函数的使用规则。所有函数都是以 XYYZZZ 的形 式命名,对于某些函数,没有第六个字符,只是 XYYZZ 的形式。 第一个字母 X 代表以下的数据类型: S D C Z 注: 在新版 LAPACK 中含有使用重复迭代法的函数 DSGESV 和 ZCDESV。 头 2 个字母表示使用的精度: DS REAL,单精度实数 DOUBLE PRECISION,双精度实数 COMPLEX,单精度复数 COMPLEX*16 或 DOUBLE COMPLEX 输入数据是 double 双精度,算法使用单精度 输入数据是 complex*16,算法使用 complex 单精度复数 ZC 接下面两个字母 YY 代表数组的类型。 bidiagonal,双对角矩阵 diagonal,对角矩阵 general band,一般带状矩阵 general (i.e., unsymmetric, in some cases rectangular),一般情形(即非对称, 在有些情形下为矩形) general matrices, generalized problem (i.e., a pair of general matrices),一般矩 阵,广义问题(即一对一般矩阵) general tridiagonal,一般三对角矩阵 (complex) Hermitian band,(复数)厄尔米特带状阵 (complex) Hermitian,(复数)厄尔米特矩阵 upper Hessenberg matrix, generalized problem (i.e a Hessenberg and a triangular matrix),上海森伯格矩阵,广义问题(即一个海森伯格矩阵和一 个三角矩阵) (complex) Hermitian, packed storage,(复数)压缩储存的厄尔米特矩阵 upper Hessenberg,上海森博格矩阵 (real) orthogonal, packed storage,(实数)压缩储存的正交阵 (real) orthogonal,(实数)正交阵 symmetric or Hermitian positive definite band,对称或厄尔米特正定带状矩 阵 symmetric or Hermitian positive definite,对称或厄尔米特正定矩阵 symmetric or Hermitian positive definite, packed storage,压缩储存的对称或 厄尔米特正定矩阵 symmetric or Hermitian positive definite tridiagonal,对称或厄尔米特正定三 对角阵 (real) symmetric band,(实数)对称带状阵 symmetric, packed storage,压缩储存的对称阵 (real) symmetric tridiagonal,(实数)对称三对角阵 BD DI GB GE GG GT HB HE HG HP HS OP OR PB PO PP PT SB SP ST
SY TB TG TP TR TZ UN UP symmetric,对称阵 triangular band,三角形带状矩阵 triangular matrices, generalized problem (i.e., a pair of triangular matrices),三 角形矩阵,广义问题(即一对三角形阵) triangular, packed storage,压缩储存的三角形阵 triangular (or in some cases quasi-triangular),三角形阵(在某些情形下为类 三角形阵) trapezoidal,梯形阵 (complex) unitary,(复数)酉矩阵 (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 ++ ){
分享到:
收藏