- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Please ask this in the MKL forum. It isn't a compiler question.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page