Tree Reductions and Floating-Point Summation Order

Summing an array on a parallel machine is where the difference between a computation that is correct and one that is defined by a specification becomes visible. Floating-point addition is not associative: $a + (b + c)$ and $(a + b) + c$ can differ in the last bit, and over a million terms those differences accumulate. A tree reduction changes the order of every addition, so it produces a different result from a sequential loop over the same data. Both facts hold at once: the tree is faster, and it produces a different number. The design question is which result the caller needs.

The dependency chain is the cost

A sequential sum over $n$ terms is $n$ dependent additions. On a CPU, each add has a latency of a few cycles and the loop runs at one add per latency period. On a GPU the chain is worse in relative terms because the machine’s whole purpose is to issue many independent operations per cycle, and a chain of dependent adds exposes none of that parallelism.

The tree form breaks the chain. Summing pairs, then pairs of pairs, gives depth $\lceil \log_2 n \rceil$, which is 20 for a million terms, and each level is fully parallel. A block of 1024 threads reduces in ten levels: five shuffles inside each warp, then five more over the 32 warp totals. The per-thread sum of its grid-stride elements is a separate chain of length $n / 1024$, and that chain is usually the dominant term for large $n$.

Within a warp there are no barriers and no shared memory: the levels are five shuffle-and-add steps. Across warps, the partial sums are written to a small shared array, a barrier separates the write from the read, and one warp reduces them. The barrier is the expensive part of that stage, and its cost is set by the slowest warp in the block, as described in block barriers and divergent synchronization. Across blocks, either a second kernel reduces the per-block partials, or atomics accumulate them, or the last block to finish does the final sum.

The order of operations, concretely

Consider summing a warp’s 32 values with __shfl_down_sync. The first step adds, for each lane $k$, the value from lane $k+16$ to its own. After five steps lane 0 holds the sum. Which additions were performed, and in what order?

// Warp-level sum: five levels, no shared memory, no barrier.
float warp_sum(float v) {
    for (int off = 16; off > 0; off >>= 1)
        v += __shfl_down_sync(0xffffffffu, v, off);
    return v;   // valid in lane 0; other lanes hold partial sums
}

The result is a balanced binary tree in which element $i$ is added to element $i+16$ at the first level, $i+8$ at the second, and so on. Every element participates in the same number of additions, and the operands at each node are always one element and the partial sum of the block of elements immediately adjacent to it. This is the pairwise summation order, and it is the best-behaved of the simple orders for error growth: the error bound is $O(\log n \cdot \epsilon)$ rather than $O(n \cdot \epsilon)$ for the naive sequential sum.

The constants matter for reproducibility, not just accuracy. The exact grouping depends on the block size and on how the per-thread partials were formed, so a reduction launched with 256 threads and one launched with 1024 threads produce different bit patterns from the same input. The sum changes if the element-to-thread mapping changes at all, whether through the grid size, the number of elements each thread accumulates, or the order in which a single thread walks its own elements. This is inherent to the parallel form, not a bug in it: the specification of + does not define an order, so any order is conforming. Programs that need a reproducible sum from an unordered reduction need a different approach — a fixed reduction tree, integer or fixed-point accumulation, or a compensated sum — and no amount of care in the kernel makes the tree itself deterministic across launch shapes.

Error growth and when it matters

A sequential sum of $n$ positive terms has a worst-case error that grows like $n\epsilon$ times the sum of absolute values, with $\epsilon$ the machine epsilon of the format, about $1.2 \times 10^{-7}$ for single precision and $2.2 \times 10^{-16}$ for double. For a million single-precision terms the coefficient is $10^6 \times 1.2 \times 10^{-7} \approx 0.12$, so the bound is a substantial fraction of the sum of absolute values — enough that a naive float sum of a million values cannot be trusted to be accurate to even a few digits. Pairwise summation’s $\log_2 n$ factor reduces the coefficient by about 50000 at $n = 10^6$.

The worst case is dominated by cancelling terms: a sequence of large positive and large negative values whose true sum is small. Catastrophic cancellation is not fixed by the order, because every order loses the low bits when a large and a small value are combined. What the tree fixes is the accumulation of rounding error across many additions, not the loss of significance in any individual one. A sum that is small compared with its terms is inaccurate under every order; the remedy is a compensated algorithm, not a different tree.

The practical rules are few. Single precision is adequate when the terms are comparable in magnitude and the count is modest, and inadequate when the dynamic range of the terms spans many orders of magnitude. Accumulating partial sums in double from single-precision inputs costs a factor of two in shared-memory and register pressure and removes most of the concern. The compiler’s default is not to reassociate floating-point operations, so a kernel written with a sequential accumulation will genuinely be sequential unless the compiler is asked otherwise; a tree reduction written explicitly is the only way to get tree behaviour without that permission.

Composing the stages

A complete reduction is a stack of three stages, and each has a different cost structure.

The intra-thread stage sums the elements the thread visited in its grid-stride loop. This is a dependency chain, so its cost grows with the number of elements per thread, and it is the stage that dominates when the grid is small and the problem is large — the trade between trips and parallelism set out in grid-stride loops and launch size independence. For an unordered sum this stage can be split into two or four independent accumulators that are combined at the end, which shortens the chain by the same factor.

The intra-warp stage is five shuffles, with no barrier and no shared memory. It reduces 32 values to one per warp.

The inter-warp stage takes one value per warp and reduces them. If the block has 8 warps, it is 8 values in shared memory reduced by one warp — or a log-depth sequence of 3 shuffle steps over the 8 lanes if the warp reducing them is arranged to hold them in lanes 0–7. The shared array is small, and the barrier costs more than the arithmetic it guards.

Across blocks, the cheapest mechanism depends on the size of the partial array. A few hundred partials can be accumulated with an atomic add on a single location, which serialises the additions and therefore also makes the sum order depend on arrival order — non-deterministic across runs. A second kernel with one block reduces the partials in a fixed order, which is deterministic given the first kernel’s output. The single-pass alternative, where the last block to arrive performs the final reduction, is the same idea folded into one launch.

The last block does the last stage

The single-pass form needs to know when every block has finished. A counter in global memory is incremented by each block with a device-scope atomic after its partial is written; the block that observes the increment reaching the grid size reads all the partials and reduces them. The read of the partials must be ordered after the atomic, which requires an acquire-release pair rather than a plain counter.

Two details make this work in practice. The counter must be reset before each launch, which means zeroing it on the host or in a preceding kernel; a counter that is reused across launches without resetting will fire the final stage at the wrong time. And the block that performs the final reduction must not do so inside a region where its threads have already passed a barrier that other threads of the same block have not — the same rule as any other barrier, and the reason the final stage is usually a separate small section after all threads of the last block have arrived at a block barrier.

Integer and exact sums

None of the above applies to integer sums, and this is worth stating because integer reductions are common. Integer addition is associative and exact modulo the width of the type, so any order gives the same result, and the only concern is overflow. A sum of 32-bit signed integers overflows once the running total passes about 2.1 billion, which for a sum of byte values takes only about eight million elements. Widening to 64 bits for the accumulator removes the problem, and the reduction becomes a plain tree with no ordering concerns at all.

The same is true of floating-point sums that stay inside the exact range of the format — a sum of small integers represented in double is exact until it exceeds $2^{53}$ — which is why some algorithms can use a floating-point reduction safely and others cannot. The distinction is not the format but whether the sum of partial results exceeds the precision of the representation at any point.

The inter-warp and cross-block stages also interact with the number of blocks. A reduction over 4096 blocks has 4096 partials, so the final stage is itself a reduction of moderate size and needs to be a proper kernel rather than a single atomic. The threshold where one becomes the other is a function of the atomic rate on the device, and it is the kind of limit that has to be measured rather than assumed — which is why the two-kernel form, whose cost is one extra launch, is the safe default for a library routine and the single-pass form is the optimisation for a hot path.