在 Fortran 中使用 ZGETRI 的逆矩阵错误

Question

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

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 编译的）。生成文件：

#$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 $<

我编译的时候没有错误。这是我的输出：

 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应该是身份

Answer 1

我建议您不要使用像 REAL(4) 或 COMPLEX(16) 这样带有文字数字的种类表示法。

首先，它丑陋且不便携。

其次，复杂变量可能会造成混淆。

此处您将变量定义为 COMPLEX(16)，但 ZGETRI 和所有其他 LAPACK Z 例程都需要 COMPLEX*16。这些不相同。

COMPLEX*16 是具有 REAL*8 个分量的复数的非标准表示法。 REAL*8 是 8 字节实数的非标准表示法，通常等同于 DOUBLE PRECISION.

COMPLEX(16) 是一个具有两个 REAL(16) 分量的复数，前提是存在这种复数。在提供 REAL(16) 的编译器中，这个实数是四倍精度，而不是双精度。

因此，您有效地传递了 32 字节的复杂变量，而不是 16 字节的复杂变量。

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

开始

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

和

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

在 Fortran 中使用 ZGETRI 的逆矩阵错误

Wrong inverse matrix using ZGETRI in Fortran

fortran

matrix

lapack

intel-fortran