LU 分解で連立一次方程式を解く（Fortran バージョン）
どうも，筆者です．前回に続き，Fortran で LU 分解を行うプログラムを作成してみた．LU 分解については，以前の記事を参照すること．
workspacememory.hatenablog.com
LU 分解のプログラム
需要あるのだろうか．これは，「luDecompSolver.f」というファイル名で保存している．
!----------------------------------------------------------------------! !---------------------- LU Decomposition ----------------------! !----------------------------------------------------------------------! !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12 ! 1 2 3 4 5 6 7 subroutine luDecomp(n,pivot,matrix) implicit none ! === arguments === integer*4 :: n integer*4 :: pivot(n) double precision :: matrix(n,n) ! n: rank, scalar ! pivot: pivot vector, pivot(row) ! matrix: coefficient matrix, matrix(row,col) ! --- local variables --- integer*4 :: row,col,k,ip,iptmp double precision :: tmp,valmax,weight,mineps mineps=1e-10 do k=1,n pivot(k)=k end do do k=1,n-1 valmax=abs(matrix(k,k)) ip=k do row=k+1,n tmp=abs(matrix(row,k)) if (valmax.lt.tmp) then valmax=tmp ip=row end if end do ! the diagonal component is less than mineps if (valmax.lt.mineps) then write(*,100) 'diagonal component is less than',mineps 100 format(a,1X,G15.7) return else if (ip.ne.k) then do col=k,n tmp=matrix(ip,col) matrix(ip,col)=matrix(k,col) matrix(k,col)=tmp end do iptmp=pivot(ip) pivot(ip)=pivot(k) pivot(k)=iptmp do col=1,k-1 tmp=matrix(k,col) matrix(k,col)=matrix(ip,col) matrix(ip,col)=tmp end do end if do row=k+1,n weight=matrix(row,k)/matrix(k,k) matrix(row,k)=weight do col=k+1,n matrix(row,col) * =matrix(row,col)-weight*matrix(k,col) end do end do end do return end !----------------- last line of LU Decomposition ------------------! ! 1 2 3 4 5 6 7 !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12 !----------------------------------------------------------------------! !---------------------- LU Solver ----------------------! !----------------------------------------------------------------------! !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12 ! 1 2 3 4 5 6 7 subroutine luSolver(n,pivot,matrix,b) implicit none ! === arguments === integer*4 :: n integer*4 :: pivot(n) double precision :: matrix(n,n),b(n) ! n: rank, scalar ! pivot: pivot vector, pivot(row) ! matrix: coefficient matrix, matrix(row,col) ! b: right side vector, b(row) ! --- local variables --- integer*4 :: row,id double precision :: tmp,y(n) ! ===================== ! forkward substitution ! ===================== y(1)=b(pivot(1)) do row=2,n id=row-1 y(row)=b(pivot(row))-sum(matrix(row,1:id)*y(1:id)) end do ! ===================== ! backward substitution ! ===================== b(n)=y(n)/matrix(n,n) do row=n-1,1,-1 id=row+1 tmp=y(row)-sum(matrix(row,id:n)*b(id:n)) b(row)=tmp/matrix(row,row) end do end !----------------- last line of LU Solver ------------------! ! 1 2 3 4 5 6 7 !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
main routine
いつも通り，main routine のコードも載せておく．ちょっと調べて，動的配列を使ってみた．この動的配列は Fortran90 から導入されたらしい．エラー処理等を入れたので，goto 文とかがあって読みにくい．「main_lu.f」として保存した．
!----------------------------------------------------------------------! !--------------------------- main code ---------------------------! !----------------------------------------------------------------------! !23456789|123456789|123456789|123456789|123456789|123456789|123456789|12 ! 1 2 3 4 5 6 7 implicit none character :: fmtstr*32,argv*16,fname*16 integer*4 :: n,i,j integer*4,allocatable :: pivot(:) double precision,allocatable :: matrix(:,:),b(:) ! === check arguments === if (iargc().ne.1) then call getarg(0,argv) write(*,'("Usage: ", a, "[data file name]")') argv stop ! === read arguments === else call getarg(1,argv) read(argv,'(a16)') fname end if ! === read data from file === open(10,file=fname,status='old',err=1) goto 100 1 write(*,'("Error: ", a, "cannot open")') fname stop 100 continue read(10,*)n allocate(pivot(n),matrix(n,n),b(n)) do i=1,n read(10,*,end=200) (matrix(i,j),j=1,n),b(i) end do 200 continue close(unit=10) ! =========================== write(fmtstr,'("("I0"(G15.7,1X),G15.7)")') n do i=1,n write(*,fmtstr) (matrix(i,j),j=1,n),b(i) end do ! execute LU decomposition call luDecomp(n,pivot,matrix) ! solve call luSolver(n,pivot,matrix,b) write(*,'(/ "answer:")') write(*,'(G15.7)') (b(i),i=1,n) deallocate(pivot,matrix,b) stop end
読み込むファイルは，
[次元数（n: rank）] [係数行列 (1,1) 成分] [係数行列 (1,2) 成分] ... [係数行列 (1,n) 成分] [右辺ベクトル (1)] [係数行列 (2,1) 成分] [係数行列 (2,2) 成分] ... [係数行列 (2,n) 成分] [右辺ベクトル (2)] ... ... ... ... ... ... ... ... [係数行列 (n,1) 成分] [係数行列 (n,2) 成分] ... [係数行列 (n,n) 成分] [右辺ベクトル (n)]
という形式にする．前回の例だと，
3 2.0 1.0 -2.0 1.0 1.0 1.0 -1.0 4.0 1.0 -2.0 3.0 -1.0
となる．
Makefile
前回と同様に，Makefile も載せておく．
FC := gfortran PROG := solveLU SRCS := luDecompSolver.f main_lu.f OBJS := $(SRCS:%.f=%.o) all: $(PROG) $(PROG): $(OBJS) $(FC) -o $@ $^ %.o: %.f $(FC) -c -O2 $< .PHONY: clean clean: rm -f $(OBJS) $(PROG)
実行方法
Fortran のプログラムを上から順に，「luDecompSolver.f」，「main_lu.f」として保存したとする．また，Makefile を「Makefile」として保存したとする．その段階で，
[hoge@local]$ ls Makefile luDecompSolver.f main_lu.f
となっているはずである．まず，make コマンドを利用してコンパイルする．
[hoge@local]$ make gfortran -c -O2 luDecompSolver.f gfortran -c -O2 main_lu.f gfortran -o solveLU luDecompSolver.o main_lu.o
「solveLU」という実行ファイルが出来上がる．もし，make コマンドがなければ，上の実行結果（gfortran から始まる部分）をコピペして端末で実行すればよい．後は，データファイル input.dat を用意して，
[hoge@local]$ ./solveLU input.dat
とすれば，計算結果が出力される．また，実行ファイル，オブジェクトファイルの削除は，
[hoge@local]$ make clean rm -f luDecompSolver.o main_lu.o solveLU
とすればよい．
