GPU-Specific Algorithms, in depth: SIMT constraints, work and span, and the reduce, histogram, compaction and radix-sort toolkit

The execution model as an algorithm designer sees it

A CUDA kernel is launched as a grid of thread blocks. Each block runs on one streaming multiprocessor (SM), its threads can share a fast on-chip shared memory, and they can wait for each other at __syncthreads(). Blocks cannot wait for each other inside a kernel unless you use special cooperative launches; in ordinary code the only global barrier is the end of the kernel. Inside a block, threads are executed in warps of 32 on NVIDIA hardware. Code should read the built-in warpSize rather than hard-code it if it may run elsewhere; AMD's wavefronts are 64 wide on many parts.

Three consequences drive every design below. Divergence: when lanes of a warp take different branches, the warp runs both paths with lanes masked off, so a data-dependent if can halve throughput. Coalescing: when the 32 lanes load 32 consecutive words, the hardware serves them with a few wide transactions; when they load scattered addresses, each needs its own. Latency hiding: a global load takes hundreds of cycles, and the SM hides that by switching among many resident warps, so you need enough independent warps (occupancy) and enough independent work per warp. The SIMT execution model and coalesced access articles cover the hardware side in detail.

Cost: work, span and bytes moved

Two numbers describe a parallel algorithm. Work W is the total number of operations; span S is the length of the longest chain of dependent steps. Brent's bound says p processors finish in at most W/p + S steps. On a GPU, p is huge, so the rule of thumb is: keep W within a constant factor of the best sequential algorithm (work-efficient), and make S polylogarithmic. A naive parallel algorithm that does O(n log n) work for an O(n) problem looks fine on paper and loses to a good CPU loop in practice, because the GPU is usually limited by memory bandwidth, not by arithmetic.

That last point is the second rule. Most primitives in this article do a handful of operations per element loaded, so their speed is set by how many bytes they move. Count passes over memory: a design that reads the input once and writes the output once beats a design with a lower operation count that makes three passes. This is why libraries fuse steps aggressively and why a good GPU sort minimises the number of digit passes.

Reduction, the template for everything else

Summing an array shows the whole method in miniature. The obvious tree, where thread i adds element i + stride in shared memory, works, but it synchronizes the block at every level and leaves most threads idle. The standard shape is hierarchical. Each thread first sums many elements sequentially with a grid-stride loop, which is work-efficient and coalesced because consecutive lanes read consecutive addresses. Then warps combine in registers with shuffle instructions that need no shared memory and no barrier. Only the 32 or fewer warp totals go through shared memory.

__inline__ __device__ float warp_sum(float v) {
    for (int off = 16; off > 0; off >>= 1)          // 32 lanes: 5 shuffle steps
        v += __shfl_down_sync(0xffffffff, v, off);
    return v;                                        // lane 0 holds the warp total
}

// Launch with blockDim.x a multiple of 32 (at most 1024).
__global__ void sum_kernel(const float* __restrict__ x, float* out, int n) {
    __shared__ float warp_totals[32];
    float v = 0.0f;
    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n;
         i += gridDim.x * blockDim.x)
        v += x[i];                                   // grid-stride, coalesced
    v = warp_sum(v);
    int lane = threadIdx.x % 32, warp = threadIdx.x / 32;
    if (lane == 0) warp_totals[warp] = v;
    __syncthreads();
    if (warp == 0) {
        v = (lane < blockDim.x / 32) ? warp_totals[lane] : 0.0f;
        v = warp_sum(v);
        if (lane == 0) atomicAdd(out, v);            // see the determinism note
    }
}

The final atomicAdd makes the kernel a single launch, but floating-point addition is not associative and the order in which blocks finish varies, so the last bits of the result can change between runs. For reproducible training metrics or tests, write one partial per block and reduce those in a second, fixed-order pass. Size the grid to a small multiple of the SM count rather than one block per element range; the grid-stride loop then absorbs any n.

Histograms and privatization

A histogram is a reduction with many outputs and a data-dependent address, so the naive version, one global atomicAdd per element, collapses when many elements fall into the same bin: the atomics to that address serialize. The fix is privatization. Each block counts into its own copy in shared memory, where atomics are much cheaper, and merges into the global histogram once at the end, touching each bin once per block instead of once per element.

__global__ void hist256(const unsigned char* in, unsigned int* hist, int n) {
    __shared__ unsigned int local[256];
    for (int b = threadIdx.x; b < 256; b += blockDim.x) local[b] = 0;
    __syncthreads();
    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n;
         i += gridDim.x * blockDim.x)
        atomicAdd(&local[in[i]], 1u);
    __syncthreads();
    for (int b = threadIdx.x; b < 256; b += blockDim.x)
        if (local[b]) atomicAdd(&hist[b], local[b]);
}

Test the skewed case explicitly. A uniform input hides contention, and a constant input (every element in one bin) is the worst case. If the bin count is too large for shared memory, sort or partition first, which brings in the next primitive.

Scan-built primitives: compaction and radix sort

Scan (prefix sum) is the primitive that turns "each thread decides independently" into "each thread knows where to write". Stream compaction keeps the elements that pass a predicate, densely packed. Every thread writes a flag, an exclusive scan of the flags gives each survivor its output index, and a scatter writes it there. LSD radix sort repeats the same pattern per digit: histogram the digit, exclusive-scan the counts to get each digit's starting offset, then scatter every key to its offset plus its stable rank among equal digits. The NumPy code below is a CPU reference model of that data flow, with one line per GPU kernel, and is useful as a test oracle. It matched np.sort and boolean masking on 50 random inputs each.

import numpy as np

def compact(x, keep):
    flags = keep.astype(np.int64)                    # map
    pos = np.cumsum(flags) - flags                   # exclusive scan
    out = np.empty(int(flags.sum()), dtype=x.dtype)
    out[pos[keep]] = x[keep]                         # scatter survivors
    return out

def radix_sort_lsd(keys, bits=32, r=4):
    keys = np.asarray(keys, dtype=np.uint32)
    for shift in range(0, bits, r):                  # one pass per r-bit digit
        digit = (keys >> shift) & ((1 << r) - 1)
        counts = np.bincount(digit, minlength=1 << r)    # 1. histogram
        base = np.cumsum(counts) - counts                # 2. exclusive scan
        rank = np.empty(len(keys), dtype=np.int64)       # 3. stable rank in digit
        for d in range(1 << r):
            m = digit == d
            rank[m] = np.arange(m.sum())
        out = np.empty_like(keys)
        out[base[digit] + rank] = keys                   # 4. scatter
        keys = out
    return keys

For example, compact([3, -1, 4, -1, 5, -9, 2, 6], x > 0) gives [3, 4, 5, 2, 6]. Sorting the 8-bit keys 0x23, 0x11, 0x2A, 0x13, 0x21 takes two 4-bit passes. The low-digit pass orders them 0x11, 0x21, 0x23, 0x13, 0x2A, and the stable high-digit pass yields 0x11, 0x13, 0x21, 0x23, 0x2A. Stability of each pass is what makes LSD correct. On the GPU the rank step is itself a per-block scan, and production sorts use wider digits and fuse the scatter with local sorting in shared memory so that writes stay coalesced.

Irregular work: sparse matrices and graphs

Regular primitives are easy; the hard GPU algorithms are irregular. Sparse matrix-vector multiply in CSR format is the standard example. Assigning one thread per row is coalesced for the vector of row offsets but not for the matrix values, and a single very long row stalls a whole warp while its neighbours idle. Assigning one warp per row fixes coalescing but wastes lanes on short rows. Graphs with power-law degree have both kinds of row at once. Load-balanced methods such as merge-path SpMV (Merrill and Garland) instead split the combined list of nonzeros and row ends into equal chunks per thread, so work is balanced regardless of row lengths.

Graph traversal follows the same arc. A level-synchronous BFS keeps a frontier array, expands it in parallel, and uses compaction to build the next frontier. When the frontier grows to a large fraction of the graph, switching to a bottom-up step, in which unvisited vertices look for any visited parent, saves most of the edge checks. The general lesson: when work per item varies, partition by work, not by item.

Use the libraries

Writing these kernels teaches the model; shipping them is usually a mistake. CUB provides tuned block-, warp- and device-wide primitives, and Thrust wraps them in an STL-like interface (thrust::reduce, thrust::sort, thrust::copy_if). CUB's device-wide calls use a two-call pattern: the first call with a null buffer only reports how much temporary storage is needed.

void*  d_temp = nullptr;
size_t temp_bytes = 0;
cub::DeviceRadixSort::SortKeys(d_temp, temp_bytes, d_keys_in, d_keys_out, n);
cudaMalloc(&d_temp, temp_bytes);                     // or reuse a pooled buffer
cub::DeviceRadixSort::SortKeys(d_temp, temp_bytes, d_keys_in, d_keys_out, n);

In training frameworks the same primitives are hidden inside operators: top-k sampling is a partial sort, embedding-bag backward is a sort followed by a segmented reduction, and mixture-of-experts token routing is a histogram, a scan and a scatter. Knowing which primitive an operator becomes tells you why it is slow.

Operational guidance

• Profile first. Compare achieved memory throughput against the device peak; a bandwidth-bound primitive at a small fraction of peak usually has uncoalesced accesses or too few resident warps.
• Check occupancy limits from registers and shared memory per block; the occupancy article shows how to read them.
• Keep a CPU reference implementation, like the NumPy models above, and diff against it in tests over random, sorted, constant and empty inputs.
• Decide on determinism per kernel and document it. Atomics on floats are not reproducible; integer atomics are.
• Batch small problems. Launch overhead dominates tiny kernels, so merge work or use CUDA graphs.

Failure modes

• Divergent barrier. __syncthreads() inside a branch that not every thread of the block reaches is undefined behaviour and commonly hangs.
• Implicit warp synchrony. Code that assumes lanes run in lockstep without _sync intrinsics or __syncwarp() broke when independent thread scheduling arrived. Always pass the lane mask.
• Bank conflicts. Strided shared-memory access maps many lanes to one bank and serializes them. Pad arrays by one word or change the indexing.
• Integer overflow in indices. blockIdx.x * blockDim.x in 32-bit arithmetic overflows past about two billion elements. Use 64-bit indices for large tensors.
• Hot-bin atomics. Skewed keys serialize global atomics. Privatize, as the histogram above does.

Trade-offs

What to do next

1. Run the NumPy reference models, then write the CUDA reduction and diff its output against np.sum in float64.
2. Benchmark the atomic and two-pass endings, and measure run-to-run variation in the last bits.
3. Build the privatized histogram and time it on uniform and constant inputs.
4. Replace your hand kernels with CUB or Thrust calls and compare throughput.
5. Profile one slow operator in your training job and name the primitive it reduces to.