code/day27-memory-model/memory_model.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 27: fences, scopes and atomics in the CUDA memory model.
//
// Four kernels in two pairs. sumLastBlockRacy and sumLastBlockFenced are the
// same last-block-standing reduction and differ by one line, a
// __threadfence() between the store of a block's partial sum and the atomic
// that announces the store. tallyBlocksDeviceScope and tallyBlocksBlockScope
// are the same one-line counter and differ by the scope suffix on the
// atomic.
//
// Nothing here is timed. The subject is ordering and scope, so the file has
// no cudaEvent in it and the output carries no number a stopwatch produced.
//
// The two undefined kernels are reported and never gated. A race is allowed
// to produce the right answer, and a test demanding that it produce the
// wrong one would be asserting that undefined behaviour is dependable. The
// gates are on the two defined kernels, which have to be exact on every run.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o memory_model memory_model.cu
// Run: ./memory_model
//
// 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 <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. The built-in `warpSize` is a run-time
// value, so it cannot size an array or appear in a static_assert; this
// constant can do both.
constexpr int kWarpSize = 32;
// The course default block size, eight warps. The block reduction below
// halves the stride, so it also has to be a power of two.
constexpr int kThreadsPerBlock = 256;
// A fixed grid rather than one derived from the element count, because
// gridDim.x is the number this lesson is about: it is the ticket count, the
// length of the partials buffer, and the size of the window in which the
// race can happen. 1024 blocks against 40 SMs on the verification card means
// most blocks retire long before the last one starts reading.
constexpr int kBlocks = 1024;
// Every element is 1.0f, so the reference sum is the element count and every
// partial sum is a whole number well under 2^24. That is what lets this
// program compare exactly instead of within a tolerance, and an exact
// comparison is the only signal it has: a stale partial read costs a whole
// block's worth of elements, not a last bit.
constexpr size_t kElems = 4ull * 1024ull * 1024ull + 611ull;
// One launch proves nothing about a race. The number worth printing is how
// many launches out of many disagreed.
constexpr int kRepeats = 100;
static_assert(kThreadsPerBlock % kWarpSize == 0,
"block size must be a whole number of warps");
static_assert((kThreadsPerBlock & (kThreadsPerBlock - 1)) == 0,
"reduceBlock halves the stride, so the block size must be a "
"power of two");
static_assert(kThreadsPerBlock <= 1024,
"1024 threads per block is the hardware ceiling at every "
"compute capability this course targets");
static_assert(kBlocks > 1,
"the last-block pattern needs a block that is not the only "
"block, or there is no hand-off to get wrong");
static_assert(kElems < (1ull << 24),
"float represents every integer below 2^24 exactly, which is "
"what makes the comparison in this program exact");
static_assert(kElems >= static_cast<size_t>(kBlocks) * kThreadsPerBlock,
"every block must own at least one element, or it stores no "
"partial sum and the last block reads a cell nobody filled");
static_assert(kElems % (static_cast<size_t>(kBlocks) * kThreadsPerBlock) != 0,
"the grid-stride loop's last pass must be partial, so the loop "
"condition runs as a bounds check on every launch");
static_assert(kRepeats >= 2,
"a racy kernel that ran once has told you nothing");
// Reduces one value per thread into tile[0] and returns it to every thread.
//
// Every thread reaches every barrier here. The caller has already turned its
// out-of-range elements into contributions of zero, so there is no guard and
// no return anywhere above a __syncthreads().
//
// Launch assumption: exactly kThreadsPerBlock threads per block, and
// kThreadsPerBlock a power of two, both pinned by the static_asserts above.
// The caller may reuse `tile` afterwards only across a barrier, because the
// read of tile[0] below is the last access of this call.
__device__ float reduceBlock(float* tile, float value) {
const unsigned int tid = threadIdx.x;
tile[tid] = value;
__syncthreads();
for (unsigned int half = kThreadsPerBlock / 2u; half > 0u; half /= 2u) {
if (tid < half) {
tile[tid] += tile[tid + half];
}
__syncthreads();
}
return tile[0];
}
// Sums `in` in one launch and gets the hand-off between blocks wrong.
//
// One thread walks the array with a grid-stride loop and adds what it
// touches. The block reduces those values in shared memory, thread 0 stores
// the block's partial sum, and thread 0 then takes a ticket from an atomic
// counter. The block holding the last ticket sums all gridDim.x partial sums
// and writes the answer.
//
// Memory: one warp's 32 lanes take 32 consecutive elements on every pass of
// the grid-stride loop, so each pass is one contiguous 128-byte span per
// warp. The store to `partials` is one float per block, and the last block's
// read of `partials` is one contiguous span per warp as well.
//
// No pointer here carries __restrict__, and that is deliberate rather than an
// oversight: `partials` is written by one block and read by another through
// the same parameter, so the exclusive-access promise is not the promise this
// kernel wants to hand the compiler. CUDA-CODE-STYLE.md allows __restrict__
// on every pointer in a signature or on none, and this is the none case.
//
// Launch assumption: exactly kThreadsPerBlock threads per block, gridDim.x
// equal to the length of `partials`, and `counter` zeroed before the launch.
// Nothing separates the store to `partials` from the atomicAdd that announces
// it, which is the bug this lesson exists to show. Do not copy it.
__global__ void sumLastBlockRacy(const float* in, float* partials, float* total,
unsigned int* counter, size_t n) {
__shared__ float tile[kThreadsPerBlock];
__shared__ bool isLastBlockDone;
const unsigned int tid = threadIdx.x;
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float sum = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
sum += in[i];
}
const float blockSum = reduceBlock(tile, sum);
// snippet: racy-tail
if (tid == 0) {
partials[blockIdx.x] = blockSum;
// The __threadfence() belongs on this line.
const unsigned int ticket = atomicAdd(counter, 1u);
isLastBlockDone = (ticket == gridDim.x - 1u);
}
__syncthreads();
// end snippet: racy-tail
// isLastBlockDone is the same value for every thread in the block, so
// the barriers inside reduceBlock below are reached by all of them or by
// none of them. Day 14 is where that rule comes from.
if (isLastBlockDone) {
float last = 0.0f;
for (unsigned int b = tid; b < gridDim.x; b += blockDim.x) {
last += partials[b];
}
const float grandSum = reduceBlock(tile, last);
if (tid == 0) {
total[0] = grandSum;
}
}
}
// The same kernel with the fence put in. One line, and it is the whole
// difference between the two rows of the first table.
//
// One thread does what it does in sumLastBlockRacy. The fence sits between
// the store that publishes this block's partial sum and the atomic that
// announces the store is there, so no other block can see the ticket before
// it can see the number the ticket is promising.
//
// Memory: identical to sumLastBlockRacy. Same loads, same stores, same
// addresses, same access pattern.
//
// Launch assumption: identical to sumLastBlockRacy.
__global__ void sumLastBlockFenced(const float* in, float* partials,
float* total, unsigned int* counter,
size_t n) {
__shared__ float tile[kThreadsPerBlock];
__shared__ bool isLastBlockDone;
const unsigned int tid = threadIdx.x;
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float sum = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
sum += in[i];
}
const float blockSum = reduceBlock(tile, sum);
// snippet: fenced-tail
if (tid == 0) {
partials[blockIdx.x] = blockSum;
__threadfence();
const unsigned int ticket = atomicAdd(counter, 1u);
isLastBlockDone = (ticket == gridDim.x - 1u);
}
__syncthreads();
// end snippet: fenced-tail
if (isLastBlockDone) {
float last = 0.0f;
for (unsigned int b = tid; b < gridDim.x; b += blockDim.x) {
last += partials[b];
}
const float grandSum = reduceBlock(tile, last);
if (tid == 0) {
total[0] = grandSum;
}
}
}
// Adds 1 to one global counter from every thread, at device scope.
//
// One thread does one atomic add. The counter has to finish at exactly
// gridDim.x * blockDim.x, on every run, on every card.
//
// Memory: every lane of every warp hits the same address, which is the
// worst contention case there is. Day 26 measures what that costs; this
// kernel only checks the number.
//
// Launch assumption: `tally` zeroed before the launch, and the product of the
// launch dimensions inside the range of unsigned int.
// snippet: tally-device
__global__ void tallyBlocksDeviceScope(unsigned int* tally) {
atomicAdd(tally, 1u);
}
// end snippet: tally-device
// The same line with a scope suffix on it, which is undefined across blocks.
//
// One thread does one atomic add, atomic at cuda::thread_scope_block. The 256
// threads inside a block cannot lose an increment to each other. The 1024
// blocks racing on one device-scope address can, and the programming guide
// calls that a data race with undefined behaviour.
//
// Memory: identical to tallyBlocksDeviceScope. One address, every thread.
//
// Launch assumption: identical to tallyBlocksDeviceScope. This kernel is
// reported and never gated, and its worst outcome is a wrong number rather
// than a wedged device, which is why it is safe to run at all.
// snippet: tally-block
__global__ void tallyBlocksBlockScope(unsigned int* tally) {
atomicAdd_block(tally, 1u);
}
// end snippet: tally-block
// CPU reference for both sum kernels. Accumulates in double even though the
// kernels accumulate in float, per the house rule, though on this input both
// are exact. Written for obvious correctness, not speed. It never allocates;
// the caller owns the buffer.
static double sumCpu(const float* x, size_t n) {
double total = 0.0;
for (size_t i = 0; i < n; ++i) {
total += static_cast<double>(x[i]);
}
return total;
}
// Runs one of the two sum kernels kRepeats times, prints its row, and returns
// how many of those runs produced a total that was not the reference.
//
// All three device buffers are reset before every launch. Zeroing `partials`
// is the part that matters: left alone it would still hold the previous run's
// correct values, so a stale read would return very nearly the right number
// and hide the thing this program exists to show.
//
// The comparison is `!=` on floats rather than a tolerance because every
// value in this program is a whole number below 2^24, where float is exact.
// A difference here is an ordering bug, not a rounding one.
static int runSum(bool fenced, const float* d_in, float* d_partials,
float* d_total, unsigned int* d_counter, float want) {
const char* name = fenced ? "sumLastBlockFenced" : "sumLastBlockRacy";
int wrong = 0;
int firstBadRun = -1;
float firstBadTotal = 0.0f;
for (int run = 0; run < kRepeats; ++run) {
CUDA_CHECK(cudaMemset(d_partials, 0, kBlocks * sizeof(float)));
CUDA_CHECK(cudaMemset(d_total, 0, sizeof(float)));
CUDA_CHECK(cudaMemset(d_counter, 0, sizeof(unsigned int)));
if (fenced) {
sumLastBlockFenced<<<kBlocks, kThreadsPerBlock>>>(
d_in, d_partials, d_total, d_counter, kElems);
} else {
sumLastBlockRacy<<<kBlocks, kThreadsPerBlock>>>(
d_in, d_partials, d_total, d_counter, kElems);
}
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
float got = 0.0f;
CUDA_CHECK(
cudaMemcpy(&got, d_total, sizeof(float), cudaMemcpyDeviceToHost));
if (got != want) {
if (wrong == 0) {
firstBadRun = run;
firstBadTotal = got;
}
++wrong;
}
}
if (wrong == 0) {
std::printf("%-20s %5d %6d %13s %14s\n", name, kRepeats, wrong, "-",
"-");
} else {
std::printf("%-20s %5d %6d %13d %14.9g\n", name, kRepeats, wrong,
firstBadRun, firstBadTotal);
}
return wrong;
}
// Runs one of the two tally kernels kRepeats times, prints its row, and
// returns how many runs did not land on the exact expected count.
//
// The lowest count seen across the runs is printed as well as the number of
// wrong runs, because an atomic that is not atomic at the scope it is used at
// undercounts rather than failing, and how far it undercounts is the
// interesting half.
static int runTally(bool deviceScope, unsigned int* d_tally) {
const char* name =
deviceScope ? "tallyBlocksDeviceScope" : "tallyBlocksBlockScope";
const unsigned int want =
static_cast<unsigned int>(kBlocks) * kThreadsPerBlock;
int wrong = 0;
unsigned int lowest = want;
for (int run = 0; run < kRepeats; ++run) {
CUDA_CHECK(cudaMemset(d_tally, 0, sizeof(unsigned int)));
if (deviceScope) {
tallyBlocksDeviceScope<<<kBlocks, kThreadsPerBlock>>>(d_tally);
} else {
tallyBlocksBlockScope<<<kBlocks, kThreadsPerBlock>>>(d_tally);
}
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
unsigned int got = 0;
CUDA_CHECK(cudaMemcpy(&got, d_tally, sizeof(unsigned int),
cudaMemcpyDeviceToHost));
if (got != want) {
++wrong;
}
if (got < lowest) {
lowest = got;
}
}
std::printf("%-24s %5d %6d %13u %10u\n", name, kRepeats, wrong, lowest,
want);
return wrong;
}
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);
const size_t bytes = kElems * sizeof(float);
std::printf("n = %zu floats, every one of them 1.0f\n", kElems);
std::printf("%d blocks x %d threads, %d launches per kernel\n", kBlocks,
kThreadsPerBlock, kRepeats);
std::printf("no timing anywhere in this program\n");
std::vector<float> h_in(kElems, 1.0f);
const float want = static_cast<float>(sumCpu(h_in.data(), kElems));
float* d_in = nullptr;
float* d_partials = nullptr;
float* d_total = nullptr;
unsigned int* d_counter = nullptr;
unsigned int* d_tally = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d_partials, kBlocks * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_total, sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_counter, sizeof(unsigned int)));
CUDA_CHECK(cudaMalloc(&d_tally, sizeof(unsigned int)));
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), bytes, cudaMemcpyHostToDevice));
// One status variable and one exit. Every path below falls through to the
// frees at the bottom, including the failing ones, because a lesson that
// returns early past its own cudaFree teaches the opposite of day 5.
int status = EXIT_SUCCESS;
std::printf("\nreference total %.9g\n", want);
std::printf("\n%-20s %5s %6s %13s %14s\n", "kernel", "runs", "wrong",
"first bad run", "total there");
const int racyWrong =
runSum(false, d_in, d_partials, d_total, d_counter, want);
const int fencedWrong =
runSum(true, d_in, d_partials, d_total, d_counter, want);
std::printf("\n%-24s %5s %6s %13s %10s\n", "tally kernel", "runs",
"wrong", "lowest count", "want");
const int deviceWrong = runTally(true, d_tally);
const int blockWrong = runTally(false, d_tally);
// The racy kernel's count and the block-scope tally's count are reported
// and never gated, for the reason at the top of this file. Printing them
// beside the gated rows is the whole result.
std::printf(
"\nreported, not gated: sumLastBlockRacy wrong on %d of %d "
"runs, tallyBlocksBlockScope on %d of %d\n",
racyWrong, kRepeats, blockWrong, kRepeats);
// The gates, and they are real branches returning EXIT_FAILURE rather
// than assert(), because CI builds Release, Release defines NDEBUG, and
// NDEBUG deletes assert() out of the build that matters.
if (fencedWrong != 0) {
std::fprintf(stderr,
"sumLastBlockFenced disagreed with the reference on %d "
"of %d runs; the fence did not hold\n",
fencedWrong, kRepeats);
status = EXIT_FAILURE;
}
if (deviceWrong != 0) {
std::fprintf(stderr,
"tallyBlocksDeviceScope lost an increment on %d of %d "
"runs; a device-scope atomic must not\n",
deviceWrong, kRepeats);
status = EXIT_FAILURE;
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_partials));
CUDA_CHECK(cudaFree(d_total));
CUDA_CHECK(cudaFree(d_counter));
CUDA_CHECK(cudaFree(d_tally));
if (status != EXIT_SUCCESS) {
std::fprintf(stderr, "a gated kernel produced a wrong answer\n");
return status;
}
std::printf("\nboth gated kernels were exact on every run\n");
return EXIT_SUCCESS;
}