Gmres.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <boost/none.hpp>
7 #include <boost/optional.hpp>
8 #include <boost/optional/optional_io.hpp>
9 #include <cstddef>
10 
13 #include "Informer/Verbosity.hpp"
14 #include "NumericalAlgorithms/Convergence/Criteria.hpp"
15 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
16 #include "NumericalAlgorithms/Convergence/Reason.hpp"
18 #include "Options/Options.hpp"
19 #include "Utilities/Gsl.hpp"
20 
21 namespace LinearSolver {
22 namespace gmres::detail {
23 
24 // Perform an Arnoldi orthogonalization to find a new `operand` that is
25 // orthogonal to all vectors in `basis_history`. Appends a new column to the
26 // `orthogonalization_history` that holds the inner product of the intermediate
27 // `operand` with each vector in the `basis_history` and itself.
28 template <typename VarsType>
29 void arnoldi_orthogonalize(
30  const gsl::not_null<VarsType*> operand,
31  const gsl::not_null<DenseMatrix<double>*> orthogonalization_history,
32  const std::vector<VarsType>& basis_history,
33  const size_t iteration) noexcept {
34  // Resize matrix and make sure the new entries that are not being filled below
35  // are zero.
36  orthogonalization_history->resize(iteration + 2, iteration + 1);
37  for (size_t j = 0; j < iteration; ++j) {
38  (*orthogonalization_history)(iteration + 1, j) = 0.;
39  }
40  // Arnoldi orthogonalization
41  for (size_t j = 0; j < iteration + 1; ++j) {
42  const double orthogonalization = inner_product(basis_history[j], *operand);
43  (*orthogonalization_history)(j, iteration) = orthogonalization;
44  *operand -= orthogonalization * basis_history[j];
45  }
46  (*orthogonalization_history)(iteration + 1, iteration) =
47  sqrt(inner_product(*operand, *operand));
48  // Avoid an FPE if the new operand norm is exactly zero. In that case the
49  // problem is solved and the algorithm will terminate (see Proposition 9.3 in
50  // \cite Saad2003). Since there will be no next iteration we don't need to
51  // normalize the operand.
52  if (UNLIKELY((*orthogonalization_history)(iteration + 1, iteration) == 0.)) {
53  return;
54  }
55  *operand /= (*orthogonalization_history)(iteration + 1, iteration);
56 }
57 
58 // Solve the linear least-squares problem `||beta - H * y||` for `y`, where `H`
59 // is the Hessenberg matrix given by `orthogonalization_history` and `beta` is
60 // the vector `(initial_residual, 0, 0, ...)` by updating the QR decomposition
61 // of `H` from the previous iteration with a Givens rotation.
62 void solve_minimal_residual(
63  gsl::not_null<DenseMatrix<double>*> orthogonalization_history,
64  gsl::not_null<DenseVector<double>*> residual_history,
65  gsl::not_null<DenseVector<double>*> givens_sine_history,
66  gsl::not_null<DenseVector<double>*> givens_cosine_history,
67  size_t iteration) noexcept;
68 
69 // Find the vector that minimizes the residual by inverting the upper
70 // triangular matrix obtained above.
71 DenseVector<double> minimal_residual_vector(
72  const DenseMatrix<double>& orthogonalization_history,
73  const DenseVector<double>& residual_history) noexcept;
74 
75 } // namespace gmres::detail
76 
77 namespace Serial {
78 
79 template <typename VarsType>
81  VarsType operator()(const VarsType& arg) const noexcept { return arg; }
82 };
83 
85  void operator()(const Convergence::HasConverged& /*has_converged*/) const
86  noexcept {}
87 };
88 
89 /*!
90  * \brief A serial GMRES iterative solver for nonsymmetric linear systems of
91  * equations.
92  *
93  * This is an iterative algorithm to solve general linear equations \f$Ax=b\f$
94  * where \f$A\f$ is a linear operator. See \cite Saad2003, chapter 6.5 for a
95  * description of the GMRES algorithm and Algorithm 9.6 for this implementation.
96  * It is matrix-free, which means the operator \f$A\f$ needs not be provided
97  * explicity as a matrix but only the operator action \f$A(x)\f$ must be
98  * provided for an argument \f$x\f$.
99  *
100  * The GMRES algorithm does not require the operator \f$A\f$ to be symmetric or
101  * positive-definite. Note that other algorithms such as conjugate gradients may
102  * be more efficient for symmetric positive-definite operators.
103  *
104  * \par Convergence:
105  * Given a set of \f$N_A\f$ equations (e.g. through an \f$N_A\times N_A\f$
106  * matrix) the GMRES algorithm will converge to numerical precision in at most
107  * \f$N_A\f$ iterations. However, depending on the properties of the linear
108  * operator, an approximate solution can ideally be obtained in only a few
109  * iterations. See \cite Saad2003, section 6.11.4 for details on the convergence
110  * of the GMRES algorithm.
111  *
112  * \par Restarting:
113  * This implementation of the GMRES algorithm supports restarting, as detailed
114  * in \cite Saad2003, section 6.5.5. Since the GMRES algorithm iteratively
115  * builds up an orthogonal basis of the solution space the cost of each
116  * iteration increases linearly with the number of iterations. Therefore it is
117  * sometimes helpful to restart the algorithm every \f$N_\mathrm{restart}\f$
118  * iterations, discarding the set of basis vectors and starting again from the
119  * current solution estimate. This strategy can improve the performance of the
120  * solver, but note that the solver can stagnate for non-positive-definite
121  * operators and is not guaranteed to converge within \f$N_A\f$ iterations
122  * anymore. Set the `restart` argument of the constructor to
123  * \f$N_\mathrm{restart}\f$ to activate restarting, or set it to zero to
124  * deactivate restarting (default behaviour).
125  *
126  * \par Preconditioning:
127  * This implementation of the GMRES algorithm also supports preconditioning.
128  * You can provide a linear operator \f$P\f$ that approximates the inverse of
129  * the operator \f$A\f$ to accelerate the convergence of the linear solve.
130  * The algorithm is right-preconditioned, which allows the preconditioner to
131  * change in every iteration ("flexible" variant). See \cite Saad2003, sections
132  * 9.3.2 and 9.4.1 for details. This implementation follows Algorithm 9.6 in
133  * \cite Saad2003.
134  *
135  * \par Improvements:
136  * Further improvements can potentially be implemented for this algorithm, see
137  * e.g. \cite Ayachour2003.
138  *
139  * \example
140  * \snippet NumericalAlgorithms/LinearSolver/Test_Gmres.cpp gmres_example
141  */
142 template <typename VarsType>
143 struct Gmres {
144  private:
145  struct ConvergenceCriteria {
146  using type = Convergence::Criteria;
147  static constexpr Options::String help =
148  "Determine convergence of the algorithm";
149  };
150  struct Restart {
151  using type = size_t;
152  static constexpr Options::String help =
153  "Iterations to run before restarting";
154  static size_t default_value() noexcept { return 0; }
155  };
156  struct Verbosity {
157  using type = ::Verbosity;
158  static constexpr Options::String help = "Logging verbosity";
159  };
160 
161  public:
162  static constexpr Options::String help =
163  "A serial GMRES iterative solver for nonsymmetric linear systems of\n"
164  "equations Ax=b. It will converge to numerical precision in at most N_A\n"
165  "iterations, where N_A is the number of equations represented by the\n"
166  "linear operator A, but will ideally converge to a reasonable\n"
167  "approximation of the solution x in only a few iterations.\n"
168  "\n"
169  "Restarting: It is sometimes helpful to restart the algorithm every\n"
170  "N_restart iterations to speed it up. Note that it can stagnate for\n"
171  "non-positive-definite matrices and is not guaranteed to converge\n"
172  "within N_A iterations anymore when restarting is activated.\n"
173  "Activate restarting by setting the 'Restart' option to N_restart, or\n"
174  "deactivate restarting by setting it to zero (default).";
175  using options = tmpl::list<ConvergenceCriteria, Verbosity, Restart>;
176 
177  Gmres(Convergence::Criteria convergence_criteria, ::Verbosity verbosity,
178  size_t restart = 0) noexcept
179  // clang-tidy: trivially copyable
180  : convergence_criteria_(std::move(convergence_criteria)), // NOLINT
181  verbosity_(std::move(verbosity)), // NOLINT
182  restart_(restart > 0 ? restart : convergence_criteria_.max_iterations) {
183  initialize();
184  }
185 
186  Gmres() = default;
187  Gmres(const Gmres& /*rhs*/) = default;
188  Gmres& operator=(const Gmres& /*rhs*/) = default;
189  Gmres(Gmres&& /*rhs*/) = default;
190  Gmres& operator=(Gmres&& /*rhs*/) = default;
191  ~Gmres() = default;
192 
193  void initialize() noexcept {
194  orthogonalization_history_.reserve(restart_ + 1);
195  residual_history_.reserve(restart_ + 1);
196  givens_sine_history_.reserve(restart_);
197  givens_cosine_history_.reserve(restart_);
198  basis_history_.resize(restart_ + 1);
199  preconditioned_basis_history_.resize(restart_);
200  }
201 
202  const Convergence::Criteria& convergence_criteria() const noexcept {
203  return convergence_criteria_;
204  }
205 
206  void pup(PUP::er& p) noexcept { // NOLINT
207  p | convergence_criteria_;
208  p | verbosity_;
209  p | restart_;
210  if (p.isUnpacking()) {
211  initialize();
212  }
213  }
214 
215  /*!
216  * \brief Iteratively solve the problem \f$Ax=b\f$ for \f$x\f$ where \f$A\f$
217  * is the `linear_operator` and \f$b\f$ is the `source`, starting \f$x\f$ at
218  * `initial_guess`.
219  *
220  * Optionally provide a `preconditioner` (see class documentation).
221  *
222  * \return An instance of `Convergence::HasConverged` that provides
223  * information on the convergence status of the completed solve, and the
224  * approximate solution \f$x\f$.
225  */
226  template <typename LinearOperator, typename SourceType,
227  typename Preconditioner = IdentityPreconditioner<VarsType>,
228  typename IterationCallback = NoIterationCallback>
230  LinearOperator&& linear_operator, const SourceType& source,
231  const VarsType& initial_guess,
232  Preconditioner&& preconditioner = IdentityPreconditioner<VarsType>{},
233  IterationCallback&& iteration_callback = NoIterationCallback{}) const
234  noexcept;
235 
236  private:
237  Convergence::Criteria convergence_criteria_{};
238  ::Verbosity verbosity_{::Verbosity::Verbose};
239  size_t restart_{};
240 
241  // Memory buffers to avoid re-allocating memory for successive solves:
242  // The `orthogonalization_history_` is built iteratively from inner products
243  // between existing and potential basis vectors and then Givens-rotated to
244  // become upper-triangular.
245  mutable DenseMatrix<double> orthogonalization_history_{};
246  // The `residual_history_` holds the remaining residual in its last entry, and
247  // the other entries `g` "source" the minimum residual vector `y` in
248  // `R * y = g` where `R` is the upper-triangular `orthogonalization_history_`.
249  mutable DenseVector<double> residual_history_{};
250  // These represent the accumulated Givens rotations up to the current
251  // iteration.
252  mutable DenseVector<double> givens_sine_history_{};
253  mutable DenseVector<double> givens_cosine_history_{};
254  // These represent the orthogonal Krylov-subspace basis that is constructed
255  // iteratively by Arnoldi-orthogonalizing a new vector in each iteration and
256  // appending it to the `basis_history_`.
257  mutable std::vector<VarsType> basis_history_{};
258  // When a preconditioner is used it is applied to each new basis vector. The
259  // preconditioned basis is used to construct the solution when the algorithm
260  // has converged.
261  mutable std::vector<VarsType> preconditioned_basis_history_{};
262 };
263 
264 template <typename VarsType>
265 template <typename LinearOperator, typename SourceType, typename Preconditioner,
266  typename IterationCallback>
268  LinearOperator&& linear_operator, const SourceType& source,
269  const VarsType& initial_guess, Preconditioner&& preconditioner,
270  IterationCallback&& iteration_callback) const noexcept {
271  constexpr bool use_preconditioner =
272  not std::is_same_v<Preconditioner, IdentityPreconditioner<VarsType>>;
273  constexpr bool use_iteration_callback =
274  not std::is_same_v<IterationCallback, NoIterationCallback>;
275 
276  auto result = initial_guess;
277  Convergence::HasConverged has_converged{};
278  size_t iteration = 0;
279 
280  while (not has_converged) {
281  auto& initial_operand = basis_history_[0] =
282  source - linear_operator(result);
283  const double initial_residual_magnitude =
284  sqrt(inner_product(initial_operand, initial_operand));
285  has_converged = Convergence::HasConverged{convergence_criteria_, iteration,
286  initial_residual_magnitude,
287  initial_residual_magnitude};
288  if (use_iteration_callback) {
289  iteration_callback(has_converged);
290  }
291  if (UNLIKELY(has_converged)) {
292  break;
293  }
294  initial_operand /= initial_residual_magnitude;
295  residual_history_.resize(1);
296  residual_history_[0] = initial_residual_magnitude;
297  for (size_t k = 0; k < restart_; ++k) {
298  auto& operand = basis_history_[k + 1];
299  if (use_preconditioner) {
300  preconditioned_basis_history_[k] = preconditioner(basis_history_[k]);
301  }
302  operand =
303  linear_operator(use_preconditioner ? preconditioned_basis_history_[k]
304  : basis_history_[k]);
305  // Find a new orthogonal basis vector of the Krylov subspace
306  gmres::detail::arnoldi_orthogonalize(
307  make_not_null(&operand), make_not_null(&orthogonalization_history_),
308  basis_history_, k);
309  // Least-squares solve for the minimal residual
310  gmres::detail::solve_minimal_residual(
311  make_not_null(&orthogonalization_history_),
312  make_not_null(&residual_history_),
313  make_not_null(&givens_sine_history_),
314  make_not_null(&givens_cosine_history_), k);
315  ++iteration;
316  has_converged = Convergence::HasConverged{
317  convergence_criteria_, iteration, abs(residual_history_[k + 1]),
318  initial_residual_magnitude};
319  if (use_iteration_callback) {
320  iteration_callback(has_converged);
321  }
322  if (UNLIKELY(has_converged)) {
323  break;
324  }
325  }
326  // Find the vector w.r.t. the constructed orthogonal basis of the Krylov
327  // subspace that minimizes the residual
328  const auto minres = gmres::detail::minimal_residual_vector(
329  orthogonalization_history_, residual_history_);
330  // Construct the solution from the orthogonal basis and the minimal residual
331  // vector
332  for (size_t i = 0; i < minres.size(); ++i) {
333  result +=
334  minres[i] * gsl::at(use_preconditioner ? preconditioned_basis_history_
335  : basis_history_,
336  i);
337  }
338  }
339  // NOLINTNEXTLINE(hicpp-move-const-arg,performance-move-const-arg)
340  return {std::move(has_converged), std::move(result)};
341 }
342 
343 } // namespace Serial
344 } // namespace LinearSolver
Verbosity
Verbosity
Indicates how much informative output a class should output.
Definition: Verbosity.hpp:18
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:125
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
InnerProduct.hpp
std::pair
LinearSolver::Serial::Gmres::operator()
std::pair< Convergence::HasConverged, VarsType > operator()(LinearOperator &&linear_operator, const SourceType &source, const VarsType &initial_guess, Preconditioner &&preconditioner=IdentityPreconditioner< VarsType >{}, IterationCallback &&iteration_callback=NoIterationCallback{}) const noexcept
Iteratively solve the problem for where is the linear_operator and is the source,...
Definition: Gmres.hpp:267
Options.hpp
std::vector< VarsType >
Convergence::Criteria
Criteria that determine an iterative algorithm has converged.
Definition: Criteria.hpp:35
DenseVector.hpp
DenseMatrix.hpp
DenseVector< double >
Convergence::HasConverged
Signals convergence of the algorithm.
Definition: HasConverged.hpp:59
std::sqrt
T sqrt(T... args)
LinearSolver::Serial::Gmres
A serial GMRES iterative solver for nonsymmetric linear systems of equations.
Definition: Gmres.hpp:143
DenseMatrix< double >
LinearSolver::inner_product
double inner_product(const Lhs &lhs, const Rhs &rhs) noexcept
The local part of the Euclidean inner product on the vector space w.r.t. which the addition and scala...
Definition: InnerProduct.hpp:71
cstddef
LinearSolver::Serial::IdentityPreconditioner
Definition: Gmres.hpp:80
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::NoIterationCallback
Definition: Gmres.hpp:84
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
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183