Graph traversal on a GPU
Here is a frontier from a real graph. It holds ten vertices.
Give each vertex a thread, which is the mapping this course has used since
day 5, and the launch is one
warp. This gives the ten vertices one warp in total, not one
warp per vertex. That warp runs its neighbour loop 1,255 times because the
word the has 1,255 neighbours.
The word and has 1,012 neighbours. The other eight finish within 43
iterations, so their lanes sit idle for the rest. The
per-word degrees come from code/day38-bfs/word-graph.txt; the evidence
transcript prints only each level's maximum.
The frontier holds 2,452 edges between them. The warp issues 40,160 lane slots to read them.
By the end of this page you can look at a frontier and say which vertex sets the clock, and you will have measured two ways to stop it.
What a level-synchronous search actually launches
Breadth-first search visits a graph by distance from a source. Level 0 is the source. Level 1 is everything one hop away, level 2 is everything two hops away that was not already seen, and the search stops when a level is empty.
The set of vertices being expanded now is the frontier. A GPU BFS must choose how to map that frontier to threads.
The version without a frontier is the one most people write first, and it is worth writing down because it is correct:
// Every thread checks every vertex, on every level.
if (dist[v] == level) {
for (int e = rowOffsets[v]; e < rowOffsets[v + 1]; ++e) {
// relax the neighbour
}
}
Launch that once per level and the answer is right. It also reads all 3,485 vertices ten times, which is 34,850 vertex visits to expand 3,434 vertices. The extra work grows with the number of levels, so it can be large on a graph with a large diameter, such as a road network.
Keeping a frontier needs one array and one counter and removes that factor. The loop that drives it is short:
while (frontierSize > 0) {
if (levelSizes != nullptr && levels < kMaxLevels) {
levelSizes[levels] = frontierSize;
}
CUDA_CHECK(cudaMemset(d.nextSize, 0, sizeof(unsigned int)));
launchExpand(d, frontier, frontierSize, nextFrontier, levels,
useWarpMapping);
// The host has to see this count before it can size the next launch,
// so every level costs a round trip as well as a launch.
unsigned int produced = 0u;
CUDA_CHECK(cudaMemcpy(&produced, d.nextSize, sizeof(unsigned int),
cudaMemcpyDeviceToHost));
frontierSize = static_cast<int>(produced);
std::swap(frontier, nextFrontier);
++levels;
}
Read the loop condition again. The number of kernel launches is the number of levels, which is the source's eccentricity plus one and, over the worst source, the graph's diameter plus one. The number of streaming multiprocessors does not change this count.
Each launch pays its launch overhead, whether
the frontier holds one vertex or ten thousand. The host also needs the
cudaMemcpy result before it can size the next grid or stop the search, so
each level adds a synchronization point.
Here is the shipped graph. Vertices are 3,485 words from this project's own
research documents, joined when two words sit next to each other three or more
times; the source is await, one of the four words whose search needs ten
levels. Every column comes from the graph, not from a clock.
| Level | Frontier | Edges | Max degree | Busiest | Warps | Lane slots |
|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 1 | await |
1 | 32 |
| 1 | 1 | 2 | 2 | readfile |
1 | 64 |
| 2 | 1 | 4 | 4 | join |
1 | 128 |
| 3 | 3 | 13 | 8 | root |
1 | 256 |
| 4 | 10 | 2,452 | 1,255 | the |
1 | 40,160 |
| 5 | 1,677 | 25,235 | 757 | a |
53 | 196,832 |
| 6 | 1,579 | 3,851 | 33 | org |
50 | 13,856 |
| 7 | 147 | 184 | 4 | frontend |
5 | 448 |
| 8 | 14 | 15 | 2 | hwu |
1 | 64 |
| 9 | 1 | 1 | 1 | kirk |
1 | 32 |
The frontier is tiny for four levels, grows to half the graph, then shrinks. That shape is common in a network with hubs, and it means the GPU is nearly empty for six of the ten launches.
Diagram: two frontiers, two mappings. Three horizontal bands compare issued lane-slots with useful edge reads. Band 1: level 4's ten vertices with one thread per vertex. One warp runs 1,255 iterations, issuing 40,160 lane-slots for 2,452 edges.
Two lanes have degrees 1,255 and 1,012, eight have degree at most 43, and 22 lanes are empty. Band 2: the same ten vertices with one warp per vertex. Ten warps issue 2,656 lane-slots for the same 2,452 edges, and the tallest list takes 40 iterations.
Band 3: level 6 reverses the trade. One thread per vertex issues 13,856 lane-slots for 3,851 edges, while one warp per vertex issues 50,560 because many short lists each occupy a 32-lane pass.
Alt text: "At BFS level 4, one-thread-per-vertex issues 40,160 lane-slots versus 2,656 for one-warp-per-vertex; at level 6 the ordering reverses, 13,856 versus 50,560."
One thread, one vertex, and where that habit breaks
The mapping you have used all course is one thread to one element, and it has worked because the elements had the same size. Day 11's copy and day 20's convolution give every thread the same amount of work, so a warp finishes when its threads finish, all together.
A graph breaks that assumption because its elements are neighbour lists of different lengths. Within one frontier, their lengths differ by three orders of magnitude. A warp runs in lockstep, so it issues 32 times the largest degree, not the sum of its 32 degrees.
That is the lane-slot column above. Over the whole search, it totals 251,872 lane slots for 31,758 edges.
This is warp divergence from
day 22, except that the branch is a loop bound
rather than an if, so you cannot fix it by rewriting the predicate. It is
also day 37's problem. A CSR matrix with one dense row does the same thing to
a one-thread-per-row SpMV, and the answer there was the same as the answer
here.
The answer is to widen the unit of work. Give one vertex to a whole warp and have the 32 lanes stride through its neighbour list together:
const unsigned int lane = threadIdx.x % kWarpSize;
const size_t w = t / kWarpSize; // warp uniform: one warp, one vertex
if (w < static_cast<size_t>(frontierSize)) {
const int u = frontier[w];
const int last = rowOffsets[u + 1];
for (int e = rowOffsets[u] + static_cast<int>(lane); e < last;
e += kWarpSize) {
claimVertex(colIndices[e], level, dist, nextFrontier, nextSize);
}
}
the now costs 40 iterations rather than 1,255, and consecutive lanes read
consecutive entries of colIndices, so the read is
coalesced instead of 32 unrelated addresses.
The cost shifts to short lists: this launches 32 threads for every vertex,
so a level made of degree-1 vertices wastes 31 lanes on each of them. Level 6
is exactly that level.
Measuring the mapping without measuring the launch
The full program is
code/day38-bfs/bfs.cu. Four choices make
its numbers worth trusting.
Both mappings walk the same edges, and the program says how many. It prints the level table above from the graph it just loaded, so the edge counts sit on the screen next to the times rather than being taken on trust.
The accumulator checks the work. The timed walk sums neighbour ids into an
unsigned int, where overflow is defined to wrap, so the host predicts the
exact value for every vertex and checks it. A walk the compiler shortened
comes back wrong instead of coming back fast. The sum is XORed with the repeat
counter for the same reason: without that the inner sum is loop invariant and
the compiler may compute it once and multiply.
Events, and a warm-up for every kernel. Day 9 covers why a host clock around a launch measures the launch, and why the first launch of each kernel pays its own module load.
The search and one level are timed separately, because they answer different questions. This is the part the first draft of this program got wrong. It timed the whole search under each mapping and compared, which cannot work: ten launches and ten host round trips sit in that number against 31,758 edges of real work, so the thing being studied is a rounding error inside the thing being paid for. Both measurements are still in the program, and the whole-search rows are there to show you that.
The isolated measurement drops the atomics and the frontier build and repeats the neighbour walk 64 times inside the kernel, because one pass over level 5 reads 25,235 entries and finishes in less time than a launch takes.
The expansion kernel itself is the loop you would write:
__global__ void bfsExpandThreadPerVertex(const int* __restrict__ frontier,
int frontierSize,
const int* __restrict__ rowOffsets,
const int* __restrict__ colIndices,
int* __restrict__ dist, int level,
int* __restrict__ nextFrontier,
unsigned int* __restrict__ nextSize) {
const size_t t = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (t < static_cast<size_t>(frontierSize)) {
const int u = frontier[t];
for (int e = rowOffsets[u]; e < rowOffsets[u + 1]; ++e) {
claimVertex(colIndices[e], level, dist, nextFrontier, nextSize);
}
}
}
claimVertex writes the new distance with an atomicCAS and appends the
vertex with an atomicAdd on one counter, which is day 33's compaction with
an atomic standing in for the scan. Day 33's version marks a flag per vertex,
scans and scatters, which is deterministic and costs a pass over all 3,485
vertices per level whether the frontier holds 1,677 vertices or one.
Note. The isolated walk measures the neighbour walk and nothing else. It does not include the atomics, the frontier append or the distance write, so it is the ceiling on what changing the mapping can buy you, not the whole story. The whole-search rows are the other end of that range.
Results
Re-verified without meaningful drift on a Tesla T4 with driver 580.173.02
and CUDA 13.0 (V13.0.88) on 2026-09-02. Graph structure, frontiers,
correctness and the mapping-speedup pattern reproduced; timing changes were
under seven percent. The original CUDA 12.6 capture below and CUDA 13
transcript both remain in evidence.
GPU: Tesla T4 (compute capability 7.5), 40 SMs
Graph: word-graph.txt
3485 vertices, 15906 undirected edges, 31812 csr entries
max degree 1255 (the), median 3, 997 vertices of degree 1
BFS from "await" (vertex 212)
level frontier edges max deg busiest warps lane slots
----- -------- ----- ------- ------- ----- ----------
0 1 1 1 await 1 32
1 1 2 2 readfile 1 64
2 1 4 4 join 1 128
3 3 13 8 root 1 256
4 10 2452 1255 the 1 40160
5 1677 25235 757 a 53 196832
6 1579 3851 33 org 50 13856
7 147 184 4 frontend 5 448
8 14 15 2 hwu 1 64
9 1 1 1 kirk 1 32
reached 3434 of 3485 vertices in 10 levels, so 10 launches
31758 edges to walk, 251872 lane slots issued to walk them
Correctness
thread per vertex distances match on all 3485 vertices
thread per vertex frontiers held 3434 vertices, no duplicates
warp per vertex distances match on all 3485 vertices
warp per vertex frontiers held 3434 vertices, no duplicates
Whole BFS, 10 timed runs, reset and 10 host round trips included
mapping ms per bfs
------------------- ----------
thread per vertex 1.8091
warp per vertex 0.3340
One level in isolation, neighbour walk repeated 64 times
level frontier edges max deg thread ms warp ms speedup
----- -------- ----- ------- --------- ------- -------
4 10 2452 1255 1.3724 0.0881 15.58x
5 1677 25235 757 1.0068 0.0796 12.64x
6 1579 3851 33 0.1172 0.0485 2.41x
all checks passed
Level 4 is the whole lesson. A frontier of ten vertices issues 40,160 lane slots to walk 2,452 edges. One of those ten is the word "the", with degree 1255, and it is in the same warp as nine vertices of degree three or four. The warp runs until its longest row finishes, so 31 lanes wait while one works.
That is day 37's load imbalance again, and it is worse here because you cannot choose the distribution. A sparse matrix has whatever row lengths it has, but at least they are fixed.
A BFS frontier changes shape at every level: level 3 has three vertices, level 5 has 1,677, and level 9 has one. No single launch configuration is right for all of them.
Across the whole traversal, 251,872 lane slots walked 31,758 edges, which is 12.6 percent useful. GPU BFS methods aim to raise that share.
The warp mapping wins each measured case, and the margin follows the skew. The whole search takes 0.3340 ms with warps and 1.8091 ms with threads, a 5.4x difference. The isolated speedups are 15.58x on level 4, 12.64x on level 5, and 2.41x on level 6.
Level 4 has one vertex of degree 1255, while level 6 has a maximum degree of 33. The gain falls as the degree spread narrows. The whole-search gain is also smaller because both mappings pay for ten launches and host round trips.
Ten levels means ten kernel launches, because a level cannot start until the previous one has finished writing its frontier. The graph's diameter, not its size, sets the number of launches, and each one pays the overhead day 9 measured.
Run it yourself
Use Colab or a local CUDA GPU with compute capability 7.5 or
later. The build line is in the repo's README. Run the program from its own
directory because it opens word-graph.txt beside it:
nvcc -std=c++17 -O3 -arch=sm_75 -o bfs bfs.cu
Compiler Explorer cannot run this example because it needs the graph file. An embed has one editor buffer and no filesystem, while this graph uses 173 KB. The program itself allocates about 200 KiB on the device and finishes within the 20-second run cap.
Exercise
Add a third expansion kernel that maps one thread to one edge rather than one vertex, using the scan from day 32, and report the three isolated rows again with your kernel as a fourth column.
Time: 30 to 45 minutes. Submit: the kernel, and the three ratios against the thread mapping.
Check: the harness puts your kernel through the two gates the shipped mappings pass. Distances are compared against the host BFS exactly on all 3,485 vertices, and the first disagreement prints the vertex id, the word and both values. The frontier sizes then have to add up to the 3,434 vertices reached, which is what catches a vertex enqueued twice.
Hint 1
One thread per edge means every thread does exactly one unit of work, so no warp waits for a wide vertex. The only thing a thread is missing is which vertex its edge belongs to. What would have to be computed before the launch for a thread to answer that from its own index?
Hint 2
An exclusive scan of the frontier's degrees gives an array whose entry i is
the first edge owned by frontier vertex i, and it is sorted. Thread e
wants the last entry not greater than e.
Solution
Scan the degrees, launch totalEdges threads, binary search the scanned array
for each thread's vertex, then read one neighbour and claim it. The imbalance
leaves the arithmetic: the lane-slot count becomes the edge count, so level 4
issues 2,452 slots instead of 40,160. Each level requires a scan and a search,
both small, and both hitting the same
L1 and L2 cached
global memory from every lane of a warp.
Pitfalls
Your distances are right and the program is slow. Claiming a vertex with
a plain if (dist[v] < 0) { dist[v] = level + 1; } instead of an atomicCAS
gives correct distances, because every racing thread writes the same value,
and each of them also appends v. The frontier fills with duplicates and
every duplicate expands again. Check that the frontier sizes add up to the
number of vertices reached; the distance check cannot see this.
You time the whole search to compare two expansion kernels. On a graph whose levels are cheap, that number is launches and host round trips. Isolate one level, picked by edge count rather than by frontier size. Day 9 is the general form of this mistake.
You assume the biggest frontier is the slowest level. Level 5 holds 1,677 vertices to level 4's ten, and ten times the edges, and level 4 is the slower kernel under the thread mapping, 1.3724 ms against 1.0068, because one of those ten vertices has 1,255 neighbours.
You put the guard around the shuffle. The warp mapping reduces with a
warp shuffle, __shfl_down_sync under a full
0xffffffff mask, which promises all 32 lanes are there. Wrapping the ladder
in the if that guards the walk breaks that promise, and you get a wrong sum
rather than an error. Day 23 has the rule: guard the
loads, never the collective.
You expect the warp mapping to be free. It launches 32 threads per frontier vertex, so a level of degree-1 vertices carries 31 idle lanes on each and a worse occupancy story than the thread mapping. On this graph the warp mapping won every level measured, but the margin collapsed from 15.58x on the skewed level to 2.41x on the uniform one, and that collapsing margin, not the sign, is why real implementations dispatch per level.
Go deeper
- CUDA C++ Programming Guide, "C++ Language Extensions", for
atomicCASand the__shfl_down_syncmask rules: https://docs.nvidia.com/cuda/cuda-programming-guide/05-appendices/cpp-language-extensions.html (checked 2026-08-29) - CUDA C++ Best Practices Guide 10.2.1, "Coalesced Access to Global Memory": https://docs.nvidia.com/cuda/cuda-c-best-practices-guide/index.html (checked 2026-08-30)
- Programming Massively Parallel Processors, 4th edition, the graph traversal chapter: https://shop.elsevier.com/books/programming-massively-parallel-processors/hwu/978-0-323-91231-0
code/day38-bfs/tools/build_word_graph.py, the rule behind every edge in the shipped graph
Next
Day 39 replaces the hand-written frontier build with CCCL, where day 33's
compaction is one cub::DeviceScan call, so the question becomes which parts
of this program were worth writing. Day 40 is the second capstone, PageRank
over the same CSR layout, where the imbalance you just measured stops being
one level's problem and becomes every iteration's.