LU分解 Crout法

実行結果

ソース・プログラム

module m_mat
  implicit none
  integer, parameter :: kd = kind(0.0d0)
contains
  subroutine lu2(a, b, x)
    real(kd), intent(in)  :: a(:, :), b(:)
    real(kd), intent(out) :: x(:)
    real(kd), allocatable :: c(:, :)
    c = a ! F2003 /assume:realloc_lhs
    call lu(c)
    x = bkwd(c, fwd(c, b))
    return
  end subroutine lu2
  
  subroutine lu(a) ! Crout's method
    real(kd), intent(in out) :: a(:, :)
    real(kd) :: beta
    integer :: i, icol, irow
    do icol = 1, size(a, 2)
      do irow = 1, icol 
        a(irow, icol) = a(irow, icol) - sum( a(irow, :irow - 1) * a(:irow - 1, icol) )
      end do
      beta = 1.0_kd / a(icol, icol) 
      do irow = icol + 1, size(a, 1)
        a(irow, icol) = beta * ( a(irow, icol) - sum( a(irow, :icol - 1) * a(:icol - 1, icol) ) )
      end do
    end do
    return
  end subroutine lu
  
  function fwd(a, b) result(y) ! L y = b
    real(kd), intent(in) :: a(:, :), b(:)
    real(kd) :: y(size(b))
    integer :: i
    do i = 1, size(b)
      y(i) = b(i) - sum( a(i, :i - 1) * y(:i - 1) )
    end do
    return
  end function fwd

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

program ludecomp
  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( (matmul(a, x) - b)**2 ) ) / n  
  end do
  stop
end program ludecomp