From: Martin Kronbichler Date: Sat, 2 May 2020 13:52:12 +0000 (+0200) Subject: Group cells by parents in matrix-free setup X-Git-Tag: v9.2.0-rc1~79^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9e7a935c8ed5fc88b61de0fc531b3ce51eb54708;p=dealii.git Group cells by parents in matrix-free setup --- diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index ca57050fc0..f26f18b63b 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -460,11 +460,11 @@ MatrixFree::internal_reinit( initialize_dof_handlers(dof_handler, additional_data); std::vector dummy; std::vector dummy2; - task_info.collect_boundary_cells(cell_level_index.size(), - cell_level_index.size(), - VectorizedArrayType::size(), - dummy); - task_info.create_blocks_serial(dummy, 1, dummy, false, dummy, dummy2); + task_info.vectorization_length = VectorizedArrayType::size(); + task_info.n_active_cells = cell_level_index.size(); + task_info.create_blocks_serial( + dummy, 1, dummy, false, dummy, dummy, dummy2); + for (unsigned int i = 0; i < dof_info.size(); ++i) { Assert(dof_handler[i]->get_fe_collection().size() == 1, @@ -970,11 +970,10 @@ MatrixFree::initialize_indices( subdomain_boundary_cells.push_back(counter); } - const unsigned int n_lanes = VectorizedArrayType::size(); - task_info.collect_boundary_cells(cell_level_index_end_local, - n_active_cells, - n_lanes, - subdomain_boundary_cells); + const unsigned int n_lanes = VectorizedArrayType::size(); + task_info.n_active_cells = cell_level_index_end_local; + task_info.n_ghost_cells = n_active_cells - cell_level_index_end_local; + task_info.vectorization_length = n_lanes; // Finalize the creation of the ghost indices { @@ -1022,21 +1021,53 @@ MatrixFree::initialize_indices( std::vector irregular_cells; if (task_info.scheme == internal::MatrixFreeFunctions::TaskInfo::none) { - const bool strict_categories = + bool strict_categories = additional_data.cell_vectorization_categories_strict || dof_handlers.active_dof_handler == DoFHandlers::hp; unsigned int dofs_per_cell = 0; for (const auto &info : dof_info) dofs_per_cell = std::max(dofs_per_cell, info.dofs_per_cell[0]); + + // Detect cells with the same parent to make sure they get scheduled + // together in the loop, which increases data locality. + std::vector parent_relation(task_info.n_active_cells + + task_info.n_ghost_cells, + numbers::invalid_unsigned_int); + std::map, std::vector> cell_parents; + for (unsigned int c = 0; c < cell_level_index_end_local; ++c) + if (cell_level_index[c].first > 0) + { + typename Triangulation::cell_iterator cell( + dof_handlers.active_dof_handler == DoFHandlers::usual ? + &dof_handlers.dof_handler[0]->get_triangulation() : + &dof_handlers.hp_dof_handler[0]->get_triangulation(), + cell_level_index[c].first, + cell_level_index[c].second); + Assert(cell->level() > 0, ExcInternalError()); + cell_parents[std::make_pair(cell->parent()->level(), + cell->parent()->index())] + .push_back(c); + } + unsigned int position = 0; + for (const auto &it : cell_parents) + if (it.second.size() == GeometryInfo::max_children_per_cell) + { + for (auto i : it.second) + parent_relation[i] = position; + ++position; + } task_info.create_blocks_serial(subdomain_boundary_cells, dofs_per_cell, dof_info[0].cell_active_fe_index, strict_categories, + parent_relation, renumbering, irregular_cells); } else { + task_info.make_boundary_cells_divisible(subdomain_boundary_cells); + // For strategy with blocking before partitioning: reorganize the indices // in order to overlap communication in MPI with computations: Place all // cells with ghost indices into one chunk. Also reorder cells so that we diff --git a/include/deal.II/matrix_free/task_info.h b/include/deal.II/matrix_free/task_info.h index cfc972ccd9..ed9a128b02 100644 --- a/include/deal.II/matrix_free/task_info.h +++ b/include/deal.II/matrix_free/task_info.h @@ -135,14 +135,11 @@ namespace internal loop(MFWorkerInterface &worker) const; /** - * Determines the position of cells with ghosts for distributed-memory - * calculations. + * Make the number of cells which can only be treated in the + * communication overlap divisible by the vectorization length. */ void - collect_boundary_cells(const unsigned int n_active_cells, - const unsigned int n_active_and_ghost_cells, - const unsigned int vectorization_length, - std::vector &boundary_cells); + make_boundary_cells_divisible(std::vector &boundary_cells); /** * Sets up the blocks for running the cell loop based on the options @@ -166,6 +163,10 @@ namespace internal * strictly or whether it is allowed to insert lower categories into the * next high one(s). * + * @param parent_relation This data field is used to specify which cells + * have the same parent cell. Cells with the same ancestor are grouped + * together into the same batch(es) with vectorization across cells. + * * @param renumbering When leaving this function, the vector contains a * new numbering of the cells that aligns with the grouping stored in * this class. @@ -182,6 +183,7 @@ namespace internal const unsigned int dofs_per_cell, const std::vector &cell_vectorization_categories, const bool cell_vectorization_categories_strict, + const std::vector &parent_relation, std::vector & renumbering, std::vector & incompletely_filled_vectorization); diff --git a/source/matrix_free/task_info.cc b/source/matrix_free/task_info.cc index cbae4ea508..877395765d 100644 --- a/source/matrix_free/task_info.cc +++ b/source/matrix_free/task_info.cc @@ -695,16 +695,9 @@ namespace internal void - TaskInfo::collect_boundary_cells( - const unsigned int n_active_cells_in, - const unsigned int n_active_and_ghost_cells, - const unsigned int vectorization_length_in, + TaskInfo::make_boundary_cells_divisible( std::vector &boundary_cells) { - vectorization_length = vectorization_length_in; - n_active_cells = n_active_cells_in; - n_ghost_cells = n_active_and_ghost_cells - n_active_cells; - // try to make the number of boundary cells divisible by the number of // vectors in vectorization unsigned int fillup_needed = @@ -779,6 +772,7 @@ namespace internal const unsigned int dofs_per_cell, const std::vector &cell_vectorization_categories, const bool cell_vectorization_categories_strict, + const std::vector &parent_relation, std::vector & renumbering, std::vector & incompletely_filled_vectorization) { @@ -786,7 +780,6 @@ namespace internal (n_active_cells + vectorization_length - 1) / vectorization_length; const unsigned int n_ghost_slots = (n_ghost_cells + vectorization_length - 1) / vectorization_length; - const unsigned int n_boundary_cells = boundary_cells.size(); incompletely_filled_vectorization.resize(n_macro_cells + n_ghost_slots); renumbering.resize(n_active_cells + n_ghost_cells, @@ -799,28 +792,115 @@ namespace internal else partition_row_index.resize(5); - // Initially mark the cells according to the MPI ranking + int max_parent_index = -1; + for (unsigned int i : parent_relation) + if (i != numbers::invalid_unsigned_int) + max_parent_index = std::max(static_cast(i), max_parent_index); + unsigned int expected_group_size = + max_parent_index != -1 ? + std::count(parent_relation.begin(), parent_relation.end(), 0) : + 1; + std::vector cell_marked(n_active_cells + n_ghost_cells, 0); if (n_procs > 1) { - for (unsigned int i = 0; i < n_boundary_cells; ++i) - cell_marked[boundary_cells[i]] = 2; + // This lambda is used to mark the siblings (belong to the same + // parent) if a particular cell was touched in a pass as well + const auto mark_siblings = + [&](const unsigned int mark, + const std::vector &relevant_parents) { + for (unsigned int i = 0; i < n_active_cells; ++i) + if (cell_marked[i] == 0 && + parent_relation[i] != numbers::invalid_unsigned_int && + relevant_parents[parent_relation[i]]) + cell_marked[i] = mark; + }; + + // This lambda makes the cells at the processor boundary divisible + // by the vectorization length; we start the fillup with the more + // unstructured cells without a parent to increase chances that the + // cells sharing the parent get placed together + const auto fill_up_vectorization = [&](const unsigned int mark) { + unsigned int n_marked_cells = + std::count(cell_marked.begin(), + cell_marked.begin() + n_active_cells, + mark); + unsigned int n_available_cells = + std::count(cell_marked.begin(), + cell_marked.begin() + n_active_cells, + 0); + if (n_marked_cells % vectorization_length > 0 && + n_available_cells > 0) + { + unsigned int n_missing = + vectorization_length - + (n_marked_cells % vectorization_length); + for (unsigned int i = 0; i < n_active_cells; ++i) + if (cell_marked[i] == 0 && + parent_relation[i] == numbers::invalid_unsigned_int) + { + cell_marked[i] = mark; + --n_missing; + --n_available_cells; + ++n_marked_cells; + if (n_missing == 0) + break; + } + for (unsigned int i = 0; i < n_active_cells; ++i) + if (cell_marked[n_active_cells - 1 - i] == 0) + { + cell_marked[n_active_cells - 1 - i] = mark; + --n_missing; + --n_available_cells; + ++n_marked_cells; + if (n_missing == 0) + break; + } + } - Assert(boundary_cells.size() % vectorization_length == 0 || - boundary_cells.size() == n_active_cells, - ExcInternalError()); + Assert(n_marked_cells % vectorization_length == 0 || + n_available_cells == 0, + ExcInternalError("error " + std::to_string(n_marked_cells) + + " " + std::to_string(n_available_cells))); + return n_marked_cells; + }; + // Mark all cells needing data exchange as well as those belonging + // to the same parent with a special number + { + std::vector parent_at_boundary(max_parent_index + 1); + for (const unsigned int cell : boundary_cells) + { + cell_marked[cell] = 2; + if (parent_relation[cell] != numbers::invalid_unsigned_int) + parent_at_boundary[parent_relation[cell]] = true; + } + mark_siblings(2, parent_at_boundary); + } + const unsigned int n_boundary_cells = fill_up_vectorization(2); + + // Mark the cells that get placed before the cells at processor + // boundaries const unsigned int n_second_slot = ((n_active_cells - n_boundary_cells) / 2 / vectorization_length) * vectorization_length; - unsigned int count = 0; - unsigned int c = 0; - for (; c < n_active_cells && count < n_second_slot; ++c) - if (cell_marked[c] == 0) - { - cell_marked[c] = 1; - ++count; - } + unsigned int c = 0; + { + unsigned int count = 0; + std::vector parent_marked(max_parent_index + 1, false); + for (; c < n_active_cells && count < n_second_slot; ++c) + if (cell_marked[c] == 0) + { + if (parent_relation[c] != numbers::invalid_unsigned_int) + parent_marked[parent_relation[c]] = true; + cell_marked[c] = 1; + ++count; + } + mark_siblings(1, parent_marked); + fill_up_vectorization(1); + } + + // Finally, mark the remaining cells for (; c < n_active_cells; ++c) if (cell_marked[c] == 0) cell_marked[c] = 3; @@ -912,8 +992,78 @@ namespace internal // step 3: append cells according to categories for (unsigned int j = 0; j < n_categories; ++j) { - for (const unsigned int cell : renumbering_category[j]) - renumbering[counter++] = cell; + { + // Among the current category, we need to distinguish cells + // which we want to have grouped together and other cells. We + // first start by setting up these two categories + std::vector> + grouped_cells_tmp; + std::vector other_cells; + for (const unsigned int cell : renumbering_category[j]) + if (parent_relation[cell] == numbers::invalid_unsigned_int) + other_cells.push_back(cell); + else + grouped_cells_tmp.emplace_back(parent_relation[cell], cell); + + // Create a CRS data structure to identify each of the chunks + std::sort(grouped_cells_tmp.begin(), grouped_cells_tmp.end()); + std::vector crs_group(1); + for (unsigned int i = 1; i < grouped_cells_tmp.size(); ++i) + if (grouped_cells_tmp[i].first != + grouped_cells_tmp[i - 1].first) + crs_group.push_back(i); + crs_group.push_back(grouped_cells_tmp.size()); + + // Move groups that do not have the complete size (due to + // categories) to the 'other_cells' + std::vector grouped_cells; + for (unsigned int i = 0; i < crs_group.size() - 1; ++i) + if (crs_group[i + 1] - crs_group[i] < expected_group_size) + for (unsigned int j = crs_group[i]; j < crs_group[i + 1]; + ++j) + other_cells.push_back(grouped_cells_tmp[j].second); + else + for (unsigned int j = crs_group[i]; j < crs_group[i + 1]; + ++j) + grouped_cells.push_back(grouped_cells_tmp[j].second); + + // Sort the remaining cells + std::sort(other_cells.begin(), other_cells.end()); + + // Now fill in the cells from the two slots, the one with + // groups and the one without + auto regular = grouped_cells.begin(); + auto fillup = other_cells.begin(); + while (regular != grouped_cells.end() || + fillup != other_cells.end()) + { + // Case 1: Fill up until the next expected group size + while (counter % expected_group_size && + fillup != other_cells.end()) + renumbering[counter++] = *fillup++; + + // Case 2: If the start of the next group has a larger + // index than all indices we have queued from the + // irregular + if (fillup + expected_group_size <= other_cells.end() && + (regular == grouped_cells.end() || + *(fillup + expected_group_size - 1) < *regular)) + for (unsigned int j = 0; j < expected_group_size; ++j) + renumbering[counter++] = *fillup++; + + // Case 3: Add a group at once + if (regular != grouped_cells.end()) + for (unsigned int j = 0; j < expected_group_size; ++j) + renumbering[counter++] = *regular++; + + // Case 4: The groups are empty, so fill up from the other + // chunk + else + while (fillup != other_cells.end()) + renumbering[counter++] = *fillup++; + } + } + unsigned int remainder = renumbering_category[j].size() % vectorization_length; if (remainder)