LinearSolver.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <memory>
7 #include <optional>
8 #include <pup.h>
9 #include <pup_stl.h>
10 #include <type_traits>
11 #include <utility>
12 
13 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
14 #include "Options/Auto.hpp"
15 #include "Options/Options.hpp"
17 #include "Parallel/PupStlCpp17.hpp"
18 #include "Utilities/FakeVirtual.hpp"
19 #include "Utilities/Gsl.hpp"
20 #include "Utilities/Registration.hpp"
21 #include "Utilities/Requires.hpp"
22 
23 namespace LinearSolver::Serial {
24 
25 /// Registrars for linear solvers
26 namespace Registrars {}
27 
28 /*!
29  * \brief Base class for serial linear solvers that supports factory-creation.
30  *
31  * Derive linear solvers from this class so they can be factory-created. If your
32  * linear solver supports preconditioning, derive from
33  * `PreconditionedLinearSolver` instead to inherit utility that allows using any
34  * other factor-creatable linear solver as preconditioner.
35  */
36 template <typename LinearSolverRegistrars>
37 class LinearSolver : public PUP::able {
38  protected:
39  /// \cond
40  LinearSolver() = default;
41  LinearSolver(const LinearSolver&) = default;
42  LinearSolver(LinearSolver&&) = default;
43  LinearSolver& operator=(const LinearSolver&) = default;
44  LinearSolver& operator=(LinearSolver&&) = default;
45  /// \endcond
46 
47  public:
48  ~LinearSolver() override = default;
49 
50  /// \cond
51  explicit LinearSolver(CkMigrateMessage* m) noexcept;
53  /// \endcond
54 
55  using creatable_classes = Registration::registrants<LinearSolverRegistrars>;
56 
58  const = 0;
59 
60  /*!
61  * \brief Solve the linear equation \f$Ax=b\f$ where \f$A\f$ is the
62  * `linear_operator` and \f$b\f$ is the `source`.
63  *
64  * - The (approximate) solution \f$x\f$ is returned in the
65  * `initial_guess_in_solution_out` buffer, which also serves to provide an
66  * initial guess for \f$x\f$. Not all solvers take the initial guess into
67  * account, but all expect the buffer is sized correctly.
68  * - The `linear_operator` must be an invocable that takes a `VarsType` as
69  * const-ref argument and returns a `SourceType` by reference.
70  *
71  * Each solve may mutate the private state of the solver, for example to cache
72  * quantities to accelerate successive solves for the same operator. Invoke
73  * `reset` to discard these caches.
74  */
75  template <typename LinearOperator, typename VarsType, typename SourceType,
76  typename... Args>
78  gsl::not_null<VarsType*> initial_guess_in_solution_out,
79  const LinearOperator& linear_operator, const SourceType& source,
80  Args&&... args) const noexcept;
81 
82  /// Discard caches from previous solves. Use before solving a different linear
83  /// operator.
84  virtual void reset() noexcept = 0;
85 };
86 
87 /// \cond
88 template <typename LinearSolverRegistrars>
89 LinearSolver<LinearSolverRegistrars>::LinearSolver(CkMigrateMessage* m) noexcept
90  : PUP::able(m) {}
91 /// \endcond
92 
93 template <typename LinearSolverRegistrars>
94 template <typename LinearOperator, typename VarsType, typename SourceType,
95  typename... Args>
97  const gsl::not_null<VarsType*> initial_guess_in_solution_out,
98  const LinearOperator& linear_operator, const SourceType& source,
99  Args&&... args) const noexcept {
100  return call_with_dynamic_type<Convergence::HasConverged, creatable_classes>(
101  this, [&initial_guess_in_solution_out, &linear_operator, &source,
102  &args...](auto* const linear_solver) noexcept {
103  return linear_solver->solve(initial_guess_in_solution_out,
104  linear_operator, source,
105  std::forward<Args>(args)...);
106  });
107 }
108 
109 /// Indicates the linear solver uses no preconditioner. It may perform
110 /// compile-time optimization for this case.
112 
113 /*!
114  * \brief Base class for serial linear solvers that supports factory-creation
115  * and nested preconditioning.
116  *
117  * To enable support for preconditioning in your derived linear solver class,
118  * pass any type that has a `solve` and a `reset` function as the
119  * `Preconditioner` template parameter. It can also be an abstract
120  * `LinearSolver` type, which means that any other linear solver can be used as
121  * preconditioner. Pass `NoPreconditioner` to disable support for
122  * preconditioning.
123  */
124 template <typename Preconditioner, typename LinearSolverRegistrars>
125 class PreconditionedLinearSolver : public LinearSolver<LinearSolverRegistrars> {
126  private:
128 
129  public:
130  using PreconditionerType =
131  tmpl::conditional_t<std::is_abstract_v<Preconditioner>,
132  std::unique_ptr<Preconditioner>, Preconditioner>;
134  static std::string name() noexcept { return "Preconditioner"; }
135  // Support factory-creatable preconditioners by storing them as unique-ptrs
137  static constexpr Options::String help =
138  "An approximate linear solve in every iteration that helps the "
139  "algorithm converge.";
140  };
141 
142  protected:
143  PreconditionedLinearSolver() = default;
146 
148  std::optional<PreconditionerType> local_preconditioner) noexcept;
149 
151  PreconditionedLinearSolver& operator=(
152  const PreconditionedLinearSolver& rhs) noexcept;
153 
154  public:
155  ~PreconditionedLinearSolver() override = default;
156 
157  /// \cond
158  explicit PreconditionedLinearSolver(CkMigrateMessage* m) noexcept;
159  /// \endcond
160 
161  void pup(PUP::er& p) noexcept override { // NOLINT
162  PUP::able::pup(p);
163  if constexpr (not std::is_same_v<Preconditioner, NoPreconditioner>) {
164  p | preconditioner_;
165  }
166  }
167 
168  /// Whether or not a preconditioner is set
169  bool has_preconditioner() const noexcept {
170  if constexpr (not std::is_same_v<Preconditioner, NoPreconditioner>) {
171  return preconditioner_.has_value();
172  } else {
173  return false;
174  }
175  }
176 
177  // @{
178  /// Access to the preconditioner. Check `has_preconditioner()` before calling
179  /// this function. Calling this function when `has_preconditioner()` returns
180  /// `false` is an error.
181  template <
182  bool Enabled = not std::is_same_v<Preconditioner, NoPreconditioner>,
183  Requires<Enabled and
184  not std::is_same_v<Preconditioner, NoPreconditioner>> = nullptr>
185  const Preconditioner& preconditioner() const noexcept {
187  "No preconditioner is set. Please use `has_preconditioner()` to "
188  "check before trying to retrieve it.");
189  if constexpr (std::is_abstract_v<Preconditioner>) {
190  return **preconditioner_;
191  } else {
192  return *preconditioner_;
193  }
194  }
195 
196  template <
197  bool Enabled = not std::is_same_v<Preconditioner, NoPreconditioner>,
198  Requires<Enabled and
199  not std::is_same_v<Preconditioner, NoPreconditioner>> = nullptr>
200  Preconditioner& preconditioner() noexcept {
202  "No preconditioner is set. Please use `has_preconditioner()` to "
203  "check before trying to retrieve it.");
204  if constexpr (std::is_abstract_v<Preconditioner>) {
205  return **preconditioner_;
206  } else {
207  return *preconditioner_;
208  }
209  }
210 
211  // Keep the function virtual so derived classes must provide an
212  // implementation, but also provide an implementation below that derived
213  // classes can use to reset the preconditioner
214  void reset() noexcept override = 0;
215 
216  protected:
217  /// Copy the preconditioner. Useful to implement `get_clone` when the
218  /// preconditioner has an abstract type.
219  template <
220  bool Enabled = not std::is_same_v<Preconditioner, NoPreconditioner>,
221  Requires<Enabled and
222  not std::is_same_v<Preconditioner, NoPreconditioner>> = nullptr>
223  std::optional<PreconditionerType> clone_preconditioner() const noexcept {
224  if constexpr (std::is_abstract_v<Preconditioner>) {
225  return has_preconditioner()
226  ? std::optional((*preconditioner_)->get_clone())
227  : std::nullopt;
228  } else {
229  return preconditioner_;
230  }
231  }
232 
233  private:
234  // Only needed when preconditioning is enabled, but current C++ can't remove
235  // this variable at compile-time. Keeping the variable shouldn't have any
236  // noticeable overhead though.
237  std::optional<PreconditionerType> preconditioner_{};
238 };
239 
240 template <typename Preconditioner, typename LinearSolverRegistrars>
241 PreconditionedLinearSolver<Preconditioner, LinearSolverRegistrars>::
242  PreconditionedLinearSolver(
243  std::optional<PreconditionerType> local_preconditioner) noexcept
244  : preconditioner_(std::move(local_preconditioner)) {}
245 
246 // Override copy constructors so they can clone abstract preconditioners
247 template <typename Preconditioner, typename LinearSolverRegistrars>
248 PreconditionedLinearSolver<Preconditioner, LinearSolverRegistrars>::
249  PreconditionedLinearSolver(const PreconditionedLinearSolver& rhs) noexcept
250  : Base(rhs) {
251  if constexpr (not std::is_same_v<Preconditioner, NoPreconditioner>) {
252  preconditioner_ = rhs.clone_preconditioner();
253  }
254 }
255 template <typename Preconditioner, typename LinearSolverRegistrars>
256 PreconditionedLinearSolver<Preconditioner, LinearSolverRegistrars>&
257 PreconditionedLinearSolver<Preconditioner, LinearSolverRegistrars>::operator=(
258  const PreconditionedLinearSolver& rhs) noexcept {
259  Base::operator=(rhs);
260  if constexpr (not std::is_same_v<Preconditioner, NoPreconditioner>) {
261  preconditioner_ = rhs.clone_preconditioner();
262  }
263  return *this;
264 }
265 
266 /// \cond
267 template <typename Preconditioner, typename LinearSolverRegistrars>
268 PreconditionedLinearSolver<Preconditioner, LinearSolverRegistrars>::
269  PreconditionedLinearSolver(CkMigrateMessage* m) noexcept
270  : Base(m) {}
271 /// \endcond
272 
273 template <typename Preconditioner, typename LinearSolverRegistrars>
274 void PreconditionedLinearSolver<Preconditioner,
275  LinearSolverRegistrars>::reset() noexcept {
276  if constexpr (not std::is_same_v<Preconditioner, NoPreconditioner>) {
277  if (has_preconditioner()) {
278  preconditioner().reset();
279  }
280  }
281 }
282 
283 } // namespace LinearSolver::Serial
LinearSolver::Serial::LinearSolver
Base class for serial linear solvers that supports factory-creation.
Definition: LinearSolver.hpp:37
LinearSolver::Serial::NoPreconditioner
Indicates the linear solver uses no preconditioner. It may perform compile-time optimization for this...
Definition: LinearSolver.hpp:111
std::string
CharmPupable.hpp
utility
LinearSolver::Serial::PreconditionedLinearSolver::reset
void reset() noexcept override=0
Discard caches from previous solves. Use before solving a different linear operator.
Definition: LinearSolver.hpp:275
Options.hpp
LinearSolver::Serial::PreconditionedLinearSolver::PreconditionerOption
Definition: LinearSolver.hpp:133
LinearSolver::Serial::PreconditionedLinearSolver
Base class for serial linear solvers that supports factory-creation and nested preconditioning.
Definition: LinearSolver.hpp:125
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:68
LinearSolver::Serial::PreconditionedLinearSolver::has_preconditioner
bool has_preconditioner() const noexcept
Whether or not a preconditioner is set.
Definition: LinearSolver.hpp:169
LinearSolver::Serial::LinearSolver::reset
virtual void reset() noexcept=0
Discard caches from previous solves. Use before solving a different linear operator.
WRAPPED_PUPable_abstract
#define WRAPPED_PUPable_abstract(className)
Wraps the Charm++ macro, see the Charm++ documentation.
Definition: CharmPupable.hpp:41
memory
Options::Auto
A class indicating that a parsed value can be automatically computed instead of specified.
Definition: Auto.hpp:32
LinearSolver::Serial::LinearSolver::solve
Convergence::HasConverged solve(gsl::not_null< VarsType * > initial_guess_in_solution_out, const LinearOperator &linear_operator, const SourceType &source, Args &&... args) const noexcept
Solve the linear equation where is the linear_operator and is the source.
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
LinearSolver
Functionality for solving linear systems of equations.
Definition: Gmres.cpp:16
Gsl.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
LinearSolver::Serial::PreconditionedLinearSolver::clone_preconditioner
std::optional< PreconditionerType > clone_preconditioner() const noexcept
Copy the preconditioner. Useful to implement get_clone when the preconditioner has an abstract type.
Definition: LinearSolver.hpp:223
Requires.hpp
LinearSolver::Serial::PreconditionedLinearSolver::preconditioner
const Preconditioner & preconditioner() const noexcept
Access to the preconditioner. Check has_preconditioner() before calling this function....
Definition: LinearSolver.hpp:185
optional
Requires
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t
Definition: Requires.hpp:67
std::unique_ptr
Registration::registrants
tmpl::transform< tmpl::remove_duplicates< RegistrarList >, detail::registrant< tmpl::pin< tmpl::remove_duplicates< RegistrarList > >, tmpl::_1 > > registrants
Transform a list of registrars into the list of associated registrants. This is usually used to defin...
Definition: Registration.hpp:65
type_traits
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13