From: Peter Munch <peterrmuench@gmail.com>
Date: Wed, 14 Sep 2022 09:56:16 +0000 (+0200)
Subject: PreconditionChebyshev: add specialization for A and P with ranges
X-Git-Tag: v9.5.0-rc1~963^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14271%2Fhead;p=dealii.git

PreconditionChebyshev: add specialization for A and P with ranges
---

diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h
index 5998ecc035..d580e7d4b6 100644
--- a/include/deal.II/lac/precondition.h
+++ b/include/deal.II/lac/precondition.h
@@ -3046,14 +3046,21 @@ namespace internal
 
     // We need to have a separate declaration for static const members
 
+    // general case and the case that the preconditioner can work on
+    // ranges (covered by vector_updates())
     template <
       typename MatrixType,
       typename VectorType,
       typename PreconditionerType,
-      std::enable_if_t<!has_vmult_with_std_functions<MatrixType,
-                                                     VectorType,
-                                                     PreconditionerType>,
-                       int> * = nullptr>
+      std::enable_if_t<
+        !has_vmult_with_std_functions<MatrixType,
+                                      VectorType,
+                                      PreconditionerType> &&
+          !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
+                                                          VectorType> &&
+            has_vmult_with_std_functions_for_precondition<MatrixType,
+                                                          VectorType>),
+        int> * = nullptr>
     inline void
     vmult_and_update(const MatrixType &        matrix,
                      const PreconditionerType &preconditioner,
@@ -3079,6 +3086,126 @@ namespace internal
                      solution);
     }
 
+    // case that both the operator and the preconditioner can work on
+    // subranges
+    template <
+      typename MatrixType,
+      typename VectorType,
+      typename PreconditionerType,
+      std::enable_if_t<
+        !has_vmult_with_std_functions<MatrixType,
+                                      VectorType,
+                                      PreconditionerType> &&
+          (has_vmult_with_std_functions_for_precondition<PreconditionerType,
+                                                         VectorType> &&
+           has_vmult_with_std_functions_for_precondition<MatrixType,
+                                                         VectorType>),
+        int> * = nullptr>
+    inline void
+    vmult_and_update(const MatrixType &        matrix,
+                     const PreconditionerType &preconditioner,
+                     const VectorType &        rhs,
+                     const unsigned int        iteration_index,
+                     const double              factor1_,
+                     const double              factor2_,
+                     VectorType &              solution,
+                     VectorType &              solution_old,
+                     VectorType &              temp_vector1,
+                     VectorType &              temp_vector2)
+    {
+      using Number = typename VectorType::value_type;
+
+      const Number factor1        = factor1_;
+      const Number factor1_plus_1 = 1. + factor1_;
+      const Number factor2        = factor2_;
+
+      if (iteration_index == 0)
+        {
+          preconditioner.vmult(
+            solution,
+            rhs,
+            [&](const unsigned int start_range, const unsigned int end_range) {
+              // zero 'solution' before running the vmult operation
+              if (end_range > start_range)
+                std::memset(solution.begin() + start_range,
+                            0,
+                            sizeof(Number) * (end_range - start_range));
+            },
+            [&](const unsigned int start_range, const unsigned int end_range) {
+              const auto solution_ptr = solution.begin();
+
+              DEAL_II_OPENMP_SIMD_PRAGMA
+              for (std::size_t i = start_range; i < end_range; ++i)
+                solution_ptr[i] *= factor2;
+            });
+        }
+      else
+        {
+          temp_vector1.reinit(rhs, true);
+          temp_vector2.reinit(rhs, true);
+
+          // 1) compute rediduum (including operator application)
+          matrix.vmult(
+            temp_vector1,
+            solution,
+            [&](const unsigned int start_range, const unsigned int end_range) {
+              // zero 'temp_vector1' before running the vmult
+              // operation
+              if (end_range > start_range)
+                std::memset(temp_vector1.begin() + start_range,
+                            0,
+                            sizeof(Number) * (end_range - start_range));
+            },
+            [&](const unsigned int start_range, const unsigned int end_range) {
+              const auto rhs_ptr = rhs.begin();
+              const auto tmp_ptr = temp_vector1.begin();
+
+              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]);
+            });
+
+          // 2) perform vector updates (including preconditioner application)
+          preconditioner.vmult(
+            temp_vector2,
+            temp_vector1,
+            [&](const unsigned int start_range, const unsigned int end_range) {
+              // zero 'temp_vector2' before running the vmult
+              // operation
+              if (end_range > start_range)
+                std::memset(temp_vector2.begin() + start_range,
+                            0,
+                            sizeof(Number) * (end_range - start_range));
+            },
+            [&](const unsigned int start_range, const unsigned int end_range) {
+              const auto solution_ptr     = solution.begin();
+              const auto solution_old_ptr = solution_old.begin();
+              const auto tmp_ptr          = temp_vector2.begin();
+
+              if (iteration_index == 1)
+                {
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (std::size_t i = start_range; i < end_range; ++i)
+                    tmp_ptr[i] =
+                      factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
+                }
+              else
+                {
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (std::size_t i = start_range; i < end_range; ++i)
+                    tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
+                                 factor1 * solution_old_ptr[i] +
+                                 factor2 * tmp_ptr[i];
+                }
+            });
+
+          solution.swap(temp_vector2);
+          solution_old.swap(temp_vector2);
+        }
+    }
+
+    // case that the operator can work on subranges and the preconditioner
+    // is a diagonal
     template <typename MatrixType,
               typename VectorType,
               typename PreconditionerType,
diff --git a/tests/lac/precondition_chebyshev_07.cc b/tests/lac/precondition_chebyshev_07.cc
new file mode 100644
index 0000000000..a7c4044ea6
--- /dev/null
+++ b/tests/lac/precondition_chebyshev_07.cc
@@ -0,0 +1,332 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test PreconditionChebyshev for different optimization levels. The
+// test is similar to precondition_relaxation_01.
+
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+
+#include "../tests.h"
+
+template <typename VectorType>
+class MyDiagonalMatrix
+{
+public:
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    diagonal_matrix.vmult(dst, src);
+  }
+
+  void
+  Tvmult(VectorType &dst, const VectorType &src) const
+  {
+    diagonal_matrix.vmult(dst, src);
+  }
+
+  VectorType &
+  get_vector()
+  {
+    return diagonal_matrix.get_vector();
+  }
+
+private:
+  DiagonalMatrix<VectorType> diagonal_matrix;
+};
+
+template <typename VectorType>
+class MyDiagonalMatrixWithPreAndPost
+{
+public:
+  void
+  vmult(VectorType &      dst,
+        const VectorType &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
+  {
+    if (operation_before_matrix_vector_product)
+      operation_before_matrix_vector_product(0, src.size());
+
+    diagonal_matrix.vmult(dst, src);
+
+    if (operation_after_matrix_vector_product)
+      operation_after_matrix_vector_product(0, src.size());
+  }
+
+  VectorType &
+  get_vector()
+  {
+    return diagonal_matrix.get_vector();
+  }
+
+private:
+  DiagonalMatrix<VectorType> diagonal_matrix;
+};
+
+template <typename SparseMatrixType>
+class MySparseMatrix : public Subscriptor
+{
+public:
+  MySparseMatrix(const SparseMatrixType &sparse_matrix)
+    : sparse_matrix(sparse_matrix)
+  {}
+
+  template <typename VectorType>
+  void
+  vmult(VectorType &      dst,
+        const VectorType &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
+  {
+    if (operation_before_matrix_vector_product)
+      operation_before_matrix_vector_product(0, src.size());
+
+    sparse_matrix.vmult(dst, src);
+
+    if (operation_after_matrix_vector_product)
+      operation_after_matrix_vector_product(0, src.size());
+  }
+
+  types::global_dof_index
+  m() const
+  {
+    return sparse_matrix.m();
+  }
+
+  double
+  el(unsigned int i, unsigned int j) const
+  {
+    return sparse_matrix.el(i, j);
+  }
+
+private:
+  const SparseMatrixType &sparse_matrix;
+};
+
+
+template <typename PreconditionerType, typename VectorType>
+std::tuple<double, double>
+test(const PreconditionerType &preconditioner, const VectorType &src)
+{
+  VectorType dst;
+  dst.reinit(src);
+
+  dst = 1.0;
+  preconditioner.vmult(dst, src);
+  const double norm_0 = dst.l2_norm();
+
+  dst = 1.0;
+  preconditioner.step(dst, src);
+  const double norm_1 = dst.l2_norm();
+
+  return std::tuple<double, double>{norm_0, norm_1};
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  using Number     = double;
+  using VectorType = Vector<Number>;
+  using MatrixType = SparseMatrix<Number>;
+
+  const unsigned int dim    = 2;
+  const unsigned int degree = 1;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(3);
+
+  FE_Q<dim>      fe(degree);
+  QGauss<dim>    quad(degree + 1);
+  MappingQ1<dim> mapping;
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  MatrixType system_matrix;
+
+  DynamicSparsityPattern dsp(dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern(dof_handler, dsp);
+
+  SparsityPattern sparsity_pattern;
+  sparsity_pattern.copy_from(dsp);
+  system_matrix.reinit(sparsity_pattern);
+
+  MatrixCreator::create_laplace_matrix(mapping,
+                                       dof_handler,
+                                       quad,
+                                       system_matrix);
+
+  VectorType diagonal(dof_handler.n_dofs());
+  VectorType src(dof_handler.n_dofs());
+
+  for (const auto &entry : system_matrix)
+    if (entry.row() == entry.column())
+      diagonal[entry.row()] = 1.0 / entry.value();
+
+  src = 1.0;
+
+  std::vector<unsigned int> n_iterationss{1, 2, 3};
+
+  for (const auto n_iterations : n_iterationss)
+    {
+      std::vector<std::tuple<double, double>> results;
+
+      {
+        // Test PreconditionChebyshev + DiagonalMatrix and matrix with
+        // pre/post
+        using PreconditionerType = DiagonalMatrix<VectorType>;
+
+        using MyMatrixType = MySparseMatrix<MatrixType>;
+
+        MyMatrixType my_system_matrix(system_matrix);
+
+        PreconditionChebyshev<MyMatrixType, VectorType, PreconditionerType>
+          preconditioner;
+
+        PreconditionChebyshev<MyMatrixType, VectorType, PreconditionerType>::
+          AdditionalData ad;
+        ad.degree         = 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));
+      }
+
+      {
+        /// Test PreconditionChebyshev + DiagonalMatrix and matrix without
+        // pre/post
+        using PreconditionerType = DiagonalMatrix<VectorType>;
+
+        PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>
+          preconditioner;
+
+        PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
+          AdditionalData ad;
+        ad.degree         = n_iterations;
+        ad.preconditioner = std::make_shared<PreconditionerType>();
+        ad.preconditioner->get_vector() = diagonal;
+
+        preconditioner.initialize(system_matrix, ad);
+
+        results.emplace_back(test(preconditioner, src));
+      }
+
+      {
+        // Test PreconditionChebyshev + arbitrary preconditioner and matrix
+        // without pre/post
+        using PreconditionerType = MyDiagonalMatrix<VectorType>;
+
+        PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>
+          preconditioner;
+
+        PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
+          AdditionalData ad;
+        ad.degree         = n_iterations;
+        ad.preconditioner = std::make_shared<PreconditionerType>();
+        ad.preconditioner->get_vector() = diagonal;
+
+        preconditioner.initialize(system_matrix, ad);
+
+        results.emplace_back(test(preconditioner, src));
+      }
+
+      {
+        // Test PreconditionChebyshev + preconditioner with pre/post and
+        // matrix without pre/post
+        using PreconditionerType = MyDiagonalMatrixWithPreAndPost<VectorType>;
+
+        PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>
+          preconditioner;
+
+        PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
+          AdditionalData ad;
+        ad.degree         = n_iterations;
+        ad.preconditioner = std::make_shared<PreconditionerType>();
+        ad.preconditioner->get_vector() = diagonal;
+
+        preconditioner.initialize(system_matrix, ad);
+
+        results.emplace_back(test(preconditioner, src));
+      }
+
+      {
+        // Test PreconditionChebyshev + preconditioner with pre/post and
+        // matrix with pre/post
+        using PreconditionerType = MyDiagonalMatrixWithPreAndPost<VectorType>;
+
+        using MyMatrixType = MySparseMatrix<MatrixType>;
+
+        MyMatrixType my_system_matrix(system_matrix);
+
+        PreconditionChebyshev<MyMatrixType, VectorType, PreconditionerType>
+          preconditioner;
+
+        PreconditionChebyshev<MyMatrixType, VectorType, PreconditionerType>::
+          AdditionalData ad;
+        ad.degree         = 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));
+      }
+
+      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::abs(std::get<1>(a) - std::get<1>(b)) > 1e-6)
+                         return false;
+
+                       return true;
+                     }))
+        deallog << "OK!" << std::endl;
+      else
+        deallog << "ERROR!" << std::endl;
+    }
+}
diff --git a/tests/lac/precondition_chebyshev_07.output b/tests/lac/precondition_chebyshev_07.output
new file mode 100644
index 0000000000..7ca52c43a9
--- /dev/null
+++ b/tests/lac/precondition_chebyshev_07.output
@@ -0,0 +1,4 @@
+
+DEAL::OK!
+DEAL::OK!
+DEAL::OK!