PROGRAM EXERCISE39 IMPLICIT NONE ! DGESVを用いて逆行列と行列式を求める INTEGER, PARAMETER :: N = 10 !行列のサイズはここで決めているがFORMAT文の大きさはこのNとリンクしていない INTEGER :: i, j DOUBLE PRECISION, DIMENSION(1:N,1:N) :: Amat, Bmat, Amat_orig, Bmat_orig, checkmat DOUBLE PRECISION :: detA INTEGER, DIMENSION(1:N) :: IPIV INTEGER :: INFO ! 行列Aの初期化 DO i = 1, N DO j = 1, N Amat(i,j) = 10.0d0 - ABS(i-j) END DO END DO ! 行列Bの初期化(単位行列) Bmat(1:N,1:N) = 0.0d0 DO i = 1, N Bmat(i,i) = 1.0d0 END DO ! 行列AとBが正しく初期化されたか出力してチェック PRINT *, "Amatrix initialized" DO i = 1, N WRITE(*,'(10ES15.5)') (Amat(i,j),j = 1, N) END DO PRINT *, "Bmatrix initialized" DO i = 1, N WRITE(*,'(10ES15.5)') (Bmat(i,j),j = 1, N) END DO ! Amat, Bmatを保存(あとで使う場合) Amat_orig(1:N,1:N) = Amat(1:N,1:N) Bmat_orig(1:N,1:N) = Bmat(1:N,1:N) ! 一般行列の実線形方程式計算ルーチンDGESVの呼び出し (AX=B) CALL DGESV(N,N,Amat,N,IPIV,Bmat,N,INFO) ! エラー処理。正しく計算が完了するとINFOにゼロが代入されている。 IF(INFO .NE. 0) THEN WRITE(0,*) "DGESV INFO = ", INFO STOP END IF ! LU分解の結果がAmatに上書きされる PRINT *, "LU decomposition" DO i = 1, N WRITE(*,'(10ES15.5)') (Amat(i,j),j = 1, N) END DO ! 方程式の解(逆行列)がBmatに上書きされる PRINT *, "inverse of A matrix" DO i = 1, N WRITE(*,'(10ES15.5)') (Bmat(i,j), j = 1, N) END DO ! BmatがAmatの逆行列になっていることをチェック checkmat(1:N,1:N) = MATMUL(Amat_orig(1:N,1:N),Bmat(1:N,1:N)) - Bmat_orig(1:N,1:N) DO i = 1, N DO j = 1, N IF( ABS(checkmat(i,j)) .GT. 1.0D-10) THEN WRITE(0,*) i,j, checkmat(i,j) STOP "Error" END IF END DO END DO PRINT *, "Amat * Amatinv = I checked" ! Aの行列式はLU分解後のU行列の対角要素の積 detA = 1.0d0 DO i = 1, N detA = detA * Amat(i,i) END DO PRINT *, "det A = ", detA RETURN END PROGRAM EXERCISE39