From 7ddeeb16a676c699e5c50025a77f8b3efc88bd0b Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 4 May 2020 13:17:48 +0200 Subject: [PATCH] Implement better strategy for data locality detection --- .../matrix_free/matrix_free.templates.h | 4 +- include/deal.II/matrix_free/task_info.h | 15 +- source/matrix_free/task_info.cc | 500 ++++++++---------- tests/matrix_free/thread_correctness.cc | 12 +- 4 files changed, 233 insertions(+), 298 deletions(-) diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index f26f18b63b..faaac4e085 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -463,7 +463,7 @@ MatrixFree::internal_reinit( 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); + dummy, 1, false, dummy, false, dummy, dummy, dummy2); for (unsigned int i = 0; i < dof_info.size(); ++i) { @@ -1058,6 +1058,8 @@ MatrixFree::initialize_indices( } task_info.create_blocks_serial(subdomain_boundary_cells, dofs_per_cell, + dof_handlers.active_dof_handler == + DoFHandlers::hp, dof_info[0].cell_active_fe_index, strict_categories, parent_relation, diff --git a/include/deal.II/matrix_free/task_info.h b/include/deal.II/matrix_free/task_info.h index ed9a128b02..61299b88a3 100644 --- a/include/deal.II/matrix_free/task_info.h +++ b/include/deal.II/matrix_free/task_info.h @@ -145,14 +145,20 @@ namespace internal * Sets up the blocks for running the cell loop based on the options * controlled by the input arguments. * - * @param boundary_cells A list of cells that need to exchange data prior - * to performing computations. These will be given a certain id in the - * partitioning. + * @param cells_with_comm A list of cells that need to exchange data + * prior to performing computations. These will be given a certain id in + * the partitioning to make sure cell loops that overlap communication + * with communication have the ghost data ready. * * @param dofs_per_cell Gives an expected value for the number of degrees * of freedom on a cell, which is used to determine the block size for * interleaving cell and face integrals. * + * @param categories_are_hp Defines whether + * `cell_vectorization_categories` is originating from a hp adaptive + * computation with variable polynomial degree or a user-defined + * variant. + * * @param cell_vectorization_categories This set of categories defines * the cells that should be grouped together inside the lanes of a * vectorized array. This can be the polynomial degree in an hp-element @@ -179,8 +185,9 @@ namespace internal */ void create_blocks_serial( - const std::vector &boundary_cells, + const std::vector &cells_with_comm, const unsigned int dofs_per_cell, + const bool categories_are_hp, const std::vector &cell_vectorization_categories, const bool cell_vectorization_categories_strict, const std::vector &parent_relation, diff --git a/source/matrix_free/task_info.cc b/source/matrix_free/task_info.cc index 877395765d..953041fc53 100644 --- a/source/matrix_free/task_info.cc +++ b/source/matrix_free/task_info.cc @@ -768,175 +768,42 @@ namespace internal void TaskInfo::create_blocks_serial( - const std::vector &boundary_cells, + const std::vector &cells_with_comm, const unsigned int dofs_per_cell, + const bool categories_are_hp, 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) { - const unsigned int n_macro_cells = - (n_active_cells + vectorization_length - 1) / vectorization_length; - const unsigned int n_ghost_slots = - (n_ghost_cells + vectorization_length - 1) / vectorization_length; - - incompletely_filled_vectorization.resize(n_macro_cells + n_ghost_slots); - renumbering.resize(n_active_cells + n_ghost_cells, - numbers::invalid_unsigned_int); - - // Define the outer number of partitions. In the MPI case, we have three - // partitions (part before comm, part with comm, part after comm) - if (n_procs == 1) - partition_row_index.resize(3); - else - partition_row_index.resize(5); - - 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) - { - // 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(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 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; - for (; c < n_active_cells + n_ghost_cells; ++c) - if (cell_marked[c] == 0) - cell_marked[c] = 4; - } - else - std::fill(cell_marked.begin(), cell_marked.end(), 1); - - for (const unsigned char marker : cell_marked) - { - (void)marker; - Assert(marker != 0, ExcInternalError()); - } - + // Give the compiler a chance to detect that vectorization_length is a + // power of two, which allows it to replace integer divisions by shifts + unsigned int vectorization_length_bits = 0; + unsigned int my_length = vectorization_length; + while (my_length >>= 1) + ++vectorization_length_bits; + const unsigned int n_lanes = 1 << vectorization_length_bits; + + // Step 1: find tight map of categories for not taking exceeding amounts + // of memory below. Sort the new categories by the numbers in the + // old one to ensure we respect the given rules unsigned int n_categories = 1; - std::vector tight_category_map; + std::vector tight_category_map(n_active_cells, 0); if (cell_vectorization_categories.empty() == false) { AssertDimension(cell_vectorization_categories.size(), n_active_cells + n_ghost_cells); - // create a tight map of categories for not taking exceeding amounts - // of memory below. Sort the new categories by the numbers in the - // old one. - tight_category_map.resize(n_active_cells + n_ghost_cells); std::set used_categories; - for (unsigned int i = 0; i < n_active_cells + n_ghost_cells; ++i) + for (unsigned int i = 0; i < n_active_cells; ++i) used_categories.insert(cell_vectorization_categories[i]); std::vector used_categories_vector( used_categories.size()); n_categories = 0; for (const auto &it : used_categories) used_categories_vector[n_categories++] = it; - for (unsigned int i = 0; i < n_active_cells + n_ghost_cells; ++i) + for (unsigned int i = 0; i < n_active_cells; ++i) { const unsigned int index = std::lower_bound(used_categories_vector.begin(), @@ -946,163 +813,222 @@ namespace internal AssertIndexRange(index, used_categories_vector.size()); tight_category_map[i] = index; } - - // leave some more space for empty lanes - incompletely_filled_vectorization.resize( - incompletely_filled_vectorization.size() + 4 * n_categories); } - else - tight_category_map.resize(n_active_cells + n_ghost_cells, 0); - cell_partition_data.clear(); - cell_partition_data.resize(1, 0); - unsigned int counter = 0; - unsigned int n_cells = 0; + // Step 2: Sort the cells by the category. If we want to fill up the + // ranges in vectorization, promote some of the cells to a higher + // category std::vector> renumbering_category(n_categories); - for (unsigned int block = 1; block < (n_procs > 1u ? 5u : 3u); ++block) - { - // step 1: sort by category - for (unsigned int i = 0; i < n_active_cells + n_ghost_cells; ++i) - if (cell_marked[i] == block) - renumbering_category[tight_category_map[i]].push_back(i); - - // step 2: if we want to fill up the ranges in vectorization, promote - // some of the cells to a higher category - if (cell_vectorization_categories_strict == false && n_categories > 1) - for (unsigned int j = n_categories - 1; j > 0; --j) + for (unsigned int i = 0; i < n_active_cells; ++i) + renumbering_category[tight_category_map[i]].push_back(i); + + if (cell_vectorization_categories_strict == false && n_categories > 1) + for (unsigned int j = n_categories - 1; j > 0; --j) + { + unsigned int lower_index = j - 1; + while (renumbering_category[j].size() % n_lanes) { - unsigned int lower_index = j - 1; - while (renumbering_category[j].size() % vectorization_length) + while (renumbering_category[j].size() % n_lanes && + !renumbering_category[lower_index].empty()) { - while (renumbering_category[j].size() % - vectorization_length && - !renumbering_category[lower_index].empty()) - { - renumbering_category[j].push_back( - renumbering_category[lower_index].back()); - renumbering_category[lower_index].pop_back(); - } - if (lower_index == 0) - break; - else - --lower_index; + renumbering_category[j].push_back( + renumbering_category[lower_index].back()); + renumbering_category[lower_index].pop_back(); } + if (lower_index == 0) + break; + else + --lower_index; } + } - // step 3: append cells according to categories - for (unsigned int j = 0; j < n_categories; ++j) - { + // Step 3: Use the parent relation to find a good grouping of cells. To + // do this, we first put cells of each category defined above into two + // bins, those which we know can be grouped together by the given parent + // relation and those which cannot + std::vector temporary_numbering; + temporary_numbering.reserve(n_active_cells + + (n_lanes - 1) * n_categories); + const unsigned int n_cells_per_parent = + std::count(parent_relation.begin(), parent_relation.end(), 0); + std::vector category_size; + for (unsigned int j = 0; j < n_categories; ++j) + { + std::vector> grouped_cells; + std::vector other_cells; + for (const unsigned int cell : renumbering_category[j]) + if (parent_relation.empty() || + parent_relation[cell] == numbers::invalid_unsigned_int) + other_cells.push_back(cell); + else + grouped_cells.emplace_back(parent_relation[cell], cell); + + // Compute the number of cells per group + std::sort(grouped_cells.begin(), grouped_cells.end()); + std::vector n_cells_per_group; + unsigned int length = 0; + for (unsigned int i = 0; i < grouped_cells.size(); ++i, ++length) + if (i > 0 && grouped_cells[i].first != grouped_cells[i - 1].first) { - // 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); + n_cells_per_group.push_back(length); + length = 0; + } + if (length > 0) + n_cells_per_group.push_back(length); + + // Move groups that do not have the complete size (due to + // categories) to the 'other_cells'. The cells with correct group + // size are immediately appended to the temporary cell numbering + auto group_it = grouped_cells.begin(); + for (unsigned int length : n_cells_per_group) + if (length < n_cells_per_parent) + for (unsigned int j = 0; j < length; ++j) + other_cells.push_back((group_it++)->second); + else + { + // we should not have more cells in a group than in the first + // check we did above + AssertDimension(length, n_cells_per_parent); + for (unsigned int j = 0; j < length; ++j) + temporary_numbering.push_back((group_it++)->second); + } - // Sort the remaining cells - std::sort(other_cells.begin(), other_cells.end()); + // Sort the remaining cells and append them as well + std::sort(other_cells.begin(), other_cells.end()); + temporary_numbering.insert(temporary_numbering.end(), + 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++; - } - } + while (temporary_numbering.size() % n_lanes != 0) + temporary_numbering.push_back(numbers::invalid_unsigned_int); - unsigned int remainder = - renumbering_category[j].size() % vectorization_length; - if (remainder) - incompletely_filled_vectorization - [renumbering_category[j].size() / vectorization_length + - n_cells] = remainder; - const unsigned int n_my_macro_cells = - (renumbering_category[j].size() + vectorization_length - 1) / - vectorization_length; - renumbering_category[j].clear(); - - // step 4: create blocks for face integrals, make the number of - // cells divisible by 4 if possible - const unsigned int block_size = - std::max((2048U / dofs_per_cell) / 8 * 4, 2U); - if (block < 4) - for (unsigned int k = 0; k < n_my_macro_cells; k += block_size) - cell_partition_data.push_back( - n_cells + std::min(k + block_size, n_my_macro_cells)); - else - cell_partition_data.back() += n_my_macro_cells; - n_cells += n_my_macro_cells; - } - partition_row_index[block] = cell_partition_data.size() - 1; - if (block == 3 || (block == 1 && n_procs == 1)) - cell_partition_data.push_back(n_cells); + category_size.push_back(temporary_numbering.size()); + } + + // Step 4: Identify the batches with cells marked as "comm" + std::vector batch_with_comm(temporary_numbering.size() / n_lanes, + false); + std::vector temporary_numbering_inverse(n_active_cells); + for (unsigned int i = 0; i < temporary_numbering.size(); ++i) + if (temporary_numbering[i] != numbers::invalid_unsigned_int) + temporary_numbering_inverse[temporary_numbering[i]] = i; + for (const unsigned int cell : cells_with_comm) + batch_with_comm[temporary_numbering_inverse[cell] / n_lanes] = true; + + // Step 5: Sort the batches of cells by their last cell index to get + // good locality, assuming that the initial cell order is of good + // locality. In case we have hp calculations with categories, we need to + // sort also by the category. + std::vector> batch_order; + std::vector> batch_order_comm; + for (unsigned int i = 0; i < temporary_numbering.size(); i += n_lanes) + { + unsigned int max_index = 0; + for (unsigned int j = 0; j < n_lanes; ++j) + if (temporary_numbering[i + j] < numbers::invalid_unsigned_int) + max_index = std::max(temporary_numbering[i + j], max_index); + const unsigned int category_hp = + categories_are_hp ? + std::upper_bound(category_size.begin(), category_size.end(), i) - + category_size.begin() : + 0; + const std::array next{category_hp, max_index, i}; + if (batch_with_comm[i / n_lanes]) + batch_order_comm.emplace_back(next); + else + batch_order.emplace_back(next); } - if (cell_vectorization_categories_strict == true) + + std::sort(batch_order.begin(), batch_order.end()); + std::sort(batch_order_comm.begin(), batch_order_comm.end()); + + // Step 6: Put the cells with communication in the middle of the cell + // range. For the MPI case, we need three groups to enable overlap for + // communication and computation (part before comm, part with comm, part + // after comm), whereas we need one for the other case. And in each + // case, we allow for a slot of "ghosted" cells. + std::vector blocks; + if (n_procs == 1) { - Assert(n_cells >= n_macro_cells + n_ghost_slots, ExcInternalError()); + if (batch_order.empty()) + std::swap(batch_order_comm, batch_order); + Assert(batch_order_comm.empty(), ExcInternalError()); + partition_row_index.resize(3); + blocks = {0, static_cast(batch_order.size())}; } else { - AssertDimension(n_cells, n_macro_cells + n_ghost_slots); + partition_row_index.resize(5); + const unsigned int comm_begin = batch_order.size() / 2; + batch_order.insert(batch_order.begin() + comm_begin, + batch_order_comm.begin(), + batch_order_comm.end()); + const unsigned int comm_end = comm_begin + batch_order_comm.size(); + const unsigned int end = batch_order.size(); + blocks = {0, comm_begin, comm_end, end}; } - AssertDimension(cell_partition_data.back(), n_cells); - AssertDimension(counter, n_active_cells + n_ghost_cells); - incompletely_filled_vectorization.resize(cell_partition_data.back()); + // Step 7: Fill in the data by batches for the locally owned cells. + const unsigned int n_cell_batches = batch_order.size(); + const unsigned int n_ghost_batches = + (n_ghost_cells + n_lanes - 1) / n_lanes; + incompletely_filled_vectorization.resize(n_cell_batches + + n_ghost_batches); + + cell_partition_data.clear(); + cell_partition_data.resize(1, 0); + + renumbering.clear(); + renumbering.resize(n_active_cells + n_ghost_cells, + numbers::invalid_unsigned_int); + + unsigned int counter = 0; + for (unsigned int block = 0; block < blocks.size() - 1; ++block) + { + const unsigned int grain_size = + std::max((2048U / dofs_per_cell) / 8 * 4, 2U); + for (unsigned int k = blocks[block]; k < blocks[block + 1]; + k += grain_size) + cell_partition_data.push_back( + std::min(k + grain_size, blocks[block + 1])); + partition_row_index[block + 1] = cell_partition_data.size() - 1; + + // Set the numbering according to the reordered temporary one + for (unsigned int k = blocks[block]; k < blocks[block + 1]; ++k) + { + const unsigned int pos = batch_order[k][2]; + unsigned int j = 0; + for (; j < n_lanes && temporary_numbering[pos + j] != + numbers::invalid_unsigned_int; + ++j) + renumbering[counter++] = temporary_numbering[pos + j]; + if (j < n_lanes) + incompletely_filled_vectorization[k] = j; + } + } + AssertDimension(counter, n_active_cells); + + // Step 8: Treat the ghost cells + for (unsigned int cell = n_active_cells; + cell < n_active_cells + n_ghost_cells; + ++cell) + { + if (!cell_vectorization_categories.empty()) + AssertDimension(cell_vectorization_categories[cell], + cell_vectorization_categories[n_active_cells]); + renumbering[cell] = cell; + } + if (n_ghost_cells % n_lanes) + incompletely_filled_vectorization.back() = n_ghost_cells % n_lanes; + cell_partition_data.push_back(n_cell_batches + n_ghost_batches); + partition_row_index.back() = cell_partition_data.size() - 1; + +#ifdef DEBUG + std::vector renumber_cpy(renumbering); + std::sort(renumber_cpy.begin(), renumber_cpy.end()); + for (unsigned int i = 0; i < renumber_cpy.size(); ++i) + AssertDimension(i, renumber_cpy[i]); +#endif } diff --git a/tests/matrix_free/thread_correctness.cc b/tests/matrix_free/thread_correctness.cc index 4d9f6681d6..d8519f72e9 100644 --- a/tests/matrix_free/thread_correctness.cc +++ b/tests/matrix_free/thread_correctness.cc @@ -128,6 +128,8 @@ sub_test() // make 10 sweeps in order to get in some // variation to the threaded program + const double float_factor = + std::is_same::value ? 0.01 : 1.; for (unsigned int sweep = 0; sweep < 10; ++sweep) { mf_color.vmult(out_color, in_dist); @@ -135,14 +137,12 @@ sub_test() out_color -= out_dist; double diff_norm = out_color.linfty_norm(); - deallog << "Sweep " << sweep - << ", error in partition/color: " << diff_norm - << std::endl; + deallog << "Sweep " << sweep << ", error in partition/color: " + << diff_norm * float_factor << std::endl; out_partition -= out_dist; diff_norm = out_partition.linfty_norm(); - deallog << "Sweep " << sweep - << ", error in partition/partition: " << diff_norm - << std::endl; + deallog << "Sweep " << sweep << ", error in partition/partition: " + << diff_norm * float_factor << std::endl; } deallog << std::endl; } -- 2.39.5