code/day26-atomics/atomics.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 26: atomics and contention, measured.
//
// One sum, three ways, five sizes.
//
// 1. sumGlobalAtomic one atomicAdd per element into one global float.
// 2. sumSharedAtomic one atomicAdd per element into one shared float per
// block, then one global atomicAdd per block.
// 3. sumTreeAtomic day 24's shared-memory tree reduction inside the
// block, then one global atomicAdd per block.
//
// All three read every element exactly once and in the same coalesced order,
// so the bytes read are identical across the three and the only variable is
// where the additions land.
//
// Two data sets, on purpose. The timing sweep uses values in {0, 1, 2, 3},
// whose every partial sum is an integer below 2^24 and therefore exact in
// float in whatever order the additions happen. That lets the correctness
// gate be bitwise equality instead of a tolerance, and it separates
// contention from precision so the timing table measures one thing. The
// determinism section then refills the same buffer with values that are not
// exactly representable and runs each kernel kDetRuns times, which is where
// the ordering shows up as different answers. Day 68 is where that gets
// fixed rather than measured.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o atomics atomics.cu
// Run: ./atomics
//
// Verified: run on a Tesla T4 (compute capability 7.5), driver 595.84,
// CUDA 12.6 (V12.6.85), on 2026-08-30. The only output that may be
// published as this program's output is the transcript in
// evidence/run-2026-08-30.txt. See research/REVIEW-PROCESS.md.
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#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)
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarpSize = 32; // 32 on every GPU this course targets
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr int kDetRuns = 20;
// One accumulator slot per launch timeKernel makes, warm-ups included. The
// count is a property of that helper's body rather than of these two
// constants, so the slot index below is taken modulo this value, which keeps
// the write in the buffer whatever timeKernel does, and checkLaunchCount()
// then reports the disagreement instead of a run silently reusing a slot.
constexpr int kLaunchSlots = kWarmupRuns + kTimedRuns;
// The sweep. 2^12 is small enough that 16 blocks cannot fill the card and the
// launch is most of the run; 2^22 is long enough that the accumulator, not
// the launch, is what the timer sees.
constexpr size_t kSizes[] = {4096, 65536, 262144, 1048576, 4194304};
constexpr int kNumSizes = 5;
constexpr size_t kMaxElems = 4194304;
// Sweep values are i % kExactModulus, so the largest one is
// kExactModulus - 1. Determinism values are 1.0f plus a thousandth of
// i % kInexactModulus, which lands between 1.0 and 1.999 and is almost never
// exactly representable.
constexpr size_t kExactModulus = 4;
constexpr size_t kInexactModulus = 1000;
// Every integer below 2^24 is exact in float, and so is the sum of two of
// them when the result is also below 2^24.
constexpr size_t kFloatIntegerLimit = 1ull << 24;
static_assert(static_cast<int>(sizeof(kSizes) / sizeof(kSizes[0])) == kNumSizes,
"kNumSizes must match the sweep table");
static_assert(kThreadsPerBlock % kWarpSize == 0,
"block size must be a whole number of warps");
static_assert((kThreadsPerBlock & (kThreadsPerBlock - 1)) == 0,
"sumTreeAtomic halves the stride, so the block size must be a "
"power of two");
static_assert(kSizes[kNumSizes - 1] == kMaxElems,
"the input buffer is sized for the largest sweep entry");
// The claim the sweep's correctness gate rests on, checked by the compiler
// rather than asserted in the lesson's prose. If someone raises kMaxElems or
// widens the value range, this build stops instead of the gate quietly
// turning into a test of float precision.
static_assert((kExactModulus - 1) * kMaxElems < kFloatIntegerLimit,
"bitwise equality only holds while every partial sum is an "
"integer below 2^24");
// Adds one element to a single global accumulator, one atomicAdd per thread.
//
// One thread: reads in[i] and adds it to *total.
//
// One warp: the 32 loads cover 128 contiguous bytes, so the read is coalesced
// and identical to the other two kernels here. The 32 atomic adds all name
// the same address, which is the most contention a warp can create on one
// location. What the hardware actually receives may be fewer than 32 adds:
// nvcc aggregates atomics within a warp in many cases, which is the effect
// this program exists to bound.
//
// Launch assumption: gridDim.x * blockDim.x >= n. The guard covers the work,
// and there is no early return, which is the form the other two kernels need
// because they carry barriers.
// snippet: global-atomic
__global__ void sumGlobalAtomic(const float* __restrict__ in,
float* __restrict__ total, size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
atomicAdd(total, in[i]);
}
}
// end snippet
// Accumulates into one shared float per block, then one global atomicAdd.
//
// One thread: adds in[i] to the block's shared accumulator; thread 0 then
// adds the block's total to the global one.
//
// One warp: the same coalesced 128-byte read. The 32 atomic adds now name a
// shared address, and the grid's traffic to the global accumulator falls from
// n adds to one per block.
//
// Barrier rule: every thread reaches both __syncthreads(). The guards cover
// the store and the load, never the barrier, and nothing returns above one.
// The first barrier is the one the r/CUDA post in the lesson asks about. An
// atomic orders accesses to one address and does nothing else, so it will not
// publish thread 0's store to blockSum to the rest of the block.
// snippet: shared-atomic
__global__ void sumSharedAtomic(const float* __restrict__ in,
float* __restrict__ total, size_t n) {
__shared__ float blockSum;
const unsigned int tid = threadIdx.x;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
if (tid == 0) {
blockSum = 0.0f;
}
__syncthreads();
if (i < n) {
atomicAdd(&blockSum, in[i]);
}
__syncthreads();
if (tid == 0) {
atomicAdd(total, blockSum);
}
}
// end snippet
// Day 24's tree reduction inside the block, then one global atomicAdd.
//
// One thread: stages its element in shared memory, takes part in
// log2(kThreadsPerBlock) halving steps, and thread 0 adds the block's total
// to the global one.
//
// One warp: the same coalesced 128-byte read. There is no atomic anywhere
// inside the block, and the grid performs one global atomic per block, the
// same count as sumSharedAtomic. So the only difference between this kernel
// and that one is how a block folds 256 values into one.
//
// Launch assumption: exactly kThreadsPerBlock threads per block. The shared
// array is sized from that constant and the halving loop counts down from it,
// so a different block size would reduce the wrong range rather than fail.
// snippet: tree-atomic
__global__ void sumTreeAtomic(const float* __restrict__ in,
float* __restrict__ total, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
tile[tid] = (i < n) ? in[i] : 0.0f;
__syncthreads();
for (unsigned int half = kThreadsPerBlock / 2; half > 0; half /= 2) {
if (tid < half) {
tile[tid] += tile[tid + half];
}
__syncthreads();
}
if (tid == 0) {
atomicAdd(total, tile[0]);
}
}
// end snippet
// CPU reference. Written for obvious correctness, not speed: one plain loop,
// no OpenMP, no blocking, no intrinsics. It never allocates; the caller owns
// the buffer.
//
// It accumulates in double with Kahan compensation, so it produces the right
// answer rather than one of the many answers a float summation order can
// produce. Day 68 is the one lesson that inverts this and asks both sides to
// use the same order.
static double sumCpu(const float* in, size_t n) {
double total = 0.0;
double comp = 0.0;
for (size_t i = 0; i < n; ++i) {
const double y = static_cast<double>(in[i]) - comp;
const double t = total + y;
comp = (t - total) - y;
total = t;
}
return total;
}
// 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;
}
// Counts the input stream and nothing else: n elements read once. The atomic
// traffic is deliberately outside this figure, because that traffic is the
// thing the three kernels differ in, and folding it into the denominator
// would hide the comparison the table exists to make.
static double readGBs(size_t n, float ms) {
const double bytes = static_cast<double>(n) * sizeof(float);
return bytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// Reports whether the answer is bit for bit the expected one. A real branch
// returning a failure count, never an assert(): CI builds Release, Release
// defines NDEBUG, and NDEBUG deletes assert(), so the check would be missing
// from exactly the build that matters.
static int checkExact(const char* name, size_t n, float got, float want) {
if (got != want) {
std::fprintf(stderr, "%s wrong at n = %zu: got %.9g, want %.9g\n", name,
n, static_cast<double>(got), static_cast<double>(want));
return 1;
}
return 0;
}
// Reports whether timeKernel launched exactly as many times as d_totals has
// slots. It is a real branch and not an assert, for the reason above
// checkExact. A mismatch means the timing helper and this buffer have drifted
// apart, so some launch reused a slot another launch had already accumulated
// into and the timing it produced is not the timing of one clean run.
static int checkLaunchCount(const char* name, int launches) {
if (launches != kLaunchSlots) {
std::fprintf(stderr,
"%s made %d launches per timing but d_totals holds %d "
"slots; timeKernel and the buffer have drifted apart\n",
name, launches, kLaunchSlots);
return 1;
}
return 0;
}
// How many different bit patterns appear in the first `count` results.
// Compares bits rather than values, because two answers that print the same
// at nine digits can still differ, and whether they differ at all is the
// question this section asks.
static int distinctResults(const float* results, int count) {
int distinct = 0;
for (int a = 0; a < count; ++a) {
unsigned int bitsA = 0;
std::memcpy(&bitsA, &results[a], sizeof(bitsA));
bool seen = false;
for (int b = 0; b < a; ++b) {
unsigned int bitsB = 0;
std::memcpy(&bitsB, &results[b], sizeof(bitsB));
if (bitsA == bitsB) {
seen = true;
}
}
if (!seen) {
++distinct;
}
}
return distinct;
}
// Prints one determinism row: how many different answers the same kernel gave
// over the same input, and how far apart the extremes were.
static void reportSpread(const char* name, const float* results, int count) {
float lo = results[0];
float hi = results[0];
for (int r = 1; r < count; ++r) {
if (results[r] < lo) {
lo = results[r];
}
if (results[r] > hi) {
hi = results[r];
}
}
std::printf("%-16s %10d %18.9g %18.9g %14.6g\n", name,
distinctResults(results, count), static_cast<double>(lo),
static_cast<double>(hi), static_cast<double>(hi - lo));
}
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("SMs: %d, threads per block: %d\n", prop.multiProcessorCount,
kThreadsPerBlock);
// 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.
int failures = 0;
int rows = 0;
float* d_in = nullptr;
float* d_total = nullptr;
float* d_totals = nullptr;
const size_t inBytes = kMaxElems * sizeof(float);
const size_t totalsCount = static_cast<size_t>(kLaunchSlots);
const size_t totalsBytes = totalsCount * sizeof(float);
CUDA_CHECK(cudaMalloc(&d_in, inBytes));
CUDA_CHECK(cudaMalloc(&d_total, sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_totals, totalsBytes));
// Values 0, 1, 2 and 3. Every partial sum is then an integer below 2^24,
// so float addition is exact in whatever order it happens, and the gate
// below can be == rather than a tolerance nobody can check.
std::vector<float> h_in(kMaxElems);
for (size_t i = 0; i < kMaxElems; ++i) {
h_in[i] = static_cast<float>(i % kExactModulus);
}
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), inBytes, cudaMemcpyHostToDevice));
std::printf(
"\nSweep: one sum three ways, %d timed runs after %d warm-ups\n",
kTimedRuns, kWarmupRuns);
std::printf("%10s %11s %11s %11s %9s %9s %9s %9s\n", "n", "global ms",
"shared ms", "tree ms", "gl GB/s", "sh GB/s", "tr GB/s",
"fastest");
std::printf("%10s %11s %11s %11s %9s %9s %9s %9s\n", "---------",
"----------", "----------", "----------", "--------",
"--------", "--------", "--------");
for (int s = 0; s < kNumSizes; ++s) {
const size_t n = kSizes[s];
const int blocks =
static_cast<int>((n + kThreadsPerBlock - 1) / kThreadsPerBlock);
const float want = static_cast<float>(sumCpu(h_in.data(), n));
float got = 0.0f;
// One untimed launch first, checked, so a wrong answer is reported
// before any number that came from it reaches the table.
//
// Each timed launch then gets its own slot in d_totals. An
// accumulator is updated in place, and the site's timing rule says an
// in-place kernel cannot go in a timing loop without a reset between
// runs; a fresh slot per launch is that reset, and it sits outside
// the timed region because the pointer arithmetic is on the host.
CUDA_CHECK(cudaMemset(d_total, 0, sizeof(float)));
sumGlobalAtomic<<<blocks, kThreadsPerBlock>>>(d_in, d_total, n);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(&got, d_total, sizeof(float), cudaMemcpyDeviceToHost));
failures += checkExact("sumGlobalAtomic", n, got, want);
CUDA_CHECK(cudaMemset(d_totals, 0, totalsBytes));
int globalSlot = 0;
const float globalMs = timeKernel([&] {
sumGlobalAtomic<<<blocks, kThreadsPerBlock>>>(
d_in, d_totals + (globalSlot % kLaunchSlots), n);
++globalSlot;
});
failures += checkLaunchCount("sumGlobalAtomic", globalSlot);
CUDA_CHECK(cudaMemset(d_total, 0, sizeof(float)));
sumSharedAtomic<<<blocks, kThreadsPerBlock>>>(d_in, d_total, n);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(&got, d_total, sizeof(float), cudaMemcpyDeviceToHost));
failures += checkExact("sumSharedAtomic", n, got, want);
CUDA_CHECK(cudaMemset(d_totals, 0, totalsBytes));
int sharedSlot = 0;
const float sharedMs = timeKernel([&] {
sumSharedAtomic<<<blocks, kThreadsPerBlock>>>(
d_in, d_totals + (sharedSlot % kLaunchSlots), n);
++sharedSlot;
});
failures += checkLaunchCount("sumSharedAtomic", sharedSlot);
CUDA_CHECK(cudaMemset(d_total, 0, sizeof(float)));
sumTreeAtomic<<<blocks, kThreadsPerBlock>>>(d_in, d_total, n);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(&got, d_total, sizeof(float), cudaMemcpyDeviceToHost));
failures += checkExact("sumTreeAtomic", n, got, want);
CUDA_CHECK(cudaMemset(d_totals, 0, totalsBytes));
int treeSlot = 0;
const float treeMs = timeKernel([&] {
sumTreeAtomic<<<blocks, kThreadsPerBlock>>>(
d_in, d_totals + (treeSlot % kLaunchSlots), n);
++treeSlot;
});
failures += checkLaunchCount("sumTreeAtomic", treeSlot);
// The winner is reported, never gated. Which strategy is fastest at
// which size is the question the run exists to answer, so encoding an
// answer as a check would make the test unfalsifiable.
const char* fastest = "global";
if (sharedMs < globalMs && sharedMs <= treeMs) {
fastest = "shared";
} else if (treeMs < globalMs && treeMs < sharedMs) {
fastest = "tree";
}
std::printf("%10zu %11.4f %11.4f %11.4f %9.1f %9.1f %9.1f %9s\n", n,
globalMs, sharedMs, treeMs, readGBs(n, globalMs),
readGBs(n, sharedMs), readGBs(n, treeMs), fastest);
++rows;
}
std::printf(
"\nEvery row above read the same %zu bytes per element and returned\n"
"the same answer bit for bit, checked against a Kahan double\n"
"reference with == and not a tolerance. Only the accumulation\n"
"strategy changed.\n",
sizeof(float));
// The second data set. These values are not exactly representable, so the
// order the additions happen in now changes the answer. Same buffer, same
// three kernels, one line of host code different.
for (size_t i = 0; i < kMaxElems; ++i) {
h_in[i] = 1.0f + static_cast<float>(i % kInexactModulus) * 0.001f;
}
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), inBytes, cudaMemcpyHostToDevice));
const size_t detN = kMaxElems;
const int detBlocks =
static_cast<int>((detN + kThreadsPerBlock - 1) / kThreadsPerBlock);
const double detWant = sumCpu(h_in.data(), detN);
float globalRuns[kDetRuns];
float sharedRuns[kDetRuns];
float treeRuns[kDetRuns];
for (int r = 0; r < kDetRuns; ++r) {
CUDA_CHECK(cudaMemset(d_total, 0, sizeof(float)));
sumGlobalAtomic<<<detBlocks, kThreadsPerBlock>>>(d_in, d_total, detN);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(&globalRuns[r], d_total, sizeof(float),
cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaMemset(d_total, 0, sizeof(float)));
sumSharedAtomic<<<detBlocks, kThreadsPerBlock>>>(d_in, d_total, detN);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(&sharedRuns[r], d_total, sizeof(float),
cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaMemset(d_total, 0, sizeof(float)));
sumTreeAtomic<<<detBlocks, kThreadsPerBlock>>>(d_in, d_total, detN);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(&treeRuns[r], d_total, sizeof(float),
cudaMemcpyDeviceToHost));
}
std::printf("\nDeterminism at n = %zu, %d runs each, input that is not\n",
detN, kDetRuns);
std::printf("exactly representable. Kahan double reference: %.9g\n",
detWant);
std::printf("%-16s %10s %18s %18s %14s\n", "kernel", "distinct", "min",
"max", "max - min");
std::printf("%-16s %10s %18s %18s %14s\n", "----------------", "--------",
"----------------", "----------------", "------------");
reportSpread("global atomic", globalRuns, kDetRuns);
reportSpread("shared atomic", sharedRuns, kDetRuns);
reportSpread("tree + atomic", treeRuns, kDetRuns);
rows += 3;
// A wrong answer here would be a bug, but a spread is not a bug, so the
// gate is a wide tolerance around the reference and never a claim about
// how many distinct answers appeared. The bound is the site's reduction
// rule, 4 * 2^-23 * sqrt(K) relative, which is the uncorrelated-rounding
// estimate rather than the vacuous worst case.
const double detRtol =
4.0 * 1.1920929e-7 * std::sqrt(static_cast<double>(detN));
const double detTol = detRtol * std::fabs(detWant);
for (int r = 0; r < kDetRuns; ++r) {
const double diff = static_cast<double>(globalRuns[r]) - detWant;
if (std::fabs(diff) > detTol) {
std::fprintf(stderr,
"global atomic run %d is %.9g, outside %.6g of the "
"reference %.9g\n",
r, static_cast<double>(globalRuns[r]), detTol,
detWant);
++failures;
break;
}
}
// The row count is checked so the lesson's tables cannot drift from what
// the program prints. A real branch, not an assert, for the reason given
// above checkExact.
const int kExpectedRows = kNumSizes + 3;
if (rows != kExpectedRows) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's tables and "
"this program disagree\n",
rows, kExpectedRows);
++failures;
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_total));
CUDA_CHECK(cudaFree(d_totals));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}