Matrix-free evaluation: Avoid pointer aliasing

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 9c3bcfdb334aae3403439f47bcb89ba12750515c..4bbcf8de5d27b4db95a76f0ab5ce7fd934f6ba57 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -3403,6 +3403,16 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             ->component_dof_indices_offset[this->active_fe_index]
                                           [this->first_selected_component] *
           n_lanes;
+
+      std::array<typename VectorType::value_type *, n_components> src_ptrs;
+      if (n_components == 1 || this->n_fe_components == 1)
+        for (unsigned int comp = 0; comp < n_components; ++comp)
+          src_ptrs[comp] =
+            const_cast<typename VectorType::value_type *>(src[comp]->begin());
+      else
+        src_ptrs[0] =
+          const_cast<typename VectorType::value_type *>(src[0]->begin());
+
       if (n_components == 1 || this->n_fe_components == 1)
         for (unsigned int i = 0; i < dofs_per_component;
              ++i, dof_indices += n_lanes)
@@ -3410,14 +3420,19 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             operation.process_dof_gather(dof_indices,
                                          *src[comp],
                                          0,
+                                         src_ptrs[comp],
                                          values_dofs[comp][i],
                                          vector_selector);
       else
         for (unsigned int comp = 0; comp < n_components; ++comp)
           for (unsigned int i = 0; i < dofs_per_component;
                ++i, dof_indices += n_lanes)
-            operation.process_dof_gather(
-              dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
+            operation.process_dof_gather(dof_indices,
+                                         *src[0],
+                                         0,
+                                         src_ptrs[0],
+                                         values_dofs[comp][i],
+                                         vector_selector);
       return;
     }
 
@@ -3962,6 +3977,15 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
                  interleaved_contiguous_strided)
         {
+          std::array<typename VectorType::value_type *, n_components> src_ptrs;
+          if (n_components == 1 || this->n_fe_components == 1)
+            for (unsigned int comp = 0; comp < n_components; ++comp)
+              src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
+                src[comp]->begin());
+          else
+            src_ptrs[0] =
+              const_cast<typename VectorType::value_type *>(src[0]->begin());
+
           if (n_components == 1 || this->n_fe_components == 1)
             for (unsigned int i = 0; i < dofs_per_component; ++i)
               {
@@ -3969,6 +3993,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                   operation.process_dof_gather(dof_indices.data(),
                                                *src[comp],
                                                i * n_lanes,
+                                               src_ptrs[comp] + i * n_lanes,
                                                values_dofs[comp][i],
                                                vector_selector);
               }
@@ -3976,12 +4001,13 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             for (unsigned int comp = 0; comp < n_components; ++comp)
               for (unsigned int i = 0; i < dofs_per_component; ++i)
                 {
-                  operation.process_dof_gather(dof_indices.data(),
-                                               *src[0],
-                                               (comp * dofs_per_component + i) *
-                                                 n_lanes,
-                                               values_dofs[comp][i],
-                                               vector_selector);
+                  operation.process_dof_gather(
+                    dof_indices.data(),
+                    *src[0],
+                    (comp * dofs_per_component + i) * n_lanes,
+                    src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
+                    values_dofs[comp][i],
+                    vector_selector);
                 }
         }
       else
@@ -3990,6 +4016,15 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                    internal::MatrixFreeFunctions::DoFInfo::
                      IndexStorageVariants::interleaved_contiguous_mixed_strides,
                  ExcNotImplemented());
+          std::array<typename VectorType::value_type *, n_components> src_ptrs;
+          if (n_components == 1 || this->n_fe_components == 1)
+            for (unsigned int comp = 0; comp < n_components; ++comp)
+              src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
+                src[comp]->begin());
+          else
+            src_ptrs[0] =
+              const_cast<typename VectorType::value_type *>(src[0]->begin());
+
           const unsigned int *offsets =
             &dof_info.dof_indices_interleave_strides[ind][n_lanes * this->cell];
           if (n_components == 1 || this->n_fe_components == 1)
@@ -3999,6 +4034,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                   operation.process_dof_gather(dof_indices.data(),
                                                *src[comp],
                                                0,
+                                               src_ptrs[comp],
                                                values_dofs[comp][i],
                                                vector_selector);
                 DEAL_II_OPENMP_SIMD_PRAGMA
@@ -4012,6 +4048,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                   operation.process_dof_gather(dof_indices.data(),
                                                *src[0],
                                                0,
+                                               src_ptrs[0],
                                                values_dofs[comp][i],
                                                vector_selector);
                   DEAL_II_OPENMP_SIMD_PRAGMA

diff --git a/include/deal.II/matrix_free/vector_access_internal.h b/include/deal.II/matrix_free/vector_access_internal.h
index 9b5535d209ad700b81f3eb836574659536a0c913..1504b5b3d8f89fc003a4a85d1d167156898603b2 100644 (file)
--- a/include/deal.II/matrix_free/vector_access_internal.h
+++ b/include/deal.II/matrix_free/vector_access_internal.h
@@ -439,22 +439,28 @@ namespace internal
     // gather
     template <typename VectorType>
     void
-    process_dof_gather(const unsigned int  *indices,
-                       VectorType          &vec,
-                       const unsigned int   constant_offset,
-                       VectorizedArrayType &res,
+    process_dof_gather(const unsigned int              *indices,
+                       VectorType                      &vec,
+                       const unsigned int               constant_offset,
+                       typename VectorType::value_type *vec_ptr,
+                       VectorizedArrayType             &res,
                        std::integral_constant<bool, true>) const
     {
+      (void)constant_offset;
+      (void)vec;
+
 #ifdef DEBUG
       // in debug mode, run non-vectorized version because this path
       // has additional checks (e.g., regarding ghosting)
+      Assert(vec_ptr == vec.begin() + constant_offset, ExcInternalError());
       process_dof_gather(indices,
                          vec,
                          constant_offset,
+                         vec_ptr,
                          res,
                          std::integral_constant<bool, false>());
 #else
-      res.gather(vec.begin() + constant_offset, indices);
+      res.gather(vec_ptr, indices);
 #endif
     }
 
@@ -464,9 +470,10 @@ namespace internal
     // manually load the data
     template <typename VectorType>
     void
-    process_dof_gather(const unsigned int  *indices,
-                       const VectorType    &vec,
-                       const unsigned int   constant_offset,
+    process_dof_gather(const unsigned int *indices,
+                       const VectorType   &vec,
+                       const unsigned int  constant_offset,
+                       typename VectorType::value_type *,
                        VectorizedArrayType &res,
                        std::integral_constant<bool, false>) const
     {
@@ -683,21 +690,26 @@ namespace internal
     // scatter
     template <typename VectorType>
     void
-    process_dof_gather(const unsigned int  *indices,
-                       VectorType          &vec,
-                       const unsigned int   constant_offset,
-                       VectorizedArrayType &res,
+    process_dof_gather(const unsigned int              *indices,
+                       VectorType                      &vec,
+                       const unsigned int               constant_offset,
+                       typename VectorType::value_type *vec_ptr,
+                       VectorizedArrayType             &res,
                        std::integral_constant<bool, true>) const
     {
+      (void)constant_offset;
+      (void)vec_ptr;
+      (void)vec;
+
 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512
       for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
         vector_access(vec, indices[v] + constant_offset) += res[v];
 #else
       // only use gather in case there is also scatter.
       VectorizedArrayType tmp;
-      tmp.gather(vec.begin() + constant_offset, indices);
+      tmp.gather(vec_ptr, indices);
       tmp += res;
-      tmp.scatter(indices, vec.begin() + constant_offset);
+      tmp.scatter(indices, vec_ptr);
 #endif
     }
 
@@ -707,9 +719,10 @@ namespace internal
     // manually append all data
     template <typename VectorType>
     void
-    process_dof_gather(const unsigned int  *indices,
-                       VectorType          &vec,
-                       const unsigned int   constant_offset,
+    process_dof_gather(const unsigned int *indices,
+                       VectorType         &vec,
+                       const unsigned int  constant_offset,
+                       typename VectorType::value_type *,
                        VectorizedArrayType &res,
                        std::integral_constant<bool, false>) const
     {
@@ -913,22 +926,28 @@ namespace internal
 
     template <typename VectorType>
     void
-    process_dof_gather(const unsigned int  *indices,
-                       VectorType          &vec,
-                       const unsigned int   constant_offset,
-                       VectorizedArrayType &res,
+    process_dof_gather(const unsigned int              *indices,
+                       VectorType                      &vec,
+                       const unsigned int               constant_offset,
+                       typename VectorType::value_type *vec_ptr,
+                       VectorizedArrayType             &res,
                        std::integral_constant<bool, true>) const
     {
-      res.scatter(indices, vec.begin() + constant_offset);
+      (void)vec_ptr;
+      (void)constant_offset;
+      (void)vec;
+      Assert(vec_ptr == vec.begin() + constant_offset, ExcInternalError());
+      res.scatter(indices, vec_ptr);
     }
 
 
 
     template <typename VectorType>
     void
-    process_dof_gather(const unsigned int  *indices,
-                       VectorType          &vec,
-                       const unsigned int   constant_offset,
+    process_dof_gather(const unsigned int *indices,
+                       VectorType         &vec,
+                       const unsigned int  constant_offset,
+                       typename VectorType::value_type *,
                        VectorizedArrayType &res,
                        std::integral_constant<bool, false>) const
     {