在Fortran中使用ZGETRI的错误逆矩阵 [英] Wrong inverse matrix using ZGETRI in Fortran

问题描述

我正在尝试使用ZGETRI计算复杂矩阵的逆，但是 即使执行没有错误(info = 0)， 它没有给我正确的逆矩阵，我绝对有 不知道错误是从哪里来的.

I am trying to compute the inverse of a complex matrix with ZGETRI, but even if it executes without error (info = 0), it does not give me the correct inverse matrix and I have absolutely no idea where the error comes from.

PROGRAM solvelinear
implicit none
INTEGER                        :: i,j,info,lwork
INTEGER,dimension(3)        :: ipiv
COMPLEX(16), dimension(3,3) :: C,Cinv,M,LU
COMPLEX(16),allocatable :: work(:)

info=0
lwork=100
allocate(work(lwork))
ipiv=0
work=(0.d0,0.d0)

C(1,1)=(0.d0,-1.d0)
C(1,2)=(1.d0,5.d0)
C(1,3)=(2.d0,-2.d0)
C(2,1)=(4.d0,-1.d0)
C(2,2)=(2.d0,-3.d0)
C(2,3)=(-1.d0,2.d0)
C(3,1)=(1.d0,0.d0)
C(3,2)=(3.d0,-2.d0)
C(3,3)=(0.d0,1.d0)

write(*,*)"C = "
do i=1,3
   write(*,10)(C(i,j),j=1,3)
end do

!-- LU factorisation
LU=C
CALL ZGETRF(3,3,LU,3,ipiv,info)
write(*,*)'info = ',info
write(*,*)"LU = "
do i=1,3
   write(*,10)(LU(i,j),j=1,3)
end do

!-- Inversion of matrix C using the LU

Cinv=LU
CALL ZGETRI(3,Cinv,3,ipiv,work,lwork,info)
write(*,*)'info = ',info
write(*,*)"Cinv = "
do i=1,3
   write(*,10)(Cinv(i,j),j=1,3)
end do

!-- computation of C^-1 * C to check the inverse
M = matmul(Cinv,C)
write(*,*)"M = "
do i=1,3
   write(*,10)(M(i,j),j=1,3)
end do
          10 FORMAT(3('(',F20.10,',',F20.10,') '))

END PROGRAM solvelinear

我使用ifort进行编译(而我的LAPACK图书馆版本3.7.1也使用ifort进行编译).生成文件:

I compile with ifort (and my LAPACK librairies version 3.7.1 are also compiled with ifort). Makefile:

#$Id: Makefile $
.SUFFIXES: .f90 .f .c .o
FC = ifort
FFLAGS = -g -check all -zmuldefs -i8
LIBS = -L/path/to/lapack-3.7.1 -llapack -ltmglib -lrefblas
MAIN = prog.o
EXEC = xx
all:  ${MAIN} Makefile
    ${FC} ${FFLAGS} -o ${EXEC} ${MAIN} ${LIBS}
.f.o: ${MODS} Makefile
    ${FC} ${FFLAGS} -c $<
.f90.o: ${MODS} Makefile
    ${FC} ${FFLAGS} -c $<

编译时没有错误.这是我的输出:

I have no errors when compiling. Here is my output:

 C = 
(       0.00000,      -1.00000) (       1.00000,       5.00000) (       2.00000,      -2.00000) 
(       4.00000,      -1.00000) (       2.00000,      -3.00000) (      -1.00000,       2.00000) 
(       1.00000,       0.00000) (       3.00000,      -2.00000) (       0.00000,       1.00000) 
 info =                      0
 LU = 
(       4.00000,       0.00000) (       2.00000,  120470.58824) (       2.00000,      -2.00000) 
(       0.00000,       0.00000) (28003147.29412,      -3.00000) (      -1.00000,       2.00000) 
(       1.00000,       0.00000) (       3.00000,      -2.00000) (       0.00000,       1.00000) 
 info =                      0
 Cinv = 
(       0.00000,       0.00000) (      -0.00000,      -0.00000) (       2.00000,      -2.00000) 
(      -0.00000,       0.00000) (      -0.00000,      -3.00000) (      -1.00000,       2.00000) 
(      -0.00000,      -0.00000) (       3.00000,      -2.00000) (       0.00000,       1.00000) 
 M = 
(       2.00000,      -2.00000) (       2.00000,     -10.00000) (       2.00000,       2.00000) 
(      -4.00000,     -10.00000) (      -8.00000,       2.00000) (       4.00000,       2.00000) 
(      10.00000,     -10.00000) (       2.00000,     -10.00000) (      -0.00000,       8.00000) 

如果我没记错的话，M应该是身份.

And M should be the identity if I'm not wrong.

推荐答案

我建议您不要对诸如REAL(4)或COMPLEX(16)这样的文字数字使用种类符号.

I suggest you to NOT use the kind notation with literal numbers like REAL(4) or COMPLEX(16).

首先，它丑陋且不便于携带.

First, it is ugly and not portable.

第二，对于复杂的变量可能会造成混淆.

Second, it can be confusing for complex variables.

在这里，将变量定义为COMPLEX(16)，但是将ZGETRI和所有其他LAPACK Z例程定义为COMPLEX*16.这些不相同.

Here you define your variables as COMPLEX(16), but ZGETRI, and all other LAPACK Z routines, expects COMPLEX*16. These are NOT the same.

COMPLEX*16是具有REAL*8组件的复数的非标准符号. REAL*8是8字节实数的非标准符号，通常等同于DOUBLE PRECISION.

COMPLEX*16 is a non-standard notation for complex numbers with REAL*8 components. REAL*8 is a nonstandard notation for 8 byte real numbers that are normally equivalent to DOUBLE PRECISION.

COMPLEX(16)是具有两个REAL(16)成分的复数，只要存在这种种类.在提供REAL(16)的编译器中，此实数是四倍精度，而不是双精度.

COMPLEX(16) is a complex number with two REAL(16) components, provided such a kind exists. In compilers which provide REAL(16) this real is a quadruple precision, not double precision.

因此，您实际上是在传递32字节的复杂变量，而不是16字节的复杂变量.

So you are effectively passing 32-byte complex variables instead of 16-byte complex variables.

有足够的资源可以学习如何正确使用Fortran类型.您可以以

There are enough resources where to learn how to use Fortran kinds properly. You can start with

integer, parameter :: dp = kind(1.d0)

和

real(dp) :: x
complex(dp) :: c

这篇关于在Fortran中使用ZGETRI的错误逆矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！