Portable tensor product kernels: Simplify loop generic

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 7d94f515fb622c736c2d98da7fbc6036d14546d7..f34bc2ff91827ddf2b22260fdd50f412f1d30c4c 100644 (file)
--- a/include/deal.II/matrix_free/portable_tensor_product_kernels.h
+++ b/include/deal.II/matrix_free/portable_tensor_product_kernels.h
@@ -288,20 +288,22 @@ namespace Portable
 
           // This loop simply multiplies the shape function at the quadrature
           // point by the value finite element coefficient.
-          t[q_point] = 0;
-          for (int k = 0; k < n_q_points_1d; ++k)
+          const int stride_shape = dof_to_quad ? n_q_points_1d : 1;
+          const int stride       = Utilities::pow(n_q_points_1d, direction);
+          const int base_shape   = dof_to_quad ? q : (q * n_q_points_1d);
+          const int base =
+            (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);
+          Number sum =
+            shape_data[base_shape] * (in_place ? out(base) : in(base));
+          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_place ? out(source_idx) : in(source_idx));
+              sum +=
+                shape_data[base_shape + k * stride_shape] *
+                (in_place ? out(base + k * stride) : in(base + k * stride));
             }
+          t[q_point] = sum;
         });
 
       if (in_place)