code/day16-tiled-matmul/tiled_matmul.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 16: tiled matrix multiplication, and the guard it needs on a size the
// tile does not divide.
//
// Two kernels compute the same C = A * B. matmulNaive gives one thread one
// output element and reads a whole row of A and a whole column of B out of
// global memory to produce it: 2K loads per thread. matmulTiled gives one
// block one 16 x 16 tile of C and stages both inputs through shared memory,
// so each thread issues two global loads per tile step instead of two per k
// step: 2 * ceil(K / 16) loads per thread. That ratio is the lesson.
//
// Three sizes, and the third is why this file exists. 611 x 613 x 617
// divides by 16 on no axis, and the three remainders differ: 611 is 38 whole
// tiles and 3 rows, 613 is 38 and 5, 617 is 38 and 9. A kernel whose guard
// is right on one axis and wrong on another cannot hide behind a square.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o tiled_matmul tiled_matmul.cu
// Run: ./tiled_matmul
//
// Verified 2026-08-30 on a Tesla T4 (compute capability 7.5), driver
// 595.84, CUDA 12.6 (V12.6.85). Transcript: evidence/run-2026-08-30.txt
#include <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. A T4 gives a block 48 KiB by default
// and 64 KiB with an opt-in, so shared memory is nowhere near binding here.
// Day 45 is where it starts to.
constexpr int kTileDim = 16;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kRelTolerance = 1e-5f;
// Case 0 and case 1 divide by the tile on every axis. Case 2 divides on none
// of them, with a different remainder on each, so the guard is exercised in
// the row direction, the column direction and the k direction at once.
constexpr int kNumCases = 3;
constexpr size_t kCaseM[kNumCases] = {256, 512, 611};
constexpr size_t kCaseN[kNumCases] = {256, 512, 613};
constexpr size_t kCaseK[kNumCases] = {256, 512, 617};
// Every buffer is allocated once at kMaxDim x kMaxDim floats, about 1.5 MiB,
// so no allocation ever sits inside a timed region and one bound covers all
// three cases.
constexpr size_t kMaxDim = 617;
static_assert(kTileDim * kTileDim == 256,
"a 16 x 16 tile is one 256-thread block, the course default");
static_assert((kTileDim * kTileDim) % 32 == 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(kCaseM[0] % kTileDim == 0 && kCaseN[0] % kTileDim == 0 &&
kCaseK[0] % kTileDim == 0,
"case 0 must divide by the tile on every axis, or the table "
"never shows a clean size next to a ragged one");
static_assert(kCaseM[2] % kTileDim != 0 && kCaseN[2] % kTileDim != 0 &&
kCaseK[2] % kTileDim != 0,
"case 2 must divide by the tile on no axis: it is the case a "
"tiled kernel written for multiples gets wrong");
static_assert(kCaseM[0] <= kMaxDim && kCaseN[0] <= kMaxDim &&
kCaseK[0] <= kMaxDim && kCaseM[1] <= kMaxDim &&
kCaseN[1] <= kMaxDim && kCaseK[1] <= kMaxDim &&
kCaseM[2] <= kMaxDim && kCaseN[2] <= kMaxDim &&
kCaseK[2] <= kMaxDim,
"every dimension must fit the single kMaxDim allocation");
// Integer ceiling division. constexpr so one function sizes a grid at run
// time and appears 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;
}
static_assert(ceilDiv(kMaxDim, static_cast<size_t>(kTileDim)) <= 65535,
"gridDim.y stops at 65535, unlike gridDim.x which goes to "
"2^31 - 1");
// 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.
//
// Launch assumption: the grid covers C. It rounds up on both axes, so both
// halves of the guard run on the ragged case. This kernel has no barrier,
// which is the only reason the guard is allowed to wrap the whole body.
// snippet: naive-kernel
__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;
}
}
// end snippet
// The same product, one 16 x 16 tile of C per block, staged through shared
// memory. Each block walks ceil(K / 16) tile steps. On every step each thread
// loads one element of A and one of B, the block synchronises, and each
// thread does 16 multiply-adds out of 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 rows of the table differ in the number of loads and in nothing
// else. In shared memory a warp reads tileA[threadIdx.y][p], two addresses
// the halves of the warp share, and tileB[p][threadIdx.x], sixteen
// consecutive words read by two lanes each. Both are broadcasts or distinct
// banks, so neither conflicts; day 15 is where that stops being free.
//
// Ragged sizes: a thread whose element lies outside A or B writes 0.0f into
// the tile rather than skipping the load. Zero is the identity for the sum,
// and skipping is not: the tile is reused on the next step, so a skipped
// write leaves the previous step's value in shared memory and multiplies it
// into acc.
//
// Launch assumption: exactly kTileDim x kTileDim threads per block, and a
// grid that rounds up on both axes. Every thread reaches both barriers. The
// guard covers the two loads and the store and never the __syncthreads(),
// because a barrier only some of the block arrives at is undefined
// behaviour. Day 14 has the rule and the racecheck output.
// snippet: tiled-kernel
__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;
}
}
// end snippet
// CPU reference. Three nested loops in the obvious order, written for obvious
// correctness rather than speed: no blocking, no OpenMP, no intrinsics. 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. The Kahan
// compensation EXERCISE-DESIGN.md asks for past 64 terms would change nothing
// here: makeInputs picks values small enough that every partial sum is an
// integer under 2^24 and therefore exact in float and in double alike.
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 any case
// can produce is 6 * 617 = 3,702, 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 none of the six dimensions in the table,
// 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 aElems, size_t bElems) {
for (size_t e = 0; e < aElems; ++e) {
(*h_a)[e] = static_cast<float>(static_cast<int>(e % 7) - 3);
}
for (size_t e = 0; e < bElems; ++e) {
(*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 whole point: main() turns it back into a row and a column, and
// "wrong at row 608, column 0" names the ragged tile. "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 a launch with CUDA events and returns the mean milliseconds per run.
//
// This is the one template and the one lambda allowed in module 1 to 3 code.
// 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.
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. The figure is computed here from the case
// dimensions rather than per kernel, because both kernels do exactly this
// much arithmetic and a second copy is a second thing to get wrong.
static double gflops(size_t m, size_t n, size_t k, float ms) {
const double work = 2.0 * static_cast<double>(m) * static_cast<double>(n) *
static_cast<double>(k);
return work / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
static void printRow(size_t m, size_t n, size_t k, const char* kernel, float ms,
size_t loadsPerThread) {
char label[32];
std::snprintf(label, sizeof(label), "%zux%zux%zu", m, n, k);
std::printf("%-16s %-7s %11.3f %11.1f %13zu\n", label, kernel, ms,
gflops(m, n, k, ms), loadsPerThread);
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name, prop.major,
prop.minor);
std::printf(
"tile %d x %d, %d threads per block, %zu bytes of shared memory per "
"block\n\n",
kTileDim, kTileDim, kTileDim * kTileDim,
2 * sizeof(float) * kTileDim * kTileDim);
const size_t maxElems = kMaxDim * kMaxDim;
const size_t maxBytes = maxElems * sizeof(float);
std::vector<float> h_a(maxElems);
std::vector<float> h_b(maxElems);
std::vector<float> h_c(maxElems);
std::vector<float> h_want(maxElems);
float* d_a = nullptr;
float* d_b = nullptr;
float* d_c = nullptr;
CUDA_CHECK(cudaMalloc(&d_a, maxBytes));
CUDA_CHECK(cudaMalloc(&d_b, maxBytes));
CUDA_CHECK(cudaMalloc(&d_c, maxBytes));
std::printf("%-16s %-7s %11s %11s %13s\n", "case", "kernel", "time (ms)",
"GFLOP/s", "loads/thread");
std::printf("%-16s %-7s %11s %11s %13s\n", "----------------", "-------",
"-----------", "-----------", "-------------");
// Failure is recorded rather than returned, so every path falls through
// to the one cleanup block below and no cudaMalloc escapes its cudaFree.
int rows = 0;
int badCase = kNumCases;
const char* badKernel = nullptr;
size_t badIndex = 0;
for (int c = 0; c < kNumCases; ++c) {
const size_t m = kCaseM[c];
const size_t n = kCaseN[c];
const size_t k = kCaseK[c];
const size_t outElems = m * n;
const size_t outBytes = outElems * sizeof(float);
makeInputs(&h_a, &h_b, m * k, k * n);
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), m * k * sizeof(float),
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), k * n * sizeof(float),
cudaMemcpyHostToDevice));
matmulCpu(h_a.data(), h_b.data(), h_want.data(), m, n, k);
// 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(n, tileDim)),
static_cast<unsigned int>(ceilDiv(m, tileDim)));
// The output is cleared before each correctness launch so a kernel
// that skips elements is caught by the comparison rather than by the
// previous case's leftovers agreeing with it.
CUDA_CHECK(cudaMemset(d_c, 0, outBytes));
matmulNaive<<<grid, block>>>(d_a, d_b, d_c, m, n, k);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_c.data(), d_c, outBytes, cudaMemcpyDeviceToHost));
size_t bad =
firstMismatch(h_c.data(), h_want.data(), outElems, kRelTolerance);
if (bad != outElems) {
badCase = c;
badKernel = "naive";
badIndex = bad;
break;
}
const float naiveMs = timeKernel(
[&] { matmulNaive<<<grid, block>>>(d_a, d_b, d_c, m, n, k); });
printRow(m, n, k, "naive", naiveMs, 2 * k);
++rows;
CUDA_CHECK(cudaMemset(d_c, 0, outBytes));
matmulTiled<<<grid, block>>>(d_a, d_b, d_c, m, n, k);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_c.data(), d_c, outBytes, cudaMemcpyDeviceToHost));
bad = firstMismatch(h_c.data(), h_want.data(), outElems, kRelTolerance);
if (bad != outElems) {
badCase = c;
badKernel = "tiled";
badIndex = bad;
break;
}
const float tiledMs = timeKernel(
[&] { matmulTiled<<<grid, block>>>(d_a, d_b, d_c, m, n, k); });
const size_t tiledLoads = 2 * ceilDiv(k, tileDim);
printRow(m, n, k, "tiled", tiledMs, tiledLoads);
++rows;
std::printf(
"%-16s tiled is %.2fx faster and reads %.1fx fewer floats per "
"thread\n",
"", static_cast<double>(naiveMs) / tiledMs,
static_cast<double>(2 * k) / static_cast<double>(tiledLoads));
}
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_c));
if (badCase != kNumCases) {
const size_t n = kCaseN[badCase];
std::fprintf(stderr,
"%s wrong at %zu (row %zu, col %zu) on case %zux%zux%zu\n",
badKernel, badIndex, badIndex / n, badIndex % n,
kCaseM[badCase], n, kCaseK[badCase]);
return EXIT_FAILURE;
}
// 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. It exists so the lesson's table and
// this program cannot drift apart on how many rows there are.
if (rows != 2 * kNumCases) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's table and "
"this program disagree\n",
rows, 2 * kNumCases);
return EXIT_FAILURE;
}
std::printf("\nall %d cases match the CPU reference at every element\n",
kNumCases);
return EXIT_SUCCESS;
}