src/mathlib_src/GLU_Composition.cpp

説明を見る。

#include "GLU_Composition.h"
#include "GMatrix_Util.h"

#include <cmath>

using namespace GMathLib;
using namespace GMathLib::GMatrix_Util;

GLU_Composition::GLU_Composition()
: GObject(), combined_flag(GLU_Composition::LU_NOT_COMBINED), p_L(0), p_U(0), p_LU(0), p_P(0)
{
    Class_Name("GLU_Composition");
}

GLU_Composition::~GLU_Composition()
{
    // 保持している行列オブジェクトの破棄
    delete_mx_object();
}

void GLU_Composition::delete_mx_object()
{
    if(combined_flag == GLU_Composition::LU_NOT_COMBINED){
        if(p_L != 0) delete p_L;
        if(p_U != 0) delete p_U;
    }else if(combined_flag == GLU_Composition::LU_COMBINED){
        if(p_LU != 0) delete p_LU;
    }

    if(p_P != 0){
        delete p_P;
    }
}

bool GLU_Composition::LU_Composition(GMatrix& A_mx, int lu_combined_flag)
{
    
    // 与えられた行列が非正方行列の場合はエラー出力して偽を返す
    int n = A_mx.Row();
    if( n != A_mx.Column() ){
        const std::string context = "A given matrix is not a square matrix.";
        IO::GError_Output::Puts(Class_Name() + "::LU_Composition", context);
        
        return false;
    }
    
    // L, U 保持用の行列オブジェクトが既に生成されていた場合は破棄する. 
    // LU の結合フラグを更新する
    delete_mx_object();
    
    combined_flag = lu_combined_flag;
    
    // 正方行列 A と同じサイズの 行列オブジェクト L, U, P を生成
    // (GLU_Composition::LU_COMBINED が指定されている場合は LU を生成)
    if(combined_flag == GLU_Composition::LU_NOT_COMBINED){
        p_L = new GMatrix(n, n);
        p_U = new GMatrix(n, n);

        GMatrix_Util::Identity_Matrix(*p_L);
        p_U->Fill(0.0);
        
    }else if(combined_flag == GLU_Composition::LU_COMBINED){
        p_LU = new GMatrix(n, n);
    }

    p_P = new GMatrix(n, n);
    p_P->Fill(0.0);

    // lu_composition() を呼び出し, LU 分解をおこなう
    GMatrix A_lu(n, n);
    GVector pvot_vec(n);
    A_lu.Copy(A_mx);

    if(! lu_composition(A_lu, pvot_vec)){
        const std::string context = "LU composition has failed.";
        IO::GError_Output::Puts(Class_Name() + "::LU_Composition", context);
        return false;
    }

    // LU 結合形式で返ってきた LU 分解の結果を行列オブジェクト L, U に分けるか, 
    // 行列オブジェクト LU にコピーする
    if(combined_flag == GLU_Composition::LU_NOT_COMBINED){
        for(int i=0; i < n; i++){
            for(int j=0; j < i; j++){
                p_L->Set(i, j, A_lu.Get(i,j));
            }
        }

        for(int i=0; i < n; i++){
            for(int j=i; j < n; j++){
                p_U->Set(i, j, A_lu.Get(i, j));
            }
        }

    }else if(combined_flag == GLU_Composition::LU_COMBINED){
        p_LU->Copy(A_lu);
    }

    // 部分 Pivot 選択による行交換を表す置換行列を作成する
    // GVector オブジェクト pvot_vec から, GMatrix オブジェクトに作り替える. 
    for(int i=0; i < n; i++){
        p_P->Set(i, pvot_vec.Get(i), 1.0);
    }


    return true;
}

/*
 * 実際に Court's algorithm を使って LU 分解を行う関数 
 * (行列のスケーリング・部分ピボット選択あり) 
 * NUMERICAL RECIPES IN FORTRAN 77 : THE ART OF SCIENTIFIC COMPUTIONG (p 34- )のコードを利用.
 */
bool GLU_Composition::lu_composition(GMatrix& A, GVector& pv){
    int N = A.Row();
    double aamax;
    double TINY = 1.0e-20;

    // 各行の暗黙的なスケーリングを格納するベクトル
    GVector sc_v(N); 

    // 各行のスケーリング値を求める
    for(int i=0; i < N; i++){
        pv(i) = i; // ついでにピボットの初期化
        aamax = 0.0;
        
        for(int j=0; j < N; j++){
            if(std::fabs(A(i, j)) > aamax){
                aamax = std::fabs(A(i, j));
            }

            if( aamax == 0.0 ){
                const std::string context = "A given matrix is singular matrix.";
                IO::GError_Output::Puts(Class_Name() + "::lu_composition", context);
                return false;
            }
            
            // スケーリングを保存する
            sc_v(i) = 1.0 / aamax;
        }     
    }

    // クラウト法における列ループ
    double sum;
    double dum;
    int imax=0;

    for(int j=0; j < N; j++){
        
        //i=j を除く, U_ij=A_ij - int_{k=1}^{i-1} Lik Ukj (i=1,..,j) の計算
        for(int i=0; i <= j-1; i++){
            sum = A(i, j);
            for(int k=0; k <= i-1; k++){
                sum += - A(i, k) * A(k, j); 
            }
            A(i, j) = sum;
        }

        //最大のピボット要素の探索のための初期化
        aamax = 0.0;
        //i=j :        U_ij=A_ij - int_{k=1}^{i-1} Lik Ukj (i=1,..,j) の計算
        //i=j+1,..,N : U_ij= (A_ij - int_{k=1}^{i-1} Lik Ukj)/Ujj (i=1,..,j) の計算
        for(int i=j; i < N; i++){

            sum = A(i, j);
            for(int k=0; k <= j-1; k++){
                sum += - A(i, k) * A(k, j);
            }
            A(i, j) = sum;
            
            // ピボットの効果を計算し, その効果が最大な項を探す
            dum = sc_v(i) * std::fabs(sum);
            if( dum >= aamax ){
                imax = i;
                aamax = dum;
            }
        }

        // 行の相互交換が必要な場合は行う
        if( j != imax ){
            for(int k=0; k < N; k++){
                dum = A(imax, k);
                A(imax, k) = A(j, k);
                A(j, k) = dum;
            }

            // スケールファクタも入れ替える
            sc_v(imax) =  sc_v(j);

            // ピボット配列も入れ替える
            dum = pv(imax);
            pv(imax) = pv(j);
            pv(j) = dum;
        }

        // ピポットがゼロならば, 行列は singular である. 
        // ひとまず TINY をいれて回避する
        if( A(j, j) == 0.0 )
            A(j, j) = TINY;
        
        // 最後にピボット要素で割る
        if( j != N-1){
            dum = 1.0 / A(j, j);
            for(int i=j+1; i < N; i++){
                A(i, j) = A(i, j) * dum;
            }
        }
    }
    
    return true;
}

---