code/day49-roofline/roofline.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 49: a roofline built from this card's own two ceilings, with ten
// kernels on it and two byte counts for each of them.
//
// Day 30 built the analytical half of this and shipped a documented flaw:
// it divides COMPULSORY bytes, the traffic an algorithm asks global memory
// for, by elapsed time and calls the result GB/s. On a kernel with reuse
// that is a rate the memory system never moved. Its naive matmul reported
// 1008.5 percent of the same card's copy ceiling, which is impossible and
// which the evidence file says so about.
//
// This program is the other half. It measures the two ceilings and the ten
// times, and prints the compulsory byte count beside each row. The measured
// byte count comes from Nsight Compute, which is the only thing on the
// machine that knows what reached DRAM:
//
// sudo ncu --csv --metrics dram__bytes.sum,lts__t_bytes.sum,\
// l1tex__t_bytes.sum ./roofline --profile > profile/day49-bytes.csv
// python3 merge_dram.py --json run.jsonl --ncu profile/day49-bytes.csv
//
// `ncu` needs performance counters and a stock driver ships
// RmProfilingAdminOnly = 1, so a plain run returns ERR_NVGPUCTRPERM and
// `sudo ncu` is the working form. Colab and Kaggle give you neither root
// nor that parameter, which is why the report and both CSV exports ship in
// profile/ and merge_dram.py runs anywhere Python does. README.md has the
// full commands.
//
// Three modes, and argv carries nothing else:
// ./roofline the human table
// ./roofline --json one JSON object per row, for scripts/bench_ingest
// ./roofline --profile one launch per kernel and nothing else, so an
// ncu pass sees twelve launches instead of 168
//
// No vendor number appears anywhere in this file, for the reason day 30
// gives: a ceiling nobody can reach makes every percentage measured against
// it wrong by the same unknown factor.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo -o roofline roofline.cu
// Run: ./roofline
//
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <vector>
#include <cuda_runtime.h>
// The one error macro. This file is standalone, the way a Compiler Explorer
// embed is, so it carries its own verbatim copy. `err_` carries a trailing
// underscore so it cannot collide with a variable at the call site, and the
// do/while makes the macro one statement so it survives a braceless `if`.
#define CUDA_CHECK(call) \
do { \
cudaError_t err_ = (call); \
if (err_ != cudaSuccess) { \
std::fprintf(stderr, "CUDA error %s:%d: %s: %s\n", __FILE__, \
__LINE__, #call, cudaGetErrorString(err_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
// 4096 squared, because the transpose has to be square and every other
// vector kernel shares its buffers. 16,777,216 floats is 64 MiB per buffer
// and three buffers, so the program needs 192 MiB of device memory.
//
// That count is also 2^24, the last integer a float represents exactly. The
// transpose check fills the input with its own index and compares exactly,
// which only holds while every index survives the trip through a float.
constexpr size_t kSide = 4096;
constexpr size_t kElems = kSide * kSide;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kTransposeDim = 32; // one warp per row of the block
// The strided vector add walks the same permutation day 11 settled on: a 2D
// grid where x covers kElems / kStride columns and y covers kStride rows.
// Every index in [0, kElems) is visited exactly once at every stride, so the
// coalesced and strided rows ask global memory for identical byte counts and
// differ only in the address pattern. `(i * stride) % n` is not a bijection
// for a power-of-two n, which is the bug day 11 shipped and then fixed.
constexpr int kStride = 32;
// The matmul is square and divides the tile exactly, so neither kernel needs
// an edge guard. Day 16 owns the ragged case.
//
// 512 squared floats is 1 MiB per matrix and three matrices is 3 MiB, which
// fits inside this card's 4 MiB L2. That is not an accident of sizing, it is
// the whole reason day 30's naive row read 1008 percent: almost none of the
// traffic the model charges to DRAM reaches DRAM. The program prints the L2
// size so a reader on a different card can see whether the same is true
// there.
constexpr size_t kMatDim = 512;
constexpr int kTileDim = 16;
// The stencil, taken from day 34. A 2048 squared grid is 16 MiB, and the two
// buffers it needs are 32 MiB, which does not fit in a 4 MiB L2. The pair is
// here because day 34 measured tiling making this kernel slower and argued
// that the reads the tile removes were already cache hits. That argument was
// an inference from a timer. Two DRAM byte counts settle it.
constexpr int kGridDim = 2048;
constexpr int kHalo = 1;
constexpr int kStencilX = 32; // one warp covers one row of the tile
constexpr int kStencilY = 8;
constexpr int kSharedX = kStencilX + 2 * kHalo; // 34
constexpr int kSharedY = kStencilY + 2 * kHalo; // 10
constexpr int kSharedCells = kSharedX * kSharedY;
constexpr size_t kStencilCells =
static_cast<size_t>(kGridDim) * static_cast<size_t>(kGridDim);
constexpr size_t kStencilInterior =
static_cast<size_t>(kGridDim - 2) * static_cast<size_t>(kGridDim - 2);
constexpr float kStencilR = 0.2f; // under the 1/4 CFL bound of day 34
constexpr float kCold = 0.0f;
// The one kernel on the chart that sits right of the ridge point. A Horner
// chain of kPolyDegree fused multiply-adds over one loaded float is
// 2 * kPolyDegree + 1 flops for 8 bytes, so its intensity is in the tens
// while every other row here is under 4. A roofline whose points are all on
// the sloped half teaches only half the model.
//
// kPolyScale keeps |x| below one so the chain cannot run away to inf: the
// inputs are small integers in [1, 7] and 7 * 0.0625 is 0.4375.
constexpr int kPolyDegree = 256;
constexpr float kPolyScale = 0.0625f;
constexpr float kPolyAdd = 1.0f;
constexpr size_t kPolyCheck = 65536;
// The FP32 ceiling, unchanged from day 30. Eight independent chains per
// thread so an FMA never waits on the one before it.
//
// kFmaBlocksLaunchedPerSm is a launch multiplier, not a residency figure. On
// sm_75 an SM holds 16 blocks and 1024 threads at once, so at 256 threads a
// block the thread cap binds first and only 4 are ever resident together.
// Launching 32 per SM queues eight waves behind those 4, which makes the
// tail an eighth of the run rather than all of it.
constexpr int kFmaChains = 8;
constexpr int kFmaIters = 2048;
constexpr int kFmaBlocksLaunchedPerSm = 32;
constexpr float kFmaMul = 1.0000001f;
constexpr float kFmaAdd = 1.0f;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kRelTolerance = 1e-5f;
// A compiler that deleted the FMA loop would report a peak the card cannot
// reach and every point on the chart would be measured against it. One
// percent is loose on purpose: this is a liveness check on the loop, not an
// accuracy claim. The value grows with the iteration count, so a loop that
// ran half the steps comes back half the size and fails by a factor of two.
constexpr float kFmaTolerance = 0.01f;
constexpr int kNumRows = 10;
constexpr int kNumTimings = 12; // ten placed kernels plus the two ceilings
// Blocks for a one-thread-one-element launch over n elements. constexpr
// because the static_asserts below call it; they follow only from the
// constants in this file and never ask the device anything.
constexpr size_t blocksFor(size_t n) {
return (n + kThreadsPerBlock - 1) / kThreadsPerBlock;
}
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert((kThreadsPerBlock & (kThreadsPerBlock - 1)) == 0,
"the reduction halves the stride, so the block size must be a "
"power of two");
static_assert(kElems % kThreadsPerBlock == 0,
"the reduction's add count is exact only when every block is "
"full");
static_assert(kElems <= (1u << 24),
"the transpose check compares floats filled with their own "
"index, which stops being exact above 2^24");
static_assert(kElems % kStride == 0,
"the strided vector add's 2D grid covers kElems / kStride "
"columns by kStride rows, which is a permutation only when the "
"stride divides the element count");
static_assert(kSide % kTransposeDim == 0,
"both transposes tile the matrix exactly, which is why neither "
"carries a bounds check");
static_assert(kTransposeDim * kTransposeDim <= 1024,
"a transpose block is at most 1024 threads");
static_assert(kTransposeDim * (kTransposeDim + 1) * sizeof(float) <= 49152,
"the padded transpose tile must fit the 48 KiB a block gets by "
"default on sm_75, with no cudaFuncSetAttribute opt-in");
static_assert(kMatDim % kTileDim == 0,
"the tiled matmul here has no edge guard; day 16 is the lesson "
"that adds one");
static_assert(kMatDim * kMatDim <= kElems,
"the matmul reuses the vector buffers, so it has to fit in one");
static_assert(kMatDim * 7 * 2 < (1u << 24),
"the largest dot product must be exact in a float, or a "
"mismatch below could be rounding rather than an index bug");
static_assert(kStencilX * kStencilY == kThreadsPerBlock,
"the stencil's fill loop strides by the block size, and the "
"output tile is one cell per thread");
static_assert(kSharedCells * sizeof(float) <= 49152,
"the stencil apron must fit the 48 KiB a block gets by default");
static_assert(kStencilCells <= kElems,
"the stencil borrows the vector buffers");
static_assert(kPolyCheck <= kElems,
"the polynomial check reads a prefix of the output buffer");
static_assert(blocksFor(kElems) <= 2147483647u,
"gridDim.x is bounded at 2^31 - 1");
// out[i] = in[i]. One thread owns one element.
//
// Memory: consecutive threads take consecutive elements, so one warp's 32
// addresses cover 128 contiguous bytes on the read and 128 on the write.
// That is the best address pattern this card has, which is what makes this
// kernel the bandwidth ceiling rather than another row in the table.
//
// Launch assumption: gridDim.x * blockDim.x >= n. The guard never fires at
// this program's size, because kElems divides the block size, and it stays
// anyway: it is the shape every kernel from day 5 onward has.
__global__ void copyCoalesced(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
out[i] = in[i];
}
}
// out[i] = a[i] + b[i]. One thread owns one element. Day 5's kernel.
//
// Memory: three coalesced streams, two in and one out, so 12 bytes cross the
// bus for one add. That ratio is the whole reason this kernel is on the
// chart.
//
// Launch assumption: gridDim.x * blockDim.x >= n.
__global__ void vectorAdd(const float* __restrict__ a,
const float* __restrict__ b, float* __restrict__ out,
size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
out[i] = a[i] + b[i];
}
}
// The same add over the same elements, reached in a different order.
//
// Memory: a warp's 32 lanes hold addresses `stride` floats apart. At stride
// 32 that is 128 bytes per lane, so each lane owns a 32-byte sector of its
// own and a warp's load costs 32 sectors instead of 4. The compulsory byte
// count does not change by one byte, because this is a permutation of the
// same index set, which is exactly what makes the pair worth profiling.
//
// Launch assumption: a 2D grid, x over m = n / stride columns and y over
// `stride` rows, so col in [0, m) and blockIdx.y in [0, stride) covers every
// index in [0, n) exactly once.
// snippet: strided-add
__global__ void vectorAddPermuted(const float* __restrict__ a,
const float* __restrict__ b,
float* __restrict__ out, size_t m,
int stride) {
const size_t col =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (col < m) {
const size_t j = col * static_cast<size_t>(stride) + blockIdx.y;
out[j] = a[j] + b[j];
}
}
// end snippet: strided-add
// out[col * side + row] = in[row * side + col]. One thread owns one element
// and does no arithmetic on it. Day 12's kernel.
//
// Memory: the read is coalesced, 32 lanes over 128 contiguous bytes. The
// write is not: consecutive lanes write addresses `side` floats apart, so
// one warp's store touches 32 separate sectors. Zero flops either way, which
// puts this kernel and the next one at the same place on the chart.
//
// Launch assumption: a kTransposeDim square block over a grid that tiles the
// matrix exactly, which the static_assert above guarantees.
__global__ void transposeNaive(const float* __restrict__ in,
float* __restrict__ out, size_t side) {
const size_t row =
blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
const size_t col =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
out[col * side + row] = in[row * side + col];
}
// The same transpose through a shared tile. Day 13's kernel.
//
// Memory: both the read and the write are now coalesced, because the
// transposing step happens inside shared memory. The tile is padded to 33
// columns so that the column read below lands in 32 different banks; day 15
// is the lesson that measures what the missing pad costs.
//
// Launch assumption: the same exact tiling as transposeNaive, and every
// thread reaches the barrier.
__global__ void transposeTiled(const float* __restrict__ in,
float* __restrict__ out, size_t side) {
__shared__ float tile[kTransposeDim][kTransposeDim + 1];
const size_t inCol =
blockIdx.x * static_cast<size_t>(kTransposeDim) + threadIdx.x;
const size_t inRow =
blockIdx.y * static_cast<size_t>(kTransposeDim) + threadIdx.y;
tile[threadIdx.y][threadIdx.x] = in[inRow * side + inCol];
__syncthreads();
const size_t outCol =
blockIdx.y * static_cast<size_t>(kTransposeDim) + threadIdx.x;
const size_t outRow =
blockIdx.x * static_cast<size_t>(kTransposeDim) + threadIdx.y;
out[outRow * side + outCol] = tile[threadIdx.x][threadIdx.y];
}
// Sums one block's slice of `in` into out[blockIdx.x]. Day 24's version 3,
// sequential addressing, copied unchanged so the point on the chart belongs
// to a kernel the reader has already written.
//
// Memory: one coalesced read per thread, then everything happens in shared
// memory, then one float per block goes back to global. Four bytes read per
// add is the best a streaming kernel does, and it is still far left of the
// ridge point.
//
// Launch assumption: exactly kThreadsPerBlock threads per block, a power of
// two, because the shared array is sized from that constant and the halving
// loop counts down from it. Every thread reaches every barrier: the guard
// covers the load, not the __syncthreads(), and nothing returns above one.
__global__ void reduceSequentialAddressing(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const unsigned int width = kThreadsPerBlock;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
tile[tid] = (i < n) ? in[i] : 0.0f;
__syncthreads();
for (unsigned int s = width / 2; s > 0; s >>= 1) {
if (tid < s) {
tile[tid] += tile[tid + s];
}
__syncthreads();
}
if (tid == 0) {
out[blockIdx.x] = tile[0];
}
}
// c = a * b, square, one thread per output element. Day 16's baseline.
//
// Memory: every thread walks a whole row of `a` and a whole column of `b`,
// so it issues 2 * d loads for 2 * d flops. Consecutive lanes share one
// element of `a` and take consecutive elements of one row of `b`, so the
// address patterns are fine. The count is the problem: the same row of `a`
// is fetched again by every thread in the block, and at this size the whole
// working set is smaller than L2, so almost none of those fetches reach
// DRAM.
//
// Launch assumption: a kTileDim square block over a grid that tiles the
// matrix exactly, so there is no edge guard.
__global__ void matmulNaive(const float* __restrict__ a,
const float* __restrict__ b, float* __restrict__ c,
size_t d) {
const size_t row =
blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
const size_t col =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
float acc = 0.0f;
for (size_t p = 0; p < d; ++p) {
acc += a[row * d + p] * b[p * d + col];
}
c[row * d + col] = acc;
}
// The same product, one kTileDim square tile at a time. Day 16's kernel.
//
// Memory: each thread loads two floats per tile step instead of two per k
// step, so the compulsory load count falls by kTileDim and the flop count
// does not move. This is the only change in the file that moves a kernel
// right on the chart rather than up towards a ceiling.
//
// Launch assumption: the block is exactly kTileDim by kTileDim and the grid
// tiles the matrix exactly. Every thread reaches both barriers.
__global__ void matmulTiled(const float* __restrict__ a,
const float* __restrict__ b, float* __restrict__ c,
size_t d) {
__shared__ float tileA[kTileDim][kTileDim];
__shared__ float tileB[kTileDim][kTileDim];
const size_t row = blockIdx.y * static_cast<size_t>(kTileDim) + threadIdx.y;
const size_t col = blockIdx.x * static_cast<size_t>(kTileDim) + threadIdx.x;
float acc = 0.0f;
for (size_t step = 0; step < d / kTileDim; ++step) {
const size_t aCol = step * kTileDim + threadIdx.x;
const size_t bRow = step * kTileDim + threadIdx.y;
tileA[threadIdx.y][threadIdx.x] = a[row * d + aCol];
tileB[threadIdx.y][threadIdx.x] = b[bRow * d + col];
__syncthreads();
for (int p = 0; p < kTileDim; ++p) {
acc += tileA[threadIdx.y][p] * tileB[p][threadIdx.x];
}
__syncthreads();
}
c[row * d + col] = acc;
}
// One explicit step of the 2D heat equation, straight from global memory.
// Day 34's kernel.
//
// Memory: five reads and one write per interior cell. A warp holds 32
// consecutive columns of one row, so each of the five loads is a contiguous
// 128-byte span and coalesces into four sectors. Four of the five are cells
// a neighbouring thread also wants, which is the reuse the next kernel tries
// to capture.
//
// Launch assumption: a kStencilX by kStencilY block over a grid that covers
// the interior. The guard covers the store; there is no barrier here.
__global__ void heatStepGlobal(const float* __restrict__ in,
float* __restrict__ out, int n, float r) {
const int col = static_cast<int>(blockIdx.x * blockDim.x + threadIdx.x) + 1;
const int row = static_cast<int>(blockIdx.y * blockDim.y + threadIdx.y) + 1;
if (col < n - 1 && row < n - 1) {
const size_t stride = static_cast<size_t>(n);
const size_t c = static_cast<size_t>(row) * stride + col;
const float centre = in[c];
out[c] = centre + r * (in[c - 1] + in[c + 1] + in[c - stride] +
in[c + stride] - 4.0f * centre);
}
}
// The same step through a shared apron. Day 34's kernel.
//
// Memory: one block loads a 34 by 10 patch to update a 32 by 8 patch, so it
// asks global memory for 340 cells to produce 256. That is 1.33 reads per
// updated cell against five, and day 34 measured it running slower anyway.
// The fill loop strides by the block size because 256 threads have 340 cells
// to fetch, so 84 of them go round twice.
//
// Launch assumption: the same grid as heatStepGlobal, and every thread
// reaches the barrier. The guard covers the load and the store, never the
// __syncthreads().
__global__ void heatStepTiled(const float* __restrict__ in,
float* __restrict__ out, int n, float r) {
__shared__ float tile[kSharedY][kSharedX];
const int tileCol = static_cast<int>(blockIdx.x) * kStencilX + 1;
const int tileRow = static_cast<int>(blockIdx.y) * kStencilY + 1;
const int threads = static_cast<int>(blockDim.x * blockDim.y);
const int flat = static_cast<int>(threadIdx.y * blockDim.x + threadIdx.x);
for (int s = flat; s < kSharedCells; s += threads) {
const int sy = s / kSharedX;
const int sx = s - sy * kSharedX;
const int gy = tileRow - kHalo + sy;
const int gx = tileCol - kHalo + sx;
const bool inside = (gx >= 0 && gx < n && gy >= 0 && gy < n);
tile[sy][sx] =
inside ? in[static_cast<size_t>(gy) * static_cast<size_t>(n) + gx]
: kCold;
}
__syncthreads();
const int sx = static_cast<int>(threadIdx.x) + kHalo;
const int sy = static_cast<int>(threadIdx.y) + kHalo;
const int col = tileCol + static_cast<int>(threadIdx.x);
const int row = tileRow + static_cast<int>(threadIdx.y);
if (col < n - 1 && row < n - 1) {
const float centre = tile[sy][sx];
out[static_cast<size_t>(row) * static_cast<size_t>(n) + col] =
centre + r * (tile[sy][sx - 1] + tile[sy][sx + 1] +
tile[sy - 1][sx] + tile[sy + 1][sx] - 4.0f * centre);
}
}
// A Horner chain over one loaded float. One thread owns one element.
//
// Memory: two coalesced streams, one in and one out, so eight bytes buy
// 2 * kPolyDegree + 1 flops. That is the only intensity in this program
// large enough to land right of the ridge point, and it is here so the flat
// half of the roofline has a point on it.
//
// The loop count is a compile-time constant, unlike fmaPeak's, because 256
// fused multiply-adds unroll into a reasonable amount of code while 16,384
// would measure the instruction cache. The chain is dependent, so each FMA
// waits on the one before it; the parallelism comes from having 16,777,216
// elements rather than from inside a thread.
//
// Launch assumption: gridDim.x * blockDim.x >= n.
// snippet: polynomial
__global__ void polynomialHorner(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
const float x = in[i] * kPolyScale;
float acc = x;
#pragma unroll 16
for (int k = 0; k < kPolyDegree; ++k) {
acc = fmaf(acc, x, kPolyAdd);
}
out[i] = acc;
}
}
// end snippet: polynomial
// The FP32 ceiling: fused multiply-adds and nothing else. Day 30's kernel,
// unchanged, so the two lessons' ceilings are comparable.
//
// Memory: one coalesced store per thread, after the loop. The loop reads and
// writes registers only, so nothing in it is limited by bandwidth.
//
// Launch assumption: `out` holds one float per thread. `iters` arrives as a
// run-time argument on purpose. A constexpr count would let nvcc unroll all
// 2048 steps into straight-line code, and 16,384 instructions of it would
// measure the instruction cache instead of the FP32 pipelines.
__global__ void fmaPeak(float* __restrict__ out, int iters) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
float acc[kFmaChains];
#pragma unroll
for (int q = 0; q < kFmaChains; ++q) {
acc[q] = static_cast<float>(q) * 0.125f;
}
for (int k = 0; k < iters; ++k) {
#pragma unroll
for (int q = 0; q < kFmaChains; ++q) {
acc[q] = fmaf(acc[q], kFmaMul, kFmaAdd);
}
}
float sum = 0.0f;
#pragma unroll
for (int q = 0; q < kFmaChains; ++q) {
sum += acc[q];
}
out[i] = sum;
}
// CPU references. Written for obvious correctness, not speed: plain loops,
// no OpenMP, no intrinsics. None of them allocates; the caller owns every
// buffer.
static void vectorAddCpu(const float* a, const float* b, float* out, size_t n) {
for (size_t i = 0; i < n; ++i) {
out[i] = a[i] + b[i];
}
}
// Accumulates in double even though the kernel accumulates in float, because
// a reference exists to be right rather than to match bit for bit.
static double sumCpu(const float* x, size_t n) {
double total = 0.0;
for (size_t i = 0; i < n; ++i) {
total += static_cast<double>(x[i]);
}
return total;
}
static void matmulCpu(const float* a, const float* b, float* c, size_t d) {
for (size_t row = 0; row < d; ++row) {
for (size_t col = 0; col < d; ++col) {
double acc = 0.0;
for (size_t p = 0; p < d; ++p) {
acc += static_cast<double>(a[row * d + p]) *
static_cast<double>(b[p * d + col]);
}
c[row * d + col] = static_cast<float>(acc);
}
}
}
// One explicit heat step in double, with the boundary ring written as zero
// so the comparison can run over the whole buffer. Both device kernels leave
// the ring untouched, and the caller zeroes the device output before each of
// them, so a kernel that skipped interior cells fails here rather than
// passing on a buffer the previous kernel already filled.
static void heatStepCpu(const float* in, float* out, int n, float r) {
const size_t side = static_cast<size_t>(n);
for (size_t i = 0; i < side * side; ++i) {
out[i] = 0.0f;
}
for (int row = 1; row < n - 1; ++row) {
for (int col = 1; col < n - 1; ++col) {
const size_t c = static_cast<size_t>(row) * side + col;
const double centre = static_cast<double>(in[c]);
const double neighbours = static_cast<double>(in[c - 1]) +
static_cast<double>(in[c + 1]) +
static_cast<double>(in[c - side]) +
static_cast<double>(in[c + side]);
out[c] = static_cast<float>(
centre + static_cast<double>(r) * (neighbours - 4.0 * centre));
}
}
}
// The Horner chain in float, in the same order the kernel runs it, so the
// two agree to within rounding rather than to within an argument about
// association.
static void polynomialCpu(const float* in, float* out, size_t n) {
for (size_t i = 0; i < n; ++i) {
const float x = in[i] * kPolyScale;
float acc = x;
for (int k = 0; k < kPolyDegree; ++k) {
acc = std::fmaf(acc, x, kPolyAdd);
}
out[i] = acc;
}
}
// One thread's worth of the FMA loop, in the same order the kernel runs it.
// This is what proves the loop was not optimised away.
static float fmaSumCpu(int iters) {
float sum = 0.0f;
for (int q = 0; q < kFmaChains; ++q) {
float acc = static_cast<float>(q) * 0.125f;
for (int k = 0; k < iters; ++k) {
acc = std::fmaf(acc, kFmaMul, kFmaAdd);
}
sum += acc;
}
return sum;
}
// Returns the first index where got and want differ by more than the
// relative tolerance, or n if they agree everywhere. Returning the index
// rather than a bool is the point: "wrong at 512" names the block, "wrong"
// does not.
static size_t firstMismatch(const float* got, const float* want, size_t n,
float relTolerance) {
for (size_t i = 0; i < n; ++i) {
const float scale = (want[i] == 0.0f) ? 1.0f : std::fabs(want[i]);
if (std::fabs(got[i] - want[i]) > relTolerance * scale) {
return i;
}
}
return n;
}
// Returns the first index holding a value that is not finite, or n. The
// polynomial output is only compared against the host over a prefix, because
// the reference is 4.3 billion host multiply-adds at full size. This scan is
// what covers the rest: the buffer is filled with all-ones bytes first, and
// an all-ones float is a NaN, so any element the kernel did not write fails
// here.
static size_t firstNonFinite(const float* x, size_t n) {
for (size_t i = 0; i < n; ++i) {
if (!std::isfinite(x[i])) {
return i;
}
}
return n;
}
// Times one launch with CUDA events. A host-side clock around a launch
// measures the launch, not the kernel, because launches are asynchronous.
// See day 9.
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
// Warm up this kernel, not just the first kernel in the program. Lazy
// module loading has been the default since CUDA 12.2 on Linux, so the
// first launch of each kernel pays its own load.
for (int i = 0; i < kWarmupRuns; ++i) {
launch();
}
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < kTimedRuns; ++i) {
launch();
}
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError());
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / kTimedRuns;
}
int main(int argc, char** argv) {
// argv carries two flags and nothing else. Every size in this program is
// a constant at the top of the file, so the page and the program cannot
// disagree about what was measured.
const bool json = (argc > 1 && std::strcmp(argv[1], "--json") == 0);
const bool profile = (argc > 1 && std::strcmp(argv[1], "--profile") == 0);
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
const size_t bytes = kElems * sizeof(float);
const size_t vecBlocks = blocksFor(kElems);
const int vecGrid = static_cast<int>(vecBlocks);
const size_t matElems = kMatDim * kMatDim;
const size_t strideCols = kElems / static_cast<size_t>(kStride);
const int strideGridX = static_cast<int>(blocksFor(strideCols));
const int fmaBlocks = kFmaBlocksLaunchedPerSm * prop.multiProcessorCount;
const size_t fmaThreads = static_cast<size_t>(fmaBlocks) * kThreadsPerBlock;
// The stencil grid covers the interior, so it is sized from kGridDim - 2
// rather than kGridDim, and it does not divide evenly on either axis at
// 2048. That is why both stencil kernels carry a guard on the store.
const int stencilBlocksX = (kGridDim - 2 + kStencilX - 1) / kStencilX;
const int stencilBlocksY = (kGridDim - 2 + kStencilY - 1) / kStencilY;
// fmaPeak writes one float per thread into d_out, which holds kElems of
// them, and its grid is sized from the device's SM count. On a card far
// larger than anything this course targets that write would run off the
// end, so it is checked here, before any allocation exists to leak.
if (fmaThreads > kElems) {
std::fprintf(stderr,
"%d SMs asks for %zu FMA threads, more than the %zu "
"floats d_out holds; lower kFmaBlocksLaunchedPerSm\n",
prop.multiProcessorCount, fmaThreads, kElems);
return EXIT_FAILURE;
}
std::vector<float> h_a(kElems);
std::vector<float> h_b(kElems);
std::vector<float> h_out(kElems);
std::vector<float> h_want(kElems);
float* d_a = nullptr;
float* d_b = nullptr;
float* d_out = nullptr;
float* d_partials = nullptr;
CUDA_CHECK(cudaMalloc(&d_a, bytes));
CUDA_CHECK(cudaMalloc(&d_b, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
CUDA_CHECK(cudaMalloc(&d_partials, vecBlocks * sizeof(float)));
const dim3 transposeBlock(kTransposeDim, kTransposeDim);
const dim3 transposeGrid(static_cast<unsigned int>(kSide / kTransposeDim),
static_cast<unsigned int>(kSide / kTransposeDim));
const dim3 matBlock(kTileDim, kTileDim);
const dim3 matGrid(static_cast<unsigned int>(kMatDim / kTileDim),
static_cast<unsigned int>(kMatDim / kTileDim));
const dim3 strideGrid(static_cast<unsigned int>(strideGridX),
static_cast<unsigned int>(kStride));
const dim3 stencilBlock(kStencilX, kStencilY);
const dim3 stencilGrid(static_cast<unsigned int>(stencilBlocksX),
static_cast<unsigned int>(stencilBlocksY));
// Every launch in the program is written once, here, and used twice: by
// the checking and timing path below and by --profile. Two copies of a
// launch configuration is how a profiled kernel ends up running a
// different shape from the timed one, which would make the two byte
// columns describe two different programs.
auto launchCopy = [&] {
copyCoalesced<<<vecGrid, kThreadsPerBlock>>>(d_a, d_out, kElems);
};
auto launchVectorAdd = [&] {
vectorAdd<<<vecGrid, kThreadsPerBlock>>>(d_a, d_b, d_out, kElems);
};
auto launchVectorAddPermuted = [&] {
vectorAddPermuted<<<strideGrid, kThreadsPerBlock>>>(
d_a, d_b, d_out, strideCols, kStride);
};
auto launchTransposeNaive = [&] {
transposeNaive<<<transposeGrid, transposeBlock>>>(d_a, d_out, kSide);
};
auto launchTransposeTiled = [&] {
transposeTiled<<<transposeGrid, transposeBlock>>>(d_a, d_out, kSide);
};
auto launchMatmulNaive = [&] {
matmulNaive<<<matGrid, matBlock>>>(d_a, d_b, d_out, kMatDim);
};
auto launchMatmulTiled = [&] {
matmulTiled<<<matGrid, matBlock>>>(d_a, d_b, d_out, kMatDim);
};
auto launchReduce = [&] {
reduceSequentialAddressing<<<vecGrid, kThreadsPerBlock>>>(
d_a, d_partials, kElems);
};
auto launchStencilGlobal = [&] {
heatStepGlobal<<<stencilGrid, stencilBlock>>>(d_a, d_out, kGridDim,
kStencilR);
};
auto launchStencilTiled = [&] {
heatStepTiled<<<stencilGrid, stencilBlock>>>(d_a, d_out, kGridDim,
kStencilR);
};
auto launchPolynomial = [&] {
polynomialHorner<<<vecGrid, kThreadsPerBlock>>>(d_a, d_out, kElems);
};
auto launchFmaPeak = [&] {
fmaPeak<<<fmaBlocks, kThreadsPerBlock>>>(d_out, kFmaIters);
};
// Small integers held as floats, `a` in [1, 7] and `b` in [-2, 2], so
// every sum, every partial and every dot product below is exact and a
// mismatch can only be an index bug. The periods are 7 and 5 rather than
// 4 and 8 because kMatDim is 512: a period that divides 512 would make
// every row of `a` identical, and a matmul that read the wrong row would
// then pass.
for (size_t i = 0; i < kElems; ++i) {
h_a[i] = static_cast<float>(i % 7) + 1.0f;
h_b[i] = static_cast<float>(i % 5) - 2.0f;
}
if (profile) {
// One launch per kernel, in table order, and nothing else. An ncu
// pass over the normal path would replay 168 launches and produce a
// report nobody can read; this produces twelve.
//
// No check runs here and no number is printed, so this mode proves
// nothing on its own. It exists only to give the profiler the same
// twelve launch configurations the timed path measures.
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), bytes, cudaMemcpyHostToDevice));
launchCopy();
launchVectorAdd();
launchVectorAddPermuted();
launchTransposeNaive();
launchTransposeTiled();
launchMatmulNaive();
launchMatmulTiled();
launchReduce();
launchStencilGlobal();
launchStencilTiled();
launchPolynomial();
launchFmaPeak();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_out));
CUDA_CHECK(cudaFree(d_partials));
std::printf("profile mode: 12 launches, no checks, no timings\n");
return EXIT_SUCCESS;
}
// Nothing below returns early. A failed check increments `wrong` and the
// run carries on, so one bad kernel produces a full report instead of the
// first line of one, and the four allocations above are freed once, on
// the single path out of this function.
int wrong = 0;
int rows = 0;
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), bytes, cudaMemcpyHostToDevice));
// Ceiling 1: the coalesced copy.
launchCopy();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
size_t bad = firstMismatch(h_out.data(), h_a.data(), kElems, 0.0f);
if (bad != kElems) {
std::fprintf(stderr,
"copyCoalesced wrong at %zu: got %.9g, want %.9g\n", bad,
h_out[bad], h_a[bad]);
++wrong;
}
const float copyMs = timeKernel(launchCopy);
// Row 1: vector add.
vectorAddCpu(h_a.data(), h_b.data(), h_want.data(), kElems);
launchVectorAdd();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
bad = firstMismatch(h_out.data(), h_want.data(), kElems, kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr, "vectorAdd wrong at %zu: got %.9g, want %.9g\n",
bad, h_out[bad], h_want[bad]);
++wrong;
}
const float vectorAddMs = timeKernel(launchVectorAdd);
// Row 2: the same add over a stride-32 permutation, against the same
// reference. Filling the output with all-ones bytes first is what makes
// that check mean something: an all-ones float is a NaN, so an index the
// permutation never visits fails instead of quietly keeping the answer
// the previous kernel left there.
CUDA_CHECK(cudaMemset(d_out, 0xFF, bytes));
launchVectorAddPermuted();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
bad = firstMismatch(h_out.data(), h_want.data(), kElems, kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr,
"vectorAddPermuted wrong at %zu: got %.9g, want %.9g. "
"The 2D grid is not covering every index once\n",
bad, h_out[bad], h_want[bad]);
++wrong;
}
const float vectorAddPermutedMs = timeKernel(launchVectorAddPermuted);
// Rows 3 and 4: the two transposes. The reference is built from the
// index fill, so the input is refilled here and refilled back afterwards
// for the kernels that need small values. Over a buffer that repeats, a
// swapped index can land on an equal element and pass, which is why the
// transpose gets its own input.
for (size_t i = 0; i < kElems; ++i) {
h_out[i] = static_cast<float>(i);
}
CUDA_CHECK(cudaMemcpy(d_a, h_out.data(), bytes, cudaMemcpyHostToDevice));
for (size_t row = 0; row < kSide; ++row) {
for (size_t col = 0; col < kSide; ++col) {
h_want[col * kSide + row] = h_out[row * kSide + col];
}
}
CUDA_CHECK(cudaMemset(d_out, 0xFF, bytes));
launchTransposeNaive();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
bad = firstMismatch(h_out.data(), h_want.data(), kElems, 0.0f);
if (bad != kElems) {
std::fprintf(stderr,
"transposeNaive wrong at %zu (output row %zu, col %zu): "
"got %.9g, want %.9g\n",
bad, bad / kSide, bad % kSide, h_out[bad], h_want[bad]);
++wrong;
}
const float transposeNaiveMs = timeKernel(launchTransposeNaive);
CUDA_CHECK(cudaMemset(d_out, 0xFF, bytes));
launchTransposeTiled();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
bad = firstMismatch(h_out.data(), h_want.data(), kElems, 0.0f);
if (bad != kElems) {
std::fprintf(stderr,
"transposeTiled wrong at %zu (output row %zu, col %zu): "
"got %.9g, want %.9g\n",
bad, bad / kSide, bad % kSide, h_out[bad], h_want[bad]);
++wrong;
}
const float transposeTiledMs = timeKernel(launchTransposeTiled);
// Back to the small values for everything that follows.
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), bytes, cudaMemcpyHostToDevice));
// Rows 5 and 6: the two matmuls, over the first kMatDim squared floats of
// the same buffers. The CPU reference runs once and both kernels are
// compared against it, because the point of the pair is that they compute
// the same product while asking memory for very different amounts.
matmulCpu(h_a.data(), h_b.data(), h_want.data(), kMatDim);
launchMatmulNaive();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, matElems * sizeof(float),
cudaMemcpyDeviceToHost));
bad = firstMismatch(h_out.data(), h_want.data(), matElems, kRelTolerance);
if (bad != matElems) {
std::fprintf(stderr,
"matmulNaive wrong at %zu (row %zu, col %zu): got %.9g, "
"want %.9g\n",
bad, bad / kMatDim, bad % kMatDim, h_out[bad],
h_want[bad]);
++wrong;
}
const float matmulNaiveMs = timeKernel(launchMatmulNaive);
launchMatmulTiled();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, matElems * sizeof(float),
cudaMemcpyDeviceToHost));
bad = firstMismatch(h_out.data(), h_want.data(), matElems, kRelTolerance);
if (bad != matElems) {
std::fprintf(stderr,
"matmulTiled wrong at %zu (row %zu, col %zu): got %.9g, "
"want %.9g\n",
bad, bad / kMatDim, bad % kMatDim, h_out[bad],
h_want[bad]);
++wrong;
}
const float matmulTiledMs = timeKernel(launchMatmulTiled);
// Row 7: the block reduction. The kernel sums inside a block; the host
// sums the one float each block returned, in double, so the comparison is
// against an exact total rather than against a second tree.
std::vector<float> h_partials(vecBlocks);
launchReduce();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_partials.data(), d_partials,
vecBlocks * sizeof(float), cudaMemcpyDeviceToHost));
const double gpuSum = sumCpu(h_partials.data(), vecBlocks);
const double cpuSum = sumCpu(h_a.data(), kElems);
if (std::fabs(gpuSum - cpuSum) >
static_cast<double>(kRelTolerance) * std::fabs(cpuSum)) {
std::fprintf(stderr,
"reduceSequentialAddressing wrong: got %.9g, want %.9g\n",
gpuSum, cpuSum);
++wrong;
}
const float reduceMs = timeKernel(launchReduce);
// Rows 8 and 9: the two stencils, over the first kGridDim squared floats.
// Zeroing the device output before each one lets the check run over the
// whole grid, boundary ring included, because the reference writes zero
// there and neither kernel touches it.
const size_t stencilBytes = kStencilCells * sizeof(float);
heatStepCpu(h_a.data(), h_want.data(), kGridDim, kStencilR);
CUDA_CHECK(cudaMemset(d_out, 0, stencilBytes));
launchStencilGlobal();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_out.data(), d_out, stencilBytes, cudaMemcpyDeviceToHost));
bad = firstMismatch(h_out.data(), h_want.data(), kStencilCells,
kRelTolerance);
if (bad != kStencilCells) {
std::fprintf(stderr,
"heatStepGlobal wrong at %zu (row %zu, col %zu): got "
"%.9g, want %.9g\n",
bad, bad / kGridDim, bad % kGridDim, h_out[bad],
h_want[bad]);
++wrong;
}
const float stencilGlobalMs = timeKernel(launchStencilGlobal);
CUDA_CHECK(cudaMemset(d_out, 0, stencilBytes));
launchStencilTiled();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_out.data(), d_out, stencilBytes, cudaMemcpyDeviceToHost));
bad = firstMismatch(h_out.data(), h_want.data(), kStencilCells,
kRelTolerance);
if (bad != kStencilCells) {
std::fprintf(stderr,
"heatStepTiled wrong at %zu (row %zu, col %zu): got %.9g, "
"want %.9g\n",
bad, bad / kGridDim, bad % kGridDim, h_out[bad],
h_want[bad]);
++wrong;
}
const float stencilTiledMs = timeKernel(launchStencilTiled);
// Row 10: the Horner chain. The host reference is 2.2 billion
// multiply-adds at full size, so it checks a prefix and the finite scan
// covers the rest of the buffer.
CUDA_CHECK(cudaMemset(d_out, 0xFF, bytes));
launchPolynomial();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
polynomialCpu(h_a.data(), h_want.data(), kPolyCheck);
bad = firstMismatch(h_out.data(), h_want.data(), kPolyCheck, kRelTolerance);
if (bad != kPolyCheck) {
std::fprintf(stderr,
"polynomialHorner wrong at %zu: got %.9g, want %.9g. The "
"chain did not run the %d steps the flop count assumes\n",
bad, h_out[bad], h_want[bad], kPolyDegree);
++wrong;
}
bad = firstNonFinite(h_out.data(), kElems);
if (bad != kElems) {
std::fprintf(stderr,
"polynomialHorner left %zu unwritten or not finite: "
"%.9g\n",
bad, h_out[bad]);
++wrong;
}
const float polynomialMs = timeKernel(launchPolynomial);
// Ceiling 2: the FMA chain.
launchFmaPeak();
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, fmaThreads * sizeof(float),
cudaMemcpyDeviceToHost));
const float fmaWant = fmaSumCpu(kFmaIters);
for (size_t i = 0; i < fmaThreads; ++i) {
h_want[i] = fmaWant;
}
bad = firstMismatch(h_out.data(), h_want.data(), fmaThreads, kFmaTolerance);
if (bad != fmaThreads) {
std::fprintf(stderr,
"fmaPeak wrong at %zu: got %.9g, want %.9g. The loop did "
"not run the %d iterations the flop count assumes\n",
bad, h_out[bad], h_want[bad], kFmaIters);
++wrong;
}
const float fmaMs = timeKernel(launchFmaPeak);
// Every elapsed time is about to become a denominator. A zero would mean
// the event pair never separated and every figure in the table would be
// nonsense rather than merely wrong.
const float everyMs[kNumTimings] = {
copyMs, vectorAddMs, vectorAddPermutedMs, transposeNaiveMs,
transposeTiledMs, matmulNaiveMs, matmulTiledMs, reduceMs,
stencilGlobalMs, stencilTiledMs, polynomialMs, fmaMs};
for (int r = 0; r < kNumTimings; ++r) {
if (!(everyMs[r] > 0.0f)) {
std::fprintf(stderr,
"timing %d came back as %.9g ms, which cannot be "
"divided into\n",
r, everyMs[r]);
++wrong;
}
}
// The device is finished with. Freeing here rather than at each exit puts
// every cudaMalloc's matching cudaFree on one path, including the path
// taken when an answer is wrong. Nothing below touches the GPU.
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_out));
CUDA_CHECK(cudaFree(d_partials));
if (wrong != 0) {
std::fprintf(stderr,
"%d check(s) failed; no roofline is printed from a run "
"whose kernels are wrong\n",
wrong);
return EXIT_FAILURE;
}
// The two ceilings. Both are measured. Neither is a vendor figure, and
// the widget on the page refuses one for the same reason: a ceiling you
// cannot reach tells you nothing about how close you got.
// snippet: ridge
const double copyGBs = 2.0 * static_cast<double>(kElems) * sizeof(float) /
(static_cast<double>(copyMs) * 1.0e-3) / 1.0e9;
const double fmaFlops = static_cast<double>(fmaThreads) *
static_cast<double>(kFmaIters) * kFmaChains * 2.0;
const double peakGFlops =
fmaFlops / (static_cast<double>(fmaMs) * 1.0e-3) / 1.0e9;
const double ridge = peakGFlops / copyGBs;
// end snippet: ridge
// FLOP and compulsory byte counts, worked from the kernels above and
// nowhere else. Compulsory means the traffic the algorithm asks global
// memory for, which is the textbook convention and is not the traffic
// DRAM sees. Nsight Compute's dram__bytes.sum is the other column, and
// merge_dram.py is what puts the two side by side.
//
// The stencil count is 7 flops per interior cell, from the expression as
// written: three adds over the four neighbours, one multiply and one
// subtract for the centre term, one multiply by r and one final add. A
// compiler that fuses two of those into an FMA does not change the count
// this table uses, which is the model's own warning: you supply the
// numbers and no tool checks them.
const double n = static_cast<double>(kElems);
const double blocks = static_cast<double>(vecBlocks);
const double d3 = static_cast<double>(kMatDim) *
static_cast<double>(kMatDim) *
static_cast<double>(kMatDim);
const double d2 = static_cast<double>(matElems);
const double interior = static_cast<double>(kStencilInterior);
const double apron = static_cast<double>(kSharedCells) *
static_cast<double>(stencilBlocksX) *
static_cast<double>(stencilBlocksY);
const char* names[kNumRows] = {"vector add", "vector add, stride 32",
"transpose, naive", "transpose, tiled",
"matmul, naive", "matmul, tiled 16",
"block reduction", "stencil, global",
"stencil, tiled", "polynomial, 256"};
const char* symbols[kNumRows] = {"vectorAdd",
"vectorAddPermuted",
"transposeNaive",
"transposeTiled",
"matmulNaive",
"matmulTiled",
"reduceSequentialAddressing",
"heatStepGlobal",
"heatStepTiled",
"polynomialHorner"};
const int days[kNumRows] = {5, 11, 12, 13, 16, 16, 24, 34, 34, 49};
const double flops[kNumRows] = {n,
n,
0.0,
0.0,
2.0 * d3,
2.0 * d3,
n - blocks,
7.0 * interior,
7.0 * interior,
(2.0 * kPolyDegree + 1.0) * n};
// snippet: bytes-asked
const double moved[kNumRows] = {3.0 * n * sizeof(float),
3.0 * n * sizeof(float),
2.0 * n * sizeof(float),
2.0 * n * sizeof(float),
(2.0 * d3 + d2) * sizeof(float),
(2.0 * d3 / kTileDim + d2) * sizeof(float),
(n + blocks) * sizeof(float),
6.0 * interior * sizeof(float),
(apron + interior) * sizeof(float),
2.0 * n * sizeof(float)};
// end snippet: bytes-asked
const double milliseconds[kNumRows] = {
vectorAddMs, vectorAddPermutedMs, transposeNaiveMs, transposeTiledMs,
matmulNaiveMs, matmulTiledMs, reduceMs, stencilGlobalMs,
stencilTiledMs, polynomialMs};
if (json) {
std::printf(
"{\"variant\":\"ceiling-copy\",\"symbol\":\"copyCoalesced\","
"\"ms\":%.6f,\"metric\":\"gbps\",\"value\":%.3f,"
"\"unit\":\"GB/s\"}\n",
copyMs, copyGBs);
std::printf(
"{\"variant\":\"ceiling-fma\",\"symbol\":\"fmaPeak\","
"\"ms\":%.6f,\"metric\":\"gflops\",\"value\":%.3f,"
"\"unit\":\"GFLOP/s\"}\n",
fmaMs, peakGFlops);
std::printf(
"{\"variant\":\"ridge-point\",\"metric\":\"ridge\","
"\"value\":%.4f,\"unit\":\"FLOP/byte\"}\n",
ridge);
} else {
std::printf("GPU: %s (compute capability %d.%d), %d SMs, %d KiB L2\n",
prop.name, prop.major, prop.minor, prop.multiProcessorCount,
prop.l2CacheSize / 1024);
std::printf(
"%zu floats per buffer (%.0f MiB), matmul %zu cubed, tile %d, "
"stencil %d squared\n\n",
kElems, static_cast<double>(bytes) / (1024.0 * 1024.0), kMatDim,
kTileDim, kGridDim);
std::printf("The two ends of this card's roofline, both measured\n");
std::printf(" %-32s %9.3f ms %9.1f GB/s\n",
"copy, coalesced, in and out", copyMs, copyGBs);
std::printf(" %-32s %9.3f ms %9.1f GFLOP/s\n", "fma, 8 chains", fmaMs,
peakGFlops);
std::printf(" %-32s %9s %9.2f FLOP/byte\n", "ridge point", "",
ridge);
std::printf(
" Below that intensity nothing on this card can be compute\n"
" bound, whatever its arithmetic looks like.\n\n");
std::printf("Ten kernels, with the bytes the algorithm asks for\n");
std::printf("%-22s %3s %9s %8s %9s %9s %6s\n", "kernel", "day",
"FLOP/byte", "ms", "ask MB", "ask GB/s", "%ceil");
std::printf("%-22s %3s %9s %8s %9s %9s %6s\n", "----------------------",
"---", "---------", "--------", "---------", "---------",
"------");
}
int memoryBound = 0;
for (int r = 0; r < kNumRows; ++r) {
const double seconds = static_cast<double>(milliseconds[r]) * 1.0e-3;
const double intensity = flops[r] / moved[r];
const double askGBs = moved[r] / seconds / 1.0e9;
const double gflops = flops[r] / seconds / 1.0e9;
const bool memory = intensity < ridge;
const double percent =
memory ? 100.0 * askGBs / copyGBs : 100.0 * gflops / peakGFlops;
if (memory) {
++memoryBound;
}
if (json) {
std::printf(
"{\"variant\":\"%s\",\"symbol\":\"%s\",\"day\":%d,"
"\"ms\":%.6f,\"flop\":%.0f,\"bytes\":%.0f,"
"\"metric\":\"ms\",\"value\":%.6f,\"unit\":\"ms\"}\n",
names[r], symbols[r], days[r], milliseconds[r], flops[r],
moved[r], milliseconds[r]);
} else {
std::printf("%-22s %3d %9.3f %8.3f %9.1f %9.1f %6.1f\n", names[r],
days[r], intensity, milliseconds[r], moved[r] / 1.0e6,
askGBs, percent);
}
++rows;
}
if (!json) {
std::printf(
"\n%d of the %d sit left of the ridge point, so for those the\n"
"sloped ceiling is the one that binds.\n",
memoryBound, kNumRows);
std::printf(
"\nThe two right-hand columns are COMPULSORY bytes over measured\n"
"time. That is what the algorithm asked for, not what DRAM moved,\n"
"and any row over 100 percent is the tell. Run the Nsight\n"
"Compute pass in README.md, then merge_dram.py, for the\n"
"measured column.\n");
}
// The row count is checked so the lesson's table and this program cannot
// drift apart on how many rows there are. A real branch, not an assert:
// CI builds Release, Release defines NDEBUG, and NDEBUG deletes assert(),
// so the check would be missing from exactly the build that matters.
if (rows != kNumRows) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's table and "
"this program disagree\n",
rows, kNumRows);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}