LU分解 Doolittle 法

基本的にガウス消去法と同じ方式になっています。

実行結果

LU分解では精度がやや悪くなるようです。

Eliza’s father=Alfred P.DoolittleStation master=AlfeeYoung Freddy=Sherlock Holmes

ソース・プログラム

簡単に書けるのがうれしい。ほとんど実質３〜４行ｗ

module m_mat
  implicit none
  integer, parameter :: kd = kind(0.0d0)
contains
  subroutine lu(a) ! LU decomposition by Gaussian elimination without pivotting
    real(kd), intent(in out) :: a(:, :)
    integer :: icol, irow
    do icol = 1, size(a, 2) ! forward elimination
     do irow = icol + 1, size(a, 1) 
       a(irow, icol) = a(irow, icol) / a(icol, icol)
       a(irow, icol + 1:) = a(irow, icol + 1:) - a(irow, icol) * a(icol, icol + 1:)
     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 :: irow
    do irow = 1, size(b)
      y(irow) = b(irow) - sum( a(irow, :irow - 1) * y(:irow - 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 :: irow
    do irow = size(y), 1, -1
      x(irow) = ( y(irow) - sum( a(irow, irow + 1:) * x(irow + 1:) ) ) / a(irow, irow)
    end do
    return
  end function bkwd
end module m_mat

program Gauss
  use m_mat
  implicit none
  integer, parameter :: n = 500
  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 lu(a)
  x = bkwd(a, fwd(a, b))
  print *, sqrt( sum( abs(matmul(c, x) - d)**2 ) )
  stop
end program Gauss