読者です 読者をやめる 読者になる 読者になる

作業中のメモ

よく「計算機」を使って作業をする.知らなかったことを中心にまとめるつもり.

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

とすればよい.