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)
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.
70  */
71 template <typename LinearSolverRegistrars =
72  tmpl::list<Registrars::ExplicitInverse>>
73 class ExplicitInverse : public LinearSolver<LinearSolverRegistrars> {
74  private:
76 
77  public:
78  using options = tmpl::list<>;
79  static constexpr Options::String help =
80  "Build a matrix representation of the linear operator and invert it "
81  "directly. This means that the first solve has a large initialization "
82  "cost, but all subsequent solves converge immediately.";
83 
84  ExplicitInverse() = default;
85  ExplicitInverse(const ExplicitInverse& /*rhs*/) = default;
86  ExplicitInverse& operator=(const ExplicitInverse& /*rhs*/) = default;
87  ExplicitInverse(ExplicitInverse&& /*rhs*/) = default;
88  ExplicitInverse& operator=(ExplicitInverse&& /*rhs*/) = default;
89  ~ExplicitInverse() = default;
90 
91  /// \cond
92  explicit ExplicitInverse(CkMigrateMessage* m) noexcept : Base(m) {}
93  using PUP::able::register_constructor;
95  /// \endcond
96 
97  /*!
98  * \brief Solve the equation \f$Ax=b\f$ by explicitly constructing the
99  * operator matrix \f$A\f$ and its inverse. The first solve is computationally
100  * expensive and successive solves are cheap.
101  *
102  * Building a matrix representation of the `linear_operator` requires
103  * iterating over the `SourceType` (which is also the type returned by the
104  * `linear_operator`) in a consistent way. This can be non-trivial for
105  * heterogeneous data structures because it requires they define a data
106  * ordering. Specifically, the `SourceType` must have a `size()` function as
107  * well as `begin()` and `end()` iterators that point into the data. If the
108  * iterators have a `reset()` function it is used to avoid repeatedly
109  * re-creating the `begin()` iterator. The `reset()` function must not
110  * invalidate the `end()` iterator.
111  */
112  template <typename LinearOperator, typename VarsType, typename SourceType>
114  const LinearOperator& linear_operator,
115  const SourceType& source) const noexcept;
116 
117  /// Flags the operator to require re-initialization. No memory is released.
118  /// Call this function to rebuild the solver when the operator changed.
119  void reset() noexcept override {
121  }
122 
123  /// Size of the operator. The stored matrix will have `size^2` entries.
124  size_t size() const noexcept { return size_; }
125 
126  /// The matrix representation of the solver. This matrix approximates the
127  /// inverse of the subdomain operator.
128  const DenseMatrix<double>& matrix_representation() const noexcept {
129  return inverse_;
130  }
131 
132  // NOLINTNEXTLINE(google-runtime-references)
133  void pup(PUP::er& p) noexcept override {
134  p | size_;
135  p | inverse_;
136  if (p.isUnpacking() and size_ != std::numeric_limits<size_t>::max()) {
137  source_workspace_.resize(size_);
138  solution_workspace_.resize(size_);
139  }
140  }
141 
142  std::unique_ptr<Base> get_clone() const noexcept override {
143  return std::make_unique<ExplicitInverse>(*this);
144  }
145 
146  private:
147  // Caches for successive solves of the same operator
148  mutable size_t size_ = std::numeric_limits<size_t>::max();
149  mutable DenseMatrix<double, blaze::columnMajor> inverse_{};
150 
151  // Buffers to avoid re-allocating memory for applying the operator
152  mutable DenseVector<double> source_workspace_{};
153  mutable DenseVector<double> solution_workspace_{};
154 };
155 
156 template <typename LinearSolverRegistrars>
157 template <typename LinearOperator, typename VarsType, typename SourceType>
159  const gsl::not_null<VarsType*> solution,
160  const LinearOperator& linear_operator,
161  const SourceType& source) const noexcept {
162  if (UNLIKELY(size_ == std::numeric_limits<size_t>::max())) {
163  const auto& used_for_size = source;
164  size_ = used_for_size.size();
165  source_workspace_.resize(size_);
166  solution_workspace_.resize(size_);
167  inverse_.resize(size_, size_);
168  // Construct explicit matrix representation by "sniffing out" the operator,
169  // i.e. feeding it unit vectors
170  auto unit_vector = make_with_value<VarsType>(used_for_size, 0.);
171  auto result_buffer = make_with_value<SourceType>(used_for_size, 0.);
172  auto& operator_matrix = inverse_;
173  size_t i = 0;
174  // Re-using the iterators for all operator invocations
175  auto result_iterator_begin = result_buffer.begin();
176  auto result_iterator_end = result_buffer.end();
177  for (double& unit_vector_data : unit_vector) {
178  // Add a 1 at the unit vector location i
179  unit_vector_data = 1.;
180  // Invoke the operator on the unit vector
181  linear_operator(make_not_null(&result_buffer), unit_vector);
182  // Set the unit vector back to zero
183  unit_vector_data = 0.;
184  // Reset the iterator by calling its `reset` member function or by
185  // re-creating it
186  if constexpr (detail::is_reset_callable_v<decltype(
187  result_iterator_begin)>) {
188  result_iterator_begin.reset();
189  } else {
190  result_iterator_begin = result_buffer.begin();
191  result_iterator_end = result_buffer.end();
192  }
193  // Store the result in column i of the matrix
194  std::copy(result_iterator_begin, result_iterator_end,
195  column(operator_matrix, i).begin());
196  ++i;
197  }
198  // Directly invert the matrix
199  try {
200  blaze::invert(inverse_);
201  } catch (const std::invalid_argument& e) {
202  ERROR("Could not invert subdomain matrix (size " << size_
203  << "): " << e.what());
204  }
205  }
206  // Copy source into contiguous workspace. In cases where the source and
207  // solution data are already stored contiguously we might avoid the copy and
208  // the associated workspace memory. However, compared to the cost of building
209  // and storing the matrix this is likely insignificant.
210  std::copy(source.begin(), source.end(), source_workspace_.begin());
211  // Apply inverse
212  solution_workspace_ = inverse_ * source_workspace_;
213  // Reconstruct solution data from contiguous workspace
214  std::copy(solution_workspace_.begin(), solution_workspace_.end(),
215  solution->begin());
216  return {0, 0};
217 }
218 
219 /// \cond
220 template <typename LinearSolverRegistrars>
221 // NOLINTNEXTLINE
222 PUP::able::PUP_ID ExplicitInverse<LinearSolverRegistrars>::my_PUP_ID = 0;
223 /// \endcond
224 
225 } // 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:124
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
CREATE_IS_CALLABLE_V
#define CREATE_IS_CALLABLE_V(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:68
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:119
algorithm
Registration::Registrar
A template for defining a registrar.
Definition: Registration.hpp:42
DenseMatrix< double >
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:37
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:128
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:73
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: ReadSpecPiecewisePolynomial.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:31