CODEWITHSUNDEEP
fortranlapacklinear-equation
Vim154
Vim154

Reputation: 23

LAPACK inversion routine strangely mixes up all variables

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

Upvotes: 2

Views: 863

Answers (1)

credondo
credondo

Reputation: 744

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

Upvotes: 4

Related Questions