Wrong Умножение матриц с 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 в качестве ответа, поскольку я (тривиально) умножаюсь на идентификатор перед вычислением ожидаемого значения RHO, как определено в подпрограмме EXPECTATION. Вместо этого я получаю 1.2732139384274929E-313. Полное отсутствие смысла.

Что может быть? В моем реальном коде, я хочу сделать сложное умножение матриц вида:

RHONEW = UCONJ*RHO*U, 

где UCONJ и U являются линейными комбинациями матриц.

Я не получаю правильный результат даже для тривиального умножения. Возможные ошибки, которые я могу исправить, чтобы перейти к моей проблеме?

Обратите внимание, что код находится во всех CAPS, так как я получил его от моего руководителя, который написал его в фиксированной форме. Код компилируется как с расширениями .f, так и с .f90.

Answer 1

Вы определяете SIGMAX внутри подпрограммы WRONGRESULT еще раз. Таким образом, он включен в подпрограмму (снова) и затеняет значение, определенное в модуле.

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

SIGMAX в этой подпрограмме никогда не инициализируется, поэтому matmul возвращает мусор.

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) 

все это не плохо, но это может быть проще рассмотреть вымышленными размеров массива ... Но тогда вам нужно выделить и DEALLOCATE TEMP и RHONEW, так что, вероятно, не какие-либо преимущества в этом 2х2 , Но пример приведен в случае, если вы закончите с большими массивами:

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