module DPA
  implicit none
  real*8,dimension(:),allocatable        :: ev0,vec0
  real*8,parameter    :: RPRECI=2.22045D-16, RBIGST=4.49423D+307
  integer,parameter    :: IZERO=0, IONE=1, ITWO=2, ITHREE=3,IFOUR=4
  real*8,parameter     :: ZERO=0.D0, ONE=1.D0, TWO=2.D0,HALF=0.5D0
  
contains
 SUBROUTINE LA_DIAG(N,nmode,hessian0,ev,vec)!hessian -- not allocatable
  implicit none
  integer,intent(in)     ::  N,NMODE      !N IS THE DIMENSION OF DOF
  real*8,intent(in)      ::  hessian0(*)
!  real*8,intent(out)    ::  ev(*),vec(*)
  real*8,allocatable,dimension(:),intent(out)  :: ev,vec

!  real*8,optional,intent(out)  ::  vec(*)
  CHARACTER             ::  JOBZ, RANGE,UPLO
  INTEGER               ::  IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M
  INTEGER               ::  ISUPPZ(2*N), IWORK(15*N)
  DOUBLE PRECISION      ::  W(N), WORK(35*N)
  DOUBLE PRECISION      ::  Z(N,NMODE)
  REAL*8                ::  ABSTOL, VL, VU
  INTEGER               ::  I,J,K,isallocated 
  real*8,allocatable,dimension(:,:) :: A   !real hessian
  
  !
  JOBZ='V'; RANGE='I'; UPLO='U';
  IL=1; IU=nmode; M=IU-IL+1;
  ABSTOL=1.D-12;
  LDA=N; LDZ=N;
  LWORK=35*N; LIWORK=15*N
  !
  allocate(A(n,n),stat=isallocated)
  if(isallocated /= 0) call err_stop('FAILED ALLO. HESSIAN(N,N)');
  if(allocated(ev)) deallocate(ev)
  allocate(ev(n),stat=isallocated)          !reate atom link for a residue
  if(isallocated /= 0) call err_stop('CAN NOT ALLOCATE MEMORY FOR EV(N) IN DIAG');
  if(allocated(vec)) deallocate(vec)
  allocate(vec(n*nmode),stat=isallocated)          !reate atom link for a residue
  if(isallocated /= 0) call err_stop('CAN NOT ALLOCATE MEMORY FOR VEC(N*NMODE) IN DIAG');
 
  K=1
  DO I=1,N
     DO J=I,N
        A(I,J)=HESSIAN0(K)
        K=K+1
     ENDDO
  ENDDO
!  print*,"ENTER_DSYEVR> ", N,A(1,1),A(1,N), A(N,N)
  call DSYEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, &
           ABSTOL, M, W, Z, LDZ,ISUPPZ, WORK, LWORK,IWORK, LIWORK, INFO )
  !
!  print*, "CALC_DSYEVR: ", INFO
  !
  do i=1,nmode
     ev(i)=w(i)
  enddo
!  if(present(vec)) then
     k=1
     do i=1,nmode
        do j=1,n
           vec(k)=z(j,i)
           k=k+1
        enddo
     enddo
!  endif
!
  return
 end subroutine la_diag
 subroutine err_stop(info,id)
    !
    ! PROGRAM STOPS AT ERROR
    !
    implicit none
    character*(*),intent(in)  :: info
    integer,optional,intent(in)   :: id
    
    if(.not. present(id)) write(*,'(A,A)') 'ERRSTP> '//info
    if(present(id).and.id<1E08) write(*,'(A,I8)') 'ERRSTP> '//info,id
    if(present(id).and.id>=1E08) write(*,'(A,I12)') 'ERRSTP> '//info,id
    stop
    return
 end subroutine err_stop
    !
end module dpa
subroutine err_stop(info)
    !
    ! PROGRAM STOPS AT ERROR
    !
    implicit none
    character*(*),intent(in)  :: info
   
    write(*,'(A,A)') 'ERRSTP> '//info
    stop
    
    return
end subroutine err_stop
