!---------------------------------------------------------------------
!Purpose: use plain lanczos to get the groundstate energy
!---------------------------------------------------------------------
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






!---------------------------------------------------------------------
!Purpose: use simple Lanczos to tri-diagonalize a matrix H (defined 
! in the H*v function).
!---------------------------------------------------------------------
subroutine mpi_lanczos_tridiag_d(MpiComm,MatVec,vin,alanc,blanc,threshold)
  integer                                      :: MpiComm
  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),dimension(:),intent(inout)           :: vin !Nloc
  real(8),dimension(size(vin))                 :: vout,vtmp
  real(8),dimension(:),intent(inout)           :: alanc
  real(8),dimension(size(alanc)),intent(inout) :: blanc
  integer                                      :: Nitermax,Nloc,i
  integer                                      :: iter
  real(8)                                      :: a_,b_
  real(8),optional                             :: threshold
  !
  logical                                      :: mpi_master
  !
  if(MpiComm==Mpi_Comm_Null)return
  !
  mpi_master=get_master_MPI(MpiComm)
  !
  Nloc = size(vin)
  !
  if(present(threshold))threshold_=threshold
  !
  vtmp = vin
  Nitermax = size(alanc)
  a_=0d0
  b_=0d0
  vout=0d0
  do iter=1,Nitermax
     call mpi_lanczos_iteration_d(MpiComm,MatVec,iter,vtmp,vout,a_,b_)
     alanc(iter)=a_
     if(abs(b_)<threshold_)exit
     if(iter<nitermax)blanc(iter+1)=b_
  enddo
  if(iter==nitermax.AND.mpi_master)write(*,"(A)")"MPI_LANCZOS_TRIDIAG_D: reach Nitermax"
end subroutine mpi_lanczos_tridiag_d




!---------------------------------------------------------------------
!Purpose: plain homebrew lanczos iteration (no orthogonalization)
!note: the a,b variables are real, even in the complex matrix case
!to understand why check out the Gollub-Van Loan textbook.
!a it is easy: hermiticity->diag\in\RRR
!b: is fixed by requiring |b|^2 = <v,v> thus you can only fix the 
!the absolute value. A lemma shows that the phase can be chosen 
!identically zero
!MPI VERSION
!---------------------------------------------------------------------
subroutine mpi_lanczos_iteration_d(MpiComm,MatVec,iter,vin,vout,alfa,beta)
  integer                                    :: MpiComm
  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),dimension(:),intent(inout)         :: vin !Nloc
  real(8),dimension(size(vin)),intent(inout) :: vout
  real(8),dimension(size(vin))               :: tmp
  real(8),intent(inout)                      :: alfa,beta
  real(8)                                    :: atmp,btmp
  integer                                    :: iter,nloc
  real(8)                                    :: norm,norm_tmp
  !
  logical                                    :: mpi_master
  !
  if(MpiComm==Mpi_Comm_Null)return
  !
  nloc=size(vin)
  !
  mpi_master=get_master_MPI(MpiComm)
  !
  if(iter==1)then
     norm_tmp=dot_product(vin,vin)
     norm = 0d0 ; call AllReduce_MPI(MpiComm,norm_tmp,norm)
     if(mpi_master.AND.norm==0d0)stop "MPI_LANCZOS_ITERATION_D: norm = 0!!"
     vin=vin/sqrt(norm)
  else
     tmp = vin
     vin = vout/beta
     vout= -beta*tmp
  endif
  call MatVec(nloc,vin,tmp)
  vout = vout + tmp
  atmp = dot_product(vin,vout) ; alfa = 0d0; call AllReduce_MPI(MpiComm,atmp,alfa)
  vout = vout - alfa*vin
  btmp = dot_product(vout,vout); beta = 0d0; call AllReduce_MPI(MpiComm,btmp,beta)
  beta = sqrt(beta)
end subroutine mpi_lanczos_iteration_d
