subroutine diag_mat_per_fock_degen_core(fock_diag, mat_ref, listcore,ncore, n, thr_d, thr_nd, thr_deg, leigvec, reigvec, eigval)


  BEGIN_DOC
  !
  ! subroutine that diagonalizes a matrix mat_ref BY BLOCK
  !
  ! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag"
  !
  ! the elements of listcore are untouched 
  !
  ! examples : all elements having degeneracy 1 in fock_diag (i.e. not being degenerated) will be treated together
  !
  !          : all elements having degeneracy 2 in fock_diag (i.e. two elements are equal) will be treated together
  !
  !          : all elements having degeneracy 3 in fock_diag (i.e. two elements are equal) will be treated together
  !
  ! etc... the advantage is to guarentee no spurious mixing because of numerical problems. 
  !
  END_DOC

  implicit none
  integer,          intent(in)  :: n,ncore, listcore(ncore)
  double precision, intent(in)  :: fock_diag(n), mat_ref(n,n), thr_d, thr_nd, thr_deg
  double precision, intent(out) :: leigvec(n,n), reigvec(n,n), eigval(n)

  integer                       :: n_degen_list, n_degen,size_mat, i, j, k, icount, m, index_degen
  integer                       :: ii, jj, i_good, j_good, n_real
  integer                       :: icount_eigval
  logical,          allocatable :: is_ok(:)
  integer,          allocatable :: list_degen(:,:), list_same_degen(:)
  integer,          allocatable :: iorder(:), list_degen_sorted(:)
  double precision, allocatable :: leigvec_unsrtd(:,:), reigvec_unsrtd(:,:), eigval_unsrtd(:)
  double precision, allocatable :: mat_tmp(:,:), eigval_tmp(:), leigvec_tmp(:,:), reigvec_tmp(:,:)

  allocate(leigvec_unsrtd(n,n), reigvec_unsrtd(n,n), eigval_unsrtd(n))
  leigvec_unsrtd = 0.d0
  reigvec_unsrtd = 0.d0
  eigval_unsrtd  = 0.d0

  ! obtain degeneracies 
  allocate(list_degen(n,0:n))
  call give_degen_full_listcore(fock_diag, n, listcore, ncore, thr_deg, list_degen, n_degen_list)

  allocate(iorder(n_degen_list), list_degen_sorted(n_degen_list))
  do i = 1, n_degen_list
    n_degen = list_degen(i,0) 
    list_degen_sorted(i) = n_degen
    iorder(i) = i
  enddo

  ! sort by number of degeneracies 
  call isort(list_degen_sorted, iorder, n_degen_list)

  allocate(is_ok(n_degen_list))
  is_ok = .True.
  icount_eigval = 0

  ! loop over degeneracies 
  do i = 1, n_degen_list
    if(.not.is_ok(i)) cycle

    is_ok(i) = .False.
    n_degen  = list_degen_sorted(i)

    
    if(n_degen.ge.1000)then
     print*,'core orbital '
    else
     print *, ' diagonalizing for n_degen = ', n_degen
    endif

    k = 1

   ! group all the entries having the same degeneracies 
!!  do while (list_degen_sorted(i+k)==n_degen)
    do m = i+1, n_degen_list
      if(list_degen_sorted(m)==n_degen) then
        is_ok(i+k) = .False.
        k += 1
      endif
    enddo

    print *, ' number of identical degeneracies = ', k
    if(n_degen.ge.1000)then
     n_degen = 1
    endif
    size_mat = k*n_degen  
    print *, ' size_mat = ', size_mat
    allocate(mat_tmp(size_mat,size_mat), list_same_degen(size_mat))
    allocate(eigval_tmp(size_mat), leigvec_tmp(size_mat,size_mat), reigvec_tmp(size_mat,size_mat))
    ! group all the elements sharing the same degeneracy
    icount = 0
    do j = 1, k ! jth set of degeneracy
      index_degen = iorder(i+j-1)
      do m = 1, n_degen
        icount += 1
        list_same_degen(icount) = list_degen(index_degen,m)
      enddo
    enddo

    print *, ' list of elements '
    do icount = 1, size_mat
      print *, icount, list_same_degen(icount)
    enddo

    ! you copy subset of matrix elements having all the same degeneracy in mat_tmp
    do ii = 1, size_mat
      i_good = list_same_degen(ii)
      do jj = 1, size_mat
        j_good = list_same_degen(jj)
        mat_tmp(jj,ii) = mat_ref(j_good,i_good)
      enddo
    enddo

    call non_hrmt_bieig( size_mat, mat_tmp, thr_d, thr_nd &
                       , leigvec_tmp, reigvec_tmp         & 
                       , n_real, eigval_tmp )

    do ii = 1, size_mat
      icount_eigval += 1
      eigval_unsrtd(icount_eigval) = eigval_tmp(ii) ! copy eigenvalues 
      do jj = 1, size_mat ! copy the eigenvectors 
        j_good = list_same_degen(jj)
        leigvec_unsrtd(j_good,icount_eigval) = leigvec_tmp(jj,ii)
        reigvec_unsrtd(j_good,icount_eigval) = reigvec_tmp(jj,ii)
      enddo
    enddo

    deallocate(mat_tmp, list_same_degen)
    deallocate(eigval_tmp, leigvec_tmp, reigvec_tmp)
  enddo

  if(icount_eigval .ne. n) then
    print *, ' pb !! (icount_eigval.ne.n)'
    print *, ' icount_eigval,n', icount_eigval, n
    stop
  endif
 
  deallocate(iorder)
  allocate(iorder(n))
  do i = 1, n
    iorder(i) = i
  enddo
  call dsort(eigval_unsrtd, iorder, n)

  do i = 1, n
    print*,'sorted eigenvalues '
    i_good = iorder(i)
    eigval(i) = eigval_unsrtd(i)
    print*,'i,eigval(i) = ',i,eigval(i)
    do j = 1, n
      leigvec(j,i) = leigvec_unsrtd(j,i_good)
      reigvec(j,i) = reigvec_unsrtd(j,i_good)
    enddo
  enddo

  deallocate(leigvec_unsrtd, reigvec_unsrtd, eigval_unsrtd)
  deallocate(list_degen)
  deallocate(iorder, list_degen_sorted)
  deallocate(is_ok)

end

! ---

subroutine give_degen_full_listcore(A, n, listcore, ncore, thr, list_degen, n_degen_list)

  BEGIN_DOC
  ! you enter with an array A(n) and spits out all the elements degenerated up to thr
  !
  ! the elements of A(n) DON'T HAVE TO BE SORTED IN THE ENTRANCE: TOTALLY GENERAL 
  !
  ! list_degen(i,0) = number of degenerate entries 
  !
  ! list_degen(i,1) = index of the first degenerate entry
  !
  ! list_degen(i,2:list_degen(i,0)) = list of all other dengenerate entries 
  !
  ! if list_degen(i,0) == 1 it means that there is no degeneracy for that element
  !
  ! if list_degen(i,0) >= 1000 it means that it is core orbitals
  END_DOC

  implicit none

  double precision, intent(in)  :: A(n)
  double precision, intent(in)  :: thr
  integer,          intent(in)  :: n,ncore, listcore(ncore)
  integer,          intent(out) :: list_degen(n,0:n), n_degen_list
  integer                       :: i, j, icount, icheck,k
  logical, allocatable          :: is_ok(:)


  allocate(is_ok(n))
  n_degen_list = 0
  is_ok = .True.
  ! you first exclude the "core" orbitals  
  do i = 1, ncore
   j=listcore(i)
   is_ok(j) = .False.
  enddo
  do i = 1, n
    if(.not.is_ok(i)) cycle
    n_degen_list +=1
    is_ok(i) = .False.
    list_degen(n_degen_list,1) = i
    icount = 1
    do j = i+1, n
      if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j)) then
        is_ok(j) = .False.
        icount += 1
        list_degen(n_degen_list,icount) = j
      endif
    enddo

    list_degen(n_degen_list,0) = icount
  enddo
  ! you set all the core orbitals as separate entities 
  icheck = 0
  do i = 1, n_degen_list
    icheck += list_degen(i,0)
  enddo
  if(icheck.ne.(n-ncore))then
    print *, ' pb ! :: icheck.ne.n-ncore'
    print *, icheck, n-ncore
    stop
  endif
  k=1000
  do i = 1, ncore
   n_degen_list+= 1
   j=listcore(i)
   list_degen(n_degen_list,1) = i
   list_degen(n_degen_list,0) = k
   k+=1
  enddo



end

! ---