anchor.diff:      integer :: info, liwork, lrwork, lwork, n
       integer, allocatable :: iwork(:)
       real(dp), allocatable :: rwork(:)
  -    complex(dp), allocatable :: Amt(:,:), Bmt(:,:), work(:)
  +    complex(dp), allocatable :: 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)
  +    allocate(Bmt(n,n), work(lwork), rwork(lrwork), iwork(liwork))
  +    c = Am; Bmt = Bm  ! Bmt temporary overwritten by zhegvd
  +    call zhegvd(1,'V','L',n,c,n,Bmt,n,lam,work,lwork,rwork,lrwork,iwork,liwork,info)
       if (info /= 0) then
          print *, "zhegvd returned info =", info
          if (info < 0) then
  
positive.code:     integer :: info, lda, liwork, lrwork, lwork, n
      integer, allocatable :: iwork(:)
      real(dp), allocatable :: rwork(:)
      complex(dp), allocatable :: Amt(:, :), work(:)
  
      ! use LAPACK's ZHEEVD routine
      n = size(Am, 1)
      lda = max(1, n)
      lwork = 2*n + n**2
      lrwork = 1 + 5*n + 2*n**2
      liwork = 3 + 5*n
      allocate(Amt(lda,n), work(lwork), rwork(lrwork), &
           iwork(liwork))
      Amt = Am
      call ZHEEVD("V", "L", n, Amt, lda, lam, work, lwork, rwork, lrwork, &
           iwork, liwork, info)
      if (info /= 0) then
         print *, "zheevd returned info =", info
  
negative.code:        end if
         call stop_error('eigh: zheevd error')
      end if
    end subroutine zeigh_simple
  
    function dinv(Am) result(Bm)
  