ガウス消去法 枢軸変換なし

線型方程式 Ax=b をガウス消去法で解くプログラム。Fortran90 だと式をそのままプログラムに直せばいいので楽。でもPivotting を入れるとやや見苦しくなる。

並列計算機では LU分解で Crout法より Gauss 消去法の方が有利になることもあるとか。ベクトル機全盛時代や RISC Cache 重要時代が遠くへ行くにつれ計算アルゴリズムの流行も変わるからをかし。

実行結果　倍精度

今の計算機は速いから素朴なやり方でも 500*500 くらいなら楽勝。計算量は O(N^3) だから行列サイズ倍にすると目安的に約８倍くらい時間がかかることになります。

ソース・プログラム

実質１０行くらいで書けてしまえるのがうれしい Fortran90 ！

module m_mat
  implicit none
  integer, parameter :: kd = kind(0.0d0)
contains
  subroutine gaussbk(a, b, x) ! gaussian elimination without pivotting
    real(kd), intent(in out) :: a(:, :), b(:)
    real(kd), intent(   out) :: x(:)
    real(kd) :: alpha
    integer :: i, icol, irow
    do icol = 1, size(a, 2) ! forward elimination
     do irow = icol + 1, size(a, 1)
       alpha = a(irow, icol) / a(icol, icol)
       a(irow, :) = a(irow, :) - alpha * a(icol, :)
       b(irow)    = b(irow)   - alpha * b(icol)
     end do
    end do
    do irow = size(x), 1, -1 ! backward substitution
      x(irow) = ( b(irow) - sum( a(irow, irow + 1:) * x(irow + 1:) ) ) / a(irow, irow)
    end do
    return
  end subroutine gaussbk
end module m_mat

program Gauss
  use m_mat
  implicit none
  integer, parameter :: n = 200
  real(kd) :: a(n, n), b(n), x(n)
  real(kd) :: c(n, n), d(n)
  call random_seed()
  call random_number(a)
  call random_number(b)
  c = a
  d = b
  call gaussbk(a, b, x)
  print *, sqrt( sum( (matmul(c, x) - d)**2 ) ) / n  
  stop
end program Gauss