100 DAYS OF GPU Course preview
Day 35Module 4
in-technical-review

CUDA radix sort, built from the scan you already wrote

CPU and GPU sorts often use different methods. A comparison sort keeps asking "is this key bigger than that one?", and each answer can split a warp across two branches.

Early GPU sorts often used bitonic sort. NVIDIA writes, "Previous GPU-based sorting routines have primarily used variants of bitonic sort" (GPU Gems 3, chapter 39, section 39.3.3, checked 2026-08-30, linked below). Bitonic sort over this page's four million keys, padded to a power of two, takes 276 rounds of compare-and-swap.

The sort on this page never compares two keys. Its four-bit form runs eight passes, each one a count, a scan and a scatter, and the scan is the one you built on day 31. By the end you can estimate a radix sort's cost from its digit width alone.

One digit at a time, and why the order has to survive

Radix sort never looks at a whole key. It takes a slice of the bits, the digit, moves each key into the bucket for its digit, then takes the next slice up and does it again. After the last one the array is sorted, and nothing ever asked how two keys compare.

That works only if every pass is stable: keys sharing a digit come out in the order they went in. Pass two orders by the second digit, and the only thing keeping the keys that tie there in first-digit order is that pass two did not disturb them. Lose stability once and every pass under it was wasted.

Start with the smallest digit there is, one bit. Splitting a list by one bit is not a sorting problem, it is day 33's compaction with the predicate fixed to a bit of the key. Write 1 under every key whose bit is 0, exclusive-scan that, and each of those keys is holding its own destination:

"In a temporary buffer in shared memory, we set a 1 for all false sort keys (b = 0) and a 0 for all true sort keys. We then scan this buffer. ... For a sort key at index i, this address is t = i - f + totalFalses."

Every lane works that out from its own index and its own bit. There is no divergence worth the name: the two arms of that select are two arithmetic expressions, not two code paths.

Widening to four bits changes the offset data. Two buckets need one number, totalFalses, to separate them. Sixteen need sixteen starting offsets, and each thread block needs its own sixteen or two blocks writing digit 7 collide.

The four-bit version builds a histogram per block, then runs an exclusive scan over the whole histogram. This is the same scan on a larger input.

Diagram: one stable split, built from one scan. A row of 16 key boxes across the top, then three bands of 16 cells under it, each cell aligned with the key above it. Band 1, the flag: 1 under every key whose bit is 0, 0 under the rest. Caption "16 keys, 9 with the bit clear." Band 2, the exclusive scan of the flag. Caption "the last cell plus the last flag is 9, the size of the clear group." Band 3, the destination: a clear key goes to its scan value, a set key to 9 plus its index minus its scan value. Caption "16 destinations, all distinct, and keys within each bit group keep their input order." Alt text: "Sixteen keys split into 9 bit-clear and 7 bit-set keys use one exclusive scan to produce 16 distinct, stable destinations."

The digit width looks free, and it is not

On a CPU, cache size often sets the digit width. A byte per pass gives four passes and a 256-entry counter array that fits in L1. This makes eight bits a natural first choice on a GPU, but the costs differ there.

Here is the sequential version, worth having in front of you because the middle line is why the scan days came first:

// counting sort on one 4-bit digit, sequential
for (size_t i = 0; i < n; ++i) { count[key[i] & 15]++; }
for (int d = 1; d < 16; ++d) { count[d] += count[d - 1]; }
for (size_t i = n; i-- > 0;) { out[--count[key[i] & 15]] = key[i]; }

Two things break in the move. The counter array stops being one array: every block keeps a private copy, which is privatization from day 29, so the scan runs over buckets times blocks entries rather than buckets. Widen the digit to eight bits and that array grows sixteenfold, past what one block can scan, and you are building day 32's multi-block scan before you can sort anything.

And the last line runs backwards one element at a time, because a sequential counting sort gets stability free from the loop order. A GPU has no loop order. It buys stability by making each tile's local rank exact, which for a four-bit digit means four one-bit splits in shared memory, one per bit.

Count, scan, scatter

Full program in code/day35-radix-sort/radix_sort.cu. Three kernels run per pass, and this page is about what they cost.

Both rows are the same code. The one-bit run and the four-bit run differ by one runtime argument, bits. Two kernels would have let something else change alongside the digit width.

The grid is fixed at 128 blocks, and that is not tuning. The histogram is buckets times blocks slots and one block scans all of them, so 16 times 128 has to land inside 2,048. A static_assert says so and names the fix: widening the digit turns the build red instead of producing a wrong answer.

The keys are padded, not guarded. 4,194,915 is not a multiple of the tile, so the host rounds up and fills the tail with 0xFFFFFFFF, the largest key there is, with one cudaMemset. A partly filled tile would still occupy slots in the block's local sort and push the real keys to the wrong offsets, so the guard belongs before the sort sees the data.

Everything is checked against std::sort at every element before anything is timed, exactly, no tolerance. Two arrays that match element for element, one of which is sorted, hold the same keys in the same order, so one comparison covers both sortedness and a lost key.

The split is the piece worth reading twice. With bits set to 1 it is day 33 unchanged:

        // `bits` stable one-bit splits, low bit of the digit first, leave the
        // tile ordered by the whole digit. With bits = 1 this loop runs once
        // and is day 33's compaction with the predicate fixed to a bit.
        for (int b = 0; b < bits; ++b) {
            const unsigned int key = sKeys[tid];
            const unsigned int bit =
                (key >> static_cast<unsigned int>(shift + b)) & 1u;

            // The flag is 1 for the keys that belong in the low half. The
            // scan's own barriers separate this read of sKeys from the write
            // below, so the permute needs no extra barrier in front of it.
            const unsigned int inc = blockScanInclusive(sScan, bit ^ 1u);
            const unsigned int lows = sScan[blockDim.x - 1u];

            // inc - 1 is the exclusive scan: how many low keys came first.
            // tid - inc counts the high keys before this one, because every
            // earlier thread is in exactly one of the two halves.
            const unsigned int dst =
                (bit == 0u) ? (inc - 1u) : (lows + tid - inc);
            sKeys[dst] = key;
            __syncthreads();
        }

The barrier at the bottom is the one that matters: the scan's own barriers already separate this iteration's read of sKeys from its write, and this one separates the write from the next iteration's read. After the loop the tile is in digit order, so the write out is one contiguous run per digit rather than one address per lane:

        const unsigned int d = (sKeys[tid] >> shift) & mask;
        atomicAdd(&sCount[d], 1u);

        // The tile is ordered by digit now, so a digit's run begins at the
        // one thread whose left neighbour holds a different digit. `||`
        // short-circuits, so thread 0 never reads sKeys[-1].
        if (tid == 0u || d != ((sKeys[tid - 1u] >> shift) & mask)) {
            sStart[d] = tid;
        }
        __syncthreads();

        const size_t dst = static_cast<size_t>(sBase[d]) + (tid - sStart[d]);
        out[dst] = sKeys[tid];
        __syncthreads();

Two limits matter. The read side is coalesced and the write side is not, because moving keys somewhere else is the job. And the offset scan does not pad its indices, so its later steps hit the worst bank conflict case day 15 measures; at 2,048 slots once per pass that is not where the time goes, and the file says so rather than pretending otherwise.

Results

Measured. Tesla T4, driver 595.84, CUDA 12.6 (V12.6.85), built with nvcc -O3 -arch=sm_75. Captured 2026-09-01 on the project's verification node; full transcript in the page's evidence file.

GPU: Tesla T4 (compute capability 7.5)
40 SMs, 49152 bytes of shared memory per block

4194915 keys, padded to 4227072, 16.1 MiB per buffer
128 blocks x 256 threads, each owns 33024 keys (129 tiles)

Both rows sorted the same keys and matched std::sort at all 4194915 elements.

kernel          passes   time (ms)   ms/pass   Mkeys/s   x copy
-------------   ------   ---------   -------   -------   ------
copy, no sort        -       0.167         -         -     1.00
1-bit split         32      20.805     0.650     201.6   124.95
4-bit digit          8       9.671     1.209     433.8    58.08

The copy row moves 16908288 bytes in and 16908288 out, 203.1 GB/s.
Each sort pass reads the keys to count them, reads them again to
scatter them and writes them once: 48.4 MiB per pass either way.

Four bits per pass beats one bit by roughly two to one, 9.671 ms against 20.805 ms, and the reason is in the pass count rather than the pass cost. A 4-bit digit needs 8 passes where a 1-bit split needs 32, and each pass is only about twice as expensive (1.209 ms against 0.650). Fewer, fatter passes win because every pass re-reads the whole array.

That is the entire tuning knob for a radix sort: wider digits mean fewer passes and a bigger histogram. Sixteen counters per block fit comfortably in shared memory; 256 for an 8-bit digit start to compete with everything else the block wants.

Both versions matched std::sort on all 4,194,915 keys. A sort that is merely fast is worthless, and a 1-bit split is easy to get subtly wrong at the tile boundaries, so the reference check is the first thing this program does.

The sort costs 58 times a straight copy of the same data. That sounds enormous until you notice it is doing 8 full passes, each reading the keys twice and writing them once.

CUDA 13 widened the four-bit version's lead. On the same T4, the one-bit sort took 24.218 ms and the four-bit sort took 9.586 ms, a 2.53x advantage instead of CUDA 12.6's 2.15x. Both still matched std::sort at all 4,194,915 keys. The four-bit time was effectively unchanged; the drift was in the 32-pass one-bit path.

Run it yourself

Use a CUDA GPU with compute capability 7.5 or later. From the repo's README:

nvcc -std=c++17 -O3 -arch=sm_75 -o radix_sort radix_sort.cu

Skip Compiler Explorer for this one. The program allocates 49 MiB on the device and sorts four million keys on the host for its reference.

The sorts take 20.805 plus 9.671 ms per repetition across thirteen repetitions on the measured system. The host-side std::sort also runs once per input, so the full program is a poor fit for a short shared runner.

Exercise

Make it sort key-value pairs: carry a second array of payloads through the sort so a payload lands wherever its key does. Seed each payload with its key's original index, then check that within every run of equal keys the payloads come out increasing.

Time: 25 to 40 minutes. Submit: the diff against radix_sort.cu, and one sentence on which of your changes the stability check is testing.

Check: the harness compares your keys against the reference exactly, as now, then adds the pair check: wherever out[i] == out[i - 1] it requires payload[i] > payload[i - 1], and it prints the first index that fails with both payloads and the key they share. It names which of the two went wrong, because right keys with wrong payloads is a different bug.

Hint 1

Every key already travels to a computed destination twice: once inside the tile and once out to global memory. What else has to make both trips, and what breaks if it only makes the second?

Hint 2

The tile permute writes sKeys[dst] = key. Look at what dst depends on and ask whether you could recompute it afterwards, or whether the payload has to ride along on every one of the four splits.

Solution

Add a second shared array beside sKeys, permute it with the same dst, and write it out with the same dst as the key. It has to ride along on all four splits: after split two the tile sits in an order neither the input nor the output index describes, so there is nothing left to recompute from.

The check is testing the split, not the scatter. Each split preserves the input order inside its two halves, and the block walks its tiles in order. Break either and the keys still come out sorted, which is why the payload check exists: a sort can be right and a radix pass still be broken, and only the ties can tell you.

Pitfalls

Your sort drops keys and still looks sorted. Two threads computing one destination means a key overwrites another, and the output still looks sorted because everything in it came out of a sorted pass. Compare against a reference at every element, not a prefix; the same failure and the same check are on day 33.

One digit sorts and the whole key does not. That is a pass that is not stable. Sorting on a single digit works fine without stability, so a one-digit test passes and hides it, and the bug only shows once a later pass has to preserve what an earlier one arranged.

You guard the partial tile instead of padding it. A thread that skips its write still holds a slot in the block's local sort, and its rank pushes every real key after it to the wrong offset. Pad to a whole number of tiles with a sentinel that sorts last, and let a stable sort put the padding where it belongs.

You time the sort in a loop and it gets faster every run. A sort that writes over its own input is sorted after run one, and a sorted input gives the scatter perfect locality. Sort out of place, from an input the loop never touches, or you publish the cost of re-sorting sorted data.

You reach for a radix sort on floats or signed integers. It orders keys by bit pattern, which is the order you want for unsigned integers and not for the other two: negatives come after positives, and negative floats come out backwards. The fix is a bijection on the way in and its inverse on the way out, and every production sort carries one.

Go deeper

Next

Day 36 takes the other route: merging two sorted runs, where how much work a thread owns depends on the data and merge path is how you find it without asking. Day 39 replaces these three kernels with one library call and asks what writing them bought you, a fairer question once you know what it cost.