lapack 线性方程组求解函数dgesv介绍

1. DGESV函数简介
DGESV函数是lapack函数库求解线性方程组AX=B的函数，处于drive level。利用LU分解进行求解，A=PLU，P为置换矩阵，L为单位下三角矩阵，U为上三角矩阵。

2. dgesv函数参数介绍
a) 形式
DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
输入参数：
N（INT）： 线性方程组的数目，矩阵A的维数；
NRHS（INT）： 矩阵B的列数（一般为1，单列）；
A（Double）： 系数矩阵A，维数（LDA，N），通常LDA=N；
LDA（INT）：Leading Dimension of A, LDA>=max(1,N);
B（Double）： 矩阵B，维数（LDB，NRHS），通常LDB=N，NRHS=1；
LDB（INT）： Leading Dimension of B，LDB>=max(1,N);
输出结果：
A（Double）：存储LU分解的L和U矩阵，L为单位下三角阵且对角线的单位值1没有存储；
IPIV（INT）：维数为N，利用该矩阵的值可以推求置换矩阵P，starting from a unit matrix E, E矩阵的第i行与E矩阵的第IPIV（i）行进行交换。
如IPIV=[2 2 3 4]，则P矩阵的结果为：

B（Double）：存储最终的求解结果X值；

INFO（INT）：求解状态参数。INFO=0，成功执行；INFO<0，则INFO=-i，第i个参数是illegal value；INFO>0, U(i,i)=0, LU分解成功执行，但是U是奇异的，最终结果无法计算；

DGESV函数调用了Lapack Computation level的两个函数DGETRF和DGETRS，前者是进行LU分解的函数，后者是利用分解的结果进行前向和后向推求最终结果X。简单介绍下这两个函数的调用格式：

DGETRF（M，N，A，LDA，IPIV，INFO），矩阵A的维度为（M，N），其它参数同上，返回结果为A（PLU），IPIV以及INFO。

DGETRS（TRANS，N，NRHS，A，LDA，IPIV，B，LDB，INFO），求解AX=B或A'X=B，该函数只有在INFO=0时才可以进行。输入变量'TRANS'为字符型，其意义如下：'N'对应AX=B； ‘T’对应A'X=B; 'C' 对应A'X=B (conjugate transpose)，其它变量的意义同DGESV函数。

3. dgesv函数举例
http://www.nag.co.uk/lapack-ex/node5.html
最终返回结果：
B=X= [ 1.0000 -1.0000 3.0000 -5.0000];
A=PLU=[
5.2500 -2.9500 -0.9500 -3.8000
0.3429 3.8914 2.3757 0.4129
0.3010 -0.4631 -1.5139 0.2948
-0.2114 -0.3299 0.0047 0.1314]

IPIV=[2 2 3 4]; 从一个单位阵开始，单位阵的第一行和第2行交换，第二行和第2行交换，第三行和第3行交换，第四行和第4行交换，

所以：

P=[0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1];

L=[1 0 0 0; 0.3429 1 0 0; 0.3010 -0.4631 1 0; -0.2114 -.3299 0.0047 1];

U=[5.2500 -2.9500 -0.9500 -3.8000; 0 3.8914 2.3757 0.4129; 0 0 -1.5139 0.2948; 0 0 0 0.1314]

经过验证PLU=A；


4. dgesvx函数
dgesvx是expert dgesv函数，该函数可以返回更多的中间结果，如计算矩阵A的条件数，估计解的误差范围等等，相应的，函数也增加了很多参数，所以占用的内存也就更多。具体使用可以参见user book http://www.netlib.org/lapack/lug/node26.html和应用举例 http://www.nag.co.uk/lapack-ex/node6.html#DGESVX