code/day52-events/events_forkjoin.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 52: events and cross-stream dependencies.
//
// A fork-join diamond across five streams: one upstream kernel on stream 0,
// four independent branch kernels on streams 1 to 4, one join kernel back on
// stream 0. cudaEventRecord and cudaStreamWaitEvent are the only ordering
// between them. The program checks the diamond and a serial version against
// a CPU reference, times both, shows the documented error from timing with
// a cudaEventDisableTiming event, and then runs the diamond with one line
// moved (the record before the upstream launch) so the dependency silently
// vanishes and the branches race the kernel that feeds them.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo -o events_forkjoin \
// events_forkjoin.cu
// Run: ./events_forkjoin
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cuda_runtime.h>
#include <nvtx3/nvToolsExt.h>
#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)
// kGrindBlocks is 8 so the whole diamond fits the card: 8 upstream blocks
// plus 4 branches of 8 is 40, the T4's SM count. A branch kernel that filled
// the card would serialise against its siblings whatever the events say, and
// this lesson is about the events.
constexpr size_t kElems = 1u << 20; // 4 MiB per buffer
constexpr int kThreadsPerBlock = 256;
constexpr int kGrindBlocks = 8;
constexpr int kIters = 512; // fused multiply-adds per element
constexpr int kBranches = 4;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kA = 0.999f; // keeps x = a*x + b bounded near b/(1-a)
constexpr float kRelTolerance = 1e-6f;
// The perf gate needs the upstream and all four branches resident at once,
// 40 blocks. Below 24 SMs the branches queue for space and the ratio is
// reported, not gated.
constexpr int kMinSmsForGate = 24;
// x = fmaf(a, x, b), kIters times, per element. One thread grid-strides over
// the elements; a warp reads 32 consecutive floats, coalesced, four 32-byte
// sectors per load. No launch assumption: the grid is deliberately small
// (kGrindBlocks) and the loop covers the rest.
//
// fmaf keeps the arithmetic bit-identical to the CPU reference, which also
// uses fmaf: both sides round once per step, in the same order.
__global__ void grindFma(const float* __restrict__ in, float* __restrict__ out,
float a, float b, size_t n, int iters) {
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) {
float x = in[i];
for (int k = 0; k < iters; ++k) {
x = fmaf(a, x, b);
}
out[i] = x;
}
}
// out[i] = b0[i] + b1[i] + b2[i] + b3[i]. One thread owns one element; a
// warp's 32 addresses are consecutive in all five buffers. Launched with a
// full grid because it is cheap and nothing overlaps it.
__global__ void sumFour(const float* __restrict__ b0,
const float* __restrict__ b1,
const float* __restrict__ b2,
const float* __restrict__ b3, float* __restrict__ out,
size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
out[i] = b0[i] + b1[i] + b2[i] + b3[i];
}
}
// CPU reference. std::fmaf mirrors the kernel's fmaf step for step, and the
// four-way add in main's reference block keeps sumFour's order, so the
// tolerance below only has to absorb libm quirks, not reassociation.
static void grindCpu(const float* in, float* out, float a, float b, size_t n,
int iters) {
for (size_t i = 0; i < n; ++i) {
float x = in[i];
for (int k = 0; k < iters; ++k) {
x = std::fmaf(a, x, b);
}
out[i] = x;
}
}
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;
}
static size_t countMismatches(const float* got, const float* want, size_t n,
float relTolerance) {
size_t count = 0;
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) {
++count;
}
}
return count;
}
// Times a submission with CUDA events and returns the mean milliseconds per
// run. Copied verbatim from the course reference file; the lambda here
// enqueues a whole multi-stream diamond rather than one launch, which the
// timer does not care about: the stop event is recorded on the default
// stream, and every stream in this file is a blocking stream, so the stop
// waits for all 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;
}
// Everything the three submission functions share. Plain data, no methods.
struct Diamond {
cudaStream_t s[1 + kBranches]; // s[0] carries the upstream and the join
cudaEvent_t fork; // recorded on s[0] after the upstream
cudaEvent_t branchDone[kBranches];
const float* d_in;
float* d_mid;
float* d_branch[kBranches];
float* d_out;
};
// Per-branch fma offsets. Four different values so a branch reading another
// branch's output cannot pass by accident.
constexpr float kBranchB[kBranches] = {0.125f, 0.25f, 0.5f, 1.0f};
// snippet: fork-join
// The diamond. Every call here returns as soon as the work is queued; the
// ordering lives in the two event calls, not in the call order.
static void enqueueDiamond(const Diamond& d) {
nvtxRangePushA("upstream");
grindFma<<<kGrindBlocks, kThreadsPerBlock, 0, d.s[0]>>>(
d.d_in, d.d_mid, kA, 1.0f, kElems, kIters);
CUDA_CHECK(cudaEventRecord(d.fork, d.s[0]));
nvtxRangePop();
nvtxRangePushA("fork");
for (int j = 0; j < kBranches; ++j) {
CUDA_CHECK(cudaStreamWaitEvent(d.s[1 + j], d.fork, 0));
grindFma<<<kGrindBlocks, kThreadsPerBlock, 0, d.s[1 + j]>>>(
d.d_mid, d.d_branch[j], kA, kBranchB[j], kElems, kIters);
CUDA_CHECK(cudaEventRecord(d.branchDone[j], d.s[1 + j]));
}
nvtxRangePop();
nvtxRangePushA("join");
for (int j = 0; j < kBranches; ++j) {
CUDA_CHECK(cudaStreamWaitEvent(d.s[0], d.branchDone[j], 0));
}
const int blocks =
static_cast<int>((kElems + kThreadsPerBlock - 1) / kThreadsPerBlock);
sumFour<<<blocks, kThreadsPerBlock, 0, d.s[0]>>>(
d.d_branch[0], d.d_branch[1], d.d_branch[2], d.d_branch[3], d.d_out,
kElems);
nvtxRangePop();
}
// end snippet
// The same nine kernels with no events: everything on s[0], back to back.
static void enqueueSerial(const Diamond& d) {
nvtxRangePushA("serial");
grindFma<<<kGrindBlocks, kThreadsPerBlock, 0, d.s[0]>>>(
d.d_in, d.d_mid, kA, 1.0f, kElems, kIters);
for (int j = 0; j < kBranches; ++j) {
grindFma<<<kGrindBlocks, kThreadsPerBlock, 0, d.s[0]>>>(
d.d_mid, d.d_branch[j], kA, kBranchB[j], kElems, kIters);
}
const int blocks =
static_cast<int>((kElems + kThreadsPerBlock - 1) / kThreadsPerBlock);
sumFour<<<blocks, kThreadsPerBlock, 0, d.s[0]>>>(
d.d_branch[0], d.d_branch[1], d.d_branch[2], d.d_branch[3], d.d_out,
kElems);
nvtxRangePop();
}
// snippet: broken-order
// The trap: one line moved. The record now runs before the upstream launch,
// so the event captures s[0] while it is empty. Every cudaStreamWaitEvent
// below is honoured, immediately, against that empty snapshot, and the
// branches race the kernel that feeds them. No API returns an error.
static void enqueueBroken(const Diamond& d, cudaEvent_t forkTooEarly) {
nvtxRangePushA("broken-order");
CUDA_CHECK(cudaEventRecord(forkTooEarly, d.s[0])); // captures nothing
grindFma<<<kGrindBlocks, kThreadsPerBlock, 0, d.s[0]>>>(
d.d_in, d.d_mid, kA, 1.0f, kElems, kIters);
for (int j = 0; j < kBranches; ++j) {
CUDA_CHECK(cudaStreamWaitEvent(d.s[1 + j], forkTooEarly, 0));
grindFma<<<kGrindBlocks, kThreadsPerBlock, 0, d.s[1 + j]>>>(
d.d_mid, d.d_branch[j], kA, kBranchB[j], kElems, kIters);
CUDA_CHECK(cudaEventRecord(d.branchDone[j], d.s[1 + j]));
}
for (int j = 0; j < kBranches; ++j) {
CUDA_CHECK(cudaStreamWaitEvent(d.s[0], d.branchDone[j], 0));
}
const int blocks =
static_cast<int>((kElems + kThreadsPerBlock - 1) / kThreadsPerBlock);
sumFour<<<blocks, kThreadsPerBlock, 0, d.s[0]>>>(
d.d_branch[0], d.d_branch[1], d.d_branch[2], d.d_branch[3], d.d_out,
kElems);
nvtxRangePop();
}
// end snippet
static void syncAll(const Diamond& d) {
for (int j = 0; j <= kBranches; ++j) {
CUDA_CHECK(cudaStreamSynchronize(d.s[j]));
}
CUDA_CHECK(cudaGetLastError());
}
int main() {
CUDA_CHECK(cudaSetDevice(0));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
std::printf("GPU: %s (compute capability %d.%d, %d SMs)\n", prop.name,
prop.major, prop.minor, prop.multiProcessorCount);
std::printf("elems: %zu, grind blocks: %d, iters: %d, branches: %d\n",
kElems, kGrindBlocks, kIters, kBranches);
const size_t bytes = kElems * sizeof(float);
std::vector<float> h_in(kElems);
std::vector<float> h_out(kElems);
std::vector<float> h_want(kElems);
for (size_t i = 0; i < kElems; ++i) {
h_in[i] = static_cast<float>(i % 251) * 0.01f;
}
// Reference: mid = grind(in, b=1), branch j = grind(mid, kBranchB[j]),
// out = sum of the four branches, all in the kernels' own order.
{
std::vector<float> h_mid(kElems);
std::vector<float> h_branch(kElems);
grindCpu(h_in.data(), h_mid.data(), kA, 1.0f, kElems, kIters);
for (int j = 0; j < kBranches; ++j) {
grindCpu(h_mid.data(), h_branch.data(), kA, kBranchB[j], kElems,
kIters);
for (size_t i = 0; i < kElems; ++i) {
h_want[i] = (j == 0) ? h_branch[i] : h_want[i] + h_branch[i];
}
}
}
Diamond d = {};
float* d_in = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d.d_mid, bytes));
for (int j = 0; j < kBranches; ++j) {
CUDA_CHECK(cudaMalloc(&d.d_branch[j], bytes));
}
CUDA_CHECK(cudaMalloc(&d.d_out, bytes));
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), bytes, cudaMemcpyHostToDevice));
d.d_in = d_in;
for (int j = 0; j <= kBranches; ++j) {
CUDA_CHECK(cudaStreamCreate(&d.s[j]));
}
// The dependency events carry no timestamp. cudaEventDisableTiming is
// the documented fast path for cudaStreamWaitEvent; the timing events in
// timeKernel keep the default flags because they are read with
// cudaEventElapsedTime.
CUDA_CHECK(cudaEventCreateWithFlags(&d.fork, cudaEventDisableTiming));
for (int j = 0; j < kBranches; ++j) {
CUDA_CHECK(
cudaEventCreateWithFlags(&d.branchDone[j], cudaEventDisableTiming));
}
cudaEvent_t forkTooEarly;
CUDA_CHECK(cudaEventCreateWithFlags(&forkTooEarly, cudaEventDisableTiming));
// Gate failures set status and fall through to the cleanup at the
// bottom, so no return path leaks a buffer, stream or event.
int status = EXIT_SUCCESS;
// 1. Correctness of the diamond and the serial version.
enqueueDiamond(d);
syncAll(d);
CUDA_CHECK(
cudaMemcpy(h_out.data(), d.d_out, bytes, cudaMemcpyDeviceToHost));
size_t bad =
firstMismatch(h_out.data(), h_want.data(), kElems, kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr, "diamond wrong at %zu: got %.9g, want %.9g\n", bad,
h_out[bad], h_want[bad]);
status = EXIT_FAILURE;
} else {
std::printf("diamond: PASS (0 mismatches)\n");
}
if (status == EXIT_SUCCESS) {
enqueueSerial(d);
syncAll(d);
CUDA_CHECK(
cudaMemcpy(h_out.data(), d.d_out, bytes, cudaMemcpyDeviceToHost));
bad = firstMismatch(h_out.data(), h_want.data(), kElems, kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr, "serial wrong at %zu: got %.9g, want %.9g\n",
bad, h_out[bad], h_want[bad]);
status = EXIT_FAILURE;
} else {
std::printf("serial: PASS (0 mismatches)\n");
}
}
// 2. A dependency event cannot be read as a stopwatch. Documented: "If
// either event was created with the cudaEventDisableTiming flag, then
// this function will return cudaErrorInvalidResourceHandle."
if (status == EXIT_SUCCESS) {
float bogus = 0.0f;
const cudaError_t timingErr =
cudaEventElapsedTime(&bogus, d.fork, d.branchDone[0]);
if (timingErr != cudaErrorInvalidResourceHandle) {
std::fprintf(stderr,
"expected cudaErrorInvalidResourceHandle from "
"cudaEventElapsedTime, got %s\n",
cudaGetErrorName(timingErr));
status = EXIT_FAILURE;
} else {
std::printf("cudaEventElapsedTime on the fork event: %s\n",
cudaGetErrorName(timingErr));
}
// The failed call above set the last-error slot on purpose; clear it
// so the next checked call does not inherit it.
(void)cudaGetLastError();
}
// 3. Time both shapes. Same nine launches, only the ordering differs.
if (status == EXIT_SUCCESS) {
const float msDiamond = timeKernel([&] { enqueueDiamond(d); });
const float msSerial = timeKernel([&] { enqueueSerial(d); });
const float ratio = msSerial / msDiamond;
std::printf("diamond: %.3f ms, serial: %.3f ms, ratio %.2fx\n",
msDiamond, msSerial, ratio);
if (prop.multiProcessorCount >= kMinSmsForGate) {
if (ratio < 1.3f) {
std::fprintf(stderr,
"fork-join beats serial by only %.2fx on %d "
"SMs; the branches are not overlapping\n",
ratio, prop.multiProcessorCount);
status = EXIT_FAILURE;
}
} else {
std::printf("under %d SMs: ratio reported, not gated\n",
kMinSmsForGate);
}
}
// 4. The broken ordering. d_mid is zeroed first so a branch that runs
// too early reads zeros instead of the upstream result. This is a real
// data race, so the count below is an observation, not a gate: on some
// run, somewhere, it will be zero, and the program is still wrong.
if (status == EXIT_SUCCESS) {
CUDA_CHECK(cudaMemset(d.d_mid, 0, bytes));
enqueueBroken(d, forkTooEarly);
syncAll(d);
CUDA_CHECK(
cudaMemcpy(h_out.data(), d.d_out, bytes, cudaMemcpyDeviceToHost));
const size_t raced =
countMismatches(h_out.data(), h_want.data(), kElems, kRelTolerance);
std::printf(
"broken order: %zu of %zu elements wrong (race: the count "
"is not stable)\n",
raced, kElems);
}
CUDA_CHECK(cudaEventDestroy(forkTooEarly));
CUDA_CHECK(cudaEventDestroy(d.fork));
for (int j = 0; j < kBranches; ++j) {
CUDA_CHECK(cudaEventDestroy(d.branchDone[j]));
}
for (int j = 0; j <= kBranches; ++j) {
CUDA_CHECK(cudaStreamDestroy(d.s[j]));
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d.d_mid));
for (int j = 0; j < kBranches; ++j) {
CUDA_CHECK(cudaFree(d.d_branch[j]));
}
CUDA_CHECK(cudaFree(d.d_out));
return status;
}