From 648785bdb07e368dee609d8b39f135d9803b4206 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 21 Apr 2022 11:01:26 +0200 Subject: [PATCH] Add new renumbering functions for MatrixFree data locality --- include/deal.II/dofs/dof_renumbering.h | 65 ++++ source/dofs/dof_renumbering.cc | 399 +++++++++++++++++++++++++ source/dofs/dof_renumbering.inst.in | 82 +++++ 3 files changed, 546 insertions(+) diff --git a/include/deal.II/dofs/dof_renumbering.h b/include/deal.II/dofs/dof_renumbering.h index 6aa9a9a5a6..3fb0513451 100644 --- a/include/deal.II/dofs/dof_renumbering.h +++ b/include/deal.II/dofs/dof_renumbering.h @@ -220,6 +220,15 @@ DEAL_II_NAMESPACE_OPEN * (iterative or direct ones) on the numbering of the degrees of freedom. * * + *

Renumbering DoFs for faster matrix-free computations

+ * + * The MatrixFree class provides optimized algorithms for interleaving + * operations on vectors before and after the access of the vector data in the + * respective loops. The algorithm matrix_free_data_locality() makes sure + * that all unknowns with a short distance between the first and last access + * are grouped together, in order to increase the spatial data locality. + * + * *

A comparison of reordering strategies

* * As a benchmark of comparison, let us consider what the different sparsity @@ -1264,6 +1273,62 @@ namespace DoFRenumbering * @} */ + /** + * @name Numberings for better performance with the MatrixFree infrastructure + * @{ + */ + + /** + * Sort DoFs by their appearance in matrix-free loops in order to increase + * data locality. More specifically, this renumbering strategy will group + * DoFs touched only on a single group of cells (as traversed by a + * MatrixFree::cell_loop()) nearby, whereas DoFs touched far apart through + * the loop over cells are grouped separately. Obviously, also unknowns that + * are subject to ghost exchange will be treated as those with far reach, + * and will be grouped accordingly. Currently, this function only works for + * finite elements of type FE_Q as well as systems of a single FE_Q element. + * + * This function needs to be provided with an appropriate + * MatrixFree::AdditionalData structure that indicates whether level DoFs or + * active DoFs are to be renumbered through the variable + * MatrixFree::AdditionalData::mg_level, and with the options that determine + * the order cells are traversed in a matrix-free loop. + * + * Note that the argument `matrix_free` does not need to be initialized when + * calling this function. In case it is already set up, the unknowns stored + * in that data will be used. If it is not set up, the + * MatrixFree::AdditionalData is instead used to construct an object with + * the appropriate indices, which are then processed. In both cases, the + * content of MatrixFree is left untouched. However, both MatrixFree and the + * associated AffineConstraints object need to be computed again after + * changing the numbers of the degrees of freedom. + */ + template + void + matrix_free_data_locality( + DoFHandler & dof_handler, + const AffineConstraints & constraints, + const MatrixFreeType & matrix_free, + const typename MatrixFreeType::AdditionalData &matrix_free_data); + + /** + * Compute the renumbering vector needed by the matrix_free_data_locality() + * function. + * Does not perform the renumbering on the @p DoFHandler dofs but returns the + * renumbering vector. + */ + template + std::vector + compute_matrix_free_data_locality( + const DoFHandler & dof_handler, + const AffineConstraints & constraints, + const MatrixFreeType & matrix_free, + const typename MatrixFreeType::AdditionalData &matrix_free_data); + + /** + * @} + */ + /** * Exception diff --git a/source/dofs/dof_renumbering.cc b/source/dofs/dof_renumbering.cc index a371d9a270..aa571f55b6 100644 --- a/source/dofs/dof_renumbering.cc +++ b/source/dofs/dof_renumbering.cc @@ -39,6 +39,8 @@ #include #include +#include + #include #include @@ -2435,6 +2437,403 @@ namespace DoFRenumbering new_dof_indices[i] = component_to_nodal[local_index]; } } + + + + template + void + matrix_free_data_locality( + DoFHandler & dof_handler, + const AffineConstraints & constraints, + const MatrixFreeType & matrix_free, + const typename MatrixFreeType::AdditionalData &matrix_free_data) + { + const std::vector new_global_numbers = + compute_matrix_free_data_locality(dof_handler, + constraints, + matrix_free, + matrix_free_data); + if (matrix_free_data.mg_level == dealii::numbers::invalid_unsigned_int) + dof_handler.renumber_dofs(new_global_numbers); + else + dof_handler.renumber_dofs(matrix_free_data.mg_level, new_global_numbers); + } + + + + // Implementation details for matrix-free renumbering + namespace + { + std::vector> + group_dofs_by_rank_access( + const dealii::Utilities::MPI::Partitioner &partitioner) + { + // count the number of times a locally owned DoF is accessed by the + // remote ghost data + std::vector touch_count(partitioner.locally_owned_size()); + for (const auto &p : partitioner.import_indices()) + for (unsigned int i = p.first; i < p.second; ++i) + touch_count[i]++; + + // category 0: DoFs never touched by ghosts + std::vector> result(1); + for (unsigned int i = 0; i < touch_count.size(); ++i) + if (touch_count[i] == 0) + result.back().push_back(i); + + // DoFs with 1 appearance can be simply categorized by their (single) + // MPI rank, whereas we need to go an extra round for the remaining DoFs + // by collecting the owning processes by hand + std::map> + multiple_ranks_access_dof; + const std::vector> &import_targets = + partitioner.import_targets(); + auto it = partitioner.import_indices().begin(); + for (const std::pair &proc : import_targets) + { + result.emplace_back(); + unsigned int count_dofs = 0; + while (count_dofs < proc.second) + { + for (unsigned int i = it->first; i < it->second; + ++i, ++count_dofs) + { + if (touch_count[i] == 1) + result.back().push_back(i); + else + multiple_ranks_access_dof[i].push_back(proc.first); + } + ++it; + } + } + Assert(it == partitioner.import_indices().end(), ExcInternalError()); + + // Now go from the computed map on DoFs to a map on the processes for + // DoFs with multiple owners, and append the DoFs found this way to our + // global list + std::map, + std::vector, + std::function &, + const std::vector &)>> + dofs_by_rank{[](const std::vector &a, + const std::vector &b) { + if (a.size() < b.size()) + return true; + if (a.size() == b.size()) + { + for (unsigned int i = 0; i < a.size(); ++i) + if (a[i] < b[i]) + return true; + else if (a[i] > b[i]) + return false; + } + return false; + }}; + for (const auto &entry : multiple_ranks_access_dof) + dofs_by_rank[entry.second].push_back(entry.first); + + for (const auto &procs : dofs_by_rank) + result.push_back(procs.second); + + return result; + } + + + + template + std::pair, std::vector> + compute_mf_numbering(const MatrixFree &mf, + const unsigned int component) + { + IndexSet owned_dofs = + mf.get_dof_info(component).vector_partitioner->locally_owned_range(); + const unsigned int n_comp = + mf.get_dof_handler(component).get_fe().n_components(); + Assert(mf.get_dof_handler(component).get_fe().n_base_elements() == 1, + ExcNotImplemented()); + Assert(dynamic_cast *>( + &mf.get_dof_handler(component).get_fe().base_element(0)), + ExcNotImplemented("Matrix-free renumbering only works for " + "FE_Q elements")); + + const unsigned int fe_degree = + mf.get_dof_handler(component).get_fe().degree; + const unsigned int nn = fe_degree - 1; + + // Data structure used further down for the identification of various + // entities in the hierarchical numbering of unknowns + std::array, + dealii::Utilities::pow(3, dim)> + dofs_on_objects; + if (dim == 1) + { + dofs_on_objects[0] = std::make_pair(0U, 1U); + dofs_on_objects[1] = std::make_pair(2 * n_comp, nn); + dofs_on_objects[2] = std::make_pair(n_comp, 1U); + } + else if (dim == 2) + { + dofs_on_objects[0] = std::make_pair(0U, 1U); + dofs_on_objects[1] = std::make_pair(n_comp * (4 + 2 * nn), nn); + dofs_on_objects[2] = std::make_pair(n_comp, 1U); + dofs_on_objects[3] = std::make_pair(n_comp * 4, nn); + dofs_on_objects[4] = std::make_pair(n_comp * (4 + 4 * nn), nn * nn); + dofs_on_objects[5] = std::make_pair(n_comp * (4 + 1 * nn), nn); + dofs_on_objects[6] = std::make_pair(2 * n_comp, 1U); + dofs_on_objects[7] = std::make_pair(n_comp * (4 + 3 * nn), nn); + dofs_on_objects[8] = std::make_pair(3 * n_comp, 1U); + } + else if (dim == 3) + { + dofs_on_objects[0] = std::make_pair(0U, 1U); + dofs_on_objects[1] = std::make_pair(n_comp * (8 + 2 * nn), nn); + dofs_on_objects[2] = std::make_pair(n_comp, 1U); + dofs_on_objects[3] = std::make_pair(n_comp * 8, nn); + dofs_on_objects[4] = + std::make_pair(n_comp * (8 + 12 * nn + 4 * nn * nn), nn * nn); + dofs_on_objects[5] = std::make_pair(n_comp * (8 + 1 * nn), nn); + dofs_on_objects[6] = std::make_pair(n_comp * 2, 1U); + dofs_on_objects[7] = std::make_pair(n_comp * (8 + 3 * nn), nn); + dofs_on_objects[8] = std::make_pair(n_comp * 3, 1U); + dofs_on_objects[9] = std::make_pair(n_comp * (8 + 8 * nn), nn); + dofs_on_objects[10] = + std::make_pair(n_comp * (8 + 12 * nn + 2 * nn * nn), nn * nn); + dofs_on_objects[11] = std::make_pair(n_comp * (8 + 9 * nn), nn); + dofs_on_objects[12] = std::make_pair(n_comp * (8 + 12 * nn), nn * nn); + dofs_on_objects[13] = + std::make_pair(n_comp * (8 + 12 * nn + 6 * nn * nn), nn * nn * nn); + dofs_on_objects[14] = + std::make_pair(n_comp * (8 + 12 * nn + 1 * nn * nn), nn * nn); + dofs_on_objects[15] = std::make_pair(n_comp * (8 + 10 * nn), nn); + dofs_on_objects[16] = + std::make_pair(n_comp * (8 + 12 * nn + 3 * nn * nn), nn * nn); + dofs_on_objects[17] = std::make_pair(n_comp * (8 + 11 * nn), nn); + dofs_on_objects[18] = std::make_pair(n_comp * 4, 1U); + dofs_on_objects[19] = std::make_pair(n_comp * (8 + 6 * nn), nn); + dofs_on_objects[20] = std::make_pair(n_comp * 5, 1U); + dofs_on_objects[21] = std::make_pair(n_comp * (8 + 4 * nn), nn); + dofs_on_objects[22] = + std::make_pair(n_comp * (8 + 12 * nn + 5 * nn * nn), nn * nn); + dofs_on_objects[23] = std::make_pair(n_comp * (8 + 5 * nn), nn); + dofs_on_objects[24] = std::make_pair(n_comp * 6, 1U); + dofs_on_objects[25] = std::make_pair(n_comp * (8 + 7 * nn), nn); + dofs_on_objects[26] = std::make_pair(n_comp * 7, 1U); + } + + const auto renumber_func = [](const types::global_dof_index dof_index, + const IndexSet & owned_dofs, + std::vector & result, + unsigned int &counter_dof_numbers) { + if (owned_dofs.is_element(dof_index)) + { + const unsigned int local_dof_index = + owned_dofs.index_within_set(dof_index); + AssertIndexRange(local_dof_index, result.size()); + if (result[local_dof_index] == numbers::invalid_unsigned_int) + result[local_dof_index] = counter_dof_numbers++; + } + }; + + unsigned int counter_dof_numbers = 0; + std::vector local_dofs, dofs_extracted; + std::vector dof_indices( + mf.get_dof_handler(component).get_fe().dofs_per_cell); + + // We now define a lambda function that does two things: (a) determine + // DoF numbers in a way that fit with the order a MatrixFree loop + // travels through the cells (variable 'dof_numbers_mf_order'), and (b) + // determine which unknowns are only touched from within a single range + // of cells and which ones span multiple ranges (variable + // 'touch_count'). Note that this process is done by calling into + // MatrixFree::cell_loop, which gives the right level of granularity + // (when executed in serial) for the scheduled vector operations. Note + // that we pick the unconstrained indices in the hierarchical order for + // part (a) as this makes it easy to identify the DoFs on the individual + // entities, whereas we pick the numbers with constraints eliminated for + // part (b). For the latter, we finally need to sort the indices and + // remove duplicates to really only count the DoFs in each range of cell + // batches once. + std::vector dof_numbers_mf_order( + owned_dofs.n_elements(), dealii::numbers::invalid_unsigned_int); + std::vector touch_count(owned_dofs.n_elements()); + + const auto operation_on_cell_range = + [&](const MatrixFree &data, + unsigned int &, + const unsigned int &, + const std::pair &cell_range) { + local_dofs.clear(); + for (unsigned int cell = cell_range.first; cell < cell_range.second; + ++cell) + { + // part (a): assign beneficial numbers + for (unsigned int v = 0; + v < data.n_active_entries_per_cell_batch(cell); + ++v) + { + // get the indices for the dofs in cell_batch + if (mf.get_mg_level() == numbers::invalid_unsigned_int) + data.get_cell_iterator(cell, v)->get_dof_indices( + dof_indices); + else + data.get_cell_iterator(cell, v)->get_mg_dof_indices( + dof_indices); + + for (unsigned int a = 0; a < dofs_on_objects.size(); ++a) + { + const auto &r = dofs_on_objects[a]; + if (a == 10 || a == 16) + // switch order x-z for y faces in 3D to lexicographic + // layout + for (unsigned int i1 = 0; i1 < nn; ++i1) + for (unsigned int i0 = 0; i0 < nn; ++i0) + for (unsigned int c = 0; c < n_comp; ++c) + renumber_func(dof_indices[r.first + r.second * c + + i1 + i0 * nn], + owned_dofs, + dof_numbers_mf_order, + counter_dof_numbers); + else + for (unsigned int i = 0; i < r.second; ++i) + for (unsigned int c = 0; c < n_comp; ++c) + renumber_func( + dof_indices[r.first + r.second * c + i], + owned_dofs, + dof_numbers_mf_order, + counter_dof_numbers); + } + } + + // part (b): compute the touch count in the current cell batch + // range + data.get_dof_info(component).get_dof_indices_on_cell_batch( + dofs_extracted, cell); + local_dofs.insert(local_dofs.end(), + dofs_extracted.begin(), + dofs_extracted.end()); + } + + std::sort(local_dofs.begin(), local_dofs.end()); + local_dofs.erase(std::unique(local_dofs.begin(), local_dofs.end()), + local_dofs.end()); + + for (unsigned int i : local_dofs) + if (i < touch_count.size()) + touch_count[i]++; + }; + + // Finally run the matrix-free loop. + Assert(mf.get_task_info().scheme == + dealii::internal::MatrixFreeFunctions::TaskInfo::none, + ExcNotImplemented("Renumbering only available for non-threaded " + "version of MatrixFree::cell_loop")); + + mf.template cell_loop(operation_on_cell_range, + counter_dof_numbers, + counter_dof_numbers); + + AssertDimension(counter_dof_numbers, owned_dofs.n_elements()); + + return std::make_pair(dof_numbers_mf_order, touch_count); + } + + } // namespace + + + + template + std::vector + compute_matrix_free_data_locality( + const DoFHandler & dof_handler, + const AffineConstraints & constraints, + const MatrixFreeType & matrix_free, + const typename MatrixFreeType::AdditionalData &matrix_free_data) + { + typename MatrixFreeType::AdditionalData my_mf_data = matrix_free_data; + my_mf_data.initialize_mapping = false; + my_mf_data.tasks_parallel_scheme = MatrixFreeType::AdditionalData::none; + + // try to locate the `DoFHandler` within the given MatrixFree object. + unsigned int component = 0; + for (; component < matrix_free.n_components(); ++component) + if (&matrix_free.get_dof_handler(component) == &dof_handler && + matrix_free.indices_initialized()) + break; + + // if not found construct a new one + const MatrixFreeType *chosen_matrix_free = &matrix_free; + MatrixFreeType separate_matrix_free; + if (component == matrix_free.n_components()) + { + separate_matrix_free.reinit(dealii::MappingQ1(), + dof_handler, + constraints, + dealii::QGauss<1>(2), + my_mf_data); + chosen_matrix_free = &separate_matrix_free; + component = 0; + } + + const std::vector> dofs_by_rank_access = + group_dofs_by_rank_access( + *chosen_matrix_free->get_dof_info(component).vector_partitioner); + + const std::pair, std::vector> + local_numbering = compute_mf_numbering(*chosen_matrix_free, component); + + // Now construct the new numbering + const IndexSet &owned_dofs = chosen_matrix_free->get_dof_info(component) + .vector_partitioner->locally_owned_range(); + std::vector new_numbers; + new_numbers.reserve(owned_dofs.n_elements()); + + // step 1: Take all DoFs with reference only to the current MPI process + // and touched once ("perfect overlap" case). We define a custom + // comparator for std::sort to then order the unknowns by the specified + // matrix-free loop order + const std::vector &purely_local_dofs = dofs_by_rank_access[0]; + for (unsigned int i : purely_local_dofs) + if (local_numbering.second[i] == 1) + new_numbers.push_back(i); + + const auto comparator = [&](const unsigned int a, const unsigned int b) { + return (local_numbering.first[a] < local_numbering.first[b]); + }; + + std::sort(new_numbers.begin(), new_numbers.end(), comparator); + unsigned int sorted_size = new_numbers.size(); + + // step 2: Take all DoFs with reference to only the current MPI process + // and touched multiple times (more strain on caches). + for (auto i : purely_local_dofs) + if (local_numbering.second[i] > 1) + new_numbers.push_back(i); + std::sort(new_numbers.begin() + sorted_size, new_numbers.end(), comparator); + sorted_size = new_numbers.size(); + + // step 3: Get all DoFs with reference from other MPI ranks + for (unsigned int chunk = 1; chunk < dofs_by_rank_access.size(); ++chunk) + for (auto i : dofs_by_rank_access[chunk]) + new_numbers.push_back(i); + std::sort(new_numbers.begin() + sorted_size, new_numbers.end(), comparator); + sorted_size = new_numbers.size(); + + // step 4: Put all DoFs without any reference (constrained DoFs) + for (auto i : purely_local_dofs) + if (local_numbering.second[i] == 0) + new_numbers.push_back(i); + std::sort(new_numbers.begin() + sorted_size, new_numbers.end(), comparator); + + AssertDimension(new_numbers.size(), owned_dofs.n_elements()); + + std::vector new_global_numbers( + owned_dofs.n_elements()); + for (unsigned int i = 0; i < new_numbers.size(); ++i) + new_global_numbers[new_numbers[i]] = owned_dofs.nth_index_in_set(i); + + return new_global_numbers; + } + } // namespace DoFRenumbering diff --git a/source/dofs/dof_renumbering.inst.in b/source/dofs/dof_renumbering.inst.in index 3f97733f7e..72472b5c2d 100644 --- a/source/dofs/dof_renumbering.inst.in +++ b/source/dofs/dof_renumbering.inst.in @@ -235,3 +235,85 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) \} #endif } + + +for (deal_II_dimension : DIMENSIONS; + deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED) + { + namespace DoFRenumbering + \{ + template void + matrix_free_data_locality< + deal_II_dimension, + deal_II_dimension, + deal_II_scalar_vectorized::value_type, + MatrixFree>( + DoFHandler &, + const AffineConstraints &, + const MatrixFree &, + const typename MatrixFree::AdditionalData &); + + template std::vector + compute_matrix_free_data_locality< + deal_II_dimension, + deal_II_dimension, + deal_II_scalar_vectorized::value_type, + MatrixFree>( + const DoFHandler &, + const AffineConstraints &, + const MatrixFree &, + const MatrixFree::AdditionalData &); + \} + } + +for (deal_II_dimension : DIMENSIONS; + deal_II_float_vectorized : FLOAT_VECTORIZED) + { + namespace DoFRenumbering + \{ + template void + matrix_free_data_locality>( + DoFHandler &, + const AffineConstraints &, + const MatrixFree &, + const typename MatrixFree::AdditionalData &); + + template std::vector + compute_matrix_free_data_locality< + deal_II_dimension, + deal_II_dimension, + double, + MatrixFree>( + const DoFHandler &, + const AffineConstraints &, + const MatrixFree &, + const typename MatrixFree::AdditionalData &); + \} + } \ No newline at end of file -- 2.39.5