zeigh_generalized Subroutine

subroutine zeigh_generalized(Am, Bm, lam, c)
  ! solves generalized eigen value problem for all eigenvalues and eigenvectors
  ! Am must by hermitian, Bm hermitian positive definite.
  ! Only the lower triangular part of Am and Bm is used.
  complex(8), intent(in)  :: Am(:,:)   ! LHS matrix: Am c = lam Bm c
  complex(8), intent(in)  :: Bm(:,:)   ! RHS matrix: Am c = lam Bm c
  real(8), intent(out)    :: lam(:)      ! eigenvalues: Am c = lam Bm c
  complex(8), intent(out) :: c(:,:)   ! eigenvectors: Am c = lam Bm c; c(i,j) = ith component of jth vec.
  ! lapack variables
  integer                 :: info, liwork, lrwork, lwork, n
  integer, allocatable    :: iwork(:)
  real(8), allocatable    :: rwork(:)
  complex(8), allocatable :: Bmt(:,:), work(:)
  n = size(Am,1)
  call assert_shape(Am, [n, n], "eigh", "Am")
  call assert_shape(Bm, [n, n], "eigh", "Bm")
  call assert_shape(c, [n, n], "eigh", "c")
  lwork = 2*n + n**2
  lrwork = 1 + 5*N + 2*n**2
  liwork = 3 + 5*n
  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
        print *, "the", -info, "-th argument had an illegal value"
     else if (info <= n) then
        print *, "the algorithm failed to compute an eigenvalue while working"
        print *, "on the submatrix lying in rows and columns", 1d0*info/(n+1)
        print *, "through", mod(info, n+1)
     else
        print *, "The leading minor of order ", info-n, &
             "of B is not positive definite. The factorization of B could ", &
             "not be completed and no eigenvalues or eigenvectors were computed."
     end if
     stop 'deigh_generalized error: zhegvd'
  end if
end subroutine zeigh_generalized