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;
    }
};