subroutine mpi_lanczos_eigh_d(MpiComm,MatVec,Egs,Vect,Nitermax,iverbose,threshold,ncheck,vrandom)
integer :: MpiComm
interface
subroutine MatVec(Nloc,vin,vout)
integer :: Nloc
real(8) :: vin(Nloc)
real(8) :: vout(Nloc)
end subroutine MatVec
end interface
real(8) :: egs
real(8),dimension(:) :: vect !Nloc
integer :: Nitermax
!
integer :: Nloc
real(8),optional :: threshold
integer,optional :: ncheck
logical,optional :: iverbose
logical,optional :: vrandom
!
real(8),dimension(size(vect)) :: vin,vout
integer :: iter,nlanc
real(8),dimension(Nitermax+1) :: alanc,blanc
real(8),dimension(Nitermax,Nitermax) :: Z
real(8),dimension(Nitermax) :: diag,subdiag,esave
real(8) :: a_,b_,norm,diff,norm_tmp
integer :: i,ierr
logical :: vran=.true.
!
logical :: mpi_master
!
if(MpiComm==Mpi_Comm_Null)return
!
mpi_master=get_master_MPI(MpiComm)
!
Nloc = size(vect)
!
if(present(iverbose))verb=iverbose
if(present(threshold))threshold_=threshold
if(present(ncheck))ncheck_=ncheck
if(present(vrandom))vran=vrandom
!
norm_tmp=dot_product(vect,vect); norm=0d0
call AllReduce_MPI(MpiComm,norm_tmp,norm)
!
if(norm==0d0)then
if(vran)then
! call random_seed(size=nrandom)
! if(allocated(seed_random))deallocate(seed_random)
! allocate(seed_random(nrandom))
! seed_random=1234567
! call random_seed(put=seed_random)
! deallocate(seed_random)
! call random_number(vect)
call mt_random(vect)
if(verb.AND.mpi_master)write(*,*)"MPI_LANCZOS_EIGH: random initial vector generated:"
else
vect = 1d0
if(verb.AND.mpi_master)write(*,*)"MPI_LANCZOS_EIGH: unitary initial vector generated:"
endif
norm_tmp=dot_product(vect,vect); norm=0d0
call AllReduce_MPI(MpiComm,norm_tmp,norm)
vect=vect/sqrt(norm)
endif
!
!============= LANCZOS LOOP =====================
!
Vin = Vect !save input vector for Eigenvector calculation:
Vout = 0d0
alanc= 0d0
blanc= 0d0
nlanc= 0
!
lanc_loop: do iter=1,Nitermax
call mpi_lanczos_iteration_d(MpiComm,MatVec,iter,vin,vout,a_,b_)
if(abs(b_)<threshold_)exit lanc_loop
!
nlanc=nlanc+1
!
alanc(iter) = a_ ; blanc(iter+1) = b_
!
diag(1:Nlanc) = alanc(1:Nlanc)
subdiag(2:Nlanc) = blanc(2:Nlanc)
call eigh(diag(1:Nlanc),subdiag(2:Nlanc),Ev=Z(:Nlanc,:Nlanc))
!
if(nlanc >= Ncheck_)then
esave(nlanc-(Ncheck_-1))=diag(1)
if(nlanc >= (Ncheck_+1))then
diff=esave(Nlanc-(Ncheck_-1))-esave(Nlanc-(Ncheck_-1)-1)
if(verb.AND.mpi_master)write(*,*)'Iter, E0, deltaE = ',iter,diag(1),diff
if(abs(diff).le.threshold_)exit lanc_loop
endif
endif
enddo lanc_loop
if(nlanc==nitermax)print*,"LANCZOS_SIMPLE: reach Nitermax"
!
!============== END LANCZOS LOOP ======================
!
diag(1:Nlanc) = alanc(1:Nlanc)
subdiag(2:Nlanc) = blanc(2:Nlanc)
call eigh(diag(1:Nlanc),subdiag(2:Nlanc),Ev=Z(:Nlanc,:Nlanc))
!
!Get the Eigenvalues:
egs = diag(1)
!
!Get the Eigenvector:
Vin = Vect
vout= 0d0
vect= 0d0
do iter=1,nlanc
call mpi_lanczos_iteration_d(MpiComm,MatVec,iter,vin,vout,alanc(iter),blanc(iter))
vect = vect + vin*Z(iter,1)
end do
norm_tmp=dot_product(vect,vect); norm=0d0
call Allreduce_MPI(MpiComm,norm_tmp,norm)
vect=vect/sqrt(norm)
if(verb)then
call MatVec(Nloc,vect,vout)
if(mpi_master)write(*,*)"|H*v-E*v|=",sum(abs(vout-egs*vect))/Nloc
endif
Nitermax=Nlanc
!
end subroutine mpi_lanczos_eigh_d