code/day75-async-barriers/barrier_pipeline.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 75: asynchronous barriers.
//
// One producer-consumer tile loop, two warp groups, written twice:
//
// pipelineSyncthreads the same schedule held together by __syncthreads(),
// one all-stop per tile. Every warp waits at every
// boundary, whether or not it has anything to wait
// for.
// pipelineSplitBarrier four cuda::barrier objects in shared memory, the
// producer-consumer pattern from section 7.26.5 of
// the CUDA 12.6 programming guide. Producers never
// wait on "filled", consumers never wait on "ready".
//
// Both kernels stage tiles through a double-buffered shared array (day 54's
// trade, inside one block) and compute the same per-element FMA chain, so
// their outputs must match bit for bit before anything is timed.
//
// The Tesla T4 verification node compiles this file and cannot run it: the
// binary targets sm_80, and the subject is the hardware mbarrier that
// arrived with compute capability 8.0. cuda::barrier itself exists from CC
// 7.0 in software; the program gates on 8.0 because a T4 number would
// measure a software loop, not the hardware unit this day is about.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_80 -o barrier_pipeline \
// barrier_pipeline.cu
// Run: ./barrier_pipeline
//
// -arch=sm_80 embeds compute_80 PTX that JIT-compiles forward, so any CC
// 8.x, 9.0, 10.x or 12.x card runs this binary.
//
// UNVERIFIED: not yet compiled or run on real hardware. Do not publish any
// output as this program's output until it has run. See
// research/REVIEW-PROCESS.md.
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <vector>
#include <cuda/barrier>
#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)
// kElems is deliberately not a multiple of the tile size, so the ragged
// last tile exercises every guard on every run. Two float buffers at this
// size is 32 MiB, small enough for any CC 8.0 card.
constexpr size_t kElems = 4ull * 1024ull * 1024ull + 611ull;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kProducerWarps = 4; // the exercise's knob
constexpr int kProducerThreads = kProducerWarps * 32;
constexpr int kConsumerThreads = kThreadsPerBlock - kProducerThreads;
constexpr int kTileElems = 1024; // 4 KiB per buffer, two buffers
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr int kInnerIters = 64; // FMA chain length per element
constexpr float kMul = 0.999f;
constexpr float kAdd = 0.25f;
constexpr float kRelTolerance = 1e-5f;
static_assert(kProducerThreads % 32 == 0 && kProducerThreads > 0,
"the producer group must be whole warps");
static_assert(kConsumerThreads % 32 == 0 && kConsumerThreads > 0,
"the consumer group must be whole warps");
using barrier = cuda::barrier<cuda::thread_scope_block>;
// Each element costs kInnerIters dependent FMAs, computed identically on
// host (std::fmaf) and device (fmaf). Both are correctly rounded, so the
// reference match is exact and the tolerance below only covers habits.
__host__ __device__ inline float fmaChain(float v, float mul, float add) {
for (int k = 0; k < kInnerIters; ++k) {
v = fmaf(v, mul, add);
}
return v;
}
// Producer work: thread pt of the producer group copies elements pt,
// pt + kProducerThreads, ... of one tile from global into a shared buffer.
// Consecutive producer threads touch consecutive addresses, so the global
// loads coalesce and the shared stores are conflict-free.
__device__ void loadTile(const float* __restrict__ in, float* __restrict__ buf,
size_t tile, size_t n) {
const size_t base = tile * kTileElems;
for (int e = static_cast<int>(threadIdx.x); e < kTileElems;
e += kProducerThreads) {
const size_t g = base + static_cast<size_t>(e);
buf[e] = (g < n) ? in[g] : 0.0f;
}
}
// Consumer work: the same strided walk over the tile, offset into the
// consumer group, one FMA chain per element, result to global. The guard
// covers the store; out past n is never written.
__device__ void consumeTile(const float* __restrict__ buf,
float* __restrict__ out, size_t tile, size_t n) {
const size_t base = tile * kTileElems;
for (int e = static_cast<int>(threadIdx.x) - kProducerThreads;
e < kTileElems; e += kConsumerThreads) {
const size_t g = base + static_cast<size_t>(e);
if (g < n) {
out[g] = fmaChain(buf[e], kMul, kAdd);
}
}
}
// snippet: sync-loop
// The baseline. One thread group loads, the other computes, double
// buffered: while consumers drain tiles[cur], producers fill the other
// buffer with the block's next tile. One __syncthreads() per tile closes
// both hazards at once, and stops all 8 warps to do it.
//
// Launch requirement: gridDim.x must not exceed the tile count, so every
// block has a first tile and the prologue fill is never out of range.
// Every thread reaches every __syncthreads(): the tile loop's trip count
// depends only on blockIdx, and no barrier sits inside a divergent branch.
__global__ void pipelineSyncthreads(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tiles[2][kTileElems];
const bool producer = threadIdx.x < kProducerThreads;
const size_t numTiles = (n + kTileElems - 1) / kTileElems;
if (producer) {
loadTile(in, tiles[0], blockIdx.x, n);
}
__syncthreads();
size_t j = 0;
for (size_t tile = blockIdx.x; tile < numTiles; tile += gridDim.x, ++j) {
const int cur = static_cast<int>(j & 1);
if (producer) {
const size_t next = tile + gridDim.x;
if (next < numTiles) {
loadTile(in, tiles[cur ^ 1], next, n);
}
} else {
consumeTile(tiles[cur], out, tile, n);
}
__syncthreads();
}
}
// end snippet
// snippet: split-loop
// The same schedule on four cuda::barrier objects, the producer-consumer
// pattern of programming guide section 7.26.5. ready[b] means "tiles[b]
// may be overwritten"; filled[b] means "tiles[b] holds a complete tile".
// Each wait is one-sided: producers never wait on filled, consumers never
// wait on ready. arrive() returns a token that is the right to wait on
// that phase; a group that will never wait drops it, and the (void) says
// the drop is deliberate.
//
// All 256 threads participate in all four barriers, which is why each is
// initialized with kThreadsPerBlock: an arrive_and_wait and a bare arrive
// count the same. On sm_80 and newer these barriers are hardware mbarrier
// objects in shared memory.
__global__ void pipelineSplitBarrier(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tiles[2][kTileElems];
__shared__ barrier ready[2];
__shared__ barrier filled[2];
if (threadIdx.x < 2) {
init(&ready[threadIdx.x], kThreadsPerBlock);
init(&filled[threadIdx.x], kThreadsPerBlock);
}
__syncthreads(); // no thread may touch a barrier before init
const size_t numTiles = (n + kTileElems - 1) / kTileElems;
if (threadIdx.x < kProducerThreads) {
size_t j = 0;
for (size_t tile = blockIdx.x; tile < numTiles;
tile += gridDim.x, ++j) {
const int buf = static_cast<int>(j & 1);
ready[buf].arrive_and_wait(); // consumers done with tiles[buf]
loadTile(in, tiles[buf], tile, n);
(void)filled[buf].arrive(); // hand off; do not wait
}
} else {
(void)ready[0].arrive(); // both buffers start out writable
(void)ready[1].arrive();
size_t j = 0;
for (size_t tile = blockIdx.x; tile < numTiles;
tile += gridDim.x, ++j) {
const int buf = static_cast<int>(j & 1);
filled[buf].arrive_and_wait(); // producers filled tiles[buf]
consumeTile(tiles[buf], out, tile, n);
(void)ready[buf].arrive(); // hand back; do not wait
}
}
}
// end snippet
// CPU reference. The staging through shared memory must not change the
// answer, so the reference is one loop with the same chain, no tiles.
static void pipelineCpu(const float* in, float* out, size_t n) {
for (size_t i = 0; i < n; ++i) {
out[i] = fmaChain(in[i], kMul, kAdd);
}
}
// 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 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. Copied verbatim from the course timing helper: warm-up inside, so
// every kernel it times pays its lazy-loading cost before the clock starts.
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
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);
// cuda::barrier runs from CC 7.0, in software. The hardware mbarrier
// this day is about arrived with CC 8.0, and the binary targets sm_80,
// so the gate fails here with a sentence instead of ten lines down
// with cudaErrorNoKernelImageForDevice.
if (prop.major < 8) {
std::fprintf(stderr,
"the hardware mbarrier needs compute capability 8.0 or "
"newer; %s is %d.%d\n",
prop.name, prop.major, prop.minor);
return EXIT_FAILURE;
}
const size_t bytes = kElems * sizeof(float);
const size_t numTiles = (kElems + kTileElems - 1) / kTileElems;
// Two blocks per SM keeps every SM busy without oversubscribing the
// tile list; capped at the tile count so every block owns a first tile.
const size_t wantBlocks = 2ull * prop.multiProcessorCount;
const int blocks =
static_cast<int>(wantBlocks < numTiles ? wantBlocks : numTiles);
std::printf(
"n = %zu, %zu tiles of %d, %d blocks of %d threads "
"(%d producer warps, %d consumer warps)\n",
kElems, numTiles, kTileElems, blocks, kThreadsPerBlock, kProducerWarps,
kThreadsPerBlock / 32 - kProducerWarps);
std::vector<float> h_in(kElems);
std::vector<float> h_want(kElems);
std::vector<float> h_out(kElems);
std::vector<float> h_outSplit(kElems);
for (size_t i = 0; i < kElems; ++i) {
h_in[i] = static_cast<float>(i % 97) * 0.01f;
}
float* d_in = nullptr;
float* d_out = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), bytes, cudaMemcpyHostToDevice));
bool ok = true;
// Correctness first, both kernels, before anything is timed. The
// __syncthreads baseline is checked against the CPU reference; the
// split-barrier kernel must then match the baseline bit for bit,
// because the barriers change the schedule and must change nothing
// else. d_out is zeroed in between so a kernel that silently skipped
// its stores cannot pass by inheriting the baseline's answer.
pipelineSyncthreads<<<blocks, kThreadsPerBlock>>>(d_in, d_out, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
pipelineCpu(h_in.data(), h_want.data(), kElems);
const size_t bad =
firstMismatch(h_out.data(), h_want.data(), kElems, kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr, "syncthreads wrong at %zu: got %.9g, want %.9g\n",
bad, h_out[bad], h_want[bad]);
ok = false;
}
CUDA_CHECK(cudaMemset(d_out, 0, bytes));
pipelineSplitBarrier<<<blocks, kThreadsPerBlock>>>(d_in, d_out, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_outSplit.data(), d_out, bytes, cudaMemcpyDeviceToHost));
if (ok && std::memcmp(h_out.data(), h_outSplit.data(), bytes) != 0) {
std::fprintf(stderr,
"split-barrier output differs from the __syncthreads "
"baseline; the schedule must not change results\n");
ok = false;
}
if (ok) {
// Time both. Same grid, same block, same tile walk, same math; the
// only difference is what holds the two warp groups together.
const float msSync = timeKernel([&] {
pipelineSyncthreads<<<blocks, kThreadsPerBlock>>>(d_in, d_out,
kElems);
});
const float msSplit = timeKernel([&] {
pipelineSplitBarrier<<<blocks, kThreadsPerBlock>>>(d_in, d_out,
kElems);
});
std::printf("mean of %d runs, copies not included\n", kTimedRuns);
std::printf("syncthreads: %.3f ms\n", msSync);
std::printf("split barrier: %.3f ms\n", msSplit);
std::printf("split/sync: %.3f\n", msSplit / msSync);
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_out));
return ok ? EXIT_SUCCESS : EXIT_FAILURE;
}