topic Please ask this in the MKL in Intel® oneAPI Math Kernel Library

Matrix Inversion using Lapack (DGETRF/DGETRI)

Neto__Conrado — Thu, 12 Dec 2019 19:30:18 GMT

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

Please ask this in the MKL

Steve_Lionel — Tue, 17 Dec 2019 00:49:55 GMT

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

As Steve wrote, this is a

mecej4 — Tue, 17 Dec 2019 17:02:12 GMT

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

or look at the dgetrix.f

Gennady_F_Intel — Mon, 23 Dec 2019 05:20:08 GMT

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