Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <memory>
8 : #include <optional>
9 : #include <pup.h>
10 : #include <type_traits>
11 : #include <utility>
12 :
13 : #include "DataStructures/DynamicMatrix.hpp"
14 : #include "DataStructures/DynamicVector.hpp"
15 : #include "IO/Logging/Verbosity.hpp"
16 : #include "NumericalAlgorithms/Convergence/Criteria.hpp"
17 : #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
18 : #include "NumericalAlgorithms/Convergence/Reason.hpp"
19 : #include "NumericalAlgorithms/LinearSolver/InnerProduct.hpp"
20 : #include "NumericalAlgorithms/LinearSolver/LinearSolver.hpp"
21 : #include "Options/Auto.hpp"
22 : #include "Options/String.hpp"
23 : #include "Utilities/EqualWithinRoundoff.hpp"
24 : #include "Utilities/Gsl.hpp"
25 : #include "Utilities/Registration.hpp"
26 : #include "Utilities/Serialization/CharmPupable.hpp"
27 : #include "Utilities/TMPL.hpp"
28 :
29 : namespace LinearSolver {
30 1 : namespace gmres::detail {
31 :
32 : // Perform an Arnoldi orthogonalization to find a new `operand` that is
33 : // orthogonal to all vectors in `basis_history`. Appends a new column to the
34 : // `orthogonalization_history` that holds the inner product of the intermediate
35 : // `operand` with each vector in the `basis_history` and itself.
36 : template <typename VarsType>
37 : void arnoldi_orthogonalize(const gsl::not_null<VarsType*> operand,
38 : const gsl::not_null<blaze::DynamicMatrix<double>*>
39 : orthogonalization_history,
40 : const std::vector<VarsType>& basis_history,
41 : const size_t iteration) {
42 : // Resize matrix and make sure the new entries that are not being filled below
43 : // are zero.
44 : orthogonalization_history->resize(iteration + 2, iteration + 1);
45 : for (size_t j = 0; j < iteration; ++j) {
46 : (*orthogonalization_history)(iteration + 1, j) = 0.;
47 : }
48 : // Arnoldi orthogonalization
49 : for (size_t j = 0; j < iteration + 1; ++j) {
50 : const double orthogonalization = inner_product(basis_history[j], *operand);
51 : (*orthogonalization_history)(j, iteration) = orthogonalization;
52 : *operand -= orthogonalization * basis_history[j];
53 : }
54 : (*orthogonalization_history)(iteration + 1, iteration) =
55 : sqrt(inner_product(*operand, *operand));
56 : // Avoid an FPE if the new operand norm is exactly zero. In that case the
57 : // problem is solved and the algorithm will terminate (see Proposition 9.3 in
58 : // \cite Saad2003). Since there will be no next iteration we don't need to
59 : // normalize the operand.
60 : if (UNLIKELY((*orthogonalization_history)(iteration + 1, iteration) == 0.)) {
61 : return;
62 : }
63 : *operand /= (*orthogonalization_history)(iteration + 1, iteration);
64 : }
65 :
66 : // Solve the linear least-squares problem `||beta - H * y||` for `y`, where `H`
67 : // is the Hessenberg matrix given by `orthogonalization_history` and `beta` is
68 : // the vector `(initial_residual, 0, 0, ...)` by updating the QR decomposition
69 : // of `H` from the previous iteration with a Givens rotation.
70 : void solve_minimal_residual(
71 : gsl::not_null<blaze::DynamicMatrix<double>*> orthogonalization_history,
72 : gsl::not_null<blaze::DynamicVector<double>*> residual_history,
73 : gsl::not_null<blaze::DynamicVector<double>*> givens_sine_history,
74 : gsl::not_null<blaze::DynamicVector<double>*> givens_cosine_history,
75 : size_t iteration);
76 :
77 : // Find the vector that minimizes the residual by inverting the upper
78 : // triangular matrix obtained above.
79 : blaze::DynamicVector<double> minimal_residual_vector(
80 : const blaze::DynamicMatrix<double>& orthogonalization_history,
81 : const blaze::DynamicVector<double>& residual_history);
82 :
83 : } // namespace gmres::detail
84 :
85 : namespace Serial {
86 :
87 : /// Disables the iteration callback at compile-time
88 1 : struct NoIterationCallback {};
89 :
90 : /// \cond
91 : template <typename VarsType, typename Preconditioner,
92 : typename LinearSolverRegistrars>
93 : struct Gmres;
94 : /// \endcond
95 :
96 : namespace Registrars {
97 :
98 : /// Registers the `LinearSolver::Serial::Gmres` linear solver.
99 : template <typename VarsType>
100 1 : struct Gmres {
101 : template <typename LinearSolverRegistrars>
102 0 : using f = Serial::Gmres<VarsType, LinearSolver<LinearSolverRegistrars>,
103 : LinearSolverRegistrars>;
104 : };
105 : } // namespace Registrars
106 :
107 : /*!
108 : * \brief A serial GMRES iterative solver for nonsymmetric linear systems of
109 : * equations.
110 : *
111 : * This is an iterative algorithm to solve general linear equations \f$Ax=b\f$
112 : * where \f$A\f$ is a linear operator. See \cite Saad2003, chapter 6.5 for a
113 : * description of the GMRES algorithm and Algorithm 9.6 for this implementation.
114 : * It is matrix-free, which means the operator \f$A\f$ needs not be provided
115 : * explicity as a matrix but only the operator action \f$A(x)\f$ must be
116 : * provided for an argument \f$x\f$.
117 : *
118 : * The GMRES algorithm does not require the operator \f$A\f$ to be symmetric or
119 : * positive-definite. Note that other algorithms such as conjugate gradients may
120 : * be more efficient for symmetric positive-definite operators.
121 : *
122 : * \par Convergence:
123 : * Given a set of \f$N_A\f$ equations (e.g. through an \f$N_A\times N_A\f$
124 : * matrix) the GMRES algorithm will converge to numerical precision in at most
125 : * \f$N_A\f$ iterations. However, depending on the properties of the linear
126 : * operator, an approximate solution can ideally be obtained in only a few
127 : * iterations. See \cite Saad2003, section 6.11.4 for details on the convergence
128 : * of the GMRES algorithm.
129 : *
130 : * \par Restarting:
131 : * This implementation of the GMRES algorithm supports restarting, as detailed
132 : * in \cite Saad2003, section 6.5.5. Since the GMRES algorithm iteratively
133 : * builds up an orthogonal basis of the solution space the cost of each
134 : * iteration increases linearly with the number of iterations. Therefore it is
135 : * sometimes helpful to restart the algorithm every \f$N_\mathrm{restart}\f$
136 : * iterations, discarding the set of basis vectors and starting again from the
137 : * current solution estimate. This strategy can improve the performance of the
138 : * solver, but note that the solver can stagnate for non-positive-definite
139 : * operators and is not guaranteed to converge within \f$N_A\f$ iterations
140 : * anymore. Set the `restart` argument of the constructor to
141 : * \f$N_\mathrm{restart}\f$ to activate restarting, or set it to 'None' to
142 : * deactivate restarting.
143 : *
144 : * \par Preconditioning:
145 : * This implementation of the GMRES algorithm also supports preconditioning.
146 : * You can provide a linear operator \f$P\f$ that approximates the inverse of
147 : * the operator \f$A\f$ to accelerate the convergence of the linear solve.
148 : * The algorithm is right-preconditioned, which allows the preconditioner to
149 : * change in every iteration ("flexible" variant). See \cite Saad2003, sections
150 : * 9.3.2 and 9.4.1 for details. This implementation follows Algorithm 9.6 in
151 : * \cite Saad2003.
152 : *
153 : * \par Improvements:
154 : * Further improvements can potentially be implemented for this algorithm, see
155 : * e.g. \cite Ayachour2003.
156 : *
157 : * \example
158 : * \snippet NumericalAlgorithms/LinearSolver/Test_Gmres.cpp gmres_example
159 : */
160 : template <typename VarsType, typename Preconditioner = NoPreconditioner,
161 : typename LinearSolverRegistrars =
162 : tmpl::list<Registrars::Gmres<VarsType>>>
163 1 : class Gmres final : public PreconditionedLinearSolver<Preconditioner,
164 : LinearSolverRegistrars> {
165 : private:
166 0 : using Base =
167 : PreconditionedLinearSolver<Preconditioner, LinearSolverRegistrars>;
168 :
169 0 : struct ConvergenceCriteria {
170 0 : using type = Convergence::Criteria;
171 0 : static constexpr Options::String help =
172 : "Determine convergence of the algorithm";
173 : };
174 0 : struct Restart {
175 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
176 0 : static constexpr Options::String help =
177 : "Iterations to run before restarting, or 'None' to disable restarting. "
178 : "Note that the solver is not guaranteed to converge anymore if you "
179 : "enable restarting.";
180 0 : static type suggested_value() { return {}; }
181 : };
182 0 : struct Verbosity {
183 0 : using type = ::Verbosity;
184 0 : static constexpr Options::String help = "Logging verbosity";
185 : };
186 :
187 : public:
188 0 : static constexpr Options::String help =
189 : "A serial GMRES iterative solver for nonsymmetric linear systems of\n"
190 : "equations Ax=b. It will converge to numerical precision in at most N_A\n"
191 : "iterations, where N_A is the number of equations represented by the\n"
192 : "linear operator A, but will ideally converge to a reasonable\n"
193 : "approximation of the solution x in only a few iterations.\n"
194 : "\n"
195 : "Preconditioning: Specify a preconditioner to run in every GMRES "
196 : "iteration to accelerate the solve, or 'None' to disable "
197 : "preconditioning. The choice of preconditioner can be crucial to obtain "
198 : "good convergence.\n"
199 : "\n"
200 : "Restarting: It is sometimes helpful to restart the algorithm every\n"
201 : "N_restart iterations to speed it up. Note that it can stagnate for\n"
202 : "non-positive-definite matrices and is not guaranteed to converge\n"
203 : "within N_A iterations anymore when restarting is activated.\n"
204 : "Activate restarting by setting the 'Restart' option to N_restart, or\n"
205 : "deactivate restarting by setting it to 'None'.";
206 0 : using options = tmpl::flatten<tmpl::list<
207 : ConvergenceCriteria, Verbosity, Restart,
208 : tmpl::conditional_t<std::is_same_v<Preconditioner, NoPreconditioner>,
209 : tmpl::list<>, typename Base::PreconditionerOption>>>;
210 :
211 0 : Gmres(Convergence::Criteria convergence_criteria, ::Verbosity verbosity,
212 : std::optional<size_t> restart = std::nullopt,
213 : std::optional<typename Base::PreconditionerType> local_preconditioner =
214 : std::nullopt,
215 : const Options::Context& context = {});
216 :
217 0 : Gmres() = default;
218 0 : Gmres(Gmres&&) = default;
219 0 : Gmres& operator=(Gmres&&) = default;
220 0 : ~Gmres() override = default;
221 :
222 0 : Gmres(const Gmres& rhs);
223 0 : Gmres& operator=(const Gmres& rhs);
224 :
225 0 : std::unique_ptr<LinearSolver<LinearSolverRegistrars>> get_clone()
226 : const override {
227 : return std::make_unique<Gmres>(*this);
228 : }
229 :
230 : /// \cond
231 : explicit Gmres(CkMigrateMessage* m);
232 : using PUP::able::register_constructor;
233 : WRAPPED_PUPable_decl_template(Gmres); // NOLINT
234 : /// \endcond
235 :
236 0 : void initialize() {
237 : orthogonalization_history_.reserve(restart_ + 1);
238 : residual_history_.reserve(restart_ + 1);
239 : givens_sine_history_.reserve(restart_);
240 : givens_cosine_history_.reserve(restart_);
241 : basis_history_.resize(restart_ + 1);
242 : preconditioned_basis_history_.resize(restart_);
243 : }
244 :
245 0 : const Convergence::Criteria& convergence_criteria() const {
246 : return convergence_criteria_;
247 : }
248 0 : ::Verbosity verbosity() const { return verbosity_; }
249 0 : size_t restart() const { return restart_; }
250 :
251 0 : void pup(PUP::er& p) override { // NOLINT
252 : Base::pup(p);
253 : p | convergence_criteria_;
254 : p | verbosity_;
255 : p | restart_;
256 : if (p.isUnpacking()) {
257 : initialize();
258 : }
259 : }
260 :
261 : template <typename LinearOperator, typename SourceType,
262 : typename... OperatorArgs,
263 : typename IterationCallback = NoIterationCallback>
264 0 : Convergence::HasConverged solve(
265 : gsl::not_null<VarsType*> initial_guess_in_solution_out,
266 : const LinearOperator& linear_operator, const SourceType& source,
267 : const std::tuple<OperatorArgs...>& operator_args = std::tuple{},
268 : const IterationCallback& iteration_callback =
269 : NoIterationCallback{}) const;
270 :
271 0 : void reset() override {
272 : // Nothing to reset. Only call into base class to reset preconditioner.
273 : Base::reset();
274 : }
275 :
276 : private:
277 0 : Convergence::Criteria convergence_criteria_{};
278 0 : ::Verbosity verbosity_{::Verbosity::Verbose};
279 0 : size_t restart_{};
280 :
281 : // Memory buffers to avoid re-allocating memory for successive solves:
282 : // The `orthogonalization_history_` is built iteratively from inner products
283 : // between existing and potential basis vectors and then Givens-rotated to
284 : // become upper-triangular.
285 : // NOLINTNEXTLINE(spectre-mutable)
286 0 : mutable blaze::DynamicMatrix<double> orthogonalization_history_{};
287 : // The `residual_history_` holds the remaining residual in its last entry, and
288 : // the other entries `g` "source" the minimum residual vector `y` in
289 : // `R * y = g` where `R` is the upper-triangular `orthogonalization_history_`.
290 : // NOLINTNEXTLINE(spectre-mutable)
291 0 : mutable blaze::DynamicVector<double> residual_history_{};
292 : // These represent the accumulated Givens rotations up to the current
293 : // iteration.
294 : // NOLINTNEXTLINE(spectre-mutable)
295 0 : mutable blaze::DynamicVector<double> givens_sine_history_{};
296 : // NOLINTNEXTLINE(spectre-mutable)
297 0 : mutable blaze::DynamicVector<double> givens_cosine_history_{};
298 : // These represent the orthogonal Krylov-subspace basis that is constructed
299 : // iteratively by Arnoldi-orthogonalizing a new vector in each iteration and
300 : // appending it to the `basis_history_`.
301 : // NOLINTNEXTLINE(spectre-mutable)
302 0 : mutable std::vector<VarsType> basis_history_{};
303 : // When a preconditioner is used it is applied to each new basis vector. The
304 : // preconditioned basis is used to construct the solution when the algorithm
305 : // has converged.
306 : // NOLINTNEXTLINE(spectre-mutable)
307 0 : mutable std::vector<VarsType> preconditioned_basis_history_{};
308 : };
309 :
310 : template <typename VarsType, typename Preconditioner,
311 : typename LinearSolverRegistrars>
312 : Gmres<VarsType, Preconditioner, LinearSolverRegistrars>::Gmres(
313 : Convergence::Criteria convergence_criteria, ::Verbosity verbosity,
314 : std::optional<size_t> restart,
315 : std::optional<typename Base::PreconditionerType> local_preconditioner,
316 : const Options::Context& context)
317 : // clang-tidy: trivially copyable
318 : : Base(std::move(local_preconditioner)),
319 : convergence_criteria_(std::move(convergence_criteria)), // NOLINT
320 : verbosity_(std::move(verbosity)), // NOLINT
321 : restart_(restart.value_or(convergence_criteria_.max_iterations)) {
322 : if (restart_ == 0) {
323 : PARSE_ERROR(context,
324 : "Can't restart every '0' iterations. Set to a nonzero "
325 : "number, or to 'None' if you meant to disable restarting.");
326 : }
327 : initialize();
328 : }
329 :
330 : // Define copy constructors. They don't have to copy the memory buffers but
331 : // only resize them. They take care of copying the preconditioner by calling
332 : // into the base class.
333 : template <typename VarsType, typename Preconditioner,
334 : typename LinearSolverRegistrars>
335 : Gmres<VarsType, Preconditioner, LinearSolverRegistrars>::Gmres(const Gmres& rhs)
336 : : Base(rhs),
337 : convergence_criteria_(rhs.convergence_criteria_),
338 : verbosity_(rhs.verbosity_),
339 : restart_(rhs.restart_) {
340 : initialize();
341 : }
342 : template <typename VarsType, typename Preconditioner,
343 : typename LinearSolverRegistrars>
344 : Gmres<VarsType, Preconditioner, LinearSolverRegistrars>&
345 : Gmres<VarsType, Preconditioner, LinearSolverRegistrars>::operator=(
346 : const Gmres& rhs) {
347 : Base::operator=(rhs);
348 : convergence_criteria_ = rhs.convergence_criteria_;
349 : verbosity_ = rhs.verbosity_;
350 : restart_ = rhs.restart_;
351 : initialize();
352 : return *this;
353 : }
354 :
355 : /// \cond
356 : template <typename VarsType, typename Preconditioner,
357 : typename LinearSolverRegistrars>
358 : Gmres<VarsType, Preconditioner, LinearSolverRegistrars>::Gmres(
359 : CkMigrateMessage* m)
360 : : Base(m) {}
361 : /// \endcond
362 :
363 : template <typename VarsType, typename Preconditioner,
364 : typename LinearSolverRegistrars>
365 : template <typename LinearOperator, typename SourceType,
366 : typename... OperatorArgs, typename IterationCallback>
367 : Convergence::HasConverged
368 0 : Gmres<VarsType, Preconditioner, LinearSolverRegistrars>::solve(
369 : const gsl::not_null<VarsType*> initial_guess_in_solution_out,
370 : const LinearOperator& linear_operator, const SourceType& source,
371 : const std::tuple<OperatorArgs...>& operator_args,
372 : const IterationCallback& iteration_callback) const {
373 : constexpr bool use_preconditioner =
374 : not std::is_same_v<Preconditioner, NoPreconditioner>;
375 : constexpr bool use_iteration_callback =
376 : not std::is_same_v<IterationCallback, NoIterationCallback>;
377 :
378 : // Could pre-allocate memory for the basis-history vectors here. Not doing
379 : // that for now because we don't know how many iterations we'll need.
380 : // Estimating the number of iterations and pre-allocating memory is a possible
381 : // performance optimization.
382 :
383 : auto& solution = *initial_guess_in_solution_out;
384 : Convergence::HasConverged has_converged{};
385 : size_t iteration = 0;
386 :
387 : while (not has_converged) {
388 : const auto& initial_guess = *initial_guess_in_solution_out;
389 : auto& initial_operand = basis_history_[0];
390 : // Apply the linear operator to the initial guess. This can be skipped if
391 : // the initial guess is zero, because then the linear operator applied to it
392 : // is also zero.
393 : if (equal_within_roundoff(initial_guess, 0.)) {
394 : initial_operand = source;
395 : } else {
396 : std::apply(
397 : linear_operator,
398 : std::tuple_cat(std::forward_as_tuple(make_not_null(&initial_operand),
399 : initial_guess),
400 : operator_args));
401 : initial_operand *= -1.;
402 : initial_operand += source;
403 : }
404 : const double initial_residual_magnitude =
405 : sqrt(inner_product(initial_operand, initial_operand));
406 : has_converged = Convergence::HasConverged{convergence_criteria_, iteration,
407 : initial_residual_magnitude,
408 : initial_residual_magnitude};
409 : if constexpr (use_iteration_callback) {
410 : iteration_callback(has_converged);
411 : }
412 : if (UNLIKELY(has_converged)) {
413 : break;
414 : }
415 : initial_operand /= initial_residual_magnitude;
416 : residual_history_.resize(1);
417 : residual_history_[0] = initial_residual_magnitude;
418 : for (size_t k = 0; k < restart_; ++k) {
419 : auto& operand = basis_history_[k + 1];
420 : if constexpr (use_preconditioner) {
421 : if (this->has_preconditioner()) {
422 : // Begin the preconditioner at an initial guess of 0. Not all
423 : // preconditioners take the initial guess into account.
424 : preconditioned_basis_history_[k] =
425 : make_with_value<VarsType>(initial_operand, 0.);
426 : this->preconditioner().solve(
427 : make_not_null(&preconditioned_basis_history_[k]), linear_operator,
428 : basis_history_[k], operator_args);
429 : }
430 : }
431 : std::apply(linear_operator,
432 : std::tuple_cat(std::forward_as_tuple(
433 : make_not_null(&operand),
434 : this->has_preconditioner()
435 : ? preconditioned_basis_history_[k]
436 : : basis_history_[k]),
437 : operator_args));
438 : // Find a new orthogonal basis vector of the Krylov subspace
439 : gmres::detail::arnoldi_orthogonalize(
440 : make_not_null(&operand), make_not_null(&orthogonalization_history_),
441 : basis_history_, k);
442 : // Least-squares solve for the minimal residual
443 : gmres::detail::solve_minimal_residual(
444 : make_not_null(&orthogonalization_history_),
445 : make_not_null(&residual_history_),
446 : make_not_null(&givens_sine_history_),
447 : make_not_null(&givens_cosine_history_), k);
448 : ++iteration;
449 : has_converged = Convergence::HasConverged{
450 : convergence_criteria_, iteration, abs(residual_history_[k + 1]),
451 : initial_residual_magnitude};
452 : if constexpr (use_iteration_callback) {
453 : iteration_callback(has_converged);
454 : }
455 : if (UNLIKELY(has_converged)) {
456 : break;
457 : }
458 : }
459 : // Find the vector w.r.t. the constructed orthogonal basis of the Krylov
460 : // subspace that minimizes the residual
461 : const auto minres = gmres::detail::minimal_residual_vector(
462 : orthogonalization_history_, residual_history_);
463 : // Construct the solution from the orthogonal basis and the minimal residual
464 : // vector
465 : for (size_t i = 0; i < minres.size(); ++i) {
466 : solution += minres[i] * gsl::at(this->has_preconditioner()
467 : ? preconditioned_basis_history_
468 : : basis_history_,
469 : i);
470 : }
471 : }
472 : return has_converged;
473 : }
474 :
475 : /// \cond
476 : template <typename VarsType, typename Preconditioner,
477 : typename LinearSolverRegistrars>
478 : // NOLINTNEXTLINE
479 : PUP::able::PUP_ID
480 : Gmres<VarsType, Preconditioner, LinearSolverRegistrars>::my_PUP_ID = 0;
481 : /// \endcond
482 :
483 : } // namespace Serial
484 : } // namespace LinearSolver
|
