Fortran 95 source code to do LU-decomposition by Crout method.

lu_decomposition_3.f95

!************************************************************
! LU 分解（クラウト法(Crout method)）
!
!   date          name            version
!   2019.03.08    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2019 mk-mode.com All Rights Reserved.
!************************************************************
!
module const
  ! SP: 単精度(4), DP: 倍精度(8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module const

module lu
  use const
  implicit none
  private
  public :: decompose

contains
  ! LU 分解
  ! * クラウト法(Crout method)
  !
  ! :param(inout) real(8) a(:,:): 行列
  subroutine decompose(a)
    implicit none
    real(DP), intent(inout) :: a(:, :)
    integer(SP) :: i, j, k, n
    real(DP)    :: s, tmp

    n = int(sqrt(real(size(a))))
    do k = 1, n
      do i = k, n
        s = sum((/(a(i, j) * a(j, k), j=1,k-1)/))
        a(i, k) = a(i, k) - s
      end do
      if (a(1, 1) == 0.0_DP) then
        print *, "Can't divide by 0 ..."
        stop
      end if
      tmp = 1.0_DP / a(k, k)
      do j = k + 1, n
        s = sum((/(a(k, i) * a(i, j), i=1,k-1)/))
        a(k, j) = (a(k, j) - s) * tmp
      end do
    end do
  end subroutine decompose
end module lu

program lu_decomposition
  use const
  use lu
  implicit none
  character(9), parameter :: F_INP = "input.txt"
  integer(SP),  parameter :: UID   = 10
  integer(SP) :: n_row, n_col, i, j, k, n
  real(DP)    :: tmp
  real(DP), allocatable :: a(:, :)

  ! ファイル OPEN
  open (UID, file = F_INP, status = 'old')

  ! 行数・列数取得
  read (UID, *) n_row, n_col

  ! 配列用メモリ確保
  allocate(a(n_row, n_col))

  ! 行列取得
  do i = 1, n_row
    read (UID, *) a(i, :)
  end do
  call display(a)

  ! ファイル CLOSE
  close (UID)

  ! LU 分解
  call decompose(a)
  print *, "-->"
  call display(a)

  ! 配列用メモリ解放
  deallocate(a)

  stop
contains

  subroutine display(a)
    implicit none
    real(DP), intent(in) :: a(:, :)
    integer(SP)  :: i, j, n
    character(8) :: f

    n = int(sqrt(real(size(a))))
    f = "(IF8.2)"
    write (f(2:2), '(I1)') n
    do i = 1, n
      write (*, f) a(i, :)
    end do
  endsubroutine display
end program lu_decomposition