間違ったマトリックス乗算（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)と書いています。 RHOの期待値を計算する前に、私が（些細なことながら）同一性を掛けているので、サブルーチンEXPECTATIONで定義されているように、私は答えとして1を取るべきです。代わりに、私は1.2732139384274929E-313.完全な意味がありません。

何が起こっている可能性がありますか？私の実際のコードでは、私は、フォームの複雑なマトリックス乗算をしたい：

UCONJと Uは行列の線形結合である

RHONEW = UCONJ*RHO*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) 

何も悪いありませんが、この2×2のいずれかの利点、想定配列サイズを使用することを検討する方が簡単な場合があり...しかし、その後、あなたが割り当てと割り当て解除TEMPとRHONEWする必要があるので、おそらくありません。しかし、大きな配列で終わる場合の例が提供されています：

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