[線形代数] LAPACKのDGESVで連立方程式を解く

連立方程式の解法

FORTRAN で LAPACK の連立方程式のソルバ DGESV を使って連立方程式を解く。

今回解く方程式は keisanサイト様 より (解は $1, -1, 2, -2$)。

$$ Ax=B \rightarrow \begin{bmatrix} 1&1&1&1 \\ 1&1&1&-1 \\ 1&1&-1&1 \\ 1&-1&1&1 \\ \end{bmatrix} \begin{bmatrix} x \\ y \\ z \\ w \\ \end{bmatrix}= \begin{bmatrix} 0 \\ 4 \\ -4 \\ 2 \\ \end{bmatrix} $$

ソースコード

• 変数定義
• LAPACKの DGESV (詳しくはコチラ) に必要な値の決定
• DGESV を CALL (=連立方程式を解く)
• 結果の表示

PROGRAM RENRITU
    IMPLICIT NONE
    
    REAL(8) :: A(4,4), B(4,1), ANS(4,1)     ! Ax=B, x=ANS
    INTEGER :: N, NRHS, LDA, LDB, INFO      ! LAPACK
    INTEGER, ALLOCATABLE :: IPIV(:)         ! LAPACK

    ! LAPACKの手続き
    N = size(A,2)                           ! 行列(NxN)の次数N
    NRHS = 1                                ! Bの列数 (=1)
    LDA = size(A,1)                         ! Aの行数
    LDB= size(B,1)                          ! Bの行数
    ALLOCATE( IPIV(N) )                     ! P行列 

    ! 各係数の定義 (Ax=B)
    A(1,1)=1.d0 ; A(1,2)=1.d0  ; A(1,3)=1.d0  ; A(1,4)=1.d0  ; B(1,1)=0.d0
    A(2,1)=1.d0 ; A(2,2)=1.d0  ; A(2,3)=1.d0  ; A(2,4)=-1.d0 ; B(2,1)=4.d0
    A(3,1)=1.d0 ; A(3,2)=1.d0  ; A(3,3)=-1.d0 ; A(3,4)=1.d0  ; B(3,1)=-4.d0
    A(4,1)=1.d0 ; A(4,2)=-1.d0 ; A(4,3)=1.d0  ; A(4,4)=1.d0  ; B(4,1)=2.d0

    ! SOLVE
    ANS(:,:)=B(:,:)
    CALL DGESV( N, NRHS, A, LDA, IPIV, ANS, LDB, INFO )

    ! 結果表示
    IF ( INFO==0 ) WRITE(*,'(ES12.4)') ANS(:,:)

END PROGRAM RENRITU

コンパイル・実行

gfortran renritu44.f90 -llapack
./a.out

結果

  1.0000E+00
 -1.0000E+00
  2.0000E+00
 -2.0000E+00