When to stop hand-writing: Thrust, CUB and libcu++
Across five days this course built three algorithms by hand: a sum on days 24 and 25, a scan on days 31 and 32, and a radix sort on day 35. Here is the sum again.
float sum = thrust::reduce(thrust::device, d_in, d_in + n, 0.0f);
No -l flag and no package to install. nvcc puts those
headers on your include path itself, and the scan and the sort are one line
each too.
You may conclude that those five days were wasted. This page tests that claim by running your kernels and the library versions in one process. It compares their times and line counts, then says which kernels to keep.
This page also introduces libcu++, whose types return on days 74 and 75.
One repository, three libraries, already installed
CCCL is the CUDA C++ Core Libraries: Thrust, CUB and libcu++, developed in one repository and versioned together. It ships inside the toolkit, so unless you deliberately put a newer copy on the include path, the toolkit sets the version you compile with.
CUDA 12.6 carries CCCL 2.5.0; the mapping for every other release is the "Mapping to CTK Versions" table at https://github.com/NVIDIA/cccl (checked 2026-08-30). The program below prints the versions it was built against rather than asserting them.
Thrust is the top layer, shaped like the C++ standard
library. You hand it a range and it runs the whole algorithm, launches and
all. Calls include thrust::reduce, thrust::sort, and
thrust::inclusive_scan; you never set a block size.
CUB provides operations at three scopes. Device scope
(cub::DeviceScan, cub::DeviceRadixSort) is a whole
algorithm like Thrust's, but it takes raw pointers and its own scratch
buffer, and it returns a cudaError_t instead of a value. Block scope
(cub::BlockReduce, cub::BlockScan) and warp scope (cub::WarpScan) are
different: they go inside a kernel you still write.
Day 23's shuffle reduction and day 24's shared memory reduction each become one line inside your own loop body.
libcu++ is the standard library that works on both sides of the launch.
cuda::std:: is the parts of std:: that make sense in device code, and on
top of that sit the CUDA-shaped primitives: cuda::atomic_ref,
cuda::barrier and cuda::pipeline.
Diagram: the three CCCL layers and where each one is called. Three stacked bands, host code on the left of each and device code on the right, with a vertical line marking the launch boundary. Band 1, Thrust: one box on the host side,
thrust::reduce, with an arrow crossing the line to a greyed-out box labelled "kernels you never see". Caption "One call, host side, launches and all." Band 2, CUB device scope: two host boxes, the sizing call and the working call, sharing one scratch buffer, arrow crossing to the same greyed box. Caption "Two calls, because you own the scratch." Band 3, CUB block scope plus libcu++: nothing on the host side and three boxes on the device side inside one kernel you wrote,cub::BlockReduce,cuda::atomic_refand your own code between them. Caption "Inside your kernel, where a whole-array call cannot go." Alt text: "Three ways to call CCCL. Thrust is one host call. CUB at device scope is two, because you own the scratch buffer. CUB at block scope and libcu++ sit inside a kernel you write yourself."
Two beliefs that both send you the wrong way
CUB is not slow merely because it is a library. It has a tuning policy for
each supported architecture, and the compiler selects one from the target
passed to nvcc. The generated code can therefore differ across GPUs.
Using a whole-array library call for every operation also causes problems. It
cannot run inside your kernel, so a reduction within a loop must use
cub::BlockReduce rather than cub::DeviceReduce.
A whole-array call also cannot fuse with nearby work. Replacing three kernels with three library calls can add two trips through global memory. Day 48 measures that cost.
Days 24 and 25 did not try to beat CUB. They taught you to read a bandwidth number, why removing a branch can cost you more than it saves, and what launch overhead looks like when you halve the grid. That knowledge still applies when you use a library kernel.
Three algorithms, twice, in one process
Full program in code/day39-cccl/cccl.cu.
It runs each algorithm the hand-written way and the library way on the same
input in the same process, checks every result, and prints seven timed rows.
Four rules make the comparison valid.
No row writes over its own input. The tenth timed run sees the same data as the first. Hand a radix sort back its own sorted output and you have measured something else, because a sorted input turns the scatter into a copy.
Temp storage is queried and allocated outside every timed region, for the hand-written rows as well as the library ones.
The reduce rows count input bytes only. Two of the three are closed libraries and this program cannot count their internal passes, so the one figure all three share is the 64 MiB every one of them must read.
The sort rows carry no bandwidth column at all. How many bytes a sort moves is a property of the algorithm, so that column would divide two different quantities and print the answer as a ratio.
The library half of the program is this, all of it:
// thrust::device is what tells Thrust these are device pointers. Without it,
// Thrust treats a raw pointer as host memory and dereferences device memory
// on the CPU. Note also that this call returns a value to the host, so it
// synchronizes and it allocates its own scratch; the three CUB calls below do
// neither, which is a real difference in what you are buying and not a
// measurement artifact.
static float reduceThrust(const float* d_in, size_t n) {
return thrust::reduce(thrust::device, d_in, d_in + n, 0.0f);
}
// Every CUB device entry point is called twice. With d_temp null it writes
// the scratch it needs into tempBytes and does no work; with a real pointer
// it runs. Forgetting the second call is silent: no error, no output.
//
// They return cudaError_t, so the course's own CUDA_CHECK wraps them
// unchanged. tempBytes is taken by reference because CUB writes to it.
static void reduceCub(void* d_temp, size_t& tempBytes, const float* d_in,
float* d_out, size_t n) {
CUDA_CHECK(cub::DeviceReduce::Sum(d_temp, tempBytes, d_in, d_out,
static_cast<int>(n)));
}
static void scanCub(void* d_temp, size_t& tempBytes, const unsigned int* d_in,
unsigned int* d_out, size_t n) {
CUDA_CHECK(cub::DeviceScan::ExclusiveSum(d_temp, tempBytes, d_in, d_out,
static_cast<int>(n)));
}
static void sortCub(void* d_temp, size_t& tempBytes, const unsigned int* d_in,
unsigned int* d_out, size_t n) {
CUDA_CHECK(cub::DeviceRadixSort::SortKeys(d_temp, tempBytes, d_in, d_out,
static_cast<int>(n)));
}
And here is the other half of CCCL, the half you call from inside a kernel of
your own. cub::BlockReduce does the tree; cuda::atomic_ref is the C++20
atomic_ref with a thread scope attached, which is the promise you are
making about who else touches that address. Day 26 wrote this with a plain
atomicAdd; day 27 measured the scopes.
__global__ void countOddKeys(const unsigned int* __restrict__ keys,
int* __restrict__ total, size_t n) {
using BlockReduce = cub::BlockReduce<int, kThreadsPerBlock>;
__shared__ typename BlockReduce::TempStorage temp;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const int mine = (i < n) ? static_cast<int>(keys[i] & 1u) : 0;
const int blockSum = BlockReduce(temp).Sum(mine);
if (threadIdx.x == 0) {
cuda::atomic_ref<int, cuda::thread_scope_device> counter(*total);
counter.fetch_add(blockSum, cuda::memory_order_relaxed);
}
}
The source size is counted rather than guessed. These are non-blank, non-comment lines inside the file's marked regions, and the README gives the one-line command that reproduces them without a GPU.
| Algorithm | Hand-written | Library | Ratio |
|---|---|---|---|
| reduce | 40 | 3 (Thrust), 5 (CUB) | 8 to 13 |
| scan | 77 | 5 | 15 |
| sort | 42, on top of the scan's 77 | 5 | 24 |
| all three | 159 | 18 | 9 |
The 18 counts both reduce entry points, Thrust's and CUB's. Doing each algorithm once is 13 lines with Thrust's reduce or 15 with CUB's, so the one-of-each ratio is 12, not 9.
Note. The hand-written reduction in this program is day 24's version 4, the ladder's fourth rung, not day 25's eighth. Restating four more rungs here would double the file to make one row look better, and porting version 8 into it is this page's exercise. Read the reduce row against day 25's own table on the same card before drawing a conclusion about libraries.
Results
Re-verified on a Tesla T4 with driver 580.173.02 and CUDA 13.0 (V13.0.88) on 2026-09-02. CUDA's bundled Thrust and CUB moved from 2.5.0 to 3.0.1. All correctness gates and library winners reproduced, but slower hand-written reduce and scan timings increased the measured library speedups.
The original
CUDA 12.6 table below and CUDA 13 transcript both remain in evidence.
GPU: Tesla T4 (compute capability 7.5), 40 SMs
CUDA runtime 12.6 ships Thrust 2.5.0 and CUB 2.5.0
Key width from cuda::std::numeric_limits: 32 bits
CUB temp storage: reduce 3583 B, scan 70655 B, sort 70029567 B
Part 1: reduce 16776605 floats to one, 3 hand passes
implementation time (ms) GB/s vs hand
-------------------------------- ---------- ---------- --------
hand, day 24 version 4 0.349 192.2 1.00
thrust::reduce 0.277 242.2 1.26
cub::DeviceReduce::Sum 0.254 264.4 1.38
Part 2: exclusive scan of 16776605 uints, 3 levels
implementation time (ms) GB/s vs hand
-------------------------------- ---------- ---------- --------
hand, two-level Hillis-Steele 1.331 100.8 1.00
cub::DeviceScan::ExclusiveSum 0.571 235.1 2.33
Part 3: sort 16776605 distinct uint keys, 32 hand passes
implementation time (ms) vs hand
-------------------------------- ---------- --------
hand, one bit per pass 92.095 1.00
cub::DeviceRadixSort::SortKeys 4.484 20.54
No GB/s column here. How many bytes a sort moves is a property of
the algorithm, so the two rows would not be counting the same thing.
Part 4: cub::BlockReduce plus cuda::atomic_ref
odd keys counted on the device: 8388303, host says 8388303
The libraries win, and this page exists to say so plainly. Under CUDA 12.6, CUB's reduction was 1.38x faster than the hand-written day 24 version, its scan was 2.33x faster and its sort was 20.54x faster. CUDA 13.0 reproduced the winners at 2.44x, 3.29x and 20.98x.
The larger reduction and scan gaps came from slower hand-written rows while the library rows stayed close, so they are measured version-specific margins, not new algorithmic claims. The large sort gap remains mostly algorithmic: the hand sort moves one bit per pass and CUB's radix sort moves several.
The reductions in days 24 and 25 and the scans in days 31 and 32 teach what
those kernels do and why faster versions use their chosen structures. After
building them, you can read CUB's documentation and understand what
DeviceScan::ExclusiveSum does to your data.
The temp-storage line affects program design. CUB 2.5.0 asks for 3,583 bytes to reduce, 70,655 to scan and 70,029,567 to sort. CUB 3.0.1 kept reduce and sort unchanged and raised scan scratch by 256 bytes to 70,911.
A radix sort needs 70 MB of scratch for a 64 MiB input, about as much as the data itself. This can change a program's memory plan.
Use the library. Write the kernel when you have measured that the library is the bottleneck, when your operator does not fit its interface, or when you are fusing it into something bigger to avoid a round trip to memory, which is the case day 48 makes.
Run it yourself
Use a CUDA GPU with compute capability 7.5 or later. The build line is in the repo's README, and there is nothing to link:
nvcc -std=c++17 -O3 -arch=sm_75 -o cccl cccl.cu
The full program allocates about 450 MiB plus CUB's temp storage, and its
hand-written sort makes 32 passes over 64 MiB thirteen times, which is past
what Compiler Explorer's 20 second run cap will take. The 20 second
compile cap is the tighter one here, and it follows from the same design
that removes the -l flag: CCCL is header only, so a single
#include <cub/cub.cuh> costs real compiler time. An embed on this page has
to be one algorithm at a few million elements, not three at sixteen
million.
Exercise
Replace the hand-written scan inside the radix sort with
cub::DeviceScan::ExclusiveSum and time the sort again. Report the two sort
times and say, in one sentence, how much of the gap to
cub::DeviceRadixSort::SortKeys that closed.
Time: 25 to 40 minutes. Submit: the hand sort's time before and after, the fraction of the gap closed, and the sentence.
Check: the harness runs both sorts against std::sort on all 16,776,605
keys and prints the first index where either disagrees, so a swap that broke
the scatter fails instead of getting faster. It prints all three sort times
whether it passes or not, because the times are the answer and the pass is
only the gate.
Hint 1
The sort calls the scan once per bit. Before you measure anything, work out how many times the scan runs in one sort, and what fraction of the sort's work is the two kernels either side of it.
Hint 2
Compare two numbers you already have: how much faster is the library scan than yours in part 2, and how much faster is the library sort than yours in part 3? If the second is much larger than the first, the scan is not the thing you are behind on.
Solution
It closes only part of the gap. A faster scan makes each pass cheaper, but
your sort still walks the keys 32 times because it sorts one bit per pass.
cub::DeviceRadixSort sorts several bits per pass, so it needs fewer passes
and a wider histogram.
A faster scan cannot remove the hand-written sort's extra passes.
Pitfalls
Your CUB call does nothing and reports success. Every CUB device entry point is called twice: once with a null temp-storage pointer, which writes the byte count it needs and does no work, and once with a real buffer. Miss the second call and the output keeps whatever was in it, with no error anywhere.
You pass a device pointer to a Thrust algorithm with no execution policy.
Thrust's documentation is explicit: "Like the STL, Thrust permits this usage
and it will dispatch the host path of the algorithm"
(https://docs.nvidia.com/cuda/archive/12.2.0/thrust/index.html , checked
2026-08-30), and the host path then dereferences device memory on the CPU.
Pass thrust::device first, or wrap the pointer in a thrust::device_ptr.
You compare thrust::reduce with cub::DeviceReduce::Sum as if they were
the same call. Thrust's returns a value to the host, so it synchronizes and
allocates its own scratch every call. CUB's writes to a device pointer and
returns immediately. Which you want depends on whether the answer is needed
on the CPU.
You expect a floating-point scan to give the same answer twice. CUB says otherwise in the header: "Results are not deterministic for pseudo-associative operators (e.g., addition of floating-point types). Results for pseudo-associative operators may vary from run to run" (https://github.com/NVIDIA/cccl/blob/v2.5.0/cub/cub/device/device_scan.cuh , checked 2026-08-30). The scan here runs on integers for that reason.
Day 68 is the lesson on reproducibility.
You pin a CCCL older than the one your toolkit ships. The support rule is one-directional: a newer CCCL works with an older toolkit, and "CCCL is never forward compatible with the CUDA Toolkit" (https://github.com/NVIDIA/cccl , checked 2026-08-30).
You reach for cub::DeviceReduce inside a kernel. Device-scope entry
points are host functions that launch kernels. Inside a kernel the answer is
cub::BlockReduce or cub::WarpReduce, the layer most people never find.
Go deeper
- NVIDIA/CCCL, and its "Mapping to CTK Versions" table: https://github.com/NVIDIA/cccl
- Thrust documentation: https://nvidia.github.io/cccl/unstable/thrust/
- CUB documentation, block and warp scopes as well as device: https://nvidia.github.io/cccl/unstable/cub/
- libcu++ documentation, for
cuda::atomic_ref,cuda::barrierandcuda::pipeline: https://nvidia.github.io/cccl/unstable/libcudacxx/ - CUDA Programming Guide 4.9, "Asynchronous Barriers" (https://docs.nvidia.com/cuda/cuda-programming-guide/04-special-topics/async-barriers.html) and 4.10, "Pipelines" (https://docs.nvidia.com/cuda/cuda-programming-guide/04-special-topics/pipelines.html), both checked 2026-08-30
cuda-samples,cpp/2_Concepts_and_Techniques/radixSortThrust: https://github.com/NVIDIA/cuda-samples/tree/master/cpp/2_Concepts_and_Techniques/radixSortThrust
Next
Day 40 is capstone 2, PageRank on a graph. It uses your own CSR kernels first,
then a Thrust version on the same data. The two types introduced here return
on day 74, where
cuda::barrier replaces __syncthreads() because it can be waited on in
phases, and day 75, where cuda::pipeline drives
cp.async.
The hardware-accelerated form of both needs compute capability 8.0.