1

我试图在更大的代码中找到错误并将问题隔离到以下代码片段中，我将其放入单独的程序中进行故障排除。为什么矩阵“Ga”的一个条目（并且只有一个）会通过用零初始化 HEB 而改变？如果我删除这个初始化，两个“Ga”输出将是相同的。请向下滚动到底部以查看问题所在。非常感谢，这让我发疯了。我使用 ifort 进行编译。

program AC
 implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer :: n, ndim, k, j, i, l, m, o, steps
  real(dp) :: emax, omega, pi, EFermi, auev
  complex(dp) :: Grs,Gas, ACCond, tinyc, cunit, czero, cone

  complex(dp), allocatable :: GammaL(:,:)     
  complex(dp), allocatable :: GammaL_EB(:,:)  
  complex(dp), allocatable :: Gr(:,:)     
  complex(dp), allocatable :: Ga(:,:)  
  complex(dp), allocatable :: GrInv(:,:)    
  complex(dp), allocatable :: GaInv(:,:)  
  complex(dp), allocatable :: GrDiag(:,:)     
  complex(dp), allocatable :: GaDiag(:,:)  
  complex(dp), allocatable :: GrInvDiag(:,:)   
  complex(dp), allocatable :: GaInvDiag(:,:)  
  complex(dp), allocatable :: GammaR(:,:)    
  complex(dp), allocatable :: GammaR_EB(:,:)     
  complex(dp), allocatable :: R(:,:) 
  complex(dp), allocatable :: Yc(:,:)         
  complex(dp), allocatable :: Yd(:,:)        
  complex(dp), allocatable :: AnaInt(:,:)     
  complex(dp), allocatable :: H(:,:)         
  complex(dp), allocatable :: HEB(:,:)          
  complex(dp), allocatable :: HamEff(:,:)    
  complex(dp), allocatable :: EigVec(:,:)     
  complex(dp), allocatable :: InvEigVec(:,:)  
  complex(dp), allocatable :: EigVal(:)       
  complex(dp), allocatable :: ctemp(:,:)     
  complex(dp), allocatable :: S(:,:)          
  complex(dp), allocatable :: SelfL(:,:)     
  complex(dp), allocatable :: SelfR(:,:)      
  complex(dp), allocatable :: SHalf(:,:)     
  complex(dp), allocatable :: InvSHalf(:,:)   
  complex(dp), allocatable :: Abba(:,:)   
  complex(dp), allocatable :: Integrand(:,:)
  complex(dp), allocatable :: HDiag(:,:)
  complex(dp), allocatable :: unity(:,:)

!Lapack arrays and variables
  integer :: info, lwork
  complex(dp), allocatable :: work(:)       
  real(dp), allocatable :: rwork(:)    
  integer,allocatable :: ipiv(:)

!########################################################################

!Constants
    pi = 3.14159265359
    cunit = (0,1)
    czero = (0,0)
    cone = (1,0)
    tinyc = (0.0, 0.000000000001)


!System and calculation parameters
    ndim = 3 !Dimension of the Hamiltonian
    lwork = ndim*ndim
    EFermi = 0.0 
    emax = 5.0 !Energy in eV to which the admittance is to be calculated
    steps = 1000 !Number of energy steps for the admittance calculation


    allocate(Integrand(ndim,ndim))
    allocate(H(ndim,ndim))
    allocate(Yc(ndim,ndim))
    allocate(Yd(ndim,ndim))
    allocate(S(ndim,ndim))
    allocate(SelfL(ndim,ndim))
    allocate(SelfR(ndim,ndim))
    allocate(HamEff(ndim,ndim))
    allocate(GammaR(ndim,ndim))
    allocate(GammaL(ndim,ndim))
    allocate(AnaInt(ndim,ndim))
    allocate(EigVec(ndim,ndim))
    allocate(EigVal(ndim))
    allocate(InvEigVec(ndim,ndim))
    allocate(R(ndim,ndim))
    allocate(GammaL_EB(ndim,ndim))
    allocate(GammaR_EB(ndim,ndim))
    allocate(HEB(ndim,ndim))
    allocate(Ga(ndim,ndim))
    allocate(Gr(ndim,ndim))
    allocate(GaInv(ndim,ndim))
    allocate(GrInv(ndim,ndim))
    allocate(GaDiag(ndim,ndim))
    allocate(GrDiag(ndim,ndim))
    allocate(GaInvDiag(ndim,ndim))
    allocate(GrInvDiag(ndim,ndim))
    allocate(Abba(ndim,ndim))
    allocate(HDiag(ndim,ndim))
    allocate(unity(ndim,ndim))


!################################################


!Definition of the trial system

    H = reshape((/ -0.3, -0.8, -0.2, -0.8, -0.5, -0.14, -0.2, -0.14, -0.24 /), shape(H))
    SelfL = reshape((/ -0.1*cunit, -0.3*cunit, -0.0*cunit, -0.3*cunit, -0.0*cunit, -0.0*cunit, -0.0*cunit, -0.0*cunit, -0.0*cunit /), shape(SelfL))
    SelfR = reshape((/ -0.0*cunit, -0.0*cunit, -0.0*cunit, -0.0*cunit, -0.0*cunit, -0.6*cunit, -0.0*cunit, -0.6*cunit, -0.3*cunit /), shape(SelfR))

    HamEff(:,:)=(H(:,:) + SelfL(:,:) + SelfR(:,:))   
    GammaL(:,:)=(-2.0 * cmplx(aimag(SelfL(:,:)),0.0_dp,dp)) 
    GammaR(:,:)=(-2.0 * cmplx(aimag(SelfR(:,:)),0.0_dp,dp))

unity(:,:) = 0.0_dp
do i=1,ndim
unity(i,i) = cone
end do



!Calculate Eigenvalues and Eigenvectors of the Hamiltonian (Checked! Yields correct Eigenvalues for trial matrix).
!In the PDF: B = EigVec, B^(-1) = InvEigVec, Hk = EigVal

    allocate(ctemp(ndim,ndim))
    allocate(work(lwork),rwork(2*ndim))
    ctemp(:,:) = HamEff(:,:)
    call zgeev('N', 'V', ndim, ctemp, ndim, EigVal, InvEigVec, ndim, EigVec, ndim, work, lwork, rwork, 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)
    call zgetri(ndim,InvEigVec,ndim,ipiv,work,lwork,info)
    deallocate(work)
    deallocate(ipiv)


!Green Function at E=2
    GrInv(i,j) = 0.0_dp
    GaInv(i,j) = 0.0_dp
    Gr(i,j) = 0.0_dp
    Ga(i,j) = 0.0_dp

    do i=1,ndim
       do j=1,ndim
          GrInv(i,j) = (3.0_dp + tinyc)*unity(i,j) - HamEff(i,j)
       end do
    end do

    do i=1,ndim
       do j=1,ndim
          GaInv(i,j) = (3.0_dp - tinyc)*unity(i,j) - conjg(HamEff(j,i))
       end do
    end do


!Matrix Inversions

    Ga(:,:)=GaInv(:,:)
    lwork = 3*ndim
    allocate(ipiv(ndim))
    allocate(work(lwork))
    call zgetrf(ndim,ndim,Ga,ndim,ipiv,info)
    call zgetri(ndim,Ga,ndim,ipiv,work,lwork,info)
    deallocate(work)
    deallocate(ipiv)

    Gr(:,:)=GrInv(:,:)
    lwork = 3*ndim
    allocate(ipiv(ndim))
    allocate(work(lwork))
    call zgetrf(ndim,ndim,Gr,ndim,ipiv,info)
    call zgetri(ndim,Gr,ndim,ipiv,work,lwork,info)
    deallocate(work)
    deallocate(ipiv)




    GrInvDiag(i,j) = 0.0_dp
    GaInvDiag(i,j) = 0.0_dp
    GrDiag(i,j) = 0.0_dp
    GaDiag(i,j) = 0.0_dp
    do i=1,ndim
       do j=1,ndim
          GrInvDiag(i,j) = (3.0_dp + tinyc)*unity(i,j) - HDiag(i,j)
       end do
    end do

    do i=1,ndim
       do j=1,ndim
          GaInvDiag(i,j) = (3.0_dp - tinyc)*unity(i,j) - conjg(HDiag(j,i))
       end do
    end do



!Matrix Inversions

    GaDiag(:,:)=GaInvDiag(:,:)
    lwork = 3*ndim
    allocate(ipiv(ndim))
    allocate(work(lwork))
    call zgetrf(ndim,ndim,GaDiag,ndim,ipiv,info)
    call zgetri(ndim,GaDiag,ndim,ipiv,work,lwork,info)
    deallocate(work)
    deallocate(ipiv)

    GrDiag(:,:)=GrInvDiag(:,:)
    lwork = 3*ndim
    allocate(ipiv(ndim))
    allocate(work(lwork))
    call zgetrf(ndim,ndim,GrDiag,ndim,ipiv,info)
    call zgetri(ndim,GrDiag,ndim,ipiv,work,lwork,info)
    deallocate(work)
    deallocate(ipiv)



!The problem occurs here. Why does an entry of Ga get deleted by initialization of a new array HEB?

    write(*,*) Ga
    write(*,*) "--------------------"


     HEB(i,j) = 0.0_dp


    write(*,*) Ga




!#############################################

    deallocate(Integrand)
    deallocate(HEB)
    deallocate(Yc)
    deallocate(Yd)
    deallocate(HamEff)
    deallocate(GammaR)
    deallocate(GammaL)
    deallocate(Ga)
    deallocate(Gr)
    deallocate(GaInv)
    deallocate(GrInv)
    deallocate(GaDiag)
    deallocate(GrDiag)
    deallocate(GaInvDiag)
    deallocate(GrInvDiag)
    deallocate(AnaInt)
    deallocate(EigVec)
    deallocate(EigVal)
    deallocate(InvEigVec)
    deallocate(H)
    deallocate(S)
    deallocate(SelfL)
    deallocate(SelfR)
    deallocate(R)
    deallocate(GammaL_EB)
    deallocate(GammaR_EB)
    deallocate(Abba)
    deallocate(HDiag)
    deallocate(unity)
end program AC
4

1 回答 1

3

正如 ShinTakezou 指出的那样，您似乎在使用i和在循环j之外有问题。DO除非重新初始化，否则用于DO循环的变量比它们的最后一个值大 1（也就是说，i=j=4如果您要将其打印到屏幕上，您会发现）。显然，内存地址HEB(4,4)与GA(3,1).

解决此问题的方法是使用HEB=0._dpor HEB(:,:)=0._dp，您将在 print 语句之前和之后得到相同的答案。

于 2013-09-11T14:21:39.613 回答