mpi_lanczos_eigh_c Subroutine

subroutine mpi_lanczos_eigh_c(MpiComm,MatVec,Egs,Vect,Nitermax,iverbose,threshold,ncheck,vrandom)
  integer                              :: MpiComm
  interface 
     subroutine MatVec(Nloc,vin,vout)
       integer    :: Nloc
       complex(8) :: vin(Nloc)
       complex(8) :: vout(Nloc)
     end subroutine MatVec
  end interface
  real(8)                              :: egs
  complex(8),dimension(:)              :: vect
  integer                              :: Nitermax
  !
  integer                              :: Nloc
  real(8),optional                     :: threshold
  integer,optional                     :: ncheck
  logical,optional                     :: iverbose
  logical,optional                     :: vrandom
  !
  complex(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,ran(2),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)
        ! do i=1,Nloc
        !    call random_number(ran)
        !    vect(i)=dcmplx(ran(1),ran(2))
        ! enddo
        call mt_random(vect)
        if(verb.AND.mpi_master)write(*,*)"MPI_LANCZOS_EIGH: random initial vector generated:"
     else
        vect = one
        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
  vout= zero
  alanc=0d0
  blanc=0d0
  nlanc=0
  !
  lanc_loop: do iter=1,Nitermax
     call mpi_lanczos_iteration_c(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= zero
  vect= zero
  do iter=1,nlanc
     call mpi_lanczos_iteration_c(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_c