LAPACK DGETRF+DGETRI failing

Question

I'm trying to inverse a matrix via LAPACK's DGETRF and DGETRI routines, but the following code:

program Tester

    !use LapackMatrixOps
    use MathematicalResources

    implicit none

    real :: B(2, 2), A(2, 2), WORK(2)

    integer :: i, j, SIZ, IPIV(2), INFO

    SIZ=2

    do i=1, SIZ
        do j=1, SIZ
            !if(i==j) then
            !    A(i, j)=1
            !else
            !    A(i, j)=0
            !end if
            A(i, j)=rand(0)/2+0.5
            B(i, j)=0
        end do
    end do

    B=A

    call PrintMatrix(A)

    call dgetrf(size(B, 1), size(B, 2), B, size(B, 1), IPIV, INFO)
    print *, "========="
    call PrintMatrix(B)
    print *, IPIV
    call dgetri(size(B, 1), B, size(B, 1), IPIV, WORK, size(WORK), INFO);

    print *, "========="

    call PrintMatrix(B)

end program Tester

Compiled with:

CC=gfortran
CFLAGS=-O2 -W -g
LFLAGS=-L/usr/lib/lapack -llapack -lblas

TARGET=Tester

all: install clean

.PHONY: install
install:
    $(CC) *.f90 $(CFLAGS) $(LFLAGS) -o $(TARGET)

.PHONY: test
test:
    $(CC) *.f90 $(CFLAGS) $(LFLAGS) -o $(TARGET)
    #==============TEST=============#
    ./$(TARGET)

.PHONY: clean
clean:
    rm -rf ./*.o

Seems to be failing. When I run make test, I get:

gfortran *.f90 -O2 -W -g -L/usr/lib/lapack -llapack -lblas -o Tester
#==============TEST=============#
./Tester
  0.500003815      0.565768838    
  0.877802610      0.729325056    
 =========
  0.500003815       1.89445039E+28
  0.877802610       1.57469535    
           1           2
 =========
              NaN   6.86847806E-20
              NaN  -1.28451300    
pedro@pedro-300E5M-300E5L:~/Área de Trabalho/AED/openpypm2$ 

What am I doing wrong? I've tried both custom-compiled LAPACK libraries and the ATLAS library binaries, and they led to the same results. I've also tried to change certain parameters such as the length of the work array and all I've done kept leading to this.

Answer 1

You define single precision arrays but you call double precision routines. Call sgetrf and sgetri instead. Alternatively, if you want double precision, declare your variables with:

integer, parameter :: dp = kind(1.d0)
real(kind=dp) :: B(2,2), A(2,2), work(2)

Anyway, using single precision, the following code:

Program Tester

    implicit none
    real :: B(2, 2), A(2, 2), WORK(2)
    integer :: i, j, SIZ, IPIV(2), INFO

    SIZ=2

    do i=1, SIZ
        do j=1, SIZ
            A(i, j)=rand(0)/2+0.5
            B(i, j)=0
        end do
    end do

    B=A

    print *, A
    call sgetrf(size(B, 1), size(B, 2), B, size(B, 1), IPIV, INFO)
    print *, "========="
    print *, B
    print *, IPIV
    call sgetri(size(B, 1), B, size(B, 1), IPIV, WORK, size(WORK), INFO);
    print *, "========="
    print *, B

end program Tester

compiled with the following makefile:

FC =gfortran

MKL_LIB    = -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl -m64 -I${MKLROOT}/include

all: test.f90
        $(FC) -o test.exe $^ $(MKL_LIB) -I$(MKLROOT)/include

yield this ouput:

  0.500003815      0.877802610      0.565768838      0.729325056    
 =========
  0.877802610      0.569608510      0.729325056      0.150339082    
           2           2
 =========
  -5.52652836       6.65163040       4.28716612      -3.78882527   

N.B.: As noted in the comments, I have used MKL's implementation of LAPACK out of convenience. Regardless of the implementation, the precision of the routines and the variables should match. What happens whenever they don't may depend on the implementation, but is very likely to be junk in all cases.