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