15 #include "Informer/Verbosity.hpp"
16 #include "NumericalAlgorithms/Convergence/Criteria.hpp"
17 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
18 #include "NumericalAlgorithms/Convergence/Reason.hpp"
20 #include "NumericalAlgorithms/LinearSolver/LinearSolver.hpp"
21 #include "Options/Auto.hpp"
25 #include "Utilities/Registration.hpp"
29 namespace gmres::detail {
35 template <
typename VarsType>
36 void arnoldi_orthogonalize(
40 const size_t iteration) noexcept {
43 orthogonalization_history->resize(iteration + 2, iteration + 1);
44 for (
size_t j = 0; j < iteration; ++j) {
45 (*orthogonalization_history)(iteration + 1, j) = 0.;
48 for (
size_t j = 0; j < iteration + 1; ++j) {
49 const double orthogonalization =
inner_product(basis_history[j], *operand);
50 (*orthogonalization_history)(j, iteration) = orthogonalization;
51 *operand -= orthogonalization * basis_history[j];
53 (*orthogonalization_history)(iteration + 1, iteration) =
59 if (
UNLIKELY((*orthogonalization_history)(iteration + 1, iteration) == 0.)) {
62 *operand /= (*orthogonalization_history)(iteration + 1, iteration);
69 void solve_minimal_residual(
74 size_t iteration) noexcept;
90 template <
typename VarsType,
typename Preconditioner,
91 typename LinearSolverRegistrars>
95 namespace Registrars {
98 template <
typename VarsType>
100 template <
typename LinearSolverRegistrars>
102 LinearSolverRegistrars>;
160 typename LinearSolverRegistrars =
161 tmpl::list<Registrars::Gmres<VarsType>>>
163 LinearSolverRegistrars> {
168 struct ConvergenceCriteria {
171 "Determine convergence of the algorithm";
176 "Iterations to run before restarting, or 'None' to disable restarting. "
177 "Note that the solver is not guaranteed to converge anymore if you "
178 "enable restarting.";
179 static type suggested_value() noexcept {
return {}; }
188 "A serial GMRES iterative solver for nonsymmetric linear systems of\n"
189 "equations Ax=b. It will converge to numerical precision in at most N_A\n"
190 "iterations, where N_A is the number of equations represented by the\n"
191 "linear operator A, but will ideally converge to a reasonable\n"
192 "approximation of the solution x in only a few iterations.\n"
194 "Preconditioning: Specify a preconditioner to run in every GMRES "
195 "iteration to accelerate the solve, or 'None' to disable "
196 "preconditioning. The choice of preconditioner can be crucial to obtain "
197 "good convergence.\n"
199 "Restarting: It is sometimes helpful to restart the algorithm every\n"
200 "N_restart iterations to speed it up. Note that it can stagnate for\n"
201 "non-positive-definite matrices and is not guaranteed to converge\n"
202 "within N_A iterations anymore when restarting is activated.\n"
203 "Activate restarting by setting the 'Restart' option to N_restart, or\n"
204 "deactivate restarting by setting it to 'None'.";
205 using options = tmpl::flatten<tmpl::list<
207 tmpl::conditional_t<std::is_same_v<Preconditioner, NoPreconditioner>,
219 ~
Gmres()
override =
default;
226 return std::make_unique<Gmres>(*
this);
230 explicit Gmres(CkMigrateMessage* m) noexcept;
231 using PUP::able::register_constructor;
235 void initialize() noexcept {
236 orthogonalization_history_.reserve(restart_ + 1);
237 residual_history_.reserve(restart_ + 1);
238 givens_sine_history_.reserve(restart_);
239 givens_cosine_history_.reserve(restart_);
240 basis_history_.resize(restart_ + 1);
241 preconditioned_basis_history_.resize(restart_);
245 return convergence_criteria_;
247 ::Verbosity verbosity()
const noexcept {
return verbosity_; }
248 size_t restart()
const noexcept {
return restart_; }
250 void pup(PUP::er& p) noexcept
override {
252 p | convergence_criteria_;
255 if (p.isUnpacking()) {
260 template <
typename LinearOperator,
typename SourceType,
264 const LinearOperator& linear_operator,
const SourceType& source,
265 const IterationCallback& iteration_callback =
268 void reset() noexcept
override {
301 template <
typename VarsType,
typename Preconditioner,
302 typename LinearSolverRegistrars>
309 : Base(std::move(local_preconditioner)),
310 convergence_criteria_(std::move(convergence_criteria)),
311 verbosity_(std::move(verbosity)),
312 restart_(restart.value_or(convergence_criteria_.max_iterations)) {
315 "Can't restart every '0' iterations. Set to a nonzero "
316 "number, or to 'None' if you meant to disable restarting.");
324 template <
typename VarsType,
typename Preconditioner,
325 typename LinearSolverRegistrars>
326 Gmres<VarsType, Preconditioner, LinearSolverRegistrars>::Gmres(
327 const Gmres& rhs) noexcept
329 convergence_criteria_(rhs.convergence_criteria_),
330 verbosity_(rhs.verbosity_),
331 restart_(rhs.restart_) {
334 template <
typename VarsType,
typename Preconditioner,
335 typename LinearSolverRegistrars>
336 Gmres<VarsType, Preconditioner, LinearSolverRegistrars>&
337 Gmres<VarsType, Preconditioner, LinearSolverRegistrars>::operator=(
338 const Gmres& rhs) noexcept {
339 Base::operator=(rhs);
340 convergence_criteria_ = rhs.convergence_criteria_;
341 verbosity_ = rhs.verbosity_;
342 restart_ = rhs.restart_;
348 template <
typename VarsType,
typename Preconditioner,
349 typename LinearSolverRegistrars>
350 Gmres<VarsType, Preconditioner, LinearSolverRegistrars>::Gmres(
351 CkMigrateMessage* m) noexcept
355 template <
typename VarsType,
typename Preconditioner,
356 typename LinearSolverRegistrars>
357 template <
typename LinearOperator,
typename SourceType,
358 typename IterationCallback>
360 Gmres<VarsType, Preconditioner, LinearSolverRegistrars>::solve(
362 const LinearOperator& linear_operator,
const SourceType& source,
363 const IterationCallback& iteration_callback)
const noexcept {
364 constexpr
bool use_preconditioner =
365 not std::is_same_v<Preconditioner, NoPreconditioner>;
366 constexpr
bool use_iteration_callback =
367 not std::is_same_v<IterationCallback, NoIterationCallback>;
374 auto& solution = *initial_guess_in_solution_out;
376 size_t iteration = 0;
378 while (not has_converged) {
379 const auto& initial_guess = *initial_guess_in_solution_out;
380 auto& initial_operand = basis_history_[0];
381 linear_operator(
make_not_null(&initial_operand), initial_guess);
382 initial_operand *= -1.;
383 initial_operand += source;
384 const double initial_residual_magnitude =
387 initial_residual_magnitude,
388 initial_residual_magnitude};
389 if constexpr (use_iteration_callback) {
390 iteration_callback(has_converged);
395 initial_operand /= initial_residual_magnitude;
396 residual_history_.resize(1);
397 residual_history_[0] = initial_residual_magnitude;
398 for (
size_t k = 0; k < restart_; ++k) {
399 auto& operand = basis_history_[k + 1];
400 if constexpr (use_preconditioner) {
401 if (this->has_preconditioner()) {
404 preconditioned_basis_history_[k] =
405 make_with_value<VarsType>(initial_operand, 0.);
406 this->preconditioner().solve(
407 make_not_null(&preconditioned_basis_history_[k]), linear_operator,
412 this->has_preconditioner()
413 ? preconditioned_basis_history_[k]
414 : basis_history_[k]);
416 gmres::detail::arnoldi_orthogonalize(
420 gmres::detail::solve_minimal_residual(
427 convergence_criteria_, iteration, abs(residual_history_[k + 1]),
428 initial_residual_magnitude};
429 if constexpr (use_iteration_callback) {
430 iteration_callback(has_converged);
438 const auto minres = gmres::detail::minimal_residual_vector(
439 orthogonalization_history_, residual_history_);
442 for (
size_t i = 0; i < minres.size(); ++i) {
443 solution += minres[i] *
gsl::at(this->has_preconditioner()
444 ? preconditioned_basis_history_
449 return has_converged;
453 template <
typename VarsType,
typename Preconditioner,
454 typename LinearSolverRegistrars>
457 Gmres<VarsType, Preconditioner, LinearSolverRegistrars>::my_PUP_ID = 0;