anchor.diff:      end if
     end function zeigvals
   
  -  subroutine eigh_generalized(Am, Bm, lam, c)
  +  subroutine deigh_generalized(Am, Bm, lam, c)
       ! solves generalized eigen value problem for all eigenvalues and eigenvectors
       ! Am must by symmetric, Bm symmetric positive definite.
       ! Only the lower triangular part of Am and Bm is used.
  
positive.code:     integer :: info, liwork, lrwork, lwork, n
      integer, allocatable :: iwork(:)
      real(dp), allocatable :: rwork(:)
      complex(dp), allocatable :: Amt(:,:), Bmt(:,:), work(:)
  
      ! solve
      n = size(Am,1)
      lwork = 2*n + n**2
      lrwork = 1 + 5*N + 2*n**2
      liwork = 3 + 5*n
      allocate(Amt(n,n), Bmt(n,n), work(lwork), rwork(lrwork), iwork(liwork))
      Amt = Am; Bmt = Bm  ! Amt, Bmt temporaries overwritten by zhegvd
      call zhegvd(1,'V','L',n,Amt,n,Bmt,n,lam,work,lwork,rwork,lrwork,iwork,liwork,info)
      if (info /= 0) then
         print *, "zhegvd returned info =", info
         if (info < 0) then
  
negative.code:     ldvr = n
      lwork = 8*n  ! TODO: can this size be optimized? query first?
      lrwork = 2*n
      allocate(vl(ldvl,n), vr(ldvr,n), work(lwork), rwork(lrwork))
      c = A
      call zgeev('N', 'V', n, c, lda, lam, vl, ldvl, vr, ldvr, work, &
           lwork, rwork, info)
      if(info /= 0) then
         print *, "zgeev returned info = ", info
  