]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make interface to different variants more slim by moving iteration
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Mon, 9 May 2022 08:57:11 +0000 (10:57 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Mon, 9 May 2022 19:52:34 +0000 (21:52 +0200)
to separate class

doc/doxygen/references.bib
include/deal.II/lac/diagonal_matrix.h
include/deal.II/lac/solver_cg.h
tests/lac/solver.output
tests/lapack/solver_cg.output

index 3f9e730389f88f3666637387856d5496fa84d2d0..afa41bf28a9eb951a85846b91e42c28d763e0f1a 100644 (file)
@@ -1787,3 +1787,15 @@ url = {https://doi.org/10.1016/0045-7930(73)90027-3}
   doi    = {10.11588/heidok.00005743},
   url    = {https://doi.org/10.11588/heidok.00005743}
 }
+
+@article{Chronopoulos1989,
+  author = {Chronopoulos, A. T. and Gear, C. W.},
+  title = {S-step Iterative Methods for Symmetric Linear Systems},
+  journal = {Journal of Computational and Applied Mathematics},
+  volume = {25},
+  number = {2},
+  year = {1989},
+  pages = {153--168},
+  doi = {10.1016/0377-0427(89)90045-9},
+  url = {https://doi.org/10.1016/0377-0427(89)90045-9}
+}
index 3050ce331232ed59acad14f559b755711120f1da..46ffbf76f5738d8fc083eadd4a47e8ec1a96203d 100644 (file)
@@ -197,13 +197,16 @@ public:
   Tvmult_add(VectorType &dst, const VectorType &src) const;
 
   /**
-   * Apply the preconditioner only on a subrange of elements on the vector.
+   * Apply the preconditioner only to a subrange of elements of the given
+   * vector. To support this operation, the given `VectorType` template
+   * argument needs to support a method `begin()` to return the pointer to the
+   * start of the stored elements.
    */
   void
-  apply_on_subrange(const std::size_t index_of_first_unknown,
-                    const std::size_t length,
-                    const value_type *src,
-                    value_type *      dst) const;
+  apply_to_subrange(const unsigned int index_of_first_unknown,
+                    const unsigned int length,
+                    const value_type * src,
+                    value_type *       dst) const;
 
   /**
    * Initialize vector @p dst to have the same size and partition as
@@ -415,11 +418,11 @@ DiagonalMatrix<VectorType>::Tvmult_add(VectorType &      dst,
 
 template <typename VectorType>
 void
-DiagonalMatrix<VectorType>::apply_on_subrange(
-  const std::size_t index_of_first_unknown,
-  const std::size_t length,
-  const value_type *src,
-  value_type *      dst) const
+DiagonalMatrix<VectorType>::apply_to_subrange(
+  const unsigned int index_of_first_unknown,
+  const unsigned int length,
+  const value_type * src,
+  value_type *       dst) const
 {
   AssertIndexRange(index_of_first_unknown,
                    diagonal.locally_owned_elements().n_elements());
@@ -429,7 +432,7 @@ DiagonalMatrix<VectorType>::apply_on_subrange(
   const value_type *diagonal_entry = diagonal.begin() + index_of_first_unknown;
 
   DEAL_II_OPENMP_SIMD_PRAGMA
-  for (std::size_t i = 0; i < length; ++i)
+  for (unsigned int i = 0; i < length; ++i)
     dst[i] = diagonal_entry[i] * src[i];
 }
 
index 269d99acd76d95f90fcc1353b7ededc754804b98..f227e1830151d04b7ac55cfdfd235d08ebabaaaa 100644 (file)
@@ -34,6 +34,14 @@ DEAL_II_NAMESPACE_OPEN
 // forward declaration
 #ifndef DOXYGEN
 class PreconditionIdentity;
+namespace LinearAlgebra
+{
+  namespace distributed
+  {
+    template <typename, typename>
+    class Vector;
+  }
+} // namespace LinearAlgebra
 #endif
 
 
@@ -93,6 +101,71 @@ class PreconditionIdentity;
  * The solve() function of this class uses the mechanism described in the
  * Solver base class to determine convergence. This mechanism can also be used
  * to observe the progress of the iteration.
+ *
+ * <h4>Optimized operations with specific `MatrixType` argument</h4>
+ *
+ * This class enables to embed the vector updates into the matrix-vector
+ * product in case the `MatrixType` and `PreconditionerType` support such a
+ * mode of operation. To this end, the `VectorType` needs to be
+ * LinearAlgebra::distributed::Vector, the class `MatrixType` needs to provide
+ * a function with the signature
+ * @code
+ * void MatrixType::vmult(
+ *    VectorType &,
+ *    const VectorType &,
+ *    const std::function<void(const unsigned int, const unsigned int)> &,
+ *    const std::function<void(const unsigned int, const unsigned int)> &) const
+ * @endcode
+ * where the two given functions run before and after the matrix-vector
+ * product, respectively, and the `PreconditionerType` needs to provide a
+ * function with the signature
+ * @code
+ * void PreconditionerType::apply_to_subrange(unsigned int start_range,
+ *                                            unsigned int end_range,
+ *                                            const Number* src_ptr,
+ *                                            Number* dst_ptr)
+ * @endcode
+
+ * where the pointers `src_ptr` and `dst_ptr` point to the location in the
+ * vector where the operation should be applied to. The functions passed to
+ * `MatrixType::vmult` take as arguments a sub-range among the locally owned
+ * elements of the vector, defined as half-open intervals. The intervals are
+ * designed to be scheduled close to the time the matrix-vector product
+ * touches those entries in the `src` and `dst` vectors, respectively, with
+ * the requirement that
+ * <ul>
+ * <li> the matrix-vector product may only access an entry in `src` or `dst`
+ * once the `operation_before_matrix_vector_product` has been run on that
+ * vector entry; </li>
+ * <li> `operation_after_matrix_vector_product` may run on a range of entries
+ * `[i,j)` once the matrix-vector product does not access the entries `[i,j)`
+ * in `src` and `dst` any more. </li>
+ * </ul>
+ * The motivation for this function is to increase data locality and hence
+ * cache usage. For the example of a class similar to the one in the step-37
+ * tutorial program, the implementation is
+ * @code
+ * void
+ * vmult(LinearAlgebra::distributed::Vector<number> &      dst,
+ *       const LinearAlgebra::distributed::Vector<number> &src,
+ *       const std::function<void(const unsigned int, const unsigned int)>
+ *         &operation_before_matrix_vector_product,
+ *       const std::function<void(const unsigned int, const unsigned int)>
+ *         &operation_after_matrix_vector_product) const
+ * {
+ *   data.cell_loop(&LaplaceOperator::local_apply,
+ *                  this,
+ *                  dst,
+ *                  src,
+ *                  operation_before_matrix_vector_product,
+ *                  operation_after_matrix_vector_product);
+ * }
+ * @endcode
+ *
+ * In terms of the SolverCG implementation, the operation before the loop will
+ * 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.
  */
 template <typename VectorType = Vector<double>>
 class SolverCG : public SolverBase<VectorType>
@@ -432,8 +505,8 @@ namespace internal
       template <typename PreconditionerType2>
       static decltype(
         std::declval<PreconditionerType2 const>().apply_to_subrange(
-          std::size_t(),
-          std::size_t(),
+          0U,
+          0U,
           std::declval<const typename PreconditionerType2::value_type *>(),
           std::declval<typename PreconditionerType2::value_type *>()))
       detect_preconditioner(const PreconditionerType2 &);
@@ -448,238 +521,417 @@ namespace internal
                       bool>::value &&
         !std::is_same<decltype(detect_preconditioner(
                         std::declval<PreconditionerType>())),
-                      bool>::value;
+                      bool>::value &&
+        std::is_same<
+          VectorType,
+          LinearAlgebra::distributed::Vector<typename VectorType::value_type,
+                                             MemorySpace::Host>>::value;
     };
 
+
+
     // We need to have a separate declaration for static const members
     template <typename T, typename U, typename V>
     const bool supports_vmult_with_std_functions<T, U, V>::value;
 
-    /**
-     * Internal function to run one iteration of the conjugate gradient solver
-     * for standard matrix and preconditioner arguments.
-     */
-    template <typename Number,
-              typename VectorType,
+
+
+    // Internal class to run one iteration of the conjugate gradient solver
+    // for standard matrix and preconditioner arguments.
+    template <typename VectorType,
               typename MatrixType,
-              typename PreconditionerType,
-              typename std::enable_if<
-                !supports_vmult_with_std_functions<MatrixType,
-                                                   VectorType,
-                                                   PreconditionerType>::value,
-                int>::type * = nullptr>
-    void
-    do_cg_iteration(const MatrixType &        A,
-                    const PreconditionerType &preconditioner,
-                    const typename dealii::SolverCG<VectorType>::AdditionalData
-                      &                additional_data,
-                    const unsigned int iteration,
-                    VectorType &       x,
-                    VectorType &       r,
-                    VectorType &       p,
-                    VectorType &       v,
-                    VectorType &       z,
-                    Number &           r_dot_preconditioner_dot_r,
-                    Number &           alpha,
-                    Number &           beta,
-                    double &           residual_norm,
-                    const Number /*old_alpha*/,
-                    const Number /*old_beta*/)
+              typename PreconditionerType>
+    struct IterationWorkerBase
     {
-      const Number old_r_dot_preconditioner_dot_r = r_dot_preconditioner_dot_r;
+      using Number = typename VectorType::value_type;
+
+      const MatrixType &        A;
+      const PreconditionerType &preconditioner;
+      const bool                flexible;
+      VectorType &              x;
+
+      typename VectorMemory<VectorType>::Pointer r_pointer;
+      typename VectorMemory<VectorType>::Pointer p_pointer;
+      typename VectorMemory<VectorType>::Pointer v_pointer;
+      typename VectorMemory<VectorType>::Pointer z_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
+      // auxiliary vector. This naming convention is used e.g. by the
+      // description on
+      // https://en.wikipedia.org/wiki/Conjugate_gradient_method. The variable
+      // 'z' gets only used for the flexible variant of the CG method.
+      VectorType &r;
+      VectorType &p;
+      VectorType &v;
+      VectorType &z;
+
+      Number r_dot_preconditioner_dot_r;
+      Number alpha;
+      Number beta;
+      double residual_norm;
+      Number previous_alpha;
+      Number previous_beta;
+
+      IterationWorkerBase(const MatrixType &        A,
+                          const PreconditionerType &preconditioner,
+                          const bool                flexible,
+                          VectorMemory<VectorType> &memory,
+                          VectorType &              x)
+        : A(A)
+        , preconditioner(preconditioner)
+        , flexible(flexible)
+        , x(x)
+        , r_pointer(memory)
+        , p_pointer(memory)
+        , v_pointer(memory)
+        , z_pointer(memory)
+        , r(*r_pointer)
+        , p(*p_pointer)
+        , v(*v_pointer)
+        , z(*z_pointer)
+        , r_dot_preconditioner_dot_r(Number())
+        , alpha(Number())
+        , beta(Number())
+        , residual_norm(0.0)
+        , previous_alpha(Number())
+        , previous_beta(Number())
+      {}
+
+      void
+      startup(const VectorType &b)
+      {
+        // Initialize without setting the vector entries, as those would soon
+        // be overwritten anyway
+        r.reinit(x, true);
+        p.reinit(x, true);
+        v.reinit(x, true);
+        if (flexible)
+          z.reinit(x, true);
+
+        // compute residual. if vector is zero, then short-circuit the full
+        // computation
+        if (!x.all_zero())
+          {
+            A.vmult(r, x);
+            r.sadd(-1., 1., b);
+          }
+        else
+          r.equ(1., b);
 
-      if (std::is_same<PreconditionerType, PreconditionIdentity>::value ==
-          false)
-        {
-          preconditioner.vmult(v, r);
-          r_dot_preconditioner_dot_r = r * v;
-        }
-      else
-        r_dot_preconditioner_dot_r = residual_norm * residual_norm;
+        residual_norm = r.l2_norm();
+      }
 
-      const VectorType &direction =
-        std::is_same<PreconditionerType, PreconditionIdentity>::value ? r : v;
+      void
+      do_iteration(const unsigned int iteration_index)
+      {
+        const Number previous_r_dot_preconditioner_dot_r =
+          r_dot_preconditioner_dot_r;
+        previous_alpha = alpha;
+        previous_beta  = beta;
 
-      if (iteration > 1)
-        {
-          Assert(std::abs(old_r_dot_preconditioner_dot_r) != 0.,
-                 ExcDivideByZero());
-          beta = r_dot_preconditioner_dot_r / old_r_dot_preconditioner_dot_r;
-          if (additional_data.use_flexible_variant)
-            beta -= (r * z) / old_r_dot_preconditioner_dot_r;
-          p.sadd(beta, 1., direction);
-        }
-      else
-        p.equ(1., direction);
+        if (std::is_same<PreconditionerType, PreconditionIdentity>::value ==
+            false)
+          {
+            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 (additional_data.use_flexible_variant)
-        z.swap(v);
+        if (iteration_index > 1)
+          {
+            Assert(std::abs(previous_r_dot_preconditioner_dot_r) != 0.,
+                   ExcDivideByZero());
+            beta =
+              r_dot_preconditioner_dot_r / previous_r_dot_preconditioner_dot_r;
+            if (flexible)
+              beta -= (r * z) / previous_r_dot_preconditioner_dot_r;
+            p.sadd(beta, 1., direction);
+          }
+        else
+          p.equ(1., direction);
 
-      A.vmult(v, p);
+        if (flexible)
+          z.swap(v);
 
-      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;
+        A.vmult(v, p);
 
-      x.add(alpha, p);
-      residual_norm = std::sqrt(std::abs(r.add_and_dot(-alpha, v, r)));
-    }
+        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;
 
-    /**
-     * Internal function to run one iteration of the conjugate gradient solver
-     * for matrices and preconditioners that support interleaving the vector
-     * updates with the matrix-vector product.
-     */
-    template <typename Number,
-              typename VectorType,
+        x.add(alpha, p);
+        residual_norm = std::sqrt(std::abs(r.add_and_dot(-alpha, v, r)));
+      }
+
+      void
+      finalize_after_convergence(const unsigned int)
+      {}
+    };
+
+
+
+    // Actual class with the basic operation implemented in
+    // IterationWorkerBase
+    template <typename VectorType,
               typename MatrixType,
               typename PreconditionerType,
-              typename std::enable_if<
-                supports_vmult_with_std_functions<MatrixType,
-                                                  VectorType,
-                                                  PreconditionerType>::value,
-                int>::type * = nullptr>
-    void
-    do_cg_iteration(
-      const MatrixType &        A,
-      const PreconditionerType &preconditioner,
-      const typename dealii::SolverCG<VectorType>::AdditionalData &,
-      const unsigned int iteration,
-      VectorType &       x_vector,
-      VectorType &       r_vector,
-      VectorType &       p_vector,
-      VectorType &       v_vector,
-      VectorType &,
-      Number &     r_dot_preconditioner_dot_r,
-      Number &     alpha,
-      Number &     beta,
-      double &     residual_norm,
-      const Number old_alpha,
-      const Number old_beta)
+              typename = int>
+    struct IterationWorker
+      : public IterationWorkerBase<VectorType, MatrixType, PreconditionerType>
     {
-      const auto operation_before_loop = [&](const unsigned int start_range,
-                                             const unsigned int end_range) {
-        Number *                       x = x_vector.begin() + start_range;
-        Number *                       r = r_vector.begin() + start_range;
-        Number *                       p = p_vector.begin() + start_range;
-        Number *                       v = v_vector.begin() + start_range;
-        constexpr unsigned int         grain_size = 32;
-        std::array<Number, grain_size> prec_r;
-        if (iteration == 1)
-          {
-            for (unsigned int j = start_range; j < end_range; j += grain_size)
-              {
-                const unsigned int length = std::min(grain_size, end_range - j);
-                preconditioner.apply_on_subrange(j, length, r, prec_r.data());
-                DEAL_II_OPENMP_SIMD_PRAGMA
-                for (unsigned int i = 0; i < length; ++i)
-                  {
-                    p[i] = prec_r[i];
-                    v[i] = Number();
-                  }
-                p += length;
-                r += length;
-                v += length;
-              }
-          }
-        else if (iteration % 2 == 0)
-          {
-            for (unsigned int j = start_range; j < end_range; j += grain_size)
-              {
-                const unsigned int length = std::min(grain_size, end_range - j);
-                DEAL_II_OPENMP_SIMD_PRAGMA
-                for (unsigned int i = 0; i < length; ++i)
-                  r[i] -= alpha * v[i];
-                preconditioner.apply_on_subrange(j, length, r, prec_r.data());
-                DEAL_II_OPENMP_SIMD_PRAGMA
-                for (unsigned int i = 0; i < length; ++i)
-                  {
-                    p[i] = beta * p[i] + prec_r[i];
-                    v[i] = Number();
-                  }
-                p += length;
-                r += length;
-                v += length;
-              }
-          }
+      IterationWorker(const MatrixType &        A,
+                      const PreconditionerType &preconditioner,
+                      const bool                flexible,
+                      VectorMemory<VectorType> &memory,
+                      VectorType &              x)
+        : IterationWorkerBase<VectorType, MatrixType, PreconditionerType>(
+            A,
+            preconditioner,
+            flexible,
+            memory,
+            x)
+      {}
+    };
+
+
+
+    // Internal function to run one iteration of the conjugate gradient solver
+    // for matrices and preconditioners that support interleaving the vector
+    // updates with the matrix-vector product.
+    template <typename VectorType,
+              typename MatrixType,
+              typename PreconditionerType>
+    struct IterationWorker<
+      VectorType,
+      MatrixType,
+      PreconditionerType,
+      typename std::enable_if<
+        supports_vmult_with_std_functions<MatrixType,
+                                          VectorType,
+                                          PreconditionerType>::value,
+        int>::type>
+      : public IterationWorkerBase<VectorType, MatrixType, PreconditionerType>
+    {
+      static constexpr unsigned int grain_size = 32;
+
+      IterationWorker(const MatrixType &        A,
+                      const PreconditionerType &preconditioner,
+                      const bool                flexible,
+                      VectorMemory<VectorType> &memory,
+                      VectorType &              x)
+        : IterationWorkerBase<VectorType, MatrixType, PreconditionerType>(
+            A,
+            preconditioner,
+            flexible,
+            memory,
+            x)
+      {}
+
+      void
+      do_iteration(const unsigned int iteration_index)
+      {
+        using Number = typename VectorType::value_type;
+
+        const auto operation_before_loop = [&](const unsigned int start_range,
+                                               const unsigned int end_range) {
+          Number *                       x = this->x.begin() + start_range;
+          Number *                       r = this->r.begin() + start_range;
+          Number *                       p = this->p.begin() + start_range;
+          Number *                       v = this->v.begin() + start_range;
+          std::array<Number, grain_size> prec_r;
+          if (iteration_index == 1)
+            {
+              for (unsigned int j = start_range; j < end_range; j += grain_size)
+                {
+                  const unsigned int length =
+                    std::min(grain_size, end_range - j);
+                  this->preconditioner.apply_to_subrange(j,
+                                                         length,
+                                                         r,
+                                                         prec_r.data());
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int i = 0; i < length; ++i)
+                    {
+                      p[i] = prec_r[i];
+                      v[i] = Number();
+                    }
+                  p += length;
+                  r += length;
+                  v += length;
+                }
+            }
+          else if (iteration_index % 2 == 0)
+            {
+              for (unsigned int j = start_range; j < end_range; j += grain_size)
+                {
+                  const unsigned int length =
+                    std::min(grain_size, end_range - j);
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int i = 0; i < length; ++i)
+                    r[i] -= this->alpha * v[i];
+                  this->preconditioner.apply_to_subrange(j,
+                                                         length,
+                                                         r,
+                                                         prec_r.data());
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int i = 0; i < length; ++i)
+                    {
+                      p[i] = this->beta * p[i] + prec_r[i];
+                      v[i] = Number();
+                    }
+                  p += length;
+                  r += length;
+                  v += length;
+                }
+            }
+          else
+            {
+              const Number alpha_plus_previous_alpha_over_beta =
+                this->alpha + this->previous_alpha / this->previous_beta;
+              const Number previous_alpha_over_beta =
+                this->previous_alpha / this->previous_beta;
+              for (unsigned int j = start_range; j < end_range; j += grain_size)
+                {
+                  const unsigned int length =
+                    std::min(grain_size, end_range - j);
+                  this->preconditioner.apply_to_subrange(j,
+                                                         length,
+                                                         r,
+                                                         prec_r.data());
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int i = 0; i < length; ++i)
+                    {
+                      x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
+                              previous_alpha_over_beta * prec_r[i];
+                      r[i] -= this->alpha * v[i];
+                    }
+                  this->preconditioner.apply_to_subrange(j,
+                                                         length,
+                                                         r,
+                                                         prec_r.data());
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int i = 0; i < length; ++i)
+                    {
+                      p[i] = this->beta * p[i] + prec_r[i];
+                      v[i] = Number();
+                    }
+                  p += length;
+                  r += length;
+                  v += length;
+                  x += length;
+                }
+            }
+        };
+
+        std::array<Number, 7> local_sums = {};
+        const auto operation_after_loop  = [&](const unsigned int start_range,
+                                              const unsigned int end_range) {
+          const Number *                 x = this->x.begin() + start_range;
+          const Number *                 r = this->r.begin() + start_range;
+          const Number *                 p = this->p.begin() + start_range;
+          const Number *                 v = this->v.begin() + start_range;
+          std::array<Number, grain_size> prec_r;
+          std::array<Number, grain_size> prec_v;
+          for (unsigned int j = start_range; j < end_range; j += grain_size)
+            {
+              const unsigned int length = std::min(grain_size, end_range - j);
+              this->preconditioner.apply_to_subrange(j,
+                                                     length,
+                                                     r,
+                                                     prec_r.data());
+              this->preconditioner.apply_to_subrange(j,
+                                                     length,
+                                                     v,
+                                                     prec_v.data());
+              for (unsigned int i = 0; i < length; ++i)
+                {
+                  local_sums[0] += p[i] * v[i];
+                  local_sums[1] += v[i] * v[i];
+                  local_sums[2] += r[i] * v[i];
+                  local_sums[3] += r[i] * r[i];
+                  local_sums[4] += r[i] * prec_v[i];
+                  local_sums[5] += v[i] * prec_v[i];
+                  local_sums[6] += r[i] * prec_r[i];
+                }
+              p += length;
+              v += length;
+              r += length;
+            }
+        };
+
+        this->A.vmult(this->v,
+                      this->p,
+                      operation_before_loop,
+                      operation_after_loop);
+
+        Utilities::MPI::sum(dealii::ArrayView<const Number>(local_sums.data(),
+                                                            7),
+                            this->r.get_mpi_communicator(),
+                            dealii::ArrayView<Number>(local_sums.data(), 7));
+
+        this->previous_alpha = this->alpha;
+        this->previous_beta  = this->beta;
+
+        const Number p_dot_A_dot_p = local_sums[0];
+        Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero());
+
+        const Number previous_r_dot_preconditioner_dot_r = local_sums[6];
+        this->alpha = previous_r_dot_preconditioner_dot_r / p_dot_A_dot_p;
+        this->residual_norm = std::sqrt(
+          local_sums[3] +
+          this->alpha * (-2. * local_sums[2] + this->alpha * local_sums[1]));
+
+        this->r_dot_preconditioner_dot_r =
+          previous_r_dot_preconditioner_dot_r +
+          this->alpha * (-2. * local_sums[4] + this->alpha * local_sums[5]);
+
+        this->beta = this->r_dot_preconditioner_dot_r /
+                     previous_r_dot_preconditioner_dot_r;
+      }
+
+      void
+      finalize_after_convergence(const unsigned int iteration_index)
+      {
+        if (iteration_index % 2 == 1)
+          this->x.add(this->alpha, this->p);
         else
           {
-            const Number alpha_plus_alpha_old = alpha + old_alpha / old_beta;
-            const Number alpha_old_beta_old   = old_alpha / old_beta;
-            for (unsigned int j = start_range; j < end_range; j += grain_size)
+            using Number                 = typename VectorType::value_type;
+            const unsigned int end_range = this->x.locally_owned_size();
+
+            Number *     x = this->x.begin();
+            Number *     r = this->r.begin();
+            Number *     p = this->p.begin();
+            const Number alpha_plus_previous_alpha_over_beta =
+              this->alpha + this->previous_alpha / this->previous_beta;
+            const Number previous_alpha_over_beta =
+              this->previous_alpha / this->previous_beta;
+
+            std::array<Number, grain_size> prec_r;
+            for (unsigned int j = 0; j < end_range; j += grain_size)
               {
                 const unsigned int length = std::min(grain_size, end_range - j);
-                preconditioner.apply_on_subrange(j, length, r, prec_r.data());
+                this->preconditioner.apply_to_subrange(j,
+                                                       length,
+                                                       r,
+                                                       prec_r.data());
                 DEAL_II_OPENMP_SIMD_PRAGMA
                 for (unsigned int i = 0; i < length; ++i)
-                  {
-                    x[i] += alpha_plus_alpha_old * p[i] +
-                            alpha_old_beta_old * prec_r[i];
-                    r[i] -= alpha * v[i];
-                  }
-                preconditioner.apply_on_subrange(j, length, r, prec_r.data());
-                DEAL_II_OPENMP_SIMD_PRAGMA
-                for (unsigned int i = 0; i < length; ++i)
-                  {
-                    p[i] = beta * p[i] + prec_r[i];
-                    v[i] = 0.;
-                  }
-                p += length;
-                r += length;
-                v += length;
+                  x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
+                          previous_alpha_over_beta * prec_r[i];
+
                 x += length;
+                r += length;
+                p += length;
               }
           }
-      };
-
-      std::array<Number, 7> local_sums = {};
-      const auto operation_after_loop  = [&](const unsigned int start_range,
-                                            const unsigned int end_range) {
-        Number *                       x = x_vector.begin() + start_range;
-        Number *                       r = r_vector.begin() + start_range;
-        Number *                       p = p_vector.begin() + start_range;
-        Number *                       v = v_vector.begin() + start_range;
-        constexpr unsigned int         grain_size = 32;
-        std::array<Number, grain_size> prec_r;
-        std::array<Number, grain_size> prec_v;
-        for (unsigned int j = start_range; j < end_range; j += grain_size)
-          {
-            const unsigned int length = std::min(grain_size, end_range - j);
-            preconditioner.apply_on_subrange(j, length, r, prec_r.data());
-            preconditioner.apply_on_subrange(j, length, v, prec_v.data());
-            for (unsigned int i = 0; i < length; ++i)
-              {
-                local_sums[0] += p[i] * v[i];
-                local_sums[1] += v[i] * v[i];
-                local_sums[2] += r[i] * v[i];
-                local_sums[3] += r[i] * r[i];
-                local_sums[4] += r[i] * prec_v[i];
-                local_sums[5] += v[i] * prec_v[i];
-                local_sums[6] += r[i] * prec_r[i];
-              }
-          }
-      };
-
-      A.vmult(v_vector, p_vector, operation_before_loop, operation_after_loop);
-
-      Utilities::MPI::sum(dealii::ArrayView<const double>(local_sums.begin(),
-                                                          7),
-                          r_vector.get_mpi_communicator(),
-                          dealii::ArrayView<double>(local_sums.begin_raw(), 7));
-
-      const Number p_dot_A_dot_p = local_sums[0];
-      Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero());
-
-      const Number old_r_dot_preconditioner_dot_r = local_sums[6];
-      alpha         = old_r_dot_preconditioner_dot_r / p_dot_A_dot_p;
-      residual_norm = std::sqrt(local_sums[3] + 2 * local_sums[2] +
-                                alpha * alpha * local_sums[1]);
-
-      r_dot_preconditioner_dot_r =
-        old_r_dot_preconditioner_dot_r -
-        alpha * (2. * local_sums[4] - alpha * local_sums[5]);
-      beta = r_dot_preconditioner_dot_r / old_r_dot_preconditioner_dot_r;
-    }
+      }
+    };
   } // namespace SolverCG
 } // namespace internal
 
@@ -699,22 +951,6 @@ SolverCG<VectorType>::solve(const MatrixType &        A,
 
   LogStream::Prefix prefix("cg");
 
-  // Memory allocation
-  typename VectorMemory<VectorType>::Pointer r_pointer(this->memory);
-  typename VectorMemory<VectorType>::Pointer p_pointer(this->memory);
-  typename VectorMemory<VectorType>::Pointer v_pointer(this->memory);
-  typename VectorMemory<VectorType>::Pointer z_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. The variable 'z'
-  // gets only used for the flexible variant of the CG method.
-  VectorType &r = *r_pointer;
-  VectorType &p = *p_pointer;
-  VectorType &v = *v_pointer;
-  VectorType &z = *z_pointer;
-
   // Should we build the matrix for eigenvalue computations?
   const bool do_eigenvalues =
     !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
@@ -726,33 +962,16 @@ SolverCG<VectorType>::solve(const MatrixType &        A,
 
   typename VectorType::value_type eigen_beta_alpha = 0;
 
-  // resize the vectors, but do not set the values since they'd be overwritten
-  // soon anyway.
-  r.reinit(x, true);
-  p.reinit(x, true);
-  v.reinit(x, true);
-  if (determine_beta_by_flexible_formula)
-    z.reinit(x, true);
-
-  int    it                         = 0;
-  number r_dot_preconditioner_dot_r = number();
-  number beta                       = number();
-  number alpha                      = number();
-  number old_beta                   = number();
-  number old_alpha                  = number();
-
-  // compute residual. if vector is zero, then short-circuit the full
-  // computation
-  if (!x.all_zero())
-    {
-      A.vmult(r, x);
-      r.sadd(-1., 1., b);
-    }
-  else
-    r.equ(1., b);
+  int it = 0;
+
+  internal::SolverCG::
+    IterationWorker<VectorType, MatrixType, PreconditionerType>
+      worker(
+        A, preconditioner, determine_beta_by_flexible_formula, this->memory, x);
+
+  worker.startup(b);
 
-  double residual_norm = r.l2_norm();
-  solver_state         = this->iteration_status(0, residual_norm, x);
+  solver_state = this->iteration_status(0, worker.residual_norm, x);
   if (solver_state != SolverControl::iterate)
     return;
 
@@ -760,37 +979,22 @@ SolverCG<VectorType>::solve(const MatrixType &        A,
     {
       it++;
 
-      internal::SolverCG::do_cg_iteration(A,
-                                          preconditioner,
-                                          additional_data,
-                                          it,
-                                          x,
-                                          r,
-                                          p,
-                                          v,
-                                          z,
-                                          r_dot_preconditioner_dot_r,
-                                          alpha,
-                                          beta,
-                                          residual_norm,
-                                          old_alpha,
-                                          old_beta);
-
-      old_alpha = alpha;
-      old_beta = beta;
-
-      print_vectors(it, x, r, p);
+      worker.do_iteration(it);
+
+      print_vectors(it, x, worker.r, worker.p);
 
       if (it > 1)
         {
-          this->coefficients_signal(old_alpha, beta);
-          // set up the vectors containing the diagonal and the off diagonal of
-          // the projected matrix.
+          this->coefficients_signal(worker.previous_alpha, worker.beta);
+          // set up the vectors containing the diagonal and the off diagonal
+          // of the projected matrix.
           if (do_eigenvalues)
             {
-              diagonal.push_back(number(1.) / old_alpha + eigen_beta_alpha);
-              eigen_beta_alpha = beta / old_alpha;
-              offdiagonal.push_back(std::sqrt(beta) / old_alpha);
+              diagonal.push_back(number(1.) / worker.previous_alpha +
+                                 eigen_beta_alpha);
+              eigen_beta_alpha = worker.beta / worker.previous_alpha;
+              offdiagonal.push_back(std::sqrt(worker.beta) /
+                                    worker.previous_alpha);
             }
           compute_eigs_and_cond(diagonal,
                                 offdiagonal,
@@ -798,16 +1002,18 @@ SolverCG<VectorType>::solve(const MatrixType &        A,
                                 all_condition_numbers_signal);
         }
 
-      solver_state = this->iteration_status(it, residual_norm, x);
+      solver_state = this->iteration_status(it, worker.residual_norm, x);
     }
 
+  worker.finalize_after_convergence(it);
+
   compute_eigs_and_cond(diagonal,
                         offdiagonal,
                         eigenvalues_signal,
                         condition_number_signal);
 
   AssertThrow(solver_state == SolverControl::success,
-              SolverControl::NoConvergence(it, residual_norm));
+              SolverControl::NoConvergence(it, worker.residual_norm));
 }
 
 
index c40adef8e3a85f775e9bbfe199db491461d72e68..cacd06ee370ce1c3e214acecebbcef55164dfd8b 100644 (file)
@@ -63,7 +63,7 @@ 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, residual_norm)
+DEAL:sor::Exception: SolverControl::NoConvergence(it, worker.residual_norm)
 DEAL:sor:Bicgstab::Starting value 3.000
 DEAL:sor:Bicgstab::Convergence step 5 value 4.201e-18
 DEAL:sor:GMRES::Starting value 1.462
@@ -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, residual_norm)
+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
@@ -91,7 +91,7 @@ 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, residual_norm)
+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
@@ -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, residual_norm)
+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
@@ -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, residual_norm)
+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
index e7c652f327667ca6e8881b4e2914ec625420e976..e5ae166fea8848c3f0b1544afd153d3c7d1e8cbc 100644 (file)
@@ -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, residual_norm)
+DEAL:no-fail::Exception: SolverControl::NoConvergence(it, worker.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

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.