lanczos_eigh_d Subroutine

subroutine lanczos_eigh_d(MatVec,Egs,Vect,Nitermax,iverbose,threshold,ncheck,vrandom)
  interface
     subroutine MatVec(Nloc,vin,vout)
       integer                 :: Nloc
       real(8),dimension(Nloc) :: vin
       real(8),dimension(Nloc) :: vout
     end subroutine MatVec
  end interface
  real(8)                              :: egs
  real(8),dimension(:)                 :: vect
  integer                              :: Nitermax
  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_
  real(8)                              :: norm,diff
  integer                              :: i,ierr
  logical                              :: vran=.true.
  !
  if(present(iverbose))verb=iverbose
  if(present(threshold))threshold_=threshold
  if(present(ncheck))ncheck_=ncheck
  if(present(vrandom))vran=vrandom
  !
  norm=dot_product(vect,vect)
  if(norm==0d0)then
     if(vran)then
        call mt_random(vect)
        if(verb)write(*,*)"LANCZOS_EIGH_D: random initial vector generated:"
     else
        vect = 1d0
        if(verb)write(*,*)"LANCZOS_EIGH_D: unitary initial vector generated:"
     endif
     vect=vect/sqrt(dot_product(vect,vect))
  endif
  !
  !============= LANCZOS LOOP =====================
  !
  vin = vect
  vout= 0.d0
  alanc=0.d0
  blanc=0.d0
  nlanc=0
  !
  lanc_loop: do iter=1,Nitermax
     call lanczos_iteration_d(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)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 lanczos_iteration_d(MatVec,iter,vin,vout,alanc(iter),blanc(iter))
     vect = vect + vin*Z(iter,1)
  end do
  norm=sqrt(dot_product(vect,vect))
  vect=vect/norm
  !
  if(verb)then
     call MatVec(size(vect),vect,vout)
     write(*,*)"|H*v-E*v|=",sum(abs(vout-egs*vect))/size(vect)
  endif
  Nitermax=Nlanc
  !
end subroutine lanczos_eigh_d