GPU Parallel Primitives Lab

Eight numbers sit in GPU memory: 3, -1, 4, 1, 5, -9, 2, 6. A transpose can move every number to a new address. It can't answer three common questions: What is their total? How many values came before each position? Which value is largest, and where did it occur?

Those questions need parallel primitives, small algorithms that many larger kernels reuse. Reduction turns many inputs into fewer aggregates. Prefix scan returns an aggregate at every boundary. Arg reduction carries a value and its location together. Stable softmax combines reductions with exponentials without overflowing.

The GPU Kernel Performance Engineering Lab established coalescing, shared memory, barriers, exact checks, and event timing through one transpose. Keep those contracts. Here one stable tensor exposes new dependency shapes instead of a new memory layout.

Work eight values by hand

Start with sum. A serial loop visits all eight values in order and ends at 11. A balanced tree reaches the same integer through three parallel levels:

Eight leaves require seven additions either way. Tree depth is only $log_2 8=3$, provided enough workers can execute each level together. That shrinking output is reduction's signature.

Prefix scan keeps every boundary instead. For addition, an inclusive scan includes current value at index $i$. An exclusive scan stops before it:

The relation is exact for every valid $i$:

$$\operatorname{inclusive}[i]=\operatorname{exclusive}[i]+x_i$$

Argmax needs another payload. Maximum value is 6 and its index is 7, so output is pair (6, 7). If tensor contained another 6, operator would need a tie rule such as "smaller index wins." Value, index, and tie rule form one reduction state.

The runnable model below checks all four outputs and constructs stable softmax state used later. It uses Python's math module so mechanism stays visible.

Move sum into one warp

A naive GPU baseline assigns all eight additions, or all $n$ additions, to one thread. That gives a correct, easy-to-check answer but preserves serial dependency depth and leaves rest of GPU idle. Hierarchical version begins by giving each lane a local sum.

A warp is CUDA's group of 32 lanes that executes instructions together. __shfl_down_sync lets one participating lane read a register from a higher-numbered lane without routing value through shared memory. Five offsets, 16, 8, 4, 2, 1, reduce 32 lane-local values to lane 0.^{[1]Reference 1CUDA Programming Guide.https://docs.nvidia.com/cuda/cuda-programming-guide/}

The lab uses this helper:

Full mask is correct only because every lane reaches each shuffle. Tail lanes load additive identity 0 instead of exiting. CUDA defines a shuffle read from an inactive source lane as undefined, and participating lanes named by mask must execute same intrinsic with same mask.^{[1]Reference 1CUDA Programming Guide.https://docs.nvidia.com/cuda/cuda-programming-guide/}

That contract rules out a tempting edge shortcut:

On partial warp, mask names lanes that skipped call. Repair by keeping all lanes active and substituting identity for missing input, or compute one valid mask with __ballot_sync before lanes diverge and ignore shuffle results whose source lane isn't valid. Calling __activemask() after divergence can give participating lanes different masks. Mask says who participates. It doesn't turn invalid source lane into zero.

Only 20 values remain in final warp. Why is full mask safe in supplied lab but unsafe when lanes 20 through 31 return early?

Answer

Lab keeps all 32 lanes at every shuffle and gives lanes 20 through 31 value 0, so every named lane participates. Early return violates mask contract because full mask still names absent lanes. Use uniform participation with identity values or a correctly formed active mask plus valid-source guards.

Join warps inside one block

A 256-thread block contains eight warps. Each warp first reduces registers. Lane 0 of every warp writes one partial to eight-element shared array. One block-wide barrier then makes those writes visible before warp 0 reduces eight partials.

__syncthreads() waits for every thread in block and makes earlier shared and global memory accesses visible to that block. Placing it under a condition that differs across block can hang or produce unintended behavior.^{[1]Reference 1CUDA Programming Guide.https://docs.nvidia.com/cuda/cuda-programming-guide/}

Older teaching reductions begin with interleaved active threads, then replace divergent modulo branch with contiguous active lanes. They also give each thread several elements before shared tree begins. Those steps reduce branch divergence, address work, and synchronization while keeping $O(n)$ total work.^{[2]Reference 2Optimizing Parallel Reduction in CUDAhttps://developer.download.nvidia.com/compute/cuda/1.1-Beta/x86_website/projects/reduction/doc/reduction.pdf} Current NVIDIA samples still include a two-stage shared-memory reduction as runnable reference.^{[3]Reference 3CUDA Sampleshttps://github.com/NVIDIA/cuda-samples}

One block can't synchronize with every other ordinary block inside kernel. Hierarchical reduction therefore writes one partial per block, then launches another reduction pass over partial array. Each launch on same stream forms device-wide ordering boundary:

Two practical details matter:

• Each thread accumulates two or more inputs before warp shuffle, balancing parallel depth against useful work per thread.
• Accepted operator needs associative grouping. Commutativity is optional only for algorithms that preserve operand order, such as scan. CUB DeviceReduce doesn't support non-commutative reduction operators.^{[4]Reference 4cub::DeviceReducehttps://nvidia.github.io/cccl/unstable/cub/api/structcub_1_1DeviceReduce.html}

For mathematical addition, (a+b)+c equals a+(b+c). Floating-point addition only approximates that law because each addition rounds. Tree can differ from serial sum without either kernel containing race.

Make floating-point reductions deliberate

Suppose float32 values include 100000000.0, 1.0, and -100000000.0. Serial grouping can lose 1 before cancellation, while another tree preserves it. Integer fixture avoided that ambiguity so CUDA lab could demand exact equality.

Use numerical contract suited to workload:

Reproducibility and accuracy answer different questions. A fixed wrong order can reproduce. A more accurate algorithm can still choose different trees on different architectures.

Current CUB DeviceReduce documentation states default run_to_run determinism. Same input, build, launch configuration, tuning, CCCL version, and GPU select same fixed tree. Another architecture, policy, or toolkit release can change floating-point combining order.^{[4]Reference 4cub::DeviceReducehttps://nvidia.github.io/cccl/unstable/cub/api/structcub_1_1DeviceReduce.html} Pin CCCL version and verify installed API when reproducibility enters product contract.

Two float sum kernels pass abs(error) < 1e-4. One returns same bits on every run, while other has lower error against FP64 but changes last bit across runs. Which one is correct?

Answer

Both satisfy stated accuracy tolerance. First also satisfies same-run reproducibility; second doesn't. If product needs both properties, neither single observation is enough. Write separate accuracy and determinism gates, then choose operator, accumulator, and library policy that pass both.

Keep every prefix with scan

Reduction's tree can discard intermediate totals. Scan must send block prefix into every later output. A work-efficient block scan has two phases:

1. Upsweep: combine leaves into block total.
2. Downsweep: replace root with identity, then propagate left-prefix information back to leaves.

Supplied Blelloch-style kernel handles 512 items per 256-thread block. Every thread loads two values into shared memory, participates in both tree phases, and writes exclusive output. It also records block total before root becomes 0.

Arrays longer than one block need three layers:

If first output of block 1 is 0 instead of sum of block 0, local tree probably passed and uniform add failed. That symptom points at propagation, not arithmetic inside block.

Blocks of four scan [3,-1,4,1 | 5,-9,2,6]. Local exclusive outputs are [0,3,2,6 | 0,5,-4,-2]. What offset repairs second block?

Answer

First block total is 7, so scan of block totals gives offsets [0,7]. Add 7 to every local output in second block, producing [7,12,3,5]. Combined output is [0,3,2,6,7,12,3,5].

Multi-pass scan is easy to inspect but moves block metadata in extra kernels. Merrill and Garland's decoupled look-back scan lets block compute local scan, publish aggregate and status, then look backward through earlier block states until global prefix is known. Small redundant work overlaps global prefix propagation with local work. Their report describes about $n$ input reads plus $n$ output writes, matching sequential scan's asymptotic data movement.^{[5]Reference 5Single-pass Parallel Prefix Scan with Decoupled Look-backhttps://research.nvidia.com/sites/default/files/pubs/2016-03_Single-pass-Parallel-Prefix/nvr-2016-002.pdf}

CUB has used decoupled look-back for device-wide scan since early releases. DeviceScan docs expose inclusive and exclusive forms, in-place support, and temporary-storage contract. Classic stream overloads warn that pseudo-associative results may vary across runs. Newer execution-environment controls depend on version, type, and operator: current unstable docs default to not_guaranteed and permit run_to_run only for supported integral operators and floating-point cuda::std::plus.^{[6]Reference 6cub::DeviceScanhttps://nvidia.github.io/cccl/unstable/cub/api/structcub_1_1DeviceScan.html} Check pinned CCCL docs or headers before turning a determinism mode into a release gate.

Decoupled look-back isn't invitation to improvise spin protocol. Published states, memory ordering, forward progress, fallback path, and architecture behavior all need proof. Use library for standalone device scan unless fusion or data layout gives custom kernel clear reason to exist.

Treat stable softmax as reduction state

Softmax turns logits $x_i$ into probabilities:

$$p_i=\frac{e^{x_i}}{\sum_j e^{x_j}}$$

Direct exponentiation can overflow for large positive logit. Stable form subtracts row maximum $m$:

$$m=\max_j x_j, \qquad p_i=\frac{e^{x_i-m}}{\sum_j e^{x_j-m}}$$

Subtracting same constant from all logits doesn't change probabilities, while largest exponent becomes $e^0=1$.

Safe implementation can read logits once for maximum, again for denominator, then again to write normalized outputs. Online normalizer combines first two reductions. State (m, d) represents partition maximum $m$ and shifted exponential sum $d$. Merge two partitions $a$ and $b$ by:

$$m=\max(m_a,m_b)$$ $$d=d_a e^{m_a-m}+d_b e^{m_b-m}$$

This operator is associative over finite logits in exact real arithmetic, so same warp and block reduction machinery can merge chunk states. Empty partition uses identity (-inf, 0), but merge must return other operand before evaluating formula. Otherwise two empty partitions evaluate -inf - -inf and poison denominator with NaN. Python model's two early returns enforce identity law. Milakov and Gimelshein derive online update, prove invariant, and reduce softmax memory accesses by combining maximum and denominator passes.^{[7]Reference 7Online normalizer calculation for softmax.https://arxiv.org/abs/1805.02867}

Attention masks expose same edge. A masked logit contributes empty identity, not singleton state (-inf, 1). Fully masked row has no mathematical softmax distribution, so kernel contract must reject it or define explicit fallback before normalization. NaN and positive-infinity logits also need declared input contract because subtracting maximum doesn't make inf - inf defined.

Running tensor ends with m=6 and denominator near 1.578968. Python model then revisits logits only to write eight probabilities. Fusing next consumer can sometimes avoid writing full probability vector, but correctness contract expands to fused operation.

Online doesn't mean one sequential thread. Each lane can summarize local chunk into (m,d), then parallel tree merges those states. Intermediate maximum can rise, so old denominator must be rescaled before new contribution arrives.

Left chunk state is (5, 2) and right state is (6, 1). What merged denominator is measured relative to maximum 6?

Answer

Merged maximum is 6. Left denominator rescales by $e^{5-6}=e^{-1}$, while right stays at scale 1. Merged denominator is $2e^{-1}+1$, about 1.7358. Adding 2 and 1 directly would mix different exponential scales.

Carry segments and indices when output needs them

Primitive choice follows output contract, not operator name.

• Segmented reduction: reduce several variable-length ranges independently, such as one maximum per sequence in packed batch.
• Arg reduction: return both winning value and index, with explicit tie and NaN rule.
• Reduce by key: combine adjacent runs sharing key, often after grouping.
• Scan by key: reset prefix when key changes.

For argmax, define state (value, index). Combine larger value; on tie, choose smaller index. Empty segment also needs documented identity or sentinel. Without those rules, parallel tree may return different valid-looking index after scheduling or tile change.

CUB provides device-wide, block-wide, and warp-wide collectives plus segmented variants. Scope controls who invokes operation and where output is valid. A BlockReduce result is meaningful only for designated thread; every lane writing it would duplicate or race output.^{[8]Reference 8CUB Developer Overviewhttps://nvidia.github.io/cccl/unstable/cub/developer_overview.html}

Prefer library until fusion earns custom kernel

Current CUDA Core Compute Libraries (CCCL) CUB documentation layers reduction across thread, warp, block, and device scopes. Device APIs choose tuned policies by type, size, and architecture; block APIs expose reusable temporary storage and synchronization rules.^{[8]Reference 8CUB Developer Overviewhttps://nvidia.github.io/cccl/unstable/cub/developer_overview.html}

Two-phase device reduction call looks like this. First call queries temporary bytes, second launches operation:

Real code must check every CUDA return, own stream explicitly, reuse temporary allocation across hot calls, and free it after stream work completes. Newer CCCL execution-environment APIs can manage stream-ordered temporary storage, but version pin decides available surface.

Use this boundary:

Custom kernel owns more than arithmetic: architecture policy, partial warps, temporary storage, stream semantics, deterministic order, sanitizer coverage, and future retuning.

Run downloadable CUDA lab

Download complete CUDA source. It contains four paths:

Baselines are intentionally serial. They expose dependency and produce trusted integer reference, but they aren't competitive implementations.

Compile for installed GPU architecture. Attached local receipt used sm_120; command below keeps target explicit:

Program prints hardware identity, API versions, workload, warmups, iteration count, exact correctness status, CUDA-event time, and logical byte rate. Reduction rate counts input bytes; scan rate counts one input read plus one output write. It isn't physical DRAM bandwidth because hierarchical paths also move partials and offsets.

Run odd size through correctness tools before trusting power-of-two benchmark:

Compute Sanitizer distinguishes invalid addresses, shared-memory hazards, and synchronization misuse.^{[9]Reference 9Compute Sanitizerhttps://docs.nvidia.com/compute-sanitizer/ComputeSanitizer/} Numerical checks still matter because sanitizer can't prove prefix values.

Local benchmark receipt records full raw output from RTX 5070 Ti with CUDA 13.2, source hash, compiler command, warm-cache timing method, exact CPU checks, and clean sanitizer results. It doesn't claim CUB comparison, physical bandwidth, or portability. Re-run same receipt on target hardware before using any timing decision.

Test properties, not one fixture

One eight-value tensor makes mechanism visible. Property tests catch edge behavior.

Reduction properties

• Empty input returns declared identity or documented error.
• One value returns that value.
• Split then merge matches full reduction for truly associative operator.
• Permutation leaves output unchanged only when operator is commutative.
• Integer accumulator doesn't overflow tested bound.
• Float error stays within stated tolerance against higher-precision reference.

Scan properties

• Output length equals input length.
• Exclusive output at index 0 equals identity.
• exclusive[i+1] == exclusive[i] op input[i] for exact operator.
• Last inclusive output matches full reduction.
• Segment reset occurs exactly at supplied boundary.
• Sizes around block boundary pass: 0, 1, 31, 32, 33, 511, 512, 513.

Softmax properties

• For an accepted finite row with at least one unmasked logit, every output is finite and nonnegative.
• Row sums to 1 within tolerance.
• Adding constant to every logit leaves probabilities unchanged within tolerance.
• Probability argmax matches logit argmax under same tie rule.
• Extreme fixtures such as [-10000, 0, 10000] don't produce NaN or Inf.
• Masked positions contribute empty identity, and fully masked rows follow declared error or fallback contract.

Property that fails often identifies broken layer faster than output dump.

Diagnose by first broken invariant

Stop timing at first correctness failure. Fast wrong prefix can corrupt allocation offsets far from kernel that produced them.

Keep complete benchmark receipt

Primitive timing is incomplete without five fields:

1. Hardware and software: GPU, compute capability, driver, toolkit, compiler flags, source hash.
2. Workload: count, segments, value distribution, alignment, warmups, iterations, cache state.
3. Precision and algorithm: input, accumulator, output types, operator, hierarchy, deterministic policy.
4. Baseline: exact serial, CUB, framework, or prior kernel path under equivalent work.
5. Correctness: reference, tolerance, property sizes, sanitizer results.

Measure kernel execution with CUDA events in same stream. Separate allocation, first-use loading, and JIT cost unless latency objective includes them. Report distributions across process runs when decision is close. Attach profiler counters only after timing and correctness identify meaningful candidate.

No universal speedup follows from primitive name. Reduction shape, value type, segment distribution, architecture, library version, fusion, and launch overhead can reverse ranking.

Mastery check

Key concepts

• Reduction shrinks output; scan preserves one prefix per position.
• Warp shuffles move register values only among valid participating lanes.
• Block collectives need shared-memory visibility and uniform barriers.
• Device-wide collectives need inter-block propagation through more launches or proven look-back protocol.
• Stable softmax merges (maximum, shifted denominator) states.
• Masked softmax values contribute empty identity; fully masked rows need explicit contract.
• CUB is default for standard standalone primitives; custom CUDA needs fusion or semantic reason.

Evaluation rubric

• Foundational: Compute reduction, inclusive scan, exclusive scan, and argmax for eight-value tensor by hand.
• Intermediate: Explain warp, block, and device hierarchy plus mask, barrier, identity, and offset invariants.
• Advanced: Implement and benchmark custom fused primitive against CUB with numerical, determinism, sanitizer, and receipt gates.

Common pitfalls

• Treating active mask as replacement for valid-source check.
• Returning before block-wide barrier on tail data.
• Assuming floating-point addition is exactly associative.
• Forgetting to scan block totals before uniform add.
• Merging online-softmax denominators measured against different maxima.
• Calling demonstration baseline competitive without tuned library comparison.