前端之家

  1. 首页
  2. 问答

逆矩阵的结果错误

• 问答

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

我将测试代码粘贴到此处，以便您可以参考。

我在群集上运行带有外壳程序脚本的测试程序，并且该脚本也粘贴在此处。

PROGRAM TE
IMPLICIT NONE
INTEGER :: i,j,dime
DOUBLE COMPLEX,ALLOCATABLE :: a(:,:)
DOUBLE COMPLEX,ALLOCATABLE :: b(:,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 arguement 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 arguement 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,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,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,in_lda,in_work,in_lwork,in_info)

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

c = MATMUL(a,b)

WRITE (UNIT=3,FMT='(A1)') 'c'
DO i = 1,FMT=*) c(i,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,-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,2.00000000000000)
(1.00000000000000,3.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函数有问题吗？有人可以给我一些如何解决的建议吗？预先非常感谢。

panlikq 回答：逆矩阵的结果错误

（使用Mathematica的插图）

该问题与MKL无关。您的矩阵对于dime>2是奇异的。对于没有逆的矩阵，MKL不应给出正确的逆矩阵。

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

实际上，不难发现他们在dime>1的排名中排名第二。

In[3]:= col[n_]:=Table[i+I j,2}]
In[4]:= row[n_]:={Table[2-j,n}],Table[j-1,n}]}
In[5]:= Table[Norm[mat[n]-col[n].row[n]],2,20}]
Out[5]= {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]
,

您能用最新的MKL 2019版本检查情况吗？我记得MKL在getri版本2018中存在一些问题。

fortranintel-mklmatrix-inverse
本文链接：https://www.f2er.com/3137001.html

大家都在问