code/day74-cp-async/pipelined_gemm.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 74: a two-stage pipelined GEMM with cp.async.
//
// One register-tiled matmul, its shared-memory tiles filled two ways:
//
// sync every thread pulls its float4 from global memory into registers
// and stores it to shared memory, then the block barriers. Load,
// wait, compute, repeat: the memory system and the FMA units take
// turns.
// piped the same float4 travels global-to-shared through cp.async,
// never touching a register, behind a two-stage cuda::pipeline.
// While the block computes on stage s, stage s+1's tiles are in
// flight.
//
// Both kernels share one computeTile() function, so they execute the same
// floating-point operations in the same order and their outputs must be
// bit-identical. The fill is the only difference, so the fill is the only
// thing the ratio can measure.
//
// This program does not run on the course's Tesla T4 verification node.
// cp.async requires compute capability 8.0 ("Requires sm_80 or higher.",
// PTX ISA, cp.async Target ISA Notes, CUDA 12.6). cuda::memcpy_async
// itself works from CC 7.0 through a register-path fallback, but measuring
// the fallback measures nothing this lesson teaches, so the program checks
// for 8.0 and says so rather than reporting a number that means the wrong
// thing.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_80 -lineinfo -o pipelined_gemm \
// pipelined_gemm.cu
// Run: ./pipelined_gemm
//
// -arch=sm_80 embeds compute_80 PTX that JIT-compiles forward, so any
// RTX 30 (8.6), RTX 40 (8.9), RTX 50 (12.0), A100, L4 or H100 runs this
// binary as built.
//
// UNVERIFIED: compiled, never run. nvcc 12.6.2 accepts this file at both
// -arch=sm_80 and -arch=sm_75 (evidence/compile-2026-09-01.txt), but no GPU
// has executed it: the project's Tesla T4 node is CC 7.5 and the capability
// gate below refuses it. Do not publish any output as this program's output
// until a CC 8.0 run exists. See research/REVIEW-PROCESS.md.
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <vector>
#include <cooperative_groups.h>
#include <cuda/pipeline>
#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)
// The tile shapes. A 64 x 64 block tile over a K-depth of 16, computed by
// 256 threads holding a 4 x 4 micro-tile each: day 44's design one size
// down, small enough that every constant below is checkable by eye.
constexpr int kBlockM = 64;
constexpr int kBlockN = 64;
constexpr int kBlockK = 16;
constexpr int kThreadM = 4;
constexpr int kThreadN = 4;
constexpr int kThreadsPerBlock =
(kBlockM / kThreadM) * (kBlockN / kThreadN); // 256, 8 warps
// The exercise raises kStages to 3 and 4. Every static_assert below holds
// for all three values. The upper bound of 4 is the exercise's, not the
// hardware's: at 5 stages ptxas still reports 41048 bytes of shared memory
// per block, inside the 48 KiB static ceiling (compile-only capture,
// evidence/compile-2026-09-01.txt). Six stages is what the ceiling itself
// refuses, at 6 * 8192 + 40 = 49192 bytes against 49152.
constexpr int kStages = 2;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
// Correctness sizes. Every size is a whole number of block tiles in M, N
// and K, because neither kernel guards its bounds (the guard would sit in
// the inner loop, the one place this file cannot pay for it).
constexpr int kSizes[] = {512, 1024, 2048};
constexpr int kSizeCount = 3;
// FP32 accumulate over K = 2048 terms against a double reference. The
// rounding bound 4 * 2^-23 * sqrt(K) is 2.2e-5 at the largest size; 1e-4
// gives slack for a different-but-valid ordering without passing a broken
// kernel.
constexpr double kRelTolerance = 1e-4;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert(kBlockM % kThreadM == 0 && kBlockN % kThreadN == 0,
"micro-tiles must tile the block tile exactly");
static_assert((kBlockM * kBlockK) % (kThreadsPerBlock * 4) == 0 &&
(kBlockK * kBlockN) % (kThreadsPerBlock * 4) == 0,
"each thread moves whole float4s of both tiles");
static_assert(kBlockK % 4 == 0 && kBlockN % 4 == 0,
"float4 moves need rows that are whole numbers of float4s");
static_assert(kStages >= 2 && kStages <= 4,
"one stage cannot overlap; the exercise sweeps 2, 3 and 4");
static_assert(kStages * (kBlockM * kBlockK + kBlockK * kBlockN) *
sizeof(float) <=
48 * 1024,
"staged tiles must fit the 48 KiB static shared ceiling");
// Everything below assumes square matrices that the tiles divide.
constexpr bool tilesDivide(int s) {
return s % kBlockM == 0 && s % kBlockN == 0 && s % kBlockK == 0;
}
static_assert(tilesDivide(kSizes[0]) && tilesDivide(kSizes[1]) &&
tilesDivide(kSizes[2]),
"every size must be a whole number of block tiles");
// The inner product both kernels share. One thread reads its 4 + 4 tile
// values into registers and feeds 16 FMAs per K step, day 44's ratio at
// 0.5 loads per FMA. Because both kernels call exactly this function on
// tiles holding the same values, their outputs are bit-identical by
// construction, and main() checks it with memcmp rather than a tolerance.
__device__ inline void computeTile(const float* tA, const float* tB,
int threadRow, int threadCol,
float acc[kThreadM][kThreadN]) {
for (int p = 0; p < kBlockK; ++p) {
float regM[kThreadM];
float regN[kThreadN];
for (int i = 0; i < kThreadM; ++i) {
regM[i] = tA[(threadRow * kThreadM + i) * kBlockK + p];
}
for (int j = 0; j < kThreadN; ++j) {
regN[j] = tB[p * kBlockN + threadCol * kThreadN + j];
}
for (int i = 0; i < kThreadM; ++i) {
for (int j = 0; j < kThreadN; ++j) {
acc[i][j] += regM[i] * regN[j];
}
}
}
}
// The baseline. One thread owns a 4 x 4 micro-tile of C; the block owns a
// 64 x 64 tile filled one 16-deep K slice at a time.
//
// Memory: each thread moves one float4 of the A slice and one of the B
// slice per K step. Consecutive threads read consecutive 16-byte chunks,
// so a warp's loads coalesce; each float4 crosses global memory, lands in
// a register, and is stored to shared memory, and the whole block then
// waits at the barrier before any FMA issues. Load, wait, compute: the
// pattern this file exists to break.
//
// Launch: exactly kThreadsPerBlock threads, grid (size / kBlockN,
// size / kBlockM), size a whole number of tiles in all three dimensions.
__global__ void matmulTileSync(const float* __restrict__ a,
const float* __restrict__ b,
float* __restrict__ c, int size) {
__shared__ alignas(16) float tileA[kBlockM * kBlockK];
__shared__ alignas(16) float tileB[kBlockK * kBlockN];
const unsigned int tid = threadIdx.x;
const int rowBase = static_cast<int>(blockIdx.y) * kBlockM;
const int colBase = static_cast<int>(blockIdx.x) * kBlockN;
const int threadRow = static_cast<int>(tid) / (kBlockN / kThreadN);
const int threadCol = static_cast<int>(tid) % (kBlockN / kThreadN);
// A slice: 4 float4s per row, so thread tid owns row tid / 4. B slice:
// 16 float4s per row, row tid / 16. Both mappings reappear verbatim in
// the piped kernel; only the transport changes.
const int rowA = static_cast<int>(tid) / (kBlockK / 4);
const int colA = (static_cast<int>(tid) % (kBlockK / 4)) * 4;
const int rowB = static_cast<int>(tid) / (kBlockN / 4);
const int colB = (static_cast<int>(tid) % (kBlockN / 4)) * 4;
float acc[kThreadM][kThreadN] = {};
// snippet: sync-loop
for (int kt = 0; kt < size / kBlockK; ++kt) {
const int kBase = kt * kBlockK;
*reinterpret_cast<float4*>(&tileA[rowA * kBlockK + colA]) =
*reinterpret_cast<const float4*>(
&a[(rowBase + rowA) * size + kBase + colA]);
*reinterpret_cast<float4*>(&tileB[rowB * kBlockN + colB]) =
*reinterpret_cast<const float4*>(
&b[(kBase + rowB) * size + colBase + colB]);
__syncthreads();
computeTile(tileA, tileB, threadRow, threadCol, acc);
__syncthreads();
}
// end snippet
for (int i = 0; i < kThreadM; ++i) {
for (int j = 0; j < kThreadN; ++j) {
c[(rowBase + threadRow * kThreadM + i) * size + colBase +
threadCol * kThreadN + j] = acc[i][j];
}
}
}
// The same matmul with the fill routed through cp.async behind a
// two-stage cuda::pipeline. One thread still owns a 4 x 4 micro-tile and
// still moves one float4 of each slice per K step, at the same addresses
// as the sync kernel.
//
// Memory: cuda::memcpy_async of 16 aligned bytes compiles to one cp.async
// instruction, global to shared with no register in between. The inner
// for-loop keeps kStages fetches in flight: while computeTile() runs on
// stage `compute % kStages`, the copies for the next stage are already
// moving. producer_acquire blocks only when all kStages stage slots are
// full, which is what bounds the lookahead.
//
// Launch: identical to matmulTileSync. Same grid, same block, same sizes.
__global__ void matmulTilePiped(const float* __restrict__ a,
const float* __restrict__ b,
float* __restrict__ c, int size) {
// snippet: pipe-state
// cp.async's 16-byte form needs 16-byte-aligned shared addresses; a
// bare float array only promises 4.
__shared__ alignas(16) float tileA[kStages][kBlockM * kBlockK];
__shared__ alignas(16) float tileB[kStages][kBlockK * kBlockN];
__shared__
cuda::pipeline_shared_state<cuda::thread_scope::thread_scope_block, kStages>
state;
auto block = cooperative_groups::this_thread_block();
auto pipe = cuda::make_pipeline(block, &state);
// end snippet
const unsigned int tid = threadIdx.x;
const int rowBase = static_cast<int>(blockIdx.y) * kBlockM;
const int colBase = static_cast<int>(blockIdx.x) * kBlockN;
const int threadRow = static_cast<int>(tid) / (kBlockN / kThreadN);
const int threadCol = static_cast<int>(tid) % (kBlockN / kThreadN);
const int rowA = static_cast<int>(tid) / (kBlockK / 4);
const int colA = (static_cast<int>(tid) % (kBlockK / 4)) * 4;
const int rowB = static_cast<int>(tid) / (kBlockN / 4);
const int colB = (static_cast<int>(tid) % (kBlockN / 4)) * 4;
float acc[kThreadM][kThreadN] = {};
// snippet: piped-loop
const int tiles = size / kBlockK;
for (int compute = 0, fetch = 0; compute < tiles; ++compute) {
for (; fetch < tiles && fetch < compute + kStages; ++fetch) {
pipe.producer_acquire();
const int s = fetch % kStages;
const int kBase = fetch * kBlockK;
cuda::memcpy_async(&tileA[s][rowA * kBlockK + colA],
&a[(rowBase + rowA) * size + kBase + colA],
cuda::aligned_size_t<16>(sizeof(float4)), pipe);
cuda::memcpy_async(&tileB[s][rowB * kBlockN + colB],
&b[(kBase + rowB) * size + colBase + colB],
cuda::aligned_size_t<16>(sizeof(float4)), pipe);
pipe.producer_commit();
}
pipe.consumer_wait();
computeTile(tileA[compute % kStages], tileB[compute % kStages],
threadRow, threadCol, acc);
pipe.consumer_release();
}
// end snippet
for (int i = 0; i < kThreadM; ++i) {
for (int j = 0; j < kThreadN; ++j) {
c[(rowBase + threadRow * kThreadM + i) * size + colBase +
threadCol * kThreadN + j] = acc[i][j];
}
}
}
// CPU reference. Written for obvious correctness, not speed: plain loops,
// double accumulation. It never allocates; the caller owns every buffer.
static void matmulCpu(const float* a, const float* b, double* out, size_t s) {
for (size_t i = 0; i < s * s; ++i) {
out[i] = 0.0;
}
for (size_t row = 0; row < s; ++row) {
for (size_t p = 0; p < s; ++p) {
const double av = static_cast<double>(a[row * s + p]);
for (size_t col = 0; col < s; ++col) {
out[row * s + col] += av * static_cast<double>(b[p * s + col]);
}
}
}
}
// Returns the first index where got and want differ by more than the
// relative tolerance, or n if they agree everywhere.
static size_t firstMismatch(const float* got, const double* want, size_t n,
double relTolerance) {
for (size_t i = 0; i < n; ++i) {
const double scale = (want[i] == 0.0) ? 1.0 : std::fabs(want[i]);
if (std::fabs(static_cast<double>(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 carried over from the
// course's standard timing helper; copy it verbatim.
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() {
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);
// cp.async exists from CC 8.0. Below that, cuda::memcpy_async still
// compiles and runs, through registers, so a run would print a ratio
// that measures the fallback rather than the feature. Failing here,
// before any allocation, beats publishing that number.
if (prop.major < 8) {
std::fprintf(stderr,
"cp.async needs compute capability 8.0 or newer; %s is "
"%d.%d\n",
prop.name, prop.major, prop.minor);
return EXIT_FAILURE;
}
std::printf("tiles %dx%dx%d, micro-tile %dx%d, %d threads, %d stages\n",
kBlockM, kBlockN, kBlockK, kThreadM, kThreadN, kThreadsPerBlock,
kStages);
bool ok = true;
for (int sizeIdx = 0; sizeIdx < kSizeCount; ++sizeIdx) {
const int size = kSizes[sizeIdx];
const size_t elems = static_cast<size_t>(size) * size;
const size_t bytes = elems * sizeof(float);
std::vector<float> h_a(elems);
std::vector<float> h_b(elems);
std::vector<float> h_sync(elems);
std::vector<float> h_piped(elems);
std::vector<double> h_want(elems);
for (size_t i = 0; i < elems; ++i) {
h_a[i] = static_cast<float>(i % 97) * 0.01f;
h_b[i] = static_cast<float>(i % 53) * 0.02f - 0.5f;
}
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));
const dim3 grid(static_cast<unsigned int>(size / kBlockN),
static_cast<unsigned int>(size / kBlockM));
// Correctness before any timing: the sync kernel against the CPU
// reference, then the piped kernel bit-identical to the sync one.
// The two kernels share computeTile(), so any memcmp difference
// means the fill delivered different bytes, which is exactly the
// bug class this day is about.
matmulTileSync<<<grid, kThreadsPerBlock>>>(d_a, d_b, d_c, size);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_sync.data(), d_c, bytes, cudaMemcpyDeviceToHost));
matmulCpu(h_a.data(), h_b.data(), h_want.data(),
static_cast<size_t>(size));
// Per size, so a failure at 512 still reports 1024 and 2048 rather
// than silently skipping the checks that would localize the bug.
bool sizeOk = true;
const size_t bad =
firstMismatch(h_sync.data(), h_want.data(), elems, kRelTolerance);
if (bad != elems) {
std::fprintf(stderr,
"N=%d sync wrong at %zu: got %.9g, want %.9g\n", size,
bad, h_sync[bad], h_want[bad]);
sizeOk = false;
}
matmulTilePiped<<<grid, kThreadsPerBlock>>>(d_a, d_b, d_c, size);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_piped.data(), d_c, bytes, cudaMemcpyDeviceToHost));
if (sizeOk && std::memcmp(h_sync.data(), h_piped.data(), bytes) != 0) {
std::fprintf(stderr,
"N=%d piped output differs from sync; the fill must "
"not change results\n",
size);
sizeOk = false;
}
ok = ok && sizeOk;
if (sizeOk) {
const float msSync = timeKernel([&] {
matmulTileSync<<<grid, kThreadsPerBlock>>>(d_a, d_b, d_c, size);
});
const float msPiped = timeKernel([&] {
matmulTilePiped<<<grid, kThreadsPerBlock>>>(d_a, d_b, d_c,
size);
});
const double flop = 2.0 * static_cast<double>(elems) * size;
std::printf(
"N=%d sync %.3f ms (%.1f GFLOP/s) piped %.3f ms "
"(%.1f GFLOP/s) piped/sync %.3f\n",
size, msSync, flop / (msSync * 1e-3) / 1e9, msPiped,
flop / (msPiped * 1e-3) / 1e9, msPiped / msSync);
}
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_c));
}
return ok ? EXIT_SUCCESS : EXIT_FAILURE;
}