code/day31-scan-1/scan.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 31: prefix sum (scan), part 1, Hillis-Steele at block level.
//
// Four kernels over the same input at four block sizes. Each block scans its
// own tile of blockDim.x elements and nothing crosses a block boundary.
//
// scanInclusiveInPlace one shared array, one barrier a step. Racy.
// scanInclusiveNoCopy two arrays, idle threads skip the copy. Wrong.
// scanInclusiveDoubleBuffer two arrays, every thread writes every step.
// scanExclusiveDoubleBuffer the same scan, shifted one to the right.
//
// The first two are wrong on purpose. The program reports how wrong they are
// rather than asserting that they fail, because the shape of each failure is
// the lesson: the race hides inside a single warp and the missing copy does
// not. Only the last two are gated.
//
// Every element is (i % 8) + 1, so every prefix inside a tile is a whole
// number below 2^24 and a float holds it exactly. The comparison is therefore
// != with no tolerance at all: a mismatch here means a wrong set of elements
// was added, never a rounding difference. Day 68 is the lesson where the
// order does change the answer.
//
// Nothing here is timed. A block scan of a few hundred elements is smaller
// than the launch that carries it, so a clock on this program would measure
// day 9's subject rather than this one. Day 32 is where scan gets a
// bandwidth number.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o scan scan.cu
// Run: ./scan
//
// 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 <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)
// 32 on every GPU this course targets. main() reads the device's own warpSize
// back and fails if it disagrees, because the race row of the table below is
// read against this number.
constexpr int kWarpSize = 32;
// The shared tiles are sized from this, not from blockDim.x, so one build
// serves every row of the sweep. A kernel that only ever ran at one block size
// would size its tile from that constant and use half the shared memory; the
// ceiling here is that a block of 32 threads still reserves room for 1024.
constexpr int kMaxThreadsPerBlock = 1024;
// The block size the lesson's prose quotes. Every other size in the sweep is
// exercised by the same build.
constexpr int kThreadsPerBlock = 256; // 8 warps
// 2^20 plus 611. 611 is 13 x 47, so no block size in the sweep divides it and
// the bounds check in every kernel runs on every launch instead of never.
constexpr size_t kElems = 1048576ull + 611ull;
constexpr int kSweep = 4;
constexpr int kKernels = 4;
constexpr int kBlockSizes[kSweep] = {32, 64, 256, 1024};
// log2 of a power of two.
constexpr int intLog2(int n) {
int bits = 0;
while (n > 1) {
n /= 2;
++bits;
}
return bits;
}
// Additions one Hillis-Steele scan of n elements performs.
//
// Step `offset` adds for every element that has an element `offset` to its
// left, which is n - offset of them, and the offsets are 1, 2, 4 up to n / 2.
// That sums to n * log2(n) - (n - 1), which the static_assert below pins.
constexpr int hillisSteeleAdds(int n) {
int total = 0;
for (int offset = 1; offset < n; offset *= 2) {
total += n - offset;
}
return total;
}
// Additions a sequential scan of n elements performs: one per element after
// the first, whatever n is.
constexpr int sequentialAdds(int n) {
return n - 1;
}
static_assert(kThreadsPerBlock % kWarpSize == 0,
"the block size must be a whole number of warps");
static_assert((kThreadsPerBlock & (kThreadsPerBlock - 1)) == 0,
"the scan doubles its offset, so the block size must be a "
"power of two");
static_assert(kThreadsPerBlock <= kMaxThreadsPerBlock,
"the shared tiles are sized from kMaxThreadsPerBlock");
static_assert(2 * kMaxThreadsPerBlock * sizeof(float) <= 49152,
"a double-buffered tile for the largest block must fit the "
"48 KiB a block gets by default on a Tesla T4");
// The work claim this page rests on, settled at compile time rather than
// asserted in prose. Both hold at every power-of-two block size from 32 up,
// which is what the lesson's exercise sweeps, so changing kThreadsPerBlock
// still builds.
static_assert(hillisSteeleAdds(kThreadsPerBlock) ==
kThreadsPerBlock * intLog2(kThreadsPerBlock) -
(kThreadsPerBlock - 1),
"Hillis-Steele over n elements performs n log2(n) - (n - 1) "
"additions");
static_assert(hillisSteeleAdds(kThreadsPerBlock) >
sequentialAdds(kThreadsPerBlock),
"this scan does strictly more addition than a sequential one "
"at every block size above a single element, which is what "
"work-inefficient means");
// The one literal count the lesson quotes, guarded so a different block size
// does not fail the build. An assert that forbids the size the exercise asks
// for is worse than no assert at all.
static_assert(kThreadsPerBlock != 256 || hillisSteeleAdds(256) == 1793,
"at 256 threads a block the scan performs 1793 additions "
"against a sequential scan's 255");
// Version 1. In place: one shared array, one barrier per step.
//
// One thread: loads its element into the tile, then at every step adds the
// element `offset` places to its left, writing back into the array everybody
// else is still reading.
//
// One warp: the global load is 32 consecutive floats, 128 contiguous bytes,
// four 32-byte sectors. The shared reads are words `tid` and `tid - offset`,
// both consecutive across the warp, so they cover 32 distinct banks and
// nothing serialises. Shared memory is not what is wrong with this kernel.
//
// This one is wrong and it is here to be wrong. Thread `tid` reads
// tile[tid - offset] in the same step in which thread `tid - offset` writes
// it, and the only barrier sits below both. Day 14 measured the same shape of
// race and found it invisible at 32 threads a block.
//
// Launch assumption: blockDim.x is a power of two, at most
// kMaxThreadsPerBlock, and the grid covers ceil(n / blockDim.x) blocks.
// snippet: in-place
__global__ void scanInclusiveInPlace(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kMaxThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
// The guard covers the load and the store, never the barrier. Zero is the
// identity for addition, so a thread past the end of the array carries a
// value that changes nobody's prefix.
tile[tid] = (i < n) ? in[i] : 0.0f;
__syncthreads();
for (unsigned int offset = 1; offset < blockDim.x; offset *= 2) {
if (tid >= offset) {
tile[tid] += tile[tid - offset];
}
__syncthreads();
}
if (i < n) {
out[i] = tile[tid];
}
}
// end snippet
// Version 2. Double buffered, and still wrong: the threads with nothing to
// add do not copy themselves across.
//
// One thread: reads the half holding the current state and writes the other
// half, except at the steps where it has no element to its left, where it
// writes nothing at all.
//
// One warp: the same coalesced global load and the same conflict-free shared
// access as version 1. Only the destination changed.
//
// The bug is the one double buffering creates and the in-place version cannot
// have. In place, a thread with nothing to add is correct to do nothing: its
// value is already in the array. Across two halves, doing nothing leaves the
// destination holding whatever was there two steps ago, and that stale value
// is what the next step reads.
//
// Both halves are written by the load below, so the miss is a stale read
// rather than an undefined one and the mismatch count is reproducible. The
// version people actually write initialises one half, and then this same bug
// reads shared memory nobody wrote; `compute-sanitizer --tool initcheck`
// reports that and the answer is whatever the SM happened to hold.
//
// Launch assumption: as version 1.
__global__ void scanInclusiveNoCopy(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[2 * kMaxThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const unsigned int width = blockDim.x;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
const float value = (i < n) ? in[i] : 0.0f;
tile[tid] = value;
tile[width + tid] = value;
__syncthreads();
unsigned int src = 0;
for (unsigned int offset = 1; offset < width; offset *= 2) {
const unsigned int dst = 1 - src;
if (tid >= offset) {
tile[dst * width + tid] =
tile[src * width + tid] + tile[src * width + tid - offset];
}
__syncthreads();
src = dst;
}
if (i < n) {
out[i] = tile[src * width + tid];
}
}
// Version 3. Double buffered, correct.
//
// One thread: reads the half holding the current state, writes the other
// half, and writes it on every step whether or not it had anything to add.
//
// One warp: unchanged from version 1. The extra cost of this version is one
// more shared array and one shared write per idle thread per step, not a
// different memory pattern.
//
// Why one barrier a step is enough. Every read of a half is separated from
// the next write to that half by the barrier at the bottom of the loop, so no
// thread can be writing a word another thread has not finished reading. The
// in-place version needs two barriers a step to make the same promise, and
// the version above ships one.
//
// `src` always names the half holding the complete state, so it names the
// answer when the loop ends and there is no final flip to get wrong.
//
// Launch assumption: as version 1.
// snippet: double-buffer
__global__ void scanInclusiveDoubleBuffer(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[2 * kMaxThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const unsigned int width = blockDim.x;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
tile[tid] = (i < n) ? in[i] : 0.0f;
__syncthreads();
unsigned int src = 0;
for (unsigned int offset = 1; offset < width; offset *= 2) {
const unsigned int dst = 1 - src;
if (tid >= offset) {
tile[dst * width + tid] =
tile[src * width + tid] + tile[src * width + tid - offset];
} else {
tile[dst * width + tid] = tile[src * width + tid];
}
__syncthreads();
src = dst;
}
if (i < n) {
out[i] = tile[src * width + tid];
}
}
// end snippet
// Version 4. The exclusive scan, built from version 3 and one shift.
//
// One thread: runs the same scan, then reads its left neighbour's inclusive
// prefix instead of its own. Thread 0 has no neighbour and takes zero, the
// identity for addition.
//
// One warp: one extra shared read one word to the left, still consecutive
// across the warp and still conflict free.
//
// The read below needs no barrier of its own. The last iteration of the loop
// ended in a __syncthreads(), nothing writes the tile after it, and a shared
// read that races with no write is not a race.
//
// The shift comes after the scan rather than before it. Both orders produce
// the same exclusive result, and they differ in what else the block is
// holding when the kernel ends. The exercise on the lesson page is to build
// the other one and find the difference.
//
// Launch assumption: as version 1.
__global__ void scanExclusiveDoubleBuffer(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[2 * kMaxThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const unsigned int width = blockDim.x;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
tile[tid] = (i < n) ? in[i] : 0.0f;
__syncthreads();
unsigned int src = 0;
for (unsigned int offset = 1; offset < width; offset *= 2) {
const unsigned int dst = 1 - src;
if (tid >= offset) {
tile[dst * width + tid] =
tile[src * width + tid] + tile[src * width + tid - offset];
} else {
tile[dst * width + tid] = tile[src * width + tid];
}
__syncthreads();
src = dst;
}
// snippet: exclusive-shift
if (i < n) {
out[i] = (tid > 0) ? tile[src * width + tid - 1] : 0.0f;
}
// end snippet
}
static void launchKernel(int kernel, int blocks, int threads, const float* d_in,
float* d_out, size_t n) {
switch (kernel) {
case 0:
scanInclusiveInPlace<<<blocks, threads>>>(d_in, d_out, n);
break;
case 1:
scanInclusiveNoCopy<<<blocks, threads>>>(d_in, d_out, n);
break;
case 2:
scanInclusiveDoubleBuffer<<<blocks, threads>>>(d_in, d_out, n);
break;
default:
scanExclusiveDoubleBuffer<<<blocks, threads>>>(d_in, d_out, n);
break;
}
}
// CPU reference, inclusive. One scan per tile, restarting at every tile
// boundary, because that is what a block-level scan produces.
//
// Written for obvious correctness, not speed: a plain loop, no blocking, no
// intrinsics. It never allocates; the caller owns every buffer. It
// accumulates in double where the kernel accumulates in float, and on this
// input both are exact, so the cast back loses nothing.
static void scanInclusiveCpu(const float* in, float* out, size_t n,
size_t tile) {
for (size_t base = 0; base < n; base += tile) {
const size_t end = (base + tile < n) ? base + tile : n;
double running = 0.0;
for (size_t i = base; i < end; ++i) {
running += static_cast<double>(in[i]);
out[i] = static_cast<float>(running);
}
}
}
// CPU reference, exclusive. Same tiles, and each element takes the running
// total from before its own value is added, so element 0 of every tile is 0.
static void scanExclusiveCpu(const float* in, float* out, size_t n,
size_t tile) {
for (size_t base = 0; base < n; base += tile) {
const size_t end = (base + tile < n) ? base + tile : n;
double running = 0.0;
for (size_t i = base; i < end; ++i) {
out[i] = static_cast<float>(running);
running += static_cast<double>(in[i]);
}
}
}
// Counts the elements where got and want differ and reports the first one.
// Exact, with no tolerance: see the header note on why every prefix in this
// input is a whole number a float represents exactly. Returning the count as
// well as the index matters here, because how many elements a version gets
// wrong is the part that separates the two broken kernels.
static size_t countMismatches(const float* got, const float* want, size_t n,
size_t* firstBad) {
size_t bad = 0;
*firstBad = n;
for (size_t i = 0; i < n; ++i) {
if (got[i] != want[i]) {
if (bad == 0) {
*firstBad = i;
}
++bad;
}
}
return bad;
}
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), warp size %d\n", prop.name,
prop.major, prop.minor, prop.warpSize);
std::printf("Max threads per block: %d, shared memory per block: %zu B\n",
prop.maxThreadsPerBlock, prop.sharedMemPerBlock);
// 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 failures = 0;
if (prop.warpSize != kWarpSize) {
std::fprintf(stderr,
"this card reports warp size %d; the race row of the "
"table below is read against %d\n",
prop.warpSize, kWarpSize);
++failures;
}
// The work model, evaluated by the compiler rather than measured. It is
// the claim the lesson makes about Hillis-Steele and it needs no GPU.
std::printf("\nWork per block, computed at compile time\n");
std::printf("%8s %8s %12s %12s\n", "threads", "steps", "scan adds",
"sequential");
std::printf("%8s %8s %12s %12s\n", "-------", "-----", "---------",
"----------");
for (int s = 0; s < kSweep; ++s) {
const int threads = kBlockSizes[s];
std::printf("%8d %8d %12d %12d\n", threads, intLog2(threads),
hillisSteeleAdds(threads), sequentialAdds(threads));
}
const size_t bytes = kElems * sizeof(float);
std::vector<float> h_in(kElems);
for (size_t i = 0; i < kElems; ++i) {
h_in[i] = static_cast<float>(i % 8) + 1.0f;
}
std::vector<float> h_out(kElems);
std::vector<float> h_want(kElems);
float* d_in = nullptr;
float* d_out = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), bytes, cudaMemcpyHostToDevice));
const char* names[kKernels] = {"inclusive, in place", "inclusive, no copy",
"inclusive, double buffered",
"exclusive, double buffered"};
std::printf("\n%zu elements, in[i] = (i %% 8) + 1, one tile per block\n",
kElems);
std::printf("%7s %-27s %11s %11s %6s\n", "threads", "kernel", "mismatches",
"first bad", "gated");
std::printf("%7s %-27s %11s %11s %6s\n", "-------",
"---------------------------", "----------", "---------",
"-----");
int rows = 0;
int expectedRows = 0;
for (int s = 0; s < kSweep; ++s) {
const int threads = kBlockSizes[s];
if (threads > prop.maxThreadsPerBlock) {
std::printf("%7d skipped: this card caps a block at %d threads\n",
threads, prop.maxThreadsPerBlock);
continue;
}
expectedRows += kKernels;
const int blocks =
static_cast<int>((kElems + static_cast<size_t>(threads) - 1) /
static_cast<size_t>(threads));
for (int kernel = 0; kernel < kKernels; ++kernel) {
if (kernel == kKernels - 1) {
scanExclusiveCpu(h_in.data(), h_want.data(), kElems,
static_cast<size_t>(threads));
} else {
scanInclusiveCpu(h_in.data(), h_want.data(), kElems,
static_cast<size_t>(threads));
}
launchKernel(kernel, blocks, threads, d_in, d_out, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
size_t firstBad = kElems;
const size_t bad =
countMismatches(h_out.data(), h_want.data(), kElems, &firstBad);
char firstText[24];
if (bad == 0) {
std::snprintf(firstText, sizeof(firstText), "%s", "-");
} else {
std::snprintf(firstText, sizeof(firstText), "%zu", firstBad);
}
// Only the two correct kernels are gated. Asserting that the
// other two fail would make the run unfalsifiable, and one of
// them fails through a data race, whose result the standard does
// not promise from one launch to the next.
const int gated = (kernel >= 2) ? 1 : 0;
std::printf("%7d %-27s %11zu %11s %6s\n", threads, names[kernel],
bad, firstText, gated ? "yes" : "no");
++rows;
if (gated != 0 && bad != 0) {
std::fprintf(stderr,
"%s at %d threads: %zu of %zu elements wrong, "
"first at %zu, got %.9g want %.9g\n",
names[kernel], threads, bad, kElems, firstBad,
static_cast<double>(h_out[firstBad]),
static_cast<double>(h_want[firstBad]));
++failures;
}
}
}
// A real branch rather than an assert(). CI builds Release, Release
// defines NDEBUG, and NDEBUG deletes assert(), so the check would be
// missing from exactly the build that matters. This one exists so the
// lesson's table cannot quietly disagree with the program's.
if (rows != expectedRows) {
std::fprintf(stderr, "printed %d rows, expected %d\n", rows,
expectedRows);
++failures;
}
if (expectedRows == 0) {
std::fprintf(stderr, "no block size in the sweep fits this card\n");
++failures;
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_out));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
std::printf(
"\nboth double-buffered kernels matched the reference at "
"every block size\n");
return EXIT_SUCCESS;
}