LAPACK рутина инверсии странно смешивает все переменные

Question

Я программирование в Фортране и пытаюсь использовать матричный преобразователь DGETRI из пакета LaPack:

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

Но очень странно, что, кажется, возиться со всеми своими переменные. В этом очень простом примере моя матрица А, инициализированная в начале программы, изменяется с использованием DGETRI, хотя DGETRI не включает A ...

Может кто-нибудь сказать мне, что происходит? Благодаря!

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 

Это выход:

A = 
    1.0000000000000000  0.0000000000000000 
    -1.0000000000000000  0.33333333333333331 

Answer 1

Вы должны вычислить разложение LU перед вызовом DGETRI. Вы используете двойную версию подпрограмм, поэтому LU должен быть вычислен с помощью DGETRF (ZGETRF - сложная версия).

Следующий код компилирует и возвращает правильные значения

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 

ВЫВОД

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

Кстати, ваш исходный код возвращает ожидаемое неизменное значение А при компиляции с gfortran 4.8.2