LAPACK 反演例程奇怪地混合了所有变量

我正在使用 Fortran 进行编程，并尝试使用 Lapack 包中的 DGETRI 矩阵逆变器：

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

但非常奇怪的是，它似乎弄乱了我所有的变量。在这个非常简单的例子中，我的矩阵 A 在程序开始时初始化，随着 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

在调用 DGETRI 之前，您必须计算 LU 分解。 您正在使用例程的双版本，因此 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

OUTPUT

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 编译时，您的原始代码返回 A 的预期不变值