ガウスの消去法（Fortran 自由形式）

Fortran90 の自由形式に関して

どうも，筆者です．何度同じ記事書くんだ？って感じだが，どうも自分が扱っていた Fortran の記述方法は古いらしい． 無料のコンパイラがあることを含めると，Fortran 90/95 の自由形式（free format）が主流らしい．Fortran 2003/2008 とかもあるらしい．

という事で，Fortran 90/95 の自由形式で記述する練習を兼ねてプログラムを書き直した．また，配列の使い方も何となく理解してきた．もっとも，同じプログラムを違う形式で書き換えているだけなので，動作は同じである．

プログラム

まず，ガウスの消去法を行う subroutine を示す．モジュールなるものを使ってみたりした．

! ====================
! Gaussian elimination
! ====================
!23456789|123456789|123456789|123456789|123456789|123456789|123456789|123456789|123456789|123456789
!        1         2         3         4         5         6         7         8         9
module mdlGaussElim
  implicit none
  private
  public :: GaussianElimination
  real(8), parameter :: eps_default = real(1e-10, 8)

  contains

  ! =================================
  ! ガウスの消去法
  !
  ! A: n × n の係数行列，A(row, col)
  ! b: n 次元右辺ベクトル，b(row)
  ! =================================
  subroutine GaussianElimination(A, b, mineps)
    ! 引数
    real(8), intent(inout) :: A(:, :), b(:)
    real(8), intent(in), optional :: mineps

    ! 局所変数
    integer(4) :: n, i, id, pivot
    real(8) :: val, diag, weight, use_eps
    n = size(A, 1)

    ! 許容誤差の確認
    use_eps = eps_default
    if (present(mineps)) then
      use_eps = mineps
    end if

    if (n /= size(A, 2)) then
      write(*, '("A is not square matrix")')
      return
    end if

    ! ================
    ! 部分ピボット選択
    ! ================
    do i = 1, n - 1
      id = maxloc(abs(A(i:n, i)), 1)

      ! 対角要素が最大でない場合
      if (id /= 1) then
        pivot = id + i - 1

        ! 入れ替えを行う
        do id = 1, n
          val = A(i, id)
          A(i, id) = A(pivot, id)
          A(pivot, id) = val
        end do
        val = b(i)
        b(i) = b(pivot)
        b(pivot) = val
      end if

      ! ========
      ! 前進消去
      ! ========
      diag = A(i, i)

      if (abs(diag) < use_eps) then
        write(*, '("Error: diagonal component is less than", G15.7)') use_eps
        return
      end if

      do id = i + 1, n
        weight = A(id, i) / diag
        A(id, i:n) = A(id, i:n) - weight * A(i, i:n)
        b(id) = b(id) - weight * b(i)
      end do
    end do

    ! ========
    ! 後退代入
    ! ========
    b(n) = b(n) / A(n, n)
    do i = n - 1, 1, -1
      id = i + 1
      val = b(i) - sum(A(i, id:n) * b(id:n))
      b(i) = val / A(i, i)
    end do

    return
  end subroutine GaussianElimination
end module mdlGaussElim

次に，main routine を示す．

! ============
! main routine
! ============
!23456789|123456789|123456789|123456789|123456789|123456789|123456789|123456789|123456789|123456789
!        1         2         3         4         5         6         7         8         9

program main
  use mdlGaussElim
  implicit none

  character(32) :: fmtstr, fname
  integer(4) :: n, i, ioerr
  real(8), allocatable :: A(:, :), b(:)

  ! ======================
  ! === 引数のチェック ===
  ! ======================
  if (iargc() /= 1) then
    call getarg(0, fname)
    write(*, '("Usage: ", a, "[data file name]")') fname
    stop
  else
    ! ======================
    ! === 引数の読み込み ===
    ! ======================
    call getarg(1, fname)
  end if

  ! ========================
  ! === ファイルオープン ===
  ! ========================
  open(10, file = fname, status = 'old', iostat = ioerr)

  if (ioerr /= 0) then
    write(*, '("Error: ", a, " cannot open")') fname
    stop
  end if

  ! ====================================
  ! === ファイルからデータを読み込む ===
  ! ====================================
  read(10, *) n

  allocate(A(1:n, 1:n), b(1:n))
  write(fmtstr, '("(", i0, "(E20.13), 3X, E20.13)")') n

  do i = 1, n
    read(10, *, iostat = ioerr) A(i, 1:n), b(i)

    if (ioerr /= 0) exit
  end do
  close(unit = 10)

  ! ============================
  ! === ガウスの消去法を行う ===
  ! ============================
  call GaussianElimination(A, b)

  ! ============
  ! === 出力 ===
  ! ============
  write(*, '(E20.13)') (b(i), i = 1, n)

  deallocate(A, b)

  stop
end program main

そして，最後に Makefile である．

FC := gfortran
FFLAG := -O2
PROG := gaussSolver
SRCS := gaussianElimination.f90 mainRoutine.f90

all: $(PROG)

$(PROG): $(SRCS)
  $(FC) -o $@ $(FFLAG) $^

.PHONY: clean
clean:
  rm -f $(PROG) *.mod

ファイル読み込みは，以前と同様で，最初に次元数，次に係数行列，右辺ベクトルという並びのものを用意する．