Portable MatrixFree: arrange arrays for coalesced access 18328/head
authorMartin Kronbichler <martin.kronbichler@rub.de>
Fri, 4 Apr 2025 11:15:45 +0000 (13:15 +0200)
committerMartin Kronbichler <martin.kronbichler@rub.de>
Fri, 4 Apr 2025 12:52:52 +0000 (14:52 +0200)
include/deal.II/matrix_free/portable_fe_evaluation.h
include/deal.II/matrix_free/portable_matrix_free.h
include/deal.II/matrix_free/portable_matrix_free.templates.h
include/deal.II/matrix_free/tools.h

index 607cef6bb521dac1e2043ee88004b4e5890e744d..3b20bd5fbe752a4f3f52e5cf79cd1cdf9400d2fe 100644 (file)
@@ -321,7 +321,7 @@ namespace Portable
                            for (unsigned int c = 0; c < n_components_; ++c)
                              shared_data->values(i, c) =
                                src[precomputed_data->local_to_global(
-                                 cell_id, i + tensor_dofs_per_component * c)];
+                                 i + tensor_dofs_per_component * c, cell_id)];
                          });
     data->team_member.team_barrier();
 
@@ -362,7 +362,7 @@ namespace Portable
           [&](const int &i) {
             for (unsigned int c = 0; c < n_components_; ++c)
               dst[precomputed_data->local_to_global(
-                cell_id, i + tensor_dofs_per_component * c)] +=
+                i + tensor_dofs_per_component * c, cell_id)] +=
                 shared_data->values(i, c);
           });
       }
@@ -373,7 +373,7 @@ namespace Portable
           [&](const int &i) {
             for (unsigned int c = 0; c < n_components_; ++c)
               Kokkos::atomic_add(&dst[precomputed_data->local_to_global(
-                                   cell_id, i + (tensor_dofs_per_component)*c)],
+                                   i + (tensor_dofs_per_component)*c, cell_id)],
                                  shared_data->values(i, c));
           });
       }
@@ -542,13 +542,13 @@ namespace Portable
     if constexpr (n_components_ == 1)
       {
         shared_data->values(q_point, 0) =
-          val_in * precomputed_data->JxW(cell_id, q_point);
+          val_in * precomputed_data->JxW(q_point, cell_id);
       }
     else
       {
         for (unsigned int c = 0; c < n_components; ++c)
           shared_data->values(q_point, c) =
-            val_in[c] * precomputed_data->JxW(cell_id, q_point);
+            val_in[c] * precomputed_data->JxW(q_point, cell_id);
       }
   }
 
@@ -600,7 +600,7 @@ namespace Portable
             Number tmp = 0.;
             for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
               tmp +=
-                precomputed_data->inv_jacobian(cell_id, q_point, d_2, d_1) *
+                precomputed_data->inv_jacobian(q_point, cell_id, d_2, d_1) *
                 shared_data->gradients(q_point, d_2, 0);
             grad[d_1] = tmp;
           }
@@ -613,7 +613,7 @@ namespace Portable
               Number tmp = 0.;
               for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
                 tmp +=
-                  precomputed_data->inv_jacobian(cell_id, q_point, d_2, d_1) *
+                  precomputed_data->inv_jacobian(q_point, cell_id, d_2, d_1) *
                   shared_data->gradients(q_point, d_2, c);
               grad[c][d_1] = tmp;
             }
@@ -641,10 +641,10 @@ namespace Portable
             Number tmp = 0.;
             for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
               tmp +=
-                precomputed_data->inv_jacobian(cell_id, q_point, d_1, d_2) *
+                precomputed_data->inv_jacobian(q_point, cell_id, d_1, d_2) *
                 grad_in[d_2];
             shared_data->gradients(q_point, d_1, 0) =
-              tmp * precomputed_data->JxW(cell_id, q_point);
+              tmp * precomputed_data->JxW(q_point, cell_id);
           }
       }
     else
@@ -655,10 +655,10 @@ namespace Portable
               Number tmp = 0.;
               for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
                 tmp +=
-                  precomputed_data->inv_jacobian(cell_id, q_point, d_1, d_2) *
+                  precomputed_data->inv_jacobian(q_point, cell_id, d_1, d_2) *
                   grad_in[c][d_2];
               shared_data->gradients(q_point, d_1, c) =
-                tmp * precomputed_data->JxW(cell_id, q_point);
+                tmp * precomputed_data->JxW(q_point, cell_id);
             }
       }
   }
index c5742218883099989291b8c782c4d8c903bb24ec..d92dc4faccca682f67a05b0ab4ebacc64fe8e2e9 100644 (file)
@@ -404,7 +404,7 @@ namespace Portable
       get_quadrature_point(const unsigned int cell,
                            const unsigned int q_point) const
       {
-        return precomputed_data->q_points(cell, q_point);
+        return precomputed_data->q_points(q_point, cell);
       }
     };
 
@@ -912,7 +912,7 @@ namespace Portable
     get_quadrature_point(const unsigned int cell,
                          const unsigned int q_point) const
     {
-      return q_points(cell, q_point);
+      return q_points(q_point, cell);
     }
   };
 
index c347aa648bd84bd4b4ef364dc3efd18162e59c3c..f62e0f3564fd100d214f5e4453a67fb2da3026c7 100644 (file)
@@ -168,8 +168,8 @@ namespace Portable
                      MemorySpace::Default::kokkos_space>(
           Kokkos::view_alloc("local_to_global_" + std::to_string(color),
                              Kokkos::WithoutInitializing),
-          n_cells,
-          dofs_per_cell);
+          dofs_per_cell,
+          n_cells);
 
       if (update_flags & update_quadrature_points)
         data->q_points[color] =
@@ -177,24 +177,24 @@ namespace Portable
                        MemorySpace::Default::kokkos_space>(
             Kokkos::view_alloc("q_points_" + std::to_string(color),
                                Kokkos::WithoutInitializing),
-            n_cells,
-            q_points_per_cell);
+            q_points_per_cell,
+            n_cells);
 
       if (update_flags & update_JxW_values)
         data->JxW[color] =
           Kokkos::View<Number **, MemorySpace::Default::kokkos_space>(
             Kokkos::view_alloc("JxW_" + std::to_string(color),
                                Kokkos::WithoutInitializing),
-            n_cells,
-            q_points_per_cell);
+            q_points_per_cell,
+            n_cells);
 
       if (update_flags & update_gradients)
         data->inv_jacobian[color] =
           Kokkos::View<Number **[dim][dim], MemorySpace::Default::kokkos_space>(
             Kokkos::view_alloc("inv_jacobian_" + std::to_string(color),
                                Kokkos::WithoutInitializing),
-            n_cells,
-            q_points_per_cell);
+            q_points_per_cell,
+            n_cells);
 
       // Initialize to zero, i.e., unconstrained cell
       data->constraint_mask[color] =
@@ -265,7 +265,7 @@ namespace Portable
                                           cell_id_view);
 
           for (unsigned int i = 0; i < dofs_per_cell; ++i)
-            local_to_global_host(cell_id, i) = lexicographic_dof_indices[i];
+            local_to_global_host(i, cell_id) = lexicographic_dof_indices[i];
 
           fe_values.reinit(*cell);
 
@@ -273,13 +273,13 @@ namespace Portable
           if (update_flags & update_quadrature_points)
             {
               for (unsigned int i = 0; i < q_points_per_cell; ++i)
-                q_points_host(cell_id, i) = fe_values.quadrature_point(i);
+                q_points_host(i, cell_id) = fe_values.quadrature_point(i);
             }
 
           if (update_flags & update_JxW_values)
             {
               for (unsigned int i = 0; i < q_points_per_cell; ++i)
-                JxW_host(cell_id, i) = fe_values.JxW(i);
+                JxW_host(i, cell_id) = fe_values.JxW(i);
             }
 
           if (update_flags & update_gradients)
@@ -287,7 +287,7 @@ namespace Portable
               for (unsigned int i = 0; i < q_points_per_cell; ++i)
                 for (unsigned int d = 0; d < dim; ++d)
                   for (unsigned int e = 0; e < dim; ++e)
-                    inv_jacobian_host(cell_id, i, d, e) =
+                    inv_jacobian_host(i, cell_id, d, e) =
                       fe_values.inverse_jacobian(i)[d][e];
             }
         }
index a675381953eae6ca25f10c8f143adcafca706845..1e1aca39cae2bd2ab7027a9a308a90a5fc75d916 100644 (file)
@@ -1478,7 +1478,7 @@ namespace MatrixFreeTools
             Kokkos::parallel_for(
               Kokkos::TeamThreadRange(data->team_member, dofs_per_cell),
               [&](const int &i) {
-                dst[gpu_data->local_to_global(cell, i)] +=
+                dst[gpu_data->local_to_global(i, cell)] +=
                   data->shared_data->values(i % (dofs_per_cell / n_components),
                                             i / (dofs_per_cell / n_components));
               });
@@ -1488,7 +1488,7 @@ namespace MatrixFreeTools
             Kokkos::parallel_for(
               Kokkos::TeamThreadRange(data->team_member, dofs_per_cell),
               [&](const int &i) {
-                Kokkos::atomic_add(&dst[gpu_data->local_to_global(cell, i)],
+                Kokkos::atomic_add(&dst[gpu_data->local_to_global(i, cell)],
                                    data->shared_data->values(
                                      i % (dofs_per_cell / n_components),
                                      i / (dofs_per_cell / n_components)));

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.