code/day42-ncu/matmul_profile.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 42: the program the shipped Nsight Compute reports come from.
//
// The two kernels are day 16's, unchanged, at one square size. What this file
// adds is a fixed launch order, so a capture command can name the launch it
// wants without anyone counting by hand, and a printed prediction of what the
// report's memory rows should say, computed from the constants below.
//
// 512 x 512 x 512 divides by the 16-wide tile on every axis, so no bounds
// guard ever fires and the report describes the steady-state shape. Day 16
// owns the ragged case; a report of an edge tile would teach the wrong thing
// here. The size also matches day 30's roofline row, so the dram__bytes.sum
// the report gives can be put next to the byte count that day computed by
// hand and reported at 1008 percent of the card's copy ceiling.
//
// Each kernel launches 14 times: 1 correctness launch, 3 warm-ups, 10 timed.
// The reports profile the first timed launch, which is why the capture
// command passes --launch-skip 4. The program prints that command with the
// number filled in from its own constants, so README.md and this file cannot
// drift apart on which launch was captured.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo -o matmul_profile \
// matmul_profile.cu
// Run: ./matmul_profile
//
#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)
// A 16 x 16 tile is one 256-thread block, the course default, and two tiles
// of floats is 2 KiB of shared memory, nowhere near the 48 KiB a T4 gives a
// block. Day 45 is where shared memory starts to bind.
constexpr int kTileDim = 16;
constexpr int kWarpSize = 32; // 32 on every GPU this course targets
constexpr size_t kDim = 512;
// The launch order the capture command depends on.
constexpr int kCorrectnessRuns = 1;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr int kLaunchSkip = kCorrectnessRuns + kWarmupRuns;
constexpr float kRelTolerance = 1e-5f;
// What the report should say, derived from the constants above and printed
// before any counter is read. A block is 16 x 16 with threadIdx.x fastest, so
// one warp of 32 lanes is two rows of sixteen consecutive columns.
//
// matmulNaive issues 2 * K global load instructions per thread, and a warp
// issues one instruction for all 32 of its lanes, so 2 * K requests per warp.
// Across the warp, a[row * k + p] has two distinct addresses, one per row,
// K floats apart, so it touches two 32-byte sectors. b[p * n + col] is
// sixteen consecutive floats, 64 bytes, aligned because col starts on a
// multiple of the tile width, so it touches two sectors as well.
//
// matmulTiled issues 2 per tile step, so 2 * K / 16 per warp. Each of those
// loads sixteen consecutive floats on each of the warp's two rows, which is
// two runs of 64 bytes and therefore four sectors.
//
// The naive kernel wins on sectors per request and loses by eight on sectors
// moved. Both facts belong on the page: a good ratio over a huge count is
// still a huge count.
//
// The last two constants are the point of the whole day. Tiling does not cut
// the number of load instructions a thread issues; it moves them. Two loads
// per k step is two loads per k step either way. What changes is which memory
// each load asks, and a shared-memory load fetches no sector from L2.
// snippet: predicted-counts
constexpr size_t kSectorBytes = 32;
constexpr size_t kBlocks = (kDim / kTileDim) * (kDim / kTileDim);
constexpr size_t kWarps = kBlocks * ((kTileDim * kTileDim) / kWarpSize);
constexpr size_t kNaiveRequests = kWarps * 2 * kDim;
constexpr size_t kTiledRequests = kWarps * 2 * (kDim / kTileDim);
constexpr size_t kNaiveSectorsPerRequest = 2;
constexpr size_t kTiledSectorsPerRequest = 4;
constexpr size_t kNaiveSectors = kNaiveRequests * kNaiveSectorsPerRequest;
constexpr size_t kTiledSectors = kTiledRequests * kTiledSectorsPerRequest;
constexpr size_t kNaiveGlobalLoads = 2 * kDim;
constexpr size_t kTiledGlobalLoads = 2 * (kDim / kTileDim);
constexpr size_t kTiledSharedLoads = 2 * kTileDim * (kDim / kTileDim);
// A, B and C read or written once each. This is the floor on DRAM traffic,
// not a prediction of it: both kernels re-read A and B many times and the
// caches decide how much of that reaches memory.
constexpr size_t kCompulsoryBytes = 3 * kDim * kDim * sizeof(float);
// end snippet
static_assert(kTileDim * kTileDim == 256,
"a 16 x 16 tile is one 256-thread block, the course default");
static_assert((kTileDim * kTileDim) % kWarpSize == 0,
"block size must be a whole number of warps");
static_assert(2 * kTileDim * kTileDim * sizeof(float) <= 48 * 1024,
"the two tiles must fit the 48 KiB a block gets on a T4 "
"without the cudaFuncSetAttribute opt-in");
static_assert(kDim % kTileDim == 0,
"the profiled size divides by the tile on every axis, so no "
"guard fires and the report shows the steady-state shape; day "
"16 owns the ragged case");
static_assert(kTileDim * sizeof(float) == 2 * kSectorBytes,
"the sector counts above assume a warp's sixteen columns span "
"64 bytes, which is two 32-byte sectors");
static_assert(kTiledSharedLoads == kNaiveGlobalLoads,
"the tiled kernel's shared loads are counted as two per inner "
"iteration over the whole k range and the naive kernel's "
"global loads as two per k step; they are the same number, "
"which is the claim the lesson's memory chart rests on");
static_assert(kDim / kTileDim <= 65535,
"gridDim.y stops at 65535, unlike gridDim.x which goes to "
"2^31 - 1");
// Integer ceiling division. constexpr so one function sizes a grid at run
// time and can appear in a static_assert. Anything that follows only from the
// constants above has to be a static_assert rather than an assert: CI builds
// Release, Release defines NDEBUG, and NDEBUG deletes assert().
static constexpr size_t ceilDiv(size_t a, size_t b) {
return (a + b - 1) / b;
}
// C[row][col] = sum over p of A[row][p] * B[p][col]. One thread owns one
// output element and reads 2K floats from global memory to produce it.
//
// Memory: threadIdx.x is the fastest-varying index, so a warp of this block
// is two rows of sixteen consecutive columns. b[p * n + col] is therefore two
// runs of sixteen consecutive floats and a[row * k + p] is two addresses the
// halves of the warp share. Nothing here is laid out badly. The cost is how
// many times the same value comes back across the bus, not the pattern it
// comes back in, and that is the distinction this day's report makes visible.
//
// Launch assumption: the grid covers C exactly, because kDim divides the tile
// on every axis. The guard is kept anyway, so this kernel stays the one day
// 16 published rather than a trimmed copy of it.
__global__ void matmulNaive(const float* __restrict__ a,
const float* __restrict__ b, float* __restrict__ c,
size_t m, size_t n, size_t k) {
const size_t col =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t row =
blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
if (row < m && col < n) {
float acc = 0.0f;
for (size_t p = 0; p < k; ++p) {
acc += a[row * k + p] * b[p * n + col];
}
c[row * n + col] = acc;
}
}
// The same product, one 16 x 16 tile of C per block, staged through shared
// memory. Global loads per thread fall from 2K to 2 * ceil(K / 16).
//
// Memory: the global loads have the same geometry as the naive kernel's, so
// the two reports differ in how many loads issue and in nothing else about
// the address pattern. In shared memory a warp reads tileA[threadIdx.y][p],
// two addresses its halves share, and tileB[p][threadIdx.x], sixteen
// consecutive words read by two lanes each. Both broadcast or hit distinct
// banks, so neither conflicts; day 15 is where that stops being free.
//
// Launch assumption: exactly kTileDim x kTileDim threads per block. Every
// thread reaches both barriers. The guard covers the two loads and the store
// and never the __syncthreads(), because a barrier only part of a block
// arrives at is undefined behaviour. Day 14 has the rule.
__global__ void matmulTiled(const float* __restrict__ a,
const float* __restrict__ b, float* __restrict__ c,
size_t m, size_t n, size_t k) {
__shared__ float tileA[kTileDim][kTileDim];
__shared__ float tileB[kTileDim][kTileDim];
const size_t col =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t row =
blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
// ceilDiv is a host function, so the same arithmetic is written out here.
const size_t tiles = (k + kTileDim - 1) / kTileDim;
float acc = 0.0f;
for (size_t tileIdx = 0; tileIdx < tiles; ++tileIdx) {
const size_t aCol = tileIdx * kTileDim + threadIdx.x;
const size_t bRow = tileIdx * kTileDim + threadIdx.y;
tileA[threadIdx.y][threadIdx.x] =
(row < m && aCol < k) ? a[row * k + aCol] : 0.0f;
tileB[threadIdx.y][threadIdx.x] =
(bRow < k && col < n) ? b[bRow * n + col] : 0.0f;
__syncthreads();
for (int p = 0; p < kTileDim; ++p) {
acc += tileA[threadIdx.y][p] * tileB[p][threadIdx.x];
}
__syncthreads();
}
if (row < m && col < n) {
c[row * n + col] = acc;
}
}
// CPU reference. Three nested loops in the obvious order, written for obvious
// correctness rather than speed. It never allocates; the caller owns every
// buffer. It accumulates in double where the kernels accumulate in float,
// because a reference exists to be right rather than to match bit for bit.
static void matmulCpu(const float* a, const float* b, float* c, size_t m,
size_t n, size_t k) {
for (size_t row = 0; row < m; ++row) {
for (size_t col = 0; col < n; ++col) {
double total = 0.0;
for (size_t p = 0; p < k; ++p) {
total += static_cast<double>(a[row * k + p]) *
static_cast<double>(b[p * n + col]);
}
c[row * n + col] = static_cast<float>(total);
}
}
}
// Fills both inputs with small integers held as floats: A in [-3, 3] and B in
// [-2, 2]. Every product is at most 6, so the largest dot product here is
// 6 * 512 = 3,072, far below the 2^24 above which a float stops representing
// consecutive integers. Every intermediate on both processors is exact, so a
// mismatch in main() is an indexing bug and can be nothing else.
//
// The two periods, 7 and 5, divide neither dimension, so a kernel that
// transposes an index reads a different value rather than the same one back.
static void makeInputs(std::vector<float>* h_a, std::vector<float>* h_b,
size_t elems) {
for (size_t e = 0; e < elems; ++e) {
(*h_a)[e] = static_cast<float>(static_cast<int>(e % 7) - 3);
(*h_b)[e] = static_cast<float>(static_cast<int>(e % 5) - 2);
}
}
// 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: main() turns it back into a row and a column.
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 a launch with CUDA events and returns the mean milliseconds per run.
//
// Copy it verbatim; the alternative is a second copy of the event
// boilerplate, which is how a warm-up goes missing from one of them.
//
// Under ncu these numbers are meaningless, because the profiler serializes
// and replays the kernel it is collecting. Take the timings from a plain run
// and the counters from a profiled one, never both from the same run.
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;
}
// 2 * m * n * k flops for a matrix multiply: one multiply and one add per
// term of every dot product. Both kernels do exactly this much arithmetic.
static double gflops(float ms) {
const double work = 2.0 * static_cast<double>(kDim) *
static_cast<double>(kDim) * static_cast<double>(kDim);
return work / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// Theoretical occupancy from the runtime's own calculator, as a percentage.
// The report's Occupancy section prints the same figure next to the achieved
// one, and the gap between the two is what day 42 asks the reader to explain.
static double occupancyPct(int blocksPerSm, int threadsPerBlock,
int maxThreadsPerSm) {
return 100.0 * static_cast<double>(blocksPerSm) *
static_cast<double>(threadsPerBlock) /
static_cast<double>(maxThreadsPerSm);
}
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), %d SMs\n", prop.name,
prop.major, prop.minor, prop.multiProcessorCount);
std::printf(
"L2 cache %d bytes, shared memory per block %zu bytes, "
"%d threads resident per SM\n",
prop.l2CacheSize, prop.sharedMemPerBlock,
prop.maxThreadsPerMultiProcessor);
std::printf("Shape: %zu x %zu x %zu, tile %d, %zu blocks, %zu warps\n\n",
kDim, kDim, kDim, kTileDim, kBlocks, kWarps);
std::printf("Predicted from the constants, before any counter is read\n");
std::printf("%-34s %14s %14s\n", "quantity", "naive", "tiled");
std::printf("%-34s %14s %14s\n", "----------------------------------",
"--------------", "--------------");
std::printf("%-34s %14zu %14zu\n", "global loads per thread",
kNaiveGlobalLoads, kTiledGlobalLoads);
std::printf("%-34s %14zu %14zu\n", "shared loads per thread",
static_cast<size_t>(0), kTiledSharedLoads);
std::printf("%-34s %14zu %14zu\n", "load instructions per thread",
kNaiveGlobalLoads, kTiledGlobalLoads + kTiledSharedLoads);
std::printf("%-34s %14zu %14zu\n", "global load requests", kNaiveRequests,
kTiledRequests);
std::printf("%-34s %14zu %14zu\n", "sectors per request",
kNaiveSectorsPerRequest, kTiledSectorsPerRequest);
std::printf("%-34s %14zu %14zu\n", "global load sectors", kNaiveSectors,
kTiledSectors);
std::printf("%-34s %14zu %14zu\n", "bytes fetched into L1TEX",
kNaiveSectors * kSectorBytes, kTiledSectors * kSectorBytes);
std::printf("%-34s %14zu %14zu\n\n", "compulsory DRAM bytes",
kCompulsoryBytes, kCompulsoryBytes);
const size_t elems = kDim * kDim;
const size_t bytes = elems * sizeof(float);
std::vector<float> h_a(elems);
std::vector<float> h_b(elems);
std::vector<float> h_c(elems);
std::vector<float> h_want(elems);
makeInputs(&h_a, &h_b, elems);
matmulCpu(h_a.data(), h_b.data(), h_want.data(), kDim, kDim, kDim);
float* d_a = nullptr;
float* d_b = nullptr;
float* d_c = nullptr;
CUDA_CHECK(cudaMalloc(&d_a, bytes));
CUDA_CHECK(cudaMalloc(&d_b, bytes));
CUDA_CHECK(cudaMalloc(&d_c, bytes));
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), bytes, cudaMemcpyHostToDevice));
// x walks the columns of C and y walks its rows, so threadIdx.x stays the
// fastest-varying index over consecutive columns.
const size_t tileDim = static_cast<size_t>(kTileDim);
const dim3 block(kTileDim, kTileDim);
const dim3 grid(static_cast<unsigned int>(ceilDiv(kDim, tileDim)),
static_cast<unsigned int>(ceilDiv(kDim, tileDim)));
// Failure is recorded rather than returned, so every path falls through
// to the one cleanup block below and no cudaMalloc escapes its cudaFree.
const char* badKernel = nullptr;
size_t badIndex = 0;
// Launch 1 of matmulNaive. The output is cleared first so a kernel that
// skips elements is caught by the comparison rather than by leftovers.
CUDA_CHECK(cudaMemset(d_c, 0, bytes));
matmulNaive<<<grid, block>>>(d_a, d_b, d_c, kDim, kDim, kDim);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_c.data(), d_c, bytes, cudaMemcpyDeviceToHost));
size_t bad = firstMismatch(h_c.data(), h_want.data(), elems, kRelTolerance);
if (bad != elems) {
badKernel = "naive";
badIndex = bad;
}
// Launches 2 to 14 of matmulNaive: 3 warm-ups then 10 timed.
float naiveMs = 0.0f;
if (badKernel == nullptr) {
naiveMs = timeKernel([&] {
matmulNaive<<<grid, block>>>(d_a, d_b, d_c, kDim, kDim, kDim);
});
}
float tiledMs = 0.0f;
if (badKernel == nullptr) {
CUDA_CHECK(cudaMemset(d_c, 0, bytes));
matmulTiled<<<grid, block>>>(d_a, d_b, d_c, kDim, kDim, kDim);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_c.data(), d_c, bytes, cudaMemcpyDeviceToHost));
bad = firstMismatch(h_c.data(), h_want.data(), elems, kRelTolerance);
if (bad != elems) {
badKernel = "tiled";
badIndex = bad;
} else {
tiledMs = timeKernel([&] {
matmulTiled<<<grid, block>>>(d_a, d_b, d_c, kDim, kDim, kDim);
});
}
}
int naiveBlocksPerSm = 0;
int tiledBlocksPerSm = 0;
const int threadsPerBlock = kTileDim * kTileDim;
CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
&naiveBlocksPerSm, matmulNaive, threadsPerBlock, 0));
CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
&tiledBlocksPerSm, matmulTiled, threadsPerBlock, 0));
if (badKernel == nullptr) {
std::printf(
"Measured with events, mean of %d runs after %d "
"warm-ups, copies not included\n",
kTimedRuns, kWarmupRuns);
std::printf("%-8s %11s %11s %10s %22s\n", "kernel", "time (ms)",
"GFLOP/s", "blocks/SM", "theoretical occ. (%)");
std::printf("%-8s %11s %11s %10s %22s\n", "--------", "-----------",
"-----------", "----------", "----------------------");
std::printf("%-8s %11.3f %11.1f %10d %22.1f\n", "naive", naiveMs,
gflops(naiveMs), naiveBlocksPerSm,
occupancyPct(naiveBlocksPerSm, threadsPerBlock,
prop.maxThreadsPerMultiProcessor));
std::printf("%-8s %11.3f %11.1f %10d %22.1f\n", "tiled", tiledMs,
gflops(tiledMs), tiledBlocksPerSm,
occupancyPct(tiledBlocksPerSm, threadsPerBlock,
prop.maxThreadsPerMultiProcessor));
std::printf(
"\nEach kernel launched %d times: %d correctness, %d "
"warm-up, %d timed.\nThe shipped reports profile the "
"first timed launch of each, so the capture is:\n\n",
kCorrectnessRuns + kWarmupRuns + kTimedRuns, kCorrectnessRuns,
kWarmupRuns, kTimedRuns);
std::printf(
" sudo ncu --set full --kernel-name regex:matmulNaive "
"\\\n --launch-skip %d --launch-count 1 "
"--clock-control base \\\n"
" --export profile/day42-naive ./matmul_profile\n",
kLaunchSkip);
std::printf(
" sudo ncu --set full --kernel-name regex:matmulTiled "
"\\\n --launch-skip %d --launch-count 1 "
"--clock-control base \\\n"
" --export profile/day42-tiled ./matmul_profile\n",
kLaunchSkip);
}
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_c));
// 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 (badKernel != nullptr) {
std::fprintf(stderr,
"%s wrong at %zu (row %zu, col %zu); no report should "
"be captured from a kernel that does not agree with the "
"CPU reference\n",
badKernel, badIndex, badIndex / kDim, badIndex % kDim);
return EXIT_FAILURE;
}
std::printf("\nboth kernels match the CPU reference at every element\n");
return EXIT_SUCCESS;
}