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:

`ipiv`

must be an integer array of`dimension(4)`

, not a scalar.- 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.

- Why am I getting this error? <class 'TypeError'>: wrong type
- GNU Fortran - Function 'dcosd' has no IMPLICIT type
- f2py does not properly import, despite successfully compiling
- Error in Python trying to create a list of a certain length after a call to Fortran to get the length
- How to get command line arguments of unknown length in Fortran?
- implicit real - complex conversion in fortran
- Optimal way to create once, then frequently access to a large array in Fortran
- How do you iterate through an array in Fortran?
- Why Intel Fortran + Intel MPI report warn(error) when using MPI_Bcast?
- Convert c_int to default kind integer
- Program in Fortran print different results with each execution
- How does work the OpenMP "nonmonotonic:dynamic" schedule?
- GFortran error: ld: library not found for -lSystem when trying to compile
- Write unformatted (binary data) to stdout
- How to call Fortran's FINDLOC within numpy's f2py
- QR factorization in Fortran
- Distribution of for loop iterations over triangular matrix with MPI
- Why MPI_REDUCE shows different number at some array locations?
- Conditional compilation in gfortran
- Fortran with Sparse BLAS not flushing memory
- Standard conforming way to get command line arguments in FORTRAN77
- Fortran C interface for program containing non-c-interoperable derived type
- Processing a shared array by a passed index in a subroutine in a parallel loop
- Gfortran type mismatch error despite "-fallow-argument-mismatch" flag
- Removing whitespace in string
- Pointer to OpenBLAS subroutines in fortran
- Valgrind complains reading from a file
- Use FP exception traps (-ffpe-trap/-fpe0) for code linked against SIGFPE-unsafe library (libxml2)
- Reading missing data from a file
- Do most compilers optimize MATMUL(TRANSPOSE(A),B)?