Foundational kernel exercises

Harris’s fifth optimization, unrolling the last warp, is a hint the deck leaves half-explained. Once the tree has narrowed to 32 active threads they all belong to one warp, and the values they are still combining are single registers held by neighboring lanes. Routing those values through shared memory costs a store, a load, and a register to hold the address at every step. CUDA 9 replaced that with a family of warp-level primitives that move a value straight from one lane’s registers to another’s. The whole reduction becomes three lines:

#define FULL_MASK 0xffffffff
for (int offset = 16; offset > 0; offset /= 2)
    val += __shfl_down_sync(FULL_MASK, val, offset);

A warp is 32 lanes and each thread occupies one. For a thread at lane X, __shfl_down_sync(FULL_MASK, val, offset) returns the value of val held by the thread at lane X plus offset of the same warp, so halving the offset from 16 down to 1 walks the same tree the shared-memory version walks and leaves the sum in the first thread of the warp. Nothing is stored, nothing is loaded, no barrier separates the steps: the exchange happens between registers, which NVIDIA notes is more efficient than going through shared memory, where the same step needs a load, a store, and an extra register for the address. The same family covers the warp’s other collectives, __shfl_up_sync and __shfl_xor_sync for shifted and butterfly patterns, __ballot_sync for a bitmask of which lanes passed a predicate, and __match_any_sync for partitioning lanes by a value they hold.

The first argument is the part that trips people. The mask names which lanes must participate, one bit per lane ID, and the intrinsic waits until every non-exited thread named in it reaches the call. The programming guide states the contract as four rules: each calling thread must have its own bit set, each non-calling thread must have its bit set to zero, every non-exited thread in the mask must execute the intrinsic with the same mask value, and lanes may call concurrently with different masks only when those masks are disjoint. Violate it and the behavior is invalid or undefined, up to and including a kernel hang. Efficiency is best when all 32 lanes participate and the mask is 0xFFFFFFFF.

When only some lanes should participate, the mask has to be computed before the branch that separates them rather than inside it. NVIDIA’s post shows the correct shape: __ballot_sync evaluated against the branch predicate produces the membership mask, and the reduction inside the branch uses that mask. Calling __activemask() inside the branch instead looks equivalent and is not, because the execution model does not guarantee that all threads taking a branch together execute __activemask() together. The reduction then combines a subset of the lanes and returns a partial sum, a result no test on a small input is likely to catch.

The mask exists because the pre-CUDA-9 primitives did without it. __shfl, __ballot, and __any took no mask and performed no synchronization, so which threads participated was never written down anywhere, and correctness rested on implicit warp-synchronous behavior that may change from one architecture to another, from one toolkit release to another as compiler optimizations change, or even from one run to the next. The compiler and the hardware do try to reconverge threads after a divergent branch, but that reconvergence is not guaranteed. Those legacy primitives are deprecated starting in CUDA 9.0, and bracketing one with __syncwarp() calls does not rescue it: convergence is guaranteed only within an explicitly synchronous primitive, so threads that leave a __syncwarp() together may diverge again immediately.

__syncwarp() still earns its place where the data really does move through shared memory: it is __syncthreads() at warp granularity and provides a memory fence. The catch is that it must separate reads from writes, not merely appear between steps. A tree reduction written as a sequence of shmem[tid] += shmem[tid+16]; __syncwarp(); lines still races, because the model does not guarantee that all the reads of a step are performed before all its writes; the fixed version accumulates into a register, synchronizes, stores back, and synchronizes again. On Volta and later the exchange primitives may also be called inside thread-divergent branches, which is what makes it safe to put a shuffle inside a library function without knowing the caller’s control flow.

Scan is the third exercise on the list and the only one whose parallel form is not obvious. Given a list of inputs and a binary reduction operator, a prefix scan produces an output list in which each output is the reduction of the elements occurring earlier in the input. With addition as the operator it is the prefix sum, each output being the sum of the numbers before it. An inclusive scan folds the i-th input into the i-th output; an exclusive scan does not. Sequentially it is a single pass carrying a running total, n reads and n writes, which Merrill and Garland call the optimal 2n data movement. The difficulty is that every output depends on every earlier input, a serial dependence running the length of the array, and every parallel version buys its parallelism by computing something the sequential algorithm never computes.

It is worth the trouble because scan is the primitive underneath a whole class of irregular operations: adder design, linear recurrence and tridiagonal solvers, parallel allocation and queuing, list compaction and partitioning, segmented reduction. An exclusive prefix sum across a list of allocation requirements is the corresponding list of allocation offsets; run over an array of binary flags it gives the output positions for a stream compaction, since the scatter offset for a surviving item is the count of items kept before it.

The textbook parallel version, the Kogge-Stone construction and the corresponding Hillis-Steele algorithm, is recursive doubling: at each step every element adds the element a fixed distance behind it, and that distance doubles until it exceeds the array. Its depth is the minimum possible, and its work complexity is O(n log2 n) against the n-1 additions the serial algorithm performs. That inefficiency is real, but it costs least exactly where inactive processor resources cannot be scavenged for other work anyway, which describes a warp: shallow depth and simple shared-memory address computations make it the attractive strategy for SIMD widths. The programming guide’s own warp scan is that construction written with the shuffles from the previous subsection:

__global__ void scan_sub_partition_with_8_threads_kernel() {
    int laneId = threadIdx.x % 32;
    int value  = 31 - laneId; // starting value to accumulate
    // Loop to accumulate scan within my partition.
    // Scan requires log2(8) == 3 steps for 8 threads
    for (int delta = 1; delta <= 4; delta *= 2) {
        int tmp         = __shfl_up_sync(0xFFFFFFFF, value, delta, /*width=*/8); // read from laneId - delta
        int source_lane = laneId % 8 - delta;
        if (source_lane >= 0) // lanes with 'source_lane < 0' have their value unchanged
            value += tmp;
    }
    printf("Thread %d final value = %d\n", threadIdx.x, value);
}

The loop doubles delta where the reduction halved its offset, and each lane adds the value from the lane delta positions below it with __shfl_up_sync. Lanes whose source lane would be negative keep their value unchanged, which is what makes the result a scan rather than a shifted sum. The width argument of 8 splits the warp into four sub-partitions that each behave as a separate entity and scan in parallel, so the loop finishes in three steps. The guide suggests reaching for cub::WarpScan rather than hand-rolling the loop, the same advice it gives for warp reductions.

Above the warp, work complexity starts to matter, and the answer is a tree with two halves. The Brent-Kung construction, and Blelloch’s algorithm built on it, is the work-efficient strategy: 2 log2 n depth and O(n) size. Merrill and Garland describe its data flow as an hourglass, an upsweep accumulation tree exhibiting progressively less parallelism as it narrows to a single value, then a downsweep propagation tree exhibiting progressively more as it widens back out. It pays twice the depth of Kogge-Stone for linear work, the right trade the moment the input is larger than the machine and the lanes idled by the upsweep can be given something else to do.

For a device-wide scan, though, neither depth nor work complexity picks the winner, because the kernel is memory bound and the score is how many times an implementation touches memory. Merrill and Garland line the strategies up by that measure. Scan-then-propagate, the recursive high-radix dataflow in CUDPP and Thrust, incurs about 4n of global data movement. Reduce-then-scan, the three-kernel structure of an upsweep, a scan of the block aggregates, and a downsweep, incurs about 3n and makes two full passes over the input. The sequential algorithm’s 2n is the floor. Two full passes also rule those strategies out of in-place compaction, since the execution order of thread blocks in the output pass is unconstrained and separate storage is required to stop an input being overwritten before it has been read.

Chained scan reaches 2n in a single pass by giving each thread block one tile and one serial dependence: a block waits on the inclusive prefix of its predecessor, adds its own aggregate, and passes the running total along. The dependence is also the ceiling, and the paper works out how low that ceiling sits. On a Tesla C2050 with 144 GB/s of memory bandwidth, cores at 1.15 GHz, and a latency of 600 clock cycles to pass a message from one core to another, signal propagation alone limits the machine to 1.9M partitions per second; at 256 32-bit inputs per partition that is 490M inputs per second, well short of the 18 billion inputs per second the memory system could sustain for simply reading and writing the data. Hiding the signalling would take a partition of around 9,000 inputs per thread block, and partition size cannot be raised arbitrarily: a GPU core performs best when its working set fits entirely inside the on-chip register file and shared memory.

Decoupled look-back keeps the single pass and removes the serial wait. Every block still computes and records its tile aggregate independently, but instead of blocking on its immediate predecessor it inspects the descriptors of predecessors increasingly further away. Each partition owns a status descriptor with three fields: aggregate, the partition-wide reduction; inclusive_prefix, the running total across all prior partitions; and a status flag that is X while no information has been recorded, A once the aggregate has been, and P once the inclusive prefix has. A block walking backward blocks or keeps polling on X, adds the aggregate to its running exclusive prefix and steps one further back on A, and on P adds the inclusive prefix and terminates the walk, because a predecessor that knows its own inclusive prefix has already summarized everything before it. That short circuit bounds the walk: look-back is constant given a finite number of processors, which keeps overall work complexity at O(n), and a whole SIMD group can inspect a window of predecessors at once, each thread monitoring its own.

The descriptors need the discipline any lock-free structure needs. Recording an aggregate or a prefix takes three memory operations in order: update the value field, execute a memory fence, then update the status flag, so that no peer ever sees a flag claiming more than the data behind it. The fence is not free, since it stops the write to the flag being pipelined with the write to the value, and it can be eliminated wherever flag and value fit in one architectural word: for prefix sums over 32-bit integers on hardware with 64-bit loads and stores, a single 64-bit write of the pair is enough.

The payoff is that a device-wide scan stops being an operation worth routing around. Measured on Fermi, Kepler, and Maxwell parts, the CUB implementation meets or exceeds every other implementation at every problem size, and on the Kepler and Maxwell parts it matches, for large problems, the performance ceiling of memcpy, an operation that shares the same minimum I/O workload and performs no computation at all. Its harmonic-mean speedups across all problem sizes are 1.60× over StreamScan, 1.19× over Modern GPU, and 2.80× over Thrust, and the compaction primitives built on the same single pass run 4.1×, 7.1×, 3.5×, and 3.8× faster than Thrust’s select-if, partition-if, reduce-by-key, and run-length-encode.

One property comes along with it. Where the other strategies are static dataflow networks whose order of operations is fixed for a given problem size, look-back allows each processor to perform whatever redundant work it needs to avoid waiting, so the number and order of scan operators applied is not deterministic. For floating-point data that means the prefix sum across a given dataset may vary from one run to the next: nondeterminism accepted deliberately in exchange for never stalling on a predecessor.