From 8188832e55710b7128b2a0cb95a11141922477f1 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 11 Apr 2023 13:05:30 +0200 Subject: [PATCH] Alternative version with manual vectorization --- source/grid/manifold.cc | 157 +++++++++++++++------------------------- 1 file changed, 59 insertions(+), 98 deletions(-) diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index 88b2bd776f..d8cacae68a 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -15,6 +15,7 @@ #include #include +#include #include @@ -577,7 +578,7 @@ FlatManifold::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 p; @@ -630,35 +631,6 @@ FlatManifold::get_new_point( -namespace -{ - template - std::array, n_eval_points> - interpolate_points(const unsigned int n_points, - const Point * surrounding_points, - const double * weights, - const Tensor<1, spacedim> &periodicity) - { - std::array, n_eval_points> new_point; - for (unsigned int p = 0; p < n_points; ++p) - { - const Point 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 void FlatManifold::get_new_points( @@ -669,81 +641,70 @@ FlatManifold::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 *surrounding_points_start = surrounding_points.data(); - - boost::container::small_vector, 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, 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(VectorizedArray::size(), 4); + using VectorizedArrayType = VectorizedArray; + 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 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((row + i) * n_points, + (weights.size(0) - 1) * n_points); + Point point; + for (std::size_t col = 0; col < n_regular_cols; col += n_lanes) + { + std::array 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(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(&weights(row, 0), n_points)); } -- 2.39.5