0

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

请向下滚动到写语句“WTF”。

效果号 1：向 stdout 发出“ctemp”的 Write 语句的输出取决于之前发出“WTF”的 write 语句。运行代码，然后注释掉“WTF”部分并再次运行。怎么会这样？

效果号 2：应该写入标准输出的 ctemp 是在 matmul 计算中最后定义的。当我用一个全为 1 的矩阵（如当前注释掉的那样）覆盖这个结果时，输出不再依赖于早期的 WTF-Write-Statement。

我很茫然，看不到任何逻辑。到底是怎么回事？谢谢。

编辑：根据要求，我指定了我得到的不同输出。

WITH 写语句：

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 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
4

1 回答 1

3

我想我可能已经找到了答案：

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（复数）矩阵，但您只定义了对角线。如果将其初始化为零（或 set 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
于 2013-08-30T13:39:46.820 回答