From 2a436880037ba8e5c90e4e75e2cb38acf39fc15d Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 4 Apr 2025 13:15:45 +0200 Subject: [PATCH] Portable MatrixFree: arrange arrays for coalesced access --- .../matrix_free/portable_fe_evaluation.h | 22 ++++++++--------- .../matrix_free/portable_matrix_free.h | 4 ++-- .../portable_matrix_free.templates.h | 24 +++++++++---------- include/deal.II/matrix_free/tools.h | 4 ++-- 4 files changed, 27 insertions(+), 27 deletions(-) diff --git a/include/deal.II/matrix_free/portable_fe_evaluation.h b/include/deal.II/matrix_free/portable_fe_evaluation.h index 607cef6bb5..3b20bd5fbe 100644 --- a/include/deal.II/matrix_free/portable_fe_evaluation.h +++ b/include/deal.II/matrix_free/portable_fe_evaluation.h @@ -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); } } } diff --git a/include/deal.II/matrix_free/portable_matrix_free.h b/include/deal.II/matrix_free/portable_matrix_free.h index c574221888..d92dc4facc 100644 --- a/include/deal.II/matrix_free/portable_matrix_free.h +++ b/include/deal.II/matrix_free/portable_matrix_free.h @@ -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); } }; diff --git a/include/deal.II/matrix_free/portable_matrix_free.templates.h b/include/deal.II/matrix_free/portable_matrix_free.templates.h index c347aa648b..f62e0f3564 100644 --- a/include/deal.II/matrix_free/portable_matrix_free.templates.h +++ b/include/deal.II/matrix_free/portable_matrix_free.templates.h @@ -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( 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( 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]; } } diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index a675381953..1e1aca39ca 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -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))); -- 2.39.5