Function gps_reorder_unlabelled

cpp/dolfinx/graph/ordering.cpp:129–357

----------------------------------------------------------------------------- Gibbs-Poole-Stockmeyer algorithm, finding a reordering for the given graph, operating only on nodes which are yet unlabelled (indicated with -1 in the vector rlabel).

// graph, operating only on nodes which are yet unlabelled (indicated
// 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)
{
  common::Timer timer("Gibbs-Poole-Stockmeyer 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); };

  // ALGORITHM I. Finding endpoints of a pseudo-diameter.

  // A. 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;
    }
  }

  // B. Generate a level structure Lv rooted at vertex v.
  graph::AdjacencyList<int> lv = create_level_structure(graph, v);
  graph::AdjacencyList<int> lu(0);
  bool done = false;
  int u = 0;
  std::vector<int> S;
  while (!done)
  {
    // Sort final level S of Lv into increasing degree order
    auto lv_final = lv.links(lv.num_nodes() - 1);
    S.resize(lv_final.size());
    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)
    {
      graph::AdjacencyList<int> lstmp = create_level_structure(graph, s);
      if (lstmp.num_nodes() > lv.num_nodes())
      {
        // Found a deeper level structure, so restart
        v = s;
        lv = lstmp;
        done = false;
        break;
      }

      //  D. Let u be the vertex of S whose associated level structure
      //  has smallest width