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

作業中のメモ

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

ガウスの消去法(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