code/day22-divergence/divergence.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 22: warp divergence, measured.
//
// Part 1 is one kernel and four predicates. Every thread runs kIters fused
// multiply-adds down one of two paths, and under all four predicates
// exactly half of the grid's threads take each path, so the work is
// identical across rows and the only thing that changes is how the two
// halves are grouped. Two predicates split every warp down the middle. Two
// split the grid on a warp boundary, so every warp agrees with itself.
//
// Part 2 moves the same branch inside the loop, where each side is one
// instruction. That is the case the Best Practices Guide says the compiler
// may turn into predicated instructions rather than a branch, and a
// predicated region has no branch left to diverge on. A straight-line
// kernel with no branch runs first as the floor.
//
// Every row of both parts moves the same bytes: one float in and one float
// out per thread, with consecutive lanes on consecutive addresses. Global
// memory is a small constant term here, not the thing being measured.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o divergence divergence.cu
// Run: ./divergence
//
// Verified: run on a Tesla T4 (compute capability 7.5), driver 595.84,
// CUDA 12.6 (V12.6.85), on 2026-08-30. The only output that may be
// published as this program's output is the transcript in
// evidence/run-2026-08-30.txt. See research/REVIEW-PROCESS.md.
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>
#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)
// 32 on every GPU this course targets. The built-in `warpSize` is a
// run-time value, so it cannot size an array or appear in a static_assert;
// this is the compile-time copy.
constexpr unsigned int kWarpSize = 32u;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarpsPerBlock = kThreadsPerBlock / static_cast<int>(kWarpSize);
constexpr int kBlocks = 4096;
constexpr size_t kElems = static_cast<size_t>(kBlocks) * kThreadsPerBlock;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kRelTolerance = 1e-5f;
// kIters is large enough that the branch, not the launch or the two global
// accesses, is what the timer sees.
constexpr int kIters = 4096;
// The two paths. Both are one fused multiply-add per step, so a thread pays
// the same for either one and the aggregate work does not depend on the
// predicate. Every constant is a power of two or one plus a power of two,
// so it is exact in binary floating point.
//
// Both multipliers are above 1, which matters for the correctness check:
// the chain grows by a factor of about e over kIters steps on path A and
// about e squared on path B, so the answer depends on every iteration. A
// loop the compiler shortened, or a thread that took the other path, misses
// by far more than kRelTolerance.
constexpr float kInput = 1.5f;
constexpr float kMulA = 1.000244140625f; // 1 + 2^-12
constexpr float kAddA = 0.000244140625f; // 2^-12
constexpr float kMulB = 1.00048828125f; // 1 + 2^-11
constexpr float kAddB = 0.00048828125f; // 2^-11
// The four predicates the sweep runs, numbered from zero so a mode is also
// its row index, plus the constant one the no-branch kernel uses.
//
// kAlwaysA is deliberately not the next number in the run. The exercise adds
// a fifth predicate and raises kNumModes to 5, and a sentinel sitting at 4
// would be the value that predicate wants, so the sweep would run the
// no-branch case and mixNoBranch's reference would come from the reader's new
// predicate. Any value takesPathA does not name falls through to its default.
constexpr int kLaneParity = 0;
constexpr int kLaneHalf = 1;
constexpr int kWarpParity = 2;
constexpr int kBlockParity = 3;
constexpr int kNumModes = 4;
constexpr int kAlwaysA = -1;
// True when this thread takes path A.
//
// `mode` is a kernel argument, so it is the same value for every thread in
// the grid and the switch itself never diverges. `phase` is 0 for the whole
// of part 1 and the loop counter in part 2.
//
// kLaneParity and kLaneHalf both depend on the thread's position inside its
// warp, so the 32 lanes disagree and the warp splits. kWarpParity is the
// Best Practices Guide's own example of a condition that costs nothing: it
// depends only on threadIdx / warpSize, so it is constant across a warp.
// kBlockParity is constant across a whole block.
//
// __host__ __device__ so that the host reference below and the device
// kernels cannot drift apart, and constexpr so the counts underneath are
// checked by the compiler rather than claimed in a comment.
// snippet: predicate
__host__ __device__ constexpr bool takesPathA(int mode, unsigned int tid,
unsigned int block,
unsigned int phase) {
switch (mode) {
case kLaneParity:
return ((tid + phase) & 1u) == 0u;
case kLaneHalf:
return ((tid + phase) % kWarpSize) < (kWarpSize / 2u);
case kWarpParity:
return ((tid / kWarpSize + phase) & 1u) == 0u;
case kBlockParity:
return ((block + phase) & 1u) == 0u;
default:
return true;
}
}
// end snippet
// How many of a block's warps hold lanes that disagree, at phase 0. This is
// the number the whole lesson turns on, so the compiler computes it.
__host__ __device__ constexpr int splitWarpsPerBlock(int mode) {
int split = 0;
for (unsigned int w = 0; w < static_cast<unsigned int>(kWarpsPerBlock);
++w) {
const bool lane0 = takesPathA(mode, w * kWarpSize, 0u, 0u);
bool disagrees = false;
for (unsigned int lane = 1u; lane < kWarpSize; ++lane) {
if (takesPathA(mode, w * kWarpSize + lane, 0u, 0u) != lane0) {
disagrees = true;
}
}
if (disagrees) {
++split;
}
}
return split;
}
// How many threads of one block take path A, at phase 0.
__host__ __device__ constexpr int pathAThreads(int mode, unsigned int block) {
int count = 0;
for (unsigned int t = 0; t < static_cast<unsigned int>(kThreadsPerBlock);
++t) {
if (takesPathA(mode, t, block, 0u)) {
++count;
}
}
return count;
}
static_assert(kThreadsPerBlock % static_cast<int>(kWarpSize) == 0,
"block size must be a whole number of warps");
static_assert(kBlocks % 2 == 0,
"kBlockParity sends half the grid down each path only when the "
"grid holds an even number of blocks");
// The contrast the lesson rests on, checked at compile time. A predicate
// built from the lane splits all eight warps of the block whether it is
// written as a parity test or as a threshold; a predicate built from the
// warp index splits none of them.
// snippet: split-checks
static_assert(splitWarpsPerBlock(kLaneParity) == kWarpsPerBlock,
"threadIdx.x & 1 must split every warp in the block");
static_assert(splitWarpsPerBlock(kLaneHalf) == kWarpsPerBlock,
"a threshold inside the warp splits every warp too");
static_assert(splitWarpsPerBlock(kWarpParity) == 0,
"threadIdx.x / 32 is constant across a warp");
static_assert(splitWarpsPerBlock(kBlockParity) == 0,
"blockIdx.x is constant across a warp");
static_assert(splitWarpsPerBlock(kAlwaysA) == 0, "no branch, no split");
// end snippet
// And the fairness condition: the same number of threads take each path in
// every mode, so the four rows differ in grouping and in nothing else.
// snippet: fair-checks
static_assert(pathAThreads(kLaneParity, 0u) == kThreadsPerBlock / 2,
"half the block takes path A");
static_assert(pathAThreads(kLaneHalf, 0u) == kThreadsPerBlock / 2,
"half the block takes path A");
static_assert(pathAThreads(kWarpParity, 0u) == kThreadsPerBlock / 2,
"half the block takes path A. This one needs more than one warp "
"a block: at 32 threads threadIdx.x / 32 is 0 everywhere, so the "
"whole block takes path A and the warp-uniform rows stop being "
"comparable with the split ones");
static_assert(pathAThreads(kBlockParity, 0u) + pathAThreads(kBlockParity, 1u) ==
kThreadsPerBlock,
"an even and an odd block together send half their threads "
"down path A");
// end snippet
// One thread: reads one float, runs kIters fused multiply-adds down one of
// two paths, writes one float.
//
// One warp: lane L reads and writes element (base + L), so the 32 lanes
// cover 128 contiguous bytes at each end. Both accesses are coalesced and
// identical in every mode; the branch between them is the variable.
//
// Launch assumption: the grid covers n exactly, so the guard never fires
// here. It is written the way day 5 writes it because it is also the
// cheapest example of a branch nobody should worry about, and the page says
// why.
// snippet: one-branch
__global__ void mixOneBranch(const float* __restrict__ in,
float* __restrict__ out, size_t n, int mode) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
float x = (i < n) ? in[i] : 0.0f;
if (takesPathA(mode, threadIdx.x, blockIdx.x, 0u)) {
for (int k = 0; k < kIters; ++k) {
x = fmaf(x, kMulA, kAddA);
}
} else {
for (int k = 0; k < kIters; ++k) {
x = fmaf(x, kMulB, kAddB);
}
}
if (i < n) {
out[i] = x;
}
}
// end snippet
// The same arithmetic with the branch moved inside the loop, so each side
// is a single instruction and the predicate changes on every iteration.
//
// One thread: kIters iterations, one fused multiply-add each, the path
// alternating. One warp: the same two coalesced accesses as above.
//
// What this measures and what it does not: whether a short branch costs
// what a long one costs. It cannot tell you on its own whether the compiler
// emitted a branch or predicated the two sides; the README gives the
// cuobjdump line that answers that, and the timing table below implies it.
// snippet: per-iteration
__global__ void mixBranchPerIteration(const float* __restrict__ in,
float* __restrict__ out, size_t n,
int mode) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
float x = (i < n) ? in[i] : 0.0f;
for (int k = 0; k < kIters; ++k) {
if (takesPathA(mode, threadIdx.x, blockIdx.x,
static_cast<unsigned int>(k))) {
x = fmaf(x, kMulA, kAddA);
} else {
x = fmaf(x, kMulB, kAddB);
}
}
if (i < n) {
out[i] = x;
}
}
// end snippet
// The floor for part 2: the same chain with no branch anywhere in it. Every
// thread takes path A, which is what kAlwaysA means on the host side.
//
// One warp: the same two coalesced accesses as the other two kernels.
__global__ void mixNoBranch(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
float x = (i < n) ? in[i] : 0.0f;
for (int k = 0; k < kIters; ++k) {
x = fmaf(x, kMulA, kAddA);
}
if (i < n) {
out[i] = x;
}
}
// CPU reference. Written for obvious correctness, not speed: the same loop
// the kernel runs, in the same order, with the same fused multiply-add, so
// the two agree bit for bit and the tolerance below is slack rather than
// necessity.
//
// It fills 2 * kThreadsPerBlock values instead of kElems. Every predicate
// here depends on the thread's index inside its block and on whether the
// block is even or odd, so two blocks cover every thread in the grid.
static void expectedPerThread(int mode, bool perIteration, float input,
std::vector<float>& expected) {
for (unsigned int block = 0u; block < 2u; ++block) {
for (unsigned int t = 0u;
t < static_cast<unsigned int>(kThreadsPerBlock); ++t) {
float x = input;
for (int k = 0; k < kIters; ++k) {
const unsigned int phase =
perIteration ? static_cast<unsigned int>(k) : 0u;
if (takesPathA(mode, t, block, phase)) {
x = std::fmaf(x, kMulA, kAddA);
} else {
x = std::fmaf(x, kMulB, kAddB);
}
}
expected[block * kThreadsPerBlock + t] = x;
}
}
}
// Returns the first index where the device disagrees with the reference by
// more than the relative tolerance, or n if they agree everywhere.
// Returning the index rather than a bool is the point: "wrong at 33" names
// the lane, "wrong" does not.
static size_t firstMismatch(const float* got,
const std::vector<float>& expected, size_t n,
float relTolerance) {
for (size_t i = 0; i < n; ++i) {
const size_t tid = i % static_cast<size_t>(kThreadsPerBlock);
const size_t blockParity =
(i / static_cast<size_t>(kThreadsPerBlock)) & 1u;
const float want = expected[blockParity * kThreadsPerBlock + tid];
const float scale = (want == 0.0f) ? 1.0f : std::fabs(want);
if (std::fabs(got[i] - want) > relTolerance * scale) {
return i;
}
}
return n;
}
// Times a launch with CUDA events and returns the mean milliseconds per
// run. A host clock around a launch measures the launch, because launches
// are asynchronous. Day 9 takes that apart. Copy this helper 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);
std::printf("%d SMs, warp size %d, %d blocks x %d threads = %zu threads\n",
prop.multiProcessorCount, prop.warpSize, kBlocks,
kThreadsPerBlock, kElems);
const size_t bytes = kElems * sizeof(float);
std::printf("%d multiply-adds per thread, %zu bytes moved per launch\n",
kIters, 2 * bytes);
// Every failure below records itself and falls through to the one
// cleanup block at the bottom, so no path can return with device memory
// still allocated.
int failures = 0;
int rows = 0;
float* d_in = nullptr;
float* d_out = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
// One input value for the whole grid. The buffer exists so the compiler
// cannot see the value and fold the chain away; the check below is
// about which path each thread took and how many steps it ran, not
// about indexing, which is day 4's subject.
std::vector<float> h_in(kElems, kInput);
std::vector<float> h_out(kElems);
std::vector<float> expected(2 * kThreadsPerBlock);
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), bytes, cudaMemcpyHostToDevice));
// The four predicate constants are 0 to 3, so a mode is also its row
// index and branchMs[kWarpParity] is the warp-uniform row.
const char* labels[kNumModes] = {"threadIdx.x & 1",
"(threadIdx.x & 31) < 16",
"threadIdx.x / 32 & 1", "blockIdx.x & 1"};
float branchMs[kNumModes] = {0.0f, 0.0f, 0.0f, 0.0f};
for (int mode = 0; mode < kNumModes; ++mode) {
expectedPerThread(mode, false, kInput, expected);
// One untimed launch first, checked, so a wrong answer is reported
// before a number that came from it reaches the table.
mixOneBranch<<<kBlocks, kThreadsPerBlock>>>(d_in, d_out, kElems, mode);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
const size_t bad =
firstMismatch(h_out.data(), expected, kElems, kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr,
"mixOneBranch mode %d wrong at %zu: got %.9g\n", mode,
bad, h_out[bad]);
++failures;
}
branchMs[mode] = timeKernel([&] {
mixOneBranch<<<kBlocks, kThreadsPerBlock>>>(d_in, d_out, kElems,
mode);
});
++rows;
}
std::printf("\nPart 1: one branch, %d multiply-adds on each side\n",
kIters);
std::printf("%-24s %12s %12s %12s\n", "predicate", "warps split",
"time (ms)", "vs uniform");
std::printf("%-24s %12s %12s %12s\n", "------------------------",
"-----------", "----------", "----------");
for (int mode = 0; mode < kNumModes; ++mode) {
std::printf("%-24s %12d %12.3f %12.2f\n", labels[mode],
splitWarpsPerBlock(mode), branchMs[mode],
branchMs[mode] / branchMs[kWarpParity]);
}
std::printf(
"warps split is out of the %d warps in a block, computed "
"from the predicate\n",
kWarpsPerBlock);
// Part 2. The floor first, on the same buffers, in the same process.
expectedPerThread(kAlwaysA, false, kInput, expected);
mixNoBranch<<<kBlocks, kThreadsPerBlock>>>(d_in, d_out, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
size_t bad = firstMismatch(h_out.data(), expected, kElems, kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr, "mixNoBranch wrong at %zu: got %.9g\n", bad,
h_out[bad]);
++failures;
}
const float noBranchMs = timeKernel([&] {
mixNoBranch<<<kBlocks, kThreadsPerBlock>>>(d_in, d_out, kElems);
});
++rows;
const int shortModes[2] = {kLaneParity, kWarpParity};
const char* shortLabels[2] = {"per-iteration, lane parity",
"per-iteration, warp parity"};
float shortMs[2] = {0.0f, 0.0f};
for (int m = 0; m < 2; ++m) {
expectedPerThread(shortModes[m], true, kInput, expected);
mixBranchPerIteration<<<kBlocks, kThreadsPerBlock>>>(
d_in, d_out, kElems, shortModes[m]);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
bad = firstMismatch(h_out.data(), expected, kElems, kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr,
"mixBranchPerIteration mode %d wrong at %zu: "
"got %.9g\n",
shortModes[m], bad, h_out[bad]);
++failures;
}
shortMs[m] = timeKernel([&] {
mixBranchPerIteration<<<kBlocks, kThreadsPerBlock>>>(
d_in, d_out, kElems, shortModes[m]);
});
++rows;
}
std::printf(
"\nPart 2: the same branch inside the loop, one multiply-add "
"on each side\n");
std::printf("%-28s %12s %12s\n", "kernel", "time (ms)", "vs no branch");
std::printf("%-28s %12s %12s\n", "----------------------------",
"----------", "------------");
std::printf("%-28s %12.3f %12.2f\n", "no branch", noBranchMs, 1.0);
for (int m = 0; m < 2; ++m) {
std::printf("%-28s %12.3f %12.2f\n", shortLabels[m], shortMs[m],
shortMs[m] / noBranchMs);
}
std::printf(
"\nEvery row above moved the same %zu bytes and ran the same "
"%d\nmultiply-adds per thread. Only the branch changed.\n",
2 * bytes, kIters);
// The row count is checked so the lesson's tables cannot drift from what
// the program prints. A real branch, not an assert: CI builds Release,
// Release defines NDEBUG, and NDEBUG deletes assert(), so the check
// would be missing from exactly the build that matters.
const int kExpectedRows = kNumModes + 3;
if (rows != kExpectedRows) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's tables and "
"this program disagree\n",
rows, kExpectedRows);
++failures;
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_out));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}