code/day35-radix-sort/radix_sort.cuThis is the source used by the lesson and its recorded evidence. Compile commands and expected output live in the directory README.
// SPDX-License-Identifier: MIT
//
// Day 35: a radix sort built from a scan.
//
// Sorts the same 4,194,915 unsigned 32-bit keys twice: once with a one-bit
// digit, which takes 32 passes, and once with a four-bit digit, which takes
// 8. Both runs use the same three kernels per pass (count, scan, scatter) and
// the same stable one-bit split. Only the digit width moves.
//
// Nothing here compares two keys. Every branch a thread takes depends on its
// own index or on a bit it already holds, never on how two keys order.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o radix_sort radix_sort.cu
// Run: ./radix_sort
//
// Verified 2026-08-30 on a Tesla T4 (compute capability 7.5), driver
// 595.84, CUDA 12.6 (V12.6.85). Transcript: evidence/run-2026-08-30.txt
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cuda_runtime.h>
// The one error macro. This file is standalone, the way a Compiler Explorer
// embed is, so it carries its own verbatim copy. `err_` carries a trailing
// underscore so it cannot collide with a variable at the call site, and the
// do/while makes the macro one statement so it survives a braceless `if`.
#define CUDA_CHECK(call) \
do { \
cudaError_t err_ = (call); \
if (err_ != cudaSuccess) { \
std::fprintf(stderr, "CUDA error %s:%d: %s: %s\n", __FILE__, \
__LINE__, #call, cudaGetErrorString(err_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
// 4194304 is 2^22 and 611 is 13 x 47, so the key count is not a multiple of
// any tile and the padding path below runs on every launch.
constexpr size_t kElems = 4194304ull + 611ull;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kKeyBits = 32;
constexpr int kBitsPerPass = 4;
constexpr int kBuckets = 1 << kBitsPerPass; // 16 digit values
constexpr int kNarrowBits = 1; // the one-bit run
// The grid is fixed, not derived from the key count: each block owns one
// contiguous span and walks it tile by tile. 128 is not a tuning choice. The
// digit histogram is kBuckets * kBlocks slots and scanOffsets scans the whole
// thing in one block, so 16 x 128 has to land inside kHistSlots.
constexpr int kBlocks = 128;
constexpr int kScanThreads = 1024;
constexpr int kHistSlots = 2048;
// Every block owns kSpan keys and every span is a whole number of tiles, so
// no kernel below carries a bounds check. The guard lives in the host's
// padding instead, and it has to: a partially filled tile would still occupy
// slots in the block's local sort, and its ranks would push the real keys to
// the wrong global offsets.
constexpr size_t kChunk = static_cast<size_t>(kBlocks) * kThreadsPerBlock;
constexpr size_t kPadded = ((kElems + kChunk - 1) / kChunk) * kChunk;
constexpr size_t kSpan = kPadded / kBlocks;
constexpr unsigned int kSeed = 20260830u;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr int kVariants = 2;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert(kKeyBits % kBitsPerPass == 0,
"the digit width has to divide 32, or the last pass reads bits "
"the key does not have");
static_assert(kKeyBits % kNarrowBits == 0, "same rule for the one-bit run");
static_assert(kBuckets * kBlocks <= kHistSlots,
"the digit histogram is kBuckets * kBlocks slots and scanOffsets "
"scans all of them in one block. Widening the digit past 4 bits "
"needs day 32's multi-block scan, or a smaller grid");
static_assert(kHistSlots == 2 * kScanThreads,
"scanOffsets gives every thread exactly two slots");
static_assert((kHistSlots & (kHistSlots - 1)) == 0,
"the Blelloch scan halves and doubles its stride, so its length "
"must be a power of two");
static_assert(kScanThreads <= 1024,
"no compute capability in this course allows more than 1024 "
"threads in one block");
// Inclusive scan of one value per thread, in shared memory, Hillis-Steele.
//
// This is day 31's kernel cut down to what a split needs, made a device
// function so the scatter can call it once per bit. It costs n log2(n) adds
// against a sequential scan's n; at one tile per call that trade is invisible
// beside the global traffic the pass moves.
//
// Every thread in the block must call it, because it contains barriers. On
// return `s` holds the whole inclusive scan, so a caller that needs the block
// total reads s[blockDim.x - 1] before anything writes to `s` again.
__device__ unsigned int blockScanInclusive(unsigned int* s, unsigned int v) {
const unsigned int tid = threadIdx.x;
s[tid] = v;
__syncthreads();
for (unsigned int offset = 1u; offset < blockDim.x; offset <<= 1) {
// Read into a register, barrier, write, barrier. The guard covers the
// read and never the barrier, so a thread with nothing to read still
// arrives at both.
unsigned int add = 0u;
if (tid >= offset) {
add = s[tid - offset];
}
__syncthreads();
s[tid] += add;
__syncthreads();
}
return s[tid];
}
// Counts, for one block's span, how many keys carry each digit value.
//
// One thread walks the span with a stride of blockDim.x, so a warp's 32
// addresses are 32 consecutive keys on every iteration and the read is
// coalesced.
//
// The output layout is digit-major, hist[d * gridDim.x + blockIdx.x]. That is
// what lets a single exclusive scan of the whole array give, in each slot,
// the global position of that (digit, block) pair's first key. A block-major
// layout would need one scan per digit instead.
//
// Launch assumption: gridDim.x blocks, each owning `span` keys, and `span` a
// whole number of blockDim.x tiles.
__global__ void countDigits(const unsigned int* __restrict__ in,
unsigned int* __restrict__ hist, size_t span,
int shift, unsigned int mask) {
__shared__ unsigned int sHist[kBuckets];
const unsigned int tid = threadIdx.x;
const unsigned int buckets = mask + 1u;
for (unsigned int b = tid; b < buckets; b += blockDim.x) {
sHist[b] = 0u;
}
__syncthreads();
const size_t blockStart = blockIdx.x * span;
const size_t step = blockDim.x;
for (size_t t = tid; t < span; t += step) {
const unsigned int d = (in[blockStart + t] >> shift) & mask;
atomicAdd(&sHist[d], 1u);
}
__syncthreads();
for (unsigned int b = tid; b < buckets; b += blockDim.x) {
hist[b * gridDim.x + blockIdx.x] = sHist[b];
}
}
// Exclusive scan of the whole digit histogram, in one block.
//
// Day 32's Blelloch scan: an upsweep that builds partial sums in place, a
// cleared root, and a downsweep that turns them into exclusive prefixes. It
// costs 2n adds against Hillis-Steele's n log2(n), which is why day 32 ends
// on it.
//
// One block is enough because the array is kBuckets * kBlocks slots, not one
// slot per key. Collapsing 4 million keys into 2048 counters first is the
// whole reason the histogram exists.
//
// Launch assumption: exactly one block of kScanThreads threads, with
// kHistSlots equal to twice that and a power of two.
//
// Ceiling: the indices are not padded against bank conflicts, so the later
// upsweep steps hit day 15's worst case. At 2048 slots once per pass that
// does not show up beside the scatter. GPU Gems 3 chapter 39 gives the padded
// index if a wider digit ever makes this array large enough to matter.
// __launch_bounds__ is not decoration here. This is the only kernel in the
// file that asks for a full 1024-thread block, and a T4 has 65536 registers
// per SM, so ptxas has to keep it under 64 registers a thread or the launch
// fails with cudaErrorLaunchOutOfResources. Telling ptxas the block size lets
// it enforce that at compile time instead. Day 17 covers the trade.
__global__ void __launch_bounds__(kScanThreads) scanOffsets(
const unsigned int* __restrict__ hist, unsigned int* __restrict__ offsets) {
__shared__ unsigned int s[kHistSlots];
const unsigned int tid = threadIdx.x;
const unsigned int n = static_cast<unsigned int>(kHistSlots);
s[2u * tid] = hist[2u * tid];
s[2u * tid + 1u] = hist[2u * tid + 1u];
unsigned int offset = 1u;
for (unsigned int d = n >> 1; d > 0u; d >>= 1) {
__syncthreads();
if (tid < d) {
const unsigned int ai = offset * (2u * tid + 1u) - 1u;
const unsigned int bi = offset * (2u * tid + 2u) - 1u;
s[bi] += s[ai];
}
offset *= 2u;
}
// The root now holds the total. Clearing it is what makes the downsweep
// produce an exclusive scan rather than an inclusive one.
__syncthreads();
if (tid == 0u) {
s[n - 1u] = 0u;
}
for (unsigned int d = 1u; d < n; d *= 2u) {
offset >>= 1;
__syncthreads();
if (tid < d) {
const unsigned int ai = offset * (2u * tid + 1u) - 1u;
const unsigned int bi = offset * (2u * tid + 2u) - 1u;
const unsigned int t = s[ai];
s[ai] = s[bi];
s[bi] += t;
}
}
__syncthreads();
offsets[2u * tid] = s[2u * tid];
offsets[2u * tid + 1u] = s[2u * tid + 1u];
}
// Writes one pass's keys to their final positions.
//
// One thread owns one key of a tile. The block walks its span one tile at a
// time, in input order, and sorts each tile by its digit inside shared memory
// with `bits` stable one-bit splits before writing anything out.
//
// Memory: the read is coalesced, 32 consecutive keys per warp. The write is
// not, and cannot be, because moving keys somewhere else is the job. Sorting
// the tile first is what turns it into `buckets` contiguous runs rather than
// one scattered address per lane.
//
// Launch assumption: exactly kThreadsPerBlock threads per block, the same
// gridDim.x the histogram was counted with, and `span` a whole number of
// tiles. Offsets are unsigned int, which caps this program at 2^32 keys.
__global__ void scatterPass(const unsigned int* __restrict__ in,
unsigned int* __restrict__ out,
const unsigned int* __restrict__ offsets,
size_t span, int shift, int bits,
unsigned int mask) {
__shared__ unsigned int sKeys[kThreadsPerBlock];
__shared__ unsigned int sScan[kThreadsPerBlock];
__shared__ unsigned int sBase[kBuckets];
__shared__ unsigned int sStart[kBuckets];
__shared__ unsigned int sCount[kBuckets];
const unsigned int tid = threadIdx.x;
const unsigned int buckets = mask + 1u;
// This block's global write cursor per digit, straight from the scan. It
// advances by one tile's counts at the bottom of every iteration, which
// is what keeps the whole span stable rather than only each tile.
for (unsigned int b = tid; b < buckets; b += blockDim.x) {
sBase[b] = offsets[b * gridDim.x + blockIdx.x];
}
__syncthreads();
const size_t blockStart = blockIdx.x * span;
for (size_t tile = 0; tile < span; tile += blockDim.x) {
sKeys[tid] = in[blockStart + tile + tid];
__syncthreads();
// snippet: split
// `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();
}
// end snippet
for (unsigned int b = tid; b < buckets; b += blockDim.x) {
sCount[b] = 0u;
sStart[b] = 0u;
}
__syncthreads();
// snippet: scatter
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();
// end snippet
for (unsigned int b = tid; b < buckets; b += blockDim.x) {
sBase[b] += sCount[b];
}
__syncthreads();
}
}
// The floor. Reads every key and writes it straight back, from the same grid
// over the same buffers in the same process, so a sort row has something on
// this card to be a multiple of. No sort can beat it.
__global__ void copyKeys(const unsigned int* __restrict__ in,
unsigned int* __restrict__ out, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
out[i] = in[i];
}
}
// Fills the keys from a counter-based mix, so element i is a pure function of
// i and the seed and the same run produces the same bytes on any machine.
static void fillKeys(unsigned int* keys, size_t n, unsigned int seed) {
for (size_t i = 0; i < n; ++i) {
unsigned int x = static_cast<unsigned int>(i) + seed;
x ^= x >> 16;
x *= 2246822519u;
x ^= x >> 13;
x *= 3266489917u;
x ^= x >> 16;
keys[i] = x;
}
}
// The reference: std::sort over a copy of the input. Comparing the GPU output
// against it element for element is both checks at once, because two arrays
// that match element for element and one of which is sorted hold the same
// multiset in the same order. Written for obvious correctness, not speed.
static void sortCpu(const unsigned int* in, unsigned int* out, size_t n) {
for (size_t i = 0; i < n; ++i) {
out[i] = in[i];
}
std::sort(out, out + n);
}
// Exact, not within a tolerance. These are integer keys, there is no rounding
// to allow for, and a tolerance would hide a lost key.
static size_t firstMismatch(const unsigned int* got, const unsigned int* want,
size_t n) {
for (size_t i = 0; i < n; ++i) {
if (got[i] != want[i]) {
return i;
}
}
return n;
}
// Times a launch with CUDA events and returns the mean milliseconds per run.
//
// This is the one template and the one lambda the early modules allow. Copy
// it verbatim; the alternative is six copies of the event boilerplate, which
// is how a warm-up goes missing from one of them.
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
// Warm up this kernel, not just the first kernel in the program. Lazy
// module loading has been the default since CUDA 12.2 on Linux, so the
// first launch of each kernel pays its own load.
for (int i = 0; i < kWarmupRuns; ++i) {
launch();
}
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < kTimedRuns; ++i) {
launch();
}
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError());
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / kTimedRuns;
}
// Runs one whole sort and returns the buffer holding the result.
//
// Pass 0 reads `src` and every later pass ping-pongs between `a` and `b`, so
// `src` is never written. That is what lets the whole sort go inside
// timeKernel's loop: each run starts from the same unsorted keys and there is
// no restore inside the timed region.
//
// There is no error check between the launches on purpose. The caller checks
// once after the whole sequence, because cudaGetLastError() holds the first
// error until somebody reads it, and a synchronize in here would land inside
// the measurement.
static unsigned int* radixSort(const unsigned int* src, unsigned int* a,
unsigned int* b, unsigned int* hist,
unsigned int* offsets, int bits) {
const unsigned int mask = (1u << bits) - 1u;
const int passes = kKeyBits / bits;
unsigned int* result = a;
for (int p = 0; p < passes; ++p) {
const int shift = p * bits;
const unsigned int* in = (p == 0) ? src : ((p % 2 == 1) ? a : b);
unsigned int* out = ((p % 2) == 0) ? a : b;
countDigits<<<kBlocks, kThreadsPerBlock>>>(in, hist, kSpan, shift,
mask);
scanOffsets<<<1, kScanThreads>>>(hist, offsets);
scatterPass<<<kBlocks, kThreadsPerBlock>>>(in, out, offsets, kSpan,
shift, bits, mask);
result = out;
}
return result;
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name, prop.major,
prop.minor);
std::printf("%d SMs, %zu bytes of shared memory per block\n",
prop.multiProcessorCount, prop.sharedMemPerBlock);
if (prop.maxThreadsPerBlock < kScanThreads) {
std::fprintf(stderr,
"scanOffsets wants %d threads in one block and this "
"device allows %d\n",
kScanThreads, prop.maxThreadsPerBlock);
return EXIT_FAILURE;
}
const size_t keyBytes = kPadded * sizeof(unsigned int);
const size_t realBytes = kElems * sizeof(unsigned int);
const size_t histBytes =
static_cast<size_t>(kHistSlots) * sizeof(unsigned int);
std::printf("\n%zu keys, padded to %zu, %.1f MiB per buffer\n", kElems,
kPadded, static_cast<double>(keyBytes) / (1024.0 * 1024.0));
std::printf("%d blocks x %d threads, each owns %zu keys (%zu tiles)\n",
kBlocks, kThreadsPerBlock, kSpan,
kSpan / static_cast<size_t>(kThreadsPerBlock));
std::vector<unsigned int> h_keys(kElems);
std::vector<unsigned int> h_want(kElems);
std::vector<unsigned int> h_got(kElems);
fillKeys(h_keys.data(), kElems, kSeed);
sortCpu(h_keys.data(), h_want.data(), kElems);
unsigned int* d_src = nullptr;
unsigned int* d_a = nullptr;
unsigned int* d_b = nullptr;
unsigned int* d_hist = nullptr;
unsigned int* d_offsets = nullptr;
CUDA_CHECK(cudaMalloc(&d_src, keyBytes));
CUDA_CHECK(cudaMalloc(&d_a, keyBytes));
CUDA_CHECK(cudaMalloc(&d_b, keyBytes));
CUDA_CHECK(cudaMalloc(&d_hist, histBytes));
CUDA_CHECK(cudaMalloc(&d_offsets, histBytes));
// The pad keys are 0xFFFFFFFF, the largest key there is, and every pass is
// stable, so they end up after every real key and the first kElems entries
// of the output are the answer even if the input contains 0xFFFFFFFF too.
CUDA_CHECK(
cudaMemcpy(d_src, h_keys.data(), realBytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemset(d_src + kElems, 0xFF, keyBytes - realBytes));
// scanOffsets scans all kHistSlots slots, but the exclusive prefix at
// index i depends only on slots below i, and countDigits writes every
// slot below (mask + 1) * kBlocks on this pass. Slots above that carry
// whatever an earlier pass left and cannot reach an offset anyone reads.
// The memset is here so compute-sanitizer --tool initcheck stays quiet
// on the first pass.
CUDA_CHECK(cudaMemset(d_hist, 0, histBytes));
const int bitsList[kVariants] = {kNarrowBits, kBitsPerPass};
const char* nameList[kVariants] = {"1-bit split", "4-bit digit"};
float sortMs[kVariants] = {0.0f, 0.0f};
// Every failure below records itself and falls through to the one cleanup
// block at the bottom, so no path returns with device memory allocated.
int status = EXIT_SUCCESS;
for (int v = 0; v < kVariants; ++v) {
unsigned int* d_result =
radixSort(d_src, d_a, d_b, d_hist, d_offsets, bitsList[v]);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_got.data(), d_result, realBytes,
cudaMemcpyDeviceToHost));
const size_t bad = firstMismatch(h_got.data(), h_want.data(), kElems);
if (bad != kElems) {
std::fprintf(stderr, "%s wrong at %zu: got %u, want %u\n",
nameList[v], bad, h_got[bad], h_want[bad]);
status = EXIT_FAILURE;
}
}
if (status == EXIT_SUCCESS) {
const float copyMs = timeKernel([&] {
copyKeys<<<kBlocks, kThreadsPerBlock>>>(d_src, d_a, kPadded);
});
for (int v = 0; v < kVariants; ++v) {
const int bits = bitsList[v];
sortMs[v] = timeKernel(
[&] { radixSort(d_src, d_a, d_b, d_hist, d_offsets, bits); });
}
if (copyMs <= 0.0f) {
std::fprintf(stderr, "copy timed at %.6f ms, which cannot be\n",
static_cast<double>(copyMs));
status = EXIT_FAILURE;
}
for (int v = 0; v < kVariants; ++v) {
if (sortMs[v] <= 0.0f) {
std::fprintf(stderr, "%s timed at %.6f ms, which cannot be\n",
nameList[v], static_cast<double>(sortMs[v]));
status = EXIT_FAILURE;
}
}
if (status == EXIT_SUCCESS) {
const double copyGb = 2.0 * static_cast<double>(keyBytes) /
(static_cast<double>(copyMs) * 1.0e-3) /
1.0e9;
std::printf(
"\nBoth rows sorted the same keys and matched std::sort at "
"all %zu elements.\n\n",
kElems);
std::printf(
"kernel passes time (ms) ms/pass Mkeys/s "
"x copy\n");
std::printf(
"------------- ------ --------- ------- ------- "
"------\n");
std::printf("%-13s %s %9.3f %7s %7s %6.2f\n",
"copy, no sort", "-", static_cast<double>(copyMs), "-",
"-", 1.0);
for (int v = 0; v < kVariants; ++v) {
const int passes = kKeyBits / bitsList[v];
const double ms = static_cast<double>(sortMs[v]);
const double mkeys =
static_cast<double>(kElems) / (ms * 1.0e-3) / 1.0e6;
std::printf("%-13s %6d %9.3f %7.3f %7.1f %6.2f\n",
nameList[v], passes, ms, ms / passes, mkeys,
ms / static_cast<double>(copyMs));
}
std::printf(
"\nThe copy row moves %zu bytes in and %zu out, %.1f GB/s.\n",
keyBytes, keyBytes, copyGb);
std::printf(
"Each sort pass reads the keys to count them, reads them "
"again to\nscatter them and writes them once: %.1f MiB per "
"pass either way.\n",
3.0 * static_cast<double>(keyBytes) / (1024.0 * 1024.0));
}
}
CUDA_CHECK(cudaFree(d_src));
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_hist));
CUDA_CHECK(cudaFree(d_offsets));
return status;
}