使用lapack图书馆逆矩阵

阿土，直接在代码：

#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;
}

输出结果例如以下：



也能够使用以下的方法：

#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,以列为主序，表示的矩阵实际上是这种：

3    -1    -1
4    -2    -1
-3     2     1

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

#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中验证一下，例如以下：

>>  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的结果一致。

附辅助函数：

#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