dvdson_eigh_d Subroutine

subroutine dvdson_eigh_d(MatVec,eval,evec,Nblock,Nitermax,Tol)
  !Interface to Matrix-Vector routine:
  interface
     subroutine MatVec(Nloc,vin,vout)
       integer                 :: Nloc
       real(8),dimension(Nloc) :: vin
       real(8),dimension(Nloc) :: vout
     end subroutine MatVec
  end interface
  !Arguments
  real(8),intent(inout)            :: eval(:)![Neigen]
  real(8),intent(inout)            :: evec(:,:)![Ns,Neigen]
  integer,optional                 :: Nblock
  integer,optional                 :: Nitermax
  real(8),optional                 :: Tol
  !Dimensions:
  integer                          :: Ns
  integer                          :: Neigen
  integer                          :: Nume
  integer                          :: Lim
  integer                          :: Ilow,Ihigh
  integer                          :: Niv
  integer                          :: Ncv,Mblock
  integer                          :: MaxIter
  integer                          :: Iwrsz
  integer                          :: Iiwsz
  !Arrays:
  integer,dimension(:),allocatable :: Iselec
  real(8),dimension(:),allocatable :: Work
  integer,dimension(:),allocatable :: Iwork
  !Control Vars:
  integer                          :: i,j,ie,ival,ivec
  real(8)                          :: Crite
  real(8)                          :: Critc
  real(8)                          :: Critr
  real(8)                          :: Ortho
  logical                          :: Hiend
  integer                          :: Nloops
  integer                          :: Nmv
  integer                          :: Ierr
  !
  real(8),dimension(:),allocatable :: Adiag

  Ncv = Neigen ; if(present(Nblock))Ncv = Nblock

  Ns     = size(evec,1)
  Neigen = size(eval)
  if(any(shape(Evec)/=[Ns,Neigen]))stop "dvdson_eigh_d ERROR: Evec wrong shape"
  !
  !
  Nume   = Neigen !Distance of the eigenpair index ik farthest from the extreme
  !               !we look for the first, i.e. lowest, Neigen Eigenpairs here.
  Lim    = min(Nume + 40,Ns)
  Ilow   = 1
  Ihigh  = Neigen
  Niv    = 0                    !to be used for initial vectors
  Mblock = max(1,min(Ncv,Neigen))
  Crite  = 1d-15    ; if(present(tol))Crite=tol
  Critc  = 1d-15
  Critr  = 1d-15
  Ortho  = Critr
  MaxIter = max(200,Nume*40) ; if(present(Nitermax))MaxIter = Nitermax
  Iwrsz  = 2*Ns*Lim + Lim*Lim + (Nume+10)*Lim + Nume
  Iiwsz  = 6*Lim + Nume
  !
  allocate(Iselec(lim))
  allocate(Work(Iwrsz))
  allocate(Iwork(Iiwsz))
  Iselec = 0
  Work   = 0d0
  Iwork  = 0

  allocate(Adiag(Ns))
  if(sum(abs(Evec(:,1)))==0d0)then
     call MatVec_to_diagonal(Adiag)
  else
     Adiag = Evec(:,1)
  endif
  !
  call DVDSON(block_MatVec, Ns, Lim, Adiag,&
       Ilow, Ihigh, Iselec, NIV, Mblock,&
       Crite, Critc, Critr, ortho, Maxiter,&
       Work,Iwrsz,Iwork,Iiwsz,HIEND,NLOOPS,NMV,IERR)

  !POST PROCESSING:
  if(ierr/=0)then
     if(ierr<0)then
        write(*,*)"IERR = ",ierr
        stop "sp_dvdson_eigh: IERR denotes error in DSPEVX (k eigenpairs not converged)"
     else
        write(*,*)"IERR = ",ierr
        stop "sp_dvdson_eigh: IERR denotes some error, see documentation"
        ! if (INT( MOD(IERR,  2)/1) == 0 )stop "sp_dvdson_eigh: IERR signals N < LIM"
        ! If (INT( MOD(IERR,  4)/2) == 0 ) stop "sp_dvdson_eigh: IERR signals LIM < 1"
        ! If (INT( MOD(IERR,  8)/4) == 0  ) stop "sp_dvdson_eigh: IERR signals ISELEC(1)<1, and no range specified"
        ! If (INT( MOD(IERR, 16)/8) == 0  ) stop "sp_dvdson_eigh: IERR signals IHIGH > N (in range or ISELEC)"
        ! If (INT( MOD(IERR, 32)/16) == 0 ) stop "sp_dvdson_eigh: IERR signals IHIGH < ILOW (Invalid range)"       
        ! If (INT( MOD(IERR, 64)/32) == 0 ) stop "sp_dvdson_eigh: IERR signals NEIG >= LIM (Too many wanted)"      
        ! If (INT( MOD(IERR,128)/64) == 0 ) stop "sp_dvdson_eigh: IERR signals Probable duplication in ISELEC"     
        ! If (INT( MOD(IERR,256)/128) == 0) stop "sp_dvdson_eigh: IERR signals NUME >= LIM (max eigen very far)"   
        ! If (INT( MOD(IERR,512)/256) == 0) stop "sp_dvdson_eigh: IERR signals MBLOCK is out of bounds"   
        ! If (INT( MOD(IERR,1024)/512) == 0) stop "sp_dvdson_eigh: IERR signals IWRSZ or IIWSZ is not enough"
        ! If (INT( MOD(IERR,2048)/1024) == 0) stop "sp_dvdson_eigh: IERR signals Orthogonalization Failed"     
        ! If (INT( MOD(IERR,4096)/2048) == 0) stop "sp_dvdson_eigh: IERR signals NLOOPS > MAXITER"        
     endif
  endif

  !Retrieve the Evals,Evecs results from Work array
  do ie=1,Nume
     ivec = 1  + (ie-1)*Ns
     ival = ie + Nume*Ns 
     Eval(ie) = Work(ival)
     Evec(:,ie) = Work(ivec:ivec+Ns-1)
  enddo


contains

  subroutine block_MatVec(N,Ne,B,C)
    integer :: N,Ne
    real(8) :: B(N,Ne)
    real(8) :: C(N,Ne)
    integer :: ie
    !
    do ie=1,Ne
       call MatVec(N,B(:,ie),C(:,ie))
    enddo
    !
  end subroutine block_MatVec

  subroutine MatVec_to_diagonal(diag)
    real(8) :: diag(Ns)
    integer :: i
    real(8) :: vin(Ns),vout(Ns)
    do i=1,Ns
       vin    = 0d0
       vin(i) = 1d0
       call MatVec(Ns,Vin,Vout)
       diag(i)= Vout(i)
    enddo
  end subroutine MatVec_to_diagonal


end subroutine DVDSON_EIGH_D