PreconditionRelaxation/Chebyshev fixes

diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h
index d580e7d4b645d0c5aa8bb7b582bd8533422e2d6c..bcf4f473a7e171c06c715c8edf2e02c072f965de 100644 (file)
--- a/include/deal.II/lac/precondition.h
+++ b/include/deal.II/lac/precondition.h
@@ -1145,7 +1145,7 @@ namespace internal
 
           Assert(transposed == false, ExcNotImplemented());
 
-          A.vmult(tmp, src);
+          A.vmult(tmp, dst);
 
           preconditioner.vmult(
             dst,
@@ -1237,7 +1237,7 @@ namespace internal
 
           A.vmult(
             tmp,
-            src,
+            dst,
             [&](const unsigned int start_range, const unsigned int end_range) {
               // zero 'tmp' before running the vmult
               // operation
@@ -1267,7 +1267,10 @@ namespace internal
                 }
             });
 
-          preconditioner.vmult(dst, tmp);
+          preconditioner.vmult(dst, tmp, [](const auto, const auto) {
+            // note: `dst` vector does not have to be zeroed out
+            // since we add the result into it
+          });
         }
     }
 
@@ -3162,7 +3165,7 @@ namespace internal
 
               DEAL_II_OPENMP_SIMD_PRAGMA
               for (std::size_t i = start_range; i < end_range; ++i)
-                tmp_ptr[i] = factor2 * (rhs_ptr[i] - tmp_ptr[i]);
+                tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
             });
 
           // 2) perform vector updates (including preconditioner application)

diff --git a/tests/lac/precondition_chebyshev_07.cc b/tests/lac/precondition_chebyshev_07.cc
index a7c4044ea6f187299d9ff0da125618a198c64c9c..2e073ae9a12d580363ae5920d0f8bfe9d7102a00 100644 (file)
--- a/tests/lac/precondition_chebyshev_07.cc
+++ b/tests/lac/precondition_chebyshev_07.cc
@@ -79,7 +79,10 @@ public:
     if (operation_before_matrix_vector_product)
       operation_before_matrix_vector_product(0, src.size());
 
-    diagonal_matrix.vmult(dst, src);
+    if (operation_before_matrix_vector_product)
+      diagonal_matrix.vmult_add(dst, src);
+    else
+      diagonal_matrix.vmult(dst, src);
 
     if (operation_after_matrix_vector_product)
       operation_after_matrix_vector_product(0, src.size());
@@ -211,7 +214,7 @@ main()
       std::vector<std::tuple<double, double>> results;
 
       {
-        // Test PreconditionChebyshev + DiagonalMatrix and matrix with
+        // 1) Test PreconditionChebyshev + DiagonalMatrix and matrix with
         // pre/post
         using PreconditionerType = DiagonalMatrix<VectorType>;
 
@@ -234,7 +237,7 @@ main()
       }
 
       {
-        /// Test PreconditionChebyshev + DiagonalMatrix and matrix without
+        // 2) Test PreconditionChebyshev + DiagonalMatrix and matrix without
         // pre/post
         using PreconditionerType = DiagonalMatrix<VectorType>;
 
@@ -253,7 +256,7 @@ main()
       }
 
       {
-        // Test PreconditionChebyshev + arbitrary preconditioner and matrix
+        // 3) Test PreconditionChebyshev + arbitrary preconditioner and matrix
         // without pre/post
         using PreconditionerType = MyDiagonalMatrix<VectorType>;
 
@@ -272,7 +275,7 @@ main()
       }
 
       {
-        // Test PreconditionChebyshev + preconditioner with pre/post and
+        // 4) Test PreconditionChebyshev + preconditioner with pre/post and
         // matrix without pre/post
         using PreconditionerType = MyDiagonalMatrixWithPreAndPost<VectorType>;
 
@@ -291,7 +294,7 @@ main()
       }
 
       {
-        // Test PreconditionChebyshev + preconditioner with pre/post and
+        // 5) Test PreconditionChebyshev + preconditioner with pre/post and
         // matrix with pre/post
         using PreconditionerType = MyDiagonalMatrixWithPreAndPost<VectorType>;
 
@@ -313,18 +316,15 @@ main()
         results.emplace_back(test(preconditioner, src));
       }
 
-      if (std::equal(results.begin(),
-                     results.end(),
-                     results.begin(),
-                     [](const auto &a, const auto &b) {
-                       if (std::abs(std::get<0>(a) - std::get<0>(b)) > 1e-6)
-                         return false;
+      if (std::all_of(results.begin(), results.end(), [&](const auto &a) {
+            if (std::abs(std::get<0>(a) - std::get<0>(results[0])) > 1e-6)
+              return false;
 
-                       if (std::abs(std::get<1>(a) - std::get<1>(b)) > 1e-6)
-                         return false;
+            if (std::abs(std::get<1>(a) - std::get<1>(results[0])) > 1e-6)
+              return false;
 
-                       return true;
-                     }))
+            return true;
+          }))
         deallog << "OK!" << std::endl;
       else
         deallog << "ERROR!" << std::endl;

diff --git a/tests/lac/precondition_relaxation_01.cc b/tests/lac/precondition_relaxation_01.cc
index 8c5a83c00a7d69fa7fec0e88261cd5f5d1c00a69..197e345b8dabb20bfd2aae5aa179ebbce9a340dc 100644 (file)
--- a/tests/lac/precondition_relaxation_01.cc
+++ b/tests/lac/precondition_relaxation_01.cc
@@ -71,15 +71,20 @@ public:
   vmult(VectorType &      dst,
         const VectorType &src,
         const std::function<void(const unsigned int, const unsigned int)>
-          &operation_before_matrix_vector_product,
+          &operation_before_matrix_vector_product = {},
         const std::function<void(const unsigned int, const unsigned int)>
-          &operation_after_matrix_vector_product) const
+          &operation_after_matrix_vector_product = {}) const
   {
-    operation_before_matrix_vector_product(0, src.size());
+    if (operation_before_matrix_vector_product)
+      operation_before_matrix_vector_product(0, src.size());
 
-    diagonal_matrix.vmult(dst, src);
+    if (operation_before_matrix_vector_product)
+      diagonal_matrix.vmult_add(dst, src);
+    else
+      diagonal_matrix.vmult(dst, src);
 
-    operation_after_matrix_vector_product(0, src.size());
+    if (operation_after_matrix_vector_product)
+      operation_after_matrix_vector_product(0, src.size());
   }
 
   VectorType &
@@ -215,7 +220,7 @@ main()
         std::vector<std::tuple<double, double, double, double>> results;
 
         {
-          // Test PreconditionJacobi
+          // 0) Test PreconditionJacobi
           PreconditionJacobi<MatrixType> preconditioner;
 
           PreconditionJacobi<MatrixType>::AdditionalData ad;
@@ -228,8 +233,8 @@ main()
         }
 
         {
-          // Test PreconditionRelaxation + DiagonalMatrix: optimized path with
-          // lambdas
+          // 1) Test PreconditionRelaxation + DiagonalMatrix and matrix with
+          // pre/post
           using PreconditionerType = DiagonalMatrix<VectorType>;
 
           using MyMatrixType = MySparseMatrix<MatrixType>;
@@ -252,7 +257,8 @@ main()
         }
 
         {
-          // Test PreconditionRelaxation + DiagonalMatrix: optimized path
+          // 2) Test PreconditionRelaxation + DiagonalMatrix and matrix without
+          // pre/post
           using PreconditionerType = DiagonalMatrix<VectorType>;
 
           PreconditionRelaxation<MatrixType, PreconditionerType> preconditioner;
@@ -270,8 +276,8 @@ main()
         }
 
         {
-          // Test PreconditionRelaxation + wrapper around DiagonalMatrix:
-          // optimized path cannot be used
+          // 3) Test PreconditionRelaxation + arbitrary preconditioner and
+          // matrix without pre/post
           using PreconditionerType = MyDiagonalMatrix<VectorType>;
 
           PreconditionRelaxation<MatrixType, PreconditionerType> preconditioner;
@@ -289,8 +295,8 @@ main()
         }
 
         {
-          // Test PreconditionRelaxation + wrapper around DiagonalMatrix with
-          // pre and post: alternative optimized path is taken
+          // 4) Test PreconditionRelaxation + preconditioner with pre/post and
+          // matrix without pre/post
           using PreconditionerType = MyDiagonalMatrixWithPreAndPost<VectorType>;
 
           PreconditionRelaxation<MatrixType, PreconditionerType> preconditioner;
@@ -307,24 +313,45 @@ main()
           results.emplace_back(test(preconditioner, src, false));
         }
 
-        if (std::equal(results.begin(),
-                       results.end(),
-                       results.begin(),
-                       [](const auto &a, const auto &b) {
-                         if (std::abs(std::get<0>(a) - std::get<0>(b)) > 1e-6)
-                           return false;
+        {
+          // 5) Test PreconditionRelaxation + preconditioner with pre/post and
+          // matrix with pre/post
+          using PreconditionerType = MyDiagonalMatrixWithPreAndPost<VectorType>;
+
+          using MyMatrixType = MySparseMatrix<MatrixType>;
+
+          MyMatrixType my_system_matrix(system_matrix);
+
+          PreconditionRelaxation<MyMatrixType, PreconditionerType>
+            preconditioner;
+
+          PreconditionRelaxation<MyMatrixType,
+                                 PreconditionerType>::AdditionalData ad;
+          ad.relaxation     = relaxation;
+          ad.n_iterations   = n_iterations;
+          ad.preconditioner = std::make_shared<PreconditionerType>();
+          ad.preconditioner->get_vector() = diagonal;
+
+          preconditioner.initialize(my_system_matrix, ad);
+
+          results.emplace_back(test(preconditioner, src, false));
+        }
+
+        if (std::all_of(results.begin(), results.end(), [&](const auto &a) {
+              if (std::abs(std::get<0>(a) - std::get<0>(results[0])) > 1e-6)
+                return false;
 
-                         if (std::abs(std::get<1>(a) - std::get<1>(b)) > 1e-6)
-                           return false;
+              if (std::abs(std::get<1>(a) - std::get<1>(results[0])) > 1e-6)
+                return false;
 
-                         if (std::abs(std::get<2>(a) - std::get<2>(b)) > 1e-6)
-                           return false;
+              if (std::abs(std::get<2>(a) - std::get<2>(results[0])) > 1e-6)
+                return false;
 
-                         if (std::abs(std::get<3>(a) - std::get<3>(b)) > 1e-6)
-                           return false;
+              if (std::abs(std::get<3>(a) - std::get<3>(results[0])) > 1e-6)
+                return false;
 
-                         return true;
-                       }))
+              return true;
+            }))
           deallog << "OK!" << std::endl;
         else
           deallog << "ERROR!" << std::endl;