C++ LU Decomposition

2010/11/23
(この記事の文字数: 1470)

LU分解のプログラムをC++で作ってみました。
行列を三角行列に分解することで演算効率をあげようだなんてよく考えたものです。あまり大規模な行列演算をやった経験がないので、どれだけ効率が上がるのか実感はありませんが、LU分解はよく使われるようです。

I made LU decomposition program in C++. I have no experience for computing large amout of matrices, so I don't know how fast it work. Anyway, it is a very common algorithm in matrix computing.

以下がソースコードです。

Here's a source code.


void LUDecomposition(double **L, double **U,

                 double **mat, int row_size)

{

 int itr_i, itr_j, itr_k;

 int col, row;


 for(itr_i=0; itr_i<row_size; itr_i++){

  for(itr_j=0; itr_j<itr_i+1; itr_j++)

  {

   row = itr_i;

   col = itr_j;

   L[row][col] = mat[row][col];

   for(itr_k=0; itr_k<col; itr_k++){

    L[row][col] -= L[row][itr_k] * U[itr_k][col];

   }

   row = itr_j;

   col = itr_i;

   U[row][col] = mat[row][col];

   for(itr_k=0; itr_k<row; itr_k++){

    U[row][col] -= L[row][itr_k] * U[itr_k][col];

   }

   U[row][col] /= L[row][row];

  }

  for(itr_j=itr_i+1; itr_j<row_size; itr_j++){

   L[itr_i][itr_j]=0.0;

   U[itr_j][itr_i]=0.0;

  }

 }

}

Example

#include<iostream>

using namespace std;

void main(){
    double **L, **U, **mat;
    int i, j;
    mat = new double*[3];
    L = new double*[3];
    U = new double*[3];
    for(i=0; i<3; i++){
        mat[i] = new double[4];
        L[i] = new double[3];
        U[i] = new double[3];
    }
    mat[0][0] = 2; mat[0][1] = 3; mat[0][2] = 4; mat[0][3] = 6;
    mat[1][0] = 3; mat[1][1] = 5; mat[1][2] = 2; mat[1][3] = 5;
    mat[2][0] = 4; mat[2][1] = 3; mat[2][2] = 30;mat[2][3] = 32;
    LUDecomposition(L, U, mat, 3);
    cout<<"L"<<endl;
    for(i=0; i<3; i++){
        for(j=0; j<3; j++){
            cout<<setw(5)<<L[i][j];
        }
        cout<<endl;
    }
    cout<<"U"<<endl;
    for(i=0; i<3; i++){
        for(j=0; j<3; j++){
            cout<<setw(5)<<U[i][j];
        }
        cout<<endl;
    }
}