code/day40-pagerank/pagerank.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 40, capstone 2: PageRank on a real citation graph.
//
// Ranks 28,042 papers from OpenAlex by the PageRank of their citation graph,
// twice: once with kernels this course has already taught, and once with
// Thrust and CUB from day 39. Both are timed, both are checked against the
// same double-precision CPU reference, and the program prints the byte count
// each one moves, so the two are compared on traffic and not only on
// milliseconds.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o pagerank pagerank.cu
// Run: ./pagerank from code/day40-pagerank/, so data/ resolves
//
// 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 <cfloat>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cuda_runtime.h>
#include <cub/cub.cuh>
#include <thrust/execution_policy.h>
#include <thrust/functional.h>
#include <thrust/gather.h>
#include <thrust/transform.h>
// The one error macro, copied verbatim. `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)
// Damping is 0.85 because Brin and Page set it there and every published
// PageRank number assumes it. Fixing it is what makes one learner's ranking
// comparable to another's.
constexpr float kDamping = 0.85f;
// The correctness runs use a fixed iteration count on both the GPU and the
// CPU, so the only difference between the two answers is float against
// double and never "one of them stopped earlier". The convergence run is
// separate and reports what the tolerance actually cost in iterations.
constexpr int kCheckIterations = 100;
constexpr int kPerfIterations = 50;
constexpr int kMaxIterations = 500;
constexpr float kTolerance = 1.0e-6f;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarpSize = 32; // 32 on every GPU this course targets
constexpr int kReduceBlocks = 256; // partials written by a reduction pass
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
// Ranks sit near 1/n, so the comparison is relative with an absolute floor
// under it. The iteration is a contraction with rate kDamping, so float error
// settles at roughly sqrt(maxInDegree) * eps / (1 - kDamping) instead of
// accumulating, which is why 1e-3 is loose enough never to fire on correct
// code and tight enough to catch a reduction that drops terms.
constexpr float kRankRelTol = 1.0e-3f;
constexpr float kRankAbsTol = 1.0e-9f;
// Case `real` grades rank order, not float equality. Two correct
// implementations disagree in the last digits and agree completely on who is
// at the top, and it is the top that the page prints.
constexpr size_t kSpearmanTop = 100;
constexpr double kSpearmanMin = 0.999;
constexpr size_t kTopK = 20;
// The streaming baseline: 16 Mi floats is 64 MiB per buffer, past any L2 on
// the cards this course targets, so it measures DRAM and not cache.
constexpr size_t kCopyElems = 16ull * 1024ull * 1024ull;
static_assert(kThreadsPerBlock % kWarpSize == 0,
"block size must be a whole number of warps");
static_assert(kThreadsPerBlock <= 1024,
"1024 threads per block is the maximum on every compute "
"capability this course targets");
// Times a launch with CUDA events and returns the mean milliseconds per run.
//
// This is the one template and the one lambda the early modules allow. Copy
// it verbatim; the alternative is a dozen 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. The first launch of each kernel pays a module load cost, so
// every kernel that gets timed also gets warmed.
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;
}
// ----------------------------------------------------------------- kernels
// Sums one value per thread across the whole block and returns the total in
// thread 0. Warp shuffles first (day 23), then one word per warp through
// shared memory, then one more shuffle pass over those words (day 24).
//
// Correct at every block size that is a whole number of warps up to 1024,
// including a single-warp block: `warps` is then 1, lane 0 reads the only
// slot and the other 31 lanes contribute zero.
//
// The __shared__ array is one static allocation per kernel that inlines this
// function, so a kernel may call it once. Two calls with no barrier between
// them would alias.
__device__ float blockReduceSum(float value) {
__shared__ float warpSums[kWarpSize];
const unsigned int lane = threadIdx.x % kWarpSize;
const unsigned int warp = threadIdx.x / kWarpSize;
for (int offset = kWarpSize / 2; offset > 0; offset /= 2) {
value += __shfl_down_sync(0xffffffffu, value, offset);
}
if (lane == 0) {
warpSums[warp] = value;
}
__syncthreads();
const unsigned int warps = (blockDim.x + kWarpSize - 1) / kWarpSize;
value = (threadIdx.x < warps) ? warpSums[threadIdx.x] : 0.0f;
if (warp == 0) {
for (int offset = kWarpSize / 2; offset > 0; offset /= 2) {
value += __shfl_down_sync(0xffffffffu, value, offset);
}
}
return value;
}
// out[i] = in[i]. One thread owns one element and then strides forward.
//
// Memory: consecutive threads take consecutive elements, so one warp's 32
// addresses cover 128 contiguous bytes on the read and on the write. This is
// the ceiling every other row of the table is measured against.
__global__ void copyFloats(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
out[i] = in[i];
}
}
// out[i] = value. Resets the rank vector before a run.
//
// Memory: one contiguous store per lane, no reads.
__global__ void fillFloat(float* __restrict__ out, float value, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
out[i] = value;
}
}
// contrib[u] = x[u] / outdeg[u], with invOut holding the reciprocal so the
// kernel multiplies instead of dividing, and holding zero for a node with no
// out-edges so its mass never leaves through this path.
//
// Memory: three contiguous streams, all coalesced. Doing this once per
// iteration is what lets the gather inside the SpMV touch one array instead
// of two, which halves the irregular traffic, the part that actually costs.
__global__ void scaleByOutDegree(const float* __restrict__ x,
const float* __restrict__ invOut,
float* __restrict__ contrib, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
contrib[i] = x[i] * invOut[i];
}
}
// y[row] = sum of contrib[u] over the citers u of row. One thread owns one
// row of the CSR.
//
// Memory: rowPtr and y are coalesced and everything else is not. The 32 lanes
// of a warp sit on 32 different rows, so at any instant they read 32
// unrelated positions of colIdx, and each of those indexes an unrelated
// position of contrib. Rows also differ in length, so the lanes finish at
// different times and the warp runs at the speed of its longest row.
//
// Launch assumption: gridDim.x * blockDim.x >= n.
__global__ void spmvCsrScalar(const uint32_t* __restrict__ rowPtr,
const uint32_t* __restrict__ colIdx,
const float* __restrict__ contrib,
float* __restrict__ y, size_t n) {
const size_t row =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (row < n) {
const uint32_t begin = rowPtr[row];
const uint32_t end = rowPtr[row + 1];
float sum = 0.0f;
for (uint32_t k = begin; k < end; ++k) {
sum += contrib[colIdx[k]];
}
y[row] = sum;
}
}
// The same product with one warp per row and a shuffle reduction.
//
// Memory: the 32 lanes walk one row together, so their colIdx reads are 32
// consecutive uint32 values, which is one 128-byte span. The gather into
// contrib is still scattered, and that is the part no rearrangement fixes.
//
// Launch assumption: blockDim.x is a whole number of warps. The loop strides
// over warps rather than over rows, so every lane of a warp is on the same
// row at the same time, which is what makes the full-mask shuffle legal.
__global__ void spmvCsrWarp(const uint32_t* __restrict__ rowPtr,
const uint32_t* __restrict__ colIdx,
const float* __restrict__ contrib,
float* __restrict__ y, size_t n) {
// snippet: warp-row-loop
const unsigned int lane = threadIdx.x % kWarpSize;
const size_t warpsPerBlock = blockDim.x / kWarpSize;
const size_t warpsTotal = gridDim.x * warpsPerBlock;
const size_t first = blockIdx.x * warpsPerBlock + threadIdx.x / kWarpSize;
for (size_t row = first; row < n; row += warpsTotal) {
const uint32_t begin = rowPtr[row];
const uint32_t end = rowPtr[row + 1];
float sum = 0.0f;
for (uint32_t k = begin + lane; k < end; k += kWarpSize) {
sum += contrib[colIdx[k]];
}
for (int offset = kWarpSize / 2; offset > 0; offset /= 2) {
sum += __shfl_down_sync(0xffffffffu, sum, offset);
}
if (lane == 0) {
y[row] = sum;
}
}
// end snippet
}
// Sums x[i] over the nodes with no out-edge into one float per block.
//
// Memory: two contiguous streams. The mask is one byte per node, so a warp
// reads 32 bytes of it, which is exactly one sector.
//
// Writes: gridDim.x floats, one per block, which is the budget this stage is
// graded on. An atomicAdd per thread would write n times to one address.
__global__ void sumDanglingPartial(const float* __restrict__ x,
const unsigned char* __restrict__ dangling,
float* __restrict__ partials, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float acc = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
acc += dangling[i] ? x[i] : 0.0f;
}
const float total = blockReduceSum(acc);
if (threadIdx.x == 0) {
partials[blockIdx.x] = total;
}
}
// The L1 distance between two rank vectors, as one float per block. Same
// shape as the kernel above, so the two are comparable.
//
// Memory: two coalesced streams. Writes: gridDim.x floats.
__global__ void sumAbsDiffPartial(const float* __restrict__ x,
const float* __restrict__ xPrev,
float* __restrict__ partials, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float acc = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
acc += fabsf(x[i] - xPrev[i]);
}
const float total = blockReduceSum(acc);
if (threadIdx.x == 0) {
partials[blockIdx.x] = total;
}
}
// The same L1 distance written the way it gets written first: one atomicAdd
// per thread into one global float.
//
// Memory: the reads are identical to the kernel above. The writes are not.
// Every thread in the grid targets one address, so the hardware serialises
// them and this kernel measures contention rather than bandwidth. It exists
// to be measured against the tree reduction, which is day 26.
__global__ void sumAbsDiffAtomic(const float* __restrict__ x,
const float* __restrict__ xPrev,
float* __restrict__ out, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
atomicAdd(out, fabsf(x[i] - xPrev[i]));
}
}
// Second pass of a two-pass reduction: one block sums the partials.
//
// Memory: one short contiguous read. Writes: one float.
//
// Launch assumption: exactly one block, so blockReduceSum's result in thread
// 0 is the whole total.
__global__ void sumPartials(const float* __restrict__ partials,
float* __restrict__ out, int count) {
float acc = 0.0f;
for (int i = static_cast<int>(threadIdx.x); i < count;
i += static_cast<int>(blockDim.x)) {
acc += partials[i];
}
const float total = blockReduceSum(acc);
if (threadIdx.x == 0) {
out[0] = total;
}
}
// xNext[v] = (1 - d)/n + d * (y[v] + danglingMass/n).
//
// Memory: two contiguous streams plus one broadcast read. Every thread reads
// danglingMass[0], which the hardware serves once for the whole warp.
//
// danglingMass is a device pointer on purpose. Bringing that one float to the
// host to fold it into a kernel argument would put a synchronize in the
// middle of every iteration, and this loop runs fifty of them.
__global__ void combineRanks(const float* __restrict__ y,
const float* __restrict__ danglingMass,
float* __restrict__ xNext, size_t n,
float damping) {
const float invN = 1.0f / static_cast<float>(n);
const float share =
(1.0f - damping) * invN + damping * danglingMass[0] * invN;
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
xNext[i] = share + damping * y[i];
}
}
// -------------------------------------------------------- Thrust functors
// Three function objects, because Thrust takes the operation as a value and
// none of these three is in <thrust/functional.h>. They are the whole of the
// C++ that the library version adds over the hand-written one.
struct MaskedValue {
__host__ __device__ float operator()(float value,
unsigned char keep) const {
return keep ? value : 0.0f;
}
};
struct AbsDiff {
__host__ __device__ float operator()(float a, float b) const {
return fabsf(a - b);
}
};
// Reads the dangling mass through a device pointer for the same reason
// combineRanks does: to keep the host out of the iteration.
struct Combine {
const float* danglingMass;
float damping;
float invN;
__host__ __device__ float operator()(float yv) const {
return (1.0f - damping) * invN +
damping * (yv + danglingMass[0] * invN);
}
};
// ------------------------------------------------------------- host graph
// A graph in the orientation PageRank walks: for node v, the list of nodes
// that cite it. The data file ships out-edges, so buildGraph transposes.
struct Graph {
size_t n;
size_t nnz;
std::vector<uint32_t> rowPtr; // n + 1
std::vector<uint32_t> colIdx; // nnz, the citing nodes
std::vector<uint32_t> outDeg; // n
};
// Transposes an out-edge list into CSR over in-edges with a counting sort.
//
// Counting the in-degrees is day 29's histogram, the running sum over those
// counts is day 32's exclusive scan, and the cursor pass is the scatter. All
// three run on the host, once per graph, because the capstone is about
// PageRank and not about writing a transpose.
static Graph buildGraph(const std::vector<std::vector<uint32_t>>& out) {
Graph g;
g.n = out.size();
g.outDeg.assign(g.n, 0);
std::vector<uint32_t> inDeg(g.n, 0);
size_t nnz = 0;
for (size_t u = 0; u < g.n; ++u) {
g.outDeg[u] = static_cast<uint32_t>(out[u].size());
nnz += out[u].size();
for (size_t k = 0; k < out[u].size(); ++k) {
inDeg[out[u][k]] += 1;
}
}
g.nnz = nnz;
g.rowPtr.assign(g.n + 1, 0);
for (size_t v = 0; v < g.n; ++v) {
g.rowPtr[v + 1] = g.rowPtr[v] + inDeg[v];
}
g.colIdx.assign(nnz, 0);
std::vector<uint32_t> cursor(g.rowPtr.begin(), g.rowPtr.end() - 1);
for (size_t u = 0; u < g.n; ++u) {
for (size_t k = 0; k < out[u].size(); ++k) {
g.colIdx[cursor[out[u][k]]++] = static_cast<uint32_t>(u);
}
}
return g;
}
// Six nodes, nine edges. Node 5 has no out-edge, so it is dangling and its
// mass has to be collected and redistributed every iteration. Nodes 0 and 1
// cite each other, which is the two-cycle that makes the fixed point worth
// working out by hand.
static std::vector<std::vector<uint32_t>> makeTiny() {
return {{1u, 2u}, {0u, 2u}, {3u}, {0u, 4u}, {2u, 5u}, {}};
}
// One node, no edges. Its rank is 1.0 and the loop has to terminate.
static std::vector<std::vector<uint32_t>> makeSingle() {
return {{}};
}
// A ring with a chord, so every node past the first `danglers` has
// out-degree two. The multiplier 7 is odd, so the chord never lands on the
// node it left and the graph carries no self-loop.
static std::vector<std::vector<uint32_t>> makeRing(size_t n, size_t danglers) {
std::vector<std::vector<uint32_t>> out(n);
for (size_t u = danglers; u < n; ++u) {
out[u].push_back(static_cast<uint32_t>((u + 1) % n));
out[u].push_back(static_cast<uint32_t>((u * 7 + 3) % n));
}
return out;
}
// One node cited by every other node, and citing nobody. The centre's rank
// solves a two-line system, so an off-by-one in the row pointer moves it
// somewhere a check can see.
static std::vector<std::vector<uint32_t>> makeStar(size_t n) {
std::vector<std::vector<uint32_t>> out(n);
for (size_t u = 1; u < n; ++u) {
out[u].push_back(0u);
}
return out;
}
// A path. The slowest convergence in the ladder, so a wrong tolerance shows
// up as a wrong iteration count rather than as a wrong answer.
static std::vector<std::vector<uint32_t>> makeChain(size_t n) {
std::vector<std::vector<uint32_t>> out(n);
for (size_t u = 0; u + 1 < n; ++u) {
out[u].push_back(static_cast<uint32_t>(u + 1));
}
return out;
}
// Opens the first of two paths that exists, so the program runs from the day
// directory or from the repository root. The buffer is static and shared, so
// the caller reads the name before calling again.
static std::FILE* openData(const char* name, const char** used) {
static char buffer[256];
std::snprintf(buffer, sizeof(buffer), "data/%s", name);
std::FILE* f = std::fopen(buffer, "rb");
if (f == nullptr) {
std::snprintf(buffer, sizeof(buffer), "code/day40-pagerank/data/%s",
name);
f = std::fopen(buffer, "rb");
}
*used = buffer;
return f;
}
// Reads data/graph.txt: comment lines, then a node count and an edge count,
// then one line per node holding its out-degree and that many node ids.
//
// std::fscanf over 342,000 integers is not the fastest way to read 1.8 MB
// and it is the shortest correct one. The file is read once and the run
// takes seconds, so a hand-rolled parser would buy nothing.
static bool readGraph(std::vector<std::vector<uint32_t>>* out, size_t* edges) {
const char* path = nullptr;
std::FILE* f = openData("graph.txt", &path);
if (f == nullptr) {
std::fprintf(stderr, "cannot open %s; run from code/day40-pagerank/\n",
path);
return false;
}
int c = 0;
while ((c = std::fgetc(f)) != EOF) {
if (c == '#') {
while ((c = std::fgetc(f)) != EOF && c != '\n') {
}
} else if (c != ' ' && c != '\n' && c != '\r' && c != '\t') {
std::ungetc(c, f);
break;
}
}
unsigned long long n = 0;
unsigned long long m = 0;
bool ok = std::fscanf(f, "%llu %llu", &n, &m) == 2 && n > 0;
size_t total = 0;
if (ok) {
out->clear();
out->resize(static_cast<size_t>(n));
for (size_t u = 0; ok && u < static_cast<size_t>(n); ++u) {
unsigned long long degree = 0;
ok = std::fscanf(f, "%llu", °ree) == 1;
for (size_t k = 0; ok && k < static_cast<size_t>(degree); ++k) {
unsigned long long v = 0;
ok = std::fscanf(f, "%llu", &v) == 1 && v < n;
if (ok) {
(*out)[u].push_back(static_cast<uint32_t>(v));
total += 1;
}
}
}
ok = ok && total == static_cast<size_t>(m);
}
std::fclose(f);
*edges = total;
if (!ok) {
std::fprintf(stderr, "%s is malformed\n", path);
}
return ok;
}
// Titles, as one blob of bytes plus an offset per node, so printing twenty
// of them needs no string container.
struct Titles {
std::vector<char> blob;
std::vector<size_t> offset;
};
static bool readTitles(size_t n, Titles* t) {
const char* path = nullptr;
std::FILE* f = openData("titles.tsv", &path);
if (f == nullptr) {
std::fprintf(stderr, "cannot open %s\n", path);
return false;
}
std::fseek(f, 0, SEEK_END);
const size_t size = static_cast<size_t>(std::ftell(f));
std::fseek(f, 0, SEEK_SET);
t->blob.assign(size + 1, '\0');
const size_t got = std::fread(t->blob.data(), 1, size, f);
std::fclose(f);
if (got != size || got == 0) {
std::fprintf(stderr, "short read on %s\n", path);
return false;
}
t->offset.clear();
t->offset.push_back(0);
for (size_t i = 0; i < got; ++i) {
if (t->blob[i] == '\n') {
t->blob[i] = '\0';
if (i + 1 < got) {
t->offset.push_back(i + 1);
}
}
}
if (t->offset.size() != n) {
std::fprintf(stderr, "%s has %zu lines, the graph has %zu nodes\n",
path, t->offset.size(), n);
return false;
}
return true;
}
// ---------------------------------------------------------- CPU reference
// PageRank in double, written for obvious correctness and not for speed: one
// loop over the rows, one over each row's citers, no blocking and no
// intrinsics. It is the oracle for every case in the ladder.
//
// Returns the iteration count. Stops on the L1 distance between successive
// vectors when tolerance is above zero, and otherwise runs the full count.
static int pagerankCpu(const Graph& g, double damping, double tolerance,
int maxIterations, std::vector<double>* rank) {
const double invN = 1.0 / static_cast<double>(g.n);
std::vector<double> x(g.n, invN);
std::vector<double> next(g.n, 0.0);
std::vector<double> contrib(g.n, 0.0);
int iteration = 0;
for (; iteration < maxIterations; ++iteration) {
double danglingMass = 0.0;
for (size_t u = 0; u < g.n; ++u) {
if (g.outDeg[u] == 0) {
danglingMass += x[u];
contrib[u] = 0.0;
} else {
contrib[u] = x[u] / static_cast<double>(g.outDeg[u]);
}
}
const double share =
(1.0 - damping) * invN + damping * danglingMass * invN;
double residual = 0.0;
for (size_t v = 0; v < g.n; ++v) {
double sum = 0.0;
for (uint32_t k = g.rowPtr[v]; k < g.rowPtr[v + 1]; ++k) {
sum += contrib[g.colIdx[k]];
}
next[v] = share + damping * sum;
residual += std::fabs(next[v] - x[v]);
}
x.swap(next);
if (tolerance > 0.0 && residual < tolerance) {
iteration += 1;
break;
}
}
*rank = x;
return iteration;
}
// ------------------------------------------------------------ device side
struct DeviceGraph {
uint32_t* rowPtr;
uint32_t* colIdx;
float* invOut;
unsigned char* dangling;
size_t n;
size_t nnz;
};
struct Work {
float* x;
float* xNext;
float* y;
float* contrib;
float* gathered; // nnz floats, written and read back by the library path
float* diff; // n floats, the same
float* partials; // kReduceBlocks floats
float* scalars; // [0] dangling mass, [1] residual
void* cubTemp;
size_t cubSegBytes;
size_t cubRedBytes;
};
struct PrParams {
float damping;
float tolerance; // zero disables the convergence test
int maxIterations;
bool warpPerRow;
};
static int gridFor(size_t n) {
const size_t blocks =
(n + kThreadsPerBlock - 1) / static_cast<size_t>(kThreadsPerBlock);
return static_cast<int>(blocks < 1 ? 1 : blocks);
}
// One warp per row, so the grid is sized in warps and not in threads.
static int gridForWarps(size_t n) {
const size_t warpsPerBlock = kThreadsPerBlock / kWarpSize;
const size_t blocks = (n + warpsPerBlock - 1) / warpsPerBlock;
return static_cast<int>(blocks < 1 ? 1 : blocks);
}
static void uploadGraph(const Graph& g, DeviceGraph* d) {
std::vector<float> h_invOut(g.n, 0.0f);
std::vector<unsigned char> h_dangling(g.n, 0);
for (size_t u = 0; u < g.n; ++u) {
if (g.outDeg[u] == 0) {
h_dangling[u] = 1;
} else {
h_invOut[u] = 1.0f / static_cast<float>(g.outDeg[u]);
}
}
d->n = g.n;
d->nnz = g.nnz;
// The +1 on colIdx keeps the allocation non-empty for the `single` case,
// where the graph has one node and no edges at all.
CUDA_CHECK(cudaMalloc(&d->rowPtr, (g.n + 1) * sizeof(uint32_t)));
CUDA_CHECK(cudaMalloc(&d->colIdx, (g.nnz + 1) * sizeof(uint32_t)));
CUDA_CHECK(cudaMalloc(&d->invOut, g.n * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d->dangling, g.n * sizeof(unsigned char)));
CUDA_CHECK(cudaMemcpy(d->rowPtr, g.rowPtr.data(),
(g.n + 1) * sizeof(uint32_t),
cudaMemcpyHostToDevice));
if (g.nnz > 0) {
CUDA_CHECK(cudaMemcpy(d->colIdx, g.colIdx.data(),
g.nnz * sizeof(uint32_t),
cudaMemcpyHostToDevice));
}
CUDA_CHECK(cudaMemcpy(d->invOut, h_invOut.data(), g.n * sizeof(float),
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d->dangling, h_dangling.data(),
g.n * sizeof(unsigned char), cudaMemcpyHostToDevice));
}
static void freeGraph(DeviceGraph* d) {
CUDA_CHECK(cudaFree(d->rowPtr));
CUDA_CHECK(cudaFree(d->colIdx));
CUDA_CHECK(cudaFree(d->invOut));
CUDA_CHECK(cudaFree(d->dangling));
}
// CUB asks for scratch space in two steps: call it with a null pointer to
// learn the size, then call it again with the allocation. The two reductions
// share one buffer, because they never run at the same time.
static void allocWork(const DeviceGraph& g, Work* w) {
const size_t nBytes = g.n * sizeof(float);
CUDA_CHECK(cudaMalloc(&w->x, nBytes));
CUDA_CHECK(cudaMalloc(&w->xNext, nBytes));
CUDA_CHECK(cudaMalloc(&w->y, nBytes));
CUDA_CHECK(cudaMalloc(&w->contrib, nBytes));
CUDA_CHECK(cudaMalloc(&w->diff, nBytes));
CUDA_CHECK(cudaMalloc(&w->gathered, (g.nnz + 1) * sizeof(float)));
CUDA_CHECK(cudaMalloc(&w->partials, kReduceBlocks * sizeof(float)));
CUDA_CHECK(cudaMalloc(&w->scalars, 2 * sizeof(float)));
w->cubTemp = nullptr;
w->cubSegBytes = 0;
w->cubRedBytes = 0;
CUDA_CHECK(cub::DeviceSegmentedReduce::Sum(
nullptr, w->cubSegBytes, w->gathered, w->y, static_cast<int>(g.n),
g.rowPtr, g.rowPtr + 1));
CUDA_CHECK(cub::DeviceReduce::Sum(nullptr, w->cubRedBytes, w->diff,
w->scalars, g.n));
const size_t bytes =
(w->cubSegBytes > w->cubRedBytes) ? w->cubSegBytes : w->cubRedBytes;
CUDA_CHECK(cudaMalloc(&w->cubTemp, bytes + 1));
}
static void freeWork(Work* w) {
CUDA_CHECK(cudaFree(w->x));
CUDA_CHECK(cudaFree(w->xNext));
CUDA_CHECK(cudaFree(w->y));
CUDA_CHECK(cudaFree(w->contrib));
CUDA_CHECK(cudaFree(w->diff));
CUDA_CHECK(cudaFree(w->gathered));
CUDA_CHECK(cudaFree(w->partials));
CUDA_CHECK(cudaFree(w->scalars));
CUDA_CHECK(cudaFree(w->cubTemp));
}
// One PageRank iteration with the course's own kernels. Four or five
// launches, no synchronize, and nothing copied back to the host.
// snippet: hand-iteration
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);
}
// end snippet
// The same iteration with Thrust and CUB, and no kernel of ours.
//
// The shape is forced by what the libraries offer. There is no fused
// gather-and-segmented-reduce, so the gathered edge values are written to
// memory and read straight back, and the masked rank vector is materialised
// before it is summed. Both arrays are traffic the hand-written version
// never moves, which is why the program prints both byte counts.
// snippet: library-iteration
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});
}
// end snippet
// The L1 residual, kept out of the iteration so the timed run can leave it
// out and both implementations time the same work.
static void residualHand(const Work& w, const float* x, const float* xPrev,
size_t n) {
sumAbsDiffPartial<<<kReduceBlocks, kThreadsPerBlock>>>(x, xPrev, w.partials,
n);
sumPartials<<<1, kThreadsPerBlock>>>(w.partials, w.scalars + 1,
kReduceBlocks);
}
static void residualLibrary(Work& w, const float* x, const float* xPrev,
size_t n) {
thrust::transform(thrust::device, x, x + n, xPrev, w.diff, AbsDiff());
CUDA_CHECK(cub::DeviceReduce::Sum(w.cubTemp, w.cubRedBytes, w.diff,
w.scalars + 1, n));
}
// Runs PageRank to a fixed count or to a tolerance and returns the pointer
// holding the answer, which is w.x or w.xNext depending on how many times
// the two were swapped. No synchronize: the caller either times this or
// copies the answer back, and both of those are sync points already.
static float* runPagerank(const DeviceGraph& g, Work& w, const PrParams& p,
bool useLibrary, int* iterations) {
float* x = w.x;
float* xNext = w.xNext;
fillFloat<<<gridFor(g.n), kThreadsPerBlock>>>(
x, 1.0f / static_cast<float>(g.n), g.n);
CUDA_CHECK(cudaGetLastError());
int iteration = 0;
for (; iteration < p.maxIterations; ++iteration) {
if (useLibrary) {
iterateLibrary(g, w, x, xNext, p);
} else {
iterateHand(g, w, x, xNext, p);
}
if (p.tolerance > 0.0f) {
if (useLibrary) {
residualLibrary(w, xNext, x, g.n);
} else {
residualHand(w, xNext, x, g.n);
}
}
float* swap = x;
x = xNext;
xNext = swap;
// Reading the residual is the one place the host waits for the
// device inside the loop. The timed run sets tolerance to zero and
// pays none of it, which is why the two numbers are not comparable
// and the program reports them separately.
if (p.tolerance > 0.0f) {
float residual = 0.0f;
CUDA_CHECK(cudaMemcpy(&residual, w.scalars + 1, sizeof(float),
cudaMemcpyDeviceToHost));
if (residual < p.tolerance) {
iteration += 1;
break;
}
}
}
CUDA_CHECK(cudaGetLastError());
*iterations = iteration;
return x;
}
// ---------------------------------------------------------------- checking
// The largest per-element error as a multiple of what the tolerance allows,
// so a pass is a number at or below 1.0. Returns -1 when anything is not
// finite, which is what a wrong reduction usually produces first.
static double toleranceRatio(const std::vector<float>& got,
const std::vector<double>& want) {
double worst = 0.0;
for (size_t i = 0; i < want.size(); ++i) {
if (!std::isfinite(got[i])) {
return -1.0;
}
const double allowed =
static_cast<double>(kRankAbsTol) +
static_cast<double>(kRankRelTol) * std::fabs(want[i]);
const double error = std::fabs(static_cast<double>(got[i]) - want[i]);
const double ratio = error / allowed;
if (ratio > worst) {
worst = ratio;
}
}
return worst;
}
static double sumOf(const std::vector<float>& v) {
double total = 0.0;
for (size_t i = 0; i < v.size(); ++i) {
total += static_cast<double>(v[i]);
}
return total;
}
// Every node's position in a descending ordering by score. Ties break on the
// node id so the ordering is total and the comparison is reproducible.
template <typename T>
static std::vector<size_t> positionsByScore(const std::vector<T>& score) {
std::vector<size_t> order(score.size());
for (size_t i = 0; i < order.size(); ++i) {
order[i] = i;
}
std::sort(order.begin(), order.end(), [&score](size_t a, size_t b) {
return score[a] > score[b] || (score[a] == score[b] && a < b);
});
std::vector<size_t> position(score.size(), 0);
for (size_t p = 0; p < order.size(); ++p) {
position[order[p]] = p;
}
return position;
}
// Spearman's rho over the reference's top k papers, using each paper's
// position in the full ordering on both sides. A paper the GPU pushes out of
// its own top k is penalised by how far it actually fell, which is stricter
// than re-ranking inside the subset and is the number worth reporting.
static double spearmanTopK(const std::vector<double>& want,
const std::vector<float>& got, size_t k) {
const std::vector<size_t> posWant = positionsByScore(want);
const std::vector<size_t> posGot = positionsByScore(got);
double sumSquared = 0.0;
size_t counted = 0;
for (size_t node = 0; node < want.size(); ++node) {
if (posWant[node] < k) {
const double d = static_cast<double>(posWant[node]) -
static_cast<double>(posGot[node]);
sumSquared += d * d;
counted += 1;
}
}
const double m = static_cast<double>(counted);
return 1.0 - 6.0 * sumSquared / (m * (m * m - 1.0));
}
// Runs one correctness case through all three implementations and prints one
// line each. Returns false if any of them is wrong.
static bool runCase(const char* name,
const std::vector<std::vector<uint32_t>>& outEdges) {
const Graph g = buildGraph(outEdges);
std::vector<double> want;
pagerankCpu(g, static_cast<double>(kDamping), 0.0, kCheckIterations, &want);
const double sumAllowed = 4.0 * static_cast<double>(FLT_EPSILON) *
std::sqrt(static_cast<double>(g.n)) +
1.0e-6;
DeviceGraph d;
uploadGraph(g, &d);
Work w;
allocWork(d, &w);
const char* labels[3] = {"thread per row", "warp per row", "Thrust + CUB"};
bool ok = true;
for (int variant = 0; variant < 3; ++variant) {
PrParams p;
p.damping = kDamping;
p.tolerance = 0.0f;
p.maxIterations = kCheckIterations;
p.warpPerRow = (variant == 1);
int iterations = 0;
float* answer = runPagerank(d, w, p, variant == 2, &iterations);
CUDA_CHECK(cudaDeviceSynchronize());
std::vector<float> got(g.n, 0.0f);
CUDA_CHECK(cudaMemcpy(got.data(), answer, g.n * sizeof(float),
cudaMemcpyDeviceToHost));
const double ratio = toleranceRatio(got, want);
const double total = sumOf(got);
const bool pass = ratio >= 0.0 && ratio <= 1.0 &&
std::fabs(total - 1.0) <= sumAllowed;
std::printf(
" %-11s %-15s n=%-6zu nnz=%-7zu err/tol=%6.3f "
"sum-1=%+.2e %s\n",
name, labels[variant], g.n, g.nnz, ratio, total - 1.0,
pass ? "pass" : "FAIL");
ok = ok && pass;
}
freeWork(&w);
freeGraph(&d);
return ok;
}
// ------------------------------------------------------------- the timing
// Bytes one iteration must touch, counted per kernel so the arithmetic is
// visible rather than asserted. Reads and writes both count: the number this
// feeds is effective bandwidth, not DRAM traffic, and on a graph small
// enough to sit in L2 those two are very different things.
static double handIterationBytes(size_t n, size_t nnz) {
const double fn = static_cast<double>(n);
const double fnnz = static_cast<double>(nnz);
return 12.0 * fn // scale: x in, invOut in, contrib out
+ 4.0 * (fn + 1.0) // spmv: rowPtr
+ 8.0 * fnnz // spmv: colIdx, and the gather it drives
+ 4.0 * fn // spmv: y out
+ 5.0 * fn // dangling: x and the one-byte mask
+ 8.0 * fn; // combine: y in, xNext out
}
static double libraryIterationBytes(size_t n, size_t nnz) {
const double fn = static_cast<double>(n);
const double fnnz = static_cast<double>(nnz);
return 12.0 * fn // transform: x in, invOut in, contrib out
+ 8.0 * fnnz // gather: colIdx, and the gather it drives
+ 4.0 * fnnz // gather: the edge values written out
+ 4.0 * fnnz // segmented reduce: the same values back in
+ 4.0 * (fn + 1.0) // segmented reduce: rowPtr
+ 4.0 * fn // segmented reduce: y out
+ 9.0 * fn // transform: x, the mask, the masked copy out
+ 4.0 * fn // reduce: the masked copy back in
+ 8.0 * fn; // transform: y in, xNext out
}
static double gbPerSecond(double bytes, double ms) {
return bytes / (ms * 1.0e-3) / 1.0e9;
}
static void printRow(const char* name, double ms, double bytes,
double dramGbs) {
const double gbs = gbPerSecond(bytes, ms);
std::printf("%-33s %10.4f %9.1f %9.2f\n", name, ms, gbs, gbs / dramGbs);
}
// Runs the shipped graph: correctness, convergence, the timed comparison and
// the top twenty. Returns false if anything is wrong.
static bool runReal(const Graph& g, const Titles& titles, int l2Bytes) {
std::vector<double> want;
const int cpuIterations = pagerankCpu(
g, static_cast<double>(kDamping),
static_cast<double>(kTolerance) * 1.0e-4, kMaxIterations, &want);
DeviceGraph d;
uploadGraph(g, &d);
Work w;
allocWork(d, &w);
PrParams converge;
converge.damping = kDamping;
converge.tolerance = kTolerance;
converge.maxIterations = kMaxIterations;
converge.warpPerRow = true;
int gpuIterations = 0;
float* answer = runPagerank(d, w, converge, false, &gpuIterations);
CUDA_CHECK(cudaDeviceSynchronize());
std::vector<float> got(g.n, 0.0f);
CUDA_CHECK(cudaMemcpy(got.data(), answer, g.n * sizeof(float),
cudaMemcpyDeviceToHost));
// The library path is a required deliverable, so it is checked on the
// real graph too. Comparing the run time of an implementation nobody
// graded is worth nothing.
int libIterations = 0;
float* libAnswer = runPagerank(d, w, converge, true, &libIterations);
CUDA_CHECK(cudaDeviceSynchronize());
std::vector<float> gotLibrary(g.n, 0.0f);
CUDA_CHECK(cudaMemcpy(gotLibrary.data(), libAnswer, g.n * sizeof(float),
cudaMemcpyDeviceToHost));
// Everything below sorts by rank, and std::sort with a comparator that
// sees a NaN is undefined behaviour rather than a wrong answer. A broken
// reduction produces NaN before it produces anything else, so the guard
// comes first and reports instead of crashing.
bool finite = true;
for (size_t i = 0; i < g.n; ++i) {
if (!std::isfinite(got[i]) || !std::isfinite(gotLibrary[i])) {
finite = false;
}
}
if (!finite) {
std::fprintf(stderr,
"case real: the rank vector is not finite; "
"check the reductions before reading anything "
"below\n");
freeWork(&w);
freeGraph(&d);
return false;
}
const double total = sumOf(got);
const double totalLibrary = sumOf(gotLibrary);
const double sumAllowed = 4.0 * static_cast<double>(FLT_EPSILON) *
std::sqrt(static_cast<double>(g.n));
const double rho = spearmanTopK(want, got, kSpearmanTop);
const double rhoLibrary = spearmanTopK(want, gotLibrary, kSpearmanTop);
const bool sumOk = std::fabs(total - 1.0) <= sumAllowed &&
std::fabs(totalLibrary - 1.0) <= sumAllowed;
const bool rhoOk = rho >= kSpearmanMin && rhoLibrary >= kSpearmanMin;
std::printf("\ncase real %zu nodes, %zu edges\n", g.n, g.nnz);
std::printf(" CPU reference, double, converged in %d iterations\n",
cpuIterations);
std::printf(
" our kernels, float, L1 residual under %.0e in %d "
"iterations\n",
static_cast<double>(kTolerance), gpuIterations);
std::printf(
" Thrust + CUB, float, the same tolerance in %d "
"iterations\n",
libIterations);
std::printf(
" sum of ranks minus 1: ours %+.3e, library %+.3e "
"(allowed %.3e) %s\n",
total - 1.0, totalLibrary - 1.0, sumAllowed, sumOk ? "pass" : "FAIL");
std::printf(
" Spearman rho over the reference's top %zu: ours %.6f, "
"library %.6f, need %.3f %s\n",
kSpearmanTop, rho, rhoLibrary, kSpearmanMin, rhoOk ? "pass" : "FAIL");
// The top twenty, which is the artifact a human can judge.
std::vector<size_t> order(g.n);
for (size_t i = 0; i < g.n; ++i) {
order[i] = i;
}
std::sort(order.begin(), order.end(), [&got](size_t a, size_t b) {
return got[a] > got[b] || (got[a] == got[b] && a < b);
});
std::printf("\ntop %zu by PageRank (work id, year, title)\n", kTopK);
std::printf(" rank score in-deg out-deg paper\n");
for (size_t r = 0; r < kTopK && r < g.n; ++r) {
const size_t node = order[r];
const uint32_t inDegree = g.rowPtr[node + 1] - g.rowPtr[node];
std::printf(" %4zu %.7f %6u %7u %s\n", r + 1,
static_cast<double>(got[node]), inDegree, g.outDeg[node],
&titles.blob[titles.offset[node]]);
}
// Two copy baselines. The first streams far more than any L2 holds, so
// it is the card's DRAM bandwidth. The second moves exactly what one
// PageRank iteration moves, which on this graph may be small enough that
// the cache answers most of it. The gap between the two rows is why the
// PageRank rows are not compared against the first one alone.
const double handBytes = handIterationBytes(g.n, g.nnz);
const double libBytes = libraryIterationBytes(g.n, g.nnz);
const size_t smallElems = static_cast<size_t>(handBytes / 8.0) + 1;
float* d_copySrc = nullptr;
float* d_copyDst = nullptr;
CUDA_CHECK(cudaMalloc(&d_copySrc, kCopyElems * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_copyDst, kCopyElems * sizeof(float)));
CUDA_CHECK(cudaMemset(d_copySrc, 0, kCopyElems * sizeof(float)));
const float bigCopyMs = timeKernel([&] {
copyFloats<<<gridFor(kCopyElems), kThreadsPerBlock>>>(
d_copySrc, d_copyDst, kCopyElems);
});
const float smallCopyMs = timeKernel([&] {
copyFloats<<<gridFor(smallElems), kThreadsPerBlock>>>(
d_copySrc, d_copyDst, smallElems);
});
const double bigCopyBytes =
2.0 * static_cast<double>(kCopyElems) * sizeof(float);
const double dramGbs =
gbPerSecond(bigCopyBytes, static_cast<double>(bigCopyMs));
// Individual kernels. None of these writes its own input, so a loop of
// ten launches is ten identical measurements.
const int blocks = gridFor(g.n);
const float scaleMs = timeKernel([&] {
scaleByOutDegree<<<blocks, kThreadsPerBlock>>>(w.x, d.invOut, w.contrib,
g.n);
});
const float spmvScalarMs = timeKernel([&] {
spmvCsrScalar<<<blocks, kThreadsPerBlock>>>(d.rowPtr, d.colIdx,
w.contrib, w.y, g.n);
});
const float spmvWarpMs = timeKernel([&] {
spmvCsrWarp<<<gridForWarps(g.n), kThreadsPerBlock>>>(
d.rowPtr, d.colIdx, w.contrib, w.y, g.n);
});
const float spmvLibraryMs = timeKernel([&] {
thrust::gather(thrust::device, d.colIdx, d.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),
d.rowPtr, d.rowPtr + 1));
});
const float danglingMs = timeKernel([&] {
sumDanglingPartial<<<kReduceBlocks, kThreadsPerBlock>>>(
w.x, d.dangling, w.partials, g.n);
sumPartials<<<1, kThreadsPerBlock>>>(w.partials, w.scalars,
kReduceBlocks);
});
const float residualTreeMs =
timeKernel([&] { residualHand(w, w.x, w.xNext, g.n); });
// The atomic version needs its accumulator cleared, and the clear is
// device work on the same stream, so it sits inside the measurement and
// is counted against it. It is one 4-byte write against n atomics.
const float residualAtomicMs = timeKernel([&] {
CUDA_CHECK(cudaMemsetAsync(w.scalars + 1, 0, sizeof(float)));
sumAbsDiffAtomic<<<kReduceBlocks, kThreadsPerBlock>>>(
w.x, w.xNext, w.scalars + 1, g.n);
});
// The two full runs. Fifty iterations from the same starting vector, no
// convergence test on either side, so both time the same work.
PrParams perf;
perf.damping = kDamping;
perf.tolerance = 0.0f;
perf.maxIterations = kPerfIterations;
perf.warpPerRow = true;
int ignored = 0;
const float handRunMs =
timeKernel([&] { runPagerank(d, w, perf, false, &ignored); });
const float libraryRunMs =
timeKernel([&] { runPagerank(d, w, perf, true, &ignored); });
const double fn = static_cast<double>(g.n);
const double fnnz = static_cast<double>(g.nnz);
const double spmvBytes = 4.0 * (fn + 1.0) + 8.0 * fnnz + 4.0 * fn;
const double spmvLibBytes = 4.0 * (fn + 1.0) + 16.0 * fnnz + 4.0 * fn;
std::printf(
"\ntimed on the shipped graph, mean of %d runs after %d "
"warm-ups, copies not included\n",
kTimedRuns, kWarmupRuns);
std::printf(
"one iteration moves %.0f bytes with our kernels and %.0f "
"with the library\n",
handBytes, libBytes);
std::printf("%-33s %10s %9s %9s\n", "stage", "ms", "GB/s", "x DRAM");
printRow("copy, 128 MiB (DRAM)", static_cast<double>(bigCopyMs),
bigCopyBytes, dramGbs);
printRow("copy, one iteration's bytes", static_cast<double>(smallCopyMs),
2.0 * static_cast<double>(smallElems) * sizeof(float), dramGbs);
printRow("scale by out-degree", static_cast<double>(scaleMs), 12.0 * fn,
dramGbs);
printRow("spmv, thread per row", static_cast<double>(spmvScalarMs),
spmvBytes, dramGbs);
printRow("spmv, warp per row", static_cast<double>(spmvWarpMs), spmvBytes,
dramGbs);
printRow("spmv, gather + CUB segmented", static_cast<double>(spmvLibraryMs),
spmvLibBytes, dramGbs);
printRow("dangling mass, tree reduction", static_cast<double>(danglingMs),
5.0 * fn, dramGbs);
printRow("residual, tree reduction", static_cast<double>(residualTreeMs),
8.0 * fn, dramGbs);
printRow("residual, one atomic per thread",
static_cast<double>(residualAtomicMs), 8.0 * fn, dramGbs);
const double handIterMs = static_cast<double>(handRunMs) / kPerfIterations;
const double libIterMs =
static_cast<double>(libraryRunMs) / kPerfIterations;
std::printf(
"\n%d iterations end to end, no convergence test on either "
"side\n",
kPerfIterations);
std::printf("%-33s %10s %9s %9s\n", "implementation", "ms", "ms/iter",
"GB/s");
std::printf("%-33s %10.4f %9.4f %9.1f\n", "our kernels",
static_cast<double>(handRunMs), handIterMs,
gbPerSecond(handBytes, handIterMs));
std::printf("%-33s %10.4f %9.4f %9.1f\n", "Thrust + CUB",
static_cast<double>(libraryRunMs), libIterMs,
gbPerSecond(libBytes, libIterMs));
std::printf(
"our kernels take %.2fx the library's time and move %.2fx "
"its bytes\n",
handIterMs / libIterMs, handBytes / libBytes);
std::printf(
"L2 on this card is %d bytes, and one iteration touches "
"%.0f\n",
l2Bytes, handBytes);
CUDA_CHECK(cudaFree(d_copySrc));
CUDA_CHECK(cudaFree(d_copyDst));
freeWork(&w);
freeGraph(&d);
return sumOk && rhoOk;
}
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)\n", prop.name, prop.major,
prop.minor);
std::printf("%d SMs, %d bytes of L2, %d threads per block\n",
prop.multiProcessorCount, prop.l2CacheSize,
prop.maxThreadsPerBlock);
std::vector<std::vector<uint32_t>> outEdges;
size_t edges = 0;
if (!readGraph(&outEdges, &edges)) {
return EXIT_FAILURE;
}
Titles titles;
if (!readTitles(outEdges.size(), &titles)) {
return EXIT_FAILURE;
}
std::printf("loaded %zu nodes and %zu edges from data/graph.txt\n\n",
outEdges.size(), edges);
// The ladder, smallest first. Every case runs even when an earlier one
// failed, because "passes without dangling nodes and fails with them" is
// the report that names the bug.
std::printf("correctness, %d fixed iterations on both sides\n",
kCheckIterations);
bool ok = true;
ok = runCase("tiny", makeTiny()) && ok;
ok = runCase("single", makeSingle()) && ok;
ok = runCase("nodangling", makeRing(4096, 0)) && ok;
ok = runCase("dangling", makeRing(4096, 400)) && ok;
ok = runCase("star", makeStar(10001)) && ok;
ok = runCase("chain", makeChain(8192)) && ok;
const Graph real = buildGraph(outEdges);
ok = runReal(real, titles, prop.l2CacheSize) && ok;
if (!ok) {
std::fprintf(stderr, "\nat least one case failed\n");
return EXIT_FAILURE;
}
std::printf("\nevery case passed\n");
return EXIT_SUCCESS;
}