code/day24-reduction-1/reduction.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 24: parallel reduction, versions 1 to 4.
//
// Four block-level sum kernels, one change between each neighbouring pair,
// timed on the same input in the same process:
//
// 1 reduceInterleavedDivergent stride doubles, `tid % (2 * s) == 0`
// 2 reduceInterleavedStrided same threads packed, strided shared index
// 3 reduceSequentialAddressing same threads, contiguous shared index
// 4 reduceSequentialTwoLoads version 3, plus one add during the load
//
// Every row reduces the same kElems floats down to a single float. That is
// the discipline, and it matters most for version 4: it launches half as many
// blocks, so it also leaves half as many partial sums for the next pass, and
// timing only the first launch would credit it for work it moved rather than
// removed. Day 11 is the lesson about that class of benchmark bug.
//
// Every element is 1.0f and kElems is below 2^24, so every partial sum in
// every version is a whole number that float represents exactly. Reduction
// order cannot change the answer on this input, which is why the check below
// is an exact comparison with no tolerance. Day 68 is the lesson where the
// order does change it.
//
// What this program does not check: that any version is faster than any
// other. That is the claim the run exists to test, so encoding it as a gate
// would make the test unfalsifiable.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o reduction reduction.cu
// Run: ./reduction
//
// 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 <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, and 32 shared memory banks to match.
// main() reads the device's own warpSize back and fails if it disagrees,
// because every count asserted below assumes this value.
constexpr int kWarpSize = 32;
constexpr int kBanks = 32;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr int kVersions = 4;
// 2^24 minus 611. Two properties, both load bearing. The sum is 16,776,605,
// which is under 2^24, so every partial sum in every version lands on a float
// exactly. And 611 is 13 x 47, so no block size divides the element count and
// the bounds check in every kernel runs on every launch instead of never.
constexpr size_t kElems = 16777216ull - 611ull;
// gcd(stride, kBanks) by Euclid, constexpr so the static_asserts can call it.
// Lanes L and L' of a warp share a bank when (L - L') * stride is a multiple
// of 32, so a warp reading word L * stride covers 32 / gcd(stride, 32)
// distinct banks. Day 15 measured this map on the same card.
constexpr int gcdWithBanks(int stride) {
int a = stride % kBanks;
int b = kBanks;
while (a != 0) {
const int t = b % a;
b = a;
a = t;
}
return b;
}
// How many threads the tree leaves active at step `s`, for one version.
//
// Version 1 activates the multiples of 2s and versions 2, 3 and 4 activate
// the first kThreadsPerBlock / (2s) threads. Both sets have the same size at
// every step. The versions differ in where those threads sit, not how many
// there are, which is what the two counters below take apart.
constexpr int activeThreads(int s) {
return kThreadsPerBlock / (2 * s);
}
// Warp-wide executions of the tree body in one block, summed over the whole
// tree. A warp runs the body whenever it holds at least one active thread,
// masking off the lanes that are not, so this is the count that divergence
// changes and the instruction count follows.
//
// Version 1 spreads its active threads 2s apart across the whole block, so
// every warp holds one until 2s passes the warp size. Versions 2, 3 and 4
// pack theirs into the low threads, so whole warps retire early.
constexpr int warpBodyRuns(int version) {
int total = 0;
for (int s = 1; s < kThreadsPerBlock; s *= 2) {
if (version == 1) {
const int span = 2 * s;
total += (span <= kWarpSize) ? kThreadsPerBlock / kWarpSize
: kThreadsPerBlock / span;
} else {
total += (activeThreads(s) + kWarpSize - 1) / kWarpSize;
}
}
return total;
}
// The same count, with each body run multiplied by how many ways its shared
// access splits on a bank conflict. This is the number to compare between
// versions, because a conflicted access costs the shared memory pipe one
// extra pass per way.
//
// Version 1 never conflicts: its active lanes read their own thread ids,
// which are distinct modulo 32 inside a warp. Versions 3 and 4 never conflict
// either: their active lanes read words 0 to s-1, also distinct. Version 2
// reads word 2s * lane, which puts gcd(2s, 32) lanes in every bank it uses.
constexpr int serializedBodyRuns(int version) {
if (version != 2) {
return warpBodyRuns(version);
}
int total = 0;
for (int s = 1; s < kThreadsPerBlock; s *= 2) {
const int active = activeThreads(s);
const int warps = (active + kWarpSize - 1) / kWarpSize;
const int lanes = (active < kWarpSize) ? active : kWarpSize;
const int banks = kBanks / gcdWithBanks(2 * s);
const int distinct = (lanes < banks) ? lanes : banks;
total += warps * (lanes / distinct);
}
return total;
}
// Additions one block's tree performs, which is one fewer than the elements
// it started with, in every version.
constexpr int treeAdds(int threads) {
int total = 0;
for (int s = threads / 2; s > 0; s /= 2) {
total += s;
}
return total;
}
// Elements one block consumes. Version 4 takes two per thread instead of one,
// which is the whole of the change and the reason its grid is half the size.
constexpr int elemsPerBlock(int version) {
return (version == 4) ? 2 * kThreadsPerBlock : kThreadsPerBlock;
}
static_assert(kThreadsPerBlock % kWarpSize == 0,
"block size must be a whole number of warps");
static_assert((kThreadsPerBlock & (kThreadsPerBlock - 1)) == 0,
"every tree here halves the stride, so the block size must be "
"a power of two");
static_assert(treeAdds(kThreadsPerBlock) == kThreadsPerBlock - 1,
"a tree reduction of n elements performs n - 1 additions, "
"however the steps are arranged");
// The model the timing table is checked against, settled at compile time.
//
// Read the four together and the ladder stops being four unrelated tricks.
// Version 2 cuts the body runs by a factor of four and hands all of it back
// in bank conflicts. Version 3 keeps version 2's thread set and drops the
// conflicts, so it is the first rung that removes anything. Version 4 keeps
// version 3's tree and halves how often the block pays for it.
// The two literal counts are written as implications rather than as plain
// equalities, so changing kThreadsPerBlock still builds. Day 24's exercise
// asks for exactly that, and an assert that forbids the exercise is worse
// than no assert at all. The three claims below the pair hold at every block
// size, which is the stronger statement and the one the lesson rests on.
static_assert(kThreadsPerBlock != 256 || warpBodyRuns(1) == 47,
"at 256 threads version 1 runs the tree body 47 times a block");
static_assert(kThreadsPerBlock != 256 || warpBodyRuns(2) == 12,
"at 256 threads packing the active threads cuts that to 12");
static_assert(warpBodyRuns(1) > warpBodyRuns(3) || kThreadsPerBlock <= 32,
"above one warp a block, version 1 always runs the body more "
"often than version 3");
static_assert(warpBodyRuns(2) == warpBodyRuns(3),
"versions 2 and 3 activate exactly the same threads at every "
"step, so the only difference between them is the address");
static_assert(serializedBodyRuns(1) == serializedBodyRuns(2),
"at this block size version 2 trades divergence for bank "
"conflicts one for one, which is the prediction the run judges");
static_assert(serializedBodyRuns(3) < serializedBodyRuns(2) ||
kThreadsPerBlock <= 32,
"above one warp a block, sequential addressing is the rung the "
"model says is worth the most. At 32 threads the block is one "
"warp, so there is no divergence to remove and no conflict to "
"avoid, and versions 2 and 3 do the same work");
static_assert(kThreadsPerBlock != 256 || serializedBodyRuns(2) == 47,
"at 256 threads the conflicts cost version 2 all 47 back");
// Version 1. Interleaved addressing with a divergent branch.
//
// One thread: loads one element into the tile, then at every step either adds
// its neighbour s away or sits masked.
//
// One warp: the load is 32 consecutive floats, 128 contiguous bytes, four
// 32-byte sectors. The shared reads are words `tid` and `tid + s`, which are
// distinct modulo 32 across the warp, so nothing conflicts. What costs here
// is that the active threads are the multiples of 2s, spread across every
// warp, so no warp is ever finished early.
//
// Launch assumption: exactly kThreadsPerBlock threads per block. The tile is
// sized from that constant and the loop counts up to it.
__global__ void reduceInterleavedDivergent(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const unsigned int width = kThreadsPerBlock;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
// The guard covers the load, never the barrier. Every thread in the block
// reaches every __syncthreads() below, including the threads whose
// element is past the end of the array.
tile[tid] = (i < n) ? in[i] : 0.0f;
__syncthreads();
// snippet: tree-divergent
for (unsigned int s = 1; s < width; s *= 2) {
if (tid % (2 * s) == 0) {
tile[tid] += tile[tid + s];
}
__syncthreads();
}
// end snippet
if (tid == 0) {
out[blockIdx.x] = tile[0];
}
}
// Version 2. Interleaved addressing, divergence removed.
//
// One change from version 1: the branch is on a packed index instead of on
// the thread id, so the active threads of a step are the first
// kThreadsPerBlock / (2s) of them and every warp above that retires.
//
// One warp: the global load is unchanged. The shared reads are now words
// 2s * lane and 2s * lane + s. Lane and lane + 32 / gcd(2s, 32) land in the
// same bank holding different words, so the access serializes that many ways,
// and by step 16 every active lane is in bank 0.
//
// Launch assumption: exactly kThreadsPerBlock threads per block.
__global__ void reduceInterleavedStrided(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const unsigned int width = kThreadsPerBlock;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
tile[tid] = (i < n) ? in[i] : 0.0f;
__syncthreads();
// snippet: tree-strided
for (unsigned int s = 1; s < width; s *= 2) {
const unsigned int index = 2 * s * tid;
if (index < width) {
tile[index] += tile[index + s];
}
__syncthreads();
}
// end snippet
if (tid == 0) {
out[blockIdx.x] = tile[0];
}
}
// Version 3. Sequential addressing.
//
// One change from version 2: the stride halves instead of doubling, so an
// active thread reads words `tid` and `tid + s`. The same threads are active
// at every step as in version 2, in the same numbers. Only the address moves.
//
// One warp: the global load is unchanged. The active lanes now read 32
// consecutive shared words, which is one bank each, so nothing serializes.
//
// Launch assumption: exactly kThreadsPerBlock threads per block.
__global__ void reduceSequentialAddressing(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const unsigned int width = kThreadsPerBlock;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
tile[tid] = (i < n) ? in[i] : 0.0f;
__syncthreads();
// snippet: tree-sequential
for (unsigned int s = width / 2; s > 0; s >>= 1) {
if (tid < s) {
tile[tid] += tile[tid + s];
}
__syncthreads();
}
// end snippet
if (tid == 0) {
out[blockIdx.x] = tile[0];
}
}
// Version 4. Sequential addressing, first add during the load.
//
// One change from version 3: the block covers 2 * blockDim.x elements and the
// thread adds its two before the tree starts, so the grid is half the size
// and half of the block's threads were never idle in the first step at all.
//
// One warp: two global loads, each 32 consecutive floats and 128 contiguous
// bytes, blockDim.x elements apart. They are independent, so both are in
// flight before the barrier, which is a second thing this change buys and the
// reason it is not only about the block count.
//
// The index is written as blockIdx.x * blockDim.x * 2 and not as
// 2 * (blockIdx.x * blockDim.x + tid). The second form is a partition of
// nothing: it gives a thread elements 2t and 2t + blockDim.x, so the odd
// elements are never read and a band of even ones is read twice.
//
// Launch assumption: exactly kThreadsPerBlock threads per block, and a grid
// covering ceil(n / (2 * kThreadsPerBlock)) blocks.
__global__ void reduceSequentialTwoLoads(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const unsigned int width = kThreadsPerBlock;
// snippet: two-loads
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) * 2 + tid;
float sum = (i < n) ? in[i] : 0.0f;
if (i + blockDim.x < n) {
sum += in[i + blockDim.x];
}
tile[tid] = sum;
__syncthreads();
// end snippet
for (unsigned int s = width / 2; s > 0; s >>= 1) {
if (tid < s) {
tile[tid] += tile[tid + s];
}
__syncthreads();
}
if (tid == 0) {
out[blockIdx.x] = tile[0];
}
}
static void launchVersion(int version, int blocks, const float* d_src,
float* d_dst, size_t m) {
switch (version) {
case 1:
reduceInterleavedDivergent<<<blocks, kThreadsPerBlock>>>(d_src,
d_dst, m);
break;
case 2:
reduceInterleavedStrided<<<blocks, kThreadsPerBlock>>>(d_src, d_dst,
m);
break;
case 3:
reduceSequentialAddressing<<<blocks, kThreadsPerBlock>>>(d_src,
d_dst, m);
break;
default:
reduceSequentialTwoLoads<<<blocks, kThreadsPerBlock>>>(d_src, d_dst,
m);
break;
}
}
// Runs one complete reduction, n floats down to one, and returns the buffer
// holding the answer. The two scratch buffers alternate, so no launch ever
// reads and writes the same allocation.
//
// Launches only. There is no synchronize and no error check inside, because
// this function runs inside the timed region and a synchronize there would
// measure something other than the kernels. main() checks the launches on the
// untimed call it makes first.
static const float* reduceAll(int version, const float* d_in, size_t n,
float* d_a, float* d_b, int* passes) {
const float* src = d_in;
float* dst = d_a;
size_t m = n;
int count = 0;
while (m > 1) {
const size_t per = static_cast<size_t>(elemsPerBlock(version));
const int blocks = static_cast<int>((m + per - 1) / per);
launchVersion(version, blocks, src, dst, m);
++count;
m = static_cast<size_t>(blocks);
src = dst;
dst = (dst == d_a) ? d_b : d_a;
}
*passes = count;
return src;
}
// Bytes this version reads across every pass. The passes after the first are
// under half a percent of the total, and how far under depends on the
// version, so the program counts them rather than quoting the input size and
// hoping. Writes come to the same fraction again and are not counted; the
// page says so next to the table.
static double bytesRead(int version) {
double total = 0.0;
size_t m = kElems;
while (m > 1) {
total += static_cast<double>(m) * sizeof(float);
const size_t per = static_cast<size_t>(elemsPerBlock(version));
m = (m + per - 1) / per;
}
return total;
}
// 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), %d SMs, warp size %d\n",
prop.name, prop.major, prop.minor, prop.multiProcessorCount,
prop.warpSize);
std::printf("Shared memory per block: %zu bytes\n", prop.sharedMemPerBlock);
// 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;
if (prop.warpSize != kWarpSize) {
std::fprintf(stderr,
"this card reports warp size %d; every count in this "
"program assumes %d\n",
prop.warpSize, kWarpSize);
++failures;
}
const char* names[kVersions] = {
"1 interleaved, divergent", "2 interleaved, strided",
"3 sequential addressing", "4 sequential, two loads"};
std::printf("\nModel per block, computed at compile time\n");
std::printf("%-28s %14s %12s\n", "version", "warp runs", "serialised");
std::printf("%-28s %14s %12s\n", "----------------------------",
"--------------", "------------");
for (int v = 1; v <= kVersions; ++v) {
std::printf("%-28s %14d %12d\n", names[v - 1], warpBodyRuns(v),
serializedBodyRuns(v));
}
// The input is 1.0f everywhere, so the answer is the element count and a
// wrong answer means a wrong number of elements was summed, which is
// exactly the failure the two-loads index invites. A sum is blind to
// order by construction, so no input can catch a kernel that reads the
// right elements in a different sequence, and none is claimed to.
const size_t inputBytes = kElems * sizeof(float);
const size_t scratchElems =
(kElems + kThreadsPerBlock - 1) / kThreadsPerBlock;
const size_t scratchBytes = scratchElems * sizeof(float);
const float want = static_cast<float>(kElems);
float* d_in = nullptr;
float* d_a = nullptr;
float* d_b = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, inputBytes));
CUDA_CHECK(cudaMalloc(&d_a, scratchBytes));
CUDA_CHECK(cudaMalloc(&d_b, scratchBytes));
std::vector<float> h_in(kElems, 1.0f);
CUDA_CHECK(
cudaMemcpy(d_in, h_in.data(), inputBytes, cudaMemcpyHostToDevice));
std::printf("\n%zu elements, %.0f MiB in, every element 1.0f\n", kElems,
static_cast<double>(inputBytes) / (1024.0 * 1024.0));
std::printf("%-28s %7s %12s %12s %9s\n", "version", "passes", "time (ms)",
"GB/s", "vs v1");
std::printf("%-28s %7s %12s %12s %9s\n", "----------------------------",
"------", "----------", "----------", "--------");
float versionMs[kVersions] = {0.0f, 0.0f, 0.0f, 0.0f};
int rows = 0;
for (int v = 1; v <= kVersions; ++v) {
// One untimed reduction first, checked, so a wrong answer is reported
// before a number that came from it reaches the table.
int passes = 0;
const float* d_result = reduceAll(v, d_in, kElems, d_a, d_b, &passes);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
float got = 0.0f;
CUDA_CHECK(
cudaMemcpy(&got, d_result, sizeof(float), cudaMemcpyDeviceToHost));
if (got != want) {
std::fprintf(stderr,
"version %d summed to %.9g, want %.9g; that is a "
"difference of %.9g elements\n",
v, got, want, want - got);
++failures;
}
versionMs[v - 1] = timeKernel([&] {
int timedPasses = 0;
reduceAll(v, d_in, kElems, d_a, d_b, &timedPasses);
});
++rows;
const double gbps = bytesRead(v) /
(static_cast<double>(versionMs[v - 1]) * 1.0e-3) /
1.0e9;
std::printf("%-28s %7d %12.3f %12.1f %9.2f\n", names[v - 1], passes,
versionMs[v - 1], gbps, versionMs[0] / versionMs[v - 1]);
}
// The step speedups, which are what the lesson's ladder is about. Printed
// rather than left to the reader, so the page and the program cannot
// disagree about which rung was worth the most.
std::printf("\nstep speedups\n");
int best = 1;
for (int v = 2; v <= kVersions; ++v) {
const double step = static_cast<double>(versionMs[v - 2]) /
static_cast<double>(versionMs[v - 1]);
std::printf(" v%d -> v%d %.2fx\n", v - 1, v, step);
const double bestStep = static_cast<double>(versionMs[best - 1]) /
static_cast<double>(versionMs[best]);
if (step > bestStep) {
best = v - 1;
}
}
std::printf("largest single step: v%d -> v%d\n", best, best + 1);
// The row count is checked so the lesson's table 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.
if (rows != kVersions) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's table and "
"this program disagree\n",
rows, kVersions);
++failures;
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}