Parallel reductions

01 · The puzzle

One answer, a million inputs, no shared counter

Sum an array of a million floats. Sequentially it's a for-loop. On a GPU it's genuinely hard, and the difficulty is instructive: you have a million threads but they can't all add into one variable — a single shared counter serializes into a million-step traffic jam. The answer is a tree: pairwise-add neighbors, halving the active count each round — a million values become 500k partial sums, then 250k, then… until one value remains. That's log2(1M) = 20 rounds of work instead of a million sequential adds, with the rounds themselves fully parallel.

The GPU version runs the tree in place: each round, active threads read two elements apart by the current stride and write one of them. Day 07's coalescing keeps the reads efficient, and the array shrinks by half per round.

02 · The divergence trap

Same algorithm, one branch away from 2× slower

The textbook formulation puts the condition on the output side: if (i % (2*s) == 0) a[i] += a[i + s]. Look at it through day 05's lens: within a warp of 32 consecutive threads, the even lanes are active and the odd lanes idle — every warp diverges on every round, half the lanes doing nothing while masked. Flip the condition to the thread index and the picture changes completely:

// textbook: active lanes = evens → every warp diverges
if (i % (2*s) == 0) a[i] += a[i + s];

// stride form: active lanes = a contiguous prefix
//   s ≥ 32 → every warp is fully active or fully idle
if (i < s) a[i] += a[i + s];
s >>= 1;  // 512 → 256 → … → 1

With the stride form, the active set is always a contiguous prefix of the array. Above stride 32, each warp is entirely in or entirely out — no intra-warp divergence at all. Same algorithm, same result, dramatically different SIMT economics. This single rewrite is the cleanest demonstration in all of CUDA that branch structure is a performance decision, not just a logic decision.

03 · Warp shuffles

The last 32 steps need no shared memory at all

Below stride 32 the reduction lives inside a single warp, and there's a dedicated instruction for exactly this: __shfl_down_sync exchanges values between lanes of a warp through the register file, no shared memory, no barriers:

// full-warp tree reduction; result broadcast to lane 0
for (int off = 16; off > 0; off >>= 1)
    v += __shfl_down_sync(0xffffffff, v, off);

// off=16: lane i += lane i+16 → 16 partials in lanes 0..15
// off=8, 4, 2, 1: … → the warp total sits in lane 0

A complete fast reduction combines the pieces you now own: each block grid-strides over the array (day 14), accumulating a local sum in a register; the warp shuffles collapse each warp's 32 values to one; lane 0 of each warp parks its total in shared memory; the first warp shuffles those together; the block writes one partial sum. A second tiny kernel — or one atomic add — merges the per-block partials. Two levels of tree, and the only global traffic is the original reads plus one float per block.

04 · Why this exercise matters

Reductions are everywhere, wearing disguises

"Sum an array" sounds like a toy until you notice how much of machine learning is a reduction in a trench coat: every loss, every norm, every mean/variance in batchnorm and layer norm, every softmax's max-and-sum (day 21), every attention row's denominator (day 24). The map→tree-combine pattern and its implementation details — stride shape, shuffles, per-block partials — reappear verbatim inside all of them.

Day 28 revisits this exact kernel as an optimization exercise: the naive version moves ~20% of the 4070's bandwidth, and each principled fix climbs a rung — to ~88% — with the profiler (day 18) explaining every step. If you remember one chapter's code from the Writing CUDA phase, make it this one.

Next: day 17 — streams and concurrency: overlapping compute with data transfer instead of paying for them one after the other.