modes Subroutine

subroutine modes(K, lambda, phi, tol)
  !! Ascending eigenvalues and orthonormal eigenvectors of symmetric K; eigenvalues
  !! below tol*max|lambda| (default 1e-10) are set to exactly zero
  ! not pure: stdlib eigh is not
  real(dp), intent(in) :: K(:,:)
  real(dp), intent(out) :: lambda(:), phi(:,:)
  real(dp), intent(in), optional :: tol
  real(dp) :: A(size(K,1), size(K,2)), t

  t = 1e-10_dp
  if (present(tol)) t = tol
  A = K   ! eigh overwrites its input
  call eigh(A, lambda, vectors=phi)
  where (abs(lambda) < t*maxval(abs(lambda))) lambda = 0
end subroutine modes