Why parallel sums disagree: tracing FP32 reduction order

A sum has two specifications

A reduction must produce a sum without reading invalid memory, and it must define what numerical agreement means. Those are separate obligations. Real-number addition is associative; floating-point addition is not. A correct parallel tree can differ from a serial loop without either implementation containing an indexing bug.

Begin with a deliberately ordinary two-stage design. Each block loads one bounded input segment into shared memory, substitutes zero for missing elements, and reduces a power-of-two buffer. Every thread reaches each block barrier. The block's first thread writes one partial sum; another kernel reduces the partials, repeating if necessary. NVIDIA's reduction walkthrough explains this decomposition. Its historical optimization ladder is useful, but old implicit warp-synchronous assumptions are not an implementation contract for new code.

A small rounding counterexample

Consider FP32 inputs [100000000, 1, -100000000, 1]. A left-to-right accumulation loses the first unit when it adds it to the large positive value, cancels the large terms, then retains the final unit. Its result is one. A pairwise tree first combines each large value with a unit, loses both units, and returns zero. The real-number result is two.

The CPU demo explicitly rounds every arithmetic operation through a four-byte float using Python's struct module. It compares both orders with math.fsum, and then reorders the same four values to expose the dependence on grouping. This is a numerical illustration, not an emulation of every compiler choice: fused instructions, subnormal handling, and device math settings require their own examination. Even a high-precision reference needs a documented input-conversion policy.

Optimize the collective, not its assumptions

Once the shared-memory baseline passes, one isolated change is to reduce within a warp using register shuffles. This unexecuted CUDA device helper is valid only when all 32 lanes call it with the same full mask. The block size must be a multiple of 32, and out-of-range input lanes must contribute zero rather than return early.

__device__ float sum_full_warp(float value) {
    for (int distance = 16; distance > 0; distance >>= 1)
        value += __shfl_down_sync(0xffffffffu, value, distance);
    return value; // lane zero owns the complete sum
}

The mask declares participating threads; it does not supply missing operands. An arbitrary partial mask with this unchanged shuffle tree is not a general partial-warp reduction. If only some lanes have data, keeping the whole warp participating with zero-filled inactive data is the simplest contract here. NVIDIA's warp-primitives article describes synchronized collectives and why participation must be explicit. Shared-memory communication between warps still needs an appropriate block barrier.

What would count as progress?

Test empty-input host handling, a single element, nonmultiple lengths, alternating signs, and cancellation-heavy arrays. Report absolute error and an error normalized by the sum of input magnitudes; relative error against a nearly zero result can mislead. Choose tolerances before looking at benchmark winners.

Then time the baseline and shuffle variant on the same device, dtype, sizes, and stream, including every pass in the reduction. Compare with a mature device reduction, not just the slowest teaching kernel. Record warmup, repetitions, and accumulation precision. A faster result with materially different error is a trade-off, not a free optimization. The next useful extension is higher-precision accumulation, measured as a separate numerical and performance change; no GPU latency or speedup has been established here.

Enumerate the input orders, not just the striking example

The extended test enumerates all 12 distinct permutations of the four converted FP32 inputs. Each row retains the serial result, pairwise result, high-precision reference, absolute errors, and sum of input magnitudes. The multiset is unchanged, so any result change comes from ordering and grouping rather than different data.

The values 100000000, −100000000, and 1 are exactly representable in FP32. This makes input conversion unambiguous in the selected example: the lost units arise during accumulation. A fair reduction comparison needs the same clarity for arbitrary inputs. If one implementation converts or accumulates at a different precision, numerical disagreement may be an intentional specification difference rather than a synchronization defect.