求解实数线性方程组：lapack dgesv

1. 参考来源

参考 lapack 官网上的源码，其中有 dgesv 函数的使用说明：http://www.netlib.org/lapack/lapack-3.1.1/html/dgesv.f.html
其算法就是教科书式的 LU 分解+换行。

2. 测试

然后写一个小代码测试了一下，做一个最简单的方程组：

\[\begin{bmatrix} 1 & 2 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 5 \\ 2 \end{bmatrix} \Rightarrow \begin{bmatrix} x \\ y \end{bmatrix} = \begin{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 函数放到一个文件里，配上头文件，其他场合都可以使用了。