Fortran - 行列式の計算（余因子展開による）！
Updated:
Fortran 95 で余因子展開による行列式の計算を行ってみました。
少し前に、同じことを Ruby で Array クラスを拡張する方法で実装しています。
0. 前提条件Permalink
- Debian GNU/Linux 10.3 (64bit) での作業を想定。
- GCC 9.2.0 (GFortran 9.2.0) でのコンパイルを想定。
1. 行列式の余因子展開についてPermalink
n(>1) 次正方行列 A=(aij) から第 i 行と第 j 列の成分をすべて取り除いて得られる n−1 次行列の行列式に、 (−1)i+j を掛けたものを aij の 余因子 といい、 ˜Aij で表す。すなわち、
˜Aij=(−1)i+j|a11⋯a1 j−1a1 j+1⋯a1n⋮⋮⋮⋮ai−1 1⋯ai−1 j−1ai−1 j+1⋯ai−1 nai+1 1⋯ai+1 j−1ai+1 j+1⋯ai+1 n⋮⋮⋮⋮an1⋯an j−1an j+1⋯ann|そして、次の定理が成り立つ。（証明略）
【定理】 n 次正方行列 A=(aij) に対して、
|A|=ai1˜Ai1+ai2˜Ai2+⋯+ain˜Ain (1<i≤n)|A|=a1j˜A1j+a2j˜A2j+⋯+anj˜Anj (1<j≤n)(2)を行列式｜A｜の i 行についての 余因子展開、(3)を行列式｜A｜の j 列についての 余因子展開 という。（余因子展開は、ラプラス展開や単に展開と呼ばれることもある）
2. ソースコードの作成Permalink
File: determinant.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
!****************************************************
! 行列式の計算（余因子展開）
!
! DATE AUTHOR VERSION
! 2019.12.23 mk-mode.com 1.00 新規作成
!
! Copyright(C) 2013 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 comp
use const
implicit none
private
public :: det
contains
! 行列式計算
! * 再帰的
!
! :param(in) real(8) m(:, :): 行列
! :return real(8) d: 行列式
recursive real(DP) function det(m)
implicit none
real(DP), intent(in) :: m(:, :)
integer(SP) :: j, n
det = 0.0_DP
n = size(m(1, :))
select case(n)
case (1)
det = m(1, 1)
case (2)
det = m(1, 1) * m(2, 2) &
& - m(1, 2) * m(2, 1)
case (3)
det = m(1, 1) * m(2, 2) * m(3, 3) &
& + m(1, 2) * m(2, 3) * m(3, 1) &
& + m(1, 3) * m(2, 1) * m(3, 2) &
& - m(1, 1) * m(2, 3) * m(3, 2) &
& - m(1, 2) * m(2, 1) * m(3, 3) &
& - m(1, 3) * m(2, 2) * m(3, 1)
case (4:)
det = 0
do j = 1, n
det = det + m(1, j) * cof(m, 1, j)
end do
end select
end function det
! 余因子（正方行列 A の (i, j) に対する）
!
! :param(in) m: 行列配列
! :param(in) i: 行インデックス
! :param(in) j: 列インデックス
! :param(out) c: 余因子
real(DP) function cof(m, i, j)
implicit none
real(DP), intent(in) :: m(:, :)
integer(SP), intent(in) :: i, j
integer(SP) :: n
real(DP), allocatable :: m2(:, :)
n = size(m(1, :))
allocate(m2(n - 1, n - 1))
m2(1:i-1, 1:j-1) = m(1:i-1, 1:j-1)
m2(1:i-1, j:n ) = m(1:i-1, j+1:n)
m2(i:n , 1:j-1) = m(i+1:n, 1:j-1)
m2(i:n , j:n ) = m(i+1:n, j+1:n)
cof = (-1)**(i+j) * det(m2)
deallocate(m2)
end function cof
end module comp
program determinant
use const
use comp
implicit none
integer(SP) :: n, i
real(DP) :: d
real(DP), allocatable :: m(:, :)
! データ数読み込み
read (*, *) n
! 配列用メモリ確保
allocate(m(n, n))
! データ読み込み
do i = 1, n
read (*, *) m(i, :)
end do
print *, "Matrix A ="
do i = 1, n
print *, m(i, :)
end do
! 行列式計算
d = det(m)
print *, "det(A) =", d
! 配列用メモリ解放
deallocate(m)
end program determinant
3. ソースコードのコンパイルPermalink
$ gfortran -Wall -O2 -o determinant determinant.f95
4. 動作確認Permalink
まず、以下のような入力ファイルを用意する。
（先頭行：行数、2行目以降：行列の各要素）
File: input.txt
1
2
3
4
5
6
5
1 3 5 7 9
4 2 8 6 0
9 3 7 5 1
4 0 6 8 2
3 6 9 2 5
そして、実行。
$ ./determinant < input.txt
Matrix A =
1.0000000000000000 3.0000000000000000 5.0000000000000000 7.0000000000000000 9.0000000000000000
4.0000000000000000 2.0000000000000000 8.0000000000000000 6.0000000000000000 0.0000000000000000
9.0000000000000000 3.0000000000000000 7.0000000000000000 5.0000000000000000 1.0000000000000000
4.0000000000000000 0.0000000000000000 6.0000000000000000 8.0000000000000000 2.0000000000000000
3.0000000000000000 6.0000000000000000 9.0000000000000000 2.0000000000000000 5.0000000000000000
det(A) = 2320.0000000000000
当然、以下のように実行してもよいし、
$ cat input.txt | ./determinant
以下のように実行してもよい。
$ ./determinant <<_EOF_
5
1 3 5 7 9
4 2 8 6 0
9 3 7 5 1
4 0 6 8 2
3 6 9 2 5
_EOF_
以上。
Comments