1. 参考来源
参考 lapack 官网上的源码，其中有 dgesv 函数的使用说明：http://www.netlib.org/lapack/lapack-3.1.1/html/dgesv.f.html
其算法就是教科书式的 LU 分解+换行。
2. 测试
然后写一个小代码测试了一下，做一个最简单的方程组：
[egin{bmatrix}
1 & 2 \
0 & 1
end{bmatrix}
egin{bmatrix}
x \
y
end{bmatrix}
=
egin{bmatrix}
5 \
2
end{bmatrix}
Rightarrow
egin{bmatrix}
x \
y
end{bmatrix}
=
egin{bmatrix}
1 \
2
end{bmatrix}
]
#include<iostream>
using namespace std;
extern "C" void dgesv_(int *N, int *NRHS, double *A, int *LDA, int * IPIV, double *B, int *LDB, int * INFO);
/*
* solve A x = B, A and B keep unchanged when exit
*/
void lapack_dgesv( int dim, double ** A, double * B, double * x ){
double * Atemp = new double [ dim * dim ];
for(int i=0;i<dim;i++)
for(int j=0;j<dim;j++) Atemp[ i*dim + j ] = A[j][i];
for(int i=0;i<dim;i++) x[i] = B[i];
int * IPIV = new int [ dim ];
int INFO;
dgesv_( &dim, &dim, Atemp, &dim, IPIV, x, &dim, &INFO );
delete [] Atemp; delete [] IPIV;
if( INFO == 0 ) cout<<" lapack_dgesv exited successfully "<<endl;
else if( INFO > 0 ) cout<<" lapack_dgesv error: A is singular "<<endl;
else cout<<" lapack_dgesv: the "<< -INFO <<"-th argument has an illegal value"<<endl;
}
int main(){
int dim = 2;
double Aarray[] = {1,2,0,1};
double ** A = new double * [ dim ];
for(int i=0;i<dim;i++){
A[i] = new double [ dim ];
for(int j=0;j<dim;j++) A[i][j] = Aarray[i*dim+j];
}
double * B = new double [dim]; B[0]=5; B[1]=2;
double * x = new double [dim];
lapack_dgesv( dim, A, B, x );
cout<<" The solution is : "; for(int i=0;i<dim;i++)cout<<x[i]<<","; cout<<endl;
return 0;
}
运行结果：
luyi@Swagger:~/test/dgesv$ g++ main.cpp -llapack
luyi@Swagger:~/test/dgesv$ ./a.out
lapack_dgesv exited successfully
The solution is : 1,2,
luyi@Swagger:~/test/dgesv$
3. 应用
看起来代码是对的，那就可以把代码中的 extern 声明，以及封装的 lapack_dsgev 函数放到一个文件里，配上头文件，其他场合都可以使用了。