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;
