ガウス消去法　枢軸交換付き

もともとの行列が上書きされてしまうので、中間サブルーチンを置いて作業行列などを確保します。枢軸（pivot)用のベクトルもここで用意します。

実行結果

500*500の行列の線型方程式を20回解いて残差の二乗平均を出力。おおむね数値誤差程度のまで求められています。

ソース・プログラム

f90系で書くとプログラムもスッキリして行数も減って明快になるのでうれし楽しい。

module m_mat
  implicit none
  integer, parameter :: kd = kind(0.0d0)
contains
  subroutine gaussbk2(a, b, x)
    real(kd), intent(in)  :: a(:, :), b(:)
    real(kd), intent(out) :: x(:)
    real(kd), allocatable :: c(:, :), d(:)
    c = a ! F2003 /assume:realloc_lhs or /standard-semantics
    d = b !
    call gaussbk(c, d, x)
    return
  end subroutine gaussbk2
  
  subroutine gaussbk(a, b, x) ! gaussian elimination with partial pivotting
    real(kd), intent(in out) :: a(:, :), b(:)
    real(kd), intent(   out) :: x(:)
    real(kd) :: alpha
    integer :: i, icol, irow, jcol, jrow, itmp
    integer :: ipiv(size(b))
    ipiv = [ (i, i = 1, size(ipiv)) ]
    do icol = 1, size(a, 2) ! forward elimination
      i = maxloc( abs(a(ipiv(icol:), icol)), 1 ) + icol - 1 ! choose pivot
      itmp      = ipiv(icol)
      ipiv(icol) = ipiv(i)
      ipiv(i)    = itmp 
      jcol = ipiv(icol)
      do irow = icol + 1, size(a, 1)
       jrow = ipiv(irow)
       alpha = a(jrow, icol) / a(jcol, icol)
       a(jrow, :) = a(jrow, :) - alpha * a(jcol, :)
       b(jrow)    = b(jrow)    - alpha * b(jcol)
     end do
    end do
    do i = size(x), 1, -1 ! backward substitution
      irow = ipiv(i)
      x(i) = ( b(irow) - sum( a(irow, i + 1:) * x(i + 1:) ) ) / a(irow, i)
    end do
    return
  end subroutine gaussbk
end module m_mat

program Gauss
  use m_mat
  implicit none
  integer, parameter :: n = 500
  real(kd) :: a(n, n), b(n), x(n)
  integer :: i
  call random_seed()
  do i = 1, 20
    call random_number(a)
    call random_number(b)
    call gaussbk2(a, b, x)
    print *, sqrt( sum( (matmul(a, x) - b)**2 ) ) / n  
  end do
  stop
end program Gauss

      i = maxloc( abs(a(ipiv, icol)), 1, [spread(.false., 1, icol - 1), spread(.true., 1, size(a, 2) - icol)] ) 

      i = maxloc( abs(a(ipiv, icol)), 1, mask = [ (i >= icol, i = 1, size(a, 1)) ] )