code/day68-reproducibility/reproducibility.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 68: floating point and reproducibility.
//
// Three ways to sum f(x) = x * x + 0.25f over the same 4,194,915 floats,
// run 100 times each, judged on the bit pattern of the result rather than
// on a stopwatch.
//
// sumSquaresAtomic every block adds its partial to one float with
// atomicAdd. Correct, race-free, and the order of
// the 1,024 additions is whatever the scheduler
// produced, so the bits may move run to run.
// sumSquaresKahanAtomic the same kernel with Kahan compensation on the
// per-thread accumulator. Accuracy tool, not a
// determinism tool: the atomic tail is untouched,
// so the bits may still move.
// the two-pass tree pass 1 writes one partial per block, pass 2 sums
// the partials in one block with a fixed tree. No
// atomic anywhere, every addition order pinned by
// the code, so the bits must not move. This is the
// only path the program gates.
//
// The nondeterministic kernels are reported and never gated. Ordering
// nondeterminism is legal (day 27: a relaxed atomic promises atomicity and
// nothing about order), so a check demanding that the bits move would gate
// on scheduler luck. The gates are the tree's single bit pattern across all
// 100 runs and every kernel's distance from a double Kahan reference.
//
// Nothing here is timed. Day 26 already priced these kernels; this file is
// about which bits come back.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o reproducibility \
// reproducibility.cu
// Also: nvcc -std=c++17 -O3 -arch=sm_75 -fmad=false \
// -o reproducibility_nofmad reproducibility.cu
// Run: ./reproducibility && ./reproducibility_nofmad
//
// The second build is day 47's flag pointed at day 68's question: the tree
// hash must be constant within each binary and is expected to differ
// between them, because contraction of x * x + 0.25f is the compiler's
// call and -fmad=false takes the call away.
//
#include <cmath>
#include <cstdint>
#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)
// 32 on every GPU this course targets. The built-in `warpSize` is a
// run-time value, so it cannot appear in a static_assert; this constant can.
constexpr int kWarpSize = 32;
// The course default block size, eight warps. reduceBlock halves the
// stride, so it also has to be a power of two.
constexpr int kThreadsPerBlock = 256;
// A fixed grid, never derived from the element count or the device. The
// grid shape decides which elements meet in which order, so deriving it
// from cudaOccupancyMaxActiveBlocksPerMultiprocessor or from prop
// .multiProcessorCount would make the answer a function of the card. The
// exercise changes this constant to see what that costs.
constexpr int kBlocks = 1024;
// Same element count as day 27, and the same inexact fill as day 26's
// determinism sweep: 1.0f + (i % 1000) * 0.001f is almost never
// exactly representable, so addition order can reach the answer.
constexpr size_t kElems = 4ull * 1024ull * 1024ull + 611ull;
constexpr size_t kInexactModulus = 1000;
// One run proves nothing about determinism. The claim on the page is over
// 100 runs, so the count lives here where the page can cite it.
constexpr int kRepeats = 100;
// The added constant in f(x) = x * x + 0.25f. Exactly representable, so
// the only interesting rounding in f is the multiply and the add
// themselves, which is where -fmad gets a vote.
constexpr float kQuarter = 0.25f;
// Relative distance allowed between any kernel's total and the double
// Kahan reference. Order changes the last few bits, not the neighbourhood:
// at a total near 10^7 this allows an absolute drift about 200x wider than
// the spread day 26 measured on its worst kernel.
constexpr double kRelTolerance = 1e-5;
// FNV-1a 64-bit, folded over the 100 result bit patterns in run order. One
// number per kernel per binary, so two builds can be compared with grep.
constexpr uint64_t kFnvOffset = 14695981039346656037ull;
constexpr uint64_t kFnvPrime = 1099511628211ull;
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 % kThreadsPerBlock == 0,
"pass 2 hands each thread kBlocks / kThreadsPerBlock "
"partials; a remainder would change the tree's shape");
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, "one run cannot show agreement or disagreement");
// The value each thread accumulates. The multiply next to the add is the
// point: nvcc may contract them into one FFMA with one rounding, or keep
// two roundings under -fmad=false, and the two choices give different bits
// for the same source. Day 47 measured that trade; day 68 only needs the
// fact that it exists. A bare sum of in[i] would have no multiply and
// nothing for the flag to change.
// snippet: transform
__device__ float squarePlusQuarter(float x) {
return x * x + kQuarter;
}
// end snippet
// Reduces one value per thread into tile[0] and returns it to thread 0.
// This is day 24's fixed tree: which pairs meet in which round is decided
// by tid and the halving loop, never by the scheduler, which is what makes
// every kernel built on it give the same bits on every run.
//
// Every thread reaches every barrier. The caller has already folded 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, power of
// two, pinned by the static_asserts above.
__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 f(in[i]) into *total with one atomicAdd per block.
//
// One thread walks its grid-stride slice in a fixed order, the block tree
// is fixed, and then thread 0 hands the partial to an atomicAdd. That last
// step is the deliberate day 68 nondeterminism: the 1,024 partials join
// *total in arrival order, arrival order is scheduling, and float addition
// is not associative, so the bits of *total may differ run to run. Not a
// bug in the memory-model sense. Every answer it can produce is legal.
//
// Memory: consecutive threads read consecutive elements, so one warp's 32
// loads cover 128 contiguous bytes per pass.
//
// Launch assumption: kThreadsPerBlock threads per block. *total must be
// zero on entry; the host zeroes it before every run.
__global__ void sumSquaresAtomic(const float* __restrict__ in,
float* __restrict__ total, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float acc = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
acc += squarePlusQuarter(in[i]);
}
const float blockSum = reduceBlock(tile, acc);
// snippet: atomic-tail
if (threadIdx.x == 0) {
// Deliberate (day 68): arrival order decides the addition order.
atomicAdd(total, blockSum);
}
// end snippet
}
// The same kernel with Kahan compensation on the per-thread accumulator.
//
// Kahan is an accuracy tool for a sequential sum whose order you control.
// It cannot touch the atomic tail, because the compensation term needs to
// see the running sum and an atomicAdd never shows it to you. So this
// kernel is expected to land closer to the double reference and still
// return more than one bit pattern across runs: accuracy up, determinism
// unchanged. That gap is the reason the two words are not synonyms.
//
// Memory and launch assumptions: identical to sumSquaresAtomic.
__global__ void sumSquaresKahanAtomic(const float* __restrict__ in,
float* __restrict__ total, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float acc = 0.0f;
float comp = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
const float y = squarePlusQuarter(in[i]) - comp;
const float s = acc + y;
comp = (s - acc) - y;
acc = s;
}
const float blockSum = reduceBlock(tile, acc);
if (threadIdx.x == 0) {
// Deliberate (day 68): same nondeterministic tail as above.
atomicAdd(total, blockSum);
}
}
// Pass 1 of the deterministic path: one partial per block, no atomic.
//
// Same walk and same tree as sumSquaresAtomic; the only change is where
// the partial goes. partials[blockIdx.x] is owned by exactly one block, so
// there is nothing to contend for and no order to leave to the scheduler.
//
// Memory: the store is one float per block. Launch assumption:
// kThreadsPerBlock threads, exactly kBlocks blocks, and partials holds
// kBlocks floats, all of which this launch overwrites.
__global__ void sumSquaresPerBlock(const float* __restrict__ in,
float* __restrict__ partials, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float acc = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
acc += squarePlusQuarter(in[i]);
}
partials[blockIdx.x] = reduceBlock(tile, acc);
}
// Pass 2: one block folds the kBlocks partials with the same fixed tree.
//
// Each thread sums kBlocks / kThreadsPerBlock partials in a fixed strided
// order, then the tree runs once. Every addition in both passes has its
// operands chosen by an index expression, so the whole two-pass sum is one
// fixed parenthesisation of the input and must give the same bits on every
// run of this binary on this card.
//
// Memory: one warp's 32 loads of partials are consecutive. Launch
// assumption: exactly one block of kThreadsPerBlock threads.
// snippet: tree-pass2
__global__ void sumPartials(const float* __restrict__ partials,
float* __restrict__ total) {
__shared__ float tile[kThreadsPerBlock];
float acc = 0.0f;
for (int p = static_cast<int>(threadIdx.x); p < kBlocks;
p += kThreadsPerBlock) {
acc += partials[p];
}
const float treeSum = reduceBlock(tile, acc);
if (threadIdx.x == 0) {
*total = treeSum;
}
}
// end snippet
// Double Kahan reference. Not the "true" answer either, but close enough
// that every float-ordered sum must land within kRelTolerance of it, and
// stable, so all three kernels are judged against the same number.
static double sumSquaresCpu(const float* in, size_t n) {
double acc = 0.0;
double comp = 0.0;
for (size_t i = 0; i < n; ++i) {
const double x = static_cast<double>(in[i]);
const double y = (x * x + 0.25) - comp;
const double s = acc + y;
comp = (s - acc) - y;
acc = s;
}
return acc;
}
// How many different bit patterns a run sequence contains. O(runs^2) and
// runs is 100, so nobody cares.
static int countDistinct(const uint32_t* bits, int runs) {
int distinct = 0;
for (int i = 0; i < runs; ++i) {
bool seen = false;
for (int j = 0; j < i; ++j) {
if (bits[j] == bits[i]) {
seen = true;
break;
}
}
if (!seen) {
++distinct;
}
}
return distinct;
}
// snippet: run-hash
// FNV-1a over the whole run sequence, in order. Two binaries printed the
// same hash for a kernel if and only if all 100 runs matched bit for bit,
// which is the comparison the -fmad=false build exists for.
static uint64_t hashRuns(const uint32_t* bits, int runs) {
uint64_t h = kFnvOffset;
for (int i = 0; i < runs; ++i) {
for (int b = 0; b < 4; ++b) {
h ^= (bits[i] >> (8 * b)) & 0xffu;
h *= kFnvPrime;
}
}
return h;
}
// end snippet
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);
std::printf("n = %zu floats, in[i] = 1.0f + (i %% %zu) * 0.001f\n", kElems,
kInexactModulus);
std::printf(
"summing f(x) = x * x + 0.25f, %d blocks x %d threads, "
"%d runs per kernel\n\n",
kBlocks, kThreadsPerBlock, kRepeats);
std::vector<float> h_in(kElems);
for (size_t i = 0; i < kElems; ++i) {
h_in[i] = 1.0f + static_cast<float>(i % kInexactModulus) * 0.001f;
}
const double reference = sumSquaresCpu(h_in.data(), kElems);
std::printf("reference (double, Kahan): %.6f\n\n", reference);
const size_t inBytes = kElems * sizeof(float);
float* d_in = nullptr;
float* d_partials = nullptr;
float* d_total = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, inBytes));
CUDA_CHECK(cudaMalloc(&d_partials, kBlocks * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_total, sizeof(float)));
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), inBytes, cudaMemcpyHostToDevice));
// Three kernels, kRepeats runs each, one uint32_t bit pattern per run.
// The float value is kept only long enough to score it against the
// reference; the page's subject is the bits.
uint32_t bitsAtomic[kRepeats];
uint32_t bitsKahan[kRepeats];
uint32_t bitsTree[kRepeats];
double maxRelErr[3] = {0.0, 0.0, 0.0};
for (int run = 0; run < kRepeats; ++run) {
for (int kernel = 0; kernel < 3; ++kernel) {
CUDA_CHECK(cudaMemset(d_total, 0, sizeof(float)));
if (kernel == 0) {
sumSquaresAtomic<<<kBlocks, kThreadsPerBlock>>>(d_in, d_total,
kElems);
} else if (kernel == 1) {
sumSquaresKahanAtomic<<<kBlocks, kThreadsPerBlock>>>(
d_in, d_total, kElems);
} else {
sumSquaresPerBlock<<<kBlocks, kThreadsPerBlock>>>(
d_in, d_partials, kElems);
sumPartials<<<1, kThreadsPerBlock>>>(d_partials, d_total);
}
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
float total = 0.0f;
CUDA_CHECK(cudaMemcpy(&total, d_total, sizeof(float),
cudaMemcpyDeviceToHost));
uint32_t bits = 0;
std::memcpy(&bits, &total, sizeof(bits));
if (kernel == 0) {
bitsAtomic[run] = bits;
} else if (kernel == 1) {
bitsKahan[run] = bits;
} else {
bitsTree[run] = bits;
}
const double relErr =
std::fabs(static_cast<double>(total) - reference) /
std::fabs(reference);
if (relErr > maxRelErr[kernel]) {
maxRelErr[kernel] = relErr;
}
}
}
const char* names[3] = {"sumSquaresAtomic", "sumSquaresKahanAtomic",
"two-pass tree"};
const uint32_t* allBits[3] = {bitsAtomic, bitsKahan, bitsTree};
std::printf("%-22s %4s %8s %10s %12s %18s\n", "kernel", "runs",
"distinct", "first bits", "max rel err", "run hash");
for (int kernel = 0; kernel < 3; ++kernel) {
std::printf(
"%-22s %4d %8d 0x%08x %12.3e 0x%016llx\n", names[kernel],
kRepeats, countDistinct(allBits[kernel], kRepeats),
static_cast<unsigned int>(allBits[kernel][0]), maxRelErr[kernel],
static_cast<unsigned long long>(
hashRuns(allBits[kernel], kRepeats)));
}
std::printf(
"\nreported, not gated: the two atomic kernels' distinct "
"counts.\nOrdering nondeterminism is legal, so a run where "
"they all agree\nis a result to report, not a failure.\n\n");
// The gates. Real branches, not assert(): CI builds Release, Release
// defines NDEBUG, and NDEBUG deletes an assert out of exactly the
// build that matters. Failures accumulate so the device memory is
// released in one place on every path.
int failures = 0;
const int treeDistinct = countDistinct(bitsTree, kRepeats);
if (treeDistinct != 1) {
std::fprintf(stderr,
"FAIL: two-pass tree produced %d bit patterns over %d "
"runs; a fixed-shape tree must produce exactly 1\n",
treeDistinct, kRepeats);
++failures;
}
for (int kernel = 0; kernel < 3; ++kernel) {
if (maxRelErr[kernel] > kRelTolerance) {
std::fprintf(stderr,
"FAIL: %s max relative error %.3e exceeds %.0e "
"against the double Kahan reference\n",
names[kernel], maxRelErr[kernel], kRelTolerance);
++failures;
}
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_partials));
CUDA_CHECK(cudaFree(d_total));
if (failures != 0) {
std::fprintf(stderr, "%d gate(s) failed\n", failures);
return EXIT_FAILURE;
}
std::printf(
"gates passed: tree bit pattern constant over %d runs, "
"all kernels within %.0e of the reference\n",
kRepeats, kRelTolerance);
return EXIT_SUCCESS;
}