LU分解 Crout法 部分Pivot付き

Pivot選択のためにちょっと余分な計算をしなければなりません。

実行結果

ソース・プログラム

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(:, :), y(:)
    integer :: ipiv(size(b))
    c = a ! F2003 /assume:realloc_lhs
    call lu(c, ipiv)
    y = fwd (c, b, ipiv)
    x = bkwd(c, y, ipiv)
    return
  end subroutine lu2
  
  subroutine lu(a, ipiv) ! Crout's method
    real(kd), intent(in out) :: a(:, :)
    integer, intent(out) :: ipiv(:)   
    real(kd) :: beta, wk(size(a, 1))
    integer :: i, icol, irow, jrow, jcol, itmp
    ipiv = [(i, i = 1, size(ipiv))]
    do icol = 1, size(a, 2)
      do irow = 1, icol - 1 
        jrow = ipiv(irow)
        a(jrow, icol) = a(jrow, icol) - sum( a(jrow, :irow - 1) * a(ipiv(:irow - 1), icol) )
      end do
  
      do irow = icol, size(a, 1) !pivotting axis->max|a(jrow, icol) - sum( a(jrow, :irow - 1) * a(ipiv(:irow - 1), icol) )|
        jrow = ipiv(irow) ! irow = icol
        wk(jrow) = a(jrow, icol) - sum( a(jrow, :irow - 1) * a(ipiv(:irow - 1), icol) )
      end do
      i = maxloc( abs( wk(ipiv(icol:)) ), 1 ) + icol - 1
      itmp      = ipiv(icol)
      ipiv(icol) = ipiv(i)
      ipiv(i)    = itmp 
      
      jrow = ipiv(icol) ! irow = icol
      a(jrow, icol) = a(jrow, icol) - sum( a(jrow, :icol - 1) * a(ipiv(:icol - 1), icol) )
    
      beta = 1.0_kd / a(jrow, icol) 
      do irow = icol + 1, size(a, 1)
        jrow = ipiv(irow)
        a(jrow, icol) = beta * ( a(jrow, icol) - sum( a(jrow, :icol - 1) * a(ipiv(:icol - 1), icol) ) )
      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 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