逆矩阵的错误结果

Question

我想通过 Fortran 代码计算复杂矩阵的逆矩阵以及矩阵 52 乘 52 的大小。我使用 MKL 库中的 ZGETRF 和 ZGETRI 函数来运行测试。如果复矩阵的大小是 2×2，这两个函数效果很好，因为原矩阵和逆矩阵的乘积是单位矩阵；但是，我发现如果复数矩阵的大小是 3 x 3，则结果是错误的，因为原矩阵和逆矩阵的乘法不再是单位矩阵。

我把我的测试代码贴在这里，你可以参考一下。

我在集群上用shell脚本运行测试程序，脚本也粘贴在这里。

PROGRAM TE
IMPLICIT NONE
INTEGER :: i, j, dime
DOUBLE COMPLEX, ALLOCATABLE :: a(:,:)
DOUBLE COMPLEX, ALLOCATABLE :: b(:,:)
DOUBLE COMPLEX, ALLOCATABLE :: c(:,:)

!Variables used in LU decomposition subroutine (ZGETRF)
INTEGER                       :: nu_r                 !Number of row in target matrix
INTEGER                       :: nu_c                 !Number of column in target matrix
INTEGER                       :: lu_lda               !Leading dimension of target matrix
INTEGER, ALLOCATABLE          :: lu_ipiv(:)           !Array with dimension of .GE. MIN(nu_r,nu_c)
INTEGER                       :: lu_info              !Judgement (0-successful;<0-error and infor = number, then
                                                      !'1-number'th argument has illegal value;>0 and info = 
                                                      !number - LU matrix(number,number) is zero)
!

!Variables used in inverse subroutine (ZGETRI)
INTEGER                       :: nu_o                 !Order of LU matrix
DOUBLE COMPLEX, ALLOCATABLE   :: in_work(:)           !Workspace array with dimension of .GE. MAX(1,in_lwork)
INTEGER                       :: in_lda               !Leading dimension of LU decomposed matrix .GE. MAX(1,nu_o)
INTEGER                       :: in_lwork             !Size of in_work array with value .GE. nu_o
INTEGER                       :: in_info              !Judgement (0-successful;<0-error and infor = number, then
                                                      !'1-number'th argument has illegal value;>0 and info =
                                                      !number; then, LU matrix(number,number) is zero and inversion
                                                      !cannot be accomplished)

dime = 3
ALLOCATE (a(dime,dime))
ALLOCATE (b(dime,dime))
ALLOCATE (c(dime,dime))

DO i = 1, dime, 1
   DO j = 1, dime, 1
      a(i,j) = CMPLX(DBLE(i), DBLE(j))
   END DO
END DO

b = a

OPEN (UNIT=3, FILE='te_in.dat', STATUS='UNKNOWN')

WRITE (UNIT=3, FMT='(A1)') 'a'
DO i = 1, dime, 1
   DO j = 1, dime, 1
      WRITE (UNIT=3, FMT=*) a(i,j)
   END DO
END DO
WRITE (UNIT=3, FMT=*)

!Initialising parameters for LU decomposition subroutine (ZGETRF)
nu_r = dime
nu_c = dime
lu_lda = dime
ALLOCATE (lu_ipiv(dime))
lu_ipiv = 0
!

!Initialising parameters for inverse subroutine (ZGETRI)
nu_o = dime
in_lwork = dime
ALLOCATE (in_work(in_lwork))
in_work = (0.0d0, 0.0d0)
in_lda = dime
!

CALL ZGETRF(nu_r,nu_c,b,lu_lda,lu_ipiv,lu_info)
CALL ZGETRI(nu_o,b,in_lda,lu_ipiv,in_work,in_lwork,in_info)

WRITE (UNIT=3, FMT='(A1)') 'b'
DO i = 1, dime, 1
   DO j = 1, dime, 1
      WRITE (UNIT=3, FMT=*) b(i,j)
   END DO
END DO
WRITE (UNIT=3, FMT=*)

c = MATMUL(a,b)

WRITE (UNIT=3, FMT='(A1)') 'c'
DO i = 1, dime, 1
   DO j = 1, dime, 1
      WRITE (UNIT=3, FMT=*) c(i,j)
   END DO
END DO
WRITE (UNIT=3, FMT=*)

DEALLOCATE (lu_ipiv)
DEALLOCATE (a)
DEALLOCATE (b)
DEALLOCATE (c)

CLOSE (UNIT=3)
STOP
END PROGRAM TE

下面是输出文件2×2复数矩阵的内容，正确。

a
(1.00000000000000,1.00000000000000)
(1.00000000000000,2.00000000000000)
(2.00000000000000,1.00000000000000)
(2.00000000000000,2.00000000000000)

b
(-2.00000000000000,2.00000000000000)
(2.00000000000000,-1.00000000000000)
(1.00000000000000,-2.00000000000000)
(-1.00000000000000,1.00000000000000)

c
(1.00000000000000,-2.220446049250313E-016)
(0.000000000000000E+000,2.220446049250313E-016)
(0.000000000000000E+000,4.440892098500626E-016)
(1.00000000000000,0.000000000000000E+000)

但是，程序不会为 3 x 3 复数矩阵生成正确的输出内容，如下所示。

a
(1.00000000000000,1.00000000000000)
(1.00000000000000,2.00000000000000)
(1.00000000000000,3.00000000000000)
(2.00000000000000,1.00000000000000)
(2.00000000000000,2.00000000000000)
(2.00000000000000,3.00000000000000)
(3.00000000000000,1.00000000000000)
(3.00000000000000,2.00000000000000)
(3.00000000000000,3.00000000000000)

b
(723205779577744.,262983919846453.)
(-1.446411559155488E+015,-525967839692903.)
(723205779577744.,262983919846451.)
(-1.446411559155489E+015,-525967839692904.)
(2.892823118310976E+015,1.051935679385808E+015)
(-1.446411559155487E+015,-525967839692903.)
(723205779577745.,262983919846452.)
(-1.446411559155488E+015,-525967839692905.)
(723205779577743.,262983919846452.)

c
(0.437500000000000,0.000000000000000E+000)
(1.25000000000000,0.000000000000000E+000)
(-0.250000000000000,0.000000000000000E+000)
(0.125000000000000,0.000000000000000E+000)
(0.500000000000000,0.000000000000000E+000)
(0.750000000000000,0.000000000000000E+000)
(-0.500000000000000,0.000000000000000E+000)
(0.000000000000000E+000,-1.00000000000000)
(1.50000000000000,0.500000000000000)

我还用其他更高阶的复杂矩阵测试我的代码，结果都是错误的。我在下面粘贴了我用来提交作业以运行此测试程序的 shell 脚本。

#!/bin/bash
module switch PrgEnv-cray PrgEnv-intel
#module load intel
#----------------------------------------------------------#
ifort -O3 te.f90 \
  -Wl,--start-group${MKLROOT}/lib/intel64/libmkl_intel_lp64.a \
                   ${MKLROOT}/lib/intel64/libmkl_core.a \
                   ${MKLROOT}/lib/intel64/libmkl_sequential.a \
  -Wl,--end-group-lpthread -lm –ldl

谁能告诉我我的代码有什么问题？ ZGETRF 和 ZGETRI 函数有问题吗？谁能给我一些关于如何整理它的建议？

Answer 1

（使用 Mathematica 的插图）

该问题与 MKL 无关。对于dime>2，您的矩阵是单数的。 MKL 不应该为没有逆矩阵的矩阵提供正确的逆矩阵。

In[1]:= mat[n_]:=Table[i+I j,{i,1,n},{j,1,n}]
In[2]:= Table[Det[mat[n]],{n,1,20}]
Out[2]= {1+I,-I,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}

其实不难看出他们在dime>1排名第二。

In[3]:= col[n_]:=Table[i+I j,{i,1,n},{j,1,2}]
In[4]:= row[n_]:={Table[2-j,{j,1,n}],Table[j-1,{j,1,n}]}
In[5]:= Table[Norm[mat[n]-col[n].row[n]],{n,2,20}]
Out[5]= {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}

也就是说，您的矩阵可能的rank factorization 是 C*R ：

C=[1+I 1+2I]
  [2+I 2+2I]
  [...  ...]
  [n+I n+2I]

R=[1 0 -1 ... 2-n]
  [0 1  2 ... n-1]

Answer 2

你能用最新的 MKL 2019 版本检查这个案例吗？我记得很清楚，MKL 在 getri 进入 2018 版时遇到了一些问题。