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