Lapack 科学计算包

目录

• 环境包
• BLAS 库和 CBLAS 接口
  • 安装 blas
  • 安装 cblas
• 安装 LAPACK
• 基本框架
  • 概述
  • 函数命名规则
  • 编译指令
• LAPACK 语法
  • 区分例程
  • 参数格式
• LAPACK 函数
  • degsv
  • dgels
  • dsyev
  • dgees
  • dgecon
  • dgehrd
  • dgeqrf
  • dgetrf
  • dgetri
• LAPACK 实例

 

本文整理了科学计算包 Lapack 的安装过程和使用规范。

 

环境包

需要安装 gfortran 和 cmake

sudo apt-get install gfortran

 

BLAS 库和 CBLAS 接口

安装 blas

wget http://www.netlib.org/blas/blas.tgz
tar -xzvf blas.tgz
cd BLAS-3.10.0
make
sudo cp blas_LINUX.a /usr/local/lib/libblas.a
# 将 blas_LINUX.a 复制到系统库中
# 也可以使用 ln -sf ~/BLAS-3.10.0/blas_LINUX.a /usr/local/lib/libblas.a 进行链接

 

安装 cblas

wget http://www.netlib.org/blas/blast-forum/cblas.tgz
tar -xavf cblas.tgz
cd CBLAS
cp Makefile.LINUX Makefile.in
# 修改 Makefile.in ，指出 libblas.a 的位置
# BLLIB = /usr/local/lib/libblas.a
make
sudo cp lib/cblas_LINUX.a /usr/local/lib/libcblas.a
# 将 cblas_LINUX.a 复制到系统库中

 

安装 LAPACK

由于 github 官网链接下载太慢，直接在 http://www.netlib.org/lapack/ 中寻找最新版本下载

tar -xzvf lapack-3.10.0.tar.gz 
cd lapack-3.10.0 
cp make.inc.example make.inc 
# 修改其中 BLASLIB 和 CBLASLIB 路径 
make lib
sudo cp liblapack.a /usr/local/lib/ 
sudo cp libtmglib.a /usr/local/lib/
# 将 liblapack.a 、libtmglib.a 复制到系统库中

cd TESTING
make # 编译 lapack 文件
cd LAPACKE # 进入 LAPACKE 文件夹
make # 编译 lapacke
cp include/*.h /usr/local/include/
# 复制全部头文件到系统头文件目录
cp .. # 返回上一级
cp *.a /usr/local/lib/ 
# 复制全部库文件到系统库目录

 

基本框架

概述

LAPACK API 支持两种形式：标准的 ANSI C 和标准的 FORTRAN

每个 LAPACK 例程都有四个形式：

精度	例程前缀
REAL	S
REAL DOUBLE	D
COMPLEX	C
COMPLEX DOUBLE	Z

 

函数命名规则

• LAPACK 中的每个函数名已经说明了该函数的使用规则
• 所有函数都是以 XYYZZZ 的形式命名
• 对于某些函数，没有第六个字符，只是 XYYZZ 的形式
• X 代表数据类型(S D C Z)，YY 代表数组的类型，ZZZ代表计算方法

注意：在新版中，使用重复迭代法的函数 DSGESV 和 ZCDESV 头两个字母表示使用的精度

DS 输入双精度，算法使用单精度

ZC 输入使用 DOUBLE COMPLEX，算法使用 COMPLEX 单精度

 

几种常用的数组类型

记号	类型
DI (diagonal)	对角阵
GB (general band)	一般带状矩阵
GE (general)	一般矩阵
GT (general tridiagonal)	一般三对角阵
OR (real orthogonal)	实正交阵
SB (real symmetric)	实对称带状阵
ST (real symmetric tridiagonal)	实对称三对角阵
SY (symmetric)	对称阵
TB (triangularband)	三角形带状矩阵

 

编译指令

使用 gfortran 编译，通过-l导入静态库；导入静态库的顺序不能变

gfortran test.cpp -o test -llapacke -llapack -lblas

注意：此时程序必须使用 C 语言编写

 

• 链接 C++ 库

通过添加参数来调用 C++ 相关的库和使用新标准

gfortran test.cpp -o test -llapacke -llapack -lblas -lstdc++ -std=c++11 # 链接 C++ 标准库 + C++11 标准

 

• 链接 fortran 库

如果使用 g++ 进行编译，则需要链接 fortran 库

g++ test.cpp -o test -llapacke -llapack -lblas -lgfortran

另外，需要注意的是，由于高版本的 gcc 默认会在编译参数中添加 -fPIC 参数以实现相对路径的共享，因此可能导致编译时不能正确链接库，这时就需要添加 -no-pie 参数来取消 -fPIC 参数

g++ test.cpp -o test -llapacke -llapack -lblas -lgfortran -no-pie

 

• 使用 extern 修饰

extern void c_func1(float*, int);  	 // 使用 C++ 的编译修饰规则，Fortran 无法调用
extern "C" int c_func2();            // 使用 C 的编译修饰规则，Fortran 可以调用

 

LAPACK 语法

Lapack 函数手册：

http://www.netlib.org/lapack/single/

http://www.netlib.org/lapack/double/

http://www.netlib.org/lapack/complex/

http://www.netlib.org/lapack/complex16

通过 XYYZZZ 中部分类型的组合，我们可以得到不同的函数 LAPACK_XYYZZZ

区分例程

// Linear system, solve Ax = b
LAPACKE_XYYsv(); 	// 标准求解
LAPACKE_XYYsvx(); 	// 精确求解：包括 A'x = b，条件数，误差分析等

// Linear least squares, minimize ||b - Ax||2
LAPACKE_XYYls(); 	// A满秩，QR求解

 

参数格式

// 标准求解 Ax = B A = N x N, B = N x NRHS, X = N x NRHS
lapack_int LAPACKE_dgesv( 
    int matrix_layout, 		 // 矩阵的格式
    lapack_int n,			// 方程组的行数，也即矩阵 A 的行数
    lapack_int nrhs,		// rhs: right-hand-side 右边矩阵的列数，也即 B 的列数
                          
    double* a, 				// 矩阵 A
    lapack_int lda, 		// 数组 A 的主尺寸，这是存放矩阵 A 的数组的尺寸，不小于 N
    lapack_int* ipiv,		// 长为 N 的数组，用于定义置换矩阵 P，一般初始化为0即可
                          
    double* b, 				// 矩阵 B
    lapack_int ldb 			// 数组 B 的主尺寸，这是存放矩阵 B 的数组的尺寸，不小于 N
);

matrix_layout

• LAPACK_ROW_MAJOR 按行求解（标准）
• LAPACK_COL_MAJOR 按列求解

之后介绍的 Lapack 函数中 matrix_layout 都会沿用这两个宏。

 

LAPACK 函数

degsv

求解一般实线性方程组 \(AX=B\) ，其中 \(A\in\mathbb{R}^{N\times N},\ X,B\in\mathbb{R}^{N\times NRHS}\) ，并且要求 \(A\) 可逆。求解过程使用列主元 \(LU\) 分解法，

\[A = PLU \]

其中 \(P\) 是置换阵， \(L,U\) 分别为上下三角阵。

lapack_int LAPACKE_dgesv( int matrix_layout, lapack_int n, lapack_int nrhs,
                          double* a, lapack_int lda, lapack_int* ipiv,
                          double* b, lapack_int ldb );

参数：

• N - 线性方程的个数，也就是 \(A\) 的阶数
• NRHS - 矩阵 \(B\) 的列数
• A - 矩阵 \(A\) ，返回储存 \(A\) 的 \(LU\) 分解
• LDA - 数组 \(A\) 的行维数， \(LDA \ge \max(1,N)\)
• IPIV - 存放返回维数为 \(N\) 的置换向量，行 \(i\) 被替换为行 \(IPIV(i)\)
• B - 矩阵 \(B\)
• LDB - 数组 \(B\) 的行维数， \(LDB\ge \max(1,N)\) ，注意是数组 \(B\) ，如果是单列向量，行数是 1

返回 INFO 存放计算标识：

• 0 成功退出
• < 0 若 INFO = -i ，则第 \(i\) 个变量是一个不可接受的值
• > 0 若 INFO = i ，则 \(U(i,i)\) 为零，因式分解完成，但 \(U\) 不可逆

 

dgels

求解超定或欠定实线性方程组，涉及 \(A\in\mathbb{R}^{M\times N}\) 或者它的转置，使用 \(QR\) 或 \(LQ\) 分解，它假定 \(A\) 满秩。

lapack_int LAPACKE_dgels( int matrix_layout, char trans, lapack_int m,
                          lapack_int n, lapack_int nrhs, double* a,
                          lapack_int lda, double* b, lapack_int ldb );

可选参数 TRANS ：

• 若 TRANS = 'N', m>= n ：求超定系统 \(\|B-AX\|\) 的最小解
• 若 TRANS = 'N', m<n ：求欠定系统的最小解
• 若 TRANS = 'T', m>= n ：求欠定系统 \(A^TX=B\) 的最小解
• 若 TRANS = 'T', m<n ：求超定系统 \(\|B-A^TX\|\) 的最小解

右边向量 \(b\) 和解向量 \(x\) 分别存放为 \(M\times NRHS\) 矩阵 \(B\) 和 \(N\times NRHS\) 矩阵 \(X\) 。

参数：

• TRANS - 如上
• M - 矩阵 \(A\) 的行数
• N - 矩阵 \(A\) 的列数
• NRHS - 矩阵 \(B\) 的列数
• A - 矩阵 \(A\) ，如果 \(M>=N\) ，则 \(A\) 由 \(QR\) 分解覆盖；否则 \(A\) 由 \(LQ\) 分解覆盖
• LDA - 数组 \(A\) 的行维数， \(LDA\ge\max(1,M)\)
• B - 矩阵 \(B\) ，情况比较复杂，只考虑 TRANS = 'N', m>=n ，此时第 1 到 n 行包含最小二乘向量，第 n+1 到 M 行每列的平方和给出残向量 \(B-AX\) 每列的平方和
• LDB - 数组 \(B\) 的行维数， \(LDB\ge \max(1,M,N)\)

 

dsyev

求解实对称矩阵 \(A\) 的所有特征值和对应的特征向量

lapack_int LAPACKE_dsyev( int matrix_layout, char jobz, char uplo, lapack_int n,
                          double* a, lapack_int lda, double* w );

参数：

• JOBZ - 若为 'N' ，表示只计算特征值；若为 'V' 表示计算特征值和特征向量
• UPLO - 若为 'U' ，表示存放 \(A\) 的上三角部分；若为 'L' ，表示存放 \(A\) 的下三角部分
• N - 矩阵 \(A\) 的阶数
• A - 矩阵 \(A\) ，如果 UPLO 为 'U' ，则 \(A\) 的前 \(N\times N\) 上三角部分存放在上三角部分（因为没有必要存放整个 \(A\) ），反之同理；如果 JOBZ 为 'V' ，则 \(A\) 存放正交特征向量；否则根据 UPLO 返回 \(A\) 其部分，其它部分会被摧毁
• LDA - 数组 \(A\) 的行维数， \(LDA\ge\max(1,N)\)
• W - 返回维数为 \(N\) 的向量，升序存放特征值

 

dgees

求解实非对称矩阵 \(A\in\mathbb{R}^{n\times n}\) 的特征值，实 Schur 分解，以及 Schur 向量的矩阵，给出 Schur 分解 \(A = ZTZ^T\) ；也可以选择将对角线上的特征值排序，则 \(Z\) 的第一列就是对应特征值子空间的一个正交基。

lapack_int LAPACKE_dgees( int matrix_layout, char jobvs, char sort,
                          LAPACK_D_SELECT2 select, lapack_int n, double* a,
                          lapack_int lda, lapack_int* sdim, double* wr,
                          double* wi, double* vs, lapack_int ldvs );

参数：

• JOBVS - 若为 'N' ，则不计算 Schur 向量；若为 'V' 则计算
• SORT - 若为 'N' ，则不对对角特征值排序；若为 'S' 则排序
• SELECT - 参数过于复杂，当 SORT 为 'N' ，此参数不被引用
• N - 矩阵 \(A\) 的阶数
• A - 矩阵 \(A\) ，返回 Schur 标准型
• LDA - 数组 \(A\) 的行维数， \(LDA\ge\max(1,N)\)
• SDIM - 返回值，若 SORT 为 'N' ，返回 0
• WR - 返回 \(N\) 维数组
• WI - 返回 \(N\) 维数组，它们分别存放特征值的实部和虚部，共轭特征对会以虚部为正负的顺序连续出现
• VS - 返回 \(LDVS\times N\) 数组，若 JOBVS 为 'V' ，则其中包含 Schur 向量构成的正交阵 \(Z\) ；否则不被引用
• LDVS - 数组 VS 的行维数

 

dgecon

计算一般实矩阵 \(A\) 的条件数

lapack_int LAPACKE_dgecon( int matrix_layout, char norm, lapack_int n,
                           const double* a, lapack_int lda, double anorm,
                           double* rcond );

参数：

• NORM - 若为 '1' 表示 1 范数，为 'I' 表示无穷范数
• N - 矩阵 \(A\) 的阶数
• A - 矩阵 \(A\)
• LDA - 数组 \(A\) 的行维数， \(LDA\ge\max(1,N)\)
• ANORM - 初始矩阵的范数
• RCOND - 返回 RCOND = 1/(norm(A) * norm(inv(A)))

 

dgehrd

通过正交变换 \(Q^TAQ =H\) 将一般实矩阵 \(A\) 化为上 Hessenberg 阵 \(H\) 。

lapack_int LAPACKE_dgehrd( int matrix_layout, lapack_int n, lapack_int ilo,
                           lapack_int ihi, double* a, lapack_int lda,
                           double* tau );

参数：

• N - 矩阵 \(A\) 的阶数
• ILO - 输入整型
• IHI - 输入整型；它们假设矩阵 \(A\) 已经在 1:ILO-1 行列和 IHI+1:N 行列部分上三角化。只有当已经调用了 dgebal 后才会设置它们，否则请将它们分别设为 1 和 N
• A - 矩阵 \(A\) ，返回上三角和次对角线为上 Hessenberg 化矩阵 \(H\) ，剩下的元素和 TAU 一起存放正交阵 \(Q\)
• LDA - 数组 \(A\) 的行维数， \(LDA\ge\max(1,N)\)
• TAU - 返回 \(N-1\) 维数组，存放变换阵的标量因子，具体解释如下：

/*  Further Details
*  ===============
*
*  The matrix Q is represented as a product of (ihi-ilo) elementary
*  reflectors
*
*     Q = H(ilo) H(ilo+1) . . . H(ihi-1).
*
*  Each H(i) has the form
*
*     H(i) = I - tau * v * v**T
*
*  where tau is a real scalar, and v is a real vector with
*  v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
*  exit in A(i+2:ihi,i), and tau in TAU(i).
*
*  The contents of A are illustrated by the following example, with
*  n = 7, ilo = 2 and ihi = 6:
*
*  on entry,                        on exit,
*
*  ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )
*  (     a   a   a   a   a   a )    (      a   h   h   h   h   a )
*  (     a   a   a   a   a   a )    (      h   h   h   h   h   h )
*  (     a   a   a   a   a   a )    (      v2  h   h   h   h   h )
*  (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )
*  (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )
*  (                         a )    (                          a )
*
*  where a denotes an element of the original matrix A, h denotes a
*  modified element of the upper Hessenberg matrix H, and vi denotes an
*  element of the vector defining H(i).
*/

 

dgeqrf

计算一个实 \(M\times N\) 矩阵 \(A\) 的 \(QR\) 分解

lapack_int LAPACKE_dgeqrf( int matrix_layout, lapack_int m, lapack_int n,
                           double* a, lapack_int lda, double* tau );

参数：

• M - 矩阵 \(A\) 的行数
• N - 矩阵 \(A\) 的列数
• A - 矩阵 \(A\) ，返回上三角部分为 \(R\) ，其下的部分和 TAU 存放正交阵 \(Q\)
• LDA - 数组 \(A\) 的行维数， \(LDA\ge\max(1,M)\)
• TAU - 返回标量因子 \(\min(M,N)\) 维向量，具体解释如下：

/*  Further Details
*  ===============
*
*  The matrix Q is represented as a product of elementary reflectors
*
*     Q = H(1) H(2) . . . H(k), where k = min(m,n).
*
*  Each H(i) has the form
*
*     H(i) = I - tau * v * v**T
*
*  where tau is a real scalar, and v is a real vector with
*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
*  and tau in TAU(i).
*/

 

dgetrf

计算一个实 \(M\times N\) 矩阵 \(A\) 的列主元 \(LU\) 分解

lapack_int LAPACKE_dgetrf( int matrix_layout, lapack_int m, lapack_int n,
                           double* a, lapack_int lda, lapack_int* ipiv );

参数：

• M - 矩阵 \(A\) 的行数
• N - 矩阵 \(A\) 的列数
• A - 矩阵 \(A\) ，返回存储 \(L,U\)
• LDA - 数组 \(A\) 的行维数， \(LDA\ge\max(1,M)\)
• IPIV - 返回 \(\min(M,N)\) 维向量，存放置换向量，行 \(i\) 被替换为行 \(IPIV(i)\)

 

dgetri

利用 \(LU\) 分解计算一个矩阵的逆，它先求 \(U\) 的逆，然后求解 \(A^{-1}L = U^{-1}\) 得到 \(A^{-1}\)

lapack_int LAPACKE_dgetri( int matrix_layout, lapack_int n, double* a,
                           lapack_int lda, const lapack_int* ipiv );

参数：

• N - 矩阵 \(A\) 的阶数
• A - 矩阵 \(A\) ，
• LDA - 数组 \(A\) 的行维数， \(LDA\ge\max(1,N)\)
• IPIV - 存放置换向量的 \(N\) 维数组

 

LAPACK 实例

使用标准函数求解线性方程组

利用 LAPACK_dgels求解最小二乘问题

#include <stdio.h>
#include <lapacke.h>
 
int main (int argc, const char * argv[])
{
   double a[5*3] = {1,2,3,4,5,1,3,5,2,4,1,4,2,5,3};
   double b[5*2] = {-10,12,14,16,18,-3,14,12,16,16};
   lapack_int info,m,n,lda,ldb,nrhs;
   int i,j;
 
   m = 5;
   n = 3;
   nrhs = 2;
   lda = 5;
   ldb = 5;
 
   info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',m,n,nrhs,a,lda,b,ldb);
 
   for(i=0;i<n;i++)
   {
      for(j=0;j<nrhs;j++)
      {
         printf("%lf ",b[i+ldb*j]);
      }
      printf("\n");
   }
   return(info);
}

 

最后附上我自己使用的一个求解器类

#pragma once

#include <initializer_list>
#include <iomanip>
#include <iostream>
#include <lapacke.h>

struct lapackSolver
{
    static lapack_int solve(int Row, int Col, double *A, double *b)
    {
        // 标准求解线性方程组
        lapack_int info;
        lapack_int ipiv[Col];

        info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, Col, 1, A, Col, ipiv, b, 1);
        // print(Col, b);
        return info;
    }

    static lapack_int least_square(int Row, int Col, double *A, double *b)
    {
        // 满秩求解最小二乘问题，使用 QR 分解
        lapack_int info;

        // 得到的结果存放在 b 中 Col x 1
        info = LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', Row, Col, 1, A, Row, b, 1);
        // print(Col, b);
        return info;
    }

    static lapack_int eigen(int Row, int Col, double *A, double *wr, double *wi)
    {
        // 求解实矩阵的特征值和特征向量
        LAPACK_D_SELECT2 select;
        lapack_int info, sdim;
        double vs[Row * Row];
        // A 为 Schur 标准型， wr,wi 分别存放特征值的实部和虚部
        info = LAPACKE_dgees(LAPACK_ROW_MAJOR, 'N', 'N', select, Row, A, Row, &sdim, wr, wi, vs, Row);

        // std::cout << "eigenvalues: " << std::endl;
        // for (int i = 0; i < Row; i++)
        // {
        //     std::cout << wr[i] << " + i " << wi[i] << ", ";
        // }
        // std::cout << std::endl
        //           << "Schur: " << std::endl;
        // print(Row, Col, A);
        return info;
    }

    static lapack_int eigen_sy(int Row, int Col, double *A, double *ev)
    {
        // 求解实对称矩阵的特征值和特征向量
        lapack_int info;
        // A 存放特征向量， ev 存放特征值
        info = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'U', Row, A, Row, ev);
        // std::cout << "eigenvalues: " << std::endl;
        // print(Row, ev);
        // std::cout << "eigenvectors: " << std::endl;
        // print(Row, Col, A);
        return info;
    }

    static lapack_int Hessenberg(int Row, int Col, double *A)
    {
        // 存放系数
        double tau[Row - 1];
        lapack_int info;
        // 计算矩阵的上 Hessenberg 化，返回 A 包含上 Hessenberg 部分，和下面的变换矩阵部分
        info = LAPACKE_dgehrd(LAPACK_ROW_MAJOR, Row, 1, Row, A, Row, tau);
        // std::cout << "Hessnberg: " << std::endl;
        // print(Row, Col, A);
        return info;
    }

    static lapack_int QR(int Row, int Col, double *A)
    {
        // 存放系数
        double tau[Col];
        lapack_int info;
        // 计算矩阵的 QR 分解
        info = LAPACKE_dgeqrf(LAPACK_ROW_MAJOR, Row, Col, A, Row, tau);
        // std::cout << "QR: " << std::endl;
        // print(Row, Col, A);
        return info;
    }

    static lapack_int inv(int Row, int Col, double *A)
    {
        // 求矩阵的逆
        lapack_int info, ipiv[Row];
        info = LAPACKE_dgetri(LAPACK_ROW_MAJOR, Row, A, Row, ipiv);
        // std::cout << "inv: " << std::endl;
        // print(Row, Col, A);
        return info;
    }

    // 输出向量和矩阵
    static void print(int dim, double *vec, int w = 0)
    {
        for (int i = 0; i < dim; i++)
        {
            std::cout << std::setw(w) << vec[i] << " ";
        }
        std::cout << std::endl;
    }

    static void print(int row, int col, double *mat, int w = 12)
    {
        for (int i = 0; i < row; i++)
        {
            print(col, &mat[i], w);
        }
        std::cout << std::endl;
    }
};