code/day41-nsys/pagerank_nvtx.cuThis is the source used by the lesson and its recorded evidence. Compile commands and expected output live in the directory README.
// Day 41: reading an Nsight Systems timeline.
//
// What it does: runs the day 40 PageRank iteration on a synthetic citation
// graph, three times, changing only how often the host reads the convergence
// residual. Every phase is wrapped in an NVTX range, so the same binary
// produces an anonymous timeline under `nsys profile -t cuda` and a named one
// under `nsys profile -t cuda,nvtx`.
//
// What it does not measure: DRAM bandwidth. The whole graph is a few MiB and
// fits in this card's L2, so effective GB/s here would describe the cache and
// not the memory bus. Day 40 prints the two copy rows that show it. This
// program prints milliseconds only.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo -o pagerank_nvtx \
// pagerank_nvtx.cu
// Run: ./pagerank_nvtx
#include <cfloat>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <vector>
#include <cuda_runtime.h>
#include <nvtx3/nvToolsExt.h>
// The one error macro. This file is standalone, so it carries its own
// verbatim copy. `err_` has a trailing underscore so it cannot collide with a
// variable at the call site.
#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)
// The graph is fixed here so the page and the program cannot disagree, and so
// the shipped .nsys-rep describes the same work a reader reproduces.
constexpr size_t kNodes = 28042;
constexpr float kDamping = 0.85f;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kReduceBlocks = 256; // partials the first pass writes
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
// Iteration counts. The timed loops run a fixed 50 so the three policies are
// timed on identical arithmetic; correctness runs a fixed 100 on both sides
// so the GPU and the CPU reference stop at the same place.
constexpr int kWarmupIters = 3;
constexpr int kTimedIters = 50;
constexpr int kCheckIters = 100;
// Loose on purpose. The iteration is a contraction with rate kDamping, so
// float rounding settles rather than accumulating, and 1e-3 still catches a
// reduction that drops terms.
constexpr float kRankRelTol = 1.0e-3f;
constexpr float kRankAbsTol = 1.0e-9f;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert((kThreadsPerBlock & (kThreadsPerBlock - 1)) == 0,
"blockSum halves the stride, so the block size must be a power "
"of two. 32 is the smallest value that satisfies both asserts "
"and it is a size this day's exercise invites.");
static_assert(kReduceBlocks > 0, "the first reduction pass needs a block");
static int gridFor(size_t n) {
return static_cast<int>((n + kThreadsPerBlock - 1) / kThreadsPerBlock);
}
// ---------------------------------------------------------------- kernels
// Sums one value per thread across the block and returns the total in every
// thread. The caller owns `tile`, which holds one float per thread.
//
// Every thread reaches every barrier: the guard is on the accumulate above,
// never on the __syncthreads(). Correct at 32 threads per block, where the
// halving loop runs five times and its first pass already crosses lanes
// inside a single warp.
__device__ float blockSum(float value, float* tile) {
const unsigned int tid = threadIdx.x;
tile[tid] = value;
__syncthreads();
for (unsigned int half = blockDim.x / 2; half > 0; half /= 2) {
if (tid < half) {
tile[tid] += tile[tid + half];
}
__syncthreads();
}
return tile[0];
}
// x[i] = value, for resetting the rank vector between policies.
__global__ void fillFloat(float* __restrict__ out, float value, size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
out[i] = value;
}
}
// contrib[v] = x[v] / outdeg[v]. One thread owns one node.
//
// Memory: consecutive lanes take consecutive nodes, so each of the three
// arrays is one 128-byte span per warp, four sectors, the day 11 best case.
//
// invOut is 0 on a dangling node, so a dangling node's rank contributes
// nothing here and is collected by sumDanglingPartial instead.
__global__ void scaleByOutDegree(const float* __restrict__ x,
const float* __restrict__ invOut,
float* __restrict__ contrib, size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
contrib[i] = x[i] * invOut[i];
}
}
// y[v] = sum of contrib[u] over every u that cites v. One thread owns one CSR
// row, which is day 37's thread-per-row kernel.
//
// Memory: the 32 lanes of a warp hold 32 unrelated rows, so the colIdx reads
// land in up to 32 separate 32-byte sectors, and the gather into contrib
// scatters again. This is the kernel the timeline is expected to show as the
// widest GPU bar.
//
// Launch: one thread per node, any block size.
__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 i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
const uint32_t begin = rowPtr[i];
const uint32_t end = rowPtr[i + 1];
float sum = 0.0f;
for (uint32_t e = begin; e < end; ++e) {
sum += contrib[colIdx[e]];
}
y[i] = sum;
}
}
// partials[blockIdx.x] = sum of x[v] over the dangling v in this block's
// slice. Pass one of two.
//
// Memory: a grid-stride loop, so consecutive lanes hold consecutive nodes on
// every pass and both reads coalesce.
//
// Launch: kReduceBlocks blocks, which with sumPartials makes the pair write
// kReduceBlocks + 1 floats to global memory in total.
__global__ void sumDanglingPartial(const float* __restrict__ x,
const uint32_t* __restrict__ dangling,
float* __restrict__ partials, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float sum = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
sum += (dangling[i] != 0u) ? x[i] : 0.0f;
}
const float total = blockSum(sum, tile);
if (threadIdx.x == 0) {
partials[blockIdx.x] = total;
}
}
// partials[blockIdx.x] = sum of |a[v] - b[v]| over this block's slice, the
// L1 distance between two successive rank vectors. Same shape and same
// memory pattern as sumDanglingPartial.
__global__ void sumAbsDiffPartial(const float* __restrict__ a,
const float* __restrict__ b,
float* __restrict__ partials, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float sum = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
sum += fabsf(a[i] - b[i]);
}
const float total = blockSum(sum, tile);
if (threadIdx.x == 0) {
partials[blockIdx.x] = total;
}
}
// out[0] = sum of the first `count` partials. Pass two of two, one block.
//
// Launch: exactly one block. The strided accumulate above the reduction is
// what lets one block of any size sum any number of partials.
__global__ void sumPartials(const float* __restrict__ partials,
float* __restrict__ out, int count) {
__shared__ float tile[kThreadsPerBlock];
float sum = 0.0f;
for (int k = static_cast<int>(threadIdx.x); k < count;
k += static_cast<int>(blockDim.x)) {
sum += partials[k];
}
const float total = blockSum(sum, tile);
if (threadIdx.x == 0) {
out[0] = total;
}
}
// xNext[v] = (1 - d)/n + d * (y[v] + mass[0]/n). One thread owns one node.
//
// Memory: consecutive lanes take consecutive nodes. All 32 lanes of a warp
// read the single address mass[0], which the hardware serves as one
// broadcast request, so leaving the dangling mass on the device costs nothing
// per thread and saves a device-to-host copy every iteration.
__global__ void combineRanks(const float* __restrict__ y,
const float* __restrict__ mass,
float* __restrict__ xNext, size_t n,
float damping) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
const float invN = 1.0f / static_cast<float>(n);
xNext[i] = (1.0f - damping) * invN + damping * (y[i] + mass[0] * invN);
}
}
// ------------------------------------------------------------ host: graph
struct Graph {
std::vector<uint32_t> rowPtr; // n + 1 entries, in-edge CSR
std::vector<uint32_t> colIdx; // nnz entries, the citing node
std::vector<float> invOut; // 1/outdeg, 0 on a dangling node
std::vector<uint32_t> dangling; // 1 when outdeg is 0
size_t n = 0;
size_t nnz = 0;
size_t danglingCount = 0;
};
// A counter mix, so element i is a pure function of i and the graph is the
// same on every machine and in every run. Not std::mt19937, whose output has
// changed between libstdc++ versions.
static uint64_t mix64(uint64_t z) {
z += 0x9e3779b97f4a7c15ull;
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ull;
z = (z ^ (z >> 27)) * 0x94d049bb133111ebull;
return z ^ (z >> 31);
}
// Builds an in-edge CSR graph with a skewed in-degree, which is the property
// that makes the gather in spmvCsrScalar scatter. Targets are drawn uniformly
// and then squared into range, so low-numbered nodes collect most of the
// citations, the way a real citation graph has hubs.
//
// This is not day 40's OpenAlex graph. It has the same node count and a
// similar edge count so the timeline has the same shape, and it needs no data
// file, which is what lets this program run anywhere.
static Graph buildGraph(size_t n) {
Graph g;
g.n = n;
std::vector<uint32_t> outDeg(n, 0);
std::vector<uint32_t> src;
std::vector<uint32_t> dst;
for (size_t u = 0; u < n; ++u) {
if (mix64(u) % 100u < 8u) {
continue; // cites nothing in the set: a dangling node
}
const uint32_t deg = 1u + static_cast<uint32_t>(mix64(2 * u + 1) % 24u);
for (uint32_t k = 0; k < deg; ++k) {
const uint64_t r = mix64(u * 1000003ull + k) % n;
size_t v = static_cast<size_t>((r * r) / n);
if (v == u) {
v = (v + 1) % n;
}
src.push_back(static_cast<uint32_t>(u));
dst.push_back(static_cast<uint32_t>(v));
++outDeg[u];
}
}
g.nnz = src.size();
g.rowPtr.assign(n + 1, 0);
for (size_t e = 0; e < g.nnz; ++e) {
++g.rowPtr[dst[e] + 1];
}
for (size_t v = 0; v < n; ++v) {
g.rowPtr[v + 1] += g.rowPtr[v];
}
std::vector<uint32_t> cursor(g.rowPtr.begin(), g.rowPtr.end() - 1);
g.colIdx.assign(g.nnz, 0);
for (size_t e = 0; e < g.nnz; ++e) {
g.colIdx[cursor[dst[e]]++] = src[e];
}
// invOut is a float, and the CPU reference reads this same float, so both
// sides split a node's rank the same way. It also means the total is not
// conserved to double precision: outdeg copies of float(1/outdeg) do not
// add back to 1, so the sum settles about 1e-7 away from 1 rather than
// 1e-16. That is why the gate below is 4 * FLT_EPSILON * sqrt(n) and not
// something tighter. Day 40 lands in the same place for the same reason.
g.invOut.assign(n, 0.0f);
g.dangling.assign(n, 0);
for (size_t u = 0; u < n; ++u) {
if (outDeg[u] == 0) {
g.dangling[u] = 1;
++g.danglingCount;
} else {
g.invOut[u] = 1.0f / static_cast<float>(outDeg[u]);
}
}
return g;
}
// The reference. Written for obvious correctness, not speed: plain loops, no
// blocking, and it accumulates in double where the kernels accumulate in
// float. A fixed iteration count, matching the GPU side.
static void pagerankCpu(const Graph& g, int iterations,
std::vector<double>* out) {
const size_t n = g.n;
const double invN = 1.0 / static_cast<double>(n);
std::vector<double> x(n, invN);
std::vector<double> next(n, 0.0);
std::vector<double> contrib(n, 0.0);
for (int it = 0; it < iterations; ++it) {
double mass = 0.0;
for (size_t u = 0; u < n; ++u) {
contrib[u] = x[u] * static_cast<double>(g.invOut[u]);
if (g.dangling[u] != 0u) {
mass += x[u];
}
}
for (size_t v = 0; v < n; ++v) {
double sum = 0.0;
for (uint32_t e = g.rowPtr[v]; e < g.rowPtr[v + 1]; ++e) {
sum += contrib[g.colIdx[e]];
}
next[v] = (1.0 - kDamping) * invN + kDamping * (sum + mass * invN);
}
x.swap(next);
}
*out = x;
}
// ----------------------------------------------------------- host: device
struct DeviceGraph {
uint32_t* rowPtr = nullptr;
uint32_t* colIdx = nullptr;
float* invOut = nullptr;
uint32_t* dangling = nullptr;
size_t n = 0;
};
struct Work {
float* x = nullptr;
float* xNext = nullptr;
float* contrib = nullptr;
float* y = nullptr;
float* partials = nullptr;
float* scalars = nullptr; // [0] dangling mass, [1] residual
};
// One PageRank iteration: five launches, each inside its own NVTX range.
// Without the ranges every one of these is a cudaLaunchKernel bar with the
// same name; with them the timeline reads back the phases of the algorithm.
//
// nvtxRangePushA and nvtxRangePop are header-only and cost close to nothing
// when no profiler is attached, so they stay in the shipped binary and the
// two captures on this page come from the same executable.
// snippet: iteration
static void iterate(const DeviceGraph& g, const Work& w, const float* x,
float* xNext) {
const int blocks = gridFor(g.n);
nvtxRangePushA("scale");
scaleByOutDegree<<<blocks, kThreadsPerBlock>>>(x, g.invOut, w.contrib, g.n);
nvtxRangePop();
nvtxRangePushA("spmv");
spmvCsrScalar<<<blocks, kThreadsPerBlock>>>(g.rowPtr, g.colIdx, w.contrib,
w.y, g.n);
nvtxRangePop();
nvtxRangePushA("dangling-mass");
sumDanglingPartial<<<kReduceBlocks, kThreadsPerBlock>>>(x, g.dangling,
w.partials, g.n);
sumPartials<<<1, kThreadsPerBlock>>>(w.partials, w.scalars, kReduceBlocks);
nvtxRangePop();
nvtxRangePushA("combine");
combineRanks<<<blocks, kThreadsPerBlock>>>(w.y, w.scalars, xNext, g.n,
kDamping);
nvtxRangePop();
}
// end snippet
// The convergence test, and the reason this page needs a timeline. Whether to
// stop is a host decision, so one float has to come back, and a synchronous
// device-to-host cudaMemcpy "returns only once the copy has completed", which
// means it waits for every kernel already queued.
//
// The value is computed and thrown away. That is deliberate: the three
// policies below differ only in how often this runs, so the rank vector they
// produce has to be bit-identical, and the program checks that it is.
// snippet: residual
static void residual(const Work& w, const float* x, const float* xPrev,
size_t n) {
nvtxRangePushA("residual");
sumAbsDiffPartial<<<kReduceBlocks, kThreadsPerBlock>>>(x, xPrev, w.partials,
n);
sumPartials<<<1, kThreadsPerBlock>>>(w.partials, w.scalars + 1,
kReduceBlocks);
nvtxRangePop();
float value = 0.0f;
nvtxRangePushA("copy-residual");
CUDA_CHECK(cudaMemcpy(&value, w.scalars + 1, sizeof(float),
cudaMemcpyDeviceToHost));
nvtxRangePop();
(void)value;
}
// end snippet
// Runs `iterations` iterations, reading the residual every `checkEvery` of
// them. checkEvery == 0 never reads it. Returns the buffer holding the answer,
// which is w.x or w.xNext depending on whether the count was even.
static float* runLoop(const DeviceGraph& g, const Work& w, int iterations,
int checkEvery) {
float* cur = w.x;
float* next = w.xNext;
for (int it = 0; it < iterations; ++it) {
nvtxRangePushA("iteration");
iterate(g, w, cur, next);
if (checkEvery > 0 && ((it + 1) % checkEvery) == 0) {
residual(w, next, cur, g.n);
}
nvtxRangePop();
float* tmp = cur;
cur = next;
next = tmp;
}
return cur;
}
static void resetRanks(const Work& w, size_t n) {
const float start = 1.0f / static_cast<float>(n);
fillFloat<<<gridFor(n), kThreadsPerBlock>>>(w.x, start, n);
fillFloat<<<gridFor(n), kThreadsPerBlock>>>(w.xNext, 0.0f, n);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
}
// Times one 50-iteration loop with events, after a warm-up loop that always
// reads the residual so every kernel in the program is warm whichever policy
// is timed next. Timed once rather than ten times: fifty iterations is not a
// microbenchmark, and ten repeats would bury the timeline this page is about
// under ten copies of itself.
static float timeLoop(const DeviceGraph& g, const Work& w, int checkEvery) {
resetRanks(w, g.n);
(void)runLoop(g, w, kWarmupIters, 1);
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaGetLastError());
resetRanks(w, g.n);
cudaEvent_t start;
cudaEvent_t stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
CUDA_CHECK(cudaEventRecord(start));
(void)runLoop(g, w, kTimedIters, checkEvery);
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;
}
// Times one launch with events and returns the mean milliseconds per run.
//
// This is the one function template and the one lambda this file uses. It
// warms up inside itself, which is what stops a warm-up going missing from
// one of six measurements.
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start;
cudaEvent_t stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
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;
}
// --------------------------------------------------------------- checking
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;
}
// Largest per-element error as a multiple of what the tolerance allows.
static double toleranceRatio(const std::vector<float>& got,
const std::vector<double>& want) {
double worst = 0.0;
for (size_t i = 0; i < got.size(); ++i) {
const double allowed = kRankAbsTol + kRankRelTol * std::fabs(want[i]);
const double err = std::fabs(static_cast<double>(got[i]) - want[i]);
if (err / allowed > worst) {
worst = err / allowed;
}
}
return worst;
}
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);
const Graph g = buildGraph(kNodes);
std::printf(
"synthetic citation graph: %zu nodes, %zu edges, %zu dangling\n\n", g.n,
g.nnz, g.danglingCount);
const size_t n = g.n;
const size_t floatBytes = n * sizeof(float);
const size_t flagBytes = n * sizeof(uint32_t);
DeviceGraph d;
d.n = n;
Work w;
CUDA_CHECK(cudaMalloc(&d.rowPtr, (n + 1) * sizeof(uint32_t)));
CUDA_CHECK(cudaMalloc(&d.colIdx, g.nnz * sizeof(uint32_t)));
CUDA_CHECK(cudaMalloc(&d.invOut, floatBytes));
CUDA_CHECK(cudaMalloc(&d.dangling, flagBytes));
CUDA_CHECK(cudaMalloc(&w.x, floatBytes));
CUDA_CHECK(cudaMalloc(&w.xNext, floatBytes));
CUDA_CHECK(cudaMalloc(&w.contrib, floatBytes));
CUDA_CHECK(cudaMalloc(&w.y, floatBytes));
CUDA_CHECK(cudaMalloc(&w.partials, kReduceBlocks * sizeof(float)));
CUDA_CHECK(cudaMalloc(&w.scalars, 2 * sizeof(float)));
nvtxRangePushA("upload");
CUDA_CHECK(cudaMemcpy(d.rowPtr, g.rowPtr.data(), (n + 1) * sizeof(uint32_t),
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d.colIdx, g.colIdx.data(), g.nnz * sizeof(uint32_t),
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d.invOut, g.invOut.data(), floatBytes,
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d.dangling, g.dangling.data(), flagBytes,
cudaMemcpyHostToDevice));
nvtxRangePop();
int status = EXIT_SUCCESS;
// Correctness. Both GPU policies and the CPU reference run the same fixed
// 100 iterations, because two implementations that stop at different
// points differ by more than their arithmetic does.
nvtxRangePushA("correctness");
std::vector<float> h_never(n);
std::vector<float> h_always(n);
resetRanks(w, n);
const float* answer = runLoop(d, w, kCheckIters, 0);
CUDA_CHECK(
cudaMemcpy(h_never.data(), answer, floatBytes, cudaMemcpyDeviceToHost));
resetRanks(w, n);
answer = runLoop(d, w, kCheckIters, 1);
CUDA_CHECK(cudaMemcpy(h_always.data(), answer, floatBytes,
cudaMemcpyDeviceToHost));
std::vector<double> h_ref;
pagerankCpu(g, kCheckIters, &h_ref);
nvtxRangePop();
const bool identical =
std::memcmp(h_never.data(), h_always.data(), floatBytes) == 0;
const double ratio = toleranceRatio(h_never, h_ref);
const double sumErr = sumOf(h_never) - 1.0;
const double sumAllowed = 4.0 * static_cast<double>(FLT_EPSILON) *
std::sqrt(static_cast<double>(n));
std::printf("correctness, %d fixed iterations everywhere\n", kCheckIters);
std::printf(" vs the double CPU reference err/tol = %.3f %s\n", ratio,
(ratio <= 1.0) ? "pass" : "FAIL");
std::printf(" never checked vs checked %s %s\n",
identical ? "bit-identical" : "DIFFERENT ",
identical ? "pass" : "FAIL");
std::printf(" sum of ranks minus 1 %+.3e (allowed %.3e) %s\n",
sumErr, sumAllowed,
(std::fabs(sumErr) <= sumAllowed) ? "pass" : "FAIL");
if (ratio > 1.0 || !identical || std::fabs(sumErr) > sumAllowed) {
std::fprintf(stderr, "correctness failed; not timing anything\n");
status = EXIT_FAILURE;
}
if (status == EXIT_SUCCESS) {
const int blocks = gridFor(n);
resetRanks(w, n);
nvtxRangePushA("per-kernel-timing");
const float msScale = timeKernel([&] {
scaleByOutDegree<<<blocks, kThreadsPerBlock>>>(w.x, d.invOut,
w.contrib, n);
});
const float msSpmv = timeKernel([&] {
spmvCsrScalar<<<blocks, kThreadsPerBlock>>>(d.rowPtr, d.colIdx,
w.contrib, w.y, n);
});
const float msDangling = timeKernel([&] {
sumDanglingPartial<<<kReduceBlocks, kThreadsPerBlock>>>(
w.x, d.dangling, w.partials, n);
});
const float msDiff = timeKernel([&] {
sumAbsDiffPartial<<<kReduceBlocks, kThreadsPerBlock>>>(
w.x, w.xNext, w.partials, n);
});
const float msPartials = timeKernel([&] {
sumPartials<<<1, kThreadsPerBlock>>>(w.partials, w.scalars,
kReduceBlocks);
});
const float msCombine = timeKernel([&] {
combineRanks<<<blocks, kThreadsPerBlock>>>(w.y, w.scalars, w.xNext,
n, kDamping);
});
nvtxRangePop();
std::printf("\none launch each, mean of %d runs after %d warm-ups\n",
kTimedRuns, kWarmupRuns);
std::printf("kernel ms\n");
std::printf("-------------------- ----------\n");
std::printf("scaleByOutDegree %10.4f\n",
static_cast<double>(msScale));
std::printf("spmvCsrScalar %10.4f\n",
static_cast<double>(msSpmv));
std::printf("sumDanglingPartial %10.4f\n",
static_cast<double>(msDangling));
std::printf("sumAbsDiffPartial %10.4f\n",
static_cast<double>(msDiff));
std::printf("sumPartials %10.4f\n",
static_cast<double>(msPartials));
std::printf("combineRanks %10.4f\n",
static_cast<double>(msCombine));
std::printf("one iteration, summed %10.4f\n",
static_cast<double>(msScale + msSpmv + msDangling +
msPartials + msCombine));
nvtxRangePushA("policy-never");
const float msNever = timeLoop(d, w, 0);
nvtxRangePop();
nvtxRangePushA("policy-every-10");
const float msTen = timeLoop(d, w, 10);
nvtxRangePop();
nvtxRangePushA("policy-every-step");
const float msEvery = timeLoop(d, w, 1);
nvtxRangePop();
std::printf(
"\n%d iterations, the same arithmetic. Only the convergence\n"
"test moves, and it reads four bytes when it runs.\n",
kTimedIters);
std::printf("policy total (ms) per iter (ms) D2H\n");
std::printf("-------------------- ------------ ------------- ---\n");
std::printf("never checked %12.4f %13.4f %3d\n",
static_cast<double>(msNever),
static_cast<double>(msNever) / kTimedIters, 0);
std::printf("checked every 10 %12.4f %13.4f %3d\n",
static_cast<double>(msTen),
static_cast<double>(msTen) / kTimedIters, kTimedIters / 10);
std::printf("checked every step %12.4f %13.4f %3d\n",
static_cast<double>(msEvery),
static_cast<double>(msEvery) / kTimedIters, kTimedIters);
std::printf(
"checked every step costs %.2fx never checked\n",
static_cast<double>(msEvery) / static_cast<double>(msNever));
std::printf(
"\nEvery phase above sits inside an NVTX range. Profile with:\n"
" nsys profile -t cuda -o day41-plain ./pagerank_nvtx\n"
" nsys profile -t cuda,nvtx -o day41-nvtx ./pagerank_nvtx\n"
"The first gives you cudaLaunchKernel. The second gives you\n"
"scale, spmv, dangling-mass, combine, residual, copy-residual.\n");
}
CUDA_CHECK(cudaFree(d.rowPtr));
CUDA_CHECK(cudaFree(d.colIdx));
CUDA_CHECK(cudaFree(d.invOut));
CUDA_CHECK(cudaFree(d.dangling));
CUDA_CHECK(cudaFree(w.x));
CUDA_CHECK(cudaFree(w.xNext));
CUDA_CHECK(cudaFree(w.contrib));
CUDA_CHECK(cudaFree(w.y));
CUDA_CHECK(cudaFree(w.partials));
CUDA_CHECK(cudaFree(w.scalars));
if (status == EXIT_SUCCESS) {
std::printf("\nevery case passed\n");
}
return status;
}