PageRank on a real graph in CUDA
Implement PageRank twice on a graph of 28,042 papers: first with CUDA kernels, then with Thrust and CUB. Validate both versions against one CPU reference, measure their memory traffic and run time, and justify which version you would keep.
The graph contains 2,204 papers that cite no other paper in the set. Such a paper is a dangling node, and its rank must return to the whole graph on every iteration.
Retrieve the count-and-scan pattern from days 29 and 31. You can verify both
graph counts in code/day40-pagerank/data/graph.txt with awk before the GPU
does any work.
Before reading the implementation, predict what happens if you omit that redistribution. Decide whether the program fails, whether the ranks still converge, whether their sum stays at 1, and whether the top-twenty order changes.
One iteration, three shapes of work
PageRank assigns more rank to a paper when high-ranked papers cite it. The
damping factor d is the share of rank passed through citations, and n
is the node count:
x'[v] = (1 - d)/n + d * ( sum over citers u of v of x[u]/outdeg[u]
+ danglingMass/n )
Brin and Page's PageRank paper
uses d = 0.85. This capstone fixes the same value so every run can use the
same reference.
One iteration combines sparse matrix-vector multiplication with two reductions.
The sum over citers is the CSR SpMV from day 37. The input stores out-edges, which record papers cited by each paper, while the formula needs in-edges, which record papers that cite each paper.
The host transposes the graph once with a counting sort before copying it to the GPU. That transpose uses day 29's histogram, a host prefix sum based on day 31, and a scatter.
The dangling mass uses the full-array reduction from days 24 and 25. The convergence test reduces the L1 distance, the sum of absolute differences between two successive vectors.
Now test the prediction with an invariant, a property that should hold after
each iteration. Sum the update over all v: the constant term gives 1 - d,
the SpMV term gives d times the non-dangling rank, and the dangling term
restores the rest.
The rank vector should therefore sum to 1 after every pass. This check can find a bug before a reference implementation exists.
If D is the dangling rank, omitting its redistribution changes the next
total to 1 - d * D. Keep that expression for the exercise, where you will
derive the fixed-point total.
The row mapping affects whether the col_idx reads use
memory coalescing. Predict which diagram band
will issue fewer memory transactions per lane on this graph.
Keep the dangling mass on the device
One implementation copies the dangling mass to the host and passes the value back as a kernel argument:
float mass = 0.0f;
cudaMemcpy(&mass, d_mass, sizeof(float), cudaMemcpyDeviceToHost);
combineRanks<<<blocks, threads>>>(d_y, mass, d_xNext, n, damping);
This copy synchronizes the device on every iteration. The
CUDA Runtime API synchronization rules
state that a device-to-host cudaMemcpy returns only after the copy completes.
Day 9 measures this wait.
Pass a device pointer instead and let every thread read danglingMass[0].
The 32 lanes of a warp request the same
address, so the memory system broadcasts one value to the warp.
The convergence test still copies one residual to the host because the host decides whether to stop. The timed case disables that test and runs fifty iterations, so both implementations time the same work.
Implement the loop twice
Full program in code/day40-pagerank/pagerank.cu,
with the build line, data provenance, and output guide in the README beside
it.
Use three controls when comparing the implementations.
Both implementations start from the same vector and run the same fifty iterations. CUDA events time warmed kernels, and all copies stay outside the timed region.
Every correctness case uses one double-precision CPU reference. The CPU and GPU versions run 100 iterations so a different stopping point cannot change the comparison.
Each implementation prints the bytes it moves. Effective bandwidth is bytes divided by time, so compare both traffic and milliseconds as in day 11.
Here is one iteration with the course's own kernels:
static void iterateHand(const DeviceGraph& g, const Work& w, const float* x,
float* xNext, const PrParams& p) {
const int blocks = gridFor(g.n);
scaleByOutDegree<<<blocks, kThreadsPerBlock>>>(x, g.invOut, w.contrib, g.n);
if (p.warpPerRow) {
spmvCsrWarp<<<gridForWarps(g.n), kThreadsPerBlock>>>(
g.rowPtr, g.colIdx, w.contrib, w.y, g.n);
} else {
spmvCsrScalar<<<blocks, kThreadsPerBlock>>>(g.rowPtr, g.colIdx,
w.contrib, w.y, g.n);
}
// Two passes: kReduceBlocks writes, then one. That is the budget this
// stage is graded on, and it is what an atomicAdd per thread misses.
sumDanglingPartial<<<kReduceBlocks, kThreadsPerBlock>>>(x, g.dangling,
w.partials, g.n);
sumPartials<<<1, kThreadsPerBlock>>>(w.partials, w.scalars, kReduceBlocks);
combineRanks<<<blocks, kThreadsPerBlock>>>(w.y, w.scalars, xNext, g.n,
p.damping);
}
The library version uses CCCL and contains no custom kernel:
static void iterateLibrary(const DeviceGraph& g, Work& w, const float* x,
float* xNext, const PrParams& p) {
const float invN = 1.0f / static_cast<float>(g.n);
thrust::transform(thrust::device, x, x + g.n, g.invOut, w.contrib,
thrust::multiplies<float>());
thrust::gather(thrust::device, g.colIdx, g.colIdx + g.nnz, w.contrib,
w.gathered);
CUDA_CHECK(cub::DeviceSegmentedReduce::Sum(
w.cubTemp, w.cubSegBytes, w.gathered, w.y, static_cast<int>(g.n),
g.rowPtr, g.rowPtr + 1));
thrust::transform(thrust::device, x, x + g.n, g.dangling, w.diff,
MaskedValue());
CUDA_CHECK(cub::DeviceReduce::Sum(w.cubTemp, w.cubRedBytes, w.diff,
w.scalars, g.n));
thrust::transform(thrust::device, w.y, w.y + g.n, xNext,
Combine{w.scalars, p.damping, invN});
}
Each driver has seventeen non-comment lines. The library path removes 105 non-comment lines of device code and adds 20 lines of function objects.
Predict its memory cost before looking at the results. Thrust
cannot fuse thrust::gather with the segmented reduction, so it writes
314,010 gathered floats to global memory for
CUB to read back.
The library path also materializes the masked rank vector. The custom path uses kernel fusion to combine each pair of operations in one kernel.
One hand-written iteration touches 3,437,470 bytes; one library iteration touches 6,173,886, which is 1.80 times as many.
Before measuring, predict which path is faster and whether its speed ratio will match the 1.80 traffic ratio. State which result would justify keeping 105 extra lines of device code.
Note. CUB 2.5 includes
cub::DeviceSpmv::CsrMV, but CCCL 3.0 removes it and directs users to cuSPARSE in the CCCL 3.0 migration guide. This capstone does not use an API removed by the next library version; cuSPARSE appears on day 82.
Validate before timing
Run all six correctness cases before reading the timing table. The ranks must
sum to 1 within 4 * FLT_EPSILON * sqrt(n), and Spearman's rho over the CPU
reference's top 100 papers must reach 0.999.
Spearman's rho measures agreement between two rank orders. It is a better gate than float equality because reduction order can change the last digits without changing which papers rank highest.
Then grade speed against the copy bandwidth measured in the same process. The fraction is portable across GPUs in a way that an absolute time is not.
| Tier | Requirement | 7.5 | 8.0, 8.9, 12.0 | 9.0, 10.0 |
|---|---|---|---|---|
| pass | all six cases correct | |||
| bronze | beats the thread-per-row version on the timed case | |||
| silver | write budgets met, and this fraction of the DRAM copy row | 0.22 | 0.20 | 0.15 |
| gold | the same, plus coefficient of variation under 5 percent | 0.38 | 0.35 | 0.28 |
These fractions are lower than capstone 1's
because the gather through col_idx can place each lane in a separate
32-byte sector. Most bytes in each fetched
sector then go unused, as day 11 predicts.
The 7.5 column is calibrated by the run on this page. The reference kernels reach 95.7 GB/s over 50 iterations against the 244.8 GB/s DRAM copy row, a fraction of 0.391, just above the 0.38 gold line.
This calibration follows research/CAPSTONE-SPECS.md section 0.4: a tier must
have a measured passing result before it can grade a learner.
The other columns remain uncalibrated until those cards run. The two reduction
stages may write at most gridDim.x + 1 floats to global memory in total;
this budget follows from the launch shape rather than a measured threshold.
Submit the top-twenty table and one paragraph that names the next change and its expected effect, with a number. If your GPU does not reach the tier, name the stage with the gap, cite the matching table row, and state what you would change first.
Report the GPU and test context with the measured result.
Before opening the Results section, write three predictions: which row mapping wins, how much slower the extra 1.80 times traffic makes the library path, and whether removing the dangling reduction changes the top twenty. After the run, explain each gap between prediction and measurement.
Results
Re-verified on a Tesla T4 with driver 580.173.02 and CUDA 13.0 (V13.0.88) on
2026-09-02. Every synthetic and real-graph correctness gate reproduced,
including convergence counts and top-100 rank order. End-to-end timing moved:
the hand-written implementation remained faster, but its advantage narrowed
from 5.4x to 3.7x. The original CUDA 12.6 table below and CUDA 13 transcript
both remain in evidence.
GPU: Tesla T4 (compute capability 7.5)
40 SMs, 4194304 bytes of L2, 1024 threads per block
loaded 28042 nodes and 314010 edges from data/graph.txt
correctness, 100 fixed iterations on both sides
tiny thread per row n=6 nnz=9 err/tol= 0.000 sum-1=+2.98e-08 pass
tiny warp per row n=6 nnz=9 err/tol= 0.000 sum-1=+0.00e+00 pass
tiny Thrust + CUB n=6 nnz=9 err/tol= 0.000 sum-1=+1.12e-07 pass
single thread per row n=1 nnz=0 err/tol= 0.000 sum-1=+0.00e+00 pass
single warp per row n=1 nnz=0 err/tol= 0.000 sum-1=+0.00e+00 pass
single Thrust + CUB n=1 nnz=0 err/tol= 0.000 sum-1=+0.00e+00 pass
nodangling thread per row n=4096 nnz=8192 err/tol= 0.000 sum-1=+0.00e+00 pass
nodangling warp per row n=4096 nnz=8192 err/tol= 0.000 sum-1=+0.00e+00 pass
nodangling Thrust + CUB n=4096 nnz=8192 err/tol= 0.000 sum-1=+0.00e+00 pass
dangling thread per row n=4096 nnz=7392 err/tol= 0.000 sum-1=+6.32e-09 pass
dangling warp per row n=4096 nnz=7392 err/tol= 0.000 sum-1=+6.32e-09 pass
dangling Thrust + CUB n=4096 nnz=7392 err/tol= 0.000 sum-1=+3.73e-08 pass
star thread per row n=10001 nnz=10000 err/tol= 0.216 sum-1=-2.81e-05 pass
star warp per row n=10001 nnz=10000 err/tol= 0.013 sum-1=-3.76e-07 pass
star Thrust + CUB n=10001 nnz=10000 err/tol= 0.001 sum-1=-2.13e-07 pass
chain thread per row n=8192 nnz=8191 err/tol= 0.000 sum-1=-3.16e-07 pass
chain warp per row n=8192 nnz=8191 err/tol= 0.000 sum-1=-3.16e-07 pass
chain Thrust + CUB n=8192 nnz=8191 err/tol= 0.000 sum-1=-3.16e-07 pass
case real 28042 nodes, 314010 edges
CPU reference, double, converged in 109 iterations
our kernels, float, L1 residual under 1e-06 in 53 iterations
Thrust + CUB, float, the same tolerance in 53 iterations
sum of ranks minus 1: ours +2.097e-07, library +1.966e-07 (allowed 7.985e-05) pass
Spearman rho over the reference's top 100: ours 1.000000, library 1.000000, need 0.999 pass
top 20 by PageRank (work id, year, title)
rank score in-deg out-deg paper
1 0.0075687 80 1 W2096733429 2004 Global register allocation at link time
2 0.0074066 348 2 W2002257715 1990 A set of level 3 basic linear algebra subprograms
3 0.0065285 6 8 W2012182306 1992 Experience with a software-defined machine architecture
4 0.0054581 62 1 W2152706734 1990 Supporting shared data structures on distributed memory architectures
5 0.0049975 54 1 W2085809150 1990 Run-time scheduling and execution of loops on message passing machines
6 0.0045897 68 3 W1964714157 2004 Interprocedural dependence analysis and parallelization
7 0.0037279 86 0 W2170830522 1998 Using cache memory to reduce processor-memory traffic
8 0.0036198 39 1 W2006384005 1990 Algorithm 679: A set of level 3 basic linear algebra subprograms: model implementation and test programs
9 0.0032463 538 6 W2119609467 1991 A data locality optimizing algorithm
10 0.0031577 2 0 W1878455722 2020 Advanced Architecture Computers
11 0.0030382 232 7 W2157758640 2004 Software pipelining
12 0.0027620 45 10 W2121619519 2004 Efficient instruction scheduling for a pipelined architecture
13 0.0026831 10 2 W2152835688 1991 The Scalable Processor Architecture (SPARC)
14 0.0026726 23 1 W2048712467 2004 Register windows vs. register allocation
15 0.0025360 112 2 W2152450729 1998 An evaluation of directory schemes for cache coherence
16 0.0024872 104 1 W2102890180 1991 Constant propagation with conditional branches
17 0.0024015 212 8 W2129192659 1991 Limits of instruction-level parallelism
18 0.0022408 108 10 W2102582914 2004 Improving register allocation for subscripted variables
19 0.0022309 364 3 W2098220211 1991 The cache performance and optimizations of blocked algorithms
20 0.0020716 232 0 W2163820265 1998 Lockup-free instruction fetch/prefetch cache organization
timed on the shipped graph, mean of 10 runs after 3 warm-ups, copies not included
one iteration moves 3437470 bytes with our kernels and 6173886 with the library
stage ms GB/s x DRAM
copy, 128 MiB (DRAM) 0.5483 244.8 1.00
copy, one iteration's bytes 0.0063 542.3 2.22
scale by out-degree 0.0028 118.8 0.49
spmv, thread per row 0.0332 82.4 0.34
spmv, warp per row 0.0273 100.2 0.41
spmv, gather + CUB segmented 0.1524 34.4 0.14
dangling mass, tree reduction 0.0082 17.1 0.07
residual, tree reduction 0.0051 43.6 0.18
residual, one atomic per thread 0.0583 3.8 0.02
50 iterations end to end, no convergence test on either side
implementation ms ms/iter GB/s
our kernels 1.7963 0.0359 95.7
Thrust + CUB 9.6324 0.1926 32.0
our kernels take 0.19x the library's time and move 0.56x its bytes
L2 on this card is 4194304 bytes, and one iteration touches 3437470
every case passed
All three implementations agree on six synthetic cases: a tiny graph, a single isolated node, graphs with and without dangling nodes, a star, and a chain. The dangling case checks whether the implementation redistributes rank that has no outgoing edge.
The sum-1 column reports that invariant. After 100 fixed iterations, the
values are at the 1e-8 level for the custom kernels and 1e-7 for the Thrust
and CUB version.
Float addition is not associative, so the library's reduction order produces a different rounded result. Both errors remain below the required tolerance; day 26 explains why reduction order can vary.
The custom implementation wins in both toolkits, but the margin changes with the version. Under CUDA 12.6, fifty iterations took 1.7963 ms against 9.6324 ms for Thrust and CUB, a 5.4x gap.
Under CUDA 13.0, the same work took 3.0900 and 11.4025 ms, a 3.7x gap. The custom path still moved 0.56 times as many bytes.
The per-stage table shows why the 1.80 traffic ratio does not predict the full
time ratio. spmv, gather + CUB segmented took 0.2587 ms, 5.4 times the fused
warp SpMV's 0.0475 ms.
The unfused gather adds a full global-memory round trip. It also loses the L2 residency kept by the fused kernel.
The graph has 28,042 nodes and 314,010 edges from the supplied paper data. Its skewed degree distribution exposes the row imbalance that a uniform random graph would hide.
Run it yourself
The graph and titles are 3.9 MB of text stored in the repository. The program does not download data at run time.
The command below targets sm_75. Change that target to your GPU's compute
capability when needed.
nvcc -std=c++17 -O3 -arch=sm_75 -o pagerank pagerank.cu
./pagerank
Compiler Explorer cannot run the full case. The double-precision CPU reference alone iterates 28,042 nodes to a tight tolerance, and two timed sweeps of fifty iterations follow it, which exceeds its 20-second limit.
If Compiler Explorer is your only option, cut the case list to tiny,
star, and chain. This checks correctness but cannot validate a speed tier.
For a local CUDA setup, use day 3.
Exercise
Delete the dangling stage. Take out the two reduction launches and pass a zero where the mass goes. Predict where the sum of the ranks settles before you run it, then check the top twenty against the correct version.
Time: 30 to 45 minutes. Submit: the predicted total, the measured total, whether any paper moved in the top twenty, and one sentence on why.
Check: the program prints sum-1 for every case and implementation, so
the broken version grades itself. tiny fails first at six nodes, which is
small enough to calculate by hand.
nodangling still passes because it has no dangling node. Compare those two
cases to isolate the missing redistribution term.
The per-element comparison also fails. It reports the largest error as a multiple of the tolerance, so the output shows the size of the error.
Hint 1
Add up the update over every node v and see what the total becomes. The
constant term contributes the same amount n times.
The SpMV term contributes each non-dangling node's rank exactly once because a node's mass is split across its out-edges and then reassembled. What remains?
Hint 2
You should reach S' = (1 - d) + d * (S - D), where S is the total and D
is the mass on the dangling nodes. Solve that for S at the fixed point and
note that d / (1 - d) is a large multiplier when d is 0.85.
For the second half, write both fixed points in matrix form. Both are
x = k * 1 + d * A * x for the same A, differing only in the constant k.
How does that relate the two solutions?
Solution
Summing the broken update gives S' = (1 - d) + d * (S - D). At the fixed
point S = 1 - D * d / (1 - d), and with d at 0.85 that multiplier is
17/3, so the total falls short by nearly six times the dangling mass.
The ranking does not move.
Both fixed points solve (I - d A) x = k * 1 for the same matrix A and
differ only in the constant. The correct one carries (1 - d)/n + d D / n
and the broken one carries (1 - d)/n.
Two linear systems with the same matrix and proportional right-hand sides have proportional solutions. The broken vector is the correct vector times a constant below 1, so every rank shrinks by the same factor and no pair swaps.
The printed ranking hides the bug because it stays the same. Only the total changes.
When a computation conserves a quantity, print it. This check uses a reduction you have already written and works even when the graph has no known ranking.
Pitfalls
If the ranks sum to less than 1 while the ordering looks correct, inspect the dangling term. Dropping it scales every rank by one constant, so no pair swaps.
Print the sum after every iteration. It should remain 1, and the first pass that drifts locates the failure.
PageRank reads in-edges, but citation data ships as out-edges. If you transpose in the wrong direction, the program runs and converges but ranks papers by their outgoing bibliographies.
Use the top-twenty table to diagnose this error. Compare the reported in- and out-degrees with the titles at the top.
An atomicAdd reduction is short, but every residual thread targets one
address. Even the dangling reduction sends a few thousand threads to the
same location.
That cost comes from contention, not bandwidth. The program times both versions. Day 26 has the crossover.
Do not treat a high effective-bandwidth value as proof that the gather uses DRAM well. On a graph small enough to sit in L2, a memory-bound kernel can be limited by cache rather than DRAM, so its effective bandwidth can exceed the DRAM rate.
Use the two copy rows to identify which memory level limits the run. Scale the graph until the small and large copy rows agree before assessing DRAM gather traffic.
The two drivers have the same line count. The custom path adds 105 lines of device code, while the library path moves 1.80 times as many bytes.
Use both costs when choosing an implementation. State whether maintenance or measured run time controls your choice.
A block reduction that stores one value per warp must also handle a block with one warp. Test the shuffle reduction at 32 threads per block before using 256.
Sources
- Brin and Page, "The Anatomy of a Large-Scale Hypertextual Web Search Engine", for the PageRank definition and damping factor (checked 2026-08-30)
- CCCL 3.0 migration guide, for the removal of
cub::DeviceSpmv(checked 2026-08-30) - NVIDIA
cuda-samples,cpp/2_Concepts_and_Techniques/shfl_scan, for warp-level scan and reduction examples (checked 2026-08-30) - Programming Massively Parallel Processors, fourth edition, chapters 14 and 15 on sparse matrices and graph traversal
Next
Day 41 adds NVTX ranges to this program and uses a timeline to find the slowest of its five launches. That trace tests the stage-level explanation you wrote here.