C
C  INVERSE_MATRIX.F  VER 1.2  NOVEMBER 25, 2011
C
C  BY TOM IRVINE   EMAIL: tomirvine@aol.com
C
C23456789
      PROGRAM INVERSE_MATRIX
      IMPLICIT NONE
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: A, Ainv
      INTEGER :: i,j, LWORK

      INTEGER	     INFO, LDA,	M, N
      INTEGER, ALLOCATABLE, DIMENSION(:) :: IPIV,WORK

      INTEGER DeAllocateStatus

      CHARACTER :: A_filename*20

      external DGETRF
      external DGETRI


      PRINT '(" ")'
      PRINT*,"Enter size of square matrix"
      READ*,N

      LDA = N
      LWORK = N*N
      ALLOCATE (A(LDA,N))
      ALLOCATE (Ainv(LDA,N))
      ALLOCATE (WORK(LWORK))
      ALLOCATE (IPIV(N))



      PRINT '("")'
      PRINT '( "Enter the filename name of matrix A")'
      READ(*,*)A_filename

      OPEN(UNIT=10,FILE=A_filename,STATUS='OLD',ACTION='READ')

      DO i=1,N
            READ(10,*) (A(i,j),j=1,N)
      ENDDO

      PRINT '(" ")'
      PRINT*,"Input Matrix:"

      DO i = 1, N
         WRITE(*,33)(A(i,j), j = 1, N)
      END DO
33    FORMAT(100(E10.4,1X))

      PRINT '(" ")'


C
C     DGETRF computes an LU factorization of a general M-by-N matrix A
C     using partial pivoting with row interchanges.

      M=N
      LDA=N

C  Store A in Ainv to prevent it from being overwritten by LAPACK

      Ainv = A

      CALL DGETRF( M, N, Ainv, LDA, IPIV, INFO )

      IF(INFO.EQ.0)THEN
         PRINT '(" LU decomposition successful ")'
      ENDIF
      IF(INFO.LT.0)THEN
         PRINT '(" LU decomposition:  illegal value ")'
         STOP
      ENDIF
      IF(INFO.GT.0)THEN
         WRITE(*,35)INFO,INFO
35       FORMAT( 'LU decomposition: U(',I4,',',I4,') = 0 ')
      ENDIF


      PRINT '(" ")'
      PRINT*,"LU Matrix:"
      DO i = 1, N
         WRITE(*,33)(Ainv(i,j), j = 1, N)
      END DO

C  DGETRI computes the inverse of a matrix using the LU factorization
C  computed by DGETRF.

      CALL DGETRI(N, Ainv, N, IPIV, WORK, LWORK, INFO)


      PRINT '(" ")'
      IF (info.NE.0) THEN
         stop 'Matrix inversion failed!'
      ELSE
         PRINT '(" Inverse Successful ")'
      ENDIF


      PRINT '(" ")'
      PRINT*,"Inverse Matrix:"

      DO i = 1, N
         WRITE(*,33)(Ainv(i,j), j = 1, N)
      END DO

      PRINT '(" ")'

      DEALLOCATE (A, STAT = DeAllocateStatus)
      DEALLOCATE (Ainv, STAT = DeAllocateStatus)
      DEALLOCATE (IPIV, STAT = DeAllocateStatus)
      DEALLOCATE (WORK, STAT = DeAllocateStatus)

      PRINT '(" Done ")'


      PRINT '(" ")'

      STOP
      END