Histograms are the engine behind one of the fastest ways to sort: Radix Sort, a noncomparison-based algorithm whose running time is O(N), performing a small number of linear-time passes over the data. Radix Sort builds directly on the counting sort described at the beginning of this chapter, which sorts N keys drawn from a small set of k possible values by histogramming the keys and using the bucket counts to decide where each key belongs in the output. Counting sort runs in O(N) time, but it is practical only when k is small; a 32-bit key would require a histogram of four billion buckets.
Radix Sort makes counting sort practical for wide keys by sorting on one digit at a time, where a digit is a group of b adjacent bits and so takes one of \(2^b\) values. The keys are sorted by their least significant digit, then by the next digit, and so on up to the most significant. The essential property is that each per-digit pass is stable—it preserves the relative order of keys that share the same digit—so sorting on a more significant digit never disturbs the ordering already established by the less significant digits. Once the most significant digit has been processed, the array is fully sorted.
Each pass is exactly the histogram-and-offset computation this chapter has been developing, in three steps:
Histogram the digit. For each key, extract the b-bit digit and increment the corresponding one of the \(2^b\) counters. This is the same operation as before; Listing 16-7 gives a kernel that computes it on the GPU with global memory atomics, and the privatized techniques of Section 16.3 can be used without modification.
Scan the counts. An exclusive prefix sum (Chapter 13) over the \(2^b\) counters turns each bucket’s count into its base offset: the index in the output array where that digit’s keys begin.
Scatter the keys. Walk the input once more, appending each key to its bucket and advancing that bucket’s offset. Because the walk proceeds in input order and each offset only increases, keys with equal digits land in the output in their original relative order—which is precisely what makes the pass stable.
__global__ voidRadixHistogram( int *pHist, const int *in, size_t N, int shift, int mask ){ for ( size_t i = blockIdx.x*blockDim.x + threadIdx.x; i < N; i += blockDim.x*gridDim.x ) { int digit = (in[i] & mask) >> shift; atomicAdd( &pHist[digit], 1 ); }}Listing 16-8 shows one full pass. The histogram and scan run over the \(2^b\) buckets, and the scatter loop places every key in a single stable sweep.
template<int b>voidRadixPass( int *out, const int *in, size_t N, int shift, int mask ){ const int numCounts = 1<<b; int counts[numCounts]; memset( counts, 0, sizeof(counts) ); // 1. Histogram: count the occurrences of each digit value. for ( size_t i = 0; i < N; i++ ) { int digit = (in[i] & mask) >> shift; counts[digit] += 1; } // 2. Exclusive scan: turn the counts into output base offsets. int sum = 0; for ( int i = 0; i < numCounts; i++ ) { int temp = counts[i]; counts[i] = sum; sum += temp; } // 3. Scatter: append each key to its bucket, in input order. for ( size_t i = 0; i < N; i++ ) { int digit = (in[i] & mask) >> shift; out[ counts[digit]++ ] = in[i]; }}The driver in Listing 16-9 applies one pass per digit, ping-ponging the data between two output buffers and returning whichever buffer holds the final, fully sorted result.
template<int b>int *RadixSort( int *out[2], const int *in, size_t N ){ int shift = 0; int mask = (1<<b)-1; int outIndex = 0; RadixPass<b>( out[outIndex], in, N, shift, mask ); while ( mask ) { outIndex = 1 - outIndex; shift += 1; mask <<= 1; RadixPass<b>( out[outIndex], out[1-outIndex], N, shift, mask ); } return out[outIndex];}Listings 16-8 and 16-9 keep the scan and scatter on the host, which
is enough to show the algorithm but leaves most of the work off the GPU.
The sample in radixSort.cu carries all three steps onto the GPU, sorting
one \(b\)-bit digit per pass in three
phases that ping-pong between two buffers:
Local sort. Each threadblock loads a tile of keys (one per thread) and stably sorts the tile by the current digit, using \(b\) successive one-bit splits. Each split is a block-wide exclusive scan that partitions the tile into the keys with a 0 and the keys with a 1 in the current bit (Listing 16-10). A by-product of the local sort is the tile’s per-digit histogram, accumulated in shared memory exactly as in Section 16.3—so the histogram phase uses no global atomics at all.
Global scan. The per-tile histograms form a matrix of \(2^b\) counts per tile. A single device-wide scan of that matrix (Chapter 13), laid out digit-major, hands every (tile, digit) pair the base offset of its keys in the output.
Scatter. Each block writes its already-locally-sorted keys to that base offset plus each key’s rank within its digit’s run in the tile (Listing 16-11). Because the tiles are visited in order and each tile is already stable, the global result is stable.
template<int b>__global__ voidRadixLocalSort( unsigned *sortedKeys, int *blockHist, const unsigned *in, size_t N, int shift, int numTiles ){ const int NUM_DIGITS = 1 << b; const unsigned mask = NUM_DIGITS - 1; const int tid = threadIdx.x, tile = blockIdx.x; const size_t base = (size_t) tile * RADIX_TILE; const int valid = (base + RADIX_TILE <= N) ? RADIX_TILE : (int)(N - base); __shared__ unsigned s[RADIX_TILE], sTmp[RADIX_TILE]; __shared__ int sScan[RADIX_TILE], sHist[1<<b]; // Pad the last tile with 0xFFFFFFFF so padding sorts to the end. s[tid] = (tid < valid) ? in[base+tid] : 0xFFFFFFFFu; __syncthreads(); // b one-bit stable splits leave the tile sorted by the b-bit digit. for ( int bit = 0; bit < b; bit++ ) { int flag = (s[tid] >> (shift+bit)) & 1; sScan[tid] = 1 - flag; // 1 marks a 0-bit ("false") __syncthreads(); inclusiveScanBlock( sScan ); int totalFalses = sScan[RADIX_TILE-1]; int f = sScan[tid] - (1 - flag); // exclusive falses before tid int dest = flag ? (totalFalses + tid - f) : f; __syncthreads(); sTmp[dest] = s[tid]; __syncthreads(); s[tid] = sTmp[tid]; __syncthreads(); } // Per-tile histogram in shared memory (Section 16.3): no global atomics. if ( tid < NUM_DIGITS ) sHist[tid] = 0; __syncthreads(); if ( tid < valid ) atomicAdd( &sHist[(s[tid] >> shift) & mask], 1 ); __syncthreads(); sortedKeys[base+tid] = s[tid]; if ( tid < NUM_DIGITS ) blockHist[tid*numTiles + tile] = sHist[tid];}template<int b>__global__ voidRadixScatter( unsigned *out, const unsigned *sortedKeys, const int *scanIncl, const int *blockHist, size_t N, int shift, int numTiles ){ const int NUM_DIGITS = 1 << b; const unsigned mask = NUM_DIGITS - 1; const int tid = threadIdx.x, tile = blockIdx.x; const size_t base = (size_t) tile * RADIX_TILE; const int valid = (base + RADIX_TILE <= N) ? RADIX_TILE : (int)(N - base); __shared__ int digitStart[1<<b]; unsigned key = sortedKeys[base+tid]; int d = (key >> shift) & mask; // Find the first position of each digit's run within the sorted tile. if ( tid < valid ) { int prevd = (tid > 0) ? ((sortedKeys[base+tid-1] >> shift) & mask) : -1; if ( d != prevd ) digitStart[d] = tid; } __syncthreads(); if ( tid < valid ) { int gbase = scanIncl[d*numTiles+tile] - blockHist[d*numTiles+tile]; out[gbase + (tid - digitStart[d])] = key; }}The entire sort runs without returning to the host between passes. On
an RTX 3060, the eight-pass (\(b\)=4)
configuration sorts 16 million 32-bit keys at about 550 million keys per
second—roughly 24x higher throughput than std::sort on a
single core of the same machine’s AMD Ryzen 7 7700X (about 23 million
keys per second). Among the digit widths the design allows (\(b\) must be a power of 2 less than or equal
to \(\lg k = 32\), and a tile’s \(2^b\)-entry histogram must fit in shared
memory, which leaves 1, 2, 4, and 8), \(b\)=4 is the fastest at this size: the
four-pass \(b\)=8 does half as many
passes, but must scan a histogram matrix that is 16x larger, and
overtakes \(b\)=4 only when the array
shrinks enough to be launch-bound rather than kernel-bound. Production
libraries such as Thrust and CUB are faster still—they process many keys
per thread and fuse these phases—but the microdemo makes the essential
point: a complete, fully GPU-resident sort is assembled almost entirely
from the two primitives this book has already built, the histogram and
the scan.
Because a given sort is a fixed sequence of several dozen
kernel launches—the same ones every time, however the passes and the
scan’s recursion unfold—it is a natural candidate for a CUDA graph
(Section 6.9). Capturing the sort once and replaying it with a single
cudaGraphLaunch() removes the CPU cost of reissuing every
kernel, and whether that runs faster depends entirely on problem size.
At 16 million keys, the kernels run for tens of milliseconds and the
launch overhead is lost against them, so the graph is a wash. At a few
thousand keys, each kernel instead finishes in microseconds, the
per-launch dispatch cost comes to dominate, and collapsing the whole
sequence into one launch cuts the sort’s time by about a third; on this
GPU, the crossover sits near 64 thousand keys, and the
radixSort.cu sample measures both paths across that range.
The graph does not make any kernel run faster—it removes work from the
CPU, not the device—so the benefit appears when a small, fixed launch
sequence repeats often enough that the CPU can no longer keep the GPU
fed: a batched sort of many short arrays, or a sort buried in the inner
loop of a larger computation. The digit width addresses the same cost
from the other direction: at small sizes, a wider radix and a captured
graph both speed up the sort by reducing the number of kernel
launches.