Skip to content

Instantly share code, notes, and snippets.

@komasaru
Created April 14, 2022 02:19
Show Gist options
  • Save komasaru/5971d1d02b355f4160aefedc4afaabcd to your computer and use it in GitHub Desktop.
Save komasaru/5971d1d02b355f4160aefedc4afaabcd to your computer and use it in GitHub Desktop.
Fortran95 source code to solve simultaneous equations with Gauss elimination method(pivot).
!************************************************************
! Simultaneous equations solving by Gauss-Elimination(Pivot) method
!
! DATE AUTHOR VERSION
! 2022.04.14 mk-mode.com 1.00 新規作成
!
! Copyright(C) 2022 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 gauss_elimination_pivot
use const
implicit none
private
public :: solve
contains
! ピボット選択
!
! :param(in) integer(4) n: 元数
! :param(in) integer(4) k: 対象行
! :param(inout) real(8) a(n,n+1): 係数配列
subroutine pivot(n, k, a)
implicit none
integer(SP), intent(in) :: n, k
real(DP), intent(inout) :: a(n, n + 1)
integer(SP) :: i, pv
real(DP) :: v_max, dummy(n + 1)
pv = k
v_max = abs(a(k, k))
do i = k, n
if (abs(a(i, k)) > v_max) then
pv = i
v_max = abs(a(i, k))
end if
end do
if (k /= pv) then
dummy = a(k, :)
a(k, :) = a(pv, :)
a(pv, :) = dummy
end if
end subroutine pivot
! 連立方程式を解く
!
! :param(in) integer(4) n: 元数
! :param(inout) real(8) a(n,n+1): 係数配列
subroutine solve(n, a)
implicit none
integer(SP), intent(in) :: n
real(DP), intent(inout) :: a(n, n + 1)
integer(SP) :: i, j
real(DP) :: d
! 前進消去
do j = 1, n - 1
call pivot(n, j, a) ! ピボット選択
do i = j + 1, n
d = a(i, j) / a(j, j)
a(i, j+1:n+1) = a(i, j+1:n+1) - a(j, j+1:n+1) * d
end do
end do
! 後退代入
do i = n, 1, -1
d = a(i, n + 1)
do j = i + 1, n
d = d - a(i, j) * a(j, n + 1)
end do
a(i, n + 1) = d / a(i, i)
end do
end subroutine solve
end module gauss_elimination_pivot
program simultaneous_equations
use const
use gauss_elimination_pivot
implicit none
integer(SP) :: n, i ! 元数、ループインデックス
character(20) :: f ! 書式文字列
real(DP), allocatable :: a(:,:) ! 係数配列
! 元数取得
write (*, '(A)', advance="no") "n ? "
read (*, *) n
! 配列メモリ確保
allocate(a(n, n + 1))
! 係数取得
do i = 1, n
write (*, '("row(", I0, ",1:", I0, ") ? ")', advance="no") i, n + 1
read (*, *) a(i,:)
end do
write (*, '(/A)') "Coefficients:"
write (f, '("(", I0, "(F8.2, 2X)", ")")') n + 1
do i = 1, n
write (*, f) a(i,:)
end do
! 連立方程式を解く
! (計算後 a の最右列が解)
call solve(n, a)
! 結果出力
write (*, '(A)') "Answer:"
write (f, '("(", I0, "(F8.2, 2X)", ")")') n
write (*, f) a(:, n + 1)
! 配列メモリ解放
deallocate(a)
end program simultaneous_equations
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment