作業中のメモ

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

トップ > 数学 > ガウスの消去法（Fortran バージョン）

ガウスの消去法（Fortran バージョン）

数学 プログラミング

お久しぶりです．筆者です．

アクセス解析のデータをみた

久しぶりに，ブログに戻ってきた．最初に，アクセス解析のデータをみた． すると，多くの人がガウスの消去法関連でアクセスしている事が分かった．その中でも，Fortran 言語でガウスの消去法を探している人が多い事に気付いた．

ガウスの消去法を Fortran 言語に書き直す

筆者は，ほとんど C 言語を利用するが，Fortran も少しは触れた事があるので，色々調べながら，以前書いた記事「ガウスの消去法のプログラム」を Fortran77 の固定形式で書き直す事にした． 何故，Fortran77 の固定形式かというと，筆者がそれしか読んだ事がないためである．この平成の時代に，Fortran77 の固定形式を勉強するとは思わなかったが，昔の方が書いたプログラムを使う以上，仕方ない事である．

workspacememory.hatenablog.com

書き直したもの

という事で，色々調べながら，C 言語のプログラムを Fortran で書き直した．読んだ事あるのは固定形式だが，書いた事はないので，コンパイルが通る程度にしか書けない．90/95 の表現もあるかもしれない．

!----------------------------------------------------------------------!
!----------------------   Gaussian elimination   ----------------------!
!----------------------------------------------------------------------!
!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12

      subroutine gaussSolver(n,matrix,b)
      implicit none

!     === arguments ===
      integer*4 :: n
      double precision :: matrix(n,n),b(n)
!     n:      rank, scalar
!     matrix: coefficient matrix, matrix(row,col)
!     b:      right side vector, b(row)

!     --- local variables ---
      integer*4 :: i,row,col,pivot
      double precision :: tmp,valmax,diag,scl,mineps
      mineps=1e-10

      do i=1,n-1
!       execute partial pivoting
        pivot=i
        valmax=abs(matrix(i,i))

!       search maximum value in a column
        do row=i+1,n
          tmp=abs(matrix(row,i))

          if (valmax.lt.tmp) then
            pivot=row
            valmax=tmp
          end if
        end do

!       replace
        if (pivot.ne.i) then
          do col=1,n
            tmp=matrix(i,col)
            matrix(i,col)=matrix(pivot,col)
            matrix(pivot,col)=tmp
          end do
          tmp=b(i)
          b(i)=b(pivot)
          b(pivot)=tmp
        end if
        diag=matrix(i,i)

        if (abs(diag).lt.mineps) then
          return
        end if

!       ===================
!       forward elimination
!       ===================
        do row=i+1,n
          scl=matrix(row,i)/diag

          do col=i,n
            matrix(row,col)=matrix(row,col)-matrix(i,col)*scl
          end do
          b(row)=b(row)-b(i)*scl
        end do
      end do

!     =====================
!     backward substitution
!     =====================
      b(n)=b(n)/matrix(n,n)
      do row=n-1,1,-1
        tmp=b(row)-sum(matrix(row,(row+1):n)*b((row+1):n))
        b(row)=tmp/matrix(row,row)
!------------------------------------------------------------
!     do row=n,1,-1
!        valsum=b(row)
!
!        do col=row+1,n
!          valsum=valsum-matrix(row,col)*b(col)
!        end do
!        b(row)=valsum/matrix(row,row)
!      end do
!------------------------------------------------------------

      return
      end

main routine

ついでに，main routine の部分も以下に示しておく．

!----------------------------------------------------------------------!
!---------------------------   main  code   ---------------------------!
!----------------------------------------------------------------------!
!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12

      implicit none

      character :: fmtstr*32
      integer*4 :: n,i,j
      parameter(n=3)
      double precision :: matrix(n,n),b(n)

      matrix(1,1)=2.D0
      matrix(1,2)=1.D0
      matrix(1,3)=-2.D0
      matrix(2,1)=1.D0
      matrix(2,2)=1.D0
      matrix(2,3)=-1.D0
      matrix(3,1)=1.D0
      matrix(3,2)=-2.D0
      matrix(3,3)=3.D0
      b(1)=1.D0
      b(2)=4.D0
      b(3)=-1.D0
      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

      ! output matrix and right side vector
      open(10,file='input.txt',status='replace')
      do i=1,n
        write(10,fmtstr) (matrix(i,j),j=1,n),b(i)
        matrix(i,:)=0.D0
        b(i)=0.D0
      end do
      close(unit=10)

      ! read data from file
      open(11,file='input.txt',status='old')
      do i=1,n
        read(11,*,end=100) (matrix(i,j),j=1,n),b(i)
      end do
100   continue
      close(unit=11)

      ! execute Gaussian elimination
      call gaussSolver(n,matrix,b)

      write(*,'(/ "answer:")')
      write(*,'(G15.7)') (b(i),i=1,n)

      stop

      end

いつも通り，英語が出来ないのはスルーして欲しい．無駄にファイル出力して，ファイル読み込みしているが，これは，単純に，Fortran で読み書きが出来るかどうかを確認しているだけである．たいした意味は無い．

Makefile

上のプログラムは自由に使ってもらってかまわないが，コンパイルに困る可能性もある．そのため，Makefile を作成した．Ubuntu 14.04 LTS でインストールできる gfortran で動作を確認している． Makefile を使った事がない人は，1 つのファイルに全て書いて，それをコンパイルしてもらえればよいと思う． また，この gfortran は，Fortran95 用のコンパイラなので，フリーダムな書き方をしていても文法のミスさえなければ通るみたいだ．以下に，Makefile を示す．

FC := gfortran
PROG := solveGauss
SRCS := gaussSolver.f prog_main.f
OBJS := $(SRCS:%.f=%.o)

all: $(PROG)

$(PROG): $(OBJS)
  $(FC) -o $@ $^

%.o: %.f
  $(FC) -c -O2 $<

.PHONY: clean
clean:
  rm -f $(OBJS) $(PROG)

sub routine 「gaussSolver」を「gaussSolver.f」として，main routine 「main code」を「prog_main.f」として保存し，make とすると，実行ファイル「solveGauss」が出来上がる．後は，実行すると解が得られる． make clean とすると，.o ファイルと solveGauss が削除される．

[hoge@local]$ ls
Makefile gaussSolver.f prog_main.f

[hoge@local]$ make
gfortran -c -O2 gaussSolver.f
gfortran -c -O2 prog_main.f
gfortran -o solveGauss gaussSolver.o prog_main.o

[hoge@local]$ make clean
rm -f gaussSolver.o prog_main.o solveGauss