使用環境など

この記事ではWindows10にインストールしたmingw32-gfortranを使っています。
言語は、fortran95を使っています。

LU分解とは

LU分というのは、行列Aを2つの行列LとUに分解することを言います。
これを用いることで、連立一次方程式の計算や行列式計算ができるようになります。
LU分解を使って連立一次方程式Ax=bを解くことのメリットとしては、一度AをLU分解をしたらそのAについての解をbの値を変えるだけで、解を求めることができます。

LUを用いて連立一次方程式を解く方法

先ず、
$$ A\vec{x} = \vec{b} $$
をLU分解します。つまり、以下の式が求まります。
$$ LU\vec{x} = \vec{b} $$
この式を以下のように考えます。
$$ L\vec{y} = \vec{b} $$
$$ U\vec{x} = \vec{y} $$
下の式を上の式に代入することで、解を簡単に求めることができます。

ソースコード

LU_decomposition.f95

program main
	implicit none
	real(8), allocatable, dimension(:,:) :: a,L,U
	real(8), allocatable, dimension(:) :: b,x,y
	real(8) temp, alpha, reserve, sum1, sum2
	integer i, j, k, m, n, dim, maxline

	print *, 'input the number of dimension'
	read(*,*) dim
	allocate (a(dim, dim), b(dim), L(dim,dim), U(dim,dim),x(dim),y(dim))

!	Initialize ------------------------------------------------------
	a = 0.0
	b = 0.0
	L = 0.0
	U = 0.0

!	Input -----------------------------------------------------------
!	A
	do i = 1, dim
		do j = 1, dim
			print *,"x(",i,',',j,')'
			read *,a(i,j)
		end do
	end do
!	b
	do i = 1, dim
		print *,"b(",i,')'
		read *,b(i)
	end do

!	print the input -----------------------------------------------
	print *, 'A'
	do i = 1, dim
		print *, a(i, 1:dim)
	end do

	print *, 'B'
	print *, b(1:dim)

!	Solve the problem------------------------------------------------
!	LU --------------------------------------------------------------

	do k = 1, dim
		L(k,k) = 1.0
	end do

	do i = 1, dim
		do j = 1, dim
			if ( i .le. j ) then
				sum1 = 0.0
				do m = 1, i -1
					sum1 = sum1 + L(i,m) * U(m,j)
				end do
				U(i,j) = a(i,j) - sum1
			end if
			if (i .gt. j) then
				sum2 = 0.0
				do n = 1, j-1
					sum2 = sum2 + L(i,n) * U(n,j)
				end do
				L(i,j) = (a(i,j) - sum2) / U(j,j)
			end if
		end do
	end do

!	----------------------------------------------------------------

!	print for check-------------------------------------------------
	print *, 'L'
	do i = 1,dim
		print *, L(i,1:dim)
	end do

	print *, 'U'
	do i = 1, dim
		print *, U(i,1:dim)
	end do

!   solve the problem-----------------------------------------------
!   Ly=b------------------------------------------------------------
    do i = 1, dim
        y(i) = b(i)
    end do

	do j = 1, dim - 1
        do i = j + 1, dim
          y(i) = y(i) - y(j) * L(i,j)
        end do
    end do

!   Ux=y------------------------------------------------------------
    do i = 1, dim
        x(i) = y(i)
    end do

	do j = dim, 1, -1
        x(j) = x(j) / U(j,j)
        do i = 1, j - 1
          x(i) = x(i) - U(i,j) * x(j)
        end do
      end do
!   ----------------------------------------------------------------

!   print-----------------------------------------------------------
    print *, 'x'
    do i = 1, dim
		print *, x(i)
	end do
!   ----------------------------------------------------------------

	deallocate (a,b,L,U,x,y)
end

実行結果

LUdecomposition.exe

 input the number of dimension
2
 x(           1 ,           1 )
1
 x(           1 ,           2 )
2
 x(           2 ,           1 )
3
 x(           2 ,           2 )
1
 b(           1 )
2
 b(           2 )
1
 A
   1.0000000000000000        2.0000000000000000
   3.0000000000000000        1.0000000000000000
 B
   2.0000000000000000        1.0000000000000000
 L
   1.0000000000000000        0.0000000000000000
   3.0000000000000000        1.0000000000000000
 U
   1.0000000000000000        2.0000000000000000
   0.0000000000000000       -5.0000000000000000
 x
   0.0000000000000000
   1.0000000000000000

参考文献：
「物理のかぎしっぽ」http://hooktail.org/computer/index.php?LU%CA%AC%B2%F2
「numerical algorithms group」http://www.nag-j.co.jp/fortran/