SolverCG: Switch to more common variable names

diff --git a/include/deal.II/lac/solver_cg.h b/include/deal.II/lac/solver_cg.h
index d8ea571785c10f3a727de84e6ee2aa0f71e185bf..71edd410474fb02c621ab13a324710ebae9e1f74 100644 (file)
--- 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<double>.
  *
- * @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).
  *
  *
  * <h3>Eigenvalue computation</h3>
@@ -322,19 +321,22 @@ SolverCG<VectorType>::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<VectorType>::Pointer g_pointer(this->memory);
-  typename VectorMemory<VectorType>::Pointer d_pointer(this->memory);
-  typename VectorMemory<VectorType>::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<VectorType>::Pointer r_pointer(this->memory);
+  typename VectorMemory<VectorType>::Pointer p_pointer(this->memory);
+  typename VectorMemory<VectorType>::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<VectorType>::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<PreconditionerType, PreconditionIdentity>::value ==
+          false)
         {
-          if (std::is_same<PreconditionerType, PreconditionIdentity>::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<PreconditionerType, PreconditionIdentity>::value ? r : v;
+
+      if (it > 1)
         {
-          if (std::is_same<PreconditionerType, PreconditionIdentity>::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<VectorType>::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<VectorType>::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 655497573cbec6ff88d243289e086d4b0ad44e08..c40adef8e3a85f775e9bbfe199db491461d72e68 100644 (file)
--- 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 69e9273110c3795695d0a7f5b504daa5de5d3686..e7c652f327667ca6e8881b4e2914ec625420e976 100644 (file)
--- 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