MetricChecker.h

#ifndef LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_METRICCHECKER_H
#define LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_METRICCHECKER_H
#include <Eigen/Dense>
#include <molpro/linalg/array/Span.h>
#include <molpro/linalg/itsolv/ArrayHandlers.h>
#include <molpro/linalg/itsolv/subspace/IXSpace.h>
#include <molpro/linalg/itsolv/subspace/Matrix.h>
#include <molpro/linalg/itsolv/subspace/util.h>
#include <molpro/linalg/itsolv/wrap.h>
#include <molpro/linalg/scalar_traits.h>
 
#include <cmath>
#include <format>
#include <optional>
#include <stdexcept>
#include <utility>
 
namespace molpro::linalg::itsolv::detail {
 
template <class R, class Q, class P>
class MetricChecker {
public:
  using value_type = typename array::ArrayHandler<R, R>::value_type;
  using value_type_abs = typename array::ArrayHandler<R, R>::value_type_abs;
 
  static value_type_abs hermiticity_tolerance() { return precision_scaled<value_type>(1e-6); }
 
  static value_type_abs positivity_tolerance() { return precision_scaled<value_type>(1e-8); }
 
  void reset() {
    m_reference.reset();
    m_pp_factor.reset();
  }
 
  void reset_pspace() { m_pp_factor.reset(); }
 
  void check_hermitian(const CVecRef<R>& params, const CVecRef<R>& metric_actions,
                       const subspace::IXSpace<R, Q, P>& xspace, ArrayHandlers<R, Q, P>& handlers) const {
    const Q* q = nullptr;
    const Q* metric_action_q = nullptr;
    value_type qmq = 0; // <q|M q>
    if (xspace.stores_metric_actions()) {
      if (xspace.dimensions().nQ == 0)
        return;
      q = &xspace.cparamsq().front().get();
      metric_action_q = &xspace.cmetric_actionsq().front().get();
      qmq = xspace.data.at(subspace::EqnData::M)(xspace.dimensions().oQ(), xspace.dimensions().oQ());
    } else if (m_reference) {
      q = &m_reference->q;
      metric_action_q = &m_reference->metric_action_q;
      qmq = m_reference->qmq;
    } else
      return;
    // <x|M q> and <M x|q> from one pass over q and M q
    auto left = params;
    left.insert(left.end(), metric_actions.begin(), metric_actions.end());
    const auto overlaps =
        subspace::util::overlap(left, CVecRef<Q>{std::cref(*metric_action_q), std::cref(*q)}, handlers.rq());
    for (size_t i = 0; i < params.size(); ++i) {
      const auto xmq = overlaps(i, 0);
      const auto mxq = overlaps(params.size() + i, 1);
      // the size of these elements, which the metric bounds with its values on x and q
      const auto scale = std::sqrt(std::abs(handlers.rr().dot(params[i], metric_actions[i]) * qmq));
      if (std::abs(xmq - mxq) > hermiticity_tolerance() * scale)
        throw std::runtime_error(std::format("The metric is not Hermitian: <x|M q> = {:.6e} differs from <M x|q> = "
                                             "{:.6e} for a new parameter x and a subspace vector q",
                                             double(std::abs(xmq)), double(std::abs(mxq))));
    }
  }
 
  void remember(const CVecRef<R>& params, const CVecRef<R>& metric_actions, const subspace::IXSpace<R, Q, P>& xspace,
                ArrayHandlers<R, Q, P>& handlers) {
    if (xspace.stores_metric_actions() || params.empty())
      return;
    if (m_reference) {
      handlers.qr().copy(m_reference->q, params.front());
      handlers.qr().copy(m_reference->metric_action_q, metric_actions.front());
    } else
      m_reference.emplace(Reference{handlers.qr().copy(params.front()), handlers.qr().copy(metric_actions.front())});
    m_reference->qmq = handlers.rr().dot(params.front(), metric_actions.front());
  }
 
  static void check_hermitian(const array::Span<value_type>& matrix, const size_t n) {
    // the diagonal elements, which have to be real, are compared with themselves
    for (size_t i = 0; i < n; ++i)
      for (size_t j = 0; j <= i; ++j)
        if (std::abs(matrix[i * n + j] - conjugate(matrix[j * n + i])) >
            hermiticity_tolerance() * std::sqrt(std::abs(matrix[i * n + i] * matrix[j * n + j])))
          throw std::runtime_error("The metric is not Hermitian: its matrix in the P space is not");
  }
 
  void check_positive_subspace(const subspace::Matrix<value_type>& metric, const subspace::Matrix<value_type>& overlap,
                               const size_t nP) {
    using matrix_type = Eigen::Matrix<value_type, Eigen::Dynamic, Eigen::Dynamic>;
    using row_major_type = Eigen::Matrix<value_type, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
    const auto n = Eigen::Index(metric.rows());
    const auto np = Eigen::Index(nP);
    if (np == 0 || np >= n)
      return check_positive(metric, overlap);
    const matrix_type m_rows = Eigen::Map<const row_major_type>(metric.data().data(), n, n);
    const matrix_type m = (m_rows + m_rows.adjoint()) / value_type(2);
    if (!m_pp_factor || m_pp_factor->rows() != np) {
      const Eigen::LLT<matrix_type> llt(m.topLeftCorner(np, np));
      if (llt.info() != Eigen::Success)
        return check_positive(metric, overlap);
      m_pp_factor = matrix_type(llt.matrixL());
    }
    const matrix_type y = m_pp_factor->template triangularView<Eigen::Lower>().solve(m.topRightCorner(np, n - np));
    const matrix_type schur = m.bottomRightCorner(n - np, n - np) - y.adjoint() * y;
    if (Eigen::LLT<matrix_type>(schur).info() == Eigen::Success)
      return;
    check_positive(metric, overlap);
  }
 
  static void check_positive(const subspace::Matrix<value_type>& metric, const subspace::Matrix<value_type>& overlap) {
    using matrix_type = Eigen::Matrix<value_type, Eigen::Dynamic, Eigen::Dynamic>;
    using row_major_type = Eigen::Matrix<value_type, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
    const auto n = Eigen::Index(metric.rows());
    if (n == 0)
      return;
    const matrix_type m_rows = Eigen::Map<const row_major_type>(metric.data().data(), n, n);
    const matrix_type m = (m_rows + m_rows.adjoint()) / value_type(2);
    if (Eigen::LLT<matrix_type>(m).info() == Eigen::Success)
      return;
    const matrix_type s = Eigen::Map<const row_major_type>(overlap.data().data(), n, n);
    const Eigen::SelfAdjointEigenSolver<matrix_type> s_eigen(s);
    const value_type_abs s_threshold = precision_scaled<value_type>(1e-10) * s_eigen.eigenvalues().maxCoeff();
    Eigen::Index rank = 0;
    for (Eigen::Index i = 0; i < n; ++i)
      if (s_eigen.eigenvalues()(i) > s_threshold)
        ++rank;
    if (rank == 0)
      return;
    // the eigenvalues come in ascending order, so the retained directions are the trailing ones
    const auto retained = s_eigen.eigenvectors().rightCols(rank);
    const matrix_type orthonormal =
        retained * s_eigen.eigenvalues().tail(rank).cwiseSqrt().cwiseInverse().template cast<value_type>().asDiagonal();
    const Eigen::SelfAdjointEigenSolver<matrix_type> m_eigen(orthonormal.adjoint() * m * orthonormal,
                                                             Eigen::EigenvaluesOnly);
    const value_type_abs lowest = m_eigen.eigenvalues().minCoeff();
    const value_type_abs largest = m_eigen.eigenvalues().cwiseAbs().maxCoeff();
    // the rounding errors of the elements of the metric matrix, amplified by the orthonormalisation
    const value_type_abs rounding = precision_scaled<value_type>(1e-13) * m.norm() / s_eigen.eigenvalues()(n - rank);
    if (lowest < -positivity_tolerance() * largest - rounding)
      throw std::runtime_error(std::format("The metric is not positive definite: it has the Rayleigh quotient {:.3e} "
                                           "on the subspace, against a largest one of {:.3e}",
                                           double(lowest), double(largest)));
  }
 
private:
  struct Reference {
    Q q;
    Q metric_action_q;
    value_type qmq = 0;
  };
  std::optional<Reference> m_reference;
  std::optional<Eigen::Matrix<value_type, Eigen::Dynamic, Eigen::Dynamic>> m_pp_factor;
};
 
} // namespace molpro::linalg::itsolv::detail
 
#endif // LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_METRICCHECKER_H