使用lapack库求逆矩阵

来源：互联网发布：js如何取整数编辑：程序博客网时间：2024/05/21 06:46

废话不多说，直接上代码：

[cpp] view plaincopy
#include <string>  
#include "lapacke.h"  
#include "lapack_aux.h"  
  
int main(int argc,char** argv)  
{  
    setlocale(LC_ALL,"");  
    double a[] =   
    {  
        3,-1,-1,  
        4,-2,-1,  
        -3,2,1  
    };  
    int m = 3;  
    int n = 3;  
    int lda = 3;  
    int ipiv[3];  
    int info;  
    print_matrix("a",m,n,a,lda);      
    info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR,m,n,a,lda,ipiv);  
    print_matrix("a",m,n,a,lda);  
  
    info = LAPACKE_dgetri(LAPACK_ROW_MAJOR,m,a,lda,ipiv);  
    print_matrix("a",m,n,a,lda);  
    return 0;  
}  

输出结果如下：

也可以使用下面的方法：

[cpp] view plaincopy
#include <string>  
#include "lapacke.h"  
#include "lapack_aux.h"  
  
int main(int argc,char** argv)  
{  
    setlocale(LC_ALL,"");  
    double a[] =   
    {  
        3,4,-3,  
        -1,-2,2,  
        -1,-1,1  
    };  
  
    int m = 3;  
    int n = 3;  
    int lda = 3;  
    int ipiv[3];  
    int info;  
    print_matrix("a",m,n,a,lda);      
    LAPACK_dgetrf(&m,&n,a,&lda,ipiv,&info);  
    print_matrix("a",m,n,a,lda);  
  
    double *b = new double[m]();  
    //求普通矩阵的逆矩阵  
    LAPACK_dgetri(&m,a,&lda,ipiv,b,&n,&info);  
    print_matrix("inv(a)",m,n,a,lda);  
    return 0;  
}  

输出结果如下：

这种方法的好处在于API接口的定义和对应的Fortran接口一致，比如dgetri,我们可以在双精度的函数说明(http://www.netlib.org/lapack/double/)文档中找到dgetri.f，打开这个Fortran文件，就可以知道相应的参数的含义了。

不过这里要注意存储矩阵时，两种方法之间的区别。第一种方法中，我们可以通过主序告诉lapack的接口我们的矩阵是以行为主序的，也就是在数组中，这个矩阵是按行存储的，对于一个3x3矩阵，输入的9个元素，前3个数是矩阵的第一行，紧接着是矩阵的第二行，最后是矩阵的第三行，而第二种方法中，没有主序这个参数，研究发现，Fortran默认是以列为主序的，也就是说我们在用数组输入一个3x3矩阵时，前3个数表示第1列，再3个数为第2列，最后3个数为第3列。因此在给定矩阵的时候，我们需要按列输入。因此方法2中的数组a,以列为主序，表示的矩阵实际上是这样的：

[plain] view plaincopy
 3    -1    -1  
 4    -2    -1  
-3     2     1  

这相当于把第一种方法中的主序改为LAPACK_COL_MAJOR,如下：

[cpp] view plaincopy
#include <string>  
#include "lapacke.h"  
#include "lapack_aux.h"  
  
int main(int argc,char** argv)  
{  
    setlocale(LC_ALL,"");  
    double a[] =   
    {  
        3,4,-3,  
        -1,-2,2,  
        -1,-1,1  
    };  
    int m = 3;  
    int n = 3;  
    int lda = 3;  
    int ipiv[3];  
    int info;  
    print_matrix("a",m,n,a,lda);      
    info = LAPACKE_dgetrf(LAPACK_COL_MAJOR,m,n,a,lda,ipiv);  
    print_matrix("a",m,n,a,lda);  
  
    info = LAPACKE_dgetri(LAPACK_COL_MAJOR,m,a,lda,ipiv);  
    print_matrix("a",m,n,a,lda);  
    return 0;  
}  

最后，我们在Matlab中验证一下，如下：

[plain] view plaincopy
>>  a = [3,4,-3,-1,-2,2,-1,-1,1]  
  
a =  
  
     3     4    -3    -1    -2     2    -1    -1     1  
  
>>  a = reshape(a,3,3)  
  
a =  
  
     3    -1    -1  
     4    -2    -1  
    -3     2     1  
  
>> inv(a)  
  
ans =  
  
     0     1     1  
     1     0     1  
    -2     3     2  

可见我们的计算结果好Matlab的结果一致。

附辅助函数：

[cpp] view plaincopy
#include <stdio.h>  
#include "lapack_aux.h"  
  
/* Auxiliary routine: printing a matrix */  
void print_matrix( char* desc, lapack_int m, lapack_int n, double* a, lapack_int lda )  
{  
    lapack_int i, j;  
    printf( "\n %s\n", desc );  
    for( i = 0; i < m; i++ )  
    {  
        for( j = 0; j < n; j++ ) printf( " %6.2f", a[i*lda+j] );  
        printf( "\n" );  
    }  
}  

参考文件：

http://blog.csdn.net/kevinzhangyang/article/details/6859246
http://blog.csdn.net/daiyuchao/article/details/2026173
http://blog.csdn.net/daiyuchao/article/details/2026162
http://www.cnblogs.com/xunxun1982/archive/2010/05/12/1734001.html
http://www.cnblogs.com/xunxun1982/archive/2010/05/13/1734809.html
http://hi.baidu.com/data2009/item/50bce0704cf57a14d0dcb3e8
http://blog.sina.com.cn/s/blog_40b056950100htpt.html
http://blog.csdn.net/cleverysm/article/details/1925553
http://blog.csdn.net/cleverysm/article/details/1925549
http://www.cnblogs.com/Jedimaster/archive/2008/06/22/1227656.html

0 0