LU分解 pivot 付き Doolittle 法

配列の割り当てに Fortran2008 の source オプションを使ってみました。
pivotを使うと精度が上がります。

実行結果

ソース・プログラム

module m_mat
  implicit none
  integer, parameter :: kd = kind(0.0d0)
contains
  subroutine lu2(a, b, x)
    real(kd), intent(in) :: a(:, :)
    real(kd), intent(in) :: b(:)
    real(kd), intent(out) :: x(:)
    real(kd), allocatable :: c(:, :)
    real(kd) :: y(size(x))
    integer :: ipiv(size(b))
    allocate(c, source = a) ! f2008
    call lu(c, ipiv)
    y = fwd (c, b, ipiv)
    x = bkwd(c, y, ipiv)
    return
  end subroutine lu2 

  subroutine lu(a, ipiv) ! LU decomposition by Gaussian elimination without pivotting
    real(kd), intent(in out) :: a(:, :)
    integer, intent(out) :: ipiv(:)
    real(kd) :: alpha
    integer :: i, icol, irow, jcol, jrow, itmp
    ipiv = [ (i, i = 1, size(ipiv)) ]
    do icol = 1, size(a, 2)
      i = maxloc( abs( a(ipiv(icol:), icol) ), 1 ) + icol - 1
      itmp      = ipiv(icol)
      ipiv(icol) = ipiv(i)
      ipiv(i)    = itmp 
      jcol = ipiv(icol)
      do irow = icol + 1, size(a, 1)
       jrow = ipiv(irow)
       a(jrow, icol) = a(jrow, icol) / a(jcol, icol)                              ! L
       a(jrow, icol + 1:) = a(jrow, icol + 1:) - a(jrow, icol) * a(jcol, icol + 1:) ! U
     end do
    end do
    return
  end subroutine lu
  
  function fwd(a, b, ipiv) result(y) ! L y = b
    real(kd), intent(in) :: a(:, :), b(:)
    integer, intent(in) :: ipiv(:)
    real(kd) :: y(size(b))
    integer :: irow, jrow
    do irow = 1, size(b)
      jrow = ipiv(irow)
      y(irow) = b(jrow) - sum( a(jrow, :irow - 1) * y(:irow - 1) )
    end do
    return
  end function fwd

  function bkwd(a, y, ipiv) result(x) ! U x = y
    real(kd), intent(in) :: a(:, :), y(:)
    integer, intent(in) :: ipiv(:)
    real(kd) :: x(size(y))
    integer :: irow, jrow
    do irow = size(y), 1, -1
      jrow = ipiv(irow)
      x(irow) = ( y(irow) - sum( a(jrow, irow + 1:) * x(irow + 1:) ) ) / a(jrow, irow)
    end do
    return
  end function bkwd
  
end module m_mat

program Doolittle
  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 lu2(a, b, x)
    print *, sqrt( sum( abs(matmul(a, x) - b)**2 ) )
  end do
  stop
end program Doolittle