【数値計算】LU分解（Doolittle法）C++コード付き

今回は連立一次方程式に対する直接法である、LU分解（Doolittle法）のC++コードを公開します。

今回解く連立一次方程式は

とあらわされます。ここで、 は の行列、 と は 次元のベクトルです。 と が既知で、 が未知です。具体的には

であり、これを解くと となります。プログラムでは と はテキストファイルから読み込ませています。また、LU分解の結果は一つの配列（行列）LUに格納しています。ではコードを以下に示します。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <stdio.h>
#include <time.h>
#include <stdlib.h> 

using namespace std;

//matix vector multiplication Ax=b
inline void mvm(double A[],double x[],double b[], int size)
{
	int i,j;
	
	for(i=0;i<size;i++)
	{
		b[i]=0.;
		for(j=0;j<size;j++)
		{
			b[i]+=A[i*size+j]*x[j];
		}
	}
}

//inner product//
inline double ip(double a[], double b[], int size)
{
	int i;
	
	double value=0.;
	for(i=0;i<size;i++) value+=a[i]*b[i];
	return value;
}

//input//
inline void input(double A[], double b[], int size)
{
	int i,j;
	
	ifstream fo;
	fo.open("A.txt");
	for(i=0;i<size;i++)
	{
		for(j=0;j<size;j++)
		{
			fo>>A[i*size+j];
		}
	}
	fo.close();
	
	fo.open("b.txt");
	for(i=0;i<size;i++)
	{
		fo>>b[i];
	}
	fo.close();
}

//matrix output//
inline void matout(double A[], int size)
{
	int i,j;
	
	for(i=0;i<size;i++)
	{
		for(j=0;j<size;j++)
		{
			cout<<A[i*size+j]<<" ";
		}
		cout<<endl;
	}
	cout<<endl;
}

//vector output//
inline void vecout(double b[], int size)
{
	int i;
	
	for(i=0;i<size;i++)
	{
		cout<<b[i]<<endl;
	}
	cout<<endl;
}

//LU decomposition//
inline void LUd(double A[], double LU[], int size)
{
	int i,j,k;
	double sum;
	
	for(i=0;i<size;i++)
	{
		for(j=0;j<size;j++)
		{
			if(i>j)//l_ij
			{
				sum=0.0;
				for(k=0;k<j;k++) sum+=LU[i*size+k]*LU[k*size+j];
				LU[i*size+j]=(A[i*size+j]-sum)/LU[j*size+j];
			}
			else//u_ij
			{
				sum=0.0;
				for(k=0;k<i;k++) sum+=LU[i*size+k]*LU[k*size+j];
				LU[i*size+j]=A[i*size+j]-sum;
			}
		}
	}
}

//forward substitution and backward substitution//
void inline fb(double LU[], double x[], double y[], double b[], int size)
{
	int i,j;
	double sum;
	
	y[0]=b[0];
	
	//forward//
	for(i=1;i<size;i++)
	{
		sum=0.0;
		for(j=0;j<i;j++)
		{
			sum+=LU[i*size+j]*y[j];
		}
		y[i]=b[i]-sum;
	}
	
	//backward//
	x[size-1]=y[size-1]/LU[(size-1)*size+size-1];
	
	for(i=size-2;i>=0;i--)
	{
		sum=1.0*x[i];
		for(j=i+1;j<size;j++)
		{
			sum+=LU[i*size+j]*x[j];
		}
		x[i]=(y[i]-sum)/LU[i*size+i];
	}
}

main()
{
	int size=4;
	
	double *b=new double[size];
	double *x=new double[size];
	double *A=new double[size*size];
	double *LU=new double[size*size];
	double *y=new double [size];
	
	//input//
	input(A,b,size);
	
	//LU decomposition//
	LUd(A,LU,size);
	
	//forward substitution and backward substitution//
	fb(LU,x,y,b,size);
	
	//output//
	matout(LU,size);
	vecout(x,size);
	
	delete [] b,x,A,LU,y;

    return 0;
}

参考