在 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屋！