Search code examples
fortranfortran90lapackmatrix-inversematrix-factorization

Why is this simple Fortran matrix inversion code not returning the expected value with LAPACK?

I'm trying to use the LAPACK package to compute the inverse of a matrix. For this end, I've employed the routines sgetrf and sgetri. However, when testing the code, it doesn't return the expected values for the inverse or (from what I can understand would be the expected value of) the LU decomposition. The base matrix I used for inversion was an arbitrary 4x4 matrix. It's written in Fortran 90 on Code::blocks, and uses the most recent GNU compilers and LAPACK libraries available. If it's at all relevant, I'm using Windows 10. The code is simple enough I believe it's adequate to post it here directly:

program test
    implicit none


    real, dimension(4,4)   :: R,R_inv


    integer                              :: info, ipiv
    real, dimension(4)                   :: lol


    R = reshape([real :: 1,0,0,0,          &
                         0,1.0/34.0,0,0,   &
                         0,0,1,2,          &
                         0,0,2,1]          &
                         ,shape(R), order = [2,1]                              )




R_inv = R
call sgetrf(4,4,R_inv,4,ipiv,info)

print *,info
print *, "_"
do pp = 1,4
print *, R_inv(pp,:)
end do
print *, "_"

call sgetri(2,R_inv,4,ipiv,lol,4,info)

print *,info
print *, "_"
print *, lol
print *, "_"
do pp = 1,4
print *, R_inv(pp,:)
end do
R_inv = matmul(R_inv,R)
print *, "_"
do pp = 1,4
print *, R_inv(pp,:)
end do
print *, "_"


end program

So far, here's what I noted:

  • The first info output returns 2. According to the documentation for the function sgetrf, this means that the factor U(2,2) is exactly singular. However, using a different LU factorization calculator such Matlab, I get a valid (and different) result for the factorization. Furthermore, I get a different result using an online calculator as well. This result agrees with the one obtained through Matlab.

  • The function sgetri seems to work fine.

  • This code works fine if the matrix is diagonal.

  • The LU factorized matrix seems fine for the first half of the diagonal, but becomes weird after that.

  • Using both online calculators and Matlab, the inverse of the matrix R exists and has no apparent problems. Perhaps negative values mess with the problem somehow? But I have no idea why that would be the case.

  • Switching the values from real to doubleprecision and using the respective dgetrf and dgetri routines had no effect on the result

Overall, can anyone tell what's wrong with this code? Are the assumptions about the routines called wrong? Am I misunderstanding what my results are supposed to be? If so, how would one go about calculating the inverse with the LAPACK package on Fortran?

Given I'm a bit inexperienced with Fortran, there might be issues on the way I handle the packages installation and usage. If this works on your machine, please let me know, as I might have messed up on that part.

Solution

  • You have two errors in your usage of LAPACK:

    1. ipiv must be an integer array of dimension(4), not a scalar.
    2. In your call to sgetri, the first argument (the matrix order) should be 4, not 2.

    With these changes your matrix inversion works correctly.

    Note however that your calculation is done in single precision (32-bit floating point), whereas matlab (and presumably that online LU calculator you linked to as well) by default calculates using 64-bit reals. So your code is more prone to precision issues, particularly for ill-conditioned problems.