From 310c8cd5f675ba01ae1e161ac61160270ac1e52b Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 21 Apr 2024 10:54:17 +0200 Subject: [PATCH] Portable matrix-free: Fix Kokkos index access --- .../portable_tensor_product_kernels.h | 30 +++++++++---------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/include/deal.II/matrix_free/portable_tensor_product_kernels.h b/include/deal.II/matrix_free/portable_tensor_product_kernels.h index a3c24348c1..7d94f515fb 100644 --- a/include/deal.II/matrix_free/portable_tensor_product_kernels.h +++ b/include/deal.II/matrix_free/portable_tensor_product_kernels.h @@ -132,20 +132,21 @@ namespace Portable n_q_points_1d); Kokkos::parallel_for(thread_policy, [&](const int i, const int j) { int q_point = i + j * n_q_points_1d; - t[q_point] = 0; // This loop simply multiplies the shape function at the quadrature // point by the value finite element coefficient. // FIXME check why using parallel_reduce ThreadVector is slower - for (int k = 0; k < n_q_points_1d; ++k) + const int base_shape = dof_to_quad ? j : j * n_q_points_1d; + const int stride_shape = dof_to_quad ? n_q_points_1d : 1; + const int base_in = (direction == 0) ? (n_q_points_1d * i) : i; + const int stride_in = Utilities::pow(n_q_points_1d, direction); + Number sum = shape_data[base_shape] * in(base_in); + for (int k = 1; k < n_q_points_1d; ++k) { - const unsigned int shape_idx = - dof_to_quad ? (j + k * n_q_points_1d) : (k + j * n_q_points_1d); - const unsigned int source_idx = (direction == 0) ? - (k + n_q_points_1d * i) : - (i + n_q_points_1d * k); - t[q_point] += shape_data[shape_idx] * in(source_idx); + sum += shape_data[base_shape + k * stride_shape] * + in(base_in + k * stride_in); } + t[q_point] = sum; }); if (in_place) @@ -200,20 +201,19 @@ namespace Portable // This loop simply multiplies the shape function at the quadrature // point by the value finite element coefficient. // FIXME check why using parallel_reduce ThreadVector is slower - const int base = dof_to_quad ? q : q * n_q_points_1d; - const int stride_shape = dof_to_quad ? n_q_points_1d : 1; - const Number *in_ptr = - in.data() + + const int base_shape = dof_to_quad ? q : q * n_q_points_1d; + const int stride_shape = dof_to_quad ? n_q_points_1d : 1; + const int base_in = (direction == 0 ? (n_q_points_1d * (i + n_q_points_1d * j)) : (direction == 1 ? (i + n_q_points_1d * n_q_points_1d * j) : (i + n_q_points_1d * j))); const int stride_in = Utilities::pow(n_q_points_1d, direction); - Number sum = shape_data[base] * (*in_ptr); + Number sum = shape_data[base_shape] * in(base_in); for (int k = 1; k < n_q_points_1d; ++k) { - sum += - shape_data[base + k * stride_shape] * in_ptr[k * stride_in]; + sum += shape_data[base_shape + k * stride_shape] * + in(base_in + k * stride_in); } t[q_point] = sum; }); -- 2.39.5