Matrix Inversion using Lapack (DGETRF/DGETRI)

[Neto__Conrado]

I am trying to invert the following matrix using DGETRF and DGETRI

matrix = 

  0.50E+12  0.00E+00 -0.30E+09  0.25E+12  0.00E+00  0.10E-02  0.00E+00  0.00E+00  0.00E+00  0.00E+00  0.00E+00  0.00E+00
  0.00E+00  0.12E+08  0.00E+00  0.00E+00 -0.60E+07 -0.00E+00 -0.00E+00  0.00E+00  0.00E+00  0.00E+00  0.00E+00  0.00E+00
 -0.30E+09  0.00E+00  0.48E+06  0.00E+00 -0.00E+00 -0.00E+00  0.30E+09  0.00E+00  0.00E+00  0.00E+00  0.00E+00  0.00E+00
  0.25E+12  0.00E+00  0.00E+00  0.10E+13  0.00E+00 -0.00E+00  0.25E+12  0.00E+00  0.00E+00  0.00E+00  0.00E+00  0.00E+00
  0.00E+00 -0.60E+07 -0.00E+00  0.00E+00  0.12E+08 -0.00E+00  0.00E+00 -0.60E+07 -0.00E+00 -0.00E+00  0.00E+00  0.00E+00
  0.00E+00 -0.00E+00 -0.24E+06 -0.30E+09  0.00E+00  0.60E+04  0.00E+00 -0.00E+00 -0.24E+06  0.30E+09  0.00E+00  0.00E+00
  0.00E+00 -0.00E+00  0.30E+09  0.25E+12  0.00E+00 -0.00E+00  0.10E+13  0.00E+00 -0.30E+09  0.25E+12  0.00E+00  0.00E+00
  0.00E+00  0.00E+00  0.00E+00  0.00E+00 -0.60E+07 -0.00E+00  0.00E+00  0.12E+08  0.00E+00  0.00E+00 -0.60E+07 -0.00E+00
  0.00E+00  0.00E+00  0.00E+00  0.00E+00 -0.00E+00 -0.00E+00 -0.30E+09  0.00E+00  0.48E+06  0.00E+00 -0.00E+00  0.30E+09
  0.00E+00  0.00E+00  0.00E+00  0.00E+00 -0.00E+00 -0.00E+00  0.25E+12  0.00E+00  0.00E+00  0.10E+13  0.00E+00  0.25E+12
  0.00E+00  0.00E+00  0.00E+00  0.00E+00  0.00E+00 -0.10E-02  0.00E+00 -0.60E+07 -0.00E+00  0.00E+00  0.60E+07  0.00E+00
  0.00E+00  0.00E+00  0.00E+00  0.00E+00  0.00E+00 -0.10E-02  0.00E+00 -0.00E+00  0.30E+09  0.25E+12  0.00E+00  0.50E+12

however it says the matrix is numerically singular. However I can invert the matrix using excel/matlab with no problems

Does anyone have any idea why this is happening?

in the code below, the matrix is of size MJNT3xMJNT3 however I want to invert only the first NxN terms of S or Sinv.

I tried calling 

CALL DGETRF(N, N, Sinv, N, ipiv, info)

CALL DGETRI(N, Sinv, N, ipiv, work, N, info)

it doesn't work. Then I tried a few other possibilities replacing the LDA in both DGETRF and DGETRI from N to MJNT3. and it didn't work either

I am not expert, so this is quite confusing to me. 

the whole subroutine is given below:

      SUBROUTINE INVERT(S)

      INTEGER MJNT,MJNT3
      PARAMETER (MJNT=510,MJNT3=3*MJNT)

!     DECLARING GLOBAL VARIABLES

      COMMON /CLIST21/ N,NJ3,M,NMT,NJ,NRJ,NR
      INTEGER N,NJ3,M,NMT,NJ,NRJ,NR !21    

!     DECLARING LOCAL VARIABLES
      INTEGER I,J !21  

      REAL*8 S(MJNT3,MJNT3)
      REAL*8 Sinv(MJNT3,MJNT3)
      INTEGER, DIMENSION(size(Sinv,1)) :: ipiv   ! pivot indices
      INTEGER, PARAMETER :: dp = selected_real_kind(15, 307)
      real(dp), DIMENSION(size(Sinv,1)) :: work  ! work array for LAPACK
      
      INTEGER :: info
      
      EXTERNAL DGETRF
      EXTERNAL DGETRI
      

      ! Store S in Sinv to prevent it from being overwritten by LAPACK
      DO 10 I=1,N
      DO 20 J=1,N
      Sinv(I,J) = S(I,J)
20    CONTINUE
10    CONTINUE
      
      ! DGETRF COMPUTES AN LU FACTORIZATION OF A GENERAL M-BY-N MATRIX A
      ! USING PARTIAL PIVOTING WITH ROW INTERCHANGES.
      CALL DGETRF(N, N, Sinv, N, ipiv, info)
      IF (info /= 0) THEN
          STOP 'MATRIX IS NUMERICALLY SINGULAR!'
      END IF      
      ! DGETRI COMPUTES THE INVERSE OF A MATRIX USING THE LU FACTORIZATION
      ! COMPUTED BY DGETRF.
      CALL DGETRI(N, Sinv, N, ipiv, work, N, info)
      IF (info /= 0) THEN
          STOP 'MATRIX INVERSION FAILED!'
      END IF
      
      DO 30 I=1,N
      DO 40 J=1,N
      S(I,J) = Sinv(I,J)
40    CONTINUE
30    CONTINUE
      
      PAUSE !to test
      
      RETURN
      END

[Steve_Lionel]

Please ask this in the MKL forum. It isn't a compiler question.

[mecej4]

As Steve wrote, this is a topic for the MKL forum.

Since SINV is not an N X N matrix, the 4th argument to DGETRF and the 3rd argument to DGETRI are incorrect. Since N is in a common block, you need to ensure that N is set to the correct value before the subroutine is called.

Rather than trying arbitrary values for arguments, please read the MKL/Lapack documentation and pass the correct arguments to DGETRF and DGETRI.

      program calls
      real*8 s(1530,1530)
      call invert(s)
      write(*,10)((s(i,j),j=1,12),i=1,12)
   10 format(12ES11.3)
      end program
      
      subroutine invert(s)

      integer mjnt,mjnt3
      parameter (mjnt=510,mjnt3=3*mjnt)

!     declaring global variables
      integer n  

!     declaring local variables
      integer i,j !21  

      real*8 s(mjnt3,mjnt3)
      real*8 sinv(mjnt3,mjnt3)
      integer, dimension(size(sinv,1)) :: ipiv   ! pivot indices
      integer, parameter :: dp = selected_real_kind(15, 307)
      real(dp), dimension(size(sinv,1)) :: work  ! work array for lapack
      
      integer :: info
      
      external dgetrf
      external dgetri
      
      n=12
      read(*,*)((s(i,j),j=1,n),i=1,n)
      ! store s in sinv to prevent it from being overwritten by lapack
      do 10 i=1,n
      do 20 j=1,n
      sinv(i,j) = s(i,j)
20    continue
10    continue
      
      ! dgetrf computes an lu factorization of a general m-by-n matrix a
      ! using partial pivoting with row interchanges.
      call dgetrf(n, n, sinv, mjnt3, ipiv, info)
      if (info /= 0) then
          stop 'matrix is numerically singular!'
      end if      
      ! dgetri computes the inverse of a matrix using the lu factorization
      ! computed by dgetrf.
      call dgetri(n, sinv, mjnt3, ipiv, work, n, info)
      if (info /= 0) then
          stop 'matrix inversion failed!'
      end if
      
      do 30 i=1,n
      do 40 j=1,n
      s(i,j) = sinv(i,j)
40    continue
30    continue
      
      
      return
      end

[Gennady_F_Intel]

or look at the dgetrix.f  examples from mklroot/examples/lapackf/source directory and plus this link to the MKL Developer Reference - https://software.intel.com/en-us/onemkl-developer-reference-fortran-getri