src/mathlib_src/GLinear_Homo_Eq.cpp

説明を見る。

#include "GLinear_Homo_Eq.h"
#include <iostream>
#include <sstream>
#include <cmath>
#include "../io_src/GError_Output.h"

using namespace GMathLib;
using namespace GMathLib::IO;
using namespace GMathLib::GMatrix_Util;

GLinear_Homo_Eq::GLinear_Homo_Eq()
: GObject(), p_lu(0)
{
    Class_Name("GLinear_Homo_Eq");
}

GLinear_Homo_Eq::~GLinear_Homo_Eq()
{
    if(p_lu != 0){
        delete p_lu;
    }
}

/*
 * �����s�{�b�g�I�����s���K�E�X�̏����@��p�����A���ꎟ�������\���o
 */
void GLinear_Homo_Eq::S_Gauss_Partial_Pivot(const GMatrix& matrix, const GVector& vec, GVector& answer_vec){
        GMatrix matrix_copy = matrix;
        GVector vec_copy = vec;
        
        S_Gauss_Partial_Pivot(matrix_copy, vec_copy, answer_vec);
}

void GLinear_Homo_Eq::S_Gauss_Partial_Pivot(GMatrix& A, GVector& b, GVector& x)
{
        // �^����ꂽ�s���x�N�g���̌`���K�������m�F����
        if( !check_matrix_size(A, b, x) ){
            std::string errorinfo = "Matrix or vector size is invalid !";
            GError_Output::Puts(Class_Name() + "::S_Gauss_Partial_Pivot", errorinfo);
                
            return;
        }
        
        const double TINY = 1.0e-10;
        int n = b.Size();
        int status = 0;
        int pk, max_p, pj;
        int *p = new int[n];
        double max, aik;
        double pivot, multi, sum;
                        
        //�s�{�b�g�ԍ����Z�b�g����
        for(int k = 0; k < n; k++)
            p[k] = k;
                
            for(int k = 0; k < n - 1; k++){
                
            //�v�f�̍ő�l�����ƒs�{�b�g�s�𒲂ׂ�
            pk = p[k]; max = 0.0;       max_p = k;              
            for(int i = k; i < n; i++){
                aik = A(p[i], k);

                if(std::abs(aik) > max){
                    max = std::abs(aik);
                    max_p = i;
                }
            }
                        
            //�ő�l��臒l��菬�����ꍇ�́C���[�v���o��
            if(max <= TINY){
                status = 8; break;
            }
                        
            //�s�{�b�g����
            if(max_p != k){
                p[k] = p[max_p]; p[max_p] = pk; pk = p[k];
            }
                        
            //�O�i����
            pivot = 1.0 / A(pk, k);
            for(int j = k +1; j < n; j++){
                pj = p[j];
                multi = A(pj, k) * pivot;
                                
                if(std::abs(multi) > TINY){
                    for(int i = k + 1; i < n; i++)
                        A(pj, i) = A(pj, i) - multi * A(pk, i);
                                                
                        b(pj) = b(pj) - multi * b(pk);
                }else{
                    A(pj, k) = 0.0;
                }
            }
        }
                
        //�G���[����
                
        //��i�����Ŏg�����鐔��臒l�����������Ȃ������m�F����
        if(std::abs( A(p[n-1], n -1) ) < TINY)
            status = 9;
                
        //��肪����ꍇ�̓G���[�o�͂���
        if(status != 0){
            std::ostringstream oss;
            oss << status;
            std::string errorinfo = "GAUSS, Failed! The system is singular! STATUS : " + oss.str();
            GError_Output::Puts(Class_Name() + "::S_Gauss_Partial_Pivot", errorinfo);
         }else{

            //��i�������n�߂�
            x(n-1) = b(p[n-1]) / A(p[n-1], n-1);
            for(int k = n - 2; k >= 0; k--){
                pk=p[k]; sum = b(pk);
                for(int i = k + 1; i < n; i++)
                    sum = sum - A(pk, i) * x(i);
                                
                x(k) = sum / A(pk, k);  
            }
        }
                
        delete p;                       
}


/*
 * �����s�{�b�g�I������эs��̃X�P�[�����O���s�� LU �����ɂ��A���ꎟ�������\���o
 */
void GLinear_Homo_Eq::S_LU_Partial_Pivot(GMatrix& matrix, GVector& vec, GVector& answer_vec, bool do_lu_cmp_flag)
{
    // �^����ꂽ�s���x�N�g���̌`���K�������m�F����
    if( !check_matrix_size(matrix, vec, answer_vec) ){
        std::string errorinfo = "Matrix or vector size is invalid !";
        GError_Output::Puts(Class_Name() + "::S_LU_Partial_Pivot", errorinfo);
                
        return;
    }
    
    GMatrix* p_Alu = 0;
    GMatrix* p_P = 0;
    int n = matrix.Row();

    // GLU_Composition �̃N���X�I�u�W�F�N�g���܂���������Ă��Ȃ��ꍇ�͂܂���������. 
    if(p_lu == 0){
        p_lu = new GLU_Composition();
    }

    // do_lu_cmp_flag ���^�Ȃ��, LU �����������Ȃ�
    if(do_lu_cmp_flag){
        p_lu->LU_Composition(matrix, GLU_Composition::LU_COMBINED);
    }else{
        if( p_lu->Get_LU() == 0){
            std::string errorinfo = "Specify do_lu_cmp_flag=true !";
            GError_Output::Puts(Class_Name() + "::S_LU_Partial_Pivot", errorinfo);
                
            return;
        }
    }

    // LU �����̌��ʂ���ђu���s��̎擾
    p_Alu = p_lu->Get_LU();
    p_P = p_lu->Get_P(); 
    
    // LU �����̍ۂ̍s�u���������’u���s����g����, 
    // �E�ӂ̗�x�N�g�������ւ���
    GVector b(n);
    b.Multi(*p_P, vec);
    
    // A �� L, U ���g���ĘA���������̉������߂�
    // A x = b, A=LU
    // L y = b, U x = y
    double sum;

    // �܂��O�i�����ɂ�� y �����߂�
    answer_vec(0) = b(0);
    
    for(int i=1; i < n; i++){
        sum = 0.0;
        for(int j=0; j <= i-1; j++){
            sum += (*p_Alu)(i, j) * answer_vec(j); 
        }

        answer_vec(i) = b(i) - sum;
    }
    
    // �����Č�ޏ����ɂ�� x �����߂�
    answer_vec(n-1) = answer_vec(n-1) / (*p_Alu)(n-1, n-1);

    for(int i=n-2; i >= 0; i--){
        sum = 0.0;
        for(int j=i+1; j < n; j++){
            sum += (*p_Alu)(i, j) * answer_vec(j);
        }

        answer_vec(i) = ( answer_vec(i) - sum ) / (*p_Alu)(i, i);
    }
}

bool GLinear_Homo_Eq::check_matrix_size(GMatrix& matrix, GVector& vec, GVector& answer_vec)
{
    if(matrix.Column() != matrix.Row()){         // �W���s��͐����s��ł��邩
        return false;
    }else if( matrix.Row() != vec.Row()){        // �W���s��ƒ萔�������̃x�N�g���̍s���͈�v���邩
        return false;
    }else if( matrix.Row() != answer_vec.Row()){ // �W���s��Ɠ������i�[����x�N�g���̍s���͈�v���邩
        return false;
    }else if( vec.Shape() != GVector::COLUMN_VECTOR
       || answer_vec.Shape() != GVector::COLUMN_VECTOR ){ // �x�N�g���͗�x�N�g���ł��邩
        return false;
    }else{
        return true;
    }
}

---