其他分享
首页 > 其他分享> > Lapack 科学计算包

Lapack 科学计算包

作者：互联网

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

 

环境包

需要安装 gfortran 和 cmake

sudo apt-get install gfortran

 

BLAS 库和 CBLAS 接口

BLAS（basic linear algebra subroutine） 是一系列基本线性代数运算函数的接口（interface）标准。这里的线性代数运算是指例如矢量的线性组合，矩阵乘以矢量，矩阵乘以矩阵等。接口在这里指的是诸如哪个函数名实现什么功能，有几个输入和输出变量，分别是什么。

BLAS 被广泛用于科学计算和工业界，已成为业界标准。在更高级的语言和库中，即使我们不直接使用 BLAS 接口，它们也是通过调用 BLAS 来实现的（如 Matlab 中的各种矩阵运算）。

BLAS 原本是用 Fortran 语言写的，但后来也产生了 C 语言的版本 cBLAS ， 接口与 Fortran 的略有不同（例如使用指针传递数组），但大同小异。

注意 BLAS 是一个接口的标准而不是某种具体实现（implementation），简单来说，就是不同的作者可以各自写出不同版本的 BLAS 库， 实现同样的接口和功能，但每个函数内部的算法可以不同。这些不同导致了不同版本的 BLAS 在不同机器上运行的速度也不同。

BLAS 的官网是 Netlib ，可以浏览完整的说明文档以及下载源代码。这个版本的 BLAS 被称为 reference BLAS，运行速度较慢，通常被其他版本用于衡量性能。对于 Intel CPU 的计算机，性能最高的是 Intel 的 MKL （Math Kernel Library） 中提供的 BLAS 。 MKL 虽然不是一个开源软件，但目前可以免费下载使用。如果想要免费开源的版本，可以尝试 OpenBlas 或者 ATLAS2 ，另外，无论是否使用 MKL，BLAS 的文档都推荐看 MKL 的相关页面 。

 

安装 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

 

函数命名规则

注意：在新版中，使用重复迭代法的函数 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++ 相关的库和使用新标准

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

 

如果使用 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 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 函数中 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 );

参数：

返回 INFO 存放计算标识：

 

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 ：

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

参数：

 

dsyev

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

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

参数：

 

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

参数：

 

dgecon

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

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

参数：

 

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

参数：

/*  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 );

参数：

/*  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 );

参数：

 

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

参数：

 

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

标签：LDA,lapack,Lapack,科学计算,矩阵,int,数组,double
来源： https://www.cnblogs.com/Bluemultipl/p/15915323.html