Last active
May 22, 2019 09:59
-
-
Save komasaru/306a3231e4c0466ce12d84cd22231795 to your computer and use it in GitHub Desktop.
Fortran 95 source code to do LU-decomposition by outer-product form.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
!************************************************************ | |
! LU 分解(外積形式ガウス法(outer-product form)) | |
! | |
! 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 分解 | |
! * 外積形式ガウス法(outer-product form) | |
! | |
! :param(inout) real(8) a(:,:): 行列 | |
subroutine decompose(a) | |
implicit none | |
real(DP), intent(inout) :: a(:, :) | |
integer(SP) :: i, j, k, n | |
real(DP) :: tmp | |
n = int(sqrt(real(size(a)))) | |
do k = 1, n | |
if (a(1, 1) == 0.0_DP) then | |
print *, "Can't divide by 0 ..." | |
stop | |
end if | |
tmp = 1.0_DP / a(k, k) | |
a(k+1:n, k) = a(k+1:n, k) * tmp | |
do j = k + 1, n | |
tmp = a(k, j) | |
a(k+1:n, j) = a(k+1:n, j) - a(k+1:n, k) * 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