code/day54-double-buffer/double_buffer.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 54: double buffering. A 1 GiB array streamed through a 64 MiB device
// window, twice: once serially (copy, compute, copy back, one chunk at a
// time), once pipelined (two device buffer pairs, two streams, so chunk k+1
// copies in while chunk k computes and chunk k-1 copies out).
//
// What it measures: the three component times (all-H2D, all-kernel, all-D2H)
// and the two end-to-end variants, so the page can put the naive bound (the
// sum of the components) next to the pipeline floor (the widest component)
// and next to what actually happened.
//
// What it does not measure: PCIe bandwidth in isolation. Day 53 owns that.
// The GB/s figures printed here are derived from the component passes and
// include the per-chunk enqueue overhead a real pipeline pays.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo -o double_buffer \
// double_buffer.cu
// Run: ./double_buffer (needs about 3 GiB of host RAM, 2 GiB of it
// pinned, and 256 MiB of device memory)
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <vector>
#include <cuda_runtime.h>
#include <nvtx3/nvToolsExt.h>
// The one error macro. This file is standalone, so it carries its own
// verbatim copy. `err_` has a trailing underscore so it cannot collide with a
// variable at the call site.
#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)
// 1 GiB of floats through a 64 MiB window: 16 chunks. The chunk size is the
// exercise's knob; 4 MiB and 256 MiB windows divide the total too, so the
// static_asserts below hold at every size the page invites, and at block
// size 32.
constexpr size_t kTotalElems = 256 * 1024 * 1024; // 1 GiB of float
constexpr size_t kChunkElems = 16 * 1024 * 1024; // 64 MiB window
constexpr int kChunks = static_cast<int>(kTotalElems / kChunkElems);
constexpr size_t kChunkBytes = kChunkElems * sizeof(float);
constexpr size_t kTotalBytes = kTotalElems * sizeof(float);
constexpr int kThreadsPerBlock = 256; // 8 warps
static_assert(kTotalElems % kChunkElems == 0,
"the window must divide the array");
static_assert(kChunkElems % kThreadsPerBlock == 0,
"chunk size must be a whole number of blocks");
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
constexpr int kBlocksPerChunk =
static_cast<int>(kChunkElems / kThreadsPerBlock);
// kIters chained FMAs per element give the kernel a cost the copies can
// hide. With one FMA per element the kernel stage would be a sliver and the
// pipeline would only ever overlap the two copy directions.
constexpr int kIters = 512;
constexpr float kA = 0.999f;
constexpr float kB = 0.0625f;
// The reference is the closed form of the same recurrence in double, exact
// where the kernel's 512 sequential float FMAs each round once. |kA| < 1
// makes every step contractive, so the accumulated error stays within a few
// hundred ULP; 1e-4 clears that with margin and still catches a chunk that
// was skipped, doubled or read before its copy landed.
constexpr float kRelTolerance = 1e-4f;
// One pass moves 2 GiB across PCIe and runs for hundreds of milliseconds,
// so event resolution and launch jitter are irrelevant at this scale. One
// warm pass (which also pays lazy module loading for the one kernel) and
// the mean of five is enough; the course's 3-and-10 rule is for kernels a
// thousand times shorter.
constexpr int kTimedRuns = 5;
// One thread owns one element of one chunk. A warp's 32 loads and 32 stores
// are consecutive floats, 128 contiguous bytes each way, fully coalesced.
// Launch assumption: gridDim.x * blockDim.x >= n.
// snippet: kernel
__global__ void transformChunk(const float* in, float* out, size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
float y = in[i];
for (int k = 0; k < kIters; ++k) {
y = fmaf(y, kA, kB);
}
out[i] = y;
}
}
// end snippet
// Everything a pass needs. Buffer pair s belongs to stream s and to nobody
// else; stream order is the only thing that keeps a buffer from being
// overwritten while in use, and it is enough.
struct Ctx {
float* h_in;
float* h_out;
float* d_in[2];
float* d_out[2];
cudaStream_t stream[2];
};
// The third launch argument is dynamic shared memory, none here; it exists
// only to reach the fourth, the stream, which day 51 introduced.
static void launchChunk(const float* in, float* out, cudaStream_t s) {
transformChunk<<<kBlocksPerChunk, kThreadsPerBlock, 0, s>>>(in, out,
kChunkElems);
}
// All 16 H2D copies, nothing else. If asyncEngineCount reports 2, each
// PCIe direction gets its own copy engine, which is prediction 4's bet on
// the page, not a given; either way two streams share the one H2D engine,
// so this is the serial copy-in time whatever the stream count.
static void passH2D(const Ctx& c) {
for (int k = 0; k < kChunks; ++k) {
const size_t off = static_cast<size_t>(k) * kChunkElems;
CUDA_CHECK(cudaMemcpyAsync(c.d_in[k & 1], c.h_in + off, kChunkBytes,
cudaMemcpyHostToDevice, c.stream[k & 1]));
}
}
// All 16 kernel launches on resident data, one stream, no copies. One chunk
// is 65,536 blocks, which fills the whole card, so serializing them on one
// stream costs nothing and gives the clean serial compute time.
static void passKernel(const Ctx& c) {
for (int k = 0; k < kChunks; ++k) {
launchChunk(c.d_in[0], c.d_out[0], c.stream[0]);
}
}
// All 16 D2H copies, nothing else.
static void passD2H(const Ctx& c) {
for (int k = 0; k < kChunks; ++k) {
const size_t off = static_cast<size_t>(k) * kChunkElems;
CUDA_CHECK(cudaMemcpyAsync(c.h_out + off, c.d_out[k & 1], kChunkBytes,
cudaMemcpyDeviceToHost, c.stream[k & 1]));
}
}
// The upper bound: one buffer pair, one stream. Stream order serializes the
// three stages of every chunk against each other and against the next
// chunk, so this pass costs the sum of the three component passes without
// a single blocking call.
// snippet: naive
static void passNaive(const Ctx& c) {
for (int k = 0; k < kChunks; ++k) {
const size_t off = static_cast<size_t>(k) * kChunkElems;
CUDA_CHECK(cudaMemcpyAsync(c.d_in[0], c.h_in + off, kChunkBytes,
cudaMemcpyHostToDevice, c.stream[0]));
launchChunk(c.d_in[0], c.d_out[0], c.stream[0]);
CUDA_CHECK(cudaMemcpyAsync(c.h_out + off, c.d_out[0], kChunkBytes,
cudaMemcpyDeviceToHost, c.stream[0]));
}
}
// end snippet
// The pipeline. Even chunks own buffer pair 0 and stream 0, odd chunks own
// pair 1 and stream 1. Within a stream, order still protects the buffers:
// chunk k+2 cannot overwrite d_in[k & 1] until chunk k's copy out has
// drained, because both sit in the same stream. Across streams nothing
// waits, which is the overlap.
// snippet: pipeline
static void passPipelined(const Ctx& c) {
for (int k = 0; k < kChunks; ++k) {
const int s = k & 1;
const size_t off = static_cast<size_t>(k) * kChunkElems;
CUDA_CHECK(cudaMemcpyAsync(c.d_in[s], c.h_in + off, kChunkBytes,
cudaMemcpyHostToDevice, c.stream[s]));
launchChunk(c.d_in[s], c.d_out[s], c.stream[s]);
CUDA_CHECK(cudaMemcpyAsync(c.h_out + off, c.d_out[s], kChunkBytes,
cudaMemcpyDeviceToHost, c.stream[s]));
}
}
// end snippet
// Times a whole pass with events on the legacy default stream. Both streams
// are created blocking on purpose: an event recorded on the legacy stream
// then acts as a barrier, so `start` fires before any timed work is queued
// and `stop` completes only after both streams drain. That is the sync
// point; cudaEventSynchronize(stop) is where the host waits before reading
// the clock. One warm pass covers the program's one kernel and wakes both
// copy engines. No sync inside the timed loop: a drain per pass would
// remeasure the ramp-up 5 times.
static float timePassMs(const Ctx& c, void (*pass)(const Ctx&),
const char* label) {
nvtxRangePushA(label);
pass(c); // warm-up
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
cudaEvent_t start;
cudaEvent_t stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
CUDA_CHECK(cudaEventRecord(start));
for (int r = 0; r < kTimedRuns; ++r) {
pass(c);
}
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));
nvtxRangePop();
return ms / kTimedRuns;
}
// After kIters steps of y = fma(y, a, b), y = x * a^n + b * (1 - a^n) /
// (1 - a). Computed in double from the float constants, so the only
// difference from the kernel is the kernel's per-step rounding.
static void expectedCoeffs(double* scale, double* offset) {
const double a = static_cast<double>(kA);
*scale = std::pow(a, kIters);
*offset = static_cast<double>(kB) * (1.0 - *scale) / (1.0 - a);
}
// Index of the first element outside tolerance, or n if all pass. The
// expected values sit near 25, never near zero, so a pure relative
// comparison has no cancellation trap.
static size_t firstMismatch(const float* got, const float* in, size_t n) {
double scale = 0.0;
double offset = 0.0;
expectedCoeffs(&scale, &offset);
for (size_t i = 0; i < n; ++i) {
const double want = scale * static_cast<double>(in[i]) + offset;
const double diff = std::fabs(static_cast<double>(got[i]) - want);
if (diff > static_cast<double>(kRelTolerance) * std::fabs(want)) {
return i;
}
}
return n;
}
static void freeAll(Ctx* c) {
for (int s = 0; s < 2; ++s) {
CUDA_CHECK(cudaFree(c->d_in[s]));
CUDA_CHECK(cudaFree(c->d_out[s]));
CUDA_CHECK(cudaStreamDestroy(c->stream[s]));
}
CUDA_CHECK(cudaFreeHost(c->h_in));
CUDA_CHECK(cudaFreeHost(c->h_out));
}
int main() {
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
// asyncEngineCount is an int: how many copies can move while a kernel
// runs. 2 means H2D and D2H can also overlap each other, which is what
// prediction 4 on the page bets on.
std::printf("GPU: %s (compute capability %d.%d), %d copy engine(s)\n",
prop.name, prop.major, prop.minor, prop.asyncEngineCount);
std::printf(
"n = %zu floats, %zu MiB total, %d chunks of %zu MiB, "
"%d threads/block\n\n",
kTotalElems, kTotalBytes >> 20, kChunks, kChunkBytes >> 20,
kThreadsPerBlock);
Ctx c = {};
nvtxRangePushA("setup");
// Pinned on both ends: cudaMemcpyAsync from pageable memory silently
// degrades to a staged, effectively synchronous copy and the whole
// pipeline flattens. Day 53 measured the difference; this program
// depends on it.
CUDA_CHECK(cudaMallocHost(&c.h_in, kTotalBytes));
CUDA_CHECK(cudaMallocHost(&c.h_out, kTotalBytes));
for (int s = 0; s < 2; ++s) {
CUDA_CHECK(cudaMalloc(&c.d_in[s], kChunkBytes));
CUDA_CHECK(cudaMalloc(&c.d_out[s], kChunkBytes));
CUDA_CHECK(cudaStreamCreate(&c.stream[s]));
}
// 4093 is prime, so the pattern never lines up with a chunk boundary;
// values span [-0.5, 0.5).
for (size_t i = 0; i < kTotalElems; ++i) {
c.h_in[i] = static_cast<float>(i % 4093) / 4093.0f - 0.5f;
}
nvtxRangePop();
// Correctness before any timing. The naive pass is the reference
// implementation; its output is checked against the closed form, kept,
// and then the pipelined pass must reproduce it bit for bit. A buffer
// reuse race would show up here as a mismatched chunk, not as noise.
nvtxRangePushA("verify-naive");
passNaive(c);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
size_t bad = firstMismatch(c.h_out, c.h_in, kTotalElems);
if (bad != kTotalElems) {
std::fprintf(stderr, "FAIL: naive output wrong at i=%zu, got %.9g\n",
bad, static_cast<double>(c.h_out[bad]));
freeAll(&c);
return EXIT_FAILURE;
}
std::vector<float> h_want(c.h_out, c.h_out + kTotalElems);
nvtxRangePop();
nvtxRangePushA("verify-pipeline");
passPipelined(c);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
bad = firstMismatch(c.h_out, c.h_in, kTotalElems);
if (bad != kTotalElems) {
std::fprintf(stderr,
"FAIL: pipelined output wrong at i=%zu, got %.9g\n", bad,
static_cast<double>(c.h_out[bad]));
freeAll(&c);
return EXIT_FAILURE;
}
if (std::memcmp(c.h_out, h_want.data(), kTotalBytes) != 0) {
size_t i = 0;
while (i < kTotalElems &&
std::memcmp(&c.h_out[i], &h_want[i], sizeof(float)) == 0) {
++i;
}
std::fprintf(stderr,
"FAIL: pipelined differs from naive at i=%zu "
"(chunk %zu): got %.9g, want %.9g\n",
i, i / kChunkElems, static_cast<double>(c.h_out[i]),
static_cast<double>(h_want[i]));
freeAll(&c);
return EXIT_FAILURE;
}
std::printf(
"verify: naive matches closed form, pipelined matches "
"naive bit for bit\n\n");
nvtxRangePop();
// Components first, then the two ends of the argument. Copy GB/s here
// is 1 GiB divided by the pass time: the streamed rate through the
// window, including per-chunk enqueue overhead, not a peak.
const float msH2D = timePassMs(c, passH2D, "h2d-only");
const float msKernel = timePassMs(c, passKernel, "kernel-only");
const float msD2H = timePassMs(c, passD2H, "d2h-only");
const float msNaive = timePassMs(c, passNaive, "naive");
const float msPipe = timePassMs(c, passPipelined, "double-buffer");
const double gib = static_cast<double>(kTotalBytes) / 1.0e9;
std::printf("components, mean of %d passes over the full 1 GiB\n",
kTimedRuns);
std::printf(" h2d only %10.3f ms %6.2f GB/s\n", msH2D,
gib / (msH2D / 1000.0));
std::printf(" kernel only %10.3f ms\n", msKernel);
std::printf(" d2h only %10.3f ms %6.2f GB/s\n", msD2H,
gib / (msD2H / 1000.0));
const float msSum = msH2D + msKernel + msD2H;
const float msMax = std::fmax(msH2D, std::fmax(msKernel, msD2H));
std::printf("\nnaive bound (sum of components) %10.3f ms\n", msSum);
std::printf("pipeline floor (widest component) %10.3f ms\n", msMax);
std::printf("\nmeasured naive %10.3f ms\n", msNaive);
std::printf("measured double buffer %10.3f ms\n", msPipe);
std::printf(
"\nspeedup: %.2fx of naive; edge cost above the floor: "
"%.3f ms\n",
msNaive / msPipe, msPipe - msMax);
freeAll(&c);
return EXIT_SUCCESS;
}