使用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;
}
输出结果例如以下:
watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvZ2lzZWxpdGU=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center" alt="">
也能够使用以下的方法:
#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
最新文章
- 0929mysql前缀索引如何找到合适的位数
- map的应用
- Android_ADB 常用 shell命令 和 sqlite3 简单增删改查
- 减小Delphi XE5编译出来的程序体积
- MongoDB (四) MongoDB 数据模型
- 游戏开发Camera之Cinematic Camera-深度
- vs在线工具杂烩
- Vijos P1881 闪烁的星星 (加强自己多一点。。)
- TaintDroid:智能手机监控实时隐私信息流跟踪系统(一)
- Ultimate thread group线程组和Stepping thread group线程组测试场景
- 201521123114 《Java程序设计》第6周学习总结
- 【读书笔记】【深入理解ES6】#10-改进的数组功能
- [USACO12JAN]爬山Mountain Climbing
- ThinkPHP 初探
- GT-随身调详细教程
- [原创] 扩展jquery-treegrid插件, 实现勾选功能和全删按钮.
- josn的格式化
- CF 1042F
- 【NPM】常见问题解决
- python----线程进程协程
热门文章
- cocos2dx--vs2012+lua开发环境搭建
- oracle expdp 备份脚本
- PAL相机
- Spring Boot 热部署(转)
- 【例题3-3 UVA - 401】Palindromes
- 怎么不让控制台system.out.println()打印
- 【JAVA编码专题】总结 分类: B1_JAVA 2015-02-11 15:11 290人阅读 评论(0) 收藏
- Android开发之SpannableString具体解释
- 设置secureCRT中vim的字体颜色 分类: B3_LINUX 2014-07-12 22:01 1573人阅读 评论(0) 收藏
- android,安卓get请求的提交以及我遇到的异常