From e29f37eaa0c193781d3eab0647499c9ce62d1e71 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 4 Oct 2024 12:32:50 +0200 Subject: [PATCH] SparsityTools::reorder_hierarchical: better documentation and implementation --- source/lac/sparsity_tools.cc | 235 ++++++++++++++++++++--------------- 1 file changed, 138 insertions(+), 97 deletions(-) diff --git a/source/lac/sparsity_tools.cc b/source/lac/sparsity_tools.cc index 0385f72abc..df8f841539 100644 --- a/source/lac/sparsity_tools.cc +++ b/source/lac/sparsity_tools.cc @@ -726,19 +726,36 @@ namespace SparsityTools ExcMessage( "Only valid for sparsity patterns which store all rows.")); - std::vector touched_nodes( - connectivity.n_rows(), numbers::invalid_dof_index); - std::vector row_lengths(connectivity.n_rows()); - std::set current_neighbors; - std::vector> groups; + // The algorithm below works by partitioning the rows in the + // connectivity graph, called nodes, into groups. The groups are defined + // as those nodes in immediate neighborhood of some pivot node, which we + // choose by minimal adjacency below. + + // We define two types of node categories for nodes not yet classified, + // one consisting of all nodes we've not seen at all, and one for nodes + // identified as neighbors (variable current_neighbors below) but not + // yet grouped. We use this classification in combination with an + // unsorted vector, which is much faster than keeping a sorted data + // structure (e.g. std::set) + constexpr types::global_dof_index unseen_node = + numbers::invalid_dof_index; + constexpr types::global_dof_index available_node = unseen_node - 1; + const types::global_dof_index n_nodes = connectivity.n_rows(); + std::vector touched_nodes(n_nodes, unseen_node); + + std::vector row_lengths(n_nodes); + std::vector current_neighbors; + std::vector group_starts(1); + std::vector group_indices; + group_indices.reserve(n_nodes); // First collect the number of neighbors for each node. We use this - // field to find next nodes with the minimum number of non-touched + // field to find the next node with the minimum number of non-touched // neighbors in the field n_remaining_neighbors, so we will count down // on this field. We also cache the row lengths because we need this // data frequently and getting it from the sparsity pattern is more // expensive. - for (types::global_dof_index row = 0; row < connectivity.n_rows(); ++row) + for (types::global_dof_index row = 0; row < n_nodes; ++row) { row_lengths[row] = connectivity.row_length(row); Assert(row_lengths[row] > 0, ExcInternalError()); @@ -749,139 +766,163 @@ namespace SparsityTools // graph is not connected while (true) { - // Find cell with the minimal number of neighbors (typically a - // corner node when based on FEM meshes). If no cell is left, we are + // Find node with the minimal number of neighbors (typically a + // corner node when based on FEM meshes). If no node is left, we are // done. Together with the outer while loop, this loop can possibly // be of quadratic complexity in the number of disconnected - // partitions, i.e. up to connectivity.n_rows() in the worst case, + // partitions, i.e. up to n_nodes in the worst case, // but that is not the usual use case of this loop and thus not // optimized for. - std::pair - min_neighbors(numbers::invalid_dof_index, - numbers::invalid_dof_index); - for (types::global_dof_index i = 0; i < touched_nodes.size(); ++i) - if (touched_nodes[i] == numbers::invalid_dof_index) - if (row_lengths[i] < min_neighbors.second) - { - min_neighbors = std::make_pair(i, n_remaining_neighbors[i]); - if (n_remaining_neighbors[i] <= 1) - break; - } - if (min_neighbors.first == numbers::invalid_dof_index) - break; + { + unsigned int candidate_valence = numbers::invalid_unsigned_int; + types::global_dof_index candidate_index = + numbers::invalid_dof_index; + for (types::global_dof_index i = 0; i < n_nodes; ++i) + if (touched_nodes[i] == unseen_node) + if (row_lengths[i] < candidate_valence) + { + candidate_index = i; + candidate_valence = n_remaining_neighbors[i]; + if (candidate_valence <= 1) + break; + } + if (candidate_index == numbers::invalid_dof_index) + break; - Assert(min_neighbors.second > 0, ExcInternalError()); + Assert(candidate_valence > 0, ExcInternalError()); + + current_neighbors = {candidate_index}; + touched_nodes[candidate_index] = available_node; + } - current_neighbors.clear(); - current_neighbors.insert(min_neighbors.first); - while (!current_neighbors.empty()) + while (true) { - // Find node with minimum number of untouched neighbors among the - // next set of possible neighbors - min_neighbors = std::make_pair(numbers::invalid_dof_index, - numbers::invalid_dof_index); - for (const auto current_neighbor : current_neighbors) + // Find node with minimum number of untouched neighbors among + // the next set of possible neighbors (= valence), and among the + // set of nodes with the minimal number of neighbors, choose the + // one with the largest number of touched neighbors (i.e., the + // largest row length). + // + // This loop is typically the most expensive part for large + // graphs and thus only run once. We also do some cleanup, i.e., + // the indices added to a group in the previous round need to be + // removed at this point. + unsigned int candidate_valence = numbers::invalid_unsigned_int; + types::global_dof_index candidate_index = + numbers::invalid_dof_index; + unsigned int candidate_row_length = 0; + const unsigned int loop_length = current_neighbors.size(); + unsigned int write_index = 0; + for (unsigned int i = 0; i < loop_length; ++i) { - Assert(touched_nodes[current_neighbor] == - numbers::invalid_dof_index, + const types::global_dof_index node = current_neighbors[i]; + Assert(touched_nodes[node] != unseen_node, ExcInternalError()); - if (n_remaining_neighbors[current_neighbor] < - min_neighbors.second) - min_neighbors = - std::make_pair(current_neighbor, - n_remaining_neighbors[current_neighbor]); + if (touched_nodes[node] == available_node) + { + current_neighbors[write_index] = node; + ++write_index; + if (n_remaining_neighbors[node] < candidate_valence || + (n_remaining_neighbors[node] == candidate_valence && + row_lengths[node] > candidate_row_length)) + { + candidate_index = node; + candidate_valence = n_remaining_neighbors[node]; + candidate_row_length = row_lengths[node]; + } + } } + current_neighbors.resize(write_index); + for (const types::global_dof_index node : current_neighbors) + Assert(touched_nodes[node] == available_node, + ExcInternalError()); - // Among the set of nodes with the minimal number of neighbors, - // choose the one with the largest number of touched neighbors, - // i.e., the one with the largest row length - const types::global_dof_index best_row_length = - min_neighbors.second; - for (const auto current_neighbor : current_neighbors) - if (n_remaining_neighbors[current_neighbor] == best_row_length) - if (row_lengths[current_neighbor] > min_neighbors.second) - min_neighbors = - std::make_pair(current_neighbor, - row_lengths[current_neighbor]); + // No more neighbors left -> terminate loop + if (current_neighbors.empty()) + break; // Add the pivot and all direct neighbors of the pivot node not // yet touched to the list of new entries. - groups.emplace_back(); - std::vector &next_group = groups.back(); - - next_group.push_back(min_neighbors.first); - touched_nodes[min_neighbors.first] = groups.size() - 1; - for (DynamicSparsityPattern::iterator it = - connectivity.begin(min_neighbors.first); - it != connectivity.end(min_neighbors.first); + group_indices.push_back(candidate_index); + touched_nodes[candidate_index] = group_starts.size() - 1; + const auto end_it = connectivity.end(candidate_index); + for (auto it = connectivity.begin(candidate_index); it != end_it; ++it) - if (touched_nodes[it->column()] == numbers::invalid_dof_index) + if (touched_nodes[it->column()] >= available_node) { - next_group.push_back(it->column()); - touched_nodes[it->column()] = groups.size() - 1; + group_indices.push_back(it->column()); + touched_nodes[it->column()] = group_starts.size() - 1; } - - // Add all neighbors of the current list not yet touched to the - // set of possible next pivots. The added node is no longer a - // valid neighbor (here we assume symmetry of the - // connectivity). Delete the entries of the current list from - // the set of possible next pivots. - for (const auto index : next_group) + group_starts.push_back(group_indices.size()); + + // Add all neighbors of the current list not yet seen to the set + // of possible next nodes. The added node is grouped and thus no + // longer a valid neighbor (here we assume symmetry of the + // connectivity). It will be removed from the list of neighbors + // by the code further up in the next iteration of the + // surrounding loop. + for (types::global_dof_index index = + group_starts[group_starts.size() - 2]; + index < group_starts.back(); + ++index) { - for (DynamicSparsityPattern::iterator it = - connectivity.begin(index); - it != connectivity.end(index); - ++it) + auto it = connectivity.begin(group_indices[index]); + const auto end_row = connectivity.end(group_indices[index]); + for (; it != end_row; ++it) { - if (touched_nodes[it->column()] == - numbers::invalid_dof_index) - current_neighbors.insert(it->column()); + if (touched_nodes[it->column()] == unseen_node) + { + current_neighbors.push_back(it->column()); + touched_nodes[it->column()] = available_node; + } n_remaining_neighbors[it->column()]--; } - current_neighbors.erase(index); } } } // Sanity check: for all nodes, there should not be any neighbors left - for (types::global_dof_index row = 0; row < connectivity.n_rows(); ++row) + for (types::global_dof_index row = 0; row < n_nodes; ++row) Assert(n_remaining_neighbors[row] == 0, ExcInternalError()); // If the number of groups is smaller than the number of nodes, we // continue by recursively calling this method - if (groups.size() < connectivity.n_rows()) + const unsigned int n_groups = group_starts.size() - 1; + if (n_groups < n_nodes) { // Form the connectivity of the groups - DynamicSparsityPattern connectivity_next(groups.size(), - groups.size()); - for (types::global_dof_index i = 0; i < groups.size(); ++i) - for (types::global_dof_index col = 0; col < groups[i].size(); ++col) - for (DynamicSparsityPattern::iterator it = - connectivity.begin(groups[i][col]); - it != connectivity.end(groups[i][col]); - ++it) - connectivity_next.add(i, touched_nodes[it->column()]); + DynamicSparsityPattern connectivity_next(n_groups, n_groups); + for (types::global_dof_index row = 0; row < n_groups; ++row) + for (types::global_dof_index index = group_starts[row]; + index < group_starts[row + 1]; + ++index) + { + auto it = connectivity.begin(group_indices[index]); + const auto end_it = connectivity.end(group_indices[index]); + for (; it != end_it; ++it) + connectivity_next.add(row, touched_nodes[it->column()]); + } // Recursively call the reordering - std::vector renumbering_next(groups.size()); + std::vector renumbering_next(n_groups); reorder_hierarchical(connectivity_next, renumbering_next); // Renumber the indices group by group according to the incoming // ordering for the groups - for (types::global_dof_index i = 0, count = 0; i < groups.size(); ++i) - for (types::global_dof_index col = 0; - col < groups[renumbering_next[i]].size(); - ++col, ++count) - renumbering[count] = groups[renumbering_next[i]][col]; + for (types::global_dof_index row = 0, c = 0; row < n_groups; ++row) + for (types::global_dof_index index = + group_starts[renumbering_next[row]]; + index < group_starts[renumbering_next[row] + 1]; + ++index, ++c) + renumbering[c] = group_indices[index]; } else { // All groups should have size one and no more recursion is possible, // so use the numbering of the groups - for (types::global_dof_index i = 0, count = 0; i < groups.size(); ++i) - for (types::global_dof_index col = 0; col < groups[i].size(); - ++col, ++count) - renumbering[count] = groups[i][col]; + unsigned int c = 0; + for (const types::global_dof_index i : group_indices) + renumbering[c++] = i; } } } // namespace internal -- 2.39.5