Write-Statement depends on earlier Write-Statements?

Question

Excuse me that the code I'm posting is no minimal example. When I try to further reduce the code, the situation I want to show you breaks down. I'm a beginner and simply cannot understand what is happening here.

Please scroll down to the Write-Statement "WTF".

Effect Number 1: The output of the Write statement that gives out "ctemp" to stdout is dependent on the earlier write statement giving out "WTF". Run the code, then comment out the "WTF" part and run again. How can this be?

Effect Number 2: The ctemp that is supposed to be written to stdout is last defined in the matmul-calculation. When I instead overwrite this result with a matrix full of 1's (as currently commented out), the output does NOT longer depend on an earlier WTF-Write-Statement.

I'm at a loss and cannot see any logic in this. What is going on? Thanks.

EDIT: as requested, I specify the differing output I get.

WITH Write-Statement:

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

WITHOUT Write-Statement:

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

THE COMPILER:

Intel(R) Fortran Intel(R) 64 Compiler XE for applications running on Intel(R) 64, Version 12.0.3.174 Build 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

I think I may have found the answer:

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

ctemp is a 2x2 (complex) matrix, but you have only defined the diagonals. If you initialize it to zero (or set ctemp(1,2)=0.0; ctemp(2,1)=0.0) you will get the same answer both times (and on both compilers):

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