IMSL를 이용한 역행렬 계산

---

!DEC$ OBJCOMMENT LIB:'libguide.lib'

      SUBROUTINE INV_MATRIX(A, INV_A, N)

       INCLUDE 'link_fnl_static.h'

     

       USE lin_sol_gen_int

       IMPLICIT NONE

       INTEGER N, i, j

       DOUBLEPRECISION A,INV_A, B, X, DET, INV_DET, RES

      

       DOUBLEPRECISION  ONE, ERR , NRHS

    

       PARAMETER ( ONE = 1.0 )    

    

       DIMENSION  A(N,N), INV_A(N,N),B(N,0), X(N,0), DET(2), INV_DET(2), RES(N,N)

       CALL D_lin_sol_gen(A, B, X, nrhs=0, ainv=INV_A, det=DET)

              ! Compute the determinant for the inverse matrix.

       CALL D_lin_sol_gen(INV_A, B, X, nrhs=0, ainv=A, det= INV_DET)

              ! Check residuals, A times inverse = Identity.

       RES = MATMUL(A,INV_A)

       DO i=1, N

        RES(i,i) = RES(i,i) - ONE

       ENDDO

       ERR = SUM(DABS(RES)) / SUM(DABS(A))

       IF (ERR <= DSQRT(EPSILON(ONE))) THEN

        IF (DET(1) == INV_DET(1) .AND.(DABS(DET(2)+INV_DET(2)) .LT. DABS(DET(2))*DSQRT(EPSILON(ONE)))) THEN

              WRITE (*,*) 'INVERSE MATRIX is correct.'

        END if

       END if

    

       ENDSUBROUTINE INV_MATRIX

---