2

I'm programming in fortran and trying to use the DGETRI matrix inverter from the Lapack package:

http://www.netlib.org/lapack/explore-html/df/da4/dgetri_8f.html

But very strangely it seems to be messing with all my variables. In this very simple example, my matrix A initialised at the beginning of the program changes as DGETRI is applied, even though DGETRI doesn't involve A…

Can anybody tell me what is going on? Thanks!

PROGRAM solvelinear

REAL(8), dimension(2,2)     :: A,Ainv
REAL(8),allocatable         :: work(:)
INTEGER                     :: info,lwork,i
INTEGER,dimension(2)        :: ipiv

info=0
lwork=10000
allocate(work(lwork))
work=0
ipiv=0

A = reshape((/1,-1,3,3/),(/2,2/))
Ainv = reshape((/1,-1,3,3/),(/2,2/))

CALL DGETRI(2,Ainv,2,Ipiv,work,lwork,info)

print*,"A = "
do i=1,2
  print*,A(i,:)
end do

END PROGRAM solve linear

This is the output:

 A =
   1.0000000000000000        0.0000000000000000
  -1.0000000000000000       0.33333333333333331
Vim154
  • 23
  • 3
  • 1
    Your code is not compilable and doesn't reproduce the problem (with gfortran 4.9.1 and lapack 3.5.0). Additionally, it doesn't calculate the inverse of A, because you forgot to calculate the factorization with zgetrf (which doesn't matter with respect to your specific problem). – Stefan Oct 21 '14 at 04:49
  • Thanks you. Do you know what's wrong now with my computation of the inverse? (see below). – Vim154 Oct 21 '14 at 08:12

1 Answers1

4

You have to compute the LU decomposition before calling DGETRI. You are using the double version of the routines so the LU has to be computed with DGETRF (ZGETRF is the complex version).

The following code compiles and returns the correct values

PROGRAM solvelinear
implicit none
REAL(8), dimension(3,3)     :: A,Ainv,M,LU
REAL(8),allocatable              :: work(:)
INTEGER                        :: info,lwork
INTEGER,dimension(3)        :: ipiv
INTEGER                        :: i

info=0
lwork=100
allocate(work(lwork))
work=0
ipiv=0

A = reshape((/1,-1,3,3,1,1,1,1,1,1/),(/3,3/))

!-- LU factorisation
LU = A
CALL DGETRF(3,3,LU,3,ipiv,info)

!-- Inversion of matrix A using the LU
Ainv=LU
CALL DGETRI(3,Ainv,3,Ipiv,work,lwork,info)

!-- computation of A^-1 * A to check the inverse
M = matmul(Ainv,A)

print*,"M = "
do i=1,3
  print*,M(i,:)
end do

END PROGRAM solvelinear

OUTPUT

M = 
1.0000000000000000        0.0000000000000000        0.0000000000000000     
0.0000000000000000        1.0000000000000000       -5.5511151231257827E-017
-4.4408920985006262E-016  -2.2204460492503131E-016   1.0000000000000000

By the way, your original code returns the expected unchanged value of A when compiled with gfortran 4.8.2

credondo
  • 744
  • 3
  • 10