Fortran - 連立方程式解法（ガウスの消去法（ピボット選択））！
Updated:
かつて、連立方程式を「ガウスの消去法」で解くアルゴリズムを Fortran95 で実装したことを紹介しました。
しかし、計算途中で対角成分がゼロになるケースでは計算ができませんでした。
今回はその問題を解決すべく、「ガウスの消去法（ピボット選択）」で解くアルゴリズムを実装してみました。
ちなみに、前回は Ruby で実装しました。
0. 前提条件
- Debian GNU/Linux 11.3 での作業を想定。
- GCC 11.2.0 (GFortran 11.2.0) でのコンパイルを想定。
- 連立方程式の解法（ガウスの消去法（ピボット選択））についての説明は割愛。（Web 上等で容易に確認可能）
1. ソースコードの作成
File: gauss_elimination_pivot.f95
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
!************************************************************
! 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
2. ソースコードのコンパイル
$ gfortran -Wall -O2 -o gauss_elimination_pivot gauss_elimination_pivot.f95
3. 動作確認
実行すると、元の数を問われるので入力し、エンター。
そして、1行ずつ係数（元数分）＋定数をスペースで区切って入力してエンター。
$ ./gauss_elimination_pivot
n ? 4
row(1,1:5) ? 1 2 7 6 6
row(2,1:5) ? 2 4 4 2 2
row(3,1:5) ? 1 8 5 2 12
row(4,1:5) ? 2 4 3 3 5
Coefficients:
1.00 2.00 7.00 6.00 6.00
2.00 4.00 4.00 2.00 2.00
1.00 8.00 5.00 2.00 12.00
2.00 4.00 3.00 3.00 5.00
Answer:
-3.00 2.00 -1.00 2.00
（もしくは、値のみを記述したテキストファイルを読み込ませてもよい）
ちなみに、上の例は、ピボット選択しない場合には計算ができないケース。
試しにソースコードの 63行目(call pivot...) をコメントアウトすると、以下のように計算できなくなる。
$ ./gauss_elimination_pivot
n ? 4
row(1,1:5) ? 1 2 7 6 6
row(2,1:5) ? 2 4 4 2 2
row(3,1:5) ? 1 8 5 2 12
row(4,1:5) ? 2 4 3 3 5
Coefficients:
1.00 2.00 7.00 6.00 6.00
2.00 4.00 4.00 2.00 2.00
1.00 8.00 5.00 2.00 12.00
2.00 4.00 3.00 3.00 5.00
Answer:
NaN NaN NaN NaN
以上。
Comments