subroutine eig_sym(a, m, w, z, info)
  implicit none
  real(8), intent(in)  :: a(:,:)  !(n,n)
  integer, intent(out) :: m
  real(8), intent(out) :: w(:)    !(n)
  real(8), intent(out) :: z(:,:)  !(n,n)
  integer, intent(out) :: info

  integer, parameter :: nb = 64
  character(1), parameter :: jobz = 'V'
  character(1), parameter :: range = 'I'
  character(1), parameter :: uplo = 'U'
  real(8), allocatable :: a_(:,:)
  real(8) :: abstol
  integer :: n, lda, ldz
  integer :: liwork, lwork
  integer :: il, iu
  real(8) :: vl, vu
  real(8), allocatable :: work(:)
  integer, allocatable :: iwork(:)
  integer, allocatable :: isuppz(:)
  real(8) :: dummy(1) !dummy of work
  integer :: idum(1)  !dummy of iwork

  n = size(a,1)
  allocate(a_(n,n))
  a_ = a

  il = 1
  iu = n

  lda = n
  ldz = n
  m = n
  allocate(isuppz(2*m))

  ! Use routine workspace query to get optimal workspace.
  lwork = -1
  liwork = -1
  call DSYEVR(jobz, range, uplo, &
              n, a_, lda, vl, vu, il, iu, abstol, &
              m, w, z, ldz, isuppz, &
              dummy, lwork, idum, liwork, info)

  ! Make sure that there is enough workspace for block size nb.
  lwork = max((nb+6)*n, nint(dummy(1)))
  liwork = max(10*n, idum(1))
  allocate(work(lwork), iwork(liwork))

  ! Set the absolute error tolerance for eigenvalue. With ABSTOL
  ! set to zero, the default value is used instead.
  abstol = 0.d0

  ! Solve the symmetric eigenvalue problem.
  call DSYEVR(jobz, range, uplo, &
              n, a_, lda, vl, vu, il, iu, abstol, &
              m, w, z, ldz, isuppz, &
              work, lwork, iwork, liwork, info)

  if( m < n )then
    w(m+1:n) = 0.d0
    z(:,m+1:n) = 0.d0
  endif
end subroutine eig_sym