原文链接:http://blog.csdn.net/zhangxaochen/article/details/8021375
看这里: http://is.gd/VoBVUJ
理论随便搜搜。。
分解算法如下:(其中 alpha 代表 L矩阵的元素,beta代表U矩阵的元素)
我的实现如下:(传入一个矩阵 a,函数结束时 a的数据已被破坏,变成同时存储了 L, U 两个三角阵)
void myLUdcmp(vector<vector<double>>& a){ //假定传入的vector的每行都是对齐的,不存在溢出(因为本来 vector<vector<T>> 这样的结构很可能行之间参差不齐) //假定 L矩阵的主对角线人为置为1 //列循环: for(size_t j=0; j<a[0].size(); j++){ for(size_t i=0; i<a.size(); i++){ if(i<=j){ //a[i][j] 即 beta(ij) double sum=0; for(int k=0; k<i; k++) sum+=a[i][k]*a[k][j]; a[i][j]-=sum; } if(i>j){ //a[i][j] 即 alpha(ij) double sum=0; for(int k=0; k<j; k++) sum+=a[i][k]*a[k][j]; a[i][j]=(a[i][j]-sum)/a[j][j]; } }//for i }// for j }
测试代码:
#include <vector> #include <iomanip> void myLUdcmp(vector<vector<double>>& a); void main(){ double arrA[]={3, 1, 2, 1, 2, -1, 2, 1, 2,}; vector<vector<double>> aa(3, vector<double>(3)); for(int i=0; i<3; i++){ for(int j=0; j<3; j++){ aa[i][j]=arrA[i*3+j]; } } myLUdcmp(aa); cout<<"------L is: -------\n"; for(int i=0; i<aa.size(); i++){ for(int j=0; j<i; j++){ cout<<aa[i][j]<<'\t'; } cout<<1<<endl; } cout<<endl; cout<<"------U is: -------\n"; for(int i=0; i<aa.size(); i++){ for(int j=i; j<aa[0].size(); j++){ cout<<aa[i][j]<<'\t'; } cout<<endl; } }
输出:
------L is: ------- 1 0.333333 1 0.666667 0.2 1 ------U is: ------- 3 1 2 1.66667 -1.66667 1 Press any key to continue . . .
原文链接:http://blog.csdn.net/zhangxaochen/article/details/8021375
{{OVER}}