Fix two bugs in interleaved CG method

diff --git a/include/deal.II/lac/solver_cg.h b/include/deal.II/lac/solver_cg.h
index 81f50070aa4f40ff77442de896049622b48b2973..667be6679ee023305e0dbc818cad0067db603e09 100644 (file)
--- a/include/deal.II/lac/solver_cg.h
+++ b/include/deal.II/lac/solver_cg.h
@@ -511,7 +511,6 @@ namespace internal
       Number beta;
       double residual_norm;
       Number previous_alpha;
-      Number previous_beta;
 
       IterationWorkerBase(const MatrixType &        A,
                           const PreconditionerType &preconditioner,
@@ -535,7 +534,6 @@ namespace internal
         , beta(Number())
         , residual_norm(0.0)
         , previous_alpha(Number())
-        , previous_beta(Number())
       {}
 
       void
@@ -604,8 +602,6 @@ namespace internal
 
         const Number previous_r_dot_preconditioner_dot_r =
           r_dot_preconditioner_dot_r;
-        this->previous_alpha = alpha;
-        this->previous_beta  = beta;
 
         if (std::is_same<PreconditionerType, PreconditionIdentity>::value ==
             false)
@@ -639,7 +635,9 @@ namespace internal
 
         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;
+
+        this->previous_alpha = alpha;
+        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)));
@@ -717,6 +715,9 @@ namespace internal
     {
       using Number = typename VectorType::value_type;
 
+      Number next_r_dot_preconditioner_dot_r;
+      Number previous_beta;
+
       IterationWorker(const MatrixType &        A,
                       const PreconditionerType &preconditioner,
                       const bool                flexible,
@@ -728,6 +729,8 @@ namespace internal
             flexible,
             memory,
             x)
+        , next_r_dot_preconditioner_dot_r(0.)
+        , previous_beta(0.)
       {}
 
       // This is the main iteration function, that will use some of the
@@ -735,6 +738,13 @@ namespace internal
       void
       do_iteration(const unsigned int iteration_index)
       {
+        if (iteration_index > 1)
+          {
+            previous_beta = this->beta;
+            this->beta    = next_r_dot_preconditioner_dot_r /
+                         this->r_dot_preconditioner_dot_r;
+          }
+
         std::array<VectorizedArray<Number>, 7> vectorized_sums = {};
 
         this->A.vmult(
@@ -759,24 +769,23 @@ namespace internal
                             this->r.get_mpi_communicator(),
                             dealii::ArrayView<Number>(scalar_sums.data(), 7));
 
-        this->previous_alpha = this->alpha;
-        this->previous_beta  = this->beta;
+        this->r_dot_preconditioner_dot_r = scalar_sums[6];
 
         const Number p_dot_A_dot_p = scalar_sums[0];
         Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero());
 
-        const Number previous_r_dot_preconditioner_dot_r = scalar_sums[6];
-        this->alpha = previous_r_dot_preconditioner_dot_r / p_dot_A_dot_p;
-        this->residual_norm = std::sqrt(
-          scalar_sums[3] +
-          this->alpha * (-2. * scalar_sums[2] + this->alpha * scalar_sums[1]));
+        this->previous_alpha = this->alpha;
+        this->alpha          = this->r_dot_preconditioner_dot_r / p_dot_A_dot_p;
 
-        this->r_dot_preconditioner_dot_r =
-          previous_r_dot_preconditioner_dot_r +
-          this->alpha * (-2. * scalar_sums[4] + this->alpha * scalar_sums[5]);
+        // Round-off errors near zero might yield negative values, so take
+        // the absolute value in the next two formulas
+        this->residual_norm = std::sqrt(std::abs(
+          scalar_sums[3] +
+          this->alpha * (-2. * scalar_sums[2] + this->alpha * scalar_sums[1])));
 
-        this->beta = this->r_dot_preconditioner_dot_r /
-                     previous_r_dot_preconditioner_dot_r;
+        next_r_dot_preconditioner_dot_r = std::abs(
+          this->r_dot_preconditioner_dot_r +
+          this->alpha * (-2. * scalar_sums[4] + this->alpha * scalar_sums[5]));
       }
 
       // Function that we use if the PreconditionerType implements an apply()
@@ -943,16 +952,21 @@ namespace internal
       typename std::enable_if<has_apply<PreconditionerType>, U>::type
       finalize_after_convergence(const unsigned int iteration_index)
       {
-        if (iteration_index % 2 == 1)
+        if (iteration_index % 2 == 1 || iteration_index == 2)
           this->x.add(this->alpha, this->p);
         else
           {
             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();
+            Number *x = this->x.begin();
+            Number *r = this->r.begin();
+            Number *p = this->p.begin();
+
+            // Note that we use 'beta' here rather than 'previous_beta' in the
+            // formula above, which is because the shift in beta ->
+            // previous_beta has not been applied at this stage, allowing us
+            // to recover the previous search direction
             const Number alpha_plus_previous_alpha_over_beta =
               this->alpha + this->previous_alpha / this->previous_beta;
             const Number previous_alpha_over_beta =
@@ -1147,7 +1161,7 @@ namespace internal
       typename std::enable_if<!has_apply<PreconditionerType>, U>::type
       finalize_after_convergence(const unsigned int iteration_index)
       {
-        if (iteration_index % 2 == 1)
+        if (iteration_index % 2 == 1 || iteration_index == 2)
           this->x.add(this->alpha, this->p);
         else
           {

diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h
index 402b52dc815f1a799f5bbf19a88b96d4c82f52cd..f19e2c9432a9ac0fe1a92d95a75b7f1cec9d9ad8 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.h
+++ b/include/deal.II/matrix_free/matrix_free.h
@@ -1031,7 +1031,7 @@ public:
    *
    * @param dof_handler_index_pre_post Since MatrixFree can be initialized
    * with a vector of DoFHandler objects, each of them will in general have
-   * vector sizes and thus different ranges returned to
+   * different vector sizes and thus different ranges returned to
    * `operation_before_loop` and `operation_after_loop`. Use this variable to
    * specify which one of the DoFHandler objects the index range should be
    * associated to. Defaults to the `dof_handler_index` 0.
@@ -4621,6 +4621,28 @@ namespace internal
 
 
 
+  // Apply a unit matrix operation to constrained DoFs: Default cases where we
+  // cannot detect a LinearAlgebra::distributed::Vector, we do not do
+  // anything, else we apply the constraints as a unit operation
+  template <typename VectorStruct1, typename VectorStruct2>
+  inline void
+  apply_operation_to_constrained_dofs(const std::vector<unsigned int> &,
+                                      const VectorStruct1 &,
+                                      VectorStruct2 &)
+  {}
+
+  template <typename Number>
+  inline void
+  apply_operation_to_constrained_dofs(
+    const std::vector<unsigned int> &                 constrained_dofs,
+    const LinearAlgebra::distributed::Vector<Number> &src,
+    LinearAlgebra::distributed::Vector<Number> &      dst)
+  {
+    for (const unsigned int i : constrained_dofs)
+      dst.local_element(i) = src.local_element(i);
+  }
+
+
   namespace MatrixFreeFunctions
   {
     // struct to select between a const interface and a non-const interface
@@ -4861,6 +4883,17 @@ namespace internal
     {
       if (operation_after_loop)
         {
+          // Run unit matrix operation on constrained dofs if we are at the
+          // last range
+          const std::vector<unsigned int> &partition_row_index =
+            matrix_free.get_task_info().partition_row_index;
+          if (range_index ==
+              partition_row_index[partition_row_index.size() - 2] - 1)
+            apply_operation_to_constrained_dofs(
+              matrix_free.get_constrained_dofs(dof_handler_index_pre_post),
+              src,
+              dst);
+
           const internal::MatrixFreeFunctions::DoFInfo &dof_info =
             matrix_free.get_dof_info(dof_handler_index_pre_post);
           if (range_index == numbers::invalid_unsigned_int)

diff --git a/tests/matrix_free/parallel_multigrid_interleave.cc b/tests/matrix_free/parallel_multigrid_interleave.cc
index 37fd4f7d75e553aade3507131fb80b8396f8355d..afebfe7e322c7ec0d2f829710779b1e7736fa2ed 100644 (file)
--- a/tests/matrix_free/parallel_multigrid_interleave.cc
+++ b/tests/matrix_free/parallel_multigrid_interleave.cc
@@ -153,8 +153,6 @@ public:
                    src,
                    operation_before_loop,
                    operation_after_loop);
-    for (unsigned int i : data.get_constrained_dofs())
-      dst.local_element(i) = src.local_element(i);
   }
 
   void

diff --git a/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output b/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output
index f33407f3f05f70ffa1353571435afeff1e4fb5e4..ba30099891b8d4155a9dd8ee13d3337fb21f5153 100644 (file)
--- a/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output
+++ b/tests/matrix_free/parallel_multigrid_interleave.with_p4est=true.with_lapack=true.mpirun=2.output
@@ -3,63 +3,63 @@ DEAL::Testing FE_Q<2>(1)
 DEAL::Number of degrees of freedom: 81
 DEAL:cg::Starting value 7.00000
 DEAL:cg::Convergence step 3 value 1.44113e-07
-DEAL::Number of calls to special vmult for Operator of size 4: 12
-DEAL::Number of calls to special vmult for Operator of size 9: 27
-DEAL::Number of calls to special vmult for Operator of size 25: 27
-DEAL::Number of calls to special vmult for Operator of size 81: 27
+DEAL::Number of calls to special vmult for Operator of size 4: 13
+DEAL::Number of calls to special vmult for Operator of size 9: 28
+DEAL::Number of calls to special vmult for Operator of size 25: 33
+DEAL::Number of calls to special vmult for Operator of size 81: 42
 DEAL::Testing FE_Q<2>(1)
 DEAL::Number of degrees of freedom: 289
 DEAL:cg::Starting value 15.0000
 DEAL:cg::Convergence step 3 value 1.49408e-06
-DEAL::Number of calls to special vmult for Operator of size 4: 12
-DEAL::Number of calls to special vmult for Operator of size 9: 27
-DEAL::Number of calls to special vmult for Operator of size 25: 27
-DEAL::Number of calls to special vmult for Operator of size 81: 27
-DEAL::Number of calls to special vmult for Operator of size 289: 27
+DEAL::Number of calls to special vmult for Operator of size 4: 13
+DEAL::Number of calls to special vmult for Operator of size 9: 28
+DEAL::Number of calls to special vmult for Operator of size 25: 33
+DEAL::Number of calls to special vmult for Operator of size 81: 42
+DEAL::Number of calls to special vmult for Operator of size 289: 42
 DEAL::Testing FE_Q<2>(3)
 DEAL::Number of degrees of freedom: 625
 DEAL:cg::Starting value 23.0000
-DEAL:cg::Convergence step 4 value 6.30996e-08
-DEAL::Number of calls to special vmult for Operator of size 16: 16
-DEAL::Number of calls to special vmult for Operator of size 49: 36
-DEAL::Number of calls to special vmult for Operator of size 169: 36
-DEAL::Number of calls to special vmult for Operator of size 625: 36
+DEAL:cg::Convergence step 4 value 6.60805e-08
+DEAL::Number of calls to special vmult for Operator of size 16: 19
+DEAL::Number of calls to special vmult for Operator of size 49: 48
+DEAL::Number of calls to special vmult for Operator of size 169: 51
+DEAL::Number of calls to special vmult for Operator of size 625: 51
 DEAL::Testing FE_Q<2>(3)
 DEAL::Number of degrees of freedom: 2401
 DEAL:cg::Starting value 47.0000
-DEAL:cg::Convergence step 4 value 1.99754e-07
-DEAL::Number of calls to special vmult for Operator of size 16: 16
-DEAL::Number of calls to special vmult for Operator of size 49: 36
-DEAL::Number of calls to special vmult for Operator of size 169: 36
-DEAL::Number of calls to special vmult for Operator of size 625: 36
-DEAL::Number of calls to special vmult for Operator of size 2401: 36
+DEAL:cg::Convergence step 4 value 2.29801e-07
+DEAL::Number of calls to special vmult for Operator of size 16: 19
+DEAL::Number of calls to special vmult for Operator of size 49: 48
+DEAL::Number of calls to special vmult for Operator of size 169: 51
+DEAL::Number of calls to special vmult for Operator of size 625: 51
+DEAL::Number of calls to special vmult for Operator of size 2401: 51
 DEAL::Testing FE_Q<3>(1)
 DEAL::Number of degrees of freedom: 125
 DEAL:cg::Starting value 5.19615
 DEAL:cg::Convergence step 3 value 5.16655e-09
-DEAL::Number of calls to special vmult for Operator of size 8: 12
-DEAL::Number of calls to special vmult for Operator of size 27: 27
-DEAL::Number of calls to special vmult for Operator of size 125: 27
+DEAL::Number of calls to special vmult for Operator of size 8: 13
+DEAL::Number of calls to special vmult for Operator of size 27: 28
+DEAL::Number of calls to special vmult for Operator of size 125: 34
 DEAL::Testing FE_Q<3>(1)
 DEAL::Number of degrees of freedom: 729
 DEAL:cg::Starting value 18.5203
-DEAL:cg::Convergence step 3 value 1.01773e-06
-DEAL::Number of calls to special vmult for Operator of size 8: 12
-DEAL::Number of calls to special vmult for Operator of size 27: 27
-DEAL::Number of calls to special vmult for Operator of size 125: 27
-DEAL::Number of calls to special vmult for Operator of size 729: 27
+DEAL:cg::Convergence step 3 value 1.00122e-06
+DEAL::Number of calls to special vmult for Operator of size 8: 13
+DEAL::Number of calls to special vmult for Operator of size 27: 28
+DEAL::Number of calls to special vmult for Operator of size 125: 34
+DEAL::Number of calls to special vmult for Operator of size 729: 42
 DEAL::Testing FE_Q<3>(2)
 DEAL::Number of degrees of freedom: 729
 DEAL:cg::Starting value 18.5203
-DEAL:cg::Convergence step 4 value 5.65213e-09
-DEAL::Number of calls to special vmult for Operator of size 27: 16
-DEAL::Number of calls to special vmult for Operator of size 125: 36
-DEAL::Number of calls to special vmult for Operator of size 729: 36
+DEAL:cg::Convergence step 4 value 5.65178e-09
+DEAL::Number of calls to special vmult for Operator of size 27: 17
+DEAL::Number of calls to special vmult for Operator of size 125: 43
+DEAL::Number of calls to special vmult for Operator of size 729: 51
 DEAL::Testing FE_Q<3>(2)
 DEAL::Number of degrees of freedom: 4913
 DEAL:cg::Starting value 58.0948
-DEAL:cg::Convergence step 4 value 3.37746e-08
-DEAL::Number of calls to special vmult for Operator of size 27: 16
-DEAL::Number of calls to special vmult for Operator of size 125: 36
-DEAL::Number of calls to special vmult for Operator of size 729: 36
-DEAL::Number of calls to special vmult for Operator of size 4913: 36
+DEAL:cg::Convergence step 4 value 6.26463e-08
+DEAL::Number of calls to special vmult for Operator of size 27: 17
+DEAL::Number of calls to special vmult for Operator of size 125: 43
+DEAL::Number of calls to special vmult for Operator of size 729: 51
+DEAL::Number of calls to special vmult for Operator of size 4913: 51

diff --git a/tests/matrix_free/parallel_multigrid_interleave_renumber.cc b/tests/matrix_free/parallel_multigrid_interleave_renumber.cc
index 430265140d3c8f7ab8a2818b4e897b3ac781ca50..0a26bbdfe676d3ae08074a5abc27d2f5f32590c5 100644 (file)
--- a/tests/matrix_free/parallel_multigrid_interleave_renumber.cc
+++ b/tests/matrix_free/parallel_multigrid_interleave_renumber.cc
@@ -194,8 +194,6 @@ public:
                    src,
                    operation_before_loop,
                    operation_after_loop);
-    for (unsigned int i : data.get_constrained_dofs())
-      dst.local_element(i) = src.local_element(i);
   }
 
   void

diff --git a/tests/matrix_free/parallel_multigrid_interleave_renumber.wih_p4est=true.with_lapack=true.mpirun=2.output b/tests/matrix_free/parallel_multigrid_interleave_renumber.with_p4est=true.with_lapack=true.mpirun=2.output
similarity index 83%
rename from tests/matrix_free/parallel_multigrid_interleave_renumber.wih_p4est=true.with_lapack=true.mpirun=2.output
rename to tests/matrix_free/parallel_multigrid_interleave_renumber.with_p4est=true.with_lapack=true.mpirun=2.output
index f33407f3f05f70ffa1353571435afeff1e4fb5e4..d51f7fbac3f6ba0e99999629e0c3b5c2fa84366d 100644 (file)
--- a/tests/matrix_free/parallel_multigrid_interleave_renumber.wih_p4est=true.with_lapack=true.mpirun=2.output
+++ b/tests/matrix_free/parallel_multigrid_interleave_renumber.with_p4est=true.with_lapack=true.mpirun=2.output
@@ -2,64 +2,64 @@
 DEAL::Testing FE_Q<2>(1)
 DEAL::Number of degrees of freedom: 81
 DEAL:cg::Starting value 7.00000
-DEAL:cg::Convergence step 3 value 1.44113e-07
-DEAL::Number of calls to special vmult for Operator of size 4: 12
-DEAL::Number of calls to special vmult for Operator of size 9: 27
-DEAL::Number of calls to special vmult for Operator of size 25: 27
-DEAL::Number of calls to special vmult for Operator of size 81: 27
+DEAL:cg::Convergence step 3 value 2.50859e-07
+DEAL::Number of calls to special vmult for Operator of size 4: 13
+DEAL::Number of calls to special vmult for Operator of size 9: 29
+DEAL::Number of calls to special vmult for Operator of size 25: 33
+DEAL::Number of calls to special vmult for Operator of size 81: 42
 DEAL::Testing FE_Q<2>(1)
 DEAL::Number of degrees of freedom: 289
 DEAL:cg::Starting value 15.0000
-DEAL:cg::Convergence step 3 value 1.49408e-06
-DEAL::Number of calls to special vmult for Operator of size 4: 12
-DEAL::Number of calls to special vmult for Operator of size 9: 27
-DEAL::Number of calls to special vmult for Operator of size 25: 27
-DEAL::Number of calls to special vmult for Operator of size 81: 27
-DEAL::Number of calls to special vmult for Operator of size 289: 27
+DEAL:cg::Convergence step 4 value 2.05438e-08
+DEAL::Number of calls to special vmult for Operator of size 4: 17
+DEAL::Number of calls to special vmult for Operator of size 9: 38
+DEAL::Number of calls to special vmult for Operator of size 25: 42
+DEAL::Number of calls to special vmult for Operator of size 81: 51
+DEAL::Number of calls to special vmult for Operator of size 289: 51
 DEAL::Testing FE_Q<2>(3)
 DEAL::Number of degrees of freedom: 625
 DEAL:cg::Starting value 23.0000
-DEAL:cg::Convergence step 4 value 6.30996e-08
-DEAL::Number of calls to special vmult for Operator of size 16: 16
-DEAL::Number of calls to special vmult for Operator of size 49: 36
-DEAL::Number of calls to special vmult for Operator of size 169: 36
-DEAL::Number of calls to special vmult for Operator of size 625: 36
+DEAL:cg::Convergence step 4 value 2.11999e-07
+DEAL::Number of calls to special vmult for Operator of size 16: 19
+DEAL::Number of calls to special vmult for Operator of size 49: 48
+DEAL::Number of calls to special vmult for Operator of size 169: 51
+DEAL::Number of calls to special vmult for Operator of size 625: 51
 DEAL::Testing FE_Q<2>(3)
 DEAL::Number of degrees of freedom: 2401
 DEAL:cg::Starting value 47.0000
-DEAL:cg::Convergence step 4 value 1.99754e-07
-DEAL::Number of calls to special vmult for Operator of size 16: 16
-DEAL::Number of calls to special vmult for Operator of size 49: 36
-DEAL::Number of calls to special vmult for Operator of size 169: 36
-DEAL::Number of calls to special vmult for Operator of size 625: 36
-DEAL::Number of calls to special vmult for Operator of size 2401: 36
+DEAL:cg::Convergence step 4 value 4.14038e-07
+DEAL::Number of calls to special vmult for Operator of size 16: 19
+DEAL::Number of calls to special vmult for Operator of size 49: 48
+DEAL::Number of calls to special vmult for Operator of size 169: 51
+DEAL::Number of calls to special vmult for Operator of size 625: 51
+DEAL::Number of calls to special vmult for Operator of size 2401: 51
 DEAL::Testing FE_Q<3>(1)
 DEAL::Number of degrees of freedom: 125
 DEAL:cg::Starting value 5.19615
-DEAL:cg::Convergence step 3 value 5.16655e-09
-DEAL::Number of calls to special vmult for Operator of size 8: 12
-DEAL::Number of calls to special vmult for Operator of size 27: 27
-DEAL::Number of calls to special vmult for Operator of size 125: 27
+DEAL:cg::Convergence step 3 value 4.30693e-08
+DEAL::Number of calls to special vmult for Operator of size 8: 13
+DEAL::Number of calls to special vmult for Operator of size 27: 28
+DEAL::Number of calls to special vmult for Operator of size 125: 34
 DEAL::Testing FE_Q<3>(1)
 DEAL::Number of degrees of freedom: 729
 DEAL:cg::Starting value 18.5203
-DEAL:cg::Convergence step 3 value 1.01773e-06
-DEAL::Number of calls to special vmult for Operator of size 8: 12
-DEAL::Number of calls to special vmult for Operator of size 27: 27
-DEAL::Number of calls to special vmult for Operator of size 125: 27
-DEAL::Number of calls to special vmult for Operator of size 729: 27
+DEAL:cg::Convergence step 3 value 3.62205e-07
+DEAL::Number of calls to special vmult for Operator of size 8: 13
+DEAL::Number of calls to special vmult for Operator of size 27: 28
+DEAL::Number of calls to special vmult for Operator of size 125: 34
+DEAL::Number of calls to special vmult for Operator of size 729: 42
 DEAL::Testing FE_Q<3>(2)
 DEAL::Number of degrees of freedom: 729
 DEAL:cg::Starting value 18.5203
-DEAL:cg::Convergence step 4 value 5.65213e-09
-DEAL::Number of calls to special vmult for Operator of size 27: 16
-DEAL::Number of calls to special vmult for Operator of size 125: 36
-DEAL::Number of calls to special vmult for Operator of size 729: 36
+DEAL:cg::Convergence step 4 value 1.13788e-08
+DEAL::Number of calls to special vmult for Operator of size 27: 17
+DEAL::Number of calls to special vmult for Operator of size 125: 43
+DEAL::Number of calls to special vmult for Operator of size 729: 51
 DEAL::Testing FE_Q<3>(2)
 DEAL::Number of degrees of freedom: 4913
 DEAL:cg::Starting value 58.0948
-DEAL:cg::Convergence step 4 value 3.37746e-08
-DEAL::Number of calls to special vmult for Operator of size 27: 16
-DEAL::Number of calls to special vmult for Operator of size 125: 36
-DEAL::Number of calls to special vmult for Operator of size 729: 36
-DEAL::Number of calls to special vmult for Operator of size 4913: 36
+DEAL:cg::Convergence step 4 value 9.17189e-08
+DEAL::Number of calls to special vmult for Operator of size 27: 17
+DEAL::Number of calls to special vmult for Operator of size 125: 43
+DEAL::Number of calls to special vmult for Operator of size 729: 51
+DEAL::Number of calls to special vmult for Operator of size 4913: 51