code/day30-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 30, checkpoint: the memory-bound mindset.
//
// Measures the two ends of this card's roofline, then places five kernels
// from days 5, 12, 16 and 24 on it.
//
// One program rather than five, because the two ceilings and the five points
// have to come off the same card in the same clock state. Borrowing day 11's
// copy bandwidth and day 16's matmul time would fold two buffer sizes, two
// block shapes and two thermal states into one chart, and the chart would
// still look tidy.
//
// No vendor peak appears anywhere in this file. A roofline drawn against a
// marketing number is a roofline nobody can reach, so both ceilings come from
// a kernel that ran here: a coalesced copy for bandwidth, and a chain of
// fused multiply-adds with no loads in the loop for FP32.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o roofline roofline.cu
// Run: ./roofline
//
// Verified: run on a Tesla T4 (compute capability 7.5), driver 595.84,
// CUDA 12.6 (V12.6.85), on 2026-08-30. The only output that may be
// published as this program's output is the transcript in
// evidence/run-2026-08-30.txt. See research/REVIEW-PROCESS.md.
#include <cmath>
#include <cstdio>
#include <cstdlib>
#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 the other three
// vector kernels then share 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, which is the last integer a float represents
// exactly. The transpose check below fills the input with its own index and
// compares exactly, which only works 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 matmul is square and divides the tile exactly, so no kernel here needs
// an edge guard. Day 16 owns the ragged case, where 611 by 613 by 617 divides
// by 16 on no axis and a missing guard is silently wrong.
constexpr size_t kMatDim = 512;
constexpr int kTileDim = 16;
// The FP32 ceiling. 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 of these are ever resident
// together; day 28 measures exactly that, 4 blocks per SM. Launching 32 per SM
// queues eight waves behind those 4, which is what the number buys: the tail,
// the last partial wave, is 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 = 5;
constexpr int kNumTimings = 7; // five 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, so a run-time
// check of them would look like a measurement and be nothing of the kind.
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(kSide % kTransposeDim == 0,
"the transpose grid tiles the matrix exactly, which is why its "
"kernel carries no bounds check");
static_assert(kTransposeDim * kTransposeDim <= 1024,
"a transpose block is at most 1024 threads");
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(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, and a reader
// copying this file should copy the habit.
__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, 12 bytes moved 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];
}
}
// 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 is
// what puts this kernel at the far left of the chart.
//
// Launch assumption: a kTransposeDim square block over a grid that tiles the
// matrix exactly, which the static_assert above guarantees, so no bounds
// check is needed or present.
__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];
}
// Sums one block's slice of `in` into out[blockIdx.x]. This is 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 highest intensity any streaming kernel in modules 1 to 3
// reaches, and it is still far to the left of the ridge point below.
//
// 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 both are decent
// address patterns. The problem is the count, not the pattern: the same row
// of `a` is fetched again by every thread in the block.
//
// 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 global load count falls by kTileDim and the flop count does
// not move at all. That is the only optimisation in this file that changes a
// kernel's arithmetic intensity rather than how close it gets to 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;
}
// The FP32 ceiling: fused multiply-adds and nothing else.
//
// Memory: one coalesced store per thread, after the loop. The loop reads and
// writes registers only, so nothing in it is limited by bandwidth. Eight
// accumulators per thread is well inside the T4's register budget; day 17 is
// where a per-thread array stops fitting and spills to local memory.
//
// 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.
// snippet: fma-ceiling
__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;
}
// end snippet: fma-ceiling
// 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 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;
}
// 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. 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;
}
int main() {
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 int fmaBlocks = kFmaBlocksLaunchedPerSm * prop.multiProcessorCount;
const size_t fmaThreads = static_cast<size_t>(fmaBlocks) * kThreadsPerBlock;
std::printf("GPU: %s (compute capability %d.%d), %d SMs\n", prop.name,
prop.major, prop.minor, prop.multiProcessorCount);
std::printf(
"%zu floats per buffer (%.0f MiB), matmul %zu cubed, tile %d\n\n",
kElems, static_cast<double>(bytes) / (1024.0 * 1024.0), kMatDim,
kTileDim);
// 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)));
// 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;
// Phase 1 fills the input with its own index. The transpose check needs
// distinct values: over a buffer that repeats, a swapped index can land
// on an equal element and pass. Phase 2 refills with small values that
// keep the matmul exact.
for (size_t i = 0; i < kElems; ++i) {
h_a[i] = static_cast<float>(i);
}
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), bytes, cudaMemcpyHostToDevice));
// Ceiling 1: the coalesced copy.
copyCoalesced<<<vecGrid, kThreadsPerBlock>>>(d_a, d_out, kElems);
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([&] {
copyCoalesced<<<vecGrid, kThreadsPerBlock>>>(d_a, d_out, kElems);
});
// Kernel A: the naive transpose.
const dim3 transposeBlock(kTransposeDim, kTransposeDim);
const dim3 transposeGrid(static_cast<unsigned int>(kSide / kTransposeDim),
static_cast<unsigned int>(kSide / kTransposeDim));
transposeNaive<<<transposeGrid, transposeBlock>>>(d_a, d_out, kSide);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
for (size_t row = 0; row < kSide; ++row) {
for (size_t col = 0; col < kSide; ++col) {
h_want[col * kSide + row] = h_a[row * kSide + col];
}
}
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 transposeMs = timeKernel([&] {
transposeNaive<<<transposeGrid, transposeBlock>>>(d_a, d_out, kSide);
});
// Phase 2. 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;
}
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), bytes, cudaMemcpyHostToDevice));
// Kernel B: vector add.
vectorAdd<<<vecGrid, kThreadsPerBlock>>>(d_a, d_b, d_out, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
vectorAddCpu(h_a.data(), h_b.data(), h_want.data(), kElems);
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([&] {
vectorAdd<<<vecGrid, kThreadsPerBlock>>>(d_a, d_b, d_out, kElems);
});
// Kernel C: the block reduction. The kernel sums inside a block; the host
// sums the one float each block returned, in double, so the comparison
// below is against an exact total rather than against a second tree.
std::vector<float> h_partials(vecBlocks);
reduceSequentialAddressing<<<vecGrid, kThreadsPerBlock>>>(d_a, d_partials,
kElems);
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([&] {
reduceSequentialAddressing<<<vecGrid, kThreadsPerBlock>>>(
d_a, d_partials, kElems);
});
// Kernels D and E: 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 at different intensities.
matmulCpu(h_a.data(), h_b.data(), h_want.data(), kMatDim);
const dim3 matBlock(kTileDim, kTileDim);
const dim3 matGrid(static_cast<unsigned int>(kMatDim / kTileDim),
static_cast<unsigned int>(kMatDim / kTileDim));
matmulNaive<<<matGrid, matBlock>>>(d_a, d_b, d_out, kMatDim);
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(
[&] { matmulNaive<<<matGrid, matBlock>>>(d_a, d_b, d_out, kMatDim); });
matmulTiled<<<matGrid, matBlock>>>(d_a, d_b, d_out, kMatDim);
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(
[&] { matmulTiled<<<matGrid, matBlock>>>(d_a, d_b, d_out, kMatDim); });
// Ceiling 2: the FMA chain.
fmaPeak<<<fmaBlocks, kThreadsPerBlock>>>(d_out, kFmaIters);
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(
[&] { fmaPeak<<<fmaBlocks, kThreadsPerBlock>>>(d_out, kFmaIters); });
// 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, transposeMs, vectorAddMs,
reduceMs, matmulNaiveMs, matmulTiledMs,
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;
}
// snippet: ridge
// 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.
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
char fmaLabel[64];
std::snprintf(fmaLabel, sizeof(fmaLabel), "fma, %d chains, %d iterations",
kFmaChains, kFmaIters);
std::printf("The two ends of this card's roofline, both measured\n");
std::printf(" %-34s %9.3f ms %9.1f GB/s\n", "copy, coalesced, in and out",
copyMs, copyGBs);
std::printf(" %-34s %9.3f ms %9.1f GFLOP/s\n", fmaLabel, fmaMs,
peakGFlops);
std::printf(" %-34s %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");
// FLOP and byte counts, worked from the kernels above and nowhere else.
// Bytes is the traffic the algorithm asks global memory for, which is the
// textbook convention and is not the same as the traffic DRAM sees: the
// caches serve some of it, which is why a row below can read above 100
// percent of a ceiling. The reduction's add count is exact rather than n,
// because a block of 256 does 255 adds.
// snippet: counts
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 char* names[kNumRows] = {"vector add", "transpose, naive",
"block reduction", "matmul, naive",
"matmul, tiled 16"};
const double flops[kNumRows] = {n, 0.0, n - blocks, 2.0 * d3, 2.0 * d3};
const double moved[kNumRows] = {
3.0 * n * sizeof(float), 2.0 * n * sizeof(float),
(n + blocks) * sizeof(float), (2.0 * d3 + d2) * sizeof(float),
(2.0 * d3 / kTileDim + d2) * sizeof(float)};
const double milliseconds[kNumRows] = {vectorAddMs, transposeMs, reduceMs,
matmulNaiveMs, matmulTiledMs};
// end snippet: counts
std::printf("Five kernels from days 5, 12, 16 and 24, placed on it\n");
std::printf("%-20s %9s %7s %7s %8s %9s %5s\n", "kernel", "FLOP/byte",
"ms", "GB/s", "GFLOP/s", "bound by", "%ceil");
std::printf("%-20s %9s %7s %7s %8s %9s %5s\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 gbs = moved[r] / seconds / 1.0e9;
const double gflops = flops[r] / seconds / 1.0e9;
const bool memory = intensity < ridge;
const double percent =
memory ? 100.0 * gbs / copyGBs : 100.0 * gflops / peakGFlops;
if (memory) {
++memoryBound;
}
std::printf("%-20s %9.3f %7.3f %7.1f %8.1f %9s %5.1f\n", names[r],
intensity, milliseconds[r], gbs, gflops,
memory ? "bandwidth" : "compute", percent);
++rows;
}
std::printf(
"\n%d of the %d sit left of the ridge point, so for those the sloped\n"
"ceiling is the one that binds and a faster inner loop buys nothing.\n",
memoryBound, kNumRows);
std::printf(
"FLOP/byte counts the traffic the algorithm asks for. The caches\n"
"serve part of it, so a row can read above 100 percent of a ceiling\n"
"the model says it should be under.\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;
}