Mohamad/author fix

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);
 }
 
 

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<
@@ -1086,13 +1087,16 @@ namespace internal
 
     // 1) specialized implementation with a preconditioner that accepts
     // ranges
-    template <typename MatrixType,
-              typename PreconditionerType,
-              typename VectorType,
-              std::enable_if_t<has_vmult_with_std_functions_for_precondition<
-                                 PreconditionerType,
-                                 VectorType>,
-                               int> * = nullptr>
+    template <
+      typename MatrixType,
+      typename PreconditionerType,
+      typename VectorType,
+      std::enable_if_t<
+        has_vmult_with_std_functions_for_precondition<PreconditionerType,
+                                                      VectorType> &&
+          !has_vmult_with_std_functions_for_precondition<MatrixType,
+                                                         VectorType>,
+        int> * = nullptr>
     void
     step_operations(const MatrixType &        A,
                     const PreconditionerType &preconditioner,
@@ -1172,22 +1176,116 @@ namespace internal
         }
     }
 
-    // 2) specialized implementation for inverse-diagonal preconditioner
+    // 2) specialized implementation with a preconditioner and a matrix that
+    // accepts ranges
+    template <
+      typename MatrixType,
+      typename PreconditionerType,
+      typename VectorType,
+      std::enable_if_t<
+        has_vmult_with_std_functions_for_precondition<PreconditionerType,
+                                                      VectorType> &&
+          has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
+        int> * = nullptr>
+    void
+    step_operations(const MatrixType &        A,
+                    const PreconditionerType &preconditioner,
+                    VectorType &              dst,
+                    const VectorType &        src,
+                    const double              relaxation,
+                    VectorType &              tmp,
+                    VectorType &,
+                    const unsigned int i,
+                    const bool         transposed)
+    {
+      (void)transposed;
+      using Number = typename VectorType::value_type;
+
+      if (i == 0)
+        {
+          Number *      dst_ptr = dst.begin();
+          const Number *src_ptr = src.begin();
+
+          preconditioner.vmult(
+            dst,
+            src,
+            [&](const unsigned int start_range, const unsigned int end_range) {
+              // zero 'dst' before running the vmult operation
+              if (end_range > start_range)
+                std::memset(dst.begin() + start_range,
+                            0,
+                            sizeof(Number) * (end_range - start_range));
+            },
+            [&](const unsigned int start_range, const unsigned int end_range) {
+              if (relaxation == 1.0)
+                return; // nothing to do
+
+              const auto src_ptr = src.begin();
+              const auto dst_ptr = dst.begin();
+
+              DEAL_II_OPENMP_SIMD_PRAGMA
+              for (std::size_t i = start_range; i < end_range; ++i)
+                dst_ptr[i] *= relaxation;
+            });
+        }
+      else
+        {
+          tmp.reinit(src, true);
+
+          Assert(transposed == false, ExcNotImplemented());
+
+          A.vmult(
+            tmp,
+            src,
+            [&](const unsigned int start_range, const unsigned int end_range) {
+              // zero 'tmp' before running the vmult
+              // operation
+              if (end_range > start_range)
+                std::memset(tmp.begin() + start_range,
+                            0,
+                            sizeof(Number) * (end_range - start_range));
+            },
+            [&](const unsigned int start_range, const unsigned int end_range) {
+              const auto src_ptr = src.begin();
+              const auto tmp_ptr = tmp.begin();
+
+              if (relaxation == 1.0)
+                {
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (std::size_t i = start_range; i < end_range; ++i)
+                    tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
+                }
+              else
+                {
+                  // note: we scale the residual here to be able to add into
+                  // the dst vector, which contains the solution from the last
+                  // iteration
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (std::size_t i = start_range; i < end_range; ++i)
+                    tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
+                }
+            });
+
+          preconditioner.vmult(dst, tmp);
+        }
+    }
+
+    // 3) 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)
@@ -1243,23 +1341,23 @@ namespace internal
         }
     }
 
-    // 3) specialized implementation for inverse-diagonal preconditioner and
+    // 4) specialized implementation for inverse-diagonal preconditioner and
     // 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 +2969,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 +3008,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 +3152,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 +3538,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::