Skip to content

Instantly share code, notes, and snippets.

@komasaru komasaru/sle_lu.f95
Last active Jun 4, 2019

Embed
What would you like to do?
Fortran 95 source code to solve simultaneous equations by LU-decomposition(outer-product form).
!************************************************************
! 連立1次方程式の解法 ( LU 分解(外積形式ガウス法) )
!
! date name version
! 2019.03.13 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, solve
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
! 連立方程式を解く
!
! :param(in) real(8) a(n,n): 行列 A
! :param(in) real(8) b(n): ベクトル b
! :param(out) real(8) x(n): 解ベクトル x
subroutine solve(a, b, x)
implicit none
real(DP), intent(in) :: a(:, :), b(:)
real(DP), intent(out) :: x(:)
integer(SP) :: n, i, j
real(DP) :: s
real(DP), allocatable :: y(:)
! 元数
n = size(x)
! 配列 y のメモリ確保
allocate(y(n))
! 前進代入
! * Ly = b から y を計算
do i = 1, n
s = sum((/(a(i, j) * y(j), j=1,i)/))
y(i) = b(i) - s
end do
! 後退代入
! * Ux = y から x を計算
do i = n, 1, -1
s = sum((/(a(i, j) * x(j), j=i+1,n)/))
x(i) = (y(i) - s) / a(i, i)
end do
! 配列 y のメモリ解放
deallocate(y)
end subroutine solve
end module lu
program sle_lu
use const
use lu
implicit none
character(9), parameter :: F_INP = "input.txt"
integer(SP), parameter :: UID = 10
integer(SP) :: n, i ! 元数、ループインデックス
real(DP), allocatable :: a(:,:), b(:) ! 係数配列
real(DP), allocatable :: x(:) ! 解配列
! IN ファイル OPEN
open(UID, file = F_INP, status = "old")
! 元数取得
read(UID, *) n
if (n < 1) stop
print '("n = ", I0)', n
! 配列メモリ確保
allocate(a(n, n))
allocate(b(n))
allocate(x(n))
! 行列 A 取得
do i = 1, n
read(UID, *) a(i,:)
end do
print '(A)', "A ="
call display_mtx(a)
! ベクトル B 取得
read(UID, *) b(:)
print '(A)', "b ="
call display_vec(b)
! IN ファイル CLOSE
close (UID)
! 行列 A の LU 分解
call decompose(a)
!print '(A)', "(LU) ="
!call display_mtx(a)
! 連立方程式を解く
call solve(a, b, x)
print '(A)', "x ="
call display_vec(x)
! 配列メモリ解放
deallocate(a)
deallocate(b)
deallocate(x)
contains
subroutine display_mtx(m)
implicit none
real(DP), intent(in) :: m(:, :)
integer(SP) :: n
character(20) :: f ! 書式文字列
n = size(m(1, :))
write (f, '("(", I0, "(F10.2)", ")")') n + 1
do i = 1, n
print f, m(i,:)
end do
end subroutine display_mtx
subroutine display_vec(v)
implicit none
real(DP), intent(in) :: v(:)
integer(SP) :: n
character(20) :: f ! 書式文字列
n = size(v)
write (f, '("(", I0, "(F10.2)", ")")') n
print f, v(:)
end subroutine display_vec
end program sle_lu
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.