LAPACK(3)——线性方程组求解

解线性方程组A‘x=B使用函数LAPACKE_dgesv，函数原型如下，

lapack_int LAPACKE_dgesv( int matrix_order, lapack_int n, lapack_int nrhs,
                          double* a, lapack_int lda, lapack_int* ipiv,
                          double* b, lapack_int ldb );
/*
//参数说明：
matrix_order    矩阵的顺序，两种取值，LAPACK_COL_MAJOR或是LAPACK_ROW_MAJOR。
　　　　　　　　　 　　如果是COL_MAJOR，表示a[0:N-1]是第一列，a[N:2N-1]是第二列，表明矩阵的元素存储顺序。
　　　　　　　　　　　LAPACK的优先选择是COL_MAJOR，当是ROW_MAJOR时，会先进行转置，然后进行计算，然后转置回来；
n        　　　　矩阵的size，系数矩阵应该是方阵，一个5x5的矩阵，则n=5；
nrhs        　　rhs的列数，LAPACK可以同时对一个系数矩阵，多个rhs进行求解，
　　　　　　　　　　　因为对系数矩阵只需要进行一次LU分解，多个rhs一起求解更方便；
a        　　　　系数矩阵；
lda        　　 leading dimension of a，lda>=max(1,n)；
ipiv        　　pivot的行交换记录，第i行被交换到第ipiv[i]行；
b        　　　 rhs，可以是多列。其元素存储顺序是
ldb        　　 leading dimension of b, ldb>=max(1,n)，
　　　　　　　　　　　本来看代码中应该是ldb>=max(1,nrhs)的，ldb=nrhs的时候报错了；
// 返回值：
info = 0, 程序正常运行结束。
*/

测试代码如下，参考 http://blog.sina.com.cn/s/blog_441ac3c9010000v0.html

#include <stdio.h>

//lapacke headers
#include "lapacke.h"
#include "lapacke_config.h"
#include "lapacke_utils.h"

extern lapack_int LAPACKE_dgesv( int matrix_order, lapack_int n, lapack_int nrhs,
                          double* a, lapack_int lda, lapack_int* ipiv,
                          double* b, lapack_int ldb );

int main()
{
    int N = 4;
    double A[16] = {  1,  2,  3,  1,
                      4,  2,  0,  2,
                     -2,  0, -1,  2,
                      3,  4,  2, -3};
    double B[8] = {  6,  2,  1,  8,
                     1,  2,  3,  4};
    int ipiv[4];
    int n = N;
    int nrhs = 2;
    int lda = N;
    int ldb = N;

    int info = LAPACKE_dgesv(LAPACK_COL_MAJOR,n,nrhs,A,lda,ipiv,B,ldb);
    printf("info:%d\n",info);
    if(info==0)
    {
        int i = 0;
        int j = 0;
        for(j=0;j<nrhs;j++)
        {
            printf("x%d\n",j);
            for(i=0;i<N;i++)
                printf("%.6g \t",B[i+j*N]);
            printf("\n");
        }
    }
}

程序运行输出结果如下，

info:0
x0
1     2     0     -1     
x1
1.375     0.375     0.375     -0.375 
表示的是
A = 1   4    -2    3
    2   2     0    4
    3   0    -1    2
    1   2     2   -3
B = 6   1
    2   2 
    1   3
    8   4
x = 1  1.375
    2  0.375
    0  0.375
   -1 -0.375

网上可以看到的LAPACKE的资料好像比较少，LAPACKE只是在LAPACK的基础上进行封装，接口参数有些不同，调用LAPACK函数之前进行了一系列的参数检测。所以目前可参考的就是LAPACKE的源代码。先查看LAPACKE的源代码，然后查看对应的LAPACK的函数参数，基本上就可以正确调用这些函数了。LAPACK的函数说明可以参考下面的网站 http://www.oschina.net/code/explore/OpenCV-2.1.0/3rdparty/lapack。OpenCV中使用的是CLAPACK包，CLAPACK的函数参数与LAPACK一致，不像LAPACKE，自己添加参数。