Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fortran 95 source code to do LU-decomposition by inner-product form.
!************************************************************
! LU 分解(内積形式ガウス法(inner-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 分解
! * 内積形式ガウス法(inner-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
do j = 1, k - 1
tmp = a(j, k)
a(j+1:n, k) = a(j+1:n, k) - a(j+1:n, j) * tmp
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)
a(k+1:n, k) = a(k+1:n, k) * tmp
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
You can’t perform that action at this time.