LAPACK反演例程奇怪地混合了所有变量 [英] LAPACK inversion routine strangely mixes up all variables

问题描述

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

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

但是非常奇怪的是，它似乎弄乱了我所有的变量.在这个非常简单的示例中，即使应用程序DGETRI不涉及A…，也会在应用程序开始时初始化的矩阵A随DGETRI的变化而变化.

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

这是输出:

 A =
   1.0000000000000000        0.0000000000000000
  -1.0000000000000000       0.33333333333333331

推荐答案

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

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

输出

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不变值

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

这篇关于LAPACK反演例程奇怪地混合了所有变量的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！