deigh_generalized Subroutine

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.
  real(8), intent(in)  :: Am(:,:)   ! LHS matrix: Am c = lam Bm c
  real(8), intent(in)  :: Bm(:,:)   ! RHS matrix: Am c = lam Bm c
  real(8), intent(out) :: lam(:)   ! eigenvalues: Am c = lam Bm c
  real(8), intent(out) :: c(:,:)   ! eigenvectors: Am c = lam Bm c; c(i,j) = ith component of jth vec.
  integer              :: n
  ! lapack variables
  integer              :: lwork, liwork, info
  integer, allocatable :: iwork(:)
  real(8), allocatable :: Bmt(:,:), work(:)
  ! solve
  n = size(Am,1)
  call assert_shape(Am, [n, n], "eigh", "Am")
  call assert_shape(Bm, [n, n], "eigh", "B")
  call assert_shape(c, [n, n], "eigh", "c")
  lwork = 1 + 6*n + 2*n**2
  liwork = 3 + 5*n
  allocate(Bmt(n,n), work(lwork), iwork(liwork))
  c = Am; Bmt = Bm  ! Bmt temporaries overwritten by dsygvd
  call dsygvd(1,'V','L',n,c,n,Bmt,n,lam,work,lwork,iwork,liwork,info)
  if (info /= 0) then
     print *, "dsygvd 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: dsygvd'
  end if
end subroutine deigh_generalized