Parallel BFS, in depth: level-synchronous frontiers, direction-optimizing search and what limits speedup

Why BFS resists parallelism

Sequential BFS uses one FIFO queue. Parallel BFS replaces it with a frontier: the set of vertices at distance d. All of them can be expanded at once, because nothing at distance d depends on anything else at distance d. Expanding them produces the next frontier at distance d + 1, the threads meet at a barrier, and the frontiers swap. This is level-synchronous BFS, and nearly every shared-memory, GPU and distributed implementation starts from it.

The graph should be in compressed sparse row (CSR) form: an offsets array of length n + 1 and a neighbours array of length 2m for an undirected graph, so the neighbours of u sit in nbr[off[u] .. off[u+1]). It is compact, scans sequentially within a vertex, and splits cleanly across threads. Adjacency lists of pointers waste cache lines and memory bandwidth on exactly the operation BFS repeats billions of times.

Two things make it hard. First, two threads can discover the same vertex in the same level, so claiming a vertex must be atomic or the next frontier gets duplicates. Second, the frontier on a low-diameter graph such as a social network goes from tiny to most of the graph within two or three levels, and in those big levels nearly every edge scanned leads to a vertex that is already visited. That wasted work is the target of the direction-optimizing technique.

The two step types

The top-down step gives each thread a slice of the frontier. For every neighbour it does a cheap relaxed load first, so already-visited vertices cost no atomic, then a compare-and-swap that only one thread can win. Each thread appends winners to its own buffer, so there is no contended shared queue. CSR and Bitmap are minimal helper types; the rest is standard C++20 and OpenMP.

// One top-down level. parent[v] == -1 means unvisited.
int64_t td_step(const CSR& g, std::vector<int32_t>& parent,
                const std::vector<int32_t>& frontier, std::vector<int32_t>& next) {
  int64_t scout = 0;                         // degree sum of the new frontier
  std::vector<std::vector<int32_t>> local(omp_get_max_threads());
  #pragma omp parallel reduction(+ : scout)
  {
    auto& buf = local[omp_get_thread_num()];
    #pragma omp for schedule(dynamic, 64) nowait
    for (size_t i = 0; i < frontier.size(); ++i) {
      int32_t u = frontier[i];
      for (int64_t e = g.off[u]; e < g.off[u + 1]; ++e) {
        int32_t v = g.nbr[e];
        std::atomic_ref<int32_t> pv(parent[v]);
        int32_t expected = -1;
        if (pv.load(std::memory_order_relaxed) == -1 &&
            pv.compare_exchange_strong(expected, u, std::memory_order_relaxed)) {
          buf.push_back(v);                  // exactly one thread wins v
          scout += g.off[v + 1] - g.off[v];
        }
      }
    }
  }
  next.clear();                              // a prefix sum over buffer sizes parallelises this copy
  for (auto& b : local) next.insert(next.end(), b.begin(), b.end());
  return scout;
}

// One bottom-up level over a bitmap frontier. No CAS on parent: only v's owner writes it.
int64_t bu_step(const CSR& g, std::vector<int32_t>& parent,
                const Bitmap& front, Bitmap& next) {
  int64_t awake = 0;
  next.reset();
  #pragma omp parallel for reduction(+ : awake) schedule(dynamic, 1024)
  for (int32_t v = 0; v < g.n; ++v) {
    if (parent[v] != -1) continue;
    for (int64_t e = g.off[v]; e < g.off[v + 1]; ++e) {
      int32_t u = g.nbr[e];
      if (front.get(u)) {
        parent[v] = u;
        next.set_bit_atomic(v);              // bits share words, so this one is atomic
        ++awake;
        break;                               // the early exit is the whole point
      }
    }
  }
  return awake;
}

The bottom-up step inverts the loop. Every unvisited vertex scans its own neighbours looking for any member of the current frontier, held as a bitmap with one bit per vertex. The first hit is a valid parent, so the scan stops. When the frontier covers a large share of the graph, most unvisited vertices find a parent within their first few neighbours, and the edges that top-down would have wasted on visited vertices are never touched. No atomic is needed on parent because each v is handled by exactly one thread; only the shared words of the next bitmap need atomic bit sets.

Dynamic scheduling matters in both loops. Degree distributions are skewed, so a static split leaves one thread holding a hub with a million edges while the rest idle at the barrier.

Direction-optimizing BFS

Direction-optimizing BFS, introduced by Scott Beamer, Krste Asanovic and David Patterson in 2012, runs top-down while the frontier is small and bottom-up while it is large. The switch uses two counts. Scout is the sum of degrees of the frontier, an estimate of the edges top-down would scan next. Edges-to-check is the number of edges not yet scanned. The reference implementation in the GAP Benchmark Suite, which Beamer maintains, switches to bottom-up when scout exceeds edges-to-check divided by alpha, and switches back when the frontier is shrinking and holds at most n divided by beta vertices. Its defaults are alpha = 15 and beta = 18; treat them as starting points to tune per machine and graph family.

Here is a sequential Python reference of the same control logic, useful for checking a parallel implementation's parent tree. It was fuzzed against queue BFS on 300 random graphs with alpha and beta varied, checking that every vertex's depth in the parent tree equals its BFS distance.

def bfs_hybrid(offsets, nbr, source, alpha=15, beta=18):
    n = len(offsets) - 1
    parent = [-1] * n
    parent[source] = source
    frontier, edges_to_check = [source], len(nbr)
    scout = offsets[source + 1] - offsets[source]
    while frontier:
        if scout > edges_to_check / alpha:            # frontier is heavy: go bottom-up
            in_front = set(frontier)
            while True:
                old, nxt = len(frontier), []
                for v in range(n):
                    if parent[v] == -1:
                        for i in range(offsets[v], offsets[v + 1]):
                            if nbr[i] in in_front:
                                parent[v] = nbr[i]; nxt.append(v); break
                frontier, in_front = nxt, set(nxt)
                if not frontier or (len(frontier) < old and len(frontier) <= n / beta):
                    break                             # shrinking and small: back to top-down
            scout = 1
        else:
            edges_to_check -= scout
            nxt, scout = [], 0
            for u in frontier:
                for i in range(offsets[u], offsets[u + 1]):
                    v = nbr[i]
                    if parent[v] == -1:
                        parent[v] = u; nxt.append(v)
                        scout += offsets[v + 1] - offsets[v]
            frontier = nxt
    return parent

Worked run

The measured run used a 200,000-vertex undirected graph with 1.6 million edges (3.2 million directed CSR entries), with one endpoint of each edge skewed toward low ids so a few vertices become hubs, starting from vertex 0. Every vertex was reached and the deepest was four levels away.

Plain top-down BFS examined all 3,200,000 directed edges, as it always does on a connected graph. The hybrid ran top-down for level 1 (3,564 vertices discovered) and level 2 (47,149), when scout first exceeded one fifteenth of the remaining edges. It then ran bottom-up for level 3 (146,339 vertices, three quarters of the graph) and level 4 (2,947), dropped back to top-down once the frontier shrank below n / 18, and found nothing more. It examined 606,013 edges in total, 5.3 times fewer, and produced a valid BFS tree. Wall-clock speedups on real hardware depend on memory behaviour too, but the edge count is the dominant term, which is why the technique transfers across CPUs, GPUs and clusters.

Scaling out: frontiers, GPUs and clusters

• Building the next frontier. Thread-local buffers plus a parallel prefix sum over their sizes give each thread its write offset in the shared output; the parallel prefix sum article covers the scan. Bitmaps are better when the frontier is dense.
• GPUs. The same split appears with more severe imbalance: a warp that draws a hub stalls. Merrill, Garland and Grimshaw's 2012 GPU traversal work combines per-thread, per-warp and per-block gathering by degree, and frameworks such as Gunrock implement direction switching. The GPU algorithms article explains the execution model.
• Clusters. With a 1D partition each rank owns a vertex range and sends discovered vertices to their owners every level, an all-to-all whose volume grows with edges. 2D partitioning of the adjacency matrix, studied by Buluc and Madduri, limits each exchange to a row or column of processors.
• Vertex ordering. Relabelling vertices so neighbours have nearby ids, for example by degree or a prior BFS order, improves cache and bitmap locality and can matter as much as the parallel scheme.

Operational guidance

Benchmark with TEPS on graphs shaped like yours, and report the edges examined alongside the time so you can tell algorithmic savings from hardware effects. Pin threads, allocate the CSR arrays with first-touch initialisation on the NUMA node that will scan them, and use huge pages for the neighbour array: BFS is TLB-hungry. Start a few hundred BFS runs from random sources, because one source can land in a small component and distort the average.

Keep the parent array as the output. Distances are unique, but parents are not: with several valid parents, which thread wins is a race, so the tree varies run to run. Validate with the Graph500 rules rather than by comparing parents: the source is its own parent, every tree edge exists in the graph, and depths differ by exactly one along tree edges and by at most one across any graph edge.

Failure modes

• High-diameter graphs. Road networks and meshes have thousands of levels with tiny frontiers. Barrier cost dominates and bottom-up never pays; use top-down with cheap synchronisation, or a different algorithm such as delta-stepping for weighted variants.
• Directed graphs. Bottom-up needs in-edges. Without a transposed CSR the bottom-up step silently finds the wrong parents.
• Plain reads and writes racing. A non-atomic check-then-write lets two threads claim one vertex, duplicating it in the next frontier and inflating work.
• False sharing. Per-thread counters packed into one cache line, or per-thread buffers resized in place, serialise the threads.
• Mis-tuned switch. A beta that is too large keeps bottom-up running over nearly empty frontiers, where it scans every unvisited vertex for nothing.

Trade-offs

What to do next

1. Convert your graph to CSR, and build the transpose if it is directed.
2. Implement level-synchronous top-down with CAS claims and thread-local buffers, and validate the parent tree with the Graph500 rules.
3. Add the bottom-up step and the alpha/beta switch, starting from GAP's 15 and 18, and log edges examined per level as in the run above.
4. Sweep alpha and beta on your own graphs, measure TEPS over many random sources, and keep the setting that wins on the median, not the best run.
5. Try vertex relabelling and huge pages before reaching for more threads.
6. For related patterns, read multi-source BFS and the BFS and DFS comparison.