code/day38-bfs/bfs.cuThis is the source used by the lesson and its recorded evidence. Compile commands and expected output live in the directory README.
// SPDX-License-Identifier: MIT
//
// Day 38: level-synchronous breadth-first search on a real graph.
//
// The graph is word-graph.txt beside this file. Vertices are 3485 words from
// this project's own research documents; two words share an edge when they
// sit next to each other three or more times. It is not a random graph, and
// that is the point: "the" has 1255 neighbours and the median word has 3.
// A uniform random graph of the same size would give every vertex about 9
// and the lesson would have nothing to measure.
//
// Two mappings from a frontier onto threads are compared:
// thread per vertex one thread reads one vertex's whole neighbour list
// warp per vertex 32 lanes split one vertex's neighbour list
// The first is what everybody writes first. On the level that holds "the" it
// runs one lane for 1255 iterations while the other 31 wait.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o bfs bfs.cu
// Run: ./bfs from this directory, because it opens word-graph.txt
//
// Verified 2026-08-30 on a Tesla T4 (compute capability 7.5), driver
// 595.84, CUDA 12.6 (V12.6.85). Transcript: evidence/run-2026-08-30.txt
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <utility>
#include <vector>
#include <cuda_runtime.h>
// The one error macro. This file is standalone, so it carries its own
// verbatim copy. `err_` carries a trailing underscore so it cannot collide
// with a variable at the call site, and the do/while makes the macro one
// statement so it survives a braceless `if`.
#define CUDA_CHECK(call) \
do { \
cudaError_t err_ = (call); \
if (err_ != cudaSuccess) { \
std::fprintf(stderr, "CUDA error %s:%d: %s: %s\n", __FILE__, \
__LINE__, #call, cudaGetErrorString(err_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
// kSourceWord is one of the four words whose BFS needs ten levels; every
// other source in this graph needs fewer. Ten levels means ten kernel
// launches, which is the number this lesson is about.
//
// kWalkRepeats lifts the isolated neighbour walk above launch overhead. One
// pass over the widest level reads 25235 entries, which finishes in less
// time than a launch takes, so a single pass would time the launch instead
// of the walk. Day 15's probe repeats its reads for the same reason.
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarpSize = 32; // 32 on every GPU this course targets
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr int kWalkRepeats = 64;
constexpr int kWalkMinEdges = 1000; // smaller levels are not worth timing
constexpr int kMaxWordLen = 64;
constexpr int kMaxLevels = 256;
constexpr const char* kGraphPath = "word-graph.txt";
constexpr const char* kSourceWord = "await";
// The warp kernels give one warp to one frontier vertex, so a block has to
// hold a whole number of warps. That is true at 32 threads per block as well
// as at 256, which is why it can be asserted rather than assumed.
static_assert(kThreadsPerBlock % kWarpSize == 0,
"block size must be a whole number of warps");
// CSR, the layout day 37 builds for sparse matrices. rowOffsets[v] through
// rowOffsets[v + 1] is the slice of colIndices holding v's neighbours, so a
// degree is one subtraction and a neighbour list is contiguous. The words
// are carried only so the program can name the vertex that hurts.
struct Graph {
int vertexCount;
int edgeCount;
std::vector<char> words; // vertexCount blocks of kMaxWordLen chars
std::vector<int> rowOffsets; // vertexCount + 1 entries
std::vector<int> colIndices; // 2 * edgeCount entries, both directions
};
// Every device pointer the search touches, in one bag, so runBfs takes one
// argument instead of six and cannot be handed a mismatched set.
struct DeviceGraph {
int* rowOffsets;
int* colIndices;
int* dist;
int* frontier;
int* nextFrontier;
unsigned int* nextSize;
int vertexCount;
};
static const char* wordOf(const Graph& g, int v) {
return &g.words[static_cast<size_t>(v) * kMaxWordLen];
}
static int degreeOf(const Graph& g, int v) {
return g.rowOffsets[v + 1] - g.rowOffsets[v];
}
// Reads the next line that is neither a comment nor blank. The graph file
// carries its provenance in '#' lines at the top, so the loader skips them
// rather than assuming a fixed header size.
static bool nextDataLine(std::FILE* file, char* line, int cap) {
while (std::fgets(line, cap, file) != nullptr) {
if (line[0] != '#' && line[0] != '\n' && line[0] != '\r') {
return true;
}
}
return false;
}
// Reads the edge list, then builds CSR with a counting sort: count each
// vertex's degree into rowOffsets[v + 1], prefix sum the counts in place,
// then scatter every edge twice through a cursor. That is day 32's scan and
// day 33's scatter, run once on the host, because the graph is read once.
//
// Returns false and names the first malformed line. A graph file with one
// bad edge is a wrong answer, not a warning.
static bool loadGraph(const char* path, Graph* g) {
std::FILE* file = std::fopen(path, "r");
if (file == nullptr) {
std::fprintf(stderr,
"cannot open %s: run this program from the directory "
"that holds it\n",
path);
return false;
}
char line[512];
bool ok = nextDataLine(file, line, sizeof(line)) &&
std::sscanf(line, "%d %d", &g->vertexCount, &g->edgeCount) == 2 &&
g->vertexCount > 0 && g->edgeCount > 0;
if (!ok) {
std::fprintf(stderr, "%s: bad header line\n", path);
std::fclose(file);
return false;
}
g->words.assign(static_cast<size_t>(g->vertexCount) * kMaxWordLen, '\0');
for (int v = 0; v < g->vertexCount && ok; ++v) {
ok = nextDataLine(file, line, sizeof(line));
if (!ok) {
std::fprintf(stderr, "%s: ran out of labels at %d\n", path, v);
break;
}
const size_t len = std::strcspn(line, " \t\r\n");
if (len == 0 || len >= kMaxWordLen) {
std::fprintf(stderr, "%s: label %d is empty or too long\n", path,
v);
ok = false;
break;
}
std::memcpy(&g->words[static_cast<size_t>(v) * kMaxWordLen], line, len);
}
std::vector<int> edgeFrom(g->edgeCount);
std::vector<int> edgeTo(g->edgeCount);
for (int e = 0; e < g->edgeCount && ok; ++e) {
int u = 0;
int v = 0;
ok = nextDataLine(file, line, sizeof(line)) &&
std::sscanf(line, "%d %d", &u, &v) == 2 && u >= 0 && v >= 0 &&
u < g->vertexCount && v < g->vertexCount && u != v;
if (!ok) {
std::fprintf(stderr, "%s: bad edge at index %d\n", path, e);
break;
}
edgeFrom[e] = u;
edgeTo[e] = v;
}
std::fclose(file);
if (!ok) {
return false;
}
g->rowOffsets.assign(g->vertexCount + 1, 0);
for (int e = 0; e < g->edgeCount; ++e) {
++g->rowOffsets[edgeFrom[e] + 1];
++g->rowOffsets[edgeTo[e] + 1];
}
for (int v = 0; v < g->vertexCount; ++v) {
g->rowOffsets[v + 1] += g->rowOffsets[v];
}
std::vector<int> cursor(g->rowOffsets.begin(), g->rowOffsets.end() - 1);
g->colIndices.assign(static_cast<size_t>(g->edgeCount) * 2, 0);
for (int e = 0; e < g->edgeCount; ++e) {
const int u = edgeFrom[e];
const int v = edgeTo[e];
g->colIndices[cursor[u]++] = v;
g->colIndices[cursor[v]++] = u;
}
return true;
}
static int findVertex(const Graph& g, const char* word) {
for (int v = 0; v < g.vertexCount; ++v) {
if (std::strcmp(wordOf(g, v), word) == 0) {
return v;
}
}
return -1;
}
// Textbook queue BFS, the reference every GPU result is checked against. The
// queue comes out in level order, so levelStart slices it into exactly the
// frontiers the GPU builds. The GPU holds the same vertices in a different
// order, which is why only the distances are compared.
static void bfsCpu(const Graph& g, int source, std::vector<int>* dist,
std::vector<int>* order, std::vector<int>* levelStart) {
dist->assign(g.vertexCount, -1);
order->clear();
levelStart->clear();
(*dist)[source] = 0;
order->push_back(source);
levelStart->push_back(0);
size_t head = 0;
while (head < order->size()) {
const size_t levelEnd = order->size();
levelStart->push_back(static_cast<int>(levelEnd));
while (head < levelEnd) {
const int u = (*order)[head++];
for (int e = g.rowOffsets[u]; e < g.rowOffsets[u + 1]; ++e) {
const int v = g.colIndices[e];
if ((*dist)[v] < 0) {
(*dist)[v] = (*dist)[u] + 1;
order->push_back(v);
}
}
}
}
}
// What a walk kernel must return for one vertex. Unsigned addition wraps by
// definition, so host and device agree bit for bit however many repeats run,
// and a walk the compiler shortened comes back wrong instead of fast.
static unsigned int walkChecksum(const Graph& g, int vertex, int repeats) {
unsigned int acc = 0u;
for (int r = 0; r < repeats; ++r) {
for (int e = g.rowOffsets[vertex]; e < g.rowOffsets[vertex + 1]; ++e) {
acc += static_cast<unsigned int>(g.colIndices[e]) ^
static_cast<unsigned int>(r);
}
}
return acc;
}
// Claims v for the next level if nobody else has, and appends it once.
//
// The compare-and-swap is what makes exactly one thread the owner. A plain
// `if (dist[v] < 0) { dist[v] = level + 1; }` produces the same distances,
// because every racing thread writes the same value, and then appends v once
// per racer. The distances stay right and the frontier grows, so the bug
// shows up as a slow program with a correct answer.
__device__ void claimVertex(int v, int level, int* dist, int* nextFrontier,
unsigned int* nextSize) {
if (atomicCAS(&dist[v], -1, level + 1) == -1) {
const unsigned int slot = atomicAdd(nextSize, 1u);
nextFrontier[slot] = v;
}
}
// One thread per frontier vertex. This is the mapping everybody writes
// first, and it is correct.
//
// Memory: the 32 lanes of a warp start at 32 unrelated rows of colIndices,
// so nothing about the read coalesces. Worse, the loop bound is per lane, so
// the warp runs until its widest vertex is finished and the other 31 lanes
// are switched off for the rest of it. Day 22 calls that divergence; here it
// is the entire cost model.
//
// Launch assumption: gridDim.x * blockDim.x >= frontierSize.
// snippet: expand-thread
__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);
}
}
}
// end snippet
// One warp per frontier vertex. The 32 lanes stride through one neighbour
// list together, so the warp runs ceil(degree / 32) steps instead of degree
// steps, and consecutive lanes read consecutive entries of colIndices, which
// is the coalesced pattern day 11 measured.
//
// The cost is 32 times the threads for the same frontier, so a level made of
// degree-1 vertices wastes 31 lanes on every one of them. Day 37 meets the
// same trade on a sparse matrix and calls it the same thing.
//
// Launch assumption: gridDim.x * blockDim.x >= frontierSize * 32.
__global__ void bfsExpandWarpPerVertex(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;
// snippet: expand-warp
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);
}
}
// end snippet
}
// Walks the same neighbour lists as bfsExpandThreadPerVertex and touches
// nothing else: no atomics, no frontier, no distances. It exists so that one
// level can be timed on its own, which the search itself cannot be, because
// a BFS run a second time from the same state has nothing left to discover.
//
// The XOR with the repeat counter is load bearing. Without it the inner sum
// is loop invariant, the compiler is free to compute it once and multiply by
// repeats, and the measurement becomes an empty loop.
__global__ void walkThreadPerVertex(const int* __restrict__ frontier,
int frontierSize,
const int* __restrict__ rowOffsets,
const int* __restrict__ colIndices,
int repeats,
unsigned int* __restrict__ out) {
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];
const int first = rowOffsets[u];
const int last = rowOffsets[u + 1];
unsigned int acc = 0u;
for (int r = 0; r < repeats; ++r) {
for (int e = first; e < last; ++e) {
acc += static_cast<unsigned int>(colIndices[e]) ^
static_cast<unsigned int>(r);
}
}
out[t] = acc;
}
}
// The same walk, one warp per vertex, finished with the shuffle ladder from
// day 23. Addition is commutative, so this returns the value walkChecksum
// predicts even though the lanes visit the neighbours in a different order.
//
// The guard covers the walk and never the shuffle. `w` is the same for all
// 32 lanes of a warp, so the branch is warp uniform and every lane reaches
// every __shfl_down_sync, which is what the full 0xffffffff mask promises.
__global__ void walkWarpPerVertex(const int* __restrict__ frontier,
int frontierSize,
const int* __restrict__ rowOffsets,
const int* __restrict__ colIndices,
int repeats, unsigned int* __restrict__ out) {
const size_t t = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const unsigned int lane = threadIdx.x % kWarpSize;
const size_t w = t / kWarpSize;
const bool live = w < static_cast<size_t>(frontierSize);
unsigned int acc = 0u;
if (live) {
const int u = frontier[w];
const int last = rowOffsets[u + 1];
for (int r = 0; r < repeats; ++r) {
for (int e = rowOffsets[u] + static_cast<int>(lane); e < last;
e += kWarpSize) {
acc += static_cast<unsigned int>(colIndices[e]) ^
static_cast<unsigned int>(r);
}
}
}
for (int offset = kWarpSize / 2; offset > 0; offset /= 2) {
acc += __shfl_down_sync(0xffffffffu, acc, offset);
}
if (live && lane == 0u) {
out[w] = acc;
}
}
// Times a launch with CUDA events and returns the mean milliseconds per run.
//
// This is the one template and the one lambda the course allows before
// module 8. Copy it verbatim; the alternative is six copies of the event
// boilerplate, which is how a warm-up goes missing from one of them.
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
// Warm up this kernel, not just the first kernel in the program. Lazy
// module loading has been the default since CUDA 12.2 on Linux, so the
// first launch of each kernel pays its own load.
for (int i = 0; i < kWarmupRuns; ++i) {
launch();
}
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < kTimedRuns; ++i) {
launch();
}
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError());
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / kTimedRuns;
}
// Picks the mapping and the grid for one level. Split out so the level loop
// below reads as the four things it does and nothing else.
static void launchExpand(const DeviceGraph& d, const int* frontier,
int frontierSize, int* nextFrontier, int level,
bool useWarpMapping) {
if (useWarpMapping) {
const size_t threads = static_cast<size_t>(frontierSize) * kWarpSize;
const int blocks = static_cast<int>((threads + kThreadsPerBlock - 1) /
kThreadsPerBlock);
bfsExpandWarpPerVertex<<<blocks, kThreadsPerBlock>>>(
frontier, frontierSize, d.rowOffsets, d.colIndices, d.dist, level,
nextFrontier, d.nextSize);
} else {
const int blocks =
(frontierSize + kThreadsPerBlock - 1) / kThreadsPerBlock;
bfsExpandThreadPerVertex<<<blocks, kThreadsPerBlock>>>(
frontier, frontierSize, d.rowOffsets, d.colIndices, d.dist, level,
nextFrontier, d.nextSize);
}
CUDA_CHECK(cudaGetLastError());
}
// Runs one whole search and returns the number of levels, writing each
// level's frontier size into levelSizes when that pointer is not null.
//
// Two properties of this loop are the lesson. One launch per level, so the
// launch count is a property of the graph and not of the GPU. And one copy
// back to the host per level, because the host cannot size the next grid, or
// know it has finished, without the frontier count. That copy is a
// synchronisation point sitting in the middle of the algorithm.
static int runBfs(const DeviceGraph& d, int source, bool useWarpMapping,
int* levelSizes) {
// 0xff in every byte is -1 in every int, which is the not-reached mark.
const size_t distBytes = static_cast<size_t>(d.vertexCount) * sizeof(int);
const int zero = 0;
CUDA_CHECK(cudaMemset(d.dist, 0xff, distBytes));
CUDA_CHECK(cudaMemcpy(&d.dist[source], &zero, sizeof(int),
cudaMemcpyHostToDevice));
CUDA_CHECK(
cudaMemcpy(d.frontier, &source, sizeof(int), cudaMemcpyHostToDevice));
int* frontier = d.frontier;
int* nextFrontier = d.nextFrontier;
int frontierSize = 1;
int levels = 0;
// snippet: level-loop
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;
}
// end snippet
return levels;
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
std::printf("GPU: %s (compute capability %d.%d), %d SMs\n", prop.name,
prop.major, prop.minor, prop.multiProcessorCount);
Graph g;
if (!loadGraph(kGraphPath, &g)) {
return EXIT_FAILURE;
}
const int source = findVertex(g, kSourceWord);
if (source < 0) {
std::fprintf(stderr, "%s holds no vertex named %s\n", kGraphPath,
kSourceWord);
return EXIT_FAILURE;
}
// Degree statistics first, because they are why this graph was chosen
// and everything below is a consequence of them.
std::vector<int> degrees(g.vertexCount);
int hub = 0;
int degreeOne = 0;
for (int v = 0; v < g.vertexCount; ++v) {
degrees[v] = degreeOf(g, v);
if (degrees[v] == 1) {
++degreeOne;
}
if (degrees[v] > degrees[hub]) {
hub = v;
}
}
std::vector<int> sortedDegrees(degrees);
std::sort(sortedDegrees.begin(), sortedDegrees.end());
std::printf("\nGraph: %s\n", kGraphPath);
std::printf(" %d vertices, %d undirected edges, %zu csr entries\n",
g.vertexCount, g.edgeCount, g.colIndices.size());
std::printf(" max degree %d (%s), median %d, %d vertices of degree 1\n",
degrees[hub], wordOf(g, hub),
sortedDegrees[sortedDegrees.size() / 2], degreeOne);
std::vector<int> cpuDist;
std::vector<int> order;
std::vector<int> levelStart;
bfsCpu(g, source, &cpuDist, &order, &levelStart);
const int levels = static_cast<int>(levelStart.size()) - 1;
const int reached = static_cast<int>(order.size());
// Per level: the real work is the sum of the frontier's degrees, and the
// work a thread-per-vertex launch issues is 32 lane slots per warp per
// iteration, with a warp iterating as many times as its widest vertex.
// Both numbers come from the graph, not from the clock.
std::printf("\nBFS from \"%s\" (vertex %d)\n", kSourceWord, source);
std::printf("%6s %9s %8s %9s %-12s %6s %11s\n", "level", "frontier",
"edges", "max deg", "busiest", "warps", "lane slots");
std::printf("%6s %9s %8s %9s %-12s %6s %11s\n", "-----", "--------",
"-----", "-------", "-------", "-----", "----------");
long long totalEdges = 0;
long long totalSlots = 0;
for (int lv = 0; lv < levels; ++lv) {
const int begin = levelStart[lv];
const int end = levelStart[lv + 1];
int edges = 0;
int maxDegree = 0;
int busiest = order[begin];
int warpMax = 0;
long long slots = 0;
for (int k = begin; k < end; ++k) {
const int v = order[k];
edges += degrees[v];
if (degrees[v] > maxDegree) {
maxDegree = degrees[v];
busiest = v;
}
warpMax = std::max(warpMax, degrees[v]);
if ((k - begin) % kWarpSize == kWarpSize - 1 || k == end - 1) {
slots += static_cast<long long>(kWarpSize) * warpMax;
warpMax = 0;
}
}
totalEdges += edges;
totalSlots += slots;
std::printf("%6d %9d %8d %9d %-12s %6d %11lld\n", lv, end - begin,
edges, maxDegree, wordOf(g, busiest),
(end - begin + kWarpSize - 1) / kWarpSize, slots);
}
std::printf(" reached %d of %d vertices in %d levels, so %d launches\n",
reached, g.vertexCount, levels, levels);
std::printf(" %lld edges to walk, %lld lane slots issued to walk them\n",
totalEdges, totalSlots);
// Device side. Seven allocations and one free block at the end. Every
// failure below records itself instead of returning, so nothing leaks on
// the path where the answer is wrong.
const size_t offsetBytes =
static_cast<size_t>(g.vertexCount + 1) * sizeof(int);
const size_t colBytes = g.colIndices.size() * sizeof(int);
const size_t vertexBytes = static_cast<size_t>(g.vertexCount) * sizeof(int);
DeviceGraph d;
d.rowOffsets = nullptr;
d.colIndices = nullptr;
d.dist = nullptr;
d.frontier = nullptr;
d.nextFrontier = nullptr;
d.nextSize = nullptr;
d.vertexCount = g.vertexCount;
unsigned int* d_walkOut = nullptr;
CUDA_CHECK(cudaMalloc(&d.rowOffsets, offsetBytes));
CUDA_CHECK(cudaMalloc(&d.colIndices, colBytes));
CUDA_CHECK(cudaMalloc(&d.dist, vertexBytes));
CUDA_CHECK(cudaMalloc(&d.frontier, vertexBytes));
CUDA_CHECK(cudaMalloc(&d.nextFrontier, vertexBytes));
CUDA_CHECK(cudaMalloc(&d.nextSize, sizeof(unsigned int)));
CUDA_CHECK(cudaMalloc(&d_walkOut, vertexBytes));
CUDA_CHECK(cudaMemcpy(d.rowOffsets, g.rowOffsets.data(), offsetBytes,
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d.colIndices, g.colIndices.data(), colBytes,
cudaMemcpyHostToDevice));
int failures = 0;
std::vector<int> gpuDist(g.vertexCount);
std::vector<int> levelSizes(kMaxLevels, 0);
std::printf("\nCorrectness\n");
for (int variant = 0; variant < 2; ++variant) {
const bool useWarp = variant == 1;
const char* name = useWarp ? "warp per vertex" : "thread per vertex";
const int gpuLevels = runBfs(d, source, useWarp, levelSizes.data());
CUDA_CHECK(cudaMemcpy(gpuDist.data(), d.dist, vertexBytes,
cudaMemcpyDeviceToHost));
// Distances are integers, so this comparison is exact. A tolerance
// would hide the only class of bug it can catch.
int firstBad = -1;
for (int v = 0; v < g.vertexCount && firstBad < 0; ++v) {
if (gpuDist[v] != cpuDist[v]) {
firstBad = v;
}
}
if (firstBad >= 0) {
std::fprintf(stderr,
"%s: distance wrong at vertex %d (%s): got %d, "
"want %d\n",
name, firstBad, wordOf(g, firstBad), gpuDist[firstBad],
cpuDist[firstBad]);
++failures;
} else {
std::printf(" %-18s distances match on all %d vertices\n", name,
g.vertexCount);
}
if (gpuLevels != levels) {
std::fprintf(stderr, "%s: %d levels, want %d\n", name, gpuLevels,
levels);
++failures;
}
// The frontier total is the check that catches a duplicate append.
// The distances survive that bug; this does not.
long long enqueued = 0;
for (int lv = 0; lv < gpuLevels && lv < kMaxLevels; ++lv) {
enqueued += levelSizes[lv];
}
if (enqueued != reached) {
std::fprintf(stderr,
"%s: frontiers held %lld vertices, want %d. Some "
"vertex was enqueued more than once.\n",
name, enqueued, reached);
++failures;
} else {
std::printf(" %-18s frontiers held %d vertices, no duplicates\n",
name, reached);
}
}
// Whole search, both mappings. This measurement includes one distance
// reset and one host round trip per level, because a level-synchronous
// BFS cannot run without them. CUDA-CODE-STYLE.md says never to
// synchronise inside a timed loop; that rule is about not adding a sync
// a kernel benchmark does not need, and here the sync is the algorithm.
std::printf(
"\nWhole BFS, %d timed runs, reset and %d host round trips included\n",
kTimedRuns, levels);
std::printf("%-20s %12s\n", "mapping", "ms per bfs");
std::printf("%-20s %12s\n", "-------------------", "----------");
for (int variant = 0; variant < 2; ++variant) {
const bool useWarp = variant == 1;
const float ms =
timeKernel([&] { runBfs(d, source, useWarp, nullptr); });
std::printf("%-20s %12.4f\n",
useWarp ? "warp per vertex" : "thread per vertex", ms);
if (!(ms > 0.0f)) {
std::fprintf(stderr, "whole bfs time is not positive\n");
++failures;
}
}
// One level on its own. Same frontier, same entries read, no atomics, so
// the only difference left between the two rows is which thread reads
// which neighbour.
std::printf("\nOne level in isolation, neighbour walk repeated %d times\n",
kWalkRepeats);
std::printf("%6s %9s %8s %9s %11s %11s %9s\n", "level", "frontier", "edges",
"max deg", "thread ms", "warp ms", "speedup");
std::printf("%6s %9s %8s %9s %11s %11s %9s\n", "-----", "--------", "-----",
"-------", "---------", "-------", "-------");
std::vector<unsigned int> want(g.vertexCount);
std::vector<unsigned int> got(g.vertexCount);
for (int lv = 0; lv < levels; ++lv) {
const int begin = levelStart[lv];
const int size = levelStart[lv + 1] - begin;
int edges = 0;
int maxDegree = 0;
for (int k = 0; k < size; ++k) {
edges += degrees[order[begin + k]];
maxDegree = std::max(maxDegree, degrees[order[begin + k]]);
}
if (edges < kWalkMinEdges) {
continue;
}
const size_t frontierBytes = static_cast<size_t>(size) * sizeof(int);
const size_t outBytes =
static_cast<size_t>(size) * sizeof(unsigned int);
CUDA_CHECK(cudaMemcpy(d.frontier, &order[begin], frontierBytes,
cudaMemcpyHostToDevice));
for (int k = 0; k < size; ++k) {
want[k] = walkChecksum(g, order[begin + k], kWalkRepeats);
}
const int threadBlocks =
(size + kThreadsPerBlock - 1) / kThreadsPerBlock;
const size_t warpThreads = static_cast<size_t>(size) * kWarpSize;
const int warpBlocks = static_cast<int>(
(warpThreads + kThreadsPerBlock - 1) / kThreadsPerBlock);
const float threadMs = timeKernel([&] {
walkThreadPerVertex<<<threadBlocks, kThreadsPerBlock>>>(
d.frontier, size, d.rowOffsets, d.colIndices, kWalkRepeats,
d_walkOut);
});
CUDA_CHECK(cudaMemcpy(got.data(), d_walkOut, outBytes,
cudaMemcpyDeviceToHost));
int bad = -1;
for (int k = 0; k < size && bad < 0; ++k) {
if (got[k] != want[k]) {
bad = k;
}
}
const float warpMs = timeKernel([&] {
walkWarpPerVertex<<<warpBlocks, kThreadsPerBlock>>>(
d.frontier, size, d.rowOffsets, d.colIndices, kWalkRepeats,
d_walkOut);
});
CUDA_CHECK(cudaMemcpy(got.data(), d_walkOut, outBytes,
cudaMemcpyDeviceToHost));
for (int k = 0; k < size && bad < 0; ++k) {
if (got[k] != want[k]) {
bad = k;
}
}
if (bad >= 0) {
std::fprintf(stderr,
"level %d: walk wrong at slot %d, vertex %d (%s): "
"got %u, want %u\n",
lv, bad, order[begin + bad],
wordOf(g, order[begin + bad]), got[bad], want[bad]);
++failures;
}
if (!(threadMs > 0.0f) || !(warpMs > 0.0f)) {
std::fprintf(stderr, "level %d: a walk time is not positive\n", lv);
++failures;
}
std::printf("%6d %9d %8d %9d %11.4f %11.4f %8.2fx\n", lv, size, edges,
maxDegree, threadMs, warpMs, threadMs / warpMs);
}
// One free block, reached on every path that got this far.
CUDA_CHECK(cudaFree(d.rowOffsets));
CUDA_CHECK(cudaFree(d.colIndices));
CUDA_CHECK(cudaFree(d.dist));
CUDA_CHECK(cudaFree(d.frontier));
CUDA_CHECK(cudaFree(d.nextFrontier));
CUDA_CHECK(cudaFree(d.nextSize));
CUDA_CHECK(cudaFree(d_walkOut));
if (failures > 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
std::printf("\nall checks passed\n");
return EXIT_SUCCESS;
}