Major performance improvement for re-ordering of large graphs

cpp/doc/source/graph.rst

@@ -22,6 +22,9 @@ Re-ordering
 .. doxygenfunction:: dolfinx::graph::reorder_gps
    :project: DOLFINx
 
+.. doxygenfunction:: dolfinx::graph::reorder_rcm
+   :project: DOLFINx
+
 
 Partitioning
 ------------

cpp/dolfinx/fem/DofMap.cpp

@@ -212,7 +212,7 @@ std::pair<DofMap, std::vector<std::int32_t>> DofMap::collapse(
   if (!reorder_fn)
   {
     reorder_fn = [](const graph::AdjacencyList<std::int32_t>& g)
-    { return graph::reorder_gps(g); };
+    { return graph::reorder_rcm(g); };
   }
   // Create new dofmap
   auto create_subdofmap = [](MPI_Comm comm, auto index_map_bs, auto& layout,

cpp/dolfinx/graph/ordering.cpp

@@ -1,4 +1,4 @@
-// Copyright (C) 2021 Chris Richardson
+// Copyright (C) 2021-2026 Chris Richardson and Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -13,6 +13,7 @@
 #include <dolfinx/common/log.h>
 #include <limits>
 #include <span>
+#include <thread>
 
 using namespace dolfinx;
 
@@ -89,9 +90,13 @@
 graph::AdjacencyList<int>
 create_level_structure(const graph::AdjacencyList<int>& graph, int s)
 {
-  common::Timer t("GPS: create_level_structure");
+  common::Timer t("Graph: create_level_structure");
 
-  // Note: int8 is often faster than bool
+  // Note: int8 is often faster than bool. A fresh buffer is
+  // allocated on each call (rather than a reused scratch buffer)
+  // since the allocator's bulk zero-fill is faster than resetting
+  // only the touched entries when a large fraction of the graph is
+  // typically visited, as is the case for compact/dense mesh graphs.
   std::vector<std::int8_t> labelled(graph.num_nodes(), false);
   labelled[s] = true;
 
@@ -128,7 +133,8 @@
 // with -1 in the vector rlabel).
 std::vector<std::int32_t>
 gps_reorder_unlabelled(const graph::AdjacencyList<std::int32_t>& graph,
-                       std::span<const std::int32_t> rlabel)
+                       std::span<const std::int32_t> rlabel,
+                       std::size_t max_candidates, std::size_t num_threads)
 {
   common::Timer timer("Gibbs-Poole-Stockmeyer ordering");
 
@@ -160,24 +166,65 @@
   std::vector<int> S;
   while (!done)
   {
-    // Sort final level S of Lv into increasing degree order
+    // Sort final level S of Lv into increasing degree order, capped to
+    // at most `max_candidates` (most likely peripheral) candidates.
+    // Testing every vertex in the final level costs O(|S| * (V+E)) per
+    // round; for a compact 3D mesh |S| scales as the cross-sectional
+    // area (~ N^(2/3)), which dominates run time at scale. A handful
+    // of the lowest-degree candidates often captures the same
+    // pseudo-diameter and width-minimising choice as an exhaustive
+    // search, but not always -- non-convex domains and strongly graded
+    // meshes can lose some bandwidth/profile quality with a small cap.
     auto lv_final = lv.links(lv.num_nodes() - 1);
-    S.resize(lv_final.size());
+    S.resize(std::min(lv_final.size(), max_candidates));
     std::partial_sort_copy(lv_final.begin(), lv_final.end(), S.begin(), S.end(),
                            cmp_degree);
     int w_min = std::numeric_limits<int>::max();
     done = true;
 
     // C. Generate level structures rooted at vertices s in S selected
-    // in order of increasing degree.
-    for (int s : S)
+    // in order of increasing degree. Each candidate's level structure
+    // is independent of the others, so when num_threads > 1 they are
+    // computed concurrently. The decision of which candidate to use
+    // is still made afterwards, sequentially, in the original order,
+    // so the result is identical regardless of num_threads.
+    std::vector<graph::AdjacencyList<int>> lstmps(S.size(),
+                                                  graph::AdjacencyList<int>(0));
+    std::size_t nt = std::min(num_threads, S.size());
+    if (nt <= 1)
+    {
+      for (std::size_t i = 0; i < S.size(); ++i)
+        lstmps[i] = create_level_structure(graph, S[i]);
+    }
+    else
     {
-      graph::AdjacencyList<int> lstmp = create_level_structure(graph, s);
+      std::vector<std::jthread> workers;
+      workers.reserve(nt);
+      std::size_t chunk = (S.size() + nt - 1) / nt;
+      for (std::size_t t = 0; t < nt; ++t)
+      {
+        std::size_t begin = t * chunk;
+        std::size_t end = std::min(S.size(), begin + chunk);
+        if (begin >= end)
+          break;
+        workers.emplace_back(
+            [&S, &lstmps, &graph, begin, end]()
+            {
+              for (std::size_t i = begin; i < end; ++i)
+                lstmps[i] = create_level_structure(graph, S[i]);
+            });
+      }
+    }
+
+    for (std::size_t i = 0; i < S.size(); ++i)
+    {
+      int s = S[i];
+      graph::AdjacencyList<int>& lstmp = lstmps[i];
       if (lstmp.num_nodes() > lv.num_nodes())
       {
         // Found a deeper level structure, so restart
         v = s;
-        lv = lstmp;
+        lv = std::move(lstmp);
         done = false;
         break;
       }
@@ -188,7 +235,7 @@
       {
         w_min = w;
         u = s;
-        lu = lstmp;
+        lu = std::move(lstmp);
       }
     }
   }
@@ -284,16 +331,18 @@
   rv.push_back(v);
   labelled[v] = true;
 
-  // Temporary work vectors
-  std::vector<std::int8_t> in_level;
+  // Temporary work vectors. `in_level` is sized once and kept all-false
+  // between iterations by undoing exactly the entries it set, rather
+  // than re-zeroing all n entries on every level (which is O(n * k)
+  // over the whole level structure instead of O(n)).
+  std::vector<std::int8_t> in_level(n, false);
   std::vector<int> rv_next;
   std::vector<int> nbr, nbr_next;
   std::vector<int> nrem;
 
   for (const std::vector<int>& lslevel : ls)
   {
     // Mark all nodes of the current level
-    in_level.assign(n, false);
     for (int w : lslevel)
       in_level[w] = true;
 
@@ -351,16 +400,131 @@
 
     // Insert already-labelled nodes of next level
     rv.insert(rv.end(), rv_next.begin(), rv_next.end());
+
+    // Undo the marks set for this level so in_level is all-false again
+    for (int w : lslevel)
+      in_level[w] = false;
   }
 
   return rv;
 }
+//-----------------------------------------------------------------------------
+// Reverse Cuthill-McKee algorithm, finding a reordering for the given
+// graph, operating only on nodes which are yet unlabelled (indicated
+// with -1 in the vector rlabel).
+std::vector<std::int32_t>
+rcm_reorder_unlabelled(const graph::AdjacencyList<std::int32_t>& graph,
+                       std::span<const std::int32_t> rlabel)
+{
+  common::Timer timer("Reverse Cuthill-McKee ordering");
+
+  const std::int32_t n = graph.num_nodes();
+
+  // Degree comparison function
+  auto cmp_degree = [&graph](auto a, auto b)
+  { return graph.num_links(a) < graph.num_links(b); };
+
+  // Pick an arbitrary vertex of minimal degree and call it v
+  std::int32_t v = 0;
+  std::int32_t dmin = std::numeric_limits<std::int32_t>::max();
+  for (std::int32_t i = 0; i < n; ++i)
+  {
+    if (int d = graph.num_links(i); rlabel[i] == -1 and d < dmin)
+    {
+      v = i;
+      dmin = d;
+    }
+  }
+
+  // Find a pseudo-peripheral root: repeatedly move to the minimum-degree
+  // vertex of the deepest level found so far (the classical George-Liu
+  // "double sweep"), stopping once no deeper level structure is found.
+  // Unlike Gibbs-Poole-Stockmeyer, only the single lowest-degree
+  // candidate is tried at each step -- Cuthill-McKee has no
+  // width-minimising second phase to feed a wider candidate set into,
+  // so testing more candidates here would only add cost, not quality.
+  graph::AdjacencyList<int> lv = create_level_structure(graph, v);
+  bool done = false;
+  while (!done)
+  {
+    auto lv_final = lv.links(lv.num_nodes() - 1);
+    int s = *std::ranges::min_element(lv_final, cmp_degree);
+    graph::AdjacencyList<int> lstmp = create_level_structure(graph, s);
+    if (lstmp.num_nodes() > lv.num_nodes())
+    {
+      v = s;
+      lv = std::move(lstmp);
+    }
+    else
+      done = true;
+  }
+
+  // Cuthill-McKee numbering, breadth-first from the root: each time a
+  // vertex is processed, its not-yet-discovered neighbours are appended
+  // in increasing degree order. Note this sorts each vertex's own
+  // newly-discovered neighbours as they are found, not each BFS level
+  // as a whole group -- the latter is a common but lower-quality
+  // simplification, since it discards the parent/discovery order that
+  // the standard algorithm relies on.
+  std::vector<std::int8_t> labelled(n, false);
+  std::vector<int> rv;
+  rv.reserve(n);
+  rv.push_back(v);
+  labelled[v] = true;
+
+  std::vector<int> nbr;
+  for (std::size_t current = 0; current < rv.size(); ++current)
+  {
+    nbr.clear();
+    for (int w : graph.links(rv[current]))
+    {
+      if (!labelled[w])
+      {
+        nbr.push_back(w);
+        labelled[w] = true;
+      }
+    }
+    std::ranges::sort(nbr, cmp_degree);
+    rv.insert(rv.end(), nbr.begin(), nbr.end());
+  }
+
+  // Reverse the numbering -- the "reverse" in Reverse Cuthill-McKee --
+  // which tends to reduce profile relative to plain Cuthill-McKee.
+  std::ranges::reverse(rv);
+
+  return rv;
+}
 
 } // namespace
 
 //-----------------------------------------------------------------------------
 std::vector<std::int32_t>
-graph::reorder_gps(const graph::AdjacencyList<std::int32_t>& graph)
+graph::reorder_gps(const graph::AdjacencyList<std::int32_t>& graph,
+                   std::size_t max_candidates, std::size_t num_threads)
+{
+  const std::int32_t n = graph.num_nodes();
+  std::vector<std::int32_t> r(n, -1);
+  std::vector<std::int32_t> rv;
+
+  // Repeat for each disconnected part of the graph
+  int count = 0;
+  while (count < n)
+  {
+    rv = gps_reorder_unlabelled(graph, r, max_candidates, num_threads);
+    assert(!rv.empty());
+
+    // Reverse permutation
+    for (std::int32_t q : rv)
+      r[q] = count++;
+  }
+
+  // Check all labelled
+  assert(std::find(r.begin(), r.end(), -1) == r.end());
+  return r;
+}
+//-----------------------------------------------------------------------------
+std::vector<std::int32_t>
+graph::reorder_rcm(const graph::AdjacencyList<std::int32_t>& graph)
 {
   const std::int32_t n = graph.num_nodes();
   std::vector<std::int32_t> r(n, -1);
@@ -370,7 +534,7 @@
   int count = 0;
   while (count < n)
   {
-    rv = gps_reorder_unlabelled(graph, r);
+    rv = rcm_reorder_unlabelled(graph, r);
     assert(!rv.empty());
 
     // Reverse permutation

cpp/dolfinx/graph/ordering.h

@@ -7,6 +7,7 @@
 #pragma once
 
 #include "AdjacencyList.h"
+#include <cstddef>
 #include <cstdint>
 #include <vector>
 
@@ -18,10 +19,76 @@ namespace dolfinx::graph
 /// Bandwidth and Profile of a Sparse Matrix*, SIAM Journal on Numerical
 /// Analysis, 13(2): 236-250, 1976, https://doi.org/10.1137/0713023.
 ///
+/// At each step the pseudo-diameter search (the dominant cost for
+/// dense graphs, e.g. dof-adjacency graphs from higher-order elements)
+/// tests candidate root vertices from the final level of a level
+/// structure, in increasing degree order, to find the pair of vertices
+/// that minimises level width. Testing every candidate costs
+/// O(|S| * (V+E)) per step, where |S| is the size of the final level,
+/// which can dominate run time for large graphs. `max_candidates`
+/// bounds how many of the lowest-degree candidates are tried at each
+/// step; the (up to `max_candidates`) candidates that are tried are
+/// evaluated concurrently across `num_threads` threads. The result is
+/// identical for any `num_threads`, since candidates are only
+/// evaluated concurrently, not selected concurrently: the choice of
+/// which candidate to use is still made sequentially afterwards, in
+/// the original order.
+///
+/// A small `max_candidates` (e.g. 5) is a good trade-off for simple,
+/// convex, roughly uniform-density meshes, where it matches an
+/// exhaustive search at a fraction of the cost. For non-convex domains
+/// or strongly graded meshes (e.g. boundary-layer meshes) a small cap
+/// can measurably worsen the resulting bandwidth and profile; pass a
+/// larger value if ordering quality matters more than reordering time
+/// for such meshes.
+///
+/// To recover the original (uncapped) Gibbs-Poole-Stockmeyer algorithm
+/// as published -- testing every candidate in the final level, with no
+/// truncation -- pass `max_candidates =
+/// std::numeric_limits<std::size_t>::max()` (or any value at least as
+/// large as the largest final level `S` encountered; the actual number
+/// of candidates tested is `std::min(max_candidates, S.size())`, so an
+/// oversized value is never invalid, just unnecessary).
+///
+/// @param[in] graph The graph to compute a re-ordering for
+/// @param[in] max_candidates Maximum number of final-level candidates
+/// to test at each step of the pseudo-diameter search. Pass
+/// `std::numeric_limits<std::size_t>::max()` for the original,
+/// exhaustive algorithm (see above).
+/// @param[in] num_threads Number of threads to use for the
+/// pseudo-diameter candidate search. `1` runs serially.
+/// @return Reordering array `map`, where `map[i]` is the new index of
+/// node `i`.
+std::vector<std::int32_t>
+reorder_gps(const graph::AdjacencyList<std::int32_t>& graph,
+            std::size_t max_candidates, std::size_t num_threads);
+
+/// @brief Re-order a graph using the Reverse Cuthill-McKee algorithm.
+///
+/// The algorithm is described in *Reducing the Bandwidth of Sparse
+/// Symmetric Matrices*, Proceedings of the 1969 24th National
+/// Conference, ACM, 1969, pp. 157-172,
+/// https://doi.org/10.1145/800195.805928. The pseudo-peripheral root
+/// used to start the ordering is found using the George-Liu "double
+/// sweep" heuristic (as in `reorder_gps`'s Algorithm I, but trying only
+/// the single lowest-degree candidate at each step).
+///
+/// Unlike `reorder_gps`, there is no width-minimising second phase: a
+/// single level structure is built from the pseudo-peripheral root,
+/// each level is numbered in increasing degree order, and the whole
+/// numbering is reversed (the "reverse" in Reverse Cuthill-McKee, which
+/// tends to reduce profile relative to the plain, non-reversed
+/// numbering). This makes `reorder_rcm` an O(V+E) algorithm with a much
+/// smaller constant than `reorder_gps`, at the cost of the bandwidth
+/// and profile quality that `reorder_gps`'s extra phase can provide on
+/// non-convex or strongly graded meshes -- for simple, convex,
+/// roughly uniform-density meshes the two typically produce comparable
+/// bandwidth.
+///
 /// @param[in] graph The graph to compute a re-ordering for
 /// @return Reordering array `map`, where `map[i]` is the new index of
-/// node `i`
+/// node `i`.
 std::vector<std::int32_t>
-reorder_gps(const graph::AdjacencyList<std::int32_t>& graph);
+reorder_rcm(const graph::AdjacencyList<std::int32_t>& graph);
 
 } // namespace dolfinx::graph