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

今回は連立一次方程式に対する直接法である、LU分解のC++コードを公開します。行列を上三角行列Uと下三角行列Lに分解します。一度分解してしまえば、前進代入と後退代入という簡単な操作で連立一次方程式を解くことができるのです。理論は

http://www.ced.is.utsunomiya-u.ac.jp/lecture/2011/prog/p2/kadai3/no3/lu.pdf

がわかりやすく参考になります。肝は「どのように再帰的なアルゴリズムが構成されているか」を理解することでしょう。

しくみはおいおい説明するとして、とりあえずコードをおいておきます。私も今後自身のLU分解の解説を書きます。

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

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

であり、これを解くと となります。プログラムでは と はテキストファイルから読み込ませています。また、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 pivot;
	
	for(i=0;i<size-1;i++)
	{
		//l0 pivot//
		pivot=A[i*size+i];
		LU[i*size+i]=pivot;
		
		//l1//
		for(j=i+1;j<size;j++)
		{
			LU[j*size+i]=A[j*size+i];
		}
		
		//u1//
		for(j=i+1;j<size;j++)
		{
			LU[i*size+j]=A[i*size+j]/pivot;
		}
		
		//subtracting l1u1 from A or A_jk=lj*uk//
		for(j=i+1;j<size;j++)
		{
			for(k=i+1;k<size;k++)
			{
				A[j*size+k]-=LU[j*size+i]*LU[i*size+k];
			}
		}
	}
	
	//i=n//
	i=size-1;
	LU[i*size+i]=A[i*size+i];
	
}

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

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;
}