From 87e2a81a5d1010ea94ee3a74279c9d9c81755059 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 31 Mar 2023 19:48:45 +0200 Subject: [PATCH] FEPointEvaluation::integrate: Choose better loop layout --- .../matrix_free/tensor_product_kernels.h | 94 +++++++++++-------- 1 file changed, 56 insertions(+), 38 deletions(-) diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 74b7279fe2..08d09b328e 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -3376,34 +3376,60 @@ namespace internal inline void do_apply_test_functions_xy( AlignedVector & values, - const std::vector & renumber, const ArrayView> &shapes, const std::array & test_grads_value, const int n_shapes_runtime, int & i) { - const int n_shapes = length > 0 ? length : n_shapes_runtime; - for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1) + if (length > 0) { - const Number2 test_value_y = - dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] + - test_grads_value[1] * shapes[i1][1][1]) : - test_grads_value[2]; - const Number2 test_grad_xy = dim > 1 ? - test_grads_value[0] * shapes[i1][0][1] : - test_grads_value[0]; - if (renumber.empty()) - for (int i0 = 0; i0 < n_shapes; ++i0, ++i) - { - values[i] += shapes[i0][0][0] * test_value_y; - values[i] += shapes[i0][1][0] * test_grad_xy; - } - else - for (int i0 = 0; i0 < n_shapes; ++i0, ++i) - { - values[renumber[i]] += shapes[i0][0][0] * test_value_y; - values[renumber[i]] += shapes[i0][1][0] * test_grad_xy; - } + constexpr unsigned int array_size = length > 0 ? length : 1; + std::array shape_values_x; + std::array shape_derivs_x; + for (unsigned int i = 0; i < array_size; ++i) + { + shape_values_x[i] = shapes[i][0][0]; + shape_derivs_x[i] = shapes[i][1][0]; + } + for (unsigned int i1 = 0; i1 < (dim > 1 ? length : 1); ++i1) + { + const Number2 test_value_y = + dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] + + test_grads_value[1] * shapes[i1][1][1]) : + test_grads_value[2]; + const Number2 test_grad_xy = + dim > 1 ? test_grads_value[0] * shapes[i1][0][1] : + test_grads_value[0]; + + Number2 *values_ptr = values.data() + i + i1 * length; + for (unsigned int i0 = 0; i0 < length; ++i0) + { + values_ptr[i0] += shape_values_x[i0] * test_value_y; + values_ptr[i0] += shape_derivs_x[i0] * test_grad_xy; + } + } + i += (dim > 1 ? length * length : length); + } + else + { + for (int i1 = 0; i1 < (dim > 1 ? n_shapes_runtime : 1); ++i1) + { + const Number2 test_value_y = + dim > 1 ? (test_grads_value[2] * shapes[i1][0][1] + + test_grads_value[1] * shapes[i1][1][1]) : + test_grads_value[2]; + const Number2 test_grad_xy = + dim > 1 ? test_grads_value[0] * shapes[i1][0][1] : + test_grads_value[0]; + + Number2 *values_ptr = values.data() + i + i1 * n_shapes_runtime; + for (int i0 = 0; i0 < n_shapes_runtime; ++i0) + { + values_ptr[i0] += shapes[i0][0][0] * test_value_y; + values_ptr[i0] += shapes[i0][1][0] * test_grad_xy; + } + } + i += (dim > 1 ? n_shapes_runtime * n_shapes_runtime : n_shapes_runtime); } } @@ -3419,15 +3445,12 @@ namespace internal const int n_shapes, const Number2 & value, const Tensor<1, dim, Number2> & gradient, - AlignedVector & values, - const std::vector & renumber = {}) + AlignedVector & values) { static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented"); // as in evaluate, use `int` type to produce better code in this context AssertDimension(Utilities::pow(n_shapes, dim), values.size()); - Assert(renumber.empty() || renumber.size() == values.size(), - ExcDimensionMismatch(renumber.size(), values.size())); // Implement the transpose of the function above std::array test_grads_value; @@ -3449,22 +3472,22 @@ namespace internal // cases if (n_shapes == 2) do_apply_test_functions_xy( - values, renumber, shapes, test_grads_value, n_shapes, i); + values, shapes, test_grads_value, n_shapes, i); else if (n_shapes == 3) do_apply_test_functions_xy( - values, renumber, shapes, test_grads_value, n_shapes, i); + values, shapes, test_grads_value, n_shapes, i); else if (n_shapes == 4) do_apply_test_functions_xy( - values, renumber, shapes, test_grads_value, n_shapes, i); + values, shapes, test_grads_value, n_shapes, i); else if (n_shapes == 5) do_apply_test_functions_xy( - values, renumber, shapes, test_grads_value, n_shapes, i); + values, shapes, test_grads_value, n_shapes, i); else if (n_shapes == 6) do_apply_test_functions_xy( - values, renumber, shapes, test_grads_value, n_shapes, i); + values, shapes, test_grads_value, n_shapes, i); else do_apply_test_functions_xy( - values, renumber, shapes, test_grads_value, n_shapes, i); + values, shapes, test_grads_value, n_shapes, i); } } @@ -3481,19 +3504,14 @@ namespace internal const Number2 & value, const Tensor<1, dim, Number2> & gradient, AlignedVector & values, - const Point & p, - const std::vector & renumber = {}) + const Point & p) { (void)poly; static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented"); AssertDimension(Utilities::pow(poly.size(), dim), values.size()); - Assert(renumber.empty() || renumber.size() == values.size(), - ExcDimensionMismatch(renumber.size(), values.size())); AssertDimension(poly.size(), 2); - for (unsigned int i = 0; i < renumber.size(); ++i) - AssertDimension(renumber[i], i); if (dim == 1) { -- 2.39.5