MATMUL (Fortran) 的错误矩阵乘法

Question

我在 minimalexample.f90:

中保存了以下最小示例

        MODULE FUNCTION_CONTAINER
        IMPLICIT NONE
        SAVE

        INTEGER, PARAMETER :: DBL = SELECTED_REAL_KIND(P = 15,R = 200)

        INTEGER :: DIMSYS, DIMMAT

        COMPLEX(KIND = DBL), DIMENSION(4,1) :: INSTATE_BASISSTATES

        COMPLEX(KIND = DBL), DIMENSION(2,2) :: SIGMAX

        REAL(KIND = DBL), DIMENSION(2,2) :: BASISSTATES

        REAL(KIND = DBL), DIMENSION(2,2) :: ID

        COMPLEX(KIND = DBL),DIMENSION(2,2) :: PROJECTOR

        CONTAINS

        SUBROUTINE INDEXCONVERTER(N,K,L)
        IMPLICIT NONE
        INTEGER, INTENT(IN)::N
        INTEGER, INTENT(OUT)::K,L
        INTEGER::X, REMAINDER
        X = N/DIMSYS
        REMAINDER = N - (X * DIMSYS)
        IF (REMAINDER == 0) THEN
        K = X
        L = DIMSYS
        ELSE
        K = X + 1
        L = REMAINDER
        END IF
        END SUBROUTINE INDEXCONVERTER


        SUBROUTINE DENSITYMATRIX(X,RHO)
        IMPLICIT NONE
        COMPLEX(KIND = DBL), DIMENSION(DIMMAT,1), INTENT(IN) :: X
        COMPLEX(KIND = DBL), DIMENSION(2,2),INTENT(OUT)::RHO
        INTEGER :: J, K, L
        DO J = 1, DIMMAT
        CALL INDEXCONVERTER(J,K,L)
        RHO(K,L) = X(J,1)
        END DO
        END SUBROUTINE DENSITYMATRIX

        SUBROUTINE WRONGRESULT(X,RHONEW)
        IMPLICIT NONE
        COMPLEX(KIND = DBL), DIMENSION(DIMMAT,1), INTENT(IN) :: X
        COMPLEX(KIND = DBL),DIMENSION(DIMSYS,DIMSYS),INTENT(OUT)::RHONEW
        COMPLEX(KIND = DBL), DIMENSION(DIMSYS,DIMSYS) :: RHO
        REAL(KIND = DBL), DIMENSION(DIMSYS,DIMSYS) :: ID
        COMPLEX(KIND = DBL), DIMENSION(DIMSYS,DIMSYS) :: SIGMAZ
        COMPLEX(KIND = DBL), DIMENSION(DIMSYS,DIMSYS) :: SIGMAX
        CALL DENSITYMATRIX(X,RHO)
        RHONEW = matmul(SIGMAX,rho)
        END SUBROUTINE WRONGRESULT

        SUBROUTINE EXPECTATION(X,D,ANS)
        IMPLICIT NONE
        COMPLEX(KIND = DBL), DIMENSION(DIMMAT,1), INTENT(IN) :: X
        COMPLEX(KIND = DBL), DIMENSION(DIMSYS,DIMSYS) :: RHONEW
        COMPLEX(KIND = DBL), DIMENSION(DIMSYS,DIMSYS), INTENT(IN) :: D
        REAL(KIND = DBL),INTENT(OUT)::ANS
        COMPLEX(KIND = DBL),DIMENSION(DIMSYS, DIMSYS) :: TEMP
        INTEGER :: J
        REAL(KIND = DBL)::SUMM
        SUMM = 0.0D0
        CALL WRONGRESULT(X,RHONEW)
        TEMP = MATMUL(D,RHONEW)
        DO J = 1, DIMSYS
        SUMM = SUMM + DREAL(TEMP(J,J))
        END DO
        ANS = SUMM
        END SUBROUTINE EXPECTATION

        SUBROUTINE RK(ANSWER)
        IMPLICIT NONE
        REAL(KIND = DBL), INTENT(OUT) :: ANSWER
        COMPLEX(KIND = DBL), DIMENSION(DIMMAT,1)::X
        REAL(KIND = DBL) :: X_EXPECTATION
        REAL(KIND = DBL)::T
        T = 0.0D0
        X = INSTATE_BASISSTATES
        CALL EXPECTATION(X,SIGMAX,X_EXPECTATION)
        ANSWER = X_EXPECTATION
        END SUBROUTINE RK

        END MODULE FUNCTION_CONTAINER

        PROGRAM ONE
        USE FUNCTION_CONTAINER

        IMPLICIT NONE

        REAL(KIND = DBL) :: ANS

        SIGMAX(1,1) = (0.0D0,0.0D0)
        SIGMAX(1,2) = (1.0D0,0.0D0)
        SIGMAX(2,1) = (1.0D0,0.0D0)
        SIGMAX(2,2) = (0.0D0,0.0D0)

        ID(1,1) = 1.0D0
        ID(1,2) = 0.0D0
        ID(2,1) = 0.0D0
        ID(2,2) = 1.0D0

        DIMSYS = 2
        DIMMAT = 4

        BASISSTATES = ID

        INSTATE_BASISSTATES(1,1) = (0.5D0,0.0D0)
        INSTATE_BASISSTATES(2,1) = (0.5D0,0.0D0)
        INSTATE_BASISSTATES(3,1) = (0.5D0,0.0D0)
        INSTATE_BASISSTATES(4,1) = (0.5D0,0.0D0)

        CALL RK(ANS)

        WRITE (*,*) ANS



        END PROGRAM ONE

当我运行它时，编译器 cygwin 将答案打印为 1。没关系 - 正如预期的那样。现在在子程序中出错了，我为 RHONEW. 使用不同的表达式，例如，使用 RHONEW = 2*RHO，我得到 2 作为答案。再次，正如预期的那样。

现在，我写RHONEW = matmul(id,rho)。我应该得到 1 作为答案，因为我在计算子例程 EXPECTATION 中定义的 RHO 的期望值之前（平凡地）乘以恒等式。相反，我得到 1.2732139384274929E-313. 完全是胡说八道。

这是怎么回事？在我的实际代码中，我想执行以下形式的复杂矩阵乘法：

 RHONEW = UCONJ*RHO*U,

其中 UCONJ 和 U 是矩阵的线性组合。

即使是微不足道的乘法，我也没有得到正确的结果。我可以纠正哪些可能的错误以继续解决我的问题？

请注意，代码全部大写，因为我是从我的主管那里得到的，他以固定形式编写了代码。代码使用 .f 和 .f90 扩展进行编译。

Answer 1

您再次在子例程 WRONGRESULT 中定义 SIGMAX。因此，它（再次）在子例程中限定范围，并隐藏在模块中定义的子例程。

    SUBROUTINE WRONGRESULT(X,RHONEW)
      ! ...
      COMPLEX(KIND = DBL), DIMENSION(DIMSYS,DIMSYS) :: SIGMAX
      CALL DENSITYMATRIX(X,RHO)
      RHONEW = matmul(SIGMAX,rho)
      ! ...

这个子程序中的SIGMAX从未被初始化，所以matmulreturns垃圾。

Answer 2

有时会有像“-check uninit”这样的编译器开关。使用这些通常是值得的，以及“-warn unused”和“-check bounds”。然后，一旦整个工作正常并干净，您就可以删除这些检查并放入“-O2”。一些缩进会在美学上令人愉悦。对于更大的阵列，你可以初始化整个东西...... 我认为复杂的初始化值得声明为 CMPLX，但也许 gfortran 不需要那样做？？示例：

    !... Indent
    DO J = 1, DIMMAT
      CALL INDEXCONVERTER(J,K,L)
      RHO(K,L) = X(J,1)
    END DO
    !... CMPLX ... array init
    SIGMAX(:,:) = CMPLX(1.0D0,0.0D0)
    SIGMAX(1,1) = CMPLX(0.0D0,0.0D0)
    SIGMAX(2,2) = CMPLX(0.0D0,0.0D0)
    !... Whole array init
    ID(:,:) = 0.0D0
    !... diag init
    DO I = 1, 2
      ID(I,I) = 1.0D0
    ENDDO
    !... Whole array init
    INSTATE_BASISSTATES(:,1) = CMPLX(0.5D0,0.0D0)

none 这很糟糕，但考虑使用假定的数组大小可能更容易......但是你需要分配和解除分配 TEMP 和 RHONEW，所以在这个 2x2 中可能没有任何好处。但是提供了一个示例，以防您最终得到更大的数组：

    SUBROUTINE EXPECTATION(X, D, ANS)
    IMPLICIT NONE
    COMPLEX(KIND = DBL), DIMENSION(:,1), INTENT(IN   ) :: X
    COMPLEX(KIND = DBL), DIMENSION(:,:), INTENT(IN   ) :: D
    REAL(KIND = DBL)                   , INTENT(  OUT) :: ANS
    !~~Local~~
    COMPLEX(KIND = DBL), DIMENSION(DIMSYS,DIMSYS)      :: RHONEW
    COMPLEX(KIND = DBL), DIMENSION(DIMSYS,DIMSYS)      :: TEMP
    INTEGER                                            :: J
    !REAL(KIND = DBL)                                   :: SUMM

    CALL WRONGRESULT(X,RHONEW)
    TEMP = MATMUL(D,RHONEW)
    !SUMM = 0.0D0
    DO J = LBOUND(TEMP,1), UBOUND(TEMP,1)  !DIMSYS
      !SUMM = SUMM + DREAL(TEMP(J,J))
      ANS = ANS + DREAL(TEMP(J,J))
    END DO
    !ANS = SUMM
    END SUBROUTINE EXPECTATION

MATMUL (Fortran) 的错误矩阵乘法

Wrong Matrix Multiplication with MATMUL (Fortran)

fortran

gfortran

matrix-multiplication