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.

BenBoulderite · Answer 1 · 2020-02-07T11:15:18.240

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.

You are suggesting to use MKL libraries which are part of intel suites, what if these lib are not accessible ? — R. N, Feb 06 '20 at 20:36
@R.N: He's using MKL presumably because that's what he has in his environment. However, the same issues wrt single vs double precision apply to whichever LAPACK implementation one is using. — janneb, Feb 07 '20 at 08:28
@janneb: this is correct, I used MKL's LAPACK out of convenience. I didn't expect the implementation to matter for standard routines such as `?getri` and `?getrf`. Nice to have this confirmed though. — BenBoulderite, Feb 07 '20 at 11:07

LAPACK DGETRF+DGETRI failing

1 Answers1