ExplicitInverse.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <algorithm>
7 #include <cstddef>
8 #include <vector>
9 
12 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
13 #include "NumericalAlgorithms/LinearSolver/LinearSolver.hpp"
14 #include "Options/Options.hpp"
17 #include "Utilities/TMPL.hpp"
18 #include "Utilities/TypeTraits/CreateIsCallable.hpp"
19 
20 namespace LinearSolver::Serial {
21 
22 /// \cond
23 template <typename LinearSolverRegistrars>
24 struct ExplicitInverse;
25 /// \endcond
26 
27 namespace Registrars {
28 /// Registers the `LinearSolver::Serial::ExplicitInverse` linear solver
30 } // namespace Registrars
31 
32 namespace detail {
33 CREATE_IS_CALLABLE(reset)
34 CREATE_IS_CALLABLE_V(reset)
35 } // namespace detail
36 
37 /*!
38  * \brief Linear solver that builds a matrix representation of the linear
39  * operator and inverts it directly
40  *
41  * This solver first constructs an explicit matrix representation by "sniffing
42  * out" the operator, i.e. feeding it with unit vectors, and then directly
43  * inverts the matrix. The result is an operator that solves the linear problem
44  * in a single step. This means that each element has a large initialization
45  * cost, but all successive solves converge immediately.
46  *
47  * \par Advice on using this linear solver:
48  *
49  * - This solver is entirely agnostic to the structure of the linear operator.
50  * It is usually better to implement a linear solver that is specialized for
51  * your linear operator to take advantage of its properties. For example, if
52  * the operator has a tensor-product structure, the linear solver might take
53  * advantage of that. Only use this solver if no alternatives are available
54  * and if you have verified that it speeds up your solves.
55  * - Since this linear solver stores the full inverse operator matrix it can
56  * have significant memory demands. For example, an operator representing a 3D
57  * first-order Elasticity system (9 variables) discretized on 12 grid points
58  * per dimension requires ca. 2GB of memory (per element) to store the matrix,
59  * scaling quadratically with the number of variables and with a power of 6
60  * with the number of grid points per dimension. Therefore, make sure to
61  * distribute the elements on a sufficient number of nodes to meet the memory
62  * requirements.
63  * - This linear solver can be `reset()` when the operator changes (e.g. in each
64  * nonlinear-solver iteration). However, when using this solver as
65  * preconditioner it can be advantageous to avoid the reset and the
66  * corresponding cost of re-building the matrix and its inverse if the
67  * operator only changes "a little". In that case the preconditioner solves
68  * subdomain problems only approximately, but possibly still sufficiently to
69  * provide effective preconditioning. You can toggle this behaviour with the
70  * `EnableResets` option.
71  */
72 template <typename LinearSolverRegistrars =
73  tmpl::list<Registrars::ExplicitInverse>>
74 class ExplicitInverse : public LinearSolver<LinearSolverRegistrars> {
75  private:
77 
78  public:
79  struct EnableResets {
80  using type = bool;
81  static constexpr Options::String help =
82  "Enable or disable resets. This only has an effect in cases where the "
83  "operator changes, e.g. between nonlinear-solver iterations. Only "
84  "disable resets when using this solver as preconditioner and thus "
85  "approximate solutions are desirable. Disabling resets avoids "
86  "expensive re-building of the operator, but comes at the cost of less "
87  "accurate preconditioning and thus potentially more preconditioned "
88  "iterations. Whether or not this helps convergence overall is highly "
89  "problem-dependent.";
90  };
91 
92  using options = tmpl::list<EnableResets>;
93  static constexpr Options::String help =
94  "Build a matrix representation of the linear operator and invert it "
95  "directly. This means that the first solve has a large initialization "
96  "cost, but all subsequent solves converge immediately.";
97 
98  ExplicitInverse() = default;
99  ExplicitInverse(const ExplicitInverse& /*rhs*/) = default;
100  ExplicitInverse& operator=(const ExplicitInverse& /*rhs*/) = default;
101  ExplicitInverse(ExplicitInverse&& /*rhs*/) = default;
102  ExplicitInverse& operator=(ExplicitInverse&& /*rhs*/) = default;
103  ~ExplicitInverse() = default;
104 
105  explicit ExplicitInverse(bool enable_resets) noexcept
106  : enable_resets_(enable_resets) {}
107 
108  /// \cond
109  explicit ExplicitInverse(CkMigrateMessage* m) noexcept : Base(m) {}
110  using PUP::able::register_constructor;
111  WRAPPED_PUPable_decl_template(ExplicitInverse); // NOLINT
112  /// \endcond
113 
114  /*!
115  * \brief Solve the equation \f$Ax=b\f$ by explicitly constructing the
116  * operator matrix \f$A\f$ and its inverse. The first solve is computationally
117  * expensive and successive solves are cheap.
118  *
119  * Building a matrix representation of the `linear_operator` requires
120  * iterating over the `SourceType` (which is also the type returned by the
121  * `linear_operator`) in a consistent way. This can be non-trivial for
122  * heterogeneous data structures because it requires they define a data
123  * ordering. Specifically, the `SourceType` must have a `size()` function as
124  * well as `begin()` and `end()` iterators that point into the data. If the
125  * iterators have a `reset()` function it is used to avoid repeatedly
126  * re-creating the `begin()` iterator. The `reset()` function must not
127  * invalidate the `end()` iterator.
128  */
129  template <typename LinearOperator, typename VarsType, typename SourceType>
131  const LinearOperator& linear_operator,
132  const SourceType& source) const noexcept;
133 
134  /// Flags the operator to require re-initialization. No memory is released.
135  /// Call this function to rebuild the solver when the operator changed.
136  void reset() noexcept override {
137  if (not enable_resets_) {
138  return;
139  }
141  }
142 
143  /// Size of the operator. The stored matrix will have `size^2` entries.
144  size_t size() const noexcept { return size_; }
145 
146  /// The matrix representation of the solver. This matrix approximates the
147  /// inverse of the subdomain operator.
148  const DenseMatrix<double>& matrix_representation() const noexcept {
149  return inverse_;
150  }
151 
152  // NOLINTNEXTLINE(google-runtime-references)
153  void pup(PUP::er& p) noexcept override {
154  p | enable_resets_;
155  p | size_;
156  p | inverse_;
157  if (p.isUnpacking() and size_ != std::numeric_limits<size_t>::max()) {
158  source_workspace_.resize(size_);
159  solution_workspace_.resize(size_);
160  }
161  }
162 
163  std::unique_ptr<Base> get_clone() const noexcept override {
164  return std::make_unique<ExplicitInverse>(*this);
165  }
166 
167  private:
168  bool enable_resets_{};
169 
170  // Caches for successive solves of the same operator
171  mutable size_t size_ = std::numeric_limits<size_t>::max();
172  mutable DenseMatrix<double, blaze::columnMajor> inverse_{};
173 
174  // Buffers to avoid re-allocating memory for applying the operator
175  mutable DenseVector<double> source_workspace_{};
176  mutable DenseVector<double> solution_workspace_{};
177 };
178 
179 template <typename LinearSolverRegistrars>
180 template <typename LinearOperator, typename VarsType, typename SourceType>
182  const gsl::not_null<VarsType*> solution,
183  const LinearOperator& linear_operator,
184  const SourceType& source) const noexcept {
185  if (UNLIKELY(size_ == std::numeric_limits<size_t>::max())) {
186  const auto& used_for_size = source;
187  size_ = used_for_size.size();
188  source_workspace_.resize(size_);
189  solution_workspace_.resize(size_);
190  inverse_.resize(size_, size_);
191  // Construct explicit matrix representation by "sniffing out" the operator,
192  // i.e. feeding it unit vectors
193  auto unit_vector = make_with_value<VarsType>(used_for_size, 0.);
194  auto result_buffer = make_with_value<SourceType>(used_for_size, 0.);
195  auto& operator_matrix = inverse_;
196  size_t i = 0;
197  // Re-using the iterators for all operator invocations
198  auto result_iterator_begin = result_buffer.begin();
199  auto result_iterator_end = result_buffer.end();
200  for (double& unit_vector_data : unit_vector) {
201  // Add a 1 at the unit vector location i
202  unit_vector_data = 1.;
203  // Invoke the operator on the unit vector
204  linear_operator(make_not_null(&result_buffer), unit_vector);
205  // Set the unit vector back to zero
206  unit_vector_data = 0.;
207  // Reset the iterator by calling its `reset` member function or by
208  // re-creating it
209  if constexpr (detail::is_reset_callable_v<decltype(
210  result_iterator_begin)>) {
211  result_iterator_begin.reset();
212  } else {
213  result_iterator_begin = result_buffer.begin();
214  result_iterator_end = result_buffer.end();
215  }
216  // Store the result in column i of the matrix
217  std::copy(result_iterator_begin, result_iterator_end,
218  column(operator_matrix, i).begin());
219  ++i;
220  }
221  // Directly invert the matrix
222  try {
223  blaze::invert(inverse_);
224  } catch (const std::invalid_argument& e) {
225  ERROR("Could not invert subdomain matrix (size " << size_
226  << "): " << e.what());
227  }
228  }
229  // Copy source into contiguous workspace. In cases where the source and
230  // solution data are already stored contiguously we might avoid the copy and
231  // the associated workspace memory. However, compared to the cost of building
232  // and storing the matrix this is likely insignificant.
233  std::copy(source.begin(), source.end(), source_workspace_.begin());
234  // Apply inverse
235  solution_workspace_ = inverse_ * source_workspace_;
236  // Reconstruct solution data from contiguous workspace
237  std::copy(solution_workspace_.begin(), solution_workspace_.end(),
238  solution->begin());
239  return {0, 0};
240 }
241 
242 /// \cond
243 template <typename LinearSolverRegistrars>
244 // NOLINTNEXTLINE
245 PUP::able::PUP_ID ExplicitInverse<LinearSolverRegistrars>::my_PUP_ID = 0;
246 /// \endcond
247 
248 } // namespace LinearSolver::Serial
LinearSolver::Serial::LinearSolver
Base class for serial linear solvers that supports factory-creation.
Definition: LinearSolver.hpp:37
LinearSolver::Serial::ExplicitInverse::size
size_t size() const noexcept
Size of the operator. The stored matrix will have size^2 entries.
Definition: ExplicitInverse.hpp:144
CharmPupable.hpp
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
Options.hpp
vector
LinearSolver::Serial::ExplicitInverse::solve
Convergence::HasConverged solve(gsl::not_null< VarsType * > solution, const LinearOperator &linear_operator, const SourceType &source) const noexcept
Solve the equation by explicitly constructing the operator matrix and its inverse....
DenseVector.hpp
DenseMatrix.hpp
DenseVector< double >
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:68
LinearSolver::Serial::ExplicitInverse::reset
void reset() noexcept override
Flags the operator to require re-initialization. No memory is released. Call this function to rebuild...
Definition: ExplicitInverse.hpp:136
algorithm
Registration::Registrar
A template for defining a registrar.
Definition: Registration.hpp:42
LinearSolver::Serial::ExplicitInverse::EnableResets
Definition: ExplicitInverse.hpp:79
DenseMatrix< double >
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
WRAPPED_PUPable_decl_template
#define WRAPPED_PUPable_decl_template(className)
Mark derived classes as serializable.
Definition: CharmPupable.hpp:22
LinearSolver::Serial::ExplicitInverse::matrix_representation
const DenseMatrix< double > & matrix_representation() const noexcept
The matrix representation of the solver. This matrix approximates the inverse of the subdomain operat...
Definition: ExplicitInverse.hpp:148
cstddef
MakeWithValue.hpp
std::invalid_argument
LinearSolver::Serial::ExplicitInverse
Linear solver that builds a matrix representation of the linear operator and inverts it directly.
Definition: ExplicitInverse.hpp:74
LinearSolver
Functionality for solving linear systems of equations.
Definition: ExplicitInverse.hpp:20
std::begin
T begin(T... args)
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
std::numeric_limits::max
T max(T... args)
std::unique_ptr
std::numeric_limits
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
CREATE_IS_CALLABLE
#define CREATE_IS_CALLABLE(METHOD_NAME)
Generate a type trait to check if a class has a member function that can be invoked with arguments of...
Definition: CreateIsCallable.hpp:30