Skip to content

Instantly share code, notes, and snippets.

@komasaru
Last active January 16, 2023 14:51
Show Gist options
  • Save komasaru/a5fe45b401043895eba9327dd27a1367 to your computer and use it in GitHub Desktop.
Save komasaru/a5fe45b401043895eba9327dd27a1367 to your computer and use it in GitHub Desktop.
Fortran 95 source code to do LU-decomposition by Crout method.
!************************************************************
! 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment