Use vectorized gather operation if available.

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index ea77ffa1486c2528da196602afc51664b1fa7c58..a8e8bef3e863d490a5f9872d80390e6d127c46bb 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -2765,6 +2765,25 @@ namespace internal
       res = vector_access (const_cast<const VectorType &>(vec), index);
     }
 
+    template <typename VectorType>
+    void process_dof_gather (const unsigned int      *indices,
+                             VectorType              &vec,
+                             VectorizedArray<Number> &res,
+                             internal::bool2type<true>) const
+    {
+      res.gather(vec.begin(), indices);
+    }
+
+    template <typename VectorType>
+    void process_dof_gather (const unsigned int      *indices,
+                             VectorType              &vec,
+                             VectorizedArray<Number> &res,
+                             internal::bool2type<false>) const
+    {
+      for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+        res[v] = vector_access(const_cast<const VectorType &>(vec), indices[v]);
+    }
+
     template <typename VectorType>
     void process_dof_global (const types::global_dof_index index,
                              VectorType         &vec,
@@ -2812,6 +2831,34 @@ namespace internal
       vector_access (vec, index) += res;
     }
 
+    template <typename VectorType>
+    void process_dof_gather (const unsigned int      *indices,
+                             VectorType              &vec,
+                             VectorizedArray<Number> &res,
+                             internal::bool2type<true>) const
+    {
+#if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3
+      for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+        vector_access(vec, indices[v]) += res[v];
+#else
+      // only use gather in case there is also scatter
+      VectorizedArray<Number> tmp;
+      tmp.gather(vec.begin(), indices);
+      tmp += res;
+      tmp.scatter(indices, vec.begin());
+#endif
+    }
+
+    template <typename VectorType>
+    void process_dof_gather (const unsigned int      *indices,
+                             VectorType              &vec,
+                             VectorizedArray<Number> &res,
+                             internal::bool2type<false>) const
+    {
+      for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+        vector_access(vec, indices[v]) += res[v];
+    }
+
     template <typename VectorType>
     void process_dof_global (const types::global_dof_index index,
                              VectorType         &vec,
@@ -2858,6 +2905,25 @@ namespace internal
       vector_access (vec, index) = res;
     }
 
+    template <typename VectorType>
+    void process_dof_gather (const unsigned int      *indices,
+                             VectorType              &vec,
+                             VectorizedArray<Number> &res,
+                             internal::bool2type<true>) const
+    {
+      res.scatter(indices, vec.begin());
+    }
+
+    template <typename VectorType>
+    void process_dof_gather (const unsigned int      *indices,
+                             VectorType              &vec,
+                             VectorizedArray<Number> &res,
+                             internal::bool2type<false>) const
+    {
+      for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+        vector_access(vec, indices[v]) = res[v];
+    }
+
     template <typename VectorType>
     void process_dof_global (const types::global_dof_index index,
                              VectorType         &vec,
@@ -3071,11 +3137,11 @@ FEEvaluationBase<dim,n_components_,Number>
               // vectorization loop
               AssertDimension (dof_info->end_indices(cell)-dof_indices,
                                static_cast<int>(n_local_dofs));
-              for (unsigned int j=0; j<n_local_dofs; j+=VectorizedArray<Number>::n_array_elements)
-                for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
-                  for (unsigned int comp=0; comp<n_components; ++comp)
-                    operation.process_dof (dof_indices[j+v], *src[comp],
-                                           local_data[comp][j+v]);
+              for (unsigned int j=0, ind=0; j<dofs_per_cell; ++j, ind += VectorizedArray<Number>::n_array_elements)
+                for (unsigned int comp=0; comp<n_components; ++comp)
+                  operation.process_dof_gather(dof_indices+ind,
+                                               *src[comp], values_dofs[comp][j],
+                                               dealii::internal::bool2type<types_are_equal<typename VectorType::value_type,Number>::value>());
             }
         }
 
@@ -3214,10 +3280,11 @@ FEEvaluationBase<dim,n_components_,Number>
               // vectorization loop
               AssertDimension (dof_info->end_indices(cell)-dof_indices,
                                static_cast<int>(n_local_dofs));
-              for (unsigned int j=0; j<n_local_dofs; j+=VectorizedArray<Number>::n_array_elements)
-                for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
-                  operation.process_dof (dof_indices[j+v], *src[0],
-                                         local_data[j+v]);
+              for (unsigned int comp=0, ind=0; comp<n_components; ++comp)
+                for (unsigned int j=0; j<dofs_per_cell; ++j, ind += VectorizedArray<Number>::n_array_elements)
+                  operation.process_dof_gather(dof_indices+ind,
+                                               *src[0], values_dofs[comp][j],
+                                               dealii::internal::bool2type<types_are_equal<typename VectorType::value_type,Number>::value>());
             }
         }