【ガウスの消去法(Gaussian elimination method)】Fortranによるガウスの消去法アルゴリズム

はじめに

webメインでコーディングしている僕ですが、最近、数値計算のアルゴリズムを組む機会がありました。計算式をコーディングに落としこむ作業は、とりわけwebでのものと変わりはなかったですが計算アルゴリズムを理解するのが少し難しく感じました。
頭の整理のために、記事にまとめておきます。

ガウスの消去法とは

ガウスの消去法とは、連立一次方程式を解くための多項式時間アルゴリズムのことをいいます。掃き出し法(row reduction)とも呼ばれています。大きな方程式系を系統的な方法で小さな系へ分解するアルゴリズムです。
ガウスの消去法は連立一次方程式の解法としてよく利用されていますが、そのほかに、行列の階数の計算、行列式の計算、正則行列の逆行列の計算などにも使われます。

メリット

計算時間が短いことが挙げられます。全体としては、乗算が $\frac{n^3}{3}$ 回であり、大規模な連立一次方程式にも太刀打ちできます（対してガウスジョルダン消去法は $\frac{n^3}{2}$ 回）。

2つの特徴

• 前進消去(forward elimination)
  • 文字（数字）を一つずつ消去していく操作
  • 方程式の定数倍を他の方程式に加えることで文字を順番に消していく
• 後退代入(back substitution)
  • 一成分（変数）ずつ答えを求めていく操作
  • 後ろの式から順番にたどっていき、一つずつ変数の値を求める

参考：ガウスの消去法による連立一次方程式の解き方

解法

行列を上三角行列に変形することで、変数を求めます。
つまり、行に関する基本変形を行列に可能な限り繰り返し行って行列の左下部分の成分を全て 0 にする。

行に関する基本変形には、

• 二つの行を入れ替えるもの
• ある行を0でない定数倍するもの
• ある行に他のある行の定数倍を加えるもの

の3つがあります。

プログラム

メインプログラム(main program)

main.f90

program main
  implicit none
  integer N, M
  parameter (N = 3, M = 2)
  real a(N, N), b(N, M)
  integer i, j
  integer eflag

  open(1, file = 'a.txt')
  open(2, file = 'b.txt')
  do i = 1, N
    read (1, *) (a(i, j), j = 1, N)
    read (2, *) (b(i, j), j = 1, M)
  enddo
  close(1)
  close(2)

! 元の連立方程式の左辺
  write(*, *) 'A='
  do i = 1, N
    write(*, *) (a(i, j), j = 1, N)
  enddo
  write(*, *) 'B='
  do i = 1, N
    write(*, *) (b(i, j), j = 1, M)
  enddo

  eflag = 0

  call gauss_elim(a, N, b, M, eflag)

  if (eflag == 0) then
    write(*, *) 'X satisfying AX=B is'
    do i = 1, N
      write(*, *) (b(i, j), j = 1, M)
    enddo
  else
    write(*, *) 'ERROR', eflag
  endif

  stop
end

サブルーチン(subroutine)

sub.f90

subroutine gauss_elim(a, n, b, m, eflag)
  implicit none
  integer n, m
  real a(n ,n), b(n, m)
  integer eflag

  real TINY
  parameter(TINY=1.0e-20)

  integer i, j, k
  real r

! Gaussian elimination
  do j = 1, n-1
    call ppivot(a, n, b, m, j)
    if (ABS(a(j, j)) < TINY) then
      eflag = 1
      return
    endif

    ! forward elimination
    do i = j+1, n
      r = a(i, j) / a(j, j)
      ! k = j, nでも良いが計算を少なくできるならそれに越したことはない（使わない値k=j）
      do k = j+1, n
        a(i, k) = a(i, k) - a(j, k) * r
      enddo
      do k = 1, m
        b(i, k) = b(i, k) - b(j, k) * r
      enddo
    enddo
  enddo

  if (ABS(a(n, n)) < TINY) then
    eflag = 1
    return
  endif

! Back substitution
  do j = 1, m
    do i = n, 1, -1
      r = b(i, j)
      do k = i+1, n
        r = r - a(i, k) * b(k, j)
      enddo
      b(i, j) = r / a(i, i)
    enddo
  enddo
  return
end subroutine gauss_elim

subroutine ppivot(a, n, b, m, j)
  implicit none
  integer n, m, j
  real a(n, n), b(n, m)

  integer i, k, imax
  real amax, r
! find the maximum element
  amax = ABS(a(j, j))
  imax = j
  do i = j+1, n
    r = ABS(a(i, j))
    if (r > amax) then
      amax = r
      imax = i
    endif
  enddo
! exchange rows if necessary
if (j == imax) return
do k = 1, n
  r = a(imax, k)
  a(imax, k) = a(j, k)
  a(j, k) = r
enddo
do k = 1, m
  r = b(imax, k)
  b(imax, k) = b(j, k)
  b(j, k) = r
enddo
return
end subroutine ppivot

実行

コンパイル

メインルーチンプログラムとサブルーチンプログラムを対象にしてコンパイル

$ gfortran main.f90  sub.f90

確認

a.out実行ファイルが作成されているか確認

$ ls
a.txt			b.txt			a.out*			sub.f90	        main.f90

実行

$ ./a.out

おわりに

Fortranのプログラム全体の流れ・文法の基礎に関してはこちらにまとめていますので必要に応じてご参考ください。