From: Martin Kronbichler Date: Sat, 20 Apr 2024 06:25:29 +0000 (+0200) Subject: Portable tensor product kernels: Simplify loop X-Git-Tag: v9.6.0-rc1~359^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bbfcbea72340a71efbbfdcd1617f17d3b2cae5f5;p=dealii.git Portable tensor product kernels: Simplify loop --- 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 f67730ccdc..a3c24348c1 100644 --- a/include/deal.II/matrix_free/portable_tensor_product_kernels.h +++ b/include/deal.II/matrix_free/portable_tensor_product_kernels.h @@ -196,23 +196,26 @@ namespace Portable thread_policy, [&](const int i, const int j, const int q) { const int q_point = i + j * n_q_points_1d + q * n_q_points_1d * 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 = 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() + + (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); + for (int k = 1; k < n_q_points_1d; ++k) { - const unsigned int shape_idx = - dof_to_quad ? (q + k * n_q_points_1d) : (k + q * n_q_points_1d); - const unsigned int source_idx = - (direction == 0) ? - (k + n_q_points_1d * (i + n_q_points_1d * j)) : - (direction == 1) ? - (i + n_q_points_1d * (k + n_q_points_1d * j)) : - (i + n_q_points_1d * (j + n_q_points_1d * k)); - t[q_point] += shape_data[shape_idx] * in(source_idx); + sum += + shape_data[base + k * stride_shape] * in_ptr[k * stride_in]; } + t[q_point] = sum; }); if (in_place)