template<typename Hessenberg, typename Scalar, typename Destination>
void shifted_hessenberg_solve(const Hessenberg& hess, const Scalar& shift, Destination& result)
{
  JacobiRotation<typename Hessenberg::Scalar> rotation;
  typename Hessenberg::PlainObject solver = hess;

  solver.diagonal().array() += shift;

  for (typename Hessenberg::Index k = 1; k < hess.cols(); ++k) {
    rotation.makeGivens(solver.coeff(k-1, k-1), solver.coeff(k, k-1), &solver.coeffRef(k-1, k-1));
    solver.rightCols(hess.cols() - k).applyOnTheLeft(k-1, k, rotation);
    result.applyOnTheLeft(k-1, k, rotation);
  }

  solver.template triangularView<Upper>().solveInPlace(result);
}

template<typename Hessenberg, typename Triangular>
void hessenberg_solve_triangular(const Hessenberg& hess, Triangular& result)
{
  JacobiRotation<typename Hessenberg::Scalar> rotation;
  typename Hessenberg::PlainObject solver = hess;

  for (typename Hessenberg::Index k = 1; k < hess.cols(); ++k) {
    rotation.makeGivens(solver.coeff(k-1, k-1), solver.coeff(k, k-1), &solver.coeffRef(k-1, k-1));
    solver.rightCols(hess.cols() - k).applyOnTheLeft(k-1, k, rotation);
    result.rightCols(hess.cols() - k + 1).applyOnTheLeft(k-1, k, rotation);
  }

  solver.template triangularView<Upper>().solveInPlace(result);
}

template<typename Hessenberg, typename UnitVector>
typename UnitVector::Scalar hessenberg_update_eigenvalue(const Hessenberg& matrix, const UnitVector& eigenvector)
{
  typename UnitVector::PlainObject vector = eigenvector;
  hessenberg_matrix_vector(matrix, vector);
  return eigenvector.dot(vector);
}

template<typename Hessenberg, typename UnitVector>
typename UnitVector::Scalar rayleigh_quotient_iteration(const Hessenberg& matrix, UnitVector& vector)
{
  typedef typename UnitVector::Scalar Scalar;
  Scalar eigenvalue = hessenberg_update_eigenvalue(matrix, vector);

  for (int k = 0; k < 40; ++k) {
    shifted_hessenberg_solve(matrix.template cast<Scalar>(), -eigenvalue, vector);
    vector.normalize();

    Scalar candidate = hessenberg_update_eigenvalue(matrix, vector);

    if (isApprox(eigenvalue, candidate))
      return candidate;

    eigenvalue = candidate;
  }

  return eigenvalue;
}