From: Peter Munch <peterrmuench@gmail.com>
Date: Wed, 7 Sep 2022 10:04:49 +0000 (+0200)
Subject: DiagonalMatrix: specialize for LA:d:V
X-Git-Tag: v9.5.0-rc1~986^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=76fc8de6f29385df6b60b97babf8ff0c86271022;p=dealii.git

DiagonalMatrix: specialize for LA:d:V
---

diff --git a/include/deal.II/lac/diagonal_matrix.h b/include/deal.II/lac/diagonal_matrix.h
index 0fd53c312f..4c8a0c75ac 100644
--- a/include/deal.II/lac/diagonal_matrix.h
+++ b/include/deal.II/lac/diagonal_matrix.h
@@ -24,6 +24,20 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+// forward declarations
+#ifndef DOXYGEN
+template <typename number>
+class Vector;
+namespace LinearAlgebra
+{
+  namespace distributed
+  {
+    template <typename, typename>
+    class Vector;
+  } // namespace distributed
+} // namespace LinearAlgebra
+#endif
+
 /**
  * This class represents a <i>n x n</i> diagonal matrix based on a vector of
  * size <i>n</i>. The matrix-vector products are realized by @p
@@ -386,12 +400,46 @@ DiagonalMatrix<VectorType>::add(const size_type  i,
 
 
 
+namespace internal
+{
+  namespace DiagonalMatrix
+  {
+    template <typename VectorType>
+    void
+    assign_and_scale(VectorType &      dst,
+                     const VectorType &src,
+                     const VectorType &diagonal)
+    {
+      dst = src;
+      dst.scale(diagonal);
+    }
+
+    template <typename Number>
+    void
+    assign_and_scale(
+      LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &      dst,
+      const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &src,
+      const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
+        &diagonal)
+    {
+      const auto dst_ptr      = dst.begin();
+      const auto src_ptr      = src.begin();
+      const auto diagonal_ptr = diagonal.begin();
+
+      DEAL_II_OPENMP_SIMD_PRAGMA
+      for (unsigned int i = 0; i < src.locally_owned_size(); ++i)
+        dst_ptr[i] = src_ptr[i] * diagonal_ptr[i];
+    }
+  } // namespace DiagonalMatrix
+} // namespace internal
+
+
+
 template <typename VectorType>
 void
 DiagonalMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
 {
-  dst = src;
-  dst.scale(diagonal);
+  internal::DiagonalMatrix::assign_and_scale(dst, src, diagonal);
 }
 
 
diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h
index 2972ddc408..b5702d9959 100644
--- a/include/deal.II/lac/precondition.h
+++ b/include/deal.II/lac/precondition.h
@@ -542,7 +542,8 @@ namespace internal
             typename PreconditionerType>
   constexpr bool has_vmult_with_std_functions =
     is_supported_operation<vmult_functions_t, MatrixType, VectorType> &&
-      std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value &&
+      std::is_same<PreconditionerType,
+                   dealii::DiagonalMatrix<VectorType>>::value &&
     (std::is_same<VectorType,
                   dealii::Vector<typename VectorType::value_type>>::value ||
      std::is_same<
@@ -1175,19 +1176,19 @@ namespace internal
     // 2) specialized implementation for inverse-diagonal preconditioner
     template <typename MatrixType,
               typename VectorType,
-              std::enable_if_t<
-                !IsBlockVector<VectorType>::value &&
-                  !has_vmult_with_std_functions<MatrixType,
-                                                VectorType,
-                                                DiagonalMatrix<VectorType>>,
-                VectorType> * = nullptr>
+              std::enable_if_t<!IsBlockVector<VectorType>::value &&
+                                 !has_vmult_with_std_functions<
+                                   MatrixType,
+                                   VectorType,
+                                   dealii::DiagonalMatrix<VectorType>>,
+                               VectorType> * = nullptr>
     void
-    step_operations(const MatrixType &                A,
-                    const DiagonalMatrix<VectorType> &preconditioner,
-                    VectorType &                      dst,
-                    const VectorType &                src,
-                    const double                      relaxation,
-                    VectorType &                      tmp,
+    step_operations(const MatrixType &                        A,
+                    const dealii::DiagonalMatrix<VectorType> &preconditioner,
+                    VectorType &                              dst,
+                    const VectorType &                        src,
+                    const double                              relaxation,
+                    VectorType &                              tmp,
                     VectorType &,
                     const unsigned int i,
                     const bool         transposed)
@@ -1247,19 +1248,19 @@ namespace internal
     // matrix that accepts ranges
     template <typename MatrixType,
               typename VectorType,
-              std::enable_if_t<
-                !IsBlockVector<VectorType>::value &&
-                  has_vmult_with_std_functions<MatrixType,
-                                               VectorType,
-                                               DiagonalMatrix<VectorType>>,
-                VectorType> * = nullptr>
+              std::enable_if_t<!IsBlockVector<VectorType>::value &&
+                                 has_vmult_with_std_functions<
+                                   MatrixType,
+                                   VectorType,
+                                   dealii::DiagonalMatrix<VectorType>>,
+                               VectorType> * = nullptr>
     void
-    step_operations(const MatrixType &                A,
-                    const DiagonalMatrix<VectorType> &preconditioner,
-                    VectorType &                      dst,
-                    const VectorType &                src,
-                    const double                      relaxation,
-                    VectorType &                      tmp,
+    step_operations(const MatrixType &                        A,
+                    const dealii::DiagonalMatrix<VectorType> &preconditioner,
+                    VectorType &                              dst,
+                    const VectorType &                        src,
+                    const double                              relaxation,
+                    VectorType &                              tmp,
                     VectorType &,
                     const unsigned int i,
                     const bool         transposed)
@@ -2871,15 +2872,16 @@ namespace internal
     // selection for diagonal matrix around deal.II vector
     template <typename Number>
     inline void
-    vector_updates(const ::dealii::Vector<Number> &                rhs,
-                   const DiagonalMatrix<::dealii::Vector<Number>> &jacobi,
-                   const unsigned int        iteration_index,
-                   const double              factor1,
-                   const double              factor2,
-                   ::dealii::Vector<Number> &solution_old,
-                   ::dealii::Vector<Number> &temp_vector1,
-                   ::dealii::Vector<Number> &,
-                   ::dealii::Vector<Number> &solution)
+    vector_updates(
+      const ::dealii::Vector<Number> &                        rhs,
+      const dealii::DiagonalMatrix<::dealii::Vector<Number>> &jacobi,
+      const unsigned int                                      iteration_index,
+      const double                                            factor1,
+      const double                                            factor2,
+      ::dealii::Vector<Number> &                              solution_old,
+      ::dealii::Vector<Number> &                              temp_vector1,
+      ::dealii::Vector<Number> &,
+      ::dealii::Vector<Number> &solution)
     {
       VectorUpdater<Number> upd(rhs.begin(),
                                 jacobi.get_vector().begin(),
@@ -2909,7 +2911,7 @@ namespace internal
     inline void
     vector_updates(
       const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
-      const DiagonalMatrix<
+      const dealii::DiagonalMatrix<
         LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>> &jacobi,
       const unsigned int iteration_index,
       const double       factor1,
@@ -3053,13 +3055,14 @@ namespace internal
     template <typename MatrixType, typename VectorType>
     inline void
     initialize_preconditioner(
-      const MatrixType &                           matrix,
-      std::shared_ptr<DiagonalMatrix<VectorType>> &preconditioner)
+      const MatrixType &                                   matrix,
+      std::shared_ptr<dealii::DiagonalMatrix<VectorType>> &preconditioner)
     {
       if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
         {
           if (preconditioner.get() == nullptr)
-            preconditioner = std::make_shared<DiagonalMatrix<VectorType>>();
+            preconditioner =
+              std::make_shared<dealii::DiagonalMatrix<VectorType>>();
 
           Assert(
             preconditioner->m() == 0,
@@ -3438,8 +3441,8 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
   // We do not need the second temporary vector in case we have a
   // DiagonalMatrix as preconditioner and use deal.II's own vectors
   using NumberType = typename VectorType::value_type;
-  if (std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value ==
-        false ||
+  if (std::is_same<PreconditionerType,
+                   dealii::DiagonalMatrix<VectorType>>::value == false ||
       (std::is_same<VectorType, dealii::Vector<NumberType>>::value == false &&
        ((std::is_same<VectorType,
                       LinearAlgebra::distributed::