Mesh build performance improvements (#4110)

* Small improvements

* Speed up topology

* Performance improvements

* Small improrvements

* Test fix

Author: garth-wells

cpp/dolfinx/graph/ordering.cpp

@@ -75,12 +75,14 @@ residual_graph_components(const graph::AdjacencyList<int>& graph,
 }
 //-----------------------------------------------------------------------------
 // Get the (maximum) width of a level structure
-int max_level_width(const graph::AdjacencyList<int>& levels)
+std::size_t max_level_width(const graph::AdjacencyList<int>& levels)
 {
-  int wmax = 0;
-  for (int i = 0; i < levels.num_nodes(); ++i)
-    wmax = std::max(wmax, levels.num_links(i));
-  return wmax;
+  const std::vector<std::int32_t>& offsets = levels.offsets();
+  return std::transform_reduce(
+      offsets.begin(), std::prev(offsets.end()), std::next(offsets.begin()),
+      std::size_t(0),
+      [](auto x0, auto x1) -> std::size_t { return std::max(x1, x0); },
+      [](auto x0, auto x1) -> std::size_t { return x1 - x0; });
 }
 //-----------------------------------------------------------------------------
 // Create a level structure from graph, rooted at node s
@@ -120,7 +122,6 @@ create_level_structure(const graph::AdjacencyList<int>& graph, int s)
   return graph::AdjacencyList(std::move(level_structure),
                               std::move(level_offsets));
 }
-
 //-----------------------------------------------------------------------------
 // Gibbs-Poole-Stockmeyer algorithm, finding a reordering for the given
 // graph, operating only on nodes which are yet unlabelled (indicated
@@ -134,7 +135,7 @@ gps_reorder_unlabelled(const graph::AdjacencyList<std::int32_t>& graph,
   const std::int32_t n = graph.num_nodes();
 
   // Degree comparison function
-  auto cmp_degree = [&graph](int a, int b)
+  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.
@@ -164,7 +165,6 @@ gps_reorder_unlabelled(const graph::AdjacencyList<std::int32_t>& graph,
     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;
 
@@ -258,15 +258,15 @@ gps_reorder_unlabelled(const graph::AdjacencyList<std::int32_t>& graph,
                              [](int vl, int vn) { return (vl > vn) ? vl : 0; });
 
       // Find maximum of those that did increase
-      int h0 = *std::ranges::max_element(wh);
-      int l0 = *std::ranges::max_element(wl);
+      int h0 = std::ranges::max(wh);
+      int l0 = std::ranges::max(wl);
 
       // Choose which side to use
       int side = h0 < l0 ? 0 : 1;
 
-      // If h0 == l0, then use the elements of the level pairs which arise
-      // from the rooted level structure of smaller width. If the widths are
-      // equal, use the first elements. (i.e. lvp[][0]).
+      // If h0 == l0, then use the elements of the level pairs which
+      // arise from the rooted level structure of smaller width. If the
+      // widths are equal, use the first elements. (i.e. lvp[][0]).
       if (h0 == l0)
         side = max_level_width(lu) < max_level_width(lv) ? 1 : 0;
 
@@ -336,10 +336,12 @@ gps_reorder_unlabelled(const graph::AdjacencyList<std::int32_t>& graph,
       // with lowest degree
       nrem.clear();
       for (int w : lslevel)
+      {
         if (!labelled[w])
           nrem.push_back(w);
+      }
 
-      if (nrem.size() == 0)
+      if (nrem.empty())
         break;
 
       std::ranges::sort(nrem, cmp_degree);

cpp/dolfinx/graph/partition.cpp

@@ -1,4 +1,4 @@
-// Copyright (C) 2020 Garth N. Wells
+// Copyright (C) 2020-2026 Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -8,6 +8,7 @@
 #include "AdjacencyList.h"
 #include "partitioners.h"
 #include <algorithm>
+#include <boost/sort/sort.hpp>
 #include <dolfinx/common/MPI.h>
 #include <dolfinx/common/Timer.h>
 #include <dolfinx/common/log.h>
@@ -390,12 +391,12 @@ graph::build::distribute(MPI_Comm comm, std::span<const std::int64_t> list,
           std::move(ghost_index_owner)};
 }
 //-----------------------------------------------------------------------------
-std::vector<std::int64_t>
-graph::build::compute_ghost_indices(MPI_Comm comm,
-                                    std::span<const std::int64_t> owned_indices,
-                                    std::span<const std::int64_t> ghost_indices,
-                                    std::span<const int> ghost_owners)
+std::vector<std::int64_t> graph::build::compute_ghost_indices(
+    MPI_Comm comm, std::span<const std::int64_t> owned_indices,
+    std::span<const std::int64_t> ghost_indices,
+    std::span<const int> ghost_owners, int num_threads)
 {
+  common::Timer timer("Compute ghost indices");
   spdlog::info("Compute ghost indices");
 
   // Get number of local cells determine global offset
@@ -424,7 +425,8 @@ graph::build::compute_ghost_indices(MPI_Comm comm,
 
   MPI_Comm neighbor_comm_fwd, neighbor_comm_rev;
 
-  std::vector<int> in_edges = MPI::compute_graph_edges_pcx(comm, neighbors);
+  std::vector<int> in_edges
+      = dolfinx::MPI::compute_graph_edges_pcx(comm, neighbors);
   MPI_Dist_graph_create_adjacent(comm, in_edges.size(), in_edges.data(),
                                  MPI_UNWEIGHTED, neighbors.size(),
                                  neighbors.data(), MPI_UNWEIGHTED,
@@ -487,7 +489,13 @@ graph::build::compute_ghost_indices(MPI_Comm comm,
     old_to_new.push_back(
         {idx, static_cast<std::int64_t>(offset_local + old_to_new.size())});
   }
-  std::ranges::sort(old_to_new);
+  if (num_threads > 0)
+  {
+    boost::sort::block_indirect_sort(old_to_new.begin(), old_to_new.end(),
+                                     num_threads);
+  }
+  else
+    std::ranges::sort(old_to_new);
 
   // Replace values in recv_data with new_index and send back
   std::ranges::transform(recv_data, recv_data.begin(),

cpp/dolfinx/graph/partition.h

@@ -114,12 +114,13 @@ distribute(MPI_Comm comm, std::span<const std::int64_t> list,
 /// @param[in] ghost_indices List of ghost global indices.
 /// @param[in] ghost_owners The owning rank for each entry in
 /// `ghost_indices`.
+/// @param[in] num_threads Number of threads to use.
 /// @return New global indices for the ghost indices.
 std::vector<std::int64_t>
 compute_ghost_indices(MPI_Comm comm,
                       std::span<const std::int64_t> owned_indices,
                       std::span<const std::int64_t> ghost_indices,
-                      std::span<const int> ghost_owners);
+                      std::span<const int> ghost_owners, int num_threads);
 
 /// Given an adjacency list with global, possibly non-contiguous, link
 /// indices and a local adjacency list with contiguous link indices

cpp/dolfinx/mesh/Topology.cpp

@@ -689,19 +689,38 @@ std::vector<std::array<std::int64_t, 3>> exchange_ghost_indexing(
 /// @param[in] global_to_local Sorted array of (global, local) indices.
 std::vector<std::int32_t> convert_to_local_indexing(
     std::span<const std::int64_t> g,
-    std::span<const std::pair<std::int64_t, std::int32_t>> global_to_local)
+    std::span<const std::pair<std::int64_t, std::int32_t>> global_to_local,
+    int num_threads)
 {
+  auto transform =
+      [](std::span<std::int32_t> data, std::span<const std::int64_t> g,
+         std::span<const std::pair<std::int64_t, std::int32_t>> global_to_local)
+  {
+    std::transform(g.begin(), g.end(), data.begin(),
+                   [&global_to_local](auto i)
+                   {
+                     auto it = std::ranges::lower_bound(
+                         global_to_local, i, std::ranges::less(),
+                         [](auto& e) { return e.first; });
+                     assert(it != global_to_local.end());
+                     assert(it->first == i);
+                     return it->second;
+                   });
+  };
+
   std::vector<std::int32_t> data(g.size());
-  std::transform(g.begin(), std::next(g.begin(), data.size()), data.begin(),
-                 [&global_to_local](auto i)
-                 {
-                   auto it = std::ranges::lower_bound(
-                       global_to_local, i, std::ranges::less(),
-                       [](auto& e) { return e.first; });
-                   assert(it != global_to_local.end());
-                   assert(it->first == i);
-                   return it->second;
-                 });
+  if (num_threads > 0)
+  {
+    std::vector<std::jthread> threads(num_threads);
+    for (int i = 0; i < num_threads; ++i)
+    {
+      auto [c0, c1] = dolfinx::MPI::local_range(i, g.size(), num_threads);
+      threads[i] = std::jthread(transform, std::span(data.data() + c0, c1 - c0),
+                                g.subspan(c0, c1 - c0), global_to_local);
+    }
+  }
+  else
+    transform(data, g, global_to_local);
 
   return data;
 }
@@ -733,7 +752,6 @@ build_entity_types(const std::vector<CellType>& cell_types)
   }
   return entity_types;
 }
-
 } // namespace
 
 //-----------------------------------------------------------------------------
@@ -834,9 +852,10 @@ std::shared_ptr<const common::IndexMap> Topology::index_map(int dim) const
       = this->index_maps(dim);
   if (im.empty())
   {
-    throw std::runtime_error(std::format(
-        "Missing IndexMap in Topology. Maybe you need to create_entities({}).",
-        dim));
+    throw std::runtime_error(
+        std::format("Missing IndexMap in Topology. Maybe you need to "
+                    "create_entities({}).",
+                    dim));
   }
 
   return im.at(0);
@@ -1034,7 +1053,7 @@ Topology mesh::create_topology(
     std::vector<std::span<const std::int64_t>> cells,
     std::vector<std::span<const std::int64_t>> original_cell_index,
     std::vector<std::span<const int>> ghost_owners,
-    std::span<const std::int64_t> boundary_vertices)
+    std::span<const std::int64_t> boundary_vertices, int num_threads)
 {
   common::Timer timer("Topology: create");
 
@@ -1044,6 +1063,7 @@ Topology mesh::create_topology(
 
   // Check cell data consistency and compile spans of owned and ghost
   // cells
+  common::Timer timer0("Topology: 0");
   spdlog::info("Create topology (generalised)");
   std::vector<std::int32_t> num_local_cells(cell_types.size());
   std::vector<std::span<const std::int64_t>> owned_cells;
@@ -1064,6 +1084,11 @@ Topology mesh::create_topology(
     ghost_cells.push_back(cells[i].last(ghost_owners[i].size() * num_vertices));
   }
 
+  timer0.stop();
+  timer0.flush();
+
+  common::Timer timer1("Topology: 1");
+
   // Create sets of owned and unowned vertices from the cell ownership
   // and the list of boundary vertices
   auto [owned_vertices, unowned_vertices]
@@ -1075,9 +1100,14 @@ Topology mesh::create_topology(
   const graph::AdjacencyList<int> global_vertex_to_ranks
       = determine_sharing_ranks(comm, boundary_vertices);
 
+  timer1.stop();
+  timer1.flush();
+
   // Iterate over vertices that have 'unknown' ownership, and if flagged
   // as owned by determine_sharing_ranks update ownership status
   {
+    common::Timer timer2("Topology: 2");
+
     const int mpi_rank = dolfinx::MPI::rank(comm);
     std::vector<std::int64_t> owned_shared_vertices;
     for (std::size_t i = 0; i < boundary_vertices.size(); ++i)
@@ -1102,13 +1132,15 @@ Topology mesh::create_topology(
     dolfinx::radix_sort(owned_vertices);
   }
 
+  // NOTE: This block is relatively expensive
   // Number all owned vertices, iterating over vertices cell-wise
   std::vector<std::int32_t> local_vertex_indices(owned_vertices.size(), -1);
   {
+    common::Timer timer3("Topology: 3 ");
     std::int32_t v = 0;
-    for (std::size_t i = 0; i < cell_types.size(); ++i)
+    for (std::span<const std::int64_t> cells_t : cells)
     {
-      for (auto vtx : cells[i])
+      for (auto vtx : cells_t)
       {
         if (auto it = std::ranges::lower_bound(owned_vertices, vtx);
             it != owned_vertices.end() and *it == vtx)
@@ -1121,6 +1153,8 @@ Topology mesh::create_topology(
     }
   }
 
+  common::Timer timer4("Topology: 4");
+
   // Compute the global offset for owned (local) vertex indices
   std::int64_t global_offset_v = 0;
   {
@@ -1136,29 +1170,39 @@ Topology mesh::create_topology(
     std::span cell_idx(original_cell_index[i]);
     cell_ghost_indices.push_back(graph::build::compute_ghost_indices(
         comm, cell_idx.first(num_local_cells[i]),
-        cell_idx.last(ghost_owners[i].size()), ghost_owners[i]));
+        cell_idx.last(ghost_owners[i].size()), ghost_owners[i], num_threads));
 
     // Create index maps for each cell type
     index_map_c.push_back(std::make_shared<common::IndexMap>(
         comm, num_local_cells[i], cell_ghost_indices[i], ghost_owners[i],
         static_cast<int>(dolfinx::MPI::tag::consensus_nbx) + i));
   }
 
+  timer4.stop();
+  timer4.flush();
+
   // Send and receive  ((input vertex index) -> (new global index, owner
   // rank)) data with neighbours (for vertices on 'true domain
   // boundary')
+  common::Timer timer5("Topology: 5");
+
   const std::vector<std::int64_t> unowned_vertex_data = exchange_indexing(
       comm, boundary_vertices, global_vertex_to_ranks, global_offset_v,
       owned_vertices, local_vertex_indices);
   assert(unowned_vertex_data.size() % 3 == 0);
 
+  timer5.stop();
+  timer5.flush();
+
   // Unpack received data and build array of ghost vertices and owners
   // of the ghost vertices
   std::vector<std::int64_t> ghost_vertices;
   std::vector<int> ghost_vertex_owners;
   std::vector<std::int32_t> local_vertex_indices_unowned(
       unowned_vertices.size(), -1);
   {
+    common::Timer timer6("Topology: 6");
+
     std::int32_t v = owned_vertices.size();
     for (std::size_t i = 0; i < unowned_vertex_data.size(); i += 3)
     {
@@ -1244,14 +1288,18 @@ Topology mesh::create_topology(
       { return {idx0, idx1}; });
   std::ranges::sort(global_to_local_vertices);
 
+  common::Timer timer7("Topology: 7");
   std::vector<std::vector<std::int32_t>> _cells_local_idx;
   _cells_local_idx.reserve(cells.size());
   for (std::span<const std::int64_t> c : cells)
   {
     _cells_local_idx.push_back(
-        convert_to_local_indexing(c, global_to_local_vertices));
+        convert_to_local_indexing(c, global_to_local_vertices, num_threads));
   }
 
+  timer7.stop();
+  timer7.flush();
+
   // -- Create Topology object
 
   // Determine which ranks ghost vertices that are owned by this rank.
@@ -1309,11 +1357,12 @@ Topology
 mesh::create_topology(MPI_Comm comm, std::span<const std::int64_t> cells,
                       std::span<const std::int64_t> original_cell_index,
                       std::span<const int> ghost_owners, CellType cell_type,
-                      std::span<const std::int64_t> boundary_vertices)
+                      std::span<const std::int64_t> boundary_vertices,
+                      int num_threads)
 {
   spdlog::info("Create topology (single cell type)");
   return create_topology(comm, {cell_type}, {cells}, {original_cell_index},
-                         {ghost_owners}, boundary_vertices);
+                         {ghost_owners}, boundary_vertices, num_threads);
 }
 //-----------------------------------------------------------------------------
 std::tuple<Topology, std::vector<int32_t>, std::vector<int32_t>>