Gauss-Jordan 法¶
ガウス-ジョルダン法 ( Gauss-Jordan Method )は、連立1次方程式の解くためのアルゴリズムである．解くべき連立1次方程式を
と表す．右辺行列で表現した係数と左辺ベクトル(もしくは行列)をまとめた表式を拡大係数行列と呼ぶ．ここでの Gauss-Jordan 法は、拡大係数行列に対する一連の変形操作により、解を得る手法を意味する．一連の変形操作として、 行基本変形 を用いる場合と 列基本変形 を用いる場合がある．例えば、行基本変形を用いる場合、もし行列 A が正則であれば、 
の逆行列 
を用いて、連立一次方程式は
へ変形することができる．この際に、使用する行基本変形は以下の3つである．
行列Aの i 行目と j 行目を入れ替える．
行列Aの i 行目を定数倍する．
行列Aの i 行目を j 行目に加える．
行に対する操作を列に対する操作へ置き換えると、列基本変形である．
Gauss-Jordan 法は、行基本変形によって方程式 
を 
へと変形する手法であり、連立一次方程式の解法にあたる． これは、掃き出し法 などとも呼ばれ、本質的には、中学校で習う連立方程式の解法と同じ手法である．右辺をベクトルとした拡大係数行列の場合、連立一次方程式を解くことができ、右辺を単位行列とした拡大係数行列の場合、逆行列を求めることができる．以下に、逆行列解法の数値的アルゴリズムについて解説する．
Gauss-Jordan法の数値的アルゴリズム¶
Gauss-Jordan法のアルゴリズムの種類として、 pivoting ( ピボット選択 ) の 種類によって、以下の 2つがある．
部分 pivoting ( full pivoting 行・列、いづれかの方向で pivot 探索する )
完全 pivoting ( full pivoting 行列両方向同時に pivot 探索する )
簡単のため、部分 pivoting のアルゴリズムに対して、以下に流れを示す．
行列 A を 行列 U へコピーし、連立一次方程式の解を求める場合は右辺ベクトル、逆行列を求める場合は、右辺単位行列を用意する．
i=1,2,3,...,n 行目に対して、以下の操作を繰り返す．
i行目の対角要素
を 
( j=i,i+1,...,n ) のうち、最大の値をとるものとして探索する．探索後、初期の行と異なっていれば、2つの行を交換する．i 行目を
倍して、対角要素を 1 とする．i行目以外の全ての j=1~n 行目に対して、 i 列目要素を 0 とする．つまり、i 行目を
倍して、j 行目へ加算する．
n 行目までに、全ての対角要素 がゼロとならずに計算できた場合、行列 A は正則行列であり、右辺行列が A の逆行列となる．
Gauss-Jordan法のサンプルコード¶
Gauss-Jordan 法のサンプルコードについて、以下に示す．
1! ====================================================== !
2! === Gauss Jordan Inverting Example === !
3! ====================================================== !
4program main
5 use gausJordaMod, only : gaussJordan
6 implicit none
7 character(100) , parameter :: AmatFile = 'dat/Amat.dat'
8 character(100) , parameter :: BmatFile = 'dat/Bmat.dat'
9 double precision, allocatable :: Amat(:,:), Bmat(:,:), Cmat(:,:)
10 integer :: i, j, nSize
11 integer , parameter :: lun = 50
12
13 ! ------------------------------------- !
14 ! --- [1] Data Load --- !
15 ! ------------------------------------- !
16 open(lun,file=trim(AmatFile),status='old',form='formatted')
17 read(lun,*) nSize, nSize
18 allocate( Amat(nSize,nSize), Bmat(nSize,nSize), Cmat(nSize,nSize) )
19 Amat(:,:) = 0.d0
20 Bmat(:,:) = 0.d0
21 Cmat(:,:) = 0.d0
22 do i=1, nSize
23 read(lun,*) Amat(i,:)
24 enddo
25 close(lun)
26
27 ! ------------------------------------- !
28 ! --- [2] Gauss Elimination --- !
29 ! ------------------------------------- !
30 call gaussJordan( Amat, Bmat, nSize )
31
32 ! ------------------------------------- !
33 ! --- [3] save Results --- !
34 ! ------------------------------------- !
35 open(lun,file=trim(BmatFile),status='replace',form='formatted')
36 write(lun,*) nSize
37 do i=1, nSize
38 write(lun,*) Bmat(i,:)
39 enddo
40 close(lun)
41 ! ------------------------------------- !
42 ! --- [4] confirm results --- !
43 ! ------------------------------------- !
44 do i=1, nSize
45 write(6,*) Amat(i,:), ' | ', Bmat(i,:)
46 enddo
47 Cmat = matmul( Amat, Bmat )
48 do i=1, nSize
49 write(6,*) Cmat(i,:)
50 enddo
51
52 return
53end program main
1! ====================================================== !
2! === Gauss Jordan Method :: partial pivoting ver. === !
3! ====================================================== !
4module gausJordaMod
5contains
6
7 ! ========================================================== !
8 ! === Gauss Jordan Matrix Inversion === !
9 ! ========================================================== !
10 subroutine gaussJordan( Amat, Bmat, nSize )
11 implicit none
12 integer , intent(in) :: nSize
13 double precision, intent(in) :: Amat(nSize,nSize)
14 double precision, intent(out) :: Bmat(nSize,nSize)
15 integer :: i, j, k, ipivot
16 double precision :: Dinv, Dval, vpivot, buff
17 double precision, parameter :: eps = 1.d-10
18 double precision, allocatable :: Umat(:,:)
19
20 ! ----------------------------------------- !
21 ! --- [1] Preparation --- !
22 ! ----------------------------------------- !
23 ! -- [1-1] allocate Umat -- !
24 allocate( Umat(nSize,nSize) )
25 ! -- [1-2] Copy AMat & bvec -- !
26 do j=1, nSize
27 do i=1, nSize
28 Umat(i,j) = Amat(i,j)
29 enddo
30 enddo
31 do j=1, nSize
32 do i=1, nSize
33 if ( i.eq.j ) then
34 Bmat(i,j) = 1.d0
35 else
36 Bmat(i,j) = 0.d0
37 endif
38 enddo
39 enddo
40
41 ! ----------------------------------------- !
42 ! --- [2] Upper Triangular Matrization --- !
43 ! ----------------------------------------- !
44 do k=1, nSize
45 ! -- [2-1] Partial pivot selection -- !
46 vpivot = abs( Umat(k,k) )
47 ipivot = k
48 do j=k+1, nSize
49 if ( abs( Umat(j,k) ).gt.vpivot ) then
50 vpivot = abs( Umat(j,k) )
51 ipivot = j
52 endif
53 end do
54 if ( ipivot.ne.k ) then
55 do j=k, nSize
56 buff = Umat(ipivot,j)
57 Umat(ipivot,j) = Umat(k ,j)
58 Umat(k ,j) = buff
59 enddo
60 do j=1, nSize
61 buff = Bmat(ipivot,j)
62 Bmat(ipivot,j) = Bmat(k ,j)
63 Bmat(k ,j) = buff
64 enddo
65 end if
66 ! -- [2-2] Singular Matrix Detection -- !
67 if ( abs( Umat(k,k) ).lt.eps ) then
68 write(6,*) '[gaussElimin] Amat :: Singular Matrix :: No Solution End :: @ k= ', k
69 stop
70 endif
71 ! -- [2-3] For pivot row -- !
72 ! -- a'(k,:) = a(k,:) / a(k,k) -- !
73 Dinv = 1.d0 / Umat(k,k)
74 Umat(k,k) = 1.d0
75 Umat(k,(k+1):nSize) = Dinv * Umat(k,(k+1):nSize)
76 Bmat(k, 1:nSize) = Dinv * Bmat(k, 1:nSize)
77 ! -- [2-4] For Non-pivot row -- !
78 ! -- a'(j,:) = a(j,:) - a(j,k) * a(k,:) -- !
79 do j=1, nSize
80 if ( j.ne.k ) then
81 Dval = Umat(j,k)
82 Umat(j,k:nSize) = Umat(j,k:nSize) - Dval * Umat(k,k:nSize)
83 Bmat(j,1:nSize) = Bmat(j,1:nSize) - Dval * Bmat(k,1:nSize)
84 endif
85 enddo
86 end do
87
88 return
89 end subroutine gaussJordan
90
91end module gausJordaMod
Gauss-Jordan法のサンプル実行結果¶
実行結果を以下に示す．1段目に与えた行列Aとルーチンによって得られた逆行列を示している．2段目には、得た逆行列がAの逆行列であることを示すために、 
を計算している．
