From 83d3d9a4f4260c9f1d681c44eef1d7cb461f6f68 Mon Sep 17 00:00:00 2001 From: Dario Coscia Date: Tue, 3 Sep 2024 16:46:20 +0200 Subject: [PATCH] adding additional data for solver cg Co-authored-by: Martin Kronbichler richer docstring --- include/deal.II/lac/solver_cg.h | 184 +++++++++++++++++++++++++++----- tests/lac/solver.cc | 8 ++ tests/lac/solver.output | 47 ++++++-- 3 files changed, 204 insertions(+), 35 deletions(-) diff --git a/include/deal.II/lac/solver_cg.h b/include/deal.II/lac/solver_cg.h index dca6a25b68..9c2d48d8af 100644 --- a/include/deal.II/lac/solver_cg.h +++ b/include/deal.II/lac/solver_cg.h @@ -59,7 +59,8 @@ namespace LinearAlgebra * (for the requirements on matrices and vectors in order to work with this * class, see the documentation of the Solver base class). The type of the * solution vector must be passed as template argument, and defaults to - * dealii::Vector. + * dealii::Vector. The AdditionalData structure allows to control the + * type of residual for the stopping condition. * * @note The CG method requires a symmetric preconditioner (i.e., for example, * SOR is not a possible choice). There is a variant of the solver, @@ -172,6 +173,23 @@ namespace LinearAlgebra * run the updates on the vectors according to a variant presented in * Algorithm 2.2 of @cite Chronopoulos1989 (but for a preconditioner), whereas * the operation after the loop performs a total of 7 reductions in parallel. + * + *

Preconditioned residual

+ * + * @p AdditionalData allows you to choose between using the explicit + * or implicit residual as a stopping condition for the iterative + * solver. This behavior can be overridden by using the flag + * AdditionalData::use_default_residual. A true value refers to the + * implicit residual, while false reverts + * it. The former uses the result of the matrix-vector + * product already computed in other algorithm steps to derive the residual by + * a mere vector update, whereas the latter explicitly calculates the system + * residual with an additional matrix-vector product. + * More information on explicit and implicit residual stopping + * criteria can be found + * link + * here. */ template > DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) @@ -183,13 +201,34 @@ public: */ using size_type = types::global_dof_index; + /** * Standardized data struct to pipe additional data to the solver. - * Here, it does not store anything but just exists for consistency - * with the other solver classes. */ struct AdditionalData - {}; + { + /** + * Constructor. By default, set the residual of the stopping criterion + * to the implicit residual. A true value of + * AdditionalData::use_default_residual refers to the + * implicit residual, while false reverts + * it. The former uses the result of the matrix-vector + * product already computed in other algorithm steps to derive the residual + * by a mere vector update, whereas the latter explicitly calculates the + * system residual with an additional matrix-vector product. More + * information on explicit and implicit residual stopping criteria can be + * found link + * here. + */ + explicit AdditionalData(const bool use_default_residual = true); + + /** + * Flag for the default residual that is used to measure convergence. + */ + bool use_default_residual; + }; + /** * Constructor. @@ -363,13 +402,34 @@ public: */ using size_type = types::global_dof_index; + /** * Standardized data struct to pipe additional data to the solver. - * Here, it does not store anything but just exists for consistency - * with the other solver classes. */ struct AdditionalData - {}; + { + /** + * Constructor. By default, set the residual of the stopping criterion + * to the implicit residual. A true value of + * AdditionalData::use_default_residual refers to the + * implicit residual, while false reverts + * it. The former uses the result of the matrix-vector + * product already computed in other algorithm steps to derive the residual + * by a mere vector update, whereas the latter explicitly calculates the + * system residual with an additional matrix-vector product. More + * information on explicit and implicit residual stopping criteria can be + * found link + * here. + */ + explicit AdditionalData(const bool use_default_residual = true); + + /** + * Flag for the default residual that is used to measure convergence. + */ + bool use_default_residual; + }; + /** * Constructor. @@ -384,6 +444,12 @@ public: */ SolverFlexibleCG(SolverControl &cn, const AdditionalData &data = AdditionalData()); + +protected: + /** + * Additional parameters. + */ + AdditionalData additional_data; }; @@ -394,6 +460,12 @@ public: #ifndef DOXYGEN +template +DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) +inline SolverCG::AdditionalData::AdditionalData( + const bool use_default_residual) + : use_default_residual(use_default_residual) +{} template DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) @@ -476,6 +548,7 @@ inline void SolverCG::compute_eigs_and_cond( namespace internal { + namespace SolverCG { // This base class is used to select different variants of the conjugate @@ -495,11 +568,13 @@ namespace internal const PreconditionerType &preconditioner; const bool flexible; VectorType &x; + const VectorType &b; typename VectorMemory::Pointer r_pointer; typename VectorMemory::Pointer p_pointer; typename VectorMemory::Pointer v_pointer; typename VectorMemory::Pointer z_pointer; + typename VectorMemory::Pointer explicit_r_pointer; // Define some aliases for simpler access, using the variables 'r' for // the residual b - A*x, 'p' for the search direction, and 'v' for the @@ -511,39 +586,48 @@ namespace internal VectorType &p; VectorType &v; VectorType &z; + VectorType &explicit_r; + + Number r_dot_preconditioner_dot_r; + Number alpha; + Number beta; + double residual_norm; + Number previous_alpha; + const bool use_default_residual; - Number r_dot_preconditioner_dot_r; - Number alpha; - Number beta; - double residual_norm; - Number previous_alpha; IterationWorkerBase(const MatrixType &A, const PreconditionerType &preconditioner, const bool flexible, VectorMemory &memory, - VectorType &x) + VectorType &x, + const VectorType &b, + const bool &use_default_residual) : A(A) , preconditioner(preconditioner) , flexible(flexible) , x(x) + , b(b) , r_pointer(memory) , p_pointer(memory) , v_pointer(memory) , z_pointer(memory) + , explicit_r_pointer(memory) , r(*r_pointer) , p(*p_pointer) , v(*v_pointer) , z(*z_pointer) + , explicit_r(*explicit_r_pointer) , r_dot_preconditioner_dot_r(Number()) , alpha(Number()) , beta(Number()) , residual_norm(0.0) , previous_alpha(Number()) + , use_default_residual(use_default_residual) {} void - startup(const VectorType &b) + startup() { // Initialize without setting the vector entries, as those would soon // be overwritten anyway @@ -552,6 +636,8 @@ namespace internal v.reinit(x, true); if (flexible) z.reinit(x, true); + if (!use_default_residual) + explicit_r.reinit(x, true); // compute residual. if vector is zero, then short-circuit the full // computation @@ -581,22 +667,34 @@ namespace internal using BaseClass = IterationWorkerBase; + IterationWorker(const MatrixType &A, const PreconditionerType &preconditioner, const bool flexible, VectorMemory &memory, - VectorType &x) - : BaseClass(A, preconditioner, flexible, memory, x) + VectorType &x, + const VectorType &b, + const bool &use_default_residual) + : BaseClass(A, + preconditioner, + flexible, + memory, + x, + b, + use_default_residual) {} using BaseClass::A; using BaseClass::alpha; + using BaseClass::b; using BaseClass::beta; + using BaseClass::explicit_r; using BaseClass::p; using BaseClass::preconditioner; using BaseClass::r; using BaseClass::r_dot_preconditioner_dot_r; using BaseClass::residual_norm; + using BaseClass::use_default_residual; using BaseClass::v; using BaseClass::x; using BaseClass::z; @@ -645,7 +743,23 @@ namespace internal alpha = r_dot_preconditioner_dot_r / p_dot_A_dot_p; x.add(alpha, p); - residual_norm = std::sqrt(std::abs(r.add_and_dot(-alpha, v, r))); + + // compute the residual norm with implicit residual + if (use_default_residual) + { + residual_norm = std::sqrt(std::abs(r.add_and_dot(-alpha, v, r))); + } + // compute the residual norm with the explicit residual, i.e. + // compute l2 norm of Ax - b. + else + { + // compute the residual conjugate gradient update + r.add(-alpha, v); + // compute explicit residual + A.vmult(explicit_r, x); + explicit_r.add(-1, b); + residual_norm = explicit_r.l2_norm(); + } } void @@ -727,13 +841,17 @@ namespace internal const PreconditionerType &preconditioner, const bool flexible, VectorMemory &memory, - VectorType &x) + VectorType &x, + const VectorType &b, + const bool &use_default_residual) : IterationWorkerBase( A, preconditioner, flexible, memory, - x) + x, + b, + use_default_residual) , next_r_dot_preconditioner_dot_r(0.) , previous_beta(0.) {} @@ -1296,10 +1414,15 @@ void SolverCG::solve(const MatrixType &A, internal::SolverCG:: IterationWorker - worker( - A, preconditioner, determine_beta_by_flexible_formula, this->memory, x); + worker(A, + preconditioner, + determine_beta_by_flexible_formula, + this->memory, + x, + b, + additional_data.use_default_residual); - worker.startup(b); + worker.startup(); solver_state = this->iteration_status(0, worker.residual_norm, x); if (solver_state != SolverControl::iterate) @@ -1395,29 +1518,38 @@ boost::signals2::connection SolverCG::connect_eigenvalues_slot( +template +DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) +SolverFlexibleCG::AdditionalData::AdditionalData( + const bool use_default_residual) + : use_default_residual(use_default_residual) +{} + + template DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) SolverFlexibleCG::SolverFlexibleCG(SolverControl &cn, VectorMemory &mem, - const AdditionalData &) + const AdditionalData &data) : SolverCG(cn, mem) { this->determine_beta_by_flexible_formula = true; + this->additional_data = data; } template DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector) -SolverFlexibleCG::SolverFlexibleCG(SolverControl &cn, - const AdditionalData &) +SolverFlexibleCG::SolverFlexibleCG(SolverControl &cn, + const AdditionalData &data) : SolverCG(cn) { this->determine_beta_by_flexible_formula = true; + this->additional_data = data; } - #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/tests/lac/solver.cc b/tests/lac/solver.cc index 63dee311f4..910967b037 100644 --- a/tests/lac/solver.cc +++ b/tests/lac/solver.cc @@ -89,6 +89,8 @@ main() SolverControl control(100, 1.e-3); SolverControl verbose_control(100, 1.e-3, true); SolverCG<> cg(control, mem); + SolverCG<>::AdditionalData data0(false); + SolverCG<> cg_add_data(control, mem, data0); SolverGMRES<>::AdditionalData data1(6); SolverGMRES<> gmres(control, mem, data1); SolverGMRES<>::AdditionalData data2(6, true); @@ -170,6 +172,7 @@ main() control.set_max_steps(10); check_solve(cg, A, u, f, prec_no); + check_solve(cg_add_data, A, u, f, prec_no); check_solve(bicgstab, A, u, f, prec_no); check_solve(gmres, A, u, f, prec_no); check_solve(gmresright, A, u, f, prec_no); @@ -189,6 +192,7 @@ main() rich.set_omega(1. / A.diag_element(0)); check_solve(rich, A, u, f, prec_no); check_solve(cg, A, u, f, prec_no); + check_solve(cg_add_data, A, u, f, prec_no); check_solve(bicgstab, A, u, f, prec_no); check_solve(gmres, A, u, f, prec_no); check_solve(gmresright, A, u, f, prec_no); @@ -204,6 +208,7 @@ main() rich.set_omega(1. / A.diag_element(0)); check_solve(rich, A, u, f, prec_richardson); check_solve(cg, A, u, f, prec_richardson); + check_solve(cg_add_data, A, u, f, prec_richardson); check_solve(bicgstab, A, u, f, prec_richardson); check_solve(gmres, A, u, f, prec_richardson); check_solve(gmresright, A, u, f, prec_richardson); @@ -219,6 +224,7 @@ main() check_Tsolve(rich, A, u, f, prec_ssor); check_solve(rich, A, u, f, prec_ssor); check_solve(cg, A, u, f, prec_ssor); + check_solve(cg_add_data, A, u, f, prec_ssor); check_solve(bicgstab, A, u, f, prec_ssor); check_solve(gmres, A, u, f, prec_ssor); check_solve(gmresright, A, u, f, prec_ssor); @@ -233,6 +239,7 @@ main() check_Tsolve(rich, A, u, f, prec_sor); check_solve(rich, A, u, f, prec_sor); check_solve(cg, A, u, f, prec_sor); + check_solve(cg_add_data, A, u, f, prec_sor); check_solve(bicgstab, A, u, f, prec_sor); check_solve(gmres, A, u, f, prec_sor); check_solve(gmresright, A, u, f, prec_sor); @@ -246,6 +253,7 @@ main() check_Tsolve(rich, A, u, f, prec_psor); check_solve(rich, A, u, f, prec_psor); check_solve(cg, A, u, f, prec_psor); + check_solve(cg_add_data, A, u, f, prec_psor); check_solve(bicgstab, A, u, f, prec_psor); check_solve(gmres, A, u, f, prec_psor); check_solve(gmresright, A, u, f, prec_psor); diff --git a/tests/lac/solver.output b/tests/lac/solver.output index 4dd322bb37..2a6e1eeebd 100644 --- a/tests/lac/solver.output +++ b/tests/lac/solver.output @@ -3,12 +3,14 @@ DEAL::Size 4 Unknowns 9 DEAL::SOR-diff:0.000 DEAL:no-fail:cg::Starting value 3.000 DEAL:no-fail:cg::Convergence step 3 value 4.807e-17 +DEAL:no-fail:cg::Starting value 3.000 +DEAL:no-fail:cg::Convergence step 3 value 9.679e-16 DEAL:no-fail:Bicgstab::Starting value 3.000 DEAL:no-fail:Bicgstab::Convergence step 3 value 1.104e-17 DEAL:no-fail:GMRES::Starting value 3.000 -DEAL:no-fail:GMRES::Convergence step 3 value 4.551e-16 +DEAL:no-fail:GMRES::Convergence step 3 value 3.512e-09 DEAL:no-fail:GMRES::Starting value 3.000 -DEAL:no-fail:GMRES::Convergence step 3 value 4.551e-16 +DEAL:no-fail:GMRES::Convergence step 3 value 3.512e-09 DEAL:no-fail:GMRES::Starting value 3.000 DEAL:no-fail:GMRES::Convergence step 3 value 4.414e-16 DEAL:no-fail:SQMR::Starting value 3.000 @@ -19,12 +21,14 @@ DEAL:no:Richardson::Starting value 3.000 DEAL:no:Richardson::Convergence step 24 value 0.0007118 DEAL:no:cg::Starting value 3.000 DEAL:no:cg::Convergence step 3 value 4.807e-17 +DEAL:no:cg::Starting value 3.000 +DEAL:no:cg::Convergence step 3 value 9.679e-16 DEAL:no:Bicgstab::Starting value 3.000 DEAL:no:Bicgstab::Convergence step 3 value 1.104e-17 DEAL:no:GMRES::Starting value 3.000 -DEAL:no:GMRES::Convergence step 3 value 4.551e-16 +DEAL:no:GMRES::Convergence step 3 value 3.512e-09 DEAL:no:GMRES::Starting value 3.000 -DEAL:no:GMRES::Convergence step 3 value 4.551e-16 +DEAL:no:GMRES::Convergence step 3 value 3.512e-09 DEAL:no:GMRES::Starting value 3.000 DEAL:no:GMRES::Convergence step 3 value 4.414e-16 DEAL:no:SQMR::Starting value 3.000 @@ -35,14 +39,16 @@ DEAL:rich:Richardson::Starting value 3.000 DEAL:rich:Richardson::Convergence step 42 value 0.0008696 DEAL:rich:cg::Starting value 3.000 DEAL:rich:cg::Convergence step 3 value 2.871e-16 +DEAL:rich:cg::Starting value 3.000 +DEAL:rich:cg::Convergence step 3 value 1.691e-15 DEAL:rich:Bicgstab::Starting value 3.000 DEAL:rich:Bicgstab::Convergence step 3 value 3.238e-18 DEAL:rich:GMRES::Starting value 1.800 -DEAL:rich:GMRES::Convergence step 3 value 3.294e-16 +DEAL:rich:GMRES::Convergence step 3 value 3.512e-09 DEAL:rich:GMRES::Starting value 3.000 -DEAL:rich:GMRES::Convergence step 3 value 7.148e-16 +DEAL:rich:GMRES::Convergence step 3 value 4.139e-09 DEAL:rich:GMRES::Starting value 1.800 -DEAL:rich:GMRES::Convergence step 3 value 1.934e-16 +DEAL:rich:GMRES::Convergence step 3 value 1.976e-16 DEAL:rich:SQMR::Starting value 3.000 DEAL:rich:SQMR::Convergence step 3 value 1.422e-15 DEAL:rich:FIRE::Starting value 9.000 @@ -53,6 +59,8 @@ DEAL:ssor:Richardson::Starting value 3.000 DEAL:ssor:Richardson::Convergence step 7 value 0.0006128 DEAL:ssor:cg::Starting value 3.000 DEAL:ssor:cg::Convergence step 4 value 6.018e-05 +DEAL:ssor:cg::Starting value 3.000 +DEAL:ssor:cg::Convergence step 4 value 6.018e-05 DEAL:ssor:Bicgstab::Starting value 3.000 DEAL:ssor:Bicgstab::Convergence step 2 value 0.0003670 DEAL:ssor:GMRES::Starting value 1.920 @@ -72,12 +80,15 @@ DEAL:sor:Richardson::Convergence step 7 value 0.0004339 DEAL:sor:cg::Starting value 3.000 DEAL:sor:cg::Failure step 100 value 0.2585 DEAL:sor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) +DEAL:sor:cg::Starting value 3.000 +DEAL:sor:cg::Failure step 100 value 0.2585 +DEAL:sor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:sor:Bicgstab::Starting value 3.000 DEAL:sor:Bicgstab::Convergence step 5 value 4.331e-18 DEAL:sor:GMRES::Starting value 1.462 -DEAL:sor:GMRES::Convergence step 5 value 9.441e-19 +DEAL:sor:GMRES::Convergence step 5 value 3.488e-10 DEAL:sor:GMRES::Starting value 3.000 -DEAL:sor:GMRES::Convergence step 5 value 1.393e-17 +DEAL:sor:GMRES::Convergence step 5 value 1.520e-09 DEAL:sor:GMRES::Starting value 1.462 DEAL:sor:GMRES::Convergence step 21 value 0.0007019 DEAL:sor:FIRE::Starting value 9.000 @@ -89,6 +100,9 @@ DEAL:psor:Richardson::Convergence step 8 value 0.0004237 DEAL:psor:cg::Starting value 3.000 DEAL:psor:cg::Failure step 100 value 0.1024 DEAL:psor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) +DEAL:psor:cg::Starting value 3.000 +DEAL:psor:cg::Failure step 100 value 0.1024 +DEAL:psor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:psor:Bicgstab::Starting value 3.000 DEAL:psor:Bicgstab::Convergence step 4 value 0.0007969 DEAL:psor:GMRES::Starting value 1.467 @@ -104,6 +118,9 @@ DEAL::SOR-diff:0.000 DEAL:no-fail:cg::Starting value 11.00 DEAL:no-fail:cg::Failure step 10 value 0.1496 DEAL:no-fail::Exception: SolverControl::NoConvergence(it, worker.residual_norm) +DEAL:no-fail:cg::Starting value 11.00 +DEAL:no-fail:cg::Failure step 10 value 0.1496 +DEAL:no-fail::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:no-fail:Bicgstab::Starting value 11.00 DEAL:no-fail:Bicgstab::Failure step 10 value 0.002830 DEAL:no-fail:Bicgstab::Failure step 10 value 0.001961 @@ -128,6 +145,8 @@ DEAL:no:Richardson::Failure step 100 value 0.3002 DEAL:no::Exception: SolverControl::NoConvergence(iter, last_criterion) DEAL:no:cg::Starting value 11.00 DEAL:no:cg::Convergence step 15 value 0.0005794 +DEAL:no:cg::Starting value 11.00 +DEAL:no:cg::Convergence step 15 value 0.0005794 DEAL:no:Bicgstab::Starting value 11.00 DEAL:no:Bicgstab::Convergence step 11 value 0.0006357 DEAL:no:GMRES::Starting value 11.00 @@ -146,6 +165,8 @@ DEAL:rich:Richardson::Failure step 100 value 1.219 DEAL:rich::Exception: SolverControl::NoConvergence(iter, last_criterion) DEAL:rich:cg::Starting value 11.00 DEAL:rich:cg::Convergence step 15 value 0.0005794 +DEAL:rich:cg::Starting value 11.00 +DEAL:rich:cg::Convergence step 15 value 0.0005794 DEAL:rich:Bicgstab::Starting value 11.00 DEAL:rich:Bicgstab::Convergence step 11 value 0.0006357 DEAL:rich:GMRES::Starting value 6.600 @@ -164,6 +185,8 @@ DEAL:ssor:Richardson::Starting value 11.00 DEAL:ssor:Richardson::Convergence step 48 value 0.0009502 DEAL:ssor:cg::Starting value 11.00 DEAL:ssor:cg::Convergence step 8 value 0.0005816 +DEAL:ssor:cg::Starting value 11.00 +DEAL:ssor:cg::Convergence step 8 value 0.0005816 DEAL:ssor:Bicgstab::Starting value 11.00 DEAL:ssor:Bicgstab::Convergence step 5 value 0.0002498 DEAL:ssor:GMRES::Starting value 13.27 @@ -183,6 +206,9 @@ DEAL:sor:Richardson::Convergence step 88 value 0.0009636 DEAL:sor:cg::Starting value 11.00 DEAL:sor:cg::Failure step 100 value 5.643 DEAL:sor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) +DEAL:sor:cg::Starting value 11.00 +DEAL:sor:cg::Failure step 100 value 5.643 +DEAL:sor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:sor:Bicgstab::Starting value 11.00 DEAL:sor:Bicgstab::Convergence step 14 value 0.0009987 DEAL:sor:GMRES::Starting value 7.322 @@ -200,6 +226,9 @@ DEAL:psor:Richardson::Convergence step 89 value 0.0009736 DEAL:psor:cg::Starting value 11.00 DEAL:psor:cg::Failure step 100 value 2.935 DEAL:psor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) +DEAL:psor:cg::Starting value 11.00 +DEAL:psor:cg::Failure step 100 value 2.935 +DEAL:psor::Exception: SolverControl::NoConvergence(it, worker.residual_norm) DEAL:psor:Bicgstab::Starting value 11.00 DEAL:psor:Bicgstab::Convergence step 11 value 0.0005151 DEAL:psor:GMRES::Starting value 7.345 -- 2.39.5