2

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

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
4

1 回答 1

4

在调用 DGETRI 之前，您必须计算 LU 分解。DGETRF您正在使用例程的双重版本，因此必须使用(ZGETRF是复杂版本)计算 LU 。

以下代码编译并返回正确的值

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

于 2014-10-21T08:46:19.640 回答