From 2b4a4ea671aca65d77cd97cabd3ad6d20f276be5 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 4 May 2023 11:15:27 +0200 Subject: [PATCH] Reduce overhead in calls to tensor product value function --- .../deal.II/matrix_free/fe_point_evaluation.h | 43 ++++++++------- .../matrix_free/tensor_product_kernels.h | 53 +++++++++++-------- 2 files changed, 54 insertions(+), 42 deletions(-) diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index d75b8e279d..b22152cd7c 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -1534,14 +1534,11 @@ FEPointEvaluation::do_reinit() if (fast_path && !polynomials_are_hat_functions) { const std::size_t n_batches = - n_q_points_scalar / n_lanes_internal + - (n_q_points_scalar % n_lanes_internal > 0 ? 1 : 0); + (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal; const std::size_t n_shapes = poly.size(); shapes.resize_fast(n_batches * n_shapes); for (unsigned int qb = 0; qb < n_batches; ++qb) - internal::compute_values_of_array(make_array_view(shapes, - qb * n_shapes, - n_shapes), + internal::compute_values_of_array(shapes.data() + qb * n_shapes, poly, unit_point_ptr[qb]); } @@ -1561,13 +1558,15 @@ FEPointEvaluation::evaluate_fast( if (solution_renumbered.size() != dofs_per_component) solution_renumbered.resize(dofs_per_component); for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int i = 0; i < dofs_per_component; ++i) - ETT::read_value( - solution_values[renumber[(component_in_base_element + comp) * - dofs_per_component + - i]], - comp, - solution_renumbered[i]); + { + const unsigned int *renumber_ptr = + renumber.data() + + (component_in_base_element + comp) * dofs_per_component; + for (unsigned int i = 0; i < dofs_per_component; ++i) + ETT::read_value(solution_values[renumber_ptr[i]], + comp, + solution_renumbered[i]); + } // unit gradients are currently only implemented with the fast tensor // path @@ -1586,15 +1585,15 @@ FEPointEvaluation::evaluate_fast( internal::evaluate_tensor_product_value_and_gradient_shapes< dim, scalar_value_type, - VectorizedArrayType>(make_array_view(shapes, - qb * n_shapes, - n_shapes), + VectorizedArrayType>(shapes.data() + qb * n_shapes, n_shapes, solution_renumbered); if (evaluation_flags & EvaluationFlags::values) { - for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v) + for (unsigned int v = 0; + v < stride && (stride > 1 ? q + v < n_q_points_scalar : true); + ++v) ETT::set_value(val_and_grad.first, v, values[qb * stride + v]); } if (evaluation_flags & EvaluationFlags::gradients) @@ -1603,7 +1602,9 @@ FEPointEvaluation::evaluate_fast( update_flags & update_inverse_jacobians, ExcNotInitialized()); - for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v) + for (unsigned int v = 0; + v < stride && (stride > 1 ? q + v < n_q_points_scalar : true); + ++v) { const unsigned int offset = qb * stride + v; ETT::set_gradient(val_and_grad.second, v, unit_gradients[offset]); @@ -1762,7 +1763,9 @@ FEPointEvaluation::integrate_fast( ETT::set_zero_value(values[qb], v); } - for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v) + for (unsigned int v = 0; + v < stride && (stride > 1 ? q + v < n_q_points_scalar : true); + ++v) ETT::get_value(value, v, values[qb * stride + v]); } if (integration_flags & EvaluationFlags::gradients) @@ -1778,7 +1781,9 @@ FEPointEvaluation::integrate_fast( ETT::set_zero_gradient(gradients[qb], v); } - for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v) + for (unsigned int v = 0; + v < stride && (stride > 1 ? q + v < n_q_points_scalar : true); + ++v) { const unsigned int offset = qb * stride + v; ETT::get_gradient( diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 99b69eb97f..db7732e636 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -2989,7 +2989,7 @@ namespace internal template inline void compute_values_of_array( - ArrayView> shapes, + dealii::ndarray * shapes, const std::vector> &poly, const Point & p) { @@ -3000,7 +3000,7 @@ namespace internal for (unsigned int d = 0; d < dim; ++d) point[d] = p[d]; for (int i = 0; i < n_shapes; ++i) - poly[i].values_of_array(point, 1, &shapes[i][0]); + poly[i].values_of_array(point, 1, shapes[i].data()); } @@ -3009,12 +3009,16 @@ namespace internal * Interpolate inner dimensions of tensor product shape functions. */ template - inline std::array::type, 3> - do_interpolate_xy(const std::vector & values, - const std::vector & renumber, - const ArrayView> &shapes, - const int n_shapes_runtime, - int & i) + inline +#ifndef DEBUG + DEAL_II_ALWAYS_INLINE +#endif + std::array::type, 3> + do_interpolate_xy(const std::vector & values, + const std::vector & renumber, + const dealii::ndarray *shapes, + const int n_shapes_runtime, + int & i) { const int n_shapes = length > 0 ? length : n_shapes_runtime; using Number3 = typename ProductTypeNoPoint::type; @@ -3066,10 +3070,10 @@ namespace internal typename ProductTypeNoPoint::type, Tensor<1, dim, typename ProductTypeNoPoint::type>> evaluate_tensor_product_value_and_gradient_shapes( - const ArrayView> &shapes, - const int n_shapes, - const std::vector & values, - const std::vector & renumber = {}) + const dealii::ndarray *shapes, + const int n_shapes, + const std::vector & values, + const std::vector & renumber = {}) { static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented"); @@ -3253,16 +3257,15 @@ namespace internal } else { + AssertIndexRange(poly.size(), 200); std::array, 200> shapes; - auto view = make_array_view(shapes); - - compute_values_of_array(view, poly, p); + compute_values_of_array(shapes.data(), poly, p); return evaluate_tensor_product_value_and_gradient_shapes( - view, poly.size(), values, renumber); + shapes.data(), poly.size(), values, renumber); } } @@ -3469,13 +3472,17 @@ namespace internal * Test inner dimensions of tensor product shape functions and accumulate. */ template - inline void - do_apply_test_functions_xy( - AlignedVector & values, - const ArrayView> &shapes, - const std::array & test_grads_value, - const int n_shapes_runtime, - int & i) + inline +#ifndef DEBUG + DEAL_II_ALWAYS_INLINE +#endif + void + do_apply_test_functions_xy( + AlignedVector & values, + const ArrayView> &shapes, + const std::array & test_grads_value, + const int n_shapes_runtime, + int & i) { if (length > 0) { -- 2.39.5