DavidsonCommon.h

#ifndef LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_DAVIDSONCOMMON_H
#define LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_DAVIDSONCOMMON_H
#include <molpro/linalg/itsolv/DSpaceResetter.h>
#include <molpro/linalg/itsolv/DavidsonOptions.h>
#include <molpro/linalg/itsolv/IterativeSolverTemplate.h>
#include <molpro/linalg/itsolv/Logger.h>
#include <molpro/linalg/itsolv/Options.h>
#include <molpro/linalg/itsolv/propose_rspace.h>
#include <molpro/linalg/itsolv/qspace_options.h>
#include <molpro/linalg/itsolv/rspace_options.h>
#include <molpro/linalg/itsolv/subspace/SubspaceSolverLinEig.h>
#include <molpro/linalg/itsolv/subspace/XSpace.h>
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <memory>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace molpro::linalg::itsolv {
 
template <template <class, class, class> class Solver, class R, class Q, class P>
class DavidsonCommon : public IterativeSolverTemplate<Solver, R, Q, P> {
public:
  using SolverTemplate = IterativeSolverTemplate<Solver, R, Q, P>;
  using typename SolverTemplate::value_type;
  using typename SolverTemplate::value_type_abs;
 
  bool nonlinear() const override { return false; }
 
  size_t end_iteration(const VecRef<R>& parameters, const VecRef<R>& action) override {
    auto prof = this->profiler();
    auto m_prof = prof->push("itsolv::end_iteration");
    if (this->m_working_set.empty() && !m_dspace_resetter.in_progress()) {
      // all roots have converged: the subspace is left as it is
      m_resetting_in_progress = false;
    } else if (m_dspace_resetter.do_reset(this->m_stats->iterations, this->m_xspace->dimensions())) {
      m_resetting_in_progress = true;
      // retained_solutions() can hold a truncated complex pair partner beyond the roots; with no more parameters than
      // roots, it is handed out under the index of a root
      const auto nparam = std::min(parameters.size(), this->n_roots());
      this->m_working_set =
          m_dspace_resetter.run(VecRef<R>(parameters.begin(), parameters.begin() + nparam), *this->m_xspace,
                                this->m_subspace_solver->retained_solutions(), rspace_opts.norm_thresh,
                                rspace_opts.svd_thresh, *this->m_handlers, *this->m_logger);
      this->m_solutions_outdated = true;
    } else {
      m_resetting_in_progress = false;
      prof->start("end_iteration (propose_rspace)");
      this->m_working_set =
          detail::propose_rspace(*this, parameters, action, *this->m_xspace, *this->m_subspace_solver,
                                 *this->m_handlers, *this->m_logger, rspace_opts, qspace_opts, *this->profiler().get());
      prof->stop();
    }
    this->m_stats->iterations++;
    read_handler_counts(this->m_stats, this->m_handlers);
    this->m_end_iteration_needed = false;
    return this->working_set().size();
  }
  size_t end_iteration(std::vector<R>& parameters, std::vector<R>& action) override {
    return end_iteration(wrap(parameters), wrap(action));
  }
  size_t end_iteration(R& parameters, R& actions) override {
    return end_iteration(wrap_arg(parameters), wrap_arg(actions));
  }
 
  void set_norm_thresh(value_type_abs thresh) { rspace_opts.norm_thresh = thresh; }
  value_type_abs get_norm_thresh() const { return rspace_opts.norm_thresh; }
  void set_svd_thresh(value_type_abs thresh) { rspace_opts.svd_thresh = thresh; }
  value_type_abs get_svd_thresh() const { return rspace_opts.svd_thresh; }
  void set_reset_D(size_t n) { m_dspace_resetter.set_nreset(n); }
  size_t get_reset_D() const { return m_dspace_resetter.get_nreset(); }
  void set_reset_D_maxQ_size(size_t n) { m_dspace_resetter.set_max_Qsize(n); }
  size_t get_reset_D_maxQ_size() const { return m_dspace_resetter.get_max_Qsize(); }
  void set_max_size_qspace(std::size_t n) {
    check_qspace_size_limits(get_min_size_qspace(), n);
    assign_max_size_qspace(n);
  }
  std::size_t get_max_size_qspace() const { return qspace_opts.max_size; }
  void set_min_size_qspace(std::size_t n) {
    check_qspace_size_limits(n, get_max_size_qspace());
    qspace_opts.min_size = n;
  }
  std::size_t get_min_size_qspace() const { return qspace_opts.min_size; }
 
  void set_hermiticity(bool hermitian) override {
    m_hermiticity = hermitian;
    auto xspace = std::dynamic_pointer_cast<subspace::XSpace<R, Q, P>>(this->m_xspace);
    xspace->set_hermiticity(hermitian);
    auto subspace_solver = std::dynamic_pointer_cast<subspace::SubspaceSolverLinEig<R, Q, P>>(this->m_subspace_solver);
    subspace_solver->set_hermiticity(hermitian);
  }
  bool get_hermiticity() const override { return m_hermiticity; }
 
protected:
  DavidsonCommon(const std::shared_ptr<ArrayHandlers<R, Q, P>>& handlers, const std::shared_ptr<Logger>& logger,
                 bool hermiticity)
      : SolverTemplate(std::make_shared<subspace::XSpace<R, Q, P>>(handlers, logger),
                       std::static_pointer_cast<subspace::ISubspaceSolver<R, Q, P>>(
                           std::make_shared<subspace::SubspaceSolverLinEig<R, Q, P>>(logger)),
                       handlers, std::make_shared<Statistics>(), logger) {
    set_hermiticity(hermiticity);
  }
 
  void set_davidson_options(const Options& options, const DavidsonOptions& opt) {
    const auto min_size_qspace = opt.min_size_qspace.value_or(get_min_size_qspace());
    const auto max_size_qspace = opt.max_size_qspace.value_or(get_max_size_qspace());
    check_qspace_size_limits(min_size_qspace, max_size_qspace);
    SolverTemplate::set_options(options);
    if (opt.reset_D)
      set_reset_D(opt.reset_D.value());
    if (opt.reset_D_max_Q_size)
      set_reset_D_maxQ_size(opt.reset_D_max_Q_size.value());
    if (opt.max_size_qspace)
      assign_max_size_qspace(max_size_qspace);
    if (opt.min_size_qspace)
      qspace_opts.min_size = min_size_qspace;
    if (opt.contrib_thresh)
      qspace_opts.contrib_thresh = opt.contrib_thresh.value();
    if (opt.norm_thresh)
      set_norm_thresh(opt.norm_thresh.value());
    if (opt.svd_thresh)
      set_svd_thresh(opt.svd_thresh.value());
    if (opt.hermiticity)
      set_hermiticity(opt.hermiticity.value());
  }
 
  template <class SolverOptions>
  std::shared_ptr<SolverOptions> get_davidson_options() const {
    auto opt = std::make_shared<SolverOptions>();
    opt->copy(*SolverTemplate::get_options());
    opt->reset_D = get_reset_D();
    opt->reset_D_max_Q_size = get_reset_D_maxQ_size();
    opt->max_size_qspace = get_max_size_qspace();
    opt->min_size_qspace = get_min_size_qspace();
    opt->contrib_thresh = qspace_opts.contrib_thresh;
    opt->norm_thresh = get_norm_thresh();
    opt->svd_thresh = get_svd_thresh();
    opt->hermiticity = get_hermiticity();
    return opt;
  }
 
  static void check_qspace_size_limits(std::size_t min_size, std::size_t max_size) {
    if (min_size > max_size)
      throw std::invalid_argument("minimum size of the Q space (" + std::to_string(min_size) +
                                  ") exceeds its maximum size (" + std::to_string(max_size) + ")");
  }
  void assign_max_size_qspace(std::size_t n) {
    qspace_opts.max_size = n;
    if (get_reset_D_maxQ_size() > n)
      set_reset_D_maxQ_size(n);
  }
 
  virtual std::vector<value_type> values() const = 0;
 
  void set_value_errors() override {
    const auto current_values = values();
    this->m_value_errors.assign(current_values.size(), std::numeric_limits<value_type_abs>::max());
    for (size_t i = 0; i < std::min(m_last_values.size(), current_values.size()); i++)
      this->m_value_errors[i] = std::abs(current_values[i] - m_last_values[i]);
    // resetting the D space does not change the solutions, so the change measured across it would be meaningless
    if (!m_resetting_in_progress)
      m_last_values = current_values;
  }
 
  RSpaceOptions<value_type_abs> rspace_opts;   
  QSpaceOptions qspace_opts;                   
  detail::DSpaceResetter<Q> m_dspace_resetter; 
  bool m_hermiticity = false;                  
  std::vector<value_type> m_last_values;       
  bool m_resetting_in_progress = false;        
};
 
} // namespace molpro::linalg::itsolv
 
#endif // LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_DAVIDSONCOMMON_H