矩陣的 LU 分解法（LU decomposition）

原文連結：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: -------10.333333        10.666667        0.2     1------U is: -------3       1       21.66667 -1.666671Press any key to continue . . .

原文連結：http://blog.csdn.net/zhangxaochen/article/details/8021375

{{OVER}}