QSpace.h

#ifndef LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_SUBSPACE_QSPACE_H
#define LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_SUBSPACE_QSPACE_H
#include <array>
#include <algorithm>
#include <cassert>
#include <functional>
#include <list>
 
#include <molpro/linalg/itsolv/ArrayHandlers.h>
#include <molpro/linalg/itsolv/Logger.h>
#include <molpro/linalg/itsolv/subspace/Dimensions.h>
#include <molpro/linalg/itsolv/subspace/PSpace.h>
#include <molpro/linalg/itsolv/subspace/SubspaceData.h>
#include <molpro/linalg/itsolv/wrap.h>
#include <molpro/linalg/itsolv/util.h>
#include <molpro/linalg/itsolv/helper-implementation.h>
 
namespace molpro::linalg::itsolv::subspace {
 
namespace qspace {
template <class Q>
struct QParam {
  std::unique_ptr<Q> param;         
  std::unique_ptr<Q> action;        
  size_t id;                        
  std::unique_ptr<Q> metric_action; 
};
 
template <class Q, class ForwardIt>
auto cwrap_params(ForwardIt begin, ForwardIt end) {
  auto param_action = std::array<CVecRef<Q>, 2>{};
  for (auto it = begin; it != end; ++it) {
    param_action[0].emplace_back(*it->param);
    param_action[1].emplace_back(*it->action);
  }
  return param_action;
}
 
template <class Q, class ForwardIt>
auto wrap_params(ForwardIt begin, ForwardIt end) {
  auto param_action = std::array<VecRef<Q>, 2>{};
  for (auto it = begin; it != end; ++it) {
    param_action[0].emplace_back(*it->param);
    param_action[1].emplace_back(*it->action);
  }
  return param_action;
}
} // namespace qspace
 
template <class Rt, class Qt, class Pt>
struct QSpace {
  using R = Rt;
  using Q = Qt;
  using P = Pt;
  using value_type = typename array::ArrayHandler<R, R>::value_type;
  using value_type_abs = typename array::ArrayHandler<R, R>::value_type_abs;
 
  explicit QSpace(std::shared_ptr<ArrayHandlers<R, Q, P>> handlers, std::shared_ptr<Logger> logger)
      : m_handlers(std::move(handlers)), m_logger(std::move(logger)) {}
 
  void set_logger(std::shared_ptr<Logger> logger) { m_logger = std::move(logger); }
 
  template <typename T>
  void update(const CVecRef<R>& params, const CVecRef<R>& actions, const SubspaceData<T>& qq_in,
              const SubspaceData<T>& qx_in, const SubspaceData<T>& xq_in, const Dimensions& dims,
              SubspaceData<T>& old_data, const CVecRef<R>& metric_actions = {}) {
    m_logger->trace("QSpace::update");
    assert(metric_actions.empty() || metric_actions.size() == params.size());
    auto it_begin = m_params.begin();
    for (size_t i = 0; i < params.size(); ++i) {
      m_params.emplace(
          it_begin,
          qspace::QParam<Q>{std::make_unique<Q>(m_handlers->qr().copy(params.at(i))),
                            std::make_unique<Q>(m_handlers->qr().copy(actions.at(i))), m_unique_id++,
                            metric_actions.empty() ? nullptr
                                                   : std::make_unique<Q>(m_handlers->qr().copy(metric_actions.at(i)))});
    }
    size_t nQnew = params.size();
    auto qq = qq_in;
    auto qx = qx_in;
    auto xq = xq_in;
    if (nQnew > 1) {
      auto rp = molpro::linalg::itsolv::detail::redundant_parameters(Matrix<T>(qq.at(EqnData::S)), 0, nQnew,
                                                                     precision_scaled<T>(1e-8), *m_logger);
      // New parameters occupy the first nQnew positions of m_params, in the same order as the rows/cols of qq.
      // Remove the redundant ones from the highest index down so that the remaining indices stay valid.
      std::sort(rp.begin(), rp.end(), std::greater<>());
      for (auto i : rp) {
        m_params.erase(std::next(m_params.begin(), i));
        for (auto d : square_blocks(qq)) {
          qq.at(d).remove_row_col(i, i);
          qx.at(d).remove_row(i);
          xq.at(d).remove_col(i);
        }
        if (!qq.at(EqnData::rhs).empty())
          qq.at(EqnData::rhs).remove_row(i);
      }
      nQnew -= rp.size();
    }
    const auto nXnew = dims.nX() + nQnew;
    auto data = old_data;
    for (auto d : square_blocks(qq)) {
      data[d].resize({dims.nX() + nQnew, dims.nX() + nQnew});
      data[d].slice({dims.oQ() + nQnew, dims.oQ() + nQnew}, {data[d].rows(), data[d].cols()}) =
          old_data[d].slice({dims.oQ(), dims.oQ()}, {dims.nX(), dims.nX()});
      data[d].slice({dims.oQ(), dims.oQ()}, {dims.oQ() + nQnew, dims.oQ() + nQnew}) =
          qq.at(d).slice({0, 0}, {nQnew, nQnew});
      data[d].slice({dims.oQ(), 0}, {dims.oQ() + nQnew, dims.oQ()}) = qx.at(d).slice({0, 0}, {nQnew, dims.oQ()});
      data[d].slice({dims.oQ(), dims.oQ() + nQnew}, {dims.oQ() + nQnew, nXnew}) =
          qx.at(d).slice({0, dims.oQ()}, {nQnew, dims.nX()});
      data[d].slice({0, dims.oQ()}, {dims.oQ(), dims.oQ() + nQnew}) = xq.at(d).slice({0, 0}, {dims.oQ(), nQnew});
      data[d].slice({dims.oQ() + nQnew, dims.oQ()}, {nXnew, dims.oQ() + nQnew}) =
          xq.at(d).slice({dims.oQ(), 0}, {dims.nX(), nQnew});
      data[d].slice({0, 0}, {dims.oQ(), dims.oQ()}) = old_data[d].slice({0, 0}, {dims.oQ(), dims.oQ()});
      data[d].slice({0, dims.oQ() + nQnew}, {dims.oQ(), nXnew}) =
          old_data[d].slice({0, dims.oQ()}, {dims.oQ(), dims.nX()});
      data[d].slice({dims.oQ() + nQnew, 0}, {nXnew, dims.oQ()}) =
          old_data[d].slice({dims.oQ(), 0}, {dims.nX(), dims.oQ()});
      old_data[d] = data[d];
    }
    if (!qq.at(EqnData::rhs).empty()) {
      data[EqnData::rhs].resize({dims.nX() + nQnew, dims.nRHS});
      data[EqnData::rhs].slice({dims.oQ(), 0}, {dims.oQ() + nQnew, dims.nRHS}) = qq.at(EqnData::rhs).slice();
      data[EqnData::rhs].slice({dims.oQ() + nQnew, 0}, {dims.nX() + nQnew, dims.nRHS}) =
          old_data[EqnData::rhs].slice({dims.oQ(), 0}, {dims.nX(), dims.nRHS});
      old_data[EqnData::rhs] = data[EqnData::rhs];
    }
    m_logger->data_dump("S = ", data.at(EqnData::S));
    m_logger->data_dump("H = ", data.at(EqnData::H));
    if (data.count(EqnData::M))
      m_logger->data_dump("M = ", data.at(EqnData::M));
    m_logger->data_dump("rhs = ", data.at(EqnData::rhs));
  }
 
  void clear() { m_params.clear(); }
 
  void erase(size_t i) {
    assert(m_params.size() > i);
    auto it = std::next(begin(m_params), i);
    m_params.erase(it);
  }
 
  size_t size() const { return m_params.size(); }
 
  VecRef<Q> params() {
    auto qparams = qspace::wrap_params<Q>(m_params.begin(), m_params.end());
    return qparams[0];
  }
 
  CVecRef<Q> params() const {
    auto qparams = qspace::cwrap_params<Q>(m_params.begin(), m_params.end());
    return qparams[0];
  }
 
  CVecRef<Q> cparams() const { return params(); }
 
  VecRef<Q> actions() {
    auto qparams = qspace::wrap_params<Q>(m_params.begin(), m_params.end());
    return qparams[1];
  }
 
  CVecRef<Q> actions() const {
    auto qparams = qspace::cwrap_params<Q>(m_params.begin(), m_params.end());
    return qparams[1];
  }
 
  CVecRef<Q> cactions() const { return actions(); }
 
  CVecRef<Q> cmetric_actions() const {
    auto metric_actions = CVecRef<Q>{};
    for (const auto& p : m_params)
      if (p.metric_action)
        metric_actions.emplace_back(*p.metric_action);
    return metric_actions;
  }
 
protected:
  std::shared_ptr<ArrayHandlers<R, Q, P>> m_handlers;
  std::shared_ptr<Logger> m_logger;
  size_t m_unique_id{0};                 
  std::list<qspace::QParam<Q>> m_params; 
};
 
} // namespace molpro::linalg::itsolv::subspace
 
#endif // LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_SUBSPACE_QSPACE_H