IterativeSolver.h

#ifndef LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_ITERATIVESOLVER_H
#define LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_ITERATIVESOLVER_H
#include <molpro/linalg/array/ArrayHandler.h>
#include <molpro/linalg/array/Span.h>
#include <molpro/linalg/itsolv/Options.h>
#include <molpro/linalg/itsolv/Statistics.h>
#include <molpro/linalg/itsolv/subspace/Dimensions.h>
#include <molpro/linalg/itsolv/wrap.h>
 
#include <molpro/linalg/array/DistrArray.h>
#include <molpro/linalg/array/util/Distribution.h>
 
#include <memory>
#include <ostream>
#include <vector>
 
namespace molpro::profiler {
class Profiler;
}
namespace molpro::linalg::itsolv {
 
template <typename T>
struct has_iterator {
  template <typename C>
  constexpr static std::true_type test(typename C::iterator*);
 
  template <typename>
  constexpr static std::false_type test(...);
 
  constexpr static bool value =
      std::is_same<std::true_type, decltype(test<typename std::remove_reference<T>::type>(0))>::value;
};
 
template <typename T, typename = std::enable_if_t<std::is_base_of<molpro::linalg::array::DistrArray, T>::value>>
void precondition_default(const VecRef<T>& action, const std::vector<double>& shift, const T& diagonals) {
  auto diagonals_local_buffer = diagonals.local_buffer();
  for (size_t k = 0; k < action.size(); k++) {
    auto action_local_buffer = action[k].get().local_buffer();
    auto& distribution = action[k].get().distribution();
    auto range = distribution.range(molpro::mpi::rank_global());
    for (auto i = range.first; i < range.second; i++)
      (*action_local_buffer)[i - range.first] /= ((*diagonals_local_buffer)[i - range.first] - shift[k] + 1e-15);
  }
}
 
template <class T>
void precondition_default(const VecRef<T>& action, const std::vector<double>& shift, const T& diagonals,
                          typename T::iterator* = nullptr // SFINAE
) {
  for (size_t k = 0; k < action.size(); k++) {
    auto& a = action[k].get();
    std::transform(diagonals.begin(), diagonals.end(), a.begin(), a.begin(),
                   [shift, k](const auto& first, const auto& second) { return second / (first - shift[k] + 1e-15); });
  }
}
 
template <typename T, typename = std::enable_if_t<!std::is_base_of<molpro::linalg::array::DistrArray, T>::value>,
          class = void>
void precondition_default(const VecRef<T>& action, const std::vector<double>& shift, const T& diagonals,
                          typename std::enable_if<!has_iterator<T>::value, void*>::type = nullptr // SFINAE
) {
  throw std::logic_error("Unimplemented preconditioner");
}
 
template <typename R, typename P = std::map<size_t, typename R::value_type>>
class Problem {
public:
  Problem() = default;
  virtual ~Problem() = default;
  using container_t = R;
  using value_t = typename R::value_type;
 
  virtual value_t residual(const R& parameters, R& residual) const { return 0; }
 
  virtual void action(const CVecRef<R>& parameters, const VecRef<R>& action) const { return; }
 
  virtual bool diagonals(container_t& d) const { return false; }
 
  virtual void precondition(const VecRef<R>& residual, const std::vector<value_t>& shift) const { return; }
 
  virtual void precondition(const VecRef<R>& residual, const std::vector<value_t>& shift, const R& diagonals) const {
    precondition_default(residual, shift, diagonals);
  }
 
  virtual bool RHS(R& RHS, unsigned int instance) const { return false;}
 
  virtual std::vector<double> pp_action_matrix(const std::vector<P>& pparams) const {
    if (not pparams.empty())
      throw std::logic_error("P-space unavailable: unimplemented pp_action_matrix() in Problem class");
    return std::vector<double>(0);
  }
 
  virtual void p_action(const std::vector<std::vector<value_t>>& p_coefficients, const CVecRef<P>& pparams,
                        const VecRef<container_t>& actions) const {
    if (not pparams.empty())
      throw std::logic_error("P-space unavailable: unimplemented p_action() in Problem class");
  }
 
  virtual bool test_parameters(unsigned int instance, R& parameters) const { return false; }
};
 
template <class R, class Q, class P>
class IterativeSolver {
public:
  using value_type = typename R::value_type;                          
  using scalar_type = typename array::ArrayHandler<R, Q>::value_type; 
  using value_type_abs = typename array::ArrayHandler<R, R>::value_type_abs;
  using VectorP = std::vector<value_type>; 
  using fapply_on_p_type = std::function<void(const std::vector<VectorP>&, const CVecRef<P>&, const VecRef<R>&)>;
 
  virtual ~IterativeSolver() = default;
  IterativeSolver() = default;
  IterativeSolver(const IterativeSolver<R, Q, P>&) = delete;
  IterativeSolver<R, Q, P>& operator=(const IterativeSolver<R, Q, P>&) = delete;
  IterativeSolver(IterativeSolver<R, Q, P>&&) noexcept = default;
  IterativeSolver<R, Q, P>& operator=(IterativeSolver<R, Q, P>&&) noexcept = default;
 
  virtual bool solve(const VecRef<R>& parameters, const VecRef<R>& actions, const Problem<R>& problem,
                     bool generate_initial_guess = false) = 0;
  virtual bool solve(R& parameters, R& actions, const Problem<R>& problem, bool generate_initial_guess = false) = 0;
  virtual bool solve(std::vector<R>& parameters, std::vector<R>& actions, const Problem<R>& problem,
                     bool generate_initial_guess = false) = 0;
 
  virtual int add_vector(const VecRef<R>& parameters, const VecRef<R>& actions) = 0;
 
  // FIXME this should be removed in favour of VecRef interface
  virtual int add_vector(std::vector<R>& parameters, std::vector<R>& action) = 0;
  virtual int add_vector(R& parameters, R& action, value_type value = 0) = 0;
 
  virtual size_t add_p(const CVecRef<P>& pparams, const array::Span<value_type>& pp_action_matrix,
                       const VecRef<R>& parameters, const VecRef<R>& action, fapply_on_p_type apply_p) = 0;
 
  // FIXME Is this needed?
  virtual void clearP() = 0;
 
  virtual void solution(const std::vector<int>& roots, const VecRef<R>& parameters, const VecRef<R>& residual) = 0;
 
  virtual void solution_params(const std::vector<int>& roots, const VecRef<R>& parameters) = 0;
 
  virtual size_t end_iteration(const VecRef<R>& parameters, const VecRef<R>& residual) = 0;
 
  virtual bool end_iteration_needed() = 0;
 
  virtual std::vector<size_t> suggest_p(const CVecRef<R>& solution, const CVecRef<R>& residual, size_t max_number,
                                        double threshold) = 0;
 
  virtual void solution(const std::vector<int>& roots, std::vector<R>& parameters, std::vector<R>& residual) = 0;
  virtual void solution(R& parameters, R& residual) = 0;
  virtual void solution_params(const std::vector<int>& roots, std::vector<R>& parameters) = 0;
  virtual void solution_params(R& parameters) = 0;
  virtual size_t end_iteration(std::vector<R>& parameters, std::vector<R>& action) = 0;
  virtual size_t end_iteration(R& parameters, R& action) = 0;
 
  virtual const std::vector<int>& working_set() const = 0;
  virtual std::vector<scalar_type> working_set_eigenvalues() const {
    return std::vector<scalar_type>(working_set().size(), 0);
  }
  virtual size_t n_roots() const = 0;
  virtual void set_n_roots(size_t nroots) = 0;
  virtual const std::vector<scalar_type>& errors() const = 0;
  virtual const Statistics& statistics() const = 0;
  virtual void report(std::ostream& cout, bool endl = true) const = 0;
  virtual void report() const = 0;
 
  virtual void set_convergence_threshold(double thresh) = 0;
  virtual double convergence_threshold() const = 0;
  virtual void set_convergence_threshold_value(double thresh) = 0;
  virtual double convergence_threshold_value() const = 0;
  virtual void set_verbosity(Verbosity v) = 0;
  virtual void set_verbosity(int v) = 0;
  virtual Verbosity get_verbosity() const = 0;
  virtual void set_max_iter(int n) = 0;
  virtual int get_max_iter() const = 0;
  virtual void set_max_p(int n) = 0;
  virtual int get_max_p() const = 0;
  virtual void set_p_threshold(double thresh) = 0;
  virtual double get_p_threshold() const = 0;
  virtual const subspace::Dimensions& dimensions() const = 0;
  // FIXME Missing parameters: SVD threshold
  virtual void set_options(const Options& options) = 0;
  virtual std::shared_ptr<Options> get_options() const = 0;
  virtual scalar_type value() const = 0;
  virtual bool nonlinear() const = 0;
  virtual void set_profiler(molpro::profiler::Profiler& profiler) = 0;
  virtual const std::shared_ptr<molpro::profiler::Profiler>& profiler() const = 0;
  virtual bool test_problem(const Problem<R>& problem, R& v0, R& v1, int verbosity = 0,
                            double threshold = 1e-5) const = 0;
};
 
template <class R, class Q, class P>
class LinearEigensystem : public IterativeSolver<R, Q, P> {
public:
  using typename IterativeSolver<R, Q, P>::scalar_type;
  virtual std::vector<scalar_type> eigenvalues() const = 0;
  virtual void set_hermiticity(bool hermitian) = 0;
  virtual bool get_hermiticity() const = 0;
};
 
template <class R, class Q, class P>
class LinearEquations : public IterativeSolver<R, Q, P> {
public:
  using typename IterativeSolver<R, Q, P>::scalar_type;
  virtual void add_equations(const CVecRef<R>& rhs) = 0;
  virtual void add_equations(const std::vector<R>& rhs) = 0;
  virtual void add_equations(const R& rhs) = 0;
  virtual CVecRef<Q> rhs() const = 0;
  virtual void set_hermiticity(bool hermitian) = 0;
  virtual bool get_hermiticity() const = 0;
};
 
template <class R, class Q, class P>
class Optimize : public IterativeSolver<R, Q, P> {};
 
template <class R, class Q, class P>
class NonLinearEquations : public IterativeSolver<R, Q, P> {};
 
} // namespace molpro::linalg::itsolv
 
#endif // LINEARALGEBRA_SRC_MOLPRO_LINALG_ITSOLV_ITERATIVESOLVER_H