Alternative version with manual vectorization

diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc
index 88b2bd776fad9152a3d41cb2e9cbf337bddf3153..d8cacae68aa256edfa4299b8eb19a289bc18008d 100644 (file)
--- a/source/grid/manifold.cc
+++ b/source/grid/manifold.cc
@@ -15,6 +15,7 @@
 
 #include <deal.II/base/table.h>
 #include <deal.II/base/tensor.h>
+#include <deal.II/base/vectorization.h>
 
 #include <deal.II/fe/fe_q.h>
 
@@ -577,7 +578,7 @@ FlatManifold<dim, spacedim>::get_new_point(
 {
   Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0) - 1.0) <
            1e-10,
-         ExcMessage("The weights for the individual points should sum to 1!"));
+         ExcMessage("The weights for the new point should sum to 1!"));
 
   Point<spacedim> p;
 
@@ -630,35 +631,6 @@ FlatManifold<dim, spacedim>::get_new_point(
 
 
 
-namespace
-{
-  template <int n_eval_points, int spacedim>
-  std::array<Point<spacedim>, n_eval_points>
-  interpolate_points(const unsigned int         n_points,
-                     const Point<spacedim> *    surrounding_points,
-                     const double *             weights,
-                     const Tensor<1, spacedim> &periodicity)
-  {
-    std::array<Point<spacedim>, n_eval_points> new_point;
-    for (unsigned int p = 0; p < n_points; ++p)
-      {
-        const Point<spacedim> surrounding_point = surrounding_points[p];
-        for (unsigned int i = 0; i < n_eval_points; ++i)
-          new_point[i] += surrounding_point * weights[p + i * n_points];
-      }
-
-    for (unsigned int d = 0; d < spacedim; ++d)
-      if (periodicity[d] > 0)
-        for (unsigned int i = 0; i < n_eval_points; ++i)
-          if (new_point[i][d] < 0)
-            new_point[i][d] += periodicity[d];
-
-    return new_point;
-  }
-} // namespace
-
-
-
 template <int dim, int spacedim>
 void
 FlatManifold<dim, spacedim>::get_new_points(
@@ -669,81 +641,70 @@ FlatManifold<dim, spacedim>::get_new_points(
   AssertDimension(surrounding_points.size(), weights.size(1));
   if (weights.size(0) == 0)
     return;
+  AssertDimension(new_points.size(), weights.size(0));
 
   const std::size_t n_points = surrounding_points.size();
 
-  Tensor<1, spacedim> minP = periodicity;
-  for (unsigned int d = 0; d < spacedim; ++d)
-    if (periodicity[d] > 0)
-      for (unsigned int i = 0; i < n_points; ++i)
-        {
-          minP[d] = std::min(minP[d], surrounding_points[i][d]);
-          Assert((surrounding_points[i][d] <
-                  periodicity[d] + tolerance * periodicity[d]) ||
-                   (surrounding_points[i][d] >= -tolerance * periodicity[d]),
-                 ExcPeriodicBox(d, surrounding_points[i], periodicity[i]));
-        }
-
-  // check whether periodicity shifts some of the points. Only do this if
-  // necessary to avoid memory allocation
-  const Point<spacedim> *surrounding_points_start = surrounding_points.data();
-
-  boost::container::small_vector<Point<spacedim>, 200> modified_points;
-  bool adjust_periodicity = false;
-  for (unsigned int d = 0; d < spacedim; ++d)
-    if (periodicity[d] > 0)
-      for (unsigned int i = 0; i < n_points; ++i)
-        if ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0)
-          {
-            adjust_periodicity = true;
-            break;
-          }
-  if (adjust_periodicity == true)
-    {
-      modified_points.insert(modified_points.end(),
-                             surrounding_points.begin(),
-                             surrounding_points.end());
-      for (unsigned int d = 0; d < spacedim; ++d)
-        if (periodicity[d] > 0)
-          for (unsigned int i = 0; i < n_points; ++i)
-            if ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0)
-              modified_points[i][d] -= periodicity[d];
-      surrounding_points_start = modified_points.data();
-    }
-
-  for (unsigned int r = 0; r < weights.size(0); ++r)
-    Assert(std::abs(
-             std::accumulate(&weights(r, 0), &weights(r, 0) + n_points, 0.0) -
-             1.0) < 1e-10,
-           ExcMessage("The weights each point should sum to 1!"));
-
-  // Now perform the interpolation
-  for (unsigned int row = 0; row < weights.size(0);)
+  // if there is no periodicity, use an optimized implementation with
+  // VectorizedArray, otherwise go to the get_new_point function for adjusting
+  // the domain
+  if (periodicity == Tensor<1, spacedim>())
     {
-      // For performance reasons, try to compute three points at once
-      if (row + 3 <= weights.size(0))
-        {
-          const std::array<Point<spacedim>, 3> points = interpolate_points<3>(
-            n_points, surrounding_points_start, &weights(row, 0), periodicity);
-          for (unsigned int i = 0; i < 3; ++i)
-            new_points[row + i] =
-              project_to_manifold(make_array_view(surrounding_points.begin(),
-                                                  surrounding_points.end()),
-                                  points[i]);
-          row += 3;
-        }
-      else
+      for (unsigned int row = 0; row < weights.size(0); ++row)
+        Assert(std::abs(std::accumulate(&weights(row, 0),
+                                        &weights(row, 0) + n_points,
+                                        0.0) -
+                        1.0) < 1e-10,
+               ExcMessage("The weights for each of the points should sum to "
+                          "1!"));
+
+      constexpr std::size_t n_lanes =
+        std::min<std::size_t>(VectorizedArray<double>::size(), 4);
+      using VectorizedArrayType        = VectorizedArray<double, n_lanes>;
+      const std::size_t n_regular_cols = (n_points / n_lanes) * n_lanes;
+      for (unsigned int row = 0; row < weights.size(0); row += n_lanes)
         {
-          new_points[row] =
-            project_to_manifold(make_array_view(surrounding_points.begin(),
-                                                surrounding_points.end()),
-                                interpolate_points<1>(n_points,
-                                                      surrounding_points_start,
-                                                      &weights(row, 0),
-                                                      periodicity)[0]);
-          ++row;
+          std::array<unsigned int, n_lanes> offsets;
+          // ensure to not access out of bounds, possibly duplicating some
+          // entries
+          for (std::size_t i = 0; i < n_lanes; ++i)
+            offsets[i] =
+              std::min<unsigned int>((row + i) * n_points,
+                                     (weights.size(0) - 1) * n_points);
+          Point<spacedim, VectorizedArrayType> point;
+          for (std::size_t col = 0; col < n_regular_cols; col += n_lanes)
+            {
+              std::array<VectorizedArrayType, n_lanes> vectorized_weights;
+              vectorized_load_and_transpose(n_lanes,
+                                            &weights(0, 0) + col,
+                                            offsets.data(),
+                                            vectorized_weights.data());
+              for (std::size_t i = 0; i < n_lanes; ++i)
+                point += vectorized_weights[i] * surrounding_points[col + i];
+            }
+          for (std::size_t col = n_regular_cols; col < n_points; ++col)
+            {
+              VectorizedArrayType vectorized_weights;
+              vectorized_weights.gather(&weights(0, 0) + col, offsets.data());
+              point += vectorized_weights * surrounding_points[col];
+            }
+          for (unsigned int r = row;
+               r < std::min<unsigned int>(weights.size(0), row + n_lanes);
+               ++r)
+            {
+              // unpack and project to manifold
+              for (unsigned int d = 0; d < spacedim; ++d)
+                new_points[r][d] = point[d][r - row];
+              new_points[r] =
+                project_to_manifold(surrounding_points, new_points[r]);
+            }
         }
     }
+  else
+    for (unsigned int row = 0; row < weights.size(0); ++row)
+      new_points[row] =
+        get_new_point(surrounding_points,
+                      ArrayView<const double>(&weights(row, 0), n_points));
 }