code/day39-cccl/cccl.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 39: when to stop hand-writing. Thrust, CUB and libcu++.
//
// Three algorithms this course already built by hand, each run twice in one
// process on one input: once with the kernels from the earlier day, once with
// the CCCL call that replaces them.
//
// Part 1 reduce hand (day 24 version 4) vs thrust::reduce vs
// cub::DeviceReduce::Sum
// Part 2 scan hand (two-level Hillis-Steele) vs
// cub::DeviceScan::ExclusiveSum
// Part 3 sort hand (32-pass bit split, built on the scan above) vs
// cub::DeviceRadixSort::SortKeys
// Part 4 one kernel using cub::BlockReduce, cuda::atomic_ref and
// cuda::std::numeric_limits together
//
// The hand-written reduction here is the ladder's fourth rung, not its
// eighth. Day 25 goes four rungs further and the lesson says so; porting
// day 25's version 8 into this file is day 39's exercise. Every other
// hand-written routine is the plainest correct version of its algorithm,
// which is what a reader who has just finished those days actually has.
//
// Fairness rules this file follows, because a library comparison is the
// easiest benchmark in the world to rig:
//
// - Every row reads the same input buffer and produces a checked result.
// The input is never modified, so the tenth timed run sees the same data
// as the first. A radix sort handed its own sorted output is a different
// benchmark.
// - Temp-storage queries and allocations sit outside every timed region,
// for the library rows and the hand-written rows alike.
// - The reduce rows report GB/s from the input bytes only, because that is
// the one figure all three share. Two of the three are closed libraries
// and this program cannot count their internal passes.
// - The sort rows report no GB/s at all. How many bytes a sort moves is a
// property of the algorithm, so a bandwidth column would be comparing
// two different quantities and calling it a ratio.
//
// What this program does not check: that any library row is faster than any
// hand-written row. That is the claim the run exists to test, so encoding it
// as a gate would make the test unfalsifiable.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o cccl cccl.cu
// Run: ./cccl
//
// 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>
// CCCL, all three libraries. Nothing below goes on the link line: nvcc puts
// these headers on the include path itself and every one of them is header
// only. Include what you use, per CCCL's own compatibility guidelines, rather
// than the umbrella <cub/cub.cuh>.
#include <cub/block/block_reduce.cuh>
#include <cub/device/device_radix_sort.cuh>
#include <cub/device/device_reduce.cuh>
#include <cub/device/device_scan.cuh>
#include <cub/version.cuh>
#include <cuda/atomic>
#include <cuda/std/limits>
#include <thrust/execution_policy.h>
#include <thrust/reduce.h>
#include <thrust/version.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)
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kScanThreads = 256; // one block scans this many elements
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr int kTimedRows = 7;
// 2^24 minus 611, the same count day 24 reduces, so the reduce rows here can
// be read against that day's table on the same card. Two properties carry
// weight. The sum of that many 1.0f values is 16,776,605, which is under
// 2^24, so every partial sum in every version lands on a float exactly and
// the check below needs no tolerance. And 611 is 13 x 47, so no block size
// divides it and the bounds check in every kernel runs on every launch.
constexpr size_t kElems = 16777216ull - 611ull;
// Levels of the hand-written scan. Each level divides the count by
// kScanThreads, so eight levels covers 256^8 elements and the depth is a
// runtime check rather than a promise.
constexpr int kMaxScanLevels = 8;
// Knuth's multiplicative hash constant. Multiplying by an odd number is a
// bijection modulo 2^32 and so is x ^= x >> 16, so composing them gives
// 16,776,605 distinct keys. Distinct keys make the sort check exact: there is
// only one correct output and stability cannot hide a bug.
constexpr unsigned int kHashMultiplier = 2654435761u;
// The key width, read off the type by libcu++ rather than typed as 32.
// cuda::std::numeric_limits is the same std::numeric_limits you know, usable
// in host and device code from one header, and this is the smallest honest
// use of it in the program.
constexpr int kKeyBits = cuda::std::numeric_limits<unsigned int>::digits;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert((kScanThreads & (kScanThreads - 1)) == 0,
"the block scan doubles its offset, so the scan block size must "
"be a power of two");
static_assert(kElems <= 2147483647ull,
"the CUB device entry points used here take an int item count");
static_assert(kKeyBits > 0 && kKeyBits % 2 == 0,
"the sort alternates between two buffers, so an even pass count "
"is what puts the result in the buffer the caller is handed");
// ---------------------------------------------------------------------------
// The hand-written half.
// ---------------------------------------------------------------------------
// snippet: hand-reduce
// Day 24's version 4: sequential addressing with the first add done during
// the load. One block covers 2 * blockDim.x elements.
//
// One thread: adds its two elements, then walks the halving tree in shared
// memory until thread 0 holds the block's sum.
//
// One warp: two global loads, each 32 consecutive floats and 128 contiguous
// bytes, blockDim.x elements apart. They are independent, so both are in
// flight before the barrier. The shared reads are words tid and tid + s,
// consecutive across the warp, so no bank conflicts.
//
// Launch assumption: exactly kThreadsPerBlock threads per block, and a grid
// covering ceil(n / (2 * kThreadsPerBlock)) blocks.
__global__ void reduceTwoLoads(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const unsigned int width = kThreadsPerBlock;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) * 2 + tid;
// The guard covers the loads, never the barrier. Every thread in the
// block reaches every __syncthreads() below.
float sum = (i < n) ? in[i] : 0.0f;
if (i + blockDim.x < n) {
sum += in[i + blockDim.x];
}
tile[tid] = sum;
__syncthreads();
for (unsigned int s = width / 2; s > 0; s >>= 1) {
if (tid < s) {
tile[tid] += tile[tid + s];
}
__syncthreads();
}
if (tid == 0) {
out[blockIdx.x] = tile[0];
}
}
// Runs one complete reduction, n floats down to one, and returns the buffer
// holding the answer. The two scratch buffers alternate, so no launch ever
// reads and writes the same allocation.
//
// Launches only. There is no synchronize and no error check inside, because
// this runs in the timed region and a synchronize there would measure
// something other than the kernels. main() checks the launches on the untimed
// call it makes first.
static const float* reduceAll(const float* d_in, size_t n, float* d_a,
float* d_b, int* passes) {
const size_t per = 2 * static_cast<size_t>(kThreadsPerBlock);
const float* src = d_in;
float* dst = d_a;
size_t m = n;
int count = 0;
while (m > 1) {
const int blocks = static_cast<int>((m + per - 1) / per);
reduceTwoLoads<<<blocks, kThreadsPerBlock>>>(src, dst, m);
++count;
m = static_cast<size_t>(blocks);
src = dst;
dst = (dst == d_a) ? d_b : d_a;
}
*passes = count;
return src;
}
// end snippet
// snippet: hand-scan
// One block's exclusive prefix sum over kScanThreads elements, Hillis-Steele
// in shared memory, plus the block's total written to sums[blockIdx.x].
//
// One thread: owns one element, and at step `off` adds the partial sum of the
// thread `off` places to its left.
//
// One warp: the global load is 32 consecutive words, 128 contiguous bytes.
// The shared reads are words tid and tid - off, consecutive across the warp,
// so no bank conflicts. Two tiles, not one: the step writes the buffer it is
// not reading, which is what removes the second barrier per step.
//
// In-place is safe. A block reads its whole slice into shared memory before
// it writes any of it back, and no block touches another block's range, so
// `in` and `out` may be the same array. The upper levels use that.
//
// Launch assumption: exactly kScanThreads threads per block.
__global__ void scanBlocksExclusive(const unsigned int* __restrict__ in,
unsigned int* __restrict__ out,
unsigned int* __restrict__ sums, size_t n) {
__shared__ unsigned int tile[2][kScanThreads];
const unsigned int tid = threadIdx.x;
const unsigned int width = kScanThreads;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
const unsigned int mine = (i < n) ? in[i] : 0u;
unsigned int pin = 0u;
tile[pin][tid] = mine;
__syncthreads();
for (unsigned int off = 1u; off < width; off *= 2u) {
unsigned int v = tile[pin][tid];
if (tid >= off) {
v += tile[pin][tid - off];
}
tile[pin ^ 1u][tid] = v;
__syncthreads();
pin ^= 1u;
}
// tile[pin][tid] is the inclusive scan, so subtracting this thread's own
// element turns it into the exclusive one.
if (i < n) {
out[i] = tile[pin][tid] - mine;
}
if (tid == width - 1) {
sums[blockIdx.x] = tile[pin][tid];
}
}
// Adds each block's offset back onto that block's slice. offsets[b] is the
// exclusive scan of the block totals, so block b needs exactly one word of it.
__global__ void addBlockOffsets(unsigned int* __restrict__ data,
const unsigned int* __restrict__ offsets,
size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
data[i] += offsets[blockIdx.x];
}
}
// The block-sum buffer for each level of the scan. Level 0 is the caller's
// array and owns no buffer; level L holds one word per block of level L - 1.
struct ScanLevels {
size_t count[kMaxScanLevels];
unsigned int* buffer[kMaxScanLevels];
int depth;
};
// Fills in the level counts for n elements. Returns false when n needs more
// levels than the fixed array holds, which the caller turns into a failure
// rather than a silent overrun.
static bool buildScanLevels(ScanLevels* levels, size_t n) {
levels->depth = 0;
levels->count[0] = n;
levels->buffer[0] = nullptr;
while (levels->count[levels->depth] > 1) {
if (levels->depth + 1 >= kMaxScanLevels) {
return false;
}
const size_t m = levels->count[levels->depth];
levels->count[levels->depth + 1] =
(m + static_cast<size_t>(kScanThreads) - 1) /
static_cast<size_t>(kScanThreads);
++levels->depth;
levels->buffer[levels->depth] = nullptr;
}
return true;
}
// One complete exclusive scan of levels.count[0] unsigned ints.
//
// Down the levels: scan each level's blocks and collect their totals into the
// level above. Back up: add the scanned totals of the level above onto each
// block. At 16,776,605 elements and 256 threads a block that is three levels
// down and two back up, five launches, and every element of the input is read
// twice and written twice.
//
// Launches only, for the same reason reduceAll is.
static void scanAll(const unsigned int* d_in, unsigned int* d_out,
const ScanLevels& levels) {
for (int level = 0; level < levels.depth; ++level) {
const unsigned int* src = (level == 0) ? d_in : levels.buffer[level];
unsigned int* dst = (level == 0) ? d_out : levels.buffer[level];
const size_t m = levels.count[level];
const int blocks = static_cast<int>(
(m + static_cast<size_t>(kScanThreads) - 1) / kScanThreads);
scanBlocksExclusive<<<blocks, kScanThreads>>>(
src, dst, levels.buffer[level + 1], m);
}
for (int level = levels.depth - 2; level >= 0; --level) {
unsigned int* dst = (level == 0) ? d_out : levels.buffer[level];
const size_t m = levels.count[level];
const int blocks = static_cast<int>(
(m + static_cast<size_t>(kScanThreads) - 1) / kScanThreads);
addBlockOffsets<<<blocks, kScanThreads>>>(dst, levels.buffer[level + 1],
m);
}
}
// end snippet
// snippet: hand-sort
// flags[i] is 1 when bit `bit` of keys[i] is zero. Scanning that array gives
// every element the count of zero-bit elements ahead of it, which is where it
// goes.
//
// One thread: reads one key, writes one flag.
// One warp: 32 consecutive keys in, 32 consecutive flags out, both 128
// contiguous bytes.
__global__ void markZeroBits(const unsigned int* __restrict__ keys,
unsigned int* __restrict__ flags, size_t n,
int bit) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
flags[i] = ((keys[i] >> bit) & 1u) ^ 1u;
}
}
// Moves every key to its place for this bit. Elements with a zero bit keep
// their order at the front, elements with a one bit keep theirs behind them,
// so one pass is a stable partition and thirty-two of them sort.
//
// One thread: one key read, one key written.
// One warp: the reads are three coalesced 128-byte runs. The writes are not:
// two lanes next to each other land at two different ends of the output. That
// scatter is what this kernel costs, and it is the same scatter a real radix
// sort spends its design budget avoiding.
__global__ void scatterByBit(const unsigned int* __restrict__ keys,
const unsigned int* __restrict__ flags,
const unsigned int* __restrict__ scanned,
unsigned int* __restrict__ out, size_t n) {
// Every thread reads the same two words, which the hardware serves as a
// broadcast. Computing the total here keeps the whole sort on the device:
// reading it back to the host would put a synchronize per bit inside the
// timed region. Index n - 1 is always in range, so it needs no guard.
const unsigned int totalZeros = scanned[n - 1] + flags[n - 1];
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
const size_t zerosBefore = static_cast<size_t>(scanned[i]);
const size_t dst = (flags[i] != 0u) ? zerosBefore
: static_cast<size_t>(totalZeros) +
(i - zerosBefore);
out[dst] = keys[i];
}
}
// Least-significant-bit radix sort, one bit at a time, on top of the scan
// above. Returns the buffer holding the sorted keys.
//
// d_keys is never written, so every timed run starts from the same unsorted
// input. The two scratch buffers alternate, so no launch reads and writes the
// same allocation.
//
// Launches only, for the same reason reduceAll is.
static const unsigned int* radixSortHand(const unsigned int* d_keys,
unsigned int* d_a, unsigned int* d_b,
unsigned int* d_flags,
unsigned int* d_scan,
const ScanLevels& levels) {
const size_t n = levels.count[0];
const int blocks = static_cast<int>(
(n + static_cast<size_t>(kThreadsPerBlock) - 1) / kThreadsPerBlock);
const unsigned int* src = d_keys;
unsigned int* dst = d_a;
for (int bit = 0; bit < kKeyBits; ++bit) {
markZeroBits<<<blocks, kThreadsPerBlock>>>(src, d_flags, n, bit);
scanAll(d_flags, d_scan, levels);
scatterByBit<<<blocks, kThreadsPerBlock>>>(src, d_flags, d_scan, dst,
n);
src = dst;
dst = (dst == d_a) ? d_b : d_a;
}
return src;
}
// end snippet
// ---------------------------------------------------------------------------
// The library half. Everything above, again.
// ---------------------------------------------------------------------------
// snippet: library-calls
// 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)));
}
// end snippet
// Counts the odd keys, with three CCCL layers in one kernel: cub::BlockReduce
// for the block-scope sum, cuda::atomic_ref for the device-scope accumulate,
// and cuda::std::numeric_limits above for the key width.
//
// cuda::atomic_ref is the C++20 atomic_ref with a thread scope attached. The
// scope is the promise you are making about who else touches this address:
// thread_scope_block, thread_scope_device or thread_scope_system, narrowest
// first. Day 27 is where the cost of each was measured.
//
// One thread: tests one key's low bit and contributes 0 or 1.
// One warp: 32 consecutive keys, 128 contiguous bytes, four 32-byte sectors.
//
// Launch assumption: exactly kThreadsPerBlock threads per block, because
// cub::BlockReduce is specialised on that number at compile time. Every
// thread must reach the Sum() call; it is collective and it has barriers of
// its own inside.
// snippet: libcudacxx
__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);
}
}
// end snippet
// ---------------------------------------------------------------------------
// Times a launch with CUDA events and returns the mean milliseconds per run.
// A host clock around a launch measures the launch, because launches are
// asynchronous. Day 9 takes that apart. Copy this helper verbatim.
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;
}
// GB/s from a byte count and a mean time in milliseconds.
static double bandwidth(double bytes, float ms) {
return bytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
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), %d SMs\n", prop.name,
prop.major, prop.minor, prop.multiProcessorCount);
// The version map, read from the headers and the runtime rather than
// remembered. Thrust and CUB carry the CCCL version they shipped in.
int runtimeVersion = 0;
CUDA_CHECK(cudaRuntimeGetVersion(&runtimeVersion));
std::printf("CUDA runtime %d.%d ships Thrust %d.%d.%d and CUB %d.%d.%d\n",
runtimeVersion / 1000, (runtimeVersion % 1000) / 10,
THRUST_MAJOR_VERSION, THRUST_MINOR_VERSION,
THRUST_SUBMINOR_VERSION, CUB_MAJOR_VERSION, CUB_MINOR_VERSION,
CUB_SUBMINOR_VERSION);
std::printf("Key width from cuda::std::numeric_limits: %d bits\n",
kKeyBits);
// Every failure below records itself and falls through to the one cleanup
// block at the bottom, so no path can return with device memory still
// allocated. None of them is an assert(): CI builds Release, Release
// defines NDEBUG, and NDEBUG deletes assert(), so the check would be
// missing from exactly the build that matters.
int failures = 0;
int rows = 0;
ScanLevels levels;
if (!buildScanLevels(&levels, kElems)) {
std::fprintf(stderr,
"%zu elements needs more than %d scan levels at %d "
"threads a block\n",
kElems, kMaxScanLevels, kScanThreads);
return EXIT_FAILURE;
}
const size_t floatBytes = kElems * sizeof(float);
const size_t keyBytes = kElems * sizeof(unsigned int);
const size_t reduceScratch =
(kElems + 2 * kThreadsPerBlock - 1) / (2 * kThreadsPerBlock);
const int blocks =
static_cast<int>((kElems + static_cast<size_t>(kThreadsPerBlock) - 1) /
kThreadsPerBlock);
float* d_in = nullptr;
float* d_ra = nullptr;
float* d_rb = nullptr;
float* d_reduceOut = nullptr;
unsigned int* d_keys = nullptr;
unsigned int* d_sa = nullptr;
unsigned int* d_sb = nullptr;
unsigned int* d_sortOut = nullptr;
unsigned int* d_flags = nullptr;
unsigned int* d_scan = nullptr;
int* d_oddCount = nullptr;
void* d_temp = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, floatBytes));
CUDA_CHECK(cudaMalloc(&d_ra, reduceScratch * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_rb, reduceScratch * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_reduceOut, sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_keys, keyBytes));
CUDA_CHECK(cudaMalloc(&d_sa, keyBytes));
CUDA_CHECK(cudaMalloc(&d_sb, keyBytes));
CUDA_CHECK(cudaMalloc(&d_sortOut, keyBytes));
CUDA_CHECK(cudaMalloc(&d_flags, keyBytes));
CUDA_CHECK(cudaMalloc(&d_scan, keyBytes));
CUDA_CHECK(cudaMalloc(&d_oddCount, sizeof(int)));
for (int level = 1; level <= levels.depth; ++level) {
CUDA_CHECK(cudaMalloc(&levels.buffer[level],
levels.count[level] * sizeof(unsigned int)));
}
std::vector<float> h_in(kElems, 1.0f);
std::vector<unsigned int> h_keys(kElems);
std::vector<unsigned int> h_want(kElems);
std::vector<unsigned int> h_got(kElems);
for (size_t i = 0; i < kElems; ++i) {
unsigned int x = static_cast<unsigned int>(i) * kHashMultiplier;
x ^= x >> 16;
h_keys[i] = x;
}
CUDA_CHECK(
cudaMemcpy(d_in, h_in.data(), floatBytes, cudaMemcpyHostToDevice));
CUDA_CHECK(
cudaMemcpy(d_keys, h_keys.data(), keyBytes, cudaMemcpyHostToDevice));
// Every scan row runs on an array of ones, so the correct exclusive scan
// of element i is exactly i and the check needs no reference array.
// h_got is borrowed to build it rather than a sixth 64 MiB host vector;
// nothing has been read out of it yet.
std::fill(h_got.begin(), h_got.end(), 1u);
CUDA_CHECK(
cudaMemcpy(d_flags, h_got.data(), keyBytes, cudaMemcpyHostToDevice));
// Temp storage, queried once and allocated once, outside every timed
// region. Each entry point keeps its own byte count because CUB validates
// the size it was given against the layout it wants.
size_t reduceTempBytes = 0;
size_t scanTempBytes = 0;
size_t sortTempBytes = 0;
reduceCub(nullptr, reduceTempBytes, d_in, d_reduceOut, kElems);
scanCub(nullptr, scanTempBytes, d_flags, d_scan, kElems);
sortCub(nullptr, sortTempBytes, d_keys, d_sortOut, kElems);
const size_t tempBytes =
std::max(reduceTempBytes, std::max(scanTempBytes, sortTempBytes));
CUDA_CHECK(cudaMalloc(&d_temp, tempBytes));
std::printf("CUB temp storage: reduce %zu B, scan %zu B, sort %zu B\n",
reduceTempBytes, scanTempBytes, sortTempBytes);
// -----------------------------------------------------------------------
// Part 1: reduce.
// -----------------------------------------------------------------------
const float want = static_cast<float>(kElems);
int passes = 0;
const float* d_handSum = reduceAll(d_in, kElems, d_ra, d_rb, &passes);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
float handSum = 0.0f;
CUDA_CHECK(
cudaMemcpy(&handSum, d_handSum, sizeof(float), cudaMemcpyDeviceToHost));
if (handSum != want) {
std::fprintf(stderr, "hand reduce summed to %.9g, want %.9g\n", handSum,
want);
++failures;
}
float thrustSum = reduceThrust(d_in, kElems);
if (thrustSum != want) {
std::fprintf(stderr, "thrust::reduce summed to %.9g, want %.9g\n",
thrustSum, want);
++failures;
}
reduceCub(d_temp, reduceTempBytes, d_in, d_reduceOut, kElems);
CUDA_CHECK(cudaDeviceSynchronize());
float cubSum = 0.0f;
CUDA_CHECK(cudaMemcpy(&cubSum, d_reduceOut, sizeof(float),
cudaMemcpyDeviceToHost));
if (cubSum != want) {
std::fprintf(stderr,
"cub::DeviceReduce::Sum summed to %.9g, want "
"%.9g\n",
cubSum, want);
++failures;
}
const float handReduceMs = timeKernel([&] {
int timedPasses = 0;
reduceAll(d_in, kElems, d_ra, d_rb, &timedPasses);
});
const float thrustReduceMs =
timeKernel([&] { thrustSum = reduceThrust(d_in, kElems); });
const float cubReduceMs = timeKernel(
[&] { reduceCub(d_temp, reduceTempBytes, d_in, d_reduceOut, kElems); });
std::printf("\nPart 1: reduce %zu floats to one, %d hand passes\n", kElems,
passes);
std::printf("%-32s %12s %12s %10s\n", "implementation", "time (ms)", "GB/s",
"vs hand");
std::printf("%-32s %12s %12s %10s\n", "--------------------------------",
"----------", "----------", "--------");
const double inBytes = static_cast<double>(floatBytes);
std::printf("%-32s %12.3f %12.1f %10.2f\n", "hand, day 24 version 4",
handReduceMs, bandwidth(inBytes, handReduceMs), 1.0);
std::printf("%-32s %12.3f %12.1f %10.2f\n", "thrust::reduce",
thrustReduceMs, bandwidth(inBytes, thrustReduceMs),
handReduceMs / thrustReduceMs);
std::printf("%-32s %12.3f %12.1f %10.2f\n", "cub::DeviceReduce::Sum",
cubReduceMs, bandwidth(inBytes, cubReduceMs),
handReduceMs / cubReduceMs);
rows += 3;
// -----------------------------------------------------------------------
// Part 2: exclusive scan.
// -----------------------------------------------------------------------
scanAll(d_flags, d_scan, levels);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_got.data(), d_scan, keyBytes, cudaMemcpyDeviceToHost));
for (size_t i = 0; i < kElems; ++i) {
if (h_got[i] != static_cast<unsigned int>(i)) {
std::fprintf(stderr, "hand scan wrong at %zu: got %u, want %zu\n",
i, h_got[i], i);
++failures;
break;
}
}
scanCub(d_temp, scanTempBytes, d_flags, d_scan, kElems);
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_got.data(), d_scan, keyBytes, cudaMemcpyDeviceToHost));
for (size_t i = 0; i < kElems; ++i) {
if (h_got[i] != static_cast<unsigned int>(i)) {
std::fprintf(stderr,
"cub::DeviceScan::ExclusiveSum wrong at %zu: got %u, "
"want %zu\n",
i, h_got[i], i);
++failures;
break;
}
}
const float handScanMs =
timeKernel([&] { scanAll(d_flags, d_scan, levels); });
const float cubScanMs = timeKernel(
[&] { scanCub(d_temp, scanTempBytes, d_flags, d_scan, kElems); });
std::printf("\nPart 2: exclusive scan of %zu uints, %d levels\n", kElems,
levels.depth);
std::printf("%-32s %12s %12s %10s\n", "implementation", "time (ms)", "GB/s",
"vs hand");
std::printf("%-32s %12s %12s %10s\n", "--------------------------------",
"----------", "----------", "--------");
const double scanBytes = 2.0 * static_cast<double>(keyBytes);
std::printf("%-32s %12.3f %12.1f %10.2f\n", "hand, two-level Hillis-Steele",
handScanMs, bandwidth(scanBytes, handScanMs), 1.0);
std::printf("%-32s %12.3f %12.1f %10.2f\n", "cub::DeviceScan::ExclusiveSum",
cubScanMs, bandwidth(scanBytes, cubScanMs),
handScanMs / cubScanMs);
rows += 2;
// -----------------------------------------------------------------------
// Part 3: sort.
// -----------------------------------------------------------------------
h_want = h_keys;
std::sort(h_want.begin(), h_want.end());
const unsigned int* d_handSorted =
radixSortHand(d_keys, d_sa, d_sb, d_flags, d_scan, levels);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_got.data(), d_handSorted, keyBytes,
cudaMemcpyDeviceToHost));
for (size_t i = 0; i < kElems; ++i) {
if (h_got[i] != h_want[i]) {
std::fprintf(stderr, "hand sort wrong at %zu: got %u, want %u\n", i,
h_got[i], h_want[i]);
++failures;
break;
}
}
sortCub(d_temp, sortTempBytes, d_keys, d_sortOut, kElems);
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_got.data(), d_sortOut, keyBytes, cudaMemcpyDeviceToHost));
for (size_t i = 0; i < kElems; ++i) {
if (h_got[i] != h_want[i]) {
std::fprintf(stderr,
"cub::DeviceRadixSort::SortKeys wrong at %zu: got %u, "
"want %u\n",
i, h_got[i], h_want[i]);
++failures;
break;
}
}
const float handSortMs = timeKernel(
[&] { radixSortHand(d_keys, d_sa, d_sb, d_flags, d_scan, levels); });
const float cubSortMs = timeKernel(
[&] { sortCub(d_temp, sortTempBytes, d_keys, d_sortOut, kElems); });
std::printf("\nPart 3: sort %zu distinct uint keys, %d hand passes\n",
kElems, kKeyBits);
std::printf("%-32s %12s %10s\n", "implementation", "time (ms)", "vs hand");
std::printf("%-32s %12s %10s\n", "--------------------------------",
"----------", "--------");
std::printf("%-32s %12.3f %10.2f\n", "hand, one bit per pass", handSortMs,
1.0);
std::printf("%-32s %12.3f %10.2f\n", "cub::DeviceRadixSort::SortKeys",
cubSortMs, handSortMs / cubSortMs);
rows += 2;
std::printf(
"\nNo GB/s column here. How many bytes a sort moves is a "
"property of\nthe algorithm, so the two rows would not be "
"counting the same thing.\n");
// -----------------------------------------------------------------------
// Part 4: libcu++ and CUB inside one kernel of your own.
// -----------------------------------------------------------------------
size_t wantOdd = 0;
for (size_t i = 0; i < kElems; ++i) {
wantOdd += (h_keys[i] & 1u);
}
int zero = 0;
CUDA_CHECK(
cudaMemcpy(d_oddCount, &zero, sizeof(int), cudaMemcpyHostToDevice));
countOddKeys<<<blocks, kThreadsPerBlock>>>(d_keys, d_oddCount, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
int gotOdd = 0;
CUDA_CHECK(
cudaMemcpy(&gotOdd, d_oddCount, sizeof(int), cudaMemcpyDeviceToHost));
std::printf("\nPart 4: cub::BlockReduce plus cuda::atomic_ref\n");
std::printf("odd keys counted on the device: %d, host says %zu\n", gotOdd,
wantOdd);
if (static_cast<size_t>(gotOdd) != wantOdd) {
std::fprintf(stderr, "the device count and the host count disagree\n");
++failures;
}
// The row count is checked so the lesson's tables cannot drift from what
// the program prints.
if (rows != kTimedRows) {
std::fprintf(stderr, "printed %d timed rows, expected %d\n", rows,
kTimedRows);
++failures;
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_ra));
CUDA_CHECK(cudaFree(d_rb));
CUDA_CHECK(cudaFree(d_reduceOut));
CUDA_CHECK(cudaFree(d_keys));
CUDA_CHECK(cudaFree(d_sa));
CUDA_CHECK(cudaFree(d_sb));
CUDA_CHECK(cudaFree(d_sortOut));
CUDA_CHECK(cudaFree(d_flags));
CUDA_CHECK(cudaFree(d_scan));
CUDA_CHECK(cudaFree(d_oddCount));
CUDA_CHECK(cudaFree(d_temp));
for (int level = 1; level <= levels.depth; ++level) {
CUDA_CHECK(cudaFree(levels.buffer[level]));
}
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}