[線形代数] 逆行列を求める

逆行列の計算

FORTRANでLAPACKを使って4x4行列の逆行列を求める。

今回は以下の行列の逆行列を求める。答えも載せておきます(計算サイト様)。

$$ A= \begin{bmatrix} 3&1&1&2 \\ 5&1&3&4 \\ 2&0&1&0 \\ 1&3&2&1 \\ \end{bmatrix} , \hspace{5mm} A^{-1}= \begin{bmatrix} 0.5 & -0.23 & 0.36 & -0.01 \\ 0.5 & -0.32 & -0.09 & 0.27 \\ -1 & 0.45 & 0.27 & 0.18 \\ 0 & 0.27 & -0.63 & -0.09 \\ \end{bmatrix} $$

計算

• LAPACKに必要な入力パラメタを決定
• 行列の定義
• 逆行列の計算 (LAPACK)
• 結果の表示

任意の行列は直接逆行列を出せないので、LU分解してから求める。

具体的には、LU分解 (DGETRF)してから逆行列の計算 (DGETRI)をする。

必要な入力パラメータは上記サイトを、LU分解についてはコチラを参考に。

ソースコードは以下の通り。

PROGRAM INVMAT
IMPLICIT NONE

REAL(8) :: A(4,4), invA(4,4)        ! 4x4行列
INTEGER :: M, N                     ! A の行・列
INTEGER :: LDA                      ! = m
INTEGER, ALLOCATABLE :: IPIV(:)     ! サイズ n の整数配列
INTEGER :: INFO                     ! =0 で正常終了
INTEGER :: NB=64                    ! 最大次元
INTEGER :: LWORK                    ! =n*NB
REAL(8), ALLOCATABLE :: WORK(:)     ! サイズ LWORK の実数配列


! LAPACK関係のもろもろ
m = size(A,1)                       ! 行数の取得
n = size(A,2)                       ! 列数の取得
LDA = m                             ! LDA=m
LWORK = n*NB                        ! 
ALLOCATE( IPIV(n), WORK(LWORK) )    ! アロケート


! A の決定
A(1,1)=3.d0 ; A(1,2)=1.d0 ; A(1,3)=1.d0 ; A(1,4)=2.d0
A(2,1)=5.d0 ; A(2,2)=1.d0 ; A(2,3)=3.d0 ; A(2,4)=4.d0
A(3,1)=2.d0 ; A(3,2)=0.d0 ; A(3,3)=1.d0 ; A(3,4)=0.d0
A(4,1)=1.d0 ; A(4,2)=3.d0 ; A(4,3)=2.d0 ; A(4,4)=1.d0


! 逆行列の計算 (LU分解 => 逆行列)
invA = A
    CALL DGETRF( m, n, invA, LDA, IPIV(1:min(m,n)), INFO )
    CALL DGETRI( n, invA, LDA, IPIV, WORK, LWORK, INFO )
IF (INFO == 0) WRITE(*,*) "CALCLATION FINISHED"


! A の表示
WRITE(*,*) "A = "
WRITE(*,'(4ES10.3)') A(1,:)
WRITE(*,'(4ES10.3)') A(2,:)
WRITE(*,'(4ES10.3)') A(3,:)
WRITE(*,'(4ES10.3)') A(4,:)

! 逆行列の表示
WRITE(*,*) "INVERSE A = "
WRITE(*,'(4ES10.3)') invA(1,:)
WRITE(*,'(4ES10.3)') invA(2,:)
WRITE(*,'(4ES10.3)') invA(3,:)
WRITE(*,'(4ES10.3)') invA(4,:)

END PROGRAM INVMAT

コンパイル・実行

$ gfortran invmat.f90 -llapack
$ ./a.out

結果

 CALCLATION FINISHED

 A = 
 3.000E+00    1.000E+00    1.000E+00    2.000E+00
 5.000E+00    1.000E+00    3.000E+00    4.000E+00
 2.000E+00    0.000E+00    1.000E+00    0.000E+00
 1.000E+00    3.000E+00    2.000E+00    1.000E+00

 INVERSE A = 
 5.000E-01   -2.273E-01    3.636E-01   -9.091E-02
 5.000E-01   -3.182E-01   -9.091E-02    2.727E-01
-1.000E+00    4.545E-01    2.727E-01    1.818E-01
 1.413E-16    2.727E-01   -6.364E-01   -9.091E-02

まとめ

4x4以下は計算サイト様を使えばよいが、5x5以上で他力本願は厳しい。 今回は4x4だが、拡張は簡単にできるので汎用性はありそう。