]> https://gitweb.dealii.org/ - dealii.git/commitdiff
adding additional data for solver cg 17134/head
authorDario Coscia <dariocos99@gmail.com>
Tue, 3 Sep 2024 14:46:20 +0000 (16:46 +0200)
committerDario Coscia <dariocos99@gmail.com>
Thu, 5 Sep 2024 12:37:27 +0000 (14:37 +0200)
Co-authored-by: Martin Kronbichler <martin.kronbichler@rub.de>
richer docstring

include/deal.II/lac/solver_cg.h
tests/lac/solver.cc
tests/lac/solver.output

index dca6a25b68a7dbdf1859f3a5aef44111a0f1d8ec..9c2d48d8af45b0b8546eb5801c4b235b0a410f5d 100644 (file)
@@ -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<double>.
+ * dealii::Vector<double>. 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.
+ *
+ * <h3>Preconditioned residual</h3>
+ *
+ * @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 <tt>true</tt> value refers to the
+ * implicit residual, while <tt>false</tt> 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
+ * <a
+ * href="https://en.wikipedia.org/wiki/Conjugate_gradient_method#Explicit_residual_calculation">link
+ * here</a>.
  */
 template <typename VectorType = Vector<double>>
 DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
@@ -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 <tt>true</tt> value of
+     * AdditionalData::use_default_residual refers to the
+     * implicit residual, while <tt>false</tt> 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 <a
+     * href="https://en.wikipedia.org/wiki/Conjugate_gradient_method#Explicit_residual_calculation">link
+     * here</a>.
+     */
+    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 <tt>true</tt> value of
+     * AdditionalData::use_default_residual refers to the
+     * implicit residual, while <tt>false</tt> 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 <a
+     * href="https://en.wikipedia.org/wiki/Conjugate_gradient_method#Explicit_residual_calculation">link
+     * here</a>.
+     */
+    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 <typename VectorType>
+DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
+inline SolverCG<VectorType>::AdditionalData::AdditionalData(
+  const bool use_default_residual)
+  : use_default_residual(use_default_residual)
+{}
 
 template <typename VectorType>
 DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
@@ -476,6 +548,7 @@ inline void SolverCG<VectorType>::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<VectorType>::Pointer r_pointer;
       typename VectorMemory<VectorType>::Pointer p_pointer;
       typename VectorMemory<VectorType>::Pointer v_pointer;
       typename VectorMemory<VectorType>::Pointer z_pointer;
+      typename VectorMemory<VectorType>::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<VectorType> &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<VectorType, MatrixType, PreconditionerType>;
 
+
       IterationWorker(const MatrixType         &A,
                       const PreconditionerType &preconditioner,
                       const bool                flexible,
                       VectorMemory<VectorType> &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<VectorType> &memory,
-                      VectorType               &x)
+                      VectorType               &x,
+                      const VectorType         &b,
+                      const bool               &use_default_residual)
         : IterationWorkerBase<VectorType, MatrixType, PreconditionerType>(
             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<VectorType>::solve(const MatrixType         &A,
 
   internal::SolverCG::
     IterationWorker<VectorType, MatrixType, PreconditionerType>
-      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<VectorType>::connect_eigenvalues_slot(
 
 
 
+template <typename VectorType>
+DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
+SolverFlexibleCG<VectorType>::AdditionalData::AdditionalData(
+  const bool use_default_residual)
+  : use_default_residual(use_default_residual)
+{}
+
+
 template <typename VectorType>
 DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
 SolverFlexibleCG<VectorType>::SolverFlexibleCG(SolverControl            &cn,
                                                VectorMemory<VectorType> &mem,
-                                               const AdditionalData &)
+                                               const AdditionalData     &data)
   : SolverCG<VectorType>(cn, mem)
 {
   this->determine_beta_by_flexible_formula = true;
+  this->additional_data                    = data;
 }
 
 
 
 template <typename VectorType>
 DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
-SolverFlexibleCG<VectorType>::SolverFlexibleCG(SolverControl &cn,
-                                               const AdditionalData &)
+SolverFlexibleCG<VectorType>::SolverFlexibleCG(SolverControl        &cn,
+                                               const AdditionalData &data)
   : SolverCG<VectorType>(cn)
 {
   this->determine_beta_by_flexible_formula = true;
+  this->additional_data                    = data;
 }
 
 
-
 #endif // DOXYGEN
 
 DEAL_II_NAMESPACE_CLOSE
index 63dee311f42d9c7bebfb04fb7f3b36391c95dcc9..910967b037c40d2b181f6aef7cc3c9bf0e5e80b6 100644 (file)
@@ -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);
index 4dd322bb37a372271d48995c20a4f02460304935..2a6e1eeebd6b06acbb3b0fe3f5a110ac47fa7ef2 100644 (file)
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.