写语句取决于早期的写语句？

Question

不好意思，我发布的代码不是最小的例子。当我试图进一步减少代码时，我想向您展示的情况就会崩溃。我是个初学者，根本无法理解这里发生了什么。

请向下滚动到写声明"WTF“。

效果1:将"ctemp“发送给stdout的写语句的输出取决于先前发出"WTF”的写语句。运行代码，然后注释掉"WTF“部分，然后再次运行。这怎么可能呢？

效果2:应该写入stdout的ctemp是在matmul计算中最后定义的。当我用一个满是1的矩阵覆盖这个结果时(就像目前注释掉的那样)，输出不再依赖于早期的WTF写语句。

我不知所措，看不出其中有什么逻辑。怎么一回事？谢谢。

编辑:根据要求，我指定了所获得的不同输出。

带有写语句的：

0.40000000E+01 0.00000000E+00 0.20000000E+01 0.10000000E+01 0.10000000E+01 0.20000000E+01 0.20000000E+01 0.30000000E+01

没有写语句的：

0.22583338E+01 -0.17920885E+01 0.13104573E+01 -0.21149418E+01 0.28983440E+01 0.24774309E+01 0.37416662E+01 0.47920885E+01

编译器：

Intel(R) Fortran Intel(R) 64 Compiler XE用于运行在Intel(R) 64上的应用程序，版本12.0.3.174构建20110309

program testlapack
 implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer :: n, ndim, k, j, i
  complex(dp) :: cunit, czero, cone

  complex(dp), allocatable :: H(:,:)          !Input Hamiltonian
  complex(dp), allocatable :: EigVec(:,:)     !Eigenvector matrix
  complex(dp), allocatable :: InvEigVec(:,:)  !Inverted eigenvector matrix
  complex(dp), allocatable :: EigVal(:)       !Eigenvalue vector
  complex(dp), allocatable :: ctemp(:,:)      !Temporary array
  complex(dp), allocatable :: ctemp2(:,:)      !Temporary array


!Lapack arrays and variables
  integer :: info, lwork
  complex(dp), allocatable :: work(:)       
  real(dp), allocatable :: rwork(:)    
  integer,allocatable :: ipiv(:)
  
  ndim=2
  lwork=ndim*ndim
    allocate(H(ndim,ndim))
    allocate(EigVec(ndim,ndim))
    allocate(EigVal(ndim))
    allocate(InvEigVec(ndim,ndim))
    allocate(ctemp2(ndim,ndim))


    H = reshape((/ (4,0), (1,2), (2,1), (2,3) /), shape(H))

  
   
   
     allocate(ctemp(ndim,ndim))
     ctemp(:,:) = H(:,:)
     allocate(work(lwork),rwork(2*ndim))
     call zgeev('N', 'V', ndim, ctemp, ndim, EigVal, InvEigVec, ndim, EigVec, ndim, work, lwork, rwork, info)
     if(info/=0)write(*,*) "Warning: zgeev info=", info
     deallocate(work,rwork)
     deallocate(ctemp) 
   
  
     InvEigVec(:,:)=EigVec(:,:)
     lwork = 3*ndim
     allocate(ipiv(ndim))
     allocate(work(lwork))
     call zgetrf(ndim,ndim,InvEigVec,ndim,ipiv,info)
     if(info/=0)write(*,*) "Warning: zgetrf info=", info   ! LU decomposition
     call zgetri(ndim,InvEigVec,ndim,ipiv,work,lwork,info)
     if(info/=0)write(*,*) "Warning: zgetri info=", info ! Inversion by LU decomposition (Building of InvEigVec)
     deallocate(work)
     deallocate(ipiv)



write(*,*) "WTF"
     
     allocate(ctemp(ndim,ndim))
     do i=1,ndim
        ctemp(i,i) = EigVal(i)
     end do
     ctemp2 = matmul(ctemp, InvEigVec)
     ctemp = matmul(EigVec,ctemp2)  

 !    ctemp = reshape((/ (1,1), (1,1), (1,1), (1,1) /), shape(H))

   


         do i=1, ndim
             do j=1, ndim
                write(*, '(2e17.8)', advance='NO') real(ctemp(i,j)), aimag(ctemp(i,j))
             end do
             Write(128,*)
         end do


   deallocate(H)
   deallocate(EigVal)
   deallocate(EigVec)
   deallocate(InvEigVec)
   deallocate(ctemp)
   deallocate(ctemp2)
  

end program testlapack

Answer 1

我想我可能找到了答案：

allocate(ctemp(ndim,ndim))
do i=1,ndim
   ctemp(i,i) = EigVal(i)
end do
ctemp2 = matmul(ctemp, InvEigVec)
ctemp = matmul(EigVec,ctemp2)

ctemp是一个2x2 (复)矩阵，但您只定义了对角线。如果将其初始化为零(或设置ctemp(1,2)=0.0; ctemp(2,1)=0.0)，则两次(以及在两个编译器上)都会得到相同的答案：

0.40000000E+01   0.66613381E-15
0.20000000E+01   0.10000000E+01
0.10000000E+01   0.20000000E+01
0.20000000E+01   0.30000000E+01