From 27fbed3444498ebf6ebace3e698d36bbb21480cf Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 6 May 2022 11:49:19 +0200 Subject: [PATCH] SolverCG: Switch to more common variable names --- include/deal.II/lac/solver_cg.h | 126 ++++++++++++++------------------ tests/lac/solver.output | 62 ++++++++-------- tests/lapack/solver_cg.output | 8 +- 3 files changed, 91 insertions(+), 105 deletions(-) diff --git a/include/deal.II/lac/solver_cg.h b/include/deal.II/lac/solver_cg.h index d8ea571785..71edd41047 100644 --- a/include/deal.II/lac/solver_cg.h +++ b/include/deal.II/lac/solver_cg.h @@ -52,9 +52,8 @@ class PreconditionIdentity; * solution vector must be passed as template argument, and defaults to * dealii::Vector. * - * @note This version of CG is taken from D. Braess's book "Finite Elements". - * It requires a symmetric preconditioner (i.e., for example, SOR is not a - * possible choice). + * @note The CG method requires a symmetric preconditioner (i.e., for example, + * SOR is not a possible choice). * * *

Eigenvalue computation

@@ -322,19 +321,22 @@ SolverCG::solve(const MatrixType & A, { using number = typename VectorType::value_type; - SolverControl::State conv = SolverControl::iterate; + SolverControl::State solver_state = SolverControl::iterate; LogStream::Prefix prefix("cg"); // Memory allocation - typename VectorMemory::Pointer g_pointer(this->memory); - typename VectorMemory::Pointer d_pointer(this->memory); - typename VectorMemory::Pointer h_pointer(this->memory); - - // define some aliases for simpler access - VectorType &g = *g_pointer; - VectorType &d = *d_pointer; - VectorType &h = *h_pointer; + typename VectorMemory::Pointer r_pointer(this->memory); + typename VectorMemory::Pointer p_pointer(this->memory); + typename VectorMemory::Pointer v_pointer(this->memory); + + // 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 auxiliary + // vector. This naming convention is used e.g. by the description on + // https://en.wikipedia.org/wiki/Conjugate_gradient_method + VectorType &r = *r_pointer; + VectorType &p = *p_pointer; + VectorType &v = *v_pointer; // Should we build the matrix for eigenvalue computations? const bool do_eigenvalues = @@ -349,82 +351,68 @@ SolverCG::solve(const MatrixType & A, // resize the vectors, but do not set the values since they'd be overwritten // soon anyway. - g.reinit(x, true); - d.reinit(x, true); - h.reinit(x, true); + r.reinit(x, true); + p.reinit(x, true); + v.reinit(x, true); - int it = 0; - number gh = number(); - number beta = number(); - number alpha = number(); - number old_alpha = number(); + int it = 0; + number r_dot_preconditioner_dot_r = number(); + number beta = number(); + number alpha = number(); // compute residual. if vector is zero, then short-circuit the full // computation if (!x.all_zero()) { - A.vmult(g, x); - g.add(-1., b); + A.vmult(r, x); + r.sadd(-1., 1., b); } else - g.equ(-1., b); + r.equ(1., b); - double res = g.l2_norm(); - conv = this->iteration_status(0, res, x); - if (conv != SolverControl::iterate) + double residual_norm = r.l2_norm(); + solver_state = this->iteration_status(0, residual_norm, x); + if (solver_state != SolverControl::iterate) return; - while (conv == SolverControl::iterate) + while (solver_state == SolverControl::iterate) { it++; - old_alpha = alpha; + const number old_alpha = alpha; + const number old_r_dot_preconditioner_dot_r = r_dot_preconditioner_dot_r; - if (it > 1) + if (std::is_same::value == + false) { - if (std::is_same::value == - false) - { - preconditioner.vmult(h, g); - beta = gh; - Assert(std::abs(beta) != 0., ExcDivideByZero()); - gh = g * h; - beta = gh / beta; - d.sadd(beta, -1., h); - } - else - { - beta = gh; - gh = res * res; - beta = gh / beta; - d.sadd(beta, -1., g); - } + preconditioner.vmult(v, r); + r_dot_preconditioner_dot_r = r * v; } else + r_dot_preconditioner_dot_r = residual_norm * residual_norm; + + const VectorType &direction = + std::is_same::value ? r : v; + + if (it > 1) { - if (std::is_same::value == - false) - { - preconditioner.vmult(h, g); - d.equ(-1., h); - gh = g * h; - } - else - { - d.equ(-1., g); - gh = res * res; - } + Assert(std::abs(old_r_dot_preconditioner_dot_r) != 0., + ExcDivideByZero()); + beta = r_dot_preconditioner_dot_r / old_r_dot_preconditioner_dot_r; + p.sadd(beta, 1., direction); } + else + p.equ(1., direction); - A.vmult(h, d); + A.vmult(v, p); - alpha = d * h; - Assert(std::abs(alpha) != 0., ExcDivideByZero()); - alpha = gh / alpha; + const number p_dot_A_dot_p = p * v; + Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero()); + alpha = r_dot_preconditioner_dot_r / p_dot_A_dot_p; - x.add(alpha, d); - res = std::sqrt(std::abs(g.add_and_dot(alpha, h, g))); + x.add(alpha, p); + residual_norm = std::sqrt(std::abs(r.add_and_dot(-alpha, v, r))); - print_vectors(it, x, g, d); + print_vectors(it, x, r, p); if (it > 1) { @@ -443,7 +431,7 @@ SolverCG::solve(const MatrixType & A, all_condition_numbers_signal); } - conv = this->iteration_status(it, res, x); + solver_state = this->iteration_status(it, residual_norm, x); } compute_eigs_and_cond(diagonal, @@ -451,10 +439,8 @@ SolverCG::solve(const MatrixType & A, eigenvalues_signal, condition_number_signal); - // in case of failure: throw exception - if (conv != SolverControl::success) - AssertThrow(false, SolverControl::NoConvergence(it, res)); - // otherwise exit as normal + AssertThrow(solver_state == SolverControl::success, + SolverControl::NoConvergence(it, residual_norm)); } diff --git a/tests/lac/solver.output b/tests/lac/solver.output index 655497573c..c40adef8e3 100644 --- a/tests/lac/solver.output +++ b/tests/lac/solver.output @@ -2,43 +2,43 @@ 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::Convergence step 3 value 5.917e-17 DEAL:no-fail:Bicgstab::Starting value 3.000 -DEAL:no-fail:Bicgstab::Convergence step 3 value 1.104e-17 +DEAL:no-fail:Bicgstab::Convergence step 3 value 4.768e-18 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 6.160e-16 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 6.160e-16 DEAL:no-fail:SQMR::Starting value 3.000 -DEAL:no-fail:SQMR::Convergence step 3 value 1.280e-15 +DEAL:no-fail:SQMR::Convergence step 3 value 9.088e-16 DEAL:no-fail:FIRE::Starting value 9.000 DEAL:no-fail:FIRE::Convergence step 38 value 0.0002184 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::Convergence step 3 value 5.917e-17 DEAL:no:Bicgstab::Starting value 3.000 -DEAL:no:Bicgstab::Convergence step 3 value 1.104e-17 +DEAL:no:Bicgstab::Convergence step 3 value 4.768e-18 DEAL:no:GMRES::Starting value 3.000 -DEAL:no:GMRES::Convergence step 3 value 4.551e-16 +DEAL:no:GMRES::Convergence step 3 value 6.160e-16 DEAL:no:GMRES::Starting value 3.000 -DEAL:no:GMRES::Convergence step 3 value 4.551e-16 +DEAL:no:GMRES::Convergence step 3 value 6.160e-16 DEAL:no:SQMR::Starting value 3.000 -DEAL:no:SQMR::Convergence step 3 value 1.280e-15 +DEAL:no:SQMR::Convergence step 3 value 9.088e-16 DEAL:no:FIRE::Starting value 9.000 DEAL:no:FIRE::Convergence step 38 value 0.0002184 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::Convergence step 3 value 2.009e-16 DEAL:rich:Bicgstab::Starting value 3.000 -DEAL:rich:Bicgstab::Convergence step 3 value 3.238e-18 +DEAL:rich:Bicgstab::Convergence step 3 value 1.708e-17 DEAL:rich:GMRES::Starting value 1.800 -DEAL:rich:GMRES::Convergence step 3 value 3.294e-16 +DEAL:rich:GMRES::Convergence step 3 value 2.927e-16 DEAL:rich:GMRES::Starting value 3.000 -DEAL:rich:GMRES::Convergence step 3 value 7.148e-16 +DEAL:rich:GMRES::Convergence step 3 value 6.449e-16 DEAL:rich:SQMR::Starting value 3.000 -DEAL:rich:SQMR::Convergence step 3 value 1.422e-15 +DEAL:rich:SQMR::Convergence step 3 value 1.351e-15 DEAL:rich:FIRE::Starting value 9.000 DEAL:rich:FIRE::Convergence step 32 value 0.0009708 DEAL:ssor:RichardsonT::Starting value 3.000 @@ -63,13 +63,13 @@ DEAL:sor:Richardson::Starting value 3.000 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, res) +DEAL:sor::Exception: SolverControl::NoConvergence(it, residual_norm) DEAL:sor:Bicgstab::Starting value 3.000 -DEAL:sor:Bicgstab::Convergence step 5 value 4.331e-18 +DEAL:sor:Bicgstab::Convergence step 5 value 4.201e-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 2.605e-18 DEAL:sor:GMRES::Starting value 3.000 -DEAL:sor:GMRES::Convergence step 5 value 1.393e-17 +DEAL:sor:GMRES::Convergence step 5 value 5.190e-18 DEAL:sor:FIRE::Starting value 9.000 DEAL:sor:FIRE::Convergence step 51 value 0.0004993 DEAL:psor:RichardsonT::Starting value 3.000 @@ -78,7 +78,7 @@ DEAL:psor:Richardson::Starting value 3.000 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, res) +DEAL:psor::Exception: SolverControl::NoConvergence(it, 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 @@ -91,26 +91,26 @@ DEAL::Size 12 Unknowns 121 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, res) +DEAL:no-fail::Exception: SolverControl::NoConvergence(it, 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 -DEAL:no-fail::Exception: SolverControl::NoConvergence (state.last_step, state.last_residual) +DEAL:no-fail::Exception: SolverControl::NoConvergence(state.last_step, state.last_residual) DEAL:no-fail:GMRES::Starting value 11.00 DEAL:no-fail:GMRES::Failure step 10 value 0.8414 -DEAL:no-fail::Exception: SolverControl::NoConvergence (accumulated_iterations, last_res) +DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, last_res) DEAL:no-fail:GMRES::Starting value 11.00 DEAL:no-fail:GMRES::Failure step 10 value 0.8414 -DEAL:no-fail::Exception: SolverControl::NoConvergence (accumulated_iterations, last_res) +DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, last_res) DEAL:no-fail:SQMR::Starting value 11.00 DEAL:no-fail:SQMR::Failure step 10 value 0.4215 -DEAL:no-fail::Exception: SolverControl::NoConvergence (step, state.last_residual) +DEAL:no-fail::Exception: SolverControl::NoConvergence(step, state.last_residual) DEAL:no-fail:FIRE::Starting value 121.0 DEAL:no-fail:FIRE::Failure step 50 value 2.971 -DEAL:no-fail::Exception: SolverControl::NoConvergence (iter, gradients * gradients) +DEAL:no-fail::Exception: SolverControl::NoConvergence(iter, gradients * gradients) DEAL:no:Richardson::Starting value 11.00 DEAL:no:Richardson::Failure step 100 value 0.3002 -DEAL:no::Exception: SolverControl::NoConvergence (iter, last_criterion) +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:Bicgstab::Starting value 11.00 @@ -123,10 +123,10 @@ DEAL:no:SQMR::Starting value 11.00 DEAL:no:SQMR::Convergence step 16 value 0.0002583 DEAL:no:FIRE::Starting value 121.0 DEAL:no:FIRE::Failure step 100 value 0.003317 -DEAL:no::Exception: SolverControl::NoConvergence (iter, gradients * gradients) +DEAL:no::Exception: SolverControl::NoConvergence(iter, gradients * gradients) DEAL:rich:Richardson::Starting value 11.00 DEAL:rich:Richardson::Failure step 100 value 1.219 -DEAL:rich::Exception: SolverControl::NoConvergence (iter, last_criterion) +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:Bicgstab::Starting value 11.00 @@ -161,7 +161,7 @@ DEAL:sor:Richardson::Starting value 11.00 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, res) +DEAL:sor::Exception: SolverControl::NoConvergence(it, 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 @@ -176,7 +176,7 @@ DEAL:psor:Richardson::Starting value 11.00 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, res) +DEAL:psor::Exception: SolverControl::NoConvergence(it, 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 diff --git a/tests/lapack/solver_cg.output b/tests/lapack/solver_cg.output index 69e9273110..e7c652f327 100644 --- a/tests/lapack/solver_cg.output +++ b/tests/lapack/solver_cg.output @@ -2,15 +2,15 @@ DEAL::Size 4 Unknowns 9 DEAL:no-fail:cg::Starting value 3.000 DEAL:no-fail:cg::Condition number estimate: 3.535 -DEAL:no-fail:cg::Convergence step 3 value 0 +DEAL:no-fail:cg::Convergence step 3 value 5.917e-17 DEAL:no-fail:cg::Final Eigenvalues: 1.176 4.157 DEAL:no:cg::Starting value 3.000 DEAL:no:cg::Condition number estimate: 3.535 -DEAL:no:cg::Convergence step 3 value 0 +DEAL:no:cg::Convergence step 3 value 5.917e-17 DEAL:no:cg::Final Eigenvalues: 1.176 4.157 DEAL:rich:cg::Starting value 3.000 DEAL:rich:cg::Condition number estimate: 3.535 -DEAL:rich:cg::Convergence step 3 value 0 +DEAL:rich:cg::Convergence step 3 value 2.009e-16 DEAL:rich:cg::Final Eigenvalues: 0.7056 2.494 DEAL:ssor:cg::Starting value 3.000 DEAL:ssor:cg::Condition number estimate: 1.384 @@ -29,7 +29,7 @@ DEAL:no-fail:cg::Condition number estimate: 53.54 DEAL:no-fail:cg::Condition number estimate: 54.75 DEAL:no-fail:cg::Failure step 10 value 0.1496 DEAL:no-fail:cg::Final Eigenvalues: 0.1363 0.6565 1.499 2.357 3.131 3.994 5.392 6.576 7.464 -DEAL:no-fail::Exception: SolverControl::NoConvergence (it, res) +DEAL:no-fail::Exception: SolverControl::NoConvergence(it, residual_norm) DEAL:no:cg::Starting value 11.00 DEAL:no:cg::Condition number estimate: 11.01 DEAL:no:cg::Condition number estimate: 21.10 -- 2.39.5