How CUDA CCCL works: Thrust, CUB and libcu++ in 39 graphs — article mark

How CUDA CCCL works: Thrust, CUB and libcu++ in 39 graphs

CCCL — the CUDA Core Compute Libraries — is the standard library of GPU programming: Thrust, CUB and libcu++ in one repository, NVIDIA/cccl. Almost all of it is C++ templates that compile straight into your program.

Five parts, mostly graphs: 1 the map · 2 one call followed down to the hardware · 3 the other algorithms · 4 libcu++, Python and cudax · 5 the math underneath, in plain words. Every number is read from the source at CCCL 3.6.0, commit 676b94d.

Part 1 — The map

CCCL architecture: Python, your C++ program and your own kernels enter at the top; Thrust and the JIT C library feed CUB; CUB is built on libcu++; cudax sits beside it; everything runs on the CUDA runtime, PTX and the GPU Python cuda.compute.reduce_into(…) Your C++ program (runs on the CPU) calls algorithms · launches kernels cccl.c.parallel JIT: fills in a CUB kernel, compiles it at runtime Thrust — what to compute sort · reduce · scan · transform device_vector · fancy iterators CPU backends: OpenMP · TBB · C++ Your own kernel __global__ void f(…) thousands of threads run it calls launches CUB — how to compute it fast speed-of-light building blocks Device-wide algorithms — called from the CPU DeviceReduce · DeviceScan · DeviceRadixSort · DeviceSelect · … Collectives — inside kernels BlockReduce · WarpScan · … JIT-built kernels CUDA backend or call CUB directly uses libcu++ — the C++ standard library, on the GPU cuda::std:: atomic · span · mdspan · complex · chrono · … cuda:: scoped atomics · barrier · memcpy_async · ptx:: · stream cudax (experimental) CUDASTF task graphs graph builder · places built on CUDA runtime & driver → PTX → GPU: SMs · shared memory · L2 · HBM PTX = NVIDIA's portable GPU assembly language; the driver turns it into machine code for your GPU
The whole of CCCL on one page. Three ways in at the top — Python, a C++ program, your own kernel — and they all funnel into CUB, which is built on libcu++, right above the CUDA runtime and the hardware. The C++ libraries are header-only templates compiled into your program; the Python and C layers ship as compiled libraries.

A few words carry the rest of the article. A kernel is a function the CPU launches on the GPU (the “device”). It runs as a grid of blocks; each block lands on one SM (streaming multiprocessor — an H100 has 132) and splits into warps of 32 threads, numbered lanes 0–31.

GPU execution hierarchy (grid, block, warp, thread) and memory hierarchy (registers, shared memory, L2, HBM) with H100 SXM numbers Grid — one kernel launch its blocks run in any order, on any SM block 0 block 1 block 2 … Block — up to 1,024 threads its threads can share memory and wait for each other warp 0 warp 1 warp 2 … warp 7 warp = 32 threads (lanes 0–31) running one instruction together 256 threads = 8 warps · every thread has its own registers shared memory — visible to this block only Where data lives — H100 SXM Registers one thread's own · 256 KB per SM Shared memory one block's scratchpad · ≤ 228 KB per SM L2 cache shared by the whole GPU · 50 MB HBM — global memory 80 GB at 3.35 TB/s · your arrays live here faster bigger 132 SMs (streaming multiprocessors) each runs up to 2,048 threads; every block lands on one SM
Enough hardware for the rest of the article, with numbers for an H100 SXM. A kernel launch is a grid of blocks; a block is split into warps of 32 threads that execute one instruction together. The closer data sits to a thread, the faster it is to reach.

To watch warps execute, diverge and coalesce, try the GPU SIMT simulator.

Five operations on the same four numbers: reduce, inclusive scan, exclusive scan, sort and select What CCCL computes — five operations on the same input input 3 1 4 1 reduce 9 one total of everything 3 + 1 + 4 + 1 inclusive scan 3 4 8 9 running total, up to and including me 3, 3+1, 3+1+4, … exclusive scan 0 3 4 8 running total of everything before me 0, 3, 3+1, … sort 1 1 3 4 the same values, in order select x > 1 3 4 keep only what passes a test, in order
Five words used everywhere below. A reduce gives one total; a scan (prefix sum) gives every running total — inclusive counts the item itself, exclusive stops just before it. Most of CCCL is these operations, done by thousands of threads at once.

Part 2 — One call, all the way down

The dive follows one line:

thrust::device_vector<float> v(1 << 28, 1.0f);   // 268,435,456 floats = 1 GiB on the GPU
float sum = thrust::reduce(v.begin(), v.end());  // 268,435,456 — exactly
The seven hops of thrust::reduce: iterator tag, CUDA backend, CUB size query, allocation, CUB launch, two GPU kernels, copy the result back Thrust CUB on the CPU Thrust CUB on the CPU GPU Thrust 1 thrust::reduce(v.begin(), v.end()) The iterator carries a tag: device_vector means the CUDA system. thrust/detail/reduce.inl 2 cuda_cub::reduce_n(policy, first, n, init, plus) Thrust's CUDA backend hands the whole job to CUB. thrust/system/cuda/detail/reduce.h 3 cub::DeviceReduce::Reduce(nullptr, bytes, …) Call 1 is a dry run: how many bytes of scratch memory? cub/device/device_reduce.cuh 4 temporary_array<uint8_t> tmp(policy, bytes + sizeof(T)) Thrust allocates the scratch space plus a slot for the answer. thrust/system/cuda/detail/reduce.h 5 cub::DeviceReduce::Reduce(tmp, bytes, …) Call 2: look up this GPU's tuning, size the grid, launch. cub/device/dispatch/dispatch_reduce.cuh 6 DeviceReduceKernel<<<G, 256>>> → DeviceReduceSingleTileKernel<<<1, 256>>> <<<G, 256>>> = G blocks of 256 threads; they reduce their shares, then 1 block adds the G sums. cub/device/dispatch/kernels/kernel_reduce.cuh 7 get_value(policy, result) Wait for the GPU, copy one float back to the CPU, return it. thrust/system/cuda/detail/reduce.h
Seven hops from one line to one number. The pattern in the middle is CUB's signature: it is called twice — first to ask how much scratch memory it needs, then to run — and in between the caller allocates (here, Thrust). CUB's newer overloads take an execution environment instead and allocate for you, by default with cudaMallocAsync.

The layers inside CUB

CUB's six nested layers: device-wide API, dispatch and policy, kernel and agent, block collective, warp collective, thread Device-wide API cub::DeviceReduce::Reduce the whole array · called from the CPU Dispatch + policy detail::reduce::dispatch pick tuning · size the grid · launch Kernel + agent DeviceReduceKernel · AgentReduce what one block does with its tiles Block collective BlockReduce 256 threads via shared memory Warp collective WarpReduce 32 lanes swap registers Thread ThreadReduce 16 items, added in a pairwise tree runs on the CPU runs on the GPU
CUB is layered like nesting dolls. The outer two run on the CPU; from the kernel inward everything runs on the GPU. The inner three are public API: a collective is one call that every thread of a block (or warp) makes together, so your own kernels can use exactly the pieces CUB builds its algorithms from.

Picking a tile shape

A tile is the chunk of the array one block handles at a time: threads × items per thread.

How CUB picks the tuning for a float sum: first by requested determinism, then by GPU generation DeviceReduce::Sum over floats — which tuning? cc = compute capability, the GPU's generation number (H100 = 9.0) gpu_to_gpu same bits on any GPU run_to_run (default) same bits on every run not_guaranteed order may change per run requested determinism Reproducible accumulator floats split into exponent bins, bins added in a fixed order cc ≥ 9.0: 224 threads × 13 items (Part 5 §2 shows why) Same tile + atomics blocks add their sums straight into the output with an atomic add (an add no other thread can interrupt): 1 kernel instead of 2; fastest, but the last bits may vary which GPU generation? cc ≥ 10.0 · Blackwell 512 threads × 16 items · vec 2 = 8,192 floats per tile else cc 6.0 – 9.x · e.g. H100 (9.0) 256 threads × 16 items · vec 4 = 4,096 floats per tile else older than cc 6.0 256 threads × 20 items Every leaf then goes through scale_mem_bound() to fit the element size — see Part 5 §3.
Tuning is data, not code. A policy selector maps (determinism, GPU generation, type, operator) to a tile shape. An H100 (compute capability 9.0) takes the default run_to_run branch and lands in its second box: 256 threads × 16 items (vector loads of 4 only when the input is a raw pointer). The switch itself is passed through CUB's newer overloads, e.g. cuda::execution::require(cuda::execution::determinism::gpu_to_gpu); why floats need it at all is Part 5 §2.

The numbers are measured, not guessed — tunings keep the benchmark verdict that chose them:

// ipt_17.tpb_192.ipv_2  1.128568  1.062323  1.127248  1.197567
return sm100_tuning_values{17, 192, 1 << 2};   // double sum, compute capability 10.7

Read it as: 17 items per thread, 192 threads, vector loads of 22 = 4 (ipv is a power of two). Against the default it replaced: at least 6 % faster, 13 % on average, up to 20 %.

Cutting the work into tiles

How a reduction is cut into tiles and spread over a grid of blocks, then combined by a second kernel one tile 256 × 16 = 4,096 floats tiles in 1 GiB of floats 268,435,456 ÷ 4,096 = 65,536 grid (blocks) 132 × 8 × 5 = 5,280 grid = SMs × blocks that fit on one SM × 5 (CUB's subscription factor); 8 = 2,048 threads per SM ÷ 256, if registers allow Illustrated with 4 blocks and 12 tiles — block b takes tiles b, b + 4, b + 8: tile 0 block 0 tile 1 block 1 tile 2 block 2 tile 3 block 3 tile 4 block 0 tile 5 block 1 tile 6 block 2 tile 7 block 3 tile 8 block 0 tile 9 block 1 tile 10 block 2 tile 11 block 3 + grid × tile sum 0 sum 1 sum 2 sum 3 kernel 2: one block adds them the sum kernel 1 ends with one partial sum per block — with 8 blocks per SM that is 5,280 blocks of 12–13 tiles each
Cutting one reduction up. A tile is the chunk a block handles in one go (threads × items per thread). The grid is capped at 5 × the blocks that fit on the GPU at once — the subscription factor — and each block strides through the array a tile at a time, so the input is read exactly once and every block ends with a single partial sum.

Inside one block

Inside one block: all of the block's tiles fold into 256 thread sums, then 8 warp sums, then 1 block sum the block's tiles: 4,096 floats each, one after another 256 per-thread sums 8 warp sums 1 block sum Thread level · registers per tile, each thread makes 16 striped loads and adds them into its running sum, pairwise Warp level · shuffles each warp adds its 32 sums in 5 shuffle steps (Part 5 §1); lane 0 ends up with the warp's total Block level · shared memory lane 0 of each warp writes its total; __syncthreads() (a barrier: the whole block waits); thread 0 adds the 8
One block, level by level. Threads first fold every tile the block owns into registers; only then, once per block, come the warp and block steps. Each level talks through the cheapest channel that reaches its partners: not at all, register-to-register (shuffles), then shared memory. Striped loads mean each warp instruction reads 128 neighbouring bytes; 16-byte float4 loads happen only when the input is a raw pointer, which Thrust's iterators currently are not.

Getting the data in

Reduce reads striped data straight into registers. Scan and sort need each thread’s items side by side, so they pay for a transpose:

Blocked versus striped arrangement of 16 items over 4 threads, and the warp transpose that gets both Blocked: thread t owns items 4t … 4t + 3 items 0 … 15 (4 bytes each), coloured by owning thread 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 T0 T1 T2 T3 load step 1 hits 4 addresses 16 bytes apart Striped: thread t owns t, t + 4, t + 8, t + 12 items 0 … 15 (4 bytes each), coloured by owning thread 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 T0 T1 T2 T3 load step 1 hits 4 neighbours: one transaction BLOCK_LOAD_WARP_TRANSPOSE gets both: coalesced reads, blocked data in registers read striped from global memory (coalesced) write to shared memory each warp stages its slice read back blocked each thread gets a contiguous run Blocked suits per-thread sequential work (a scan over my 16 items); striped suits the memory bus.
CUB calls these data arrangements. Striped is what memory wants — neighbouring threads touch neighbouring addresses — and blocked is what per-thread algorithms want. A sum doesn't care about order, so DeviceReduce loads striped straight into registers; scan and sort need blocked order, so they load striped and transpose through shared memory.

Part 3 — The other algorithms

A reduction only needs blocks to agree at the very end. Scan and sort need every block to know about the blocks before it — and that is where CCCL’s best ideas live.

Scan: decoupled look-back

A scan writes every running total (Part 1), so tile 7 cannot write a single output until it knows the sum of tiles 0–6.

Decoupled look-back: tile 7 sums the published aggregates of tiles 6, 5 and 4 and the inclusive prefix of tile 3 to get its own prefix X not ready yet A own sum known · PARTIAL P running total known · INCLUSIVE tile 0 P 4 tile 1 P 11 tile 2 P 13 tile 3 P 19 tile 4 A 3 tile 5 X → A 5 was still loading tile 6 A 1 tile 7 own sum 6 looking back lane 0 lane 1 · waited lane 2 lane 3 ignored: tile 3 already holds the total so far exclusive prefix of tile 7 = 1 + 5 + 3 + 19 = 28 then tile 7 publishes P = 28 + 6 = 34 so its successors can stop here What every tile does, in order 1 load its tile 2 sum it, publish A 3 look back to a P 4 publish P 5 scan, then store
Single-pass scan with decoupled look-back (Merrill & Garland). A tile publishes its own sum (A) the moment it has it, so a successor waits only for its neighbours' own sums, never for their full prefixes — that split is what "decoupled" means. It adds As until it meets a P; one warp checks 32 predecessors at a time. For 4-byte values, status and value share one 8-byte word, so a single store publishes both and no reader can see a flag without its value.

On Blackwell-class GPUs (compute capability 10.0 and up) a contiguous scan runs a newer, warp-specialized form of the same idea: separate groups of warps load, reduce, look ahead and store, pipelined through shared memory. Hopper uses it too when bit-exact float sums are requested; other GPUs and fancy iterators get the classic look-back above.

Radix sort: onesweep

Radix sort orders keys one 8-bit digit at a time, lowest digit first: count how many keys have each of the 256 digit values, turn the counts into start positions, copy every key to its slot. Part 5 §5 draws one pass.

Onesweep radix sort for 32-bit keys: one histogram kernel, one exclusive-sum kernel, then four onesweep passes that ping-pong between two buffers DeviceRadixSort, 32-bit integer keys (Hopper): 8-bit digits → 4 passes · tile = 384 × 20 keys Histogram kernel reads every key once counts all 4 digits at once → 4 × 256 counters 1 Exclusive-sum 4 blocks, one per digit counts → start offsets (Part 5 §5 shows why) 2 3 Onesweep kernel, once per digit a · claim the next tile number with atomicAdd b · rank = same-digit keys before it in the tile c · look back across tiles for each of 256 digits d · scatter: digit start + earlier tiles + rank The four passes, least-significant digit first — keys ping-pong between two buffers pass 1: bits 0 – 7 A → B pass 2: bits 8 – 15 B → A pass 3: bits 16 – 23 A → B pass 4: bits 24 – 31 B → A Traffic: 1 read for the histogram + 4 × (read + write) = 9 passes over the keys.
Onesweep (Adinets & Merrill): one kernel builds every histogram up front, and each digit pass is then a single kernel, because decoupled look-back — this time over 256 counters per tile — replaces the older three kernels per digit (count, scan, scatter). Tiles are numbered with atomicAdd in the order blocks actually start, so a block only ever waits on tiles that are already running: the look-back cannot deadlock.

Merge sort

Merge sort: each block sorts one tile, then ceil(log2(tiles)) merge passes double the sorted runs block sort sorted tile sorted tile sorted tile sorted tile sorted tile sorted tile sorted tile sorted tile merge pass 1 run of 2 tiles run of 2 tiles run of 2 tiles run of 2 tiles merge pass 2 run of 4 tiles run of 4 tiles merge pass 3 fully sorted Each merge pass = a partition kernel (merge-path binary searches: where does each block's output slice begin?) + a merge kernel. Tile = 256 threads × 17 keys = 4,352 keys; 268 M keys → 61,681 tiles → 16 merge passes.
DeviceMergeSort works for any comparison — strings, tuples, custom types — at the price of ⌈log2(tiles)⌉ passes over memory. Part 5 §5 shows the binary search that splits each merge evenly.

The whole family

CUB device algorithms grouped by the core technique each one is built on Tiles + second kernel DeviceReduce · Sum · Min · Max ArgMin · ArgMax · TransformReduce multi-pass radix sort (short keys, old GPUs) Decoupled look-back DeviceScan · ScanByKey DeviceSelect · If · Flagged · Unique DevicePartition · 2- and 3-way ReduceByKey · RunLengthEncode UniqueByKey · batched memcpy/copy Onesweep DeviceRadixSort look-back per digit Merge path DeviceMergeSort · DeviceMerge Privatized counters + atomics DeviceHistogram · Even · Range Pure streaming DeviceTransform · DeviceFor DeviceCopy::Copy (mdspan)
Six ideas cover almost everything. Decoupled look-back alone carries scan, select, partition, unique, reduce-by-key and run-length encoding — which is why Part 3 opened with it.

Part 4 — libcu++, Python and cudax

libcu++: the ground floor

What libcu++ contains: cuda::std for portable standard C++, cuda:: for GPU-only features, each with a one-line meaning cuda::std:: standard C++ that also compiles for the GPU atomic a counter many threads can update safely span · mdspan views of 1-D and N-D arrays, no copies complex · chrono · bit the usual numeric toolbox optional · tuple · expected vocabulary types algorithm · numeric sequential algorithms, inside one thread type_traits · limits compile-time facts about types cuda:: what only a GPU has atomic_ref<T, scope> an atomic, plus who must see it barrier · latch wait until a whole group arrives memcpy_async · pipeline copy the next tile while computing ptx:: one C++ function per GPU instruction stream · launch · buffer host-side plumbing, in modern C++ execution::require ask for guarantees, e.g. determinism
libcu++ is two libraries in one namespace family. cuda::std:: is the ordinary standard library, made to compile in GPU code; cuda:: adds what only a GPU has. CUB and Thrust are written in both.
Thread scopes from system down to thread, with the PTX qualifier each one compiles to thread_scope_system CPU + every GPU .sys thread_scope_device every thread on one GPU .gpu thread_scope_cluster a few blocks, Hopper+ .cluster thread_scope_block up to 1,024 threads .cta thread_scope_thread just me .cta no thread-only backend yet: libcu++ builds it as block scope C++ you write cuda::atomic_ref<int, cuda::thread_scope_device> a(x); int old = a.fetch_add(1, cuda::memory_order_relaxed); PTX it becomes atom.add.relaxed.gpu.s32 Atomic add = read, add, write as one step no other thread can cut into: 1,000 threads adding 1 always give exactly +1,000. Narrower scope = cheaper: only that group must see the update. .cta = PTX's word for a block.
A GPU has more than one meaning of "visible to others". libcu++ encodes it in the type — cuda::atomic<T, scope> — and each scope compiles to one PTX qualifier. The PTX shown is what libcu++'s inline-assembly backend writes; newer toolkits reach the same instruction through compiler builtins.
Timeline: plain loads serialise copy and compute; memcpy_async with a two-stage pipeline overlaps them Four tiles: copy (70 units) then math (50 units) each plain loads copy 0 math 0 copy 1 math 1 copy 2 math 2 copy 3 math 3 done at 480 memcpy_async + 2-stage pipeline copy 0 math 0 copy 1 math 1 copy 2 math 2 copy 3 math 3 done at 330 31 % sooner cuda::memcpy_async starts a copy and returns at once; pipeline.consumer_wait() blocks only when that stage's data is needed. On Hopper and later one bulk copy (TMA) moves a whole tile and signals a shared-memory mbarrier — CUB's DeviceTransform uses this.
Asynchronous copies turn "load, then compute" into a conveyor belt: while the math works on tile k, the copy for tile k + 1 is already in flight. With these numbers the belt finishes in 330 units instead of 480.

Python: the same kernels, compiled on demand

import cupy as cp, numpy as np
import cuda.compute

d_in, d_out = cp.array([1, 2, 3, 4, 5], dtype=np.int32), cp.empty(1, dtype=np.int32)
cuda.compute.reduce_into(d_in=d_in, d_out=d_out, num_items=5,
                         op=cuda.compute.OpKind.PLUS, h_init=np.array([0], dtype=np.int32))
# d_out[0] == 15
How a cuda.compute call is JIT-compiled: the Python operator and a CUB kernel template are compiled at runtime, linked, loaded and cached Python call cuda.compute.reduce_into( d_in=…, d_out=…, op=…, num_items=…, h_init=…) 1 Your operator a Python function? numba-cuda-mlir compiles it to LTO-IR (built-in OpKind.PLUS skips this) 2 Fill in the template DeviceReduceKernel< …your types…> as C++ source text 3 NVRTC compiles that source at runtime into LTO-IR 4 nvJitLink links kernel + operator → machine code for this GPU 5 Load and launch cuLibraryLoadData, then the same dispatch logic as C++ 6 Steps 2–5 run once per new combination of types and operator; the result is cached, so the next call goes straight to 6.
Python gets the same kernels as C++, compiled on demand. The C library (cccl.c.parallel) writes the CUB template instantiation as source text, NVRTC (NVIDIA's runtime compiler) compiles it, and nvJitLink fuses it with your operator. numba-cuda-mlir is a Python-to-GPU compiler; LTO-IR is half-compiled code that can still be optimised when it is linked.

cudax: tasks that find their own order

CUDASTF infers a task graph from declared data accesses: T1 writes X; T2 and T3 read X and run in parallel; T4 waits for both What you write — in plain program order T1 ctx.task(X.rw()) T2 ctx.task(X.read(), Y.rw()) T3 ctx.task(X.read(), Z.rw()) T4 ctx.task(Y.read(), Z.rw()) rw() = read and write read after write · write after read · write after write → an edge read after read → no edge, so T2 and T3 overlap What CUDASTF runs — a graph it built itself X X Y Z T1 T2 T3 T4 in parallel Executed on CUDA streams, or captured as one CUDA graph.
CUDASTF (in cudax) takes tasks in plain program order, reads what each one does to each piece of data, and derives the dependency graph from the classic hazards. It runs the graph on CUDA streams (in-order queues of GPU work) or hands it to the GPU in one piece as a CUDA graph.

Part 5 — The math, in simple words

Everything above rests on six ideas. Here they are one at a time, each drawn with its numbers.

1. Adding in parallel

Adding 8 numbers one at a time takes 7 rounds; as a tree the same 7 additions take 3 rounds One at a time: 7 additions, 7 rounds 1 3 + 1 = 4 2 4 + 4 = 8 3 8 + 1 = 9 4 9 + 5 = 14 5 14 + 9 = 23 6 23 + 2 = 25 7 25 + 6 = 31 As a tree: the same 7 additions, 3 rounds 3 1 4 1 5 9 2 6 4 5 14 8 9 22 31 round 1 4 at once round 2 2 at once round 3 1 addition
The whole trick of parallel reduction. Because (a + b) + c = a + (b + c), the additions can be regrouped into pairs that happen at the same time. Same work, far fewer rounds.
Rounds needed to add n numbers with a tree: 10 for a thousand, 20 for a million, 30 for a billion Rounds to add n numbers with a tree = ⌈log2 n⌉ 1 thousand 10 rounds one at a time: 999 rounds 1 million 20 rounds one at a time: 999,999 rounds 1 billion 30 rounds one at a time: 999,999,999 rounds Every round halves the pile: 1,000 → 500 → 250 → 125 → 63 → 32 → 16 → 8 → 4 → 2 → 1.
Doubling the data adds just one round, so the chain of steps that must wait for each other stays short even for a billion numbers. The time itself is set by memory (§3).
Warp shuffle reduction on 8 lanes: offsets 1, 2 and 4; lane 0 ends with the total 31 Each step, lane i adds the value held by lane i + offset — straight from its register lane 0 lane 1 lane 2 lane 3 lane 4 lane 5 lane 6 lane 7 start offset 1 offset 2 offset 4 3 1 4 1 5 9 2 6 4 5 5 6 14 11 8 6 9 11 19 17 22 17 8 6 31 28 27 23 22 17 8 6 A real warp has 32 lanes: offsets 1, 2, 4, 8, 16 → 5 steps.
shfl.sync.down lets a lane read another lane's register directly — no memory at all. Three steps collapse 8 values into lane 0; five collapse a whole warp. The grey lanes compute too, every step; their results are simply never used.

2. Floats round, so order matters

Floating-point rounding makes addition order-dependent, and the three determinism levels CUB offers float32 near 100 million can only store multiples of 8 99,999,984 99,999,992 100,000,000 100,000,008 100,000,016 100,000,001 does not exist: it rounds back Why 8: a float keeps 24 significant bits, so between 67,108,864 (226) and 134,217,728 (227) the gap between neighbours is 226 − 23 = 8. (100,000,000 + 1) − 100,000,000 = 0 (100,000,000 − 100,000,000) + 1 = 1 doubles: (0.1 + 0.2) + 0.3 = 0.6000000000000001 0.1 + (0.2 + 0.3) = 0.6 So CUB lets you choose how much the order may vary not_guaranteed atomics add block sums in arrival order fastest · last bits vary run_to_run DeviceReduce's default: one fixed tree for a given GPU and tuning → same bits on every run gpu_to_gpu RFA — reproducible floating- point accumulator: values split into exponent bins; sums inside a bin are exact, so order stops mattering: same bits on any GPU ← faster more reproducible →
Every float addition rounds, so a different grouping can give different last bits. DeviceReduce's default keeps the grouping fixed for a given GPU (DeviceScan's default is not_guaranteed); gpu_to_gpu uses a reproducible accumulator that gives the same bits everywhere, at some cost.

It bites harder than it looks: add 1.0f to a running total one at a time and it stops at 16,777,216 (224) — from there every + 1 rounds back down. The one-liner in Part 2 stays exact because the only sums that grow past 224 are whole numbers of 4,096-float tiles, and multiples of 4,096 are still exact in a float.

3. Memory is the speed limit

H100 roofline: reduce and scan sit far left of the ridge at 20 FLOP per byte, so memory bandwidth sets their speed H100: an algorithm can't run faster than the lower of the two roofs TFLOP/s 1/16 1/4 1 4 16 64 0.1 1 10 100 FLOP per byte of memory traffic (log scale) ridge: 67 ÷ 3.35 = 20 compute roof: 67 TFLOP/s (FP32) memory roof: 3.35 TB/s × FLOP per byte reduce 1 add per 4 bytes → ≤ 0.84 TFLOP/s scan 1 add per 8 bytes → ≤ 0.42 TFLOP/s big matrix multiplies live out here →
A reduction does one addition per 4-byte number it reads; the GPU could do 20 FLOPs in the time one byte arrives. So CUB's job is not faster math — it is keeping the memory system saturated.
Minimum time for 1 GiB of 32-bit data on an H100: reduce 0.32 ms, scan 0.64 ms, radix sort 2.88 ms Speed of light for 1 GiB of 32-bit values on an H100: time = bytes moved ÷ 3.35 TB/s reduce read once: 1 GiB 0.32 ms scan read + write: 2 GiB 0.64 ms radix sort, 32-bit keys 1 + 4 × (read + write) = 9 GiB 2.88 ms
"Speed of light" in CUB's vocabulary: the time the bytes alone need. A good implementation lands close to these numbers; everything in Part 2 exists to get there.
Memory coalescing: stride 1 needs 4 sectors per warp, stride 8 needs 32; efficiency falls to 12.5 percent stride 1 · coalesced — 32 lanes read words 0, 1, 2 … 31: 128 B = 4 sectors of 32 B sector 0 sector 1 sector 2 sector 3 stride 8 — lanes read words 0, 8, 16 …: each box is a whole 32 B sector with only 4 B used useful bytes ÷ bytes moved = 128 ÷ (32 × sectors touched) — stride 8: 128 ÷ 1,024 100 % stride 1 50 % stride 2 25 % stride 4 12.5 % stride 8 12.5 % stride 16 12.5 % stride 32
A warp's 32 loads are served in 32-byte sectors. Neighbouring addresses share sectors; spread them out and the same 128 useful bytes cost up to 8× the traffic. This is why CUB loads striped and uses vector loads.
Shared-memory bank conflicts: stride 8 puts 8 lanes on each of 4 banks; CUB's one-word padding spreads them over all 32 banks Shared memory has 32 banks; word i lives in bank i mod 32, and a bank serves one word per round No padding — lanes write words 0, 8, 16 … (8 items per thread): bar = lanes queued on a bank 8 rounds bank 0 8 16 24 31 CUB's padding — index + ⌊index ÷ 32⌋: every lane gets its own bank 1 round bank 0 8 16 24 31
Why BlockExchange — the transpose from "Getting the data in" — pads its shared memory: without padding, 8 lanes queue on one bank and the access takes 8 rounds; one spare word every 32 shifts each group onto fresh banks.
Little's law: bytes in flight equals bandwidth times latency; CUB's per-SM targets grow from 12 KiB on V100 to 128 KiB SM issues loads data comes back ~1 µs round trip under load — every dot is a request still in the pipe bytes in flight = bandwidth × latency H200, per SM: 36.4 bytes per ns × ~1,350 ns ≈ 48 KiB CUB's minimum bytes in flight per SM (DeviceTransform tuning table) V100 and older 12 KiB 11.2 B/ns per SM → covers 1.1 µs A100 80 GB 16 KiB 18.9 B/ns per SM → 0.87 µs Hopper · H100 / H200 48 KiB 25.4 / 36.4 B/ns → 1.9 / 1.35 µs Blackwell · B200 64 KiB 54.1 B/ns per SM → 1.21 µs cc 10.7 ("Rubin" in the source) 128 KiB Divide each target by that GPU's bandwidth per SM: 0.9 to 1.9 µs every time. The budget grows with bandwidth because latency doesn't shrink.
To stream at full bandwidth a GPU must keep this many bytes requested-but-not-yet-arrived. CUB turns it into a tuning knob: DeviceTransform grows items per thread until the target is met.
scale_mem_bound resizes the nominal tile so each thread keeps about 64 bytes in flight (32 bytes for 1-byte types) scale_mem_bound: one nominal tile (256 threads × 16 four-byte items), resized per element size items = clamp(16 × 4 ÷ size, 1, 32) threads = min(256, round_up(49,152 ÷ (size × items), 32)) → stays 256 here 49,152 bytes = 48 KiB, the shared memory every GPU guarantees one block (a different 48 KiB from Little's law) 32 items 1-byte = 32 B per thread 32 items 2-byte = 64 B per thread 16 items 4-byte = 64 B per thread 8 items 8-byte = 64 B per thread 4 items 16-byte = 64 B per thread
Tunings are written in nominal 4-byte units and scaled to the real element size, keeping each thread's bytes in flight near constant — the memory system cares about bytes, not items. 1-byte types are capped at twice the nominal item count.

4. Prefix sums

Kogge-Stone inclusive scan of 8 values in 3 steps with offsets 1, 2 and 4 Kogge–Stone scan (a warp's method): each step, every lane adds the lane offset places to its left lane 0 lane 1 lane 2 lane 3 lane 4 lane 5 lane 6 lane 7 input offset 1 offset 2 offset 4 3 1 4 1 5 9 2 6 3 4 5 5 6 14 11 8 3 4 8 9 11 19 17 22 3 4 8 9 14 23 25 31 Offset-k step: n − k lanes add → 7 + 6 + 4 = 17 additions instead of 7. Three steps instead of seven: speed bought with extra work.
A prefix sum ("scan") replaces each value with the total of everything up to it. Kogge–Stone does it in log2 n steps, which is ideal inside a warp where all 32 lanes run anyway.
Scanning 1,024 numbers: sequential 1,023 additions and rounds; Kogge-Stone 9,217 additions in 10 rounds; Blelloch 2,046 additions in 20 rounds Scanning 1,024 numbers: three textbook algorithms additions (work) rounds (time, given enough lanes) one at a time a plain loop 1,023 = n − 1 1,023 = n − 1 Kogge–Stone every lane, every step 9,217 = n·log₂n − (n − 1) 10 = log₂n Blelloch up a tree, then back down 2,046 = 2(n − 1) 20 = 2·log₂n CUB mixes them by level: one at a time inside a thread (least work), Kogge–Stone inside a warp (fewest rounds), warp-scans across a block's warps, and decoupled look-back across tiles.
No single scan wins on both axes. The cheap-work algorithm is serial; the fast one does 9× the work. CUB uses each where it is cheapest.
One CUB scan on an H100 uses four algorithms: a pairwise sum in each thread, Kogge-Stone in each warp, warp-scans across the block, decoupled look-back across tiles One DeviceScan on an H100 (floats: 128 threads × 24 items per tile) — each level does what is cheapest there Thread sums its own 24 items in registers about 2 adds per item overall · no talking Warp Kogge–Stone with shfl.sync.up 5 steps for 32 lanes · registers only Block last lane of each of 4 warps posts its total each warp adds the totals before it Grid decoupled look-back over tiles one pass · 2 bytes moved per byte …back each thread adds its warp, block and tile prefixes to its 24 running totals, then stores
The same pattern as the reduction funnel in Part 2, now for scan: serial where it is free, parallel where it is cheap, and a single global pass on top. On Blackwell-class GPUs the warp-specialized lookahead kernel replaces the bottom three rows' plumbing, with the same idea.
Stream compaction with a scan: flags for x greater than 5, exclusive scan gives each kept value its output slot Stream compaction — keep only x > 5, in order input 7 2 9 4 6 1 8 3 keep? x > 5 1 0 1 0 1 0 1 0 exclusive scan 0 1 1 2 2 3 3 4 output 7 9 6 8 slot 0 slot 1 slot 2 slot 3 Kept count = last scan value + last flag = 4 + 0 = 4. Every thread knows where to write without asking anyone.
Why scan matters so much: it turns "which items survive?" into "where does each survivor go?". DeviceSelect, DevicePartition and Unique are this picture plus decoupled look-back.
Memory traffic of three cross-tile scan strategies for 1 GiB: 3 GiB, 2 GiB but serial, and 2 GiB single pass Scanning 1 GiB across tiles: bytes moved, and the time at 3.35 TB/s reduce, then scan two kernels read read write 3 GiB · 0.96 ms chained scan one kernel, but tiles wait in a line read write 2 GiB · serial decoupled look-back one kernel, no waiting in line read write 2 GiB · 0.64 ms
Reduce-then-scan reads the input twice. Chained scan reads it once but makes every tile wait for its predecessor. Decoupled look-back gets the 2n traffic without the queue — 1.5× less data than the classic approach.

5. Sorting by digits, and merging by diagonals

One counting-sort pass: histogram of digits, exclusive scan gives start offsets, each key is scattered stably One radix pass = a counting sort on one digit (2-bit digits here, 8-bit in CUB) input keys (the small letter tells equal digits apart) 3a 1a 2a 0a 3b 1b 2b 1c ↑ after one pass — stable 0a 1a 1b 1c 2a 2b 3a 3b digit count starts at 0 1 0 1 3 1 2 2 4 3 2 6 starts at = all counts before it = an exclusive scan of the counts slot = starts at + same-digit keys to its left (1b → 1 + 1 = 2) Equal digits keep their order (1a, 1b, 1c). That stability lets each pass build on the previous one.
One pass of LSD radix sort is a counting sort: count each digit, scan the counts into starting positions, drop every key into its slot. Sorting 32-bit keys is four of these, 8 bits at a time.
Raw float bits sort negatives after positives and backwards; CUB's twiddle restores numeric order A float32 is 32 bits: s exponent · 8 bits fraction · 23 bits s = sign bit: 0 for +, 1 for − Raw float bits read as unsigned integers — the order is wrong small big +0 0000 0000 +1 3F80 0000 +2 4000 0000 −0 8000 0000 −1 BF80 0000 −2 C000 0000 After CUB's twiddle — numeric order small big −2 3FFF FFFF −1 407F FFFF −0 7FFF FFFF +0 8000 0000 +1 BF80 0000 +2 C000 0000 Twiddle: sign bit 0 → flip only the sign bit · sign bit 1 → flip every bit, so a bigger negative becomes a smaller integer. Then CUB maps −0.0 onto +0.0, so the two zeros sort as equals. Undo after sorting.
Radix sort only understands unsigned integers. IEEE floats store sign and magnitude, so CUB flips bits on the way in (and back on the way out) to make integer order equal numeric order.
Bytes moved per byte of 32-bit keys: radix sort 9 at every size; merge sort 18, 26, 34 and 38 as n grows Bytes moved per byte of 32-bit keys, if every pass streams through memory radix sort (onesweep): 1 + 4 × 2 = 9, whatever n is merge sort: 2 × (1 + ⌈log2 tiles⌉), tile = 4,352 keys 9 18 1 M keys 9 26 16.8 M 9 34 268 M 9 38 1.07 B
Radix sort's cost depends on the key's bits, not on n; merge sort pays one more pass per doubling. That is why Thrust sorts plain numbers with radix sort and keeps merge sort for custom comparisons; CUB offers both.
Merge path grid for A = 1 2 7 9 and B = 3 4 5 8; diagonals every 2 outputs give each of 4 threads exactly 2 outputs Merge A = 1 2 7 9 with B = 3 4 5 8: step down = take from A, step right = take from B 3 4 5 8 1 2 7 9 B → A ↓ d = 2 d = 4 d = 6 Walk: take the smaller head; ties go to A. Diagonal d holds the first d outputs. Binary search on it: compare A[i] with B[d−1−i]; the answer flips once — at the crossing. thread 0 → outputs 1, 2 from A, A thread 1 → outputs 3, 4 from B, B thread 2 → outputs 5, 7 from B, A thread 3 → outputs 8, 9 from B, A Each thread gets exactly 2 outputs however the values interleave — the work is split evenly before any merging starts.
Merge path (Odeh, Green et al.): merging is a walk from the top-left to the bottom-right corner. Cut the walk with evenly spaced diagonals and every thread gets the same amount of output; finding a cut costs one small binary search.

6. Histograms

HistogramEven with 4 bins over 0 to 100: 37 goes to bin 1, 99.9 to bin 3, 100 and minus 3 are dropped DeviceHistogram::HistogramEven — 4 bins over [0, 100) bin 0 bin 1 bin 2 bin 3 0 25 50 75 100 37 → bin 1 99.9 → bin 3 100 → dropped: the upper bound is exclusive −3 → dropped bin = ⌊(x − lower) × bins ÷ (upper − lower)⌋ = ⌊37 × 4 ÷ 100⌋ = ⌊1.48⌋ = 1 each block counts its samples into its own bins in shared memory then 1 atomic add per bin into the global result 268 M samples: instead of 268 M fights over 4 global counters, each block adds just 4 numbers.
A histogram is one line of arithmetic per sample. The engineering is in the counting: private per-block counters keep thousands of threads from fighting over the same global counter.

Cheat sheet

Which CCCL layer to use for which job Sort, reduce or scan a whole array — simply Thrust thrust::sort(v.begin(), v.end()) Same, but you own memory, streams, determinism CUB, device-wide cub::DeviceScan::ExclusiveSum(…) A fast block-wide sum or scan inside your kernel CUB, block and warp cub::BlockReduce<float, 256> Atomics, barriers, async copies inside a kernel libcu++ cuda::atomic_ref<int, scope> The same algorithms on Python arrays cuda.compute cuda.compute.reduce_into(…) Many dependent tasks, possibly on several GPUs cudax · CUDASTF ctx.task(x.read(), y.rw())
Pick the layer by how much control you need: Thrust for convenience, CUB when you manage memory and streams yourself or write your own kernels, libcu++ for the primitives underneath.

Sources

Questions? Let's Chat :)

gpu