code/day51-streams/streams.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 51: CUDA streams and overlap.
//
// Two kernels, four ways of launching them, one question per way: did the
// timeline overlap or not?
//
// spin-solo the small kernel alone. 8 blocks on a 40-SM card, so at
// least 32 SMs sit idle for its whole duration.
// sweep-solo the grid-filling kernel alone, four launches back to back
// in one stream.
// part 1 both in the legacy default stream. Sequential by
// definition: one stream is one ordered queue.
// part 2 small kernel in stream A, sweeps in stream B, both
// created with cudaStreamNonBlocking. This is the overlap.
// part 3 the trap. Small kernel in a stream made with plain
// cudaStreamCreate, sweeps in the legacy default stream.
// The legacy stream serialises against blocking streams,
// so the second stream buys nothing.
// part 4 two grid-filling kernels in two non-blocking streams.
// Independent, separate outputs, and still no overlap
// worth having: neither leaves the other any SMs.
//
// The timer here only corroborates. The proof is the nsys timeline, where
// part 2 shows sweepAdd bars starting and finishing inside the one spinLcg
// bar on a different stream row, and part 4 shows two lanes taking turns.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo -o streams streams.cu
// Run: ./streams
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cuda_runtime.h>
#include <nvtx3/nvToolsExt.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)
// 16 Mi elements is 64 MiB per buffer, three float buffers on the device.
// Large enough that one sweep is bandwidth-bound and visible on a timeline.
constexpr size_t kElems = 16777216;
// The correctness size. 611 is 13 x 47, so no block size divides it and the
// grid-stride loop's bounds condition does real work on every check launch.
constexpr size_t kCheckElems = 1048576 + 611;
constexpr int kThreadsPerBlock = 256; // 8 warps
// The grid-filling launch. The T4 runs 1024 resident threads per SM, so 4
// blocks of 256 fit each of the 40 SMs: 160 resident blocks fill the card.
// 320 keeps a full second wave queued, so the card stays full from the
// first cycle to nearly the last and a competing kernel finds no idle SM.
constexpr int kFullBlocks = 320;
// The deliberately small launch: 8 blocks touch at most 8 SMs and leave at
// least 32 idle. That spare capacity, not the second stream by itself, is
// what makes the overlap in part 2 possible.
constexpr int kSpinBlocks = 8;
constexpr size_t kSpinThreads =
static_cast<size_t>(kSpinBlocks) * kThreadsPerBlock;
// Each spin thread walks this many LCG steps. The steps form a dependent
// chain, so the kernel's duration is set by chain length, not thread count,
// and it can be tuned to dwarf one sweep without touching memory.
constexpr int kSpinIters = 1 << 20;
// The gate runs the CPU reference over every spin thread, so it uses a
// short chain. Correctness does not depend on the chain length; the launch
// and the arithmetic are the same.
constexpr int kCheckSpinIters = 4096;
// Four sweeps per lane, so the sweep lane's duration is comparable to the
// spin kernel's and the overlap in part 2 has something to hide.
constexpr int kSweepsPerLane = 4;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
static_assert(kThreadsPerBlock % 32 == 0,
"the block is a whole number of warps");
static_assert(kThreadsPerBlock <= 1024,
"1024 threads is the per-block maximum on every compute "
"capability this course targets");
static_assert(kCheckElems < kElems,
"the correctness pass reuses the timed buffers");
// One thread, one LCG chain, one 4-byte store at the end. The loop lives
// entirely in registers, so this kernel occupies SMs without touching
// memory bandwidth, which is what lets it share the card with a sweep.
//
// A warp's 32 final stores are 32 consecutive uints, 128 contiguous bytes.
// Launch assumption: gridDim.x * blockDim.x >= n. `iters` arrives as a
// runtime argument so the compiler cannot fold the chain at build time.
// snippet: kernels
__global__ void spinLcg(unsigned int* out, size_t n, int iters) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
unsigned int x = static_cast<unsigned int>(i);
for (int k = 0; k < iters; ++k) {
x = 1664525u * x + 1013904223u;
}
out[i] = x;
}
}
// Grid-stride sweep: out[i] = in[i] + 1.0f. Bandwidth-bound, and with
// kFullBlocks blocks it fills every SM. Consecutive threads read and write
// consecutive floats, so each warp moves 128-byte contiguous lines.
__global__ void sweepAdd(const float* in, float* out, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
out[i] = in[i] + 1.0f;
}
}
// end snippet
// Same shape as sweepAdd with a different name, so the two lanes of part 4
// are two distinct rows in `nsys stats` and cannot be confused.
__global__ void sweepScale(const float* in, float* out, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
out[i] = in[i] * 2.0f;
}
}
// CPU references. Written for obvious correctness; the caller owns every
// buffer. The inputs are small whole numbers, so + 1.0f and * 2.0f are
// exact in float and the comparisons below can demand equality.
static void spinLcgCpu(unsigned int* out, size_t n, int iters) {
for (size_t i = 0; i < n; ++i) {
unsigned int x = static_cast<unsigned int>(i);
for (int k = 0; k < iters; ++k) {
x = 1664525u * x + 1013904223u;
}
out[i] = x;
}
}
static void sweepAddCpu(const float* in, float* out, size_t n) {
for (size_t i = 0; i < n; ++i) {
out[i] = in[i] + 1.0f;
}
}
static void sweepScaleCpu(const float* in, float* out, size_t n) {
for (size_t i = 0; i < n; ++i) {
out[i] = in[i] * 2.0f;
}
}
// Returns the first index where got and want differ, or n. Exact compare:
// every operation both sides perform is exact on these inputs, so any
// difference is an index or launch bug, not rounding.
static size_t firstMismatchF(const float* got, const float* want, size_t n) {
for (size_t i = 0; i < n; ++i) {
if (got[i] != want[i]) {
return i;
}
}
return n;
}
static size_t firstMismatchU(const unsigned int* got, const unsigned int* want,
size_t n) {
for (size_t i = 0; i < n; ++i) {
if (got[i] != want[i]) {
return i;
}
}
return n;
}
// Times one closure with CUDA events recorded into the legacy default
// stream. Every closure below ends in cudaDeviceSynchronize, so the device
// is idle when `stop` is recorded: that synchronise, inside the closure, is
// where the async work is drained before the clock is read. `start` is
// recorded into an empty legacy queue, so it timestamps before any of the
// closure's launches reach the device. Warm-up runs the same closure, so
// every kernel a part times has paid its first-launch module load before
// the timed runs. Prints mean, min and max of the timed runs; returns the
// mean.
template <typename F>
static float timeRegion(const char* name, cudaEvent_t start, cudaEvent_t stop,
F&& launch) {
for (int r = 0; r < kWarmupRuns; ++r) {
launch();
}
float total = 0.0f;
float lo = 0.0f;
float hi = 0.0f;
for (int r = 0; r < kTimedRuns; ++r) {
CUDA_CHECK(cudaEventRecord(start));
launch();
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
total += ms;
lo = (r == 0 || ms < lo) ? ms : lo;
hi = (r == 0 || ms > hi) ? ms : hi;
}
CUDA_CHECK(cudaGetLastError());
const float mean = total / kTimedRuns;
std::printf("%-22s %9.3f ms (min %8.3f, max %8.3f)\n", name, mean, lo, hi);
return mean;
}
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("%d SMs, %d resident threads/SM, %d async copy engines\n",
prop.multiProcessorCount, prop.maxThreadsPerMultiProcessor,
prop.asyncEngineCount);
std::printf("full grid: %d blocks of %d; card holds %d resident\n",
kFullBlocks, kThreadsPerBlock,
prop.multiProcessorCount *
(prop.maxThreadsPerMultiProcessor / kThreadsPerBlock));
std::printf("small grid: %d blocks of %d\n\n", kSpinBlocks,
kThreadsPerBlock);
const size_t bytes = kElems * sizeof(float);
std::vector<float> h_in(kElems);
for (size_t i = 0; i < kElems; ++i) {
h_in[i] = static_cast<float>(i % 8192);
}
float* d_in = nullptr;
float* d_a = nullptr;
float* d_b = nullptr;
unsigned int* d_spin = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d_a, bytes));
CUDA_CHECK(cudaMalloc(&d_b, bytes));
CUDA_CHECK(cudaMalloc(&d_spin, kSpinThreads * sizeof(unsigned int)));
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), bytes, cudaMemcpyHostToDevice));
// Streams for the four parts. A and B are non-blocking: they do not
// synchronise with the legacy default stream, which is what part 2
// needs. C is made with plain cudaStreamCreate on purpose: it is a
// blocking stream, and part 3 exists to show what that costs.
cudaStream_t streamA = nullptr;
cudaStream_t streamB = nullptr;
cudaStream_t streamC = nullptr;
CUDA_CHECK(cudaStreamCreateWithFlags(&streamA, cudaStreamNonBlocking));
CUDA_CHECK(cudaStreamCreateWithFlags(&streamB, cudaStreamNonBlocking));
CUDA_CHECK(cudaStreamCreate(&streamC));
cudaEvent_t start = nullptr;
cudaEvent_t stop = nullptr;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
// Correctness gates, before any timing. Each kernel runs once in the
// legacy stream at the check size and must match its CPU reference
// exactly. A kernel that is fast and wrong has no business in a table.
int rc = EXIT_SUCCESS;
{
std::vector<float> h_got(kCheckElems);
std::vector<float> h_want(kCheckElems);
std::vector<unsigned int> h_gotU(kSpinThreads);
std::vector<unsigned int> h_wantU(kSpinThreads);
sweepAdd<<<kFullBlocks, kThreadsPerBlock>>>(d_in, d_a, kCheckElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_got.data(), d_a, kCheckElems * sizeof(float),
cudaMemcpyDeviceToHost));
sweepAddCpu(h_in.data(), h_want.data(), kCheckElems);
size_t bad = firstMismatchF(h_got.data(), h_want.data(), kCheckElems);
if (bad != kCheckElems) {
std::fprintf(stderr, "sweepAdd wrong at %zu: got %.9g want %.9g\n",
bad, h_got[bad], h_want[bad]);
rc = EXIT_FAILURE;
}
sweepScale<<<kFullBlocks, kThreadsPerBlock>>>(d_in, d_b, kCheckElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_got.data(), d_b, kCheckElems * sizeof(float),
cudaMemcpyDeviceToHost));
sweepScaleCpu(h_in.data(), h_want.data(), kCheckElems);
bad = firstMismatchF(h_got.data(), h_want.data(), kCheckElems);
if (bad != kCheckElems) {
std::fprintf(stderr,
"sweepScale wrong at %zu: got %.9g want %.9g\n", bad,
h_got[bad], h_want[bad]);
rc = EXIT_FAILURE;
}
spinLcg<<<kSpinBlocks, kThreadsPerBlock>>>(d_spin, kSpinThreads,
kCheckSpinIters);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_gotU.data(), d_spin,
kSpinThreads * sizeof(unsigned int),
cudaMemcpyDeviceToHost));
spinLcgCpu(h_wantU.data(), kSpinThreads, kCheckSpinIters);
bad = firstMismatchU(h_gotU.data(), h_wantU.data(), kSpinThreads);
if (bad != kSpinThreads) {
std::fprintf(stderr, "spinLcg wrong at %zu: got %u want %u\n", bad,
h_gotU[bad], h_wantU[bad]);
rc = EXIT_FAILURE;
}
}
if (rc == EXIT_SUCCESS) {
std::printf(
"correctness: all three kernels match the CPU "
"references\n\n");
// The launches below pass the stream as the fourth chevron
// argument. The third is dynamic shared memory (day 13's knob),
// zero here, and is written out only because the stream slot
// cannot be reached without it.
auto spinSolo = [&] {
nvtxRangePushA("spin-solo");
spinLcg<<<kSpinBlocks, kThreadsPerBlock>>>(d_spin, kSpinThreads,
kSpinIters);
CUDA_CHECK(cudaDeviceSynchronize());
nvtxRangePop();
};
auto sweepSolo = [&] {
nvtxRangePushA("sweep-solo");
for (int s = 0; s < kSweepsPerLane; ++s) {
sweepAdd<<<kFullBlocks, kThreadsPerBlock>>>(d_in, d_a, kElems);
}
CUDA_CHECK(cudaDeviceSynchronize());
nvtxRangePop();
};
auto scaleSolo = [&] {
nvtxRangePushA("scale-solo");
for (int s = 0; s < kSweepsPerLane; ++s) {
sweepScale<<<kFullBlocks, kThreadsPerBlock>>>(d_in, d_b,
kElems);
}
CUDA_CHECK(cudaDeviceSynchronize());
nvtxRangePop();
};
auto part1Legacy = [&] {
nvtxRangePushA("part1-legacy");
spinLcg<<<kSpinBlocks, kThreadsPerBlock>>>(d_spin, kSpinThreads,
kSpinIters);
for (int s = 0; s < kSweepsPerLane; ++s) {
sweepAdd<<<kFullBlocks, kThreadsPerBlock>>>(d_in, d_a, kElems);
}
CUDA_CHECK(cudaDeviceSynchronize());
nvtxRangePop();
};
// snippet: two-streams
auto part2TwoStreams = [&] {
nvtxRangePushA("part2-two-streams");
spinLcg<<<kSpinBlocks, kThreadsPerBlock, 0, streamA>>>(
d_spin, kSpinThreads, kSpinIters);
for (int s = 0; s < kSweepsPerLane; ++s) {
sweepAdd<<<kFullBlocks, kThreadsPerBlock, 0, streamB>>>(
d_in, d_a, kElems);
}
CUDA_CHECK(cudaDeviceSynchronize());
nvtxRangePop();
};
// end snippet
// snippet: legacy-trap
auto part3LegacyTrap = [&] {
nvtxRangePushA("part3-legacy-trap");
spinLcg<<<kSpinBlocks, kThreadsPerBlock, 0, streamC>>>(
d_spin, kSpinThreads, kSpinIters);
for (int s = 0; s < kSweepsPerLane; ++s) {
sweepAdd<<<kFullBlocks, kThreadsPerBlock>>>(d_in, d_a, kElems);
}
CUDA_CHECK(cudaDeviceSynchronize());
nvtxRangePop();
};
// end snippet
auto part4TwoFullGrids = [&] {
nvtxRangePushA("part4-two-full-grids");
for (int s = 0; s < kSweepsPerLane; ++s) {
sweepAdd<<<kFullBlocks, kThreadsPerBlock, 0, streamA>>>(
d_in, d_a, kElems);
}
for (int s = 0; s < kSweepsPerLane; ++s) {
sweepScale<<<kFullBlocks, kThreadsPerBlock, 0, streamB>>>(
d_in, d_b, kElems);
}
CUDA_CHECK(cudaDeviceSynchronize());
nvtxRangePop();
};
const float msSpin = timeRegion("spin-solo", start, stop, spinSolo);
const float msSweep = timeRegion("sweep-solo", start, stop, sweepSolo);
const float msScale = timeRegion("scale-solo", start, stop, scaleSolo);
const float msPart1 =
timeRegion("part1-legacy", start, stop, part1Legacy);
const float msPart2 =
timeRegion("part2-two-streams", start, stop, part2TwoStreams);
const float msPart3 =
timeRegion("part3-legacy-trap", start, stop, part3LegacyTrap);
const float msPart4 =
timeRegion("part4-two-full-grids", start, stop, part4TwoFullGrids);
std::printf("\n");
std::printf("part1 / part2 (overlap saving) %6.2fx\n",
msPart1 / msPart2);
std::printf("part2 / spin-solo (hiding quality) %6.2fx\n",
msPart2 / msSpin);
std::printf("part3 / part1 (trap vs sequential) %6.2fx\n",
msPart3 / msPart1);
std::printf("part4 / (sweep + scale solo) %6.2fx\n",
msPart4 / (msSweep + msScale));
}
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
CUDA_CHECK(cudaStreamDestroy(streamA));
CUDA_CHECK(cudaStreamDestroy(streamB));
CUDA_CHECK(cudaStreamDestroy(streamC));
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_spin));
return rc;
}