CODEWITHSUNDEEP
matrixfortranfortran90
Laurence_jj
Laurence_jj

Reputation: 726

LAPACK - DGETRF fails (inverting matrices)

I am trying to calculate the inverse of a matrix, but DGETRF keeps saying that the matrix is numerically singular even when the matrix isn't.

real :: testM(1:2, 1:2), workT(2)
integer :: ipivT(2), info
testM = reshape((/4,2,7,6/), shape(testM))

WRITE(*,*) "started matrix inversion"

! DGETRF computes an LU factorization of a general M-by-N matrix A
! using partial pivoting with row interchanges.
call DGETRF(2, 2, testM, 2, ipivT, info)

WRITE(*,*) "info = ", 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(2, testM, 2, ipivT, workT, 2, info)

if (info /= 0) then
   stop 'Matrix inversion failed!'
end if
WRITE(*,*) testM
WRITE(*,*) "Matrix success"

I get the error saying that the matrix is numerically singular (info = 2). The matrix however is not numerically singular and I don't know why the error is given.

In the end I want to scale the code to a (33,33) matrix. But I'm first trying to get it working for this (2,2) matrix.

Upvotes: 0

Views: 816

Answers (1)

albert
albert

Reputation: 9057

The documentation of dgetrf (http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga0019443faea08275ca60a734d0593e60.html) says:

subroutine dgetrf   (   integer     M,
        integer     N,
        double precision, dimension( lda, * )   A,
        integer     LDA,
        integer, dimension( * )     IPIV,
        integer     INFO 
    )

So a double precision variable A but you call the routine with a real variable testM.

Similar problems will probably occur for DGETRI.

Upvotes: 2

Related Questions