code/day29-histogram/histogram.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 29: a byte histogram, with and without privatization.
//
// Three kernels over the same buffer, one launch configuration:
//
// sumBytes reads every byte, no atomics at all. The floor.
// histogramGlobal one atomicAdd per byte, straight into global memory.
// histogramPrivate one atomicAdd per byte into the block's own copy in
// shared memory, then one atomic per thread to merge.
//
// Privatization does not remove the atomics. It moves most of them somewhere
// with 255 rivals instead of 262,143, and it cuts the number that reach
// global memory from kBytes to kBlocks * kBins, which is 67,109,475 down to
// 262,144 here.
//
// Three inputs, because contention is a property of the data and not of the
// kernel. The kernels never see which one they got:
//
// uniform bytes a linear congruential generator, all 256 values likely
// english text a paragraph of this project's own prose, repeated
// all one byte every byte the same value. The worst case there is.
//
// On timing a kernel that updates its output in place. CUDA-CODE-STYLE.md
// says such a kernel cannot be timed in a loop without a reset between runs,
// and that the reset would then sit inside the measurement. A histogram has
// nowhere else to put its counts, so this program breaks that rule's premise
// and has to say how. Correctness is checked first, once per kernel per
// input, from a zeroed bin array. The timed loop then never zeroes and its
// counts are never read by anything. The fourth static_assert below is what
// makes that safe: it proves the worst case cannot overflow a counter.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o histogram histogram.cu
// Run: ./histogram
//
// 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)
// 64 MiB plus 611 bytes. The 611 is day 5's number and it is here for the
// same reason: the buffer is not a multiple of the grid, so the last
// iteration of every grid-stride loop below has threads that stop and threads
// that do not, and a tail bug shows up as a wrong count rather than as
// nothing.
constexpr size_t kBytes = 64ull * 1024ull * 1024ull + 611ull;
// One bin per byte value, so a byte is its own bin index and no kernel here
// needs a mapping step. 256 bins of 4 bytes is 1 KiB, which is the number the
// whole lesson turns on: a private copy that small fits anywhere.
constexpr int kBins = 256;
constexpr int kThreadsPerBlock = 256;
// A fixed grid, walked with a grid-stride loop, rather than one block per
// 256 bytes. The block count is the knob that prices privatization: the merge
// costs kBlocks * kBins global atomics, so 1024 blocks buy a 256-fold cut in
// global atomic traffic while leaving every SM more blocks than it can hold
// at once. Doubling it halves the work per block and doubles the merge.
constexpr int kBlocks = 1024;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr int kNumInputs = 3;
// 'e' rather than 0 or 255, so a buffer full of it is a plausible file and
// not obviously a test pattern.
constexpr unsigned char kOneByteValue = 'e';
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert(kBins == 256,
"one bin per byte value, so a byte indexes its own bin and the "
"kernels need no mapping step");
static_assert(static_cast<size_t>(kBins) * sizeof(unsigned int) <= 48u * 1024u,
"the private copy must fit the 48 KiB of shared memory a Turing "
"block gets without an opt-in call to cudaFuncSetAttribute");
static_assert(static_cast<unsigned long long>(kWarmupRuns + kTimedRuns) *
kBytes <=
0xFFFFFFFFull,
"the timed loop never zeroes the bins, so all 13 runs of the "
"worst-case input, every byte the same value, have to fit in "
"one 32-bit counter");
// A paragraph of this project's own prose, repeated to fill the buffer. The
// byte frequencies are English prose's: the space is the most common byte by
// a wide margin, a couple of dozen values carry nearly all of the mass, and
// most of the 256 bins are never touched at all. Using our own words keeps
// the input clear of anyone else's licence, and the program counts the real
// distribution at run time rather than this comment guessing at it.
//
// Repetition does not change a histogram, which only ever sees the multiset
// of byte values. It also cannot help the reads: 64 MiB is far past this
// card's L2 either way.
static const char kSampleText[] =
"A histogram counts how many times each value appears. On a CPU you walk "
"the array and add one to a counter. On a GPU every thread wants to add "
"one to the same small set of counters at the same moment, so the "
"counters become the bottleneck rather than the data. Privatization "
"gives each block its own copy of the counters in shared memory, lets a "
"block contend only with itself, and merges the copies into global "
"memory once per block instead of once per element. The bin count "
"decides whether that copy fits on the chip at all, and when it does "
"not you either split the bins across several passes or keep the copies "
"in global memory and settle for a smaller cut in contention.\n";
static const char* const kInputNames[kNumInputs] = {
"uniform bytes", "english text", "all one byte"};
// Fills `out` with input number `which`. Every byte is decided by arithmetic
// with fixed constants, so two people comparing numbers are comparing the
// same bytes on the same buffer.
static void fillInput(int which, unsigned char* out, size_t n) {
if (which == 0) {
unsigned int state = 12345u;
for (size_t i = 0; i < n; ++i) {
state = state * 1664525u + 1013904223u;
// The high byte, never the low one. The low bits of a linear
// congruential generator cycle with a very short period, which
// would skew the histogram by accident and quietly turn the
// flattest of the three inputs into a lumpy one.
out[i] = static_cast<unsigned char>(state >> 24);
}
} else if (which == 1) {
const size_t len = sizeof(kSampleText) - 1;
size_t p = 0;
for (size_t i = 0; i < n; ++i) {
out[i] = static_cast<unsigned char>(kSampleText[p]);
++p;
if (p == len) {
p = 0;
}
}
} else {
for (size_t i = 0; i < n; ++i) {
out[i] = kOneByteValue;
}
}
}
// The reference. Host only, plain loop, no allocation: the caller owns both
// buffers. Written to be obviously right rather than fast.
static void histogramCpu(const unsigned char* in, unsigned int* bins,
size_t n) {
for (int binIdx = 0; binIdx < kBins; ++binIdx) {
bins[binIdx] = 0u;
}
for (size_t i = 0; i < n; ++i) {
bins[in[i]] += 1u;
}
}
// The reference for sumBytes. `unsigned long long` and not `size_t` because
// this is a checksum rather than a count or an index: 64 MiB of bytes summed
// reaches about 1.7e10, which does not fit in 32 bits on any machine.
static unsigned long long sumBytesCpu(const unsigned char* in, size_t n) {
unsigned long long total = 0ull;
for (size_t i = 0; i < n; ++i) {
total += in[i];
}
return total;
}
// The share of the input that lands in the busiest bin, and how many bins are
// touched at all. These two numbers are why the rows in the table differ from
// each other, so the program derives them from the reference histogram rather
// than leaving a reader to assume them.
static void inputSkew(const unsigned int* bins, size_t n, double* maxShare,
int* used) {
unsigned int biggest = 0u;
int touched = 0;
for (int binIdx = 0; binIdx < kBins; ++binIdx) {
if (bins[binIdx] > biggest) {
biggest = bins[binIdx];
}
if (bins[binIdx] != 0u) {
++touched;
}
}
*maxShare = 100.0 * static_cast<double>(biggest) / static_cast<double>(n);
*used = touched;
}
// One thread reads many bytes and writes one register total. No atomics, no
// shared memory, no data-dependent addresses.
//
// Memory: consecutive lanes take consecutive bytes, so one warp's 32
// addresses cover 32 contiguous bytes, which is one 32-byte sector. The one
// store per thread adds 1 MiB of writes to 64 MiB of reads.
//
// This is the floor row. Whatever the two histogram kernels cost, they cannot
// beat the time it takes to look at the bytes, and this measures that in the
// same program, on the same buffer, from the same grid. Borrowing a copy
// bandwidth from another day would fold that day's block shape, buffer size
// and clock state into the answer.
//
// Launch assumption: any grid. The loop condition is the bounds check.
// snippet: sum-kernel
__global__ void sumBytes(const unsigned char* __restrict__ in,
unsigned int* __restrict__ partials, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
const size_t t = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
unsigned int sum = 0u;
for (size_t i = t; i < n; i += step) {
sum += in[i];
}
partials[t] = sum;
}
// end snippet: sum-kernel
// One thread reads many bytes and adds one to a global counter for each.
//
// Memory: the reads are the same coalesced stream as sumBytes. The writes are
// not really writes. Every one is an atomic read-modify-write on a 1 KiB
// array the whole grid shares, so what this kernel costs is set by how many
// threads want the same counter at the same moment, which is a property of
// the input and not of this code.
//
// Launch assumption: any grid. The loop condition is the bounds check.
// snippet: global-kernel
__global__ void histogramGlobal(const unsigned char* __restrict__ in,
unsigned int* __restrict__ bins, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
atomicAdd(&bins[in[i]], 1u);
}
}
// end snippet: global-kernel
// One thread reads many bytes and adds one to its own block's copy of the
// counters, which lives in shared memory. The block then merges that copy
// into the global one with exactly one atomic per thread, because this block
// happens to have as many threads as there are bins.
//
// Memory: the reads are unchanged. The atomics move from a 1 KiB array that
// the whole grid shares to a 1 KiB array that 256 threads share, so a thread
// waits behind at most 255 rivals instead of 262,143, and the atomics that
// reach global memory fall from n to gridDim.x * kBins.
//
// The private copy is 256 words across 32 banks, so bins 32 apart share a
// bank. Day 15 is where that costs something; here it sits underneath the
// atomic conflicts and this program does not separate the two.
//
// Launch assumption: any grid and any block size. Both loops over the bins
// are strided by blockDim.x rather than guarded by `threadIdx.x < kBins`. The
// guarded form is correct only while the block has at least kBins threads,
// and it drops bins in silence the day somebody tries 128.
// snippet: private-kernel
__global__ void histogramPrivate(const unsigned char* __restrict__ in,
unsigned int* __restrict__ bins, size_t n) {
__shared__ unsigned int privateBins[kBins];
const int binStep = static_cast<int>(blockDim.x);
for (int binIdx = static_cast<int>(threadIdx.x); binIdx < kBins;
binIdx += binStep) {
privateBins[binIdx] = 0u;
}
__syncthreads();
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
atomicAdd(&privateBins[in[i]], 1u);
}
__syncthreads();
for (int binIdx = static_cast<int>(threadIdx.x); binIdx < kBins;
binIdx += binStep) {
atomicAdd(&bins[binIdx], privateBins[binIdx]);
}
}
// end snippet: private-kernel
// Times a launch with CUDA events and returns the mean milliseconds per run.
// A host-side clock around a launch measures the launch, not the kernel,
// because launches are asynchronous. Day 9 takes that apart.
//
// This is the one template and the one lambda allowed in module 1 to 3 code.
// Copy it 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;
}
// Every kernel here reads the same kBytes. The bin array is 1 KiB against 64
// MiB, so the figure below counts the input stream only and says so on the
// page. It is not the whole traffic and it is not meant to be; it is what
// makes the histogram rows comparable to the floor row. Host to device copies
// are outside every timed region.
static double bandwidthGBs(float ms) {
const double bytes = static_cast<double>(kBytes);
return bytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// Compares one finished histogram against the reference. Returns
// EXIT_SUCCESS only when all 256 bins agree and they sum to the input length.
static int checkHistogram(const char* name, const char* inputName,
const unsigned int* got, const unsigned int* want,
size_t n) {
for (int binIdx = 0; binIdx < kBins; ++binIdx) {
if (got[binIdx] != want[binIdx]) {
std::fprintf(stderr, "%s wrong on %s: bin %d holds %u, want %u\n",
name, inputName, binIdx, got[binIdx], want[binIdx]);
return EXIT_FAILURE;
}
}
// The cheap invariant, and the one that does not trust the reference:
// every byte lands in exactly one bin, so the bins sum to the number of
// bytes. A merge written `bins[b] = privateBins[b]` instead of an
// atomicAdd fails here, because then the last block to finish wins and
// the total comes out a small fraction of n.
unsigned long long total = 0ull;
for (int binIdx = 0; binIdx < kBins; ++binIdx) {
total += got[binIdx];
}
if (total != static_cast<unsigned long long>(n)) {
std::fprintf(stderr, "%s wrong on %s: bins sum to %llu, want %zu\n",
name, inputName, total, n);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
// Runs both histogram kernels once over whatever is in d_in and checks both.
// Returns EXIT_SUCCESS only when both are right.
//
// Nothing here is timed. The copies back would sit inside the measurement,
// and the first launch of each kernel still pays its own module load.
//
// The bin array is zeroed before each launch, which is the reason correctness
// lives here and not in the timed loop: this is the only place the counts are
// ever read.
static int checkHistograms(const char* inputName, const unsigned char* d_in,
unsigned int* d_bins, unsigned int* h_bins,
const unsigned int* h_want, size_t binBytes) {
int status = EXIT_SUCCESS;
CUDA_CHECK(cudaMemset(d_bins, 0, binBytes));
histogramGlobal<<<kBlocks, kThreadsPerBlock>>>(d_in, d_bins, kBytes);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_bins, d_bins, binBytes, cudaMemcpyDeviceToHost));
if (checkHistogram("histogramGlobal", inputName, h_bins, h_want, kBytes) !=
EXIT_SUCCESS) {
status = EXIT_FAILURE;
}
CUDA_CHECK(cudaMemset(d_bins, 0, binBytes));
histogramPrivate<<<kBlocks, kThreadsPerBlock>>>(d_in, d_bins, kBytes);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_bins, d_bins, binBytes, cudaMemcpyDeviceToHost));
if (checkHistogram("histogramPrivate", inputName, h_bins, h_want, kBytes) !=
EXIT_SUCCESS) {
status = EXIT_FAILURE;
}
return status;
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
const size_t binBytes = static_cast<size_t>(kBins) * sizeof(unsigned int);
const size_t threads =
static_cast<size_t>(kBlocks) * static_cast<size_t>(kThreadsPerBlock);
const size_t partialBytes = threads * sizeof(unsigned int);
std::vector<unsigned char> h_in(kBytes);
std::vector<unsigned int> h_bins(kBins);
std::vector<unsigned int> h_want(kBins);
std::vector<unsigned int> h_partials(threads);
unsigned char* d_in = nullptr;
unsigned int* d_bins = nullptr;
unsigned int* d_partials = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, kBytes));
CUDA_CHECK(cudaMalloc(&d_bins, binBytes));
CUDA_CHECK(cudaMalloc(&d_partials, partialBytes));
int status = EXIT_SUCCESS;
float floorMs = 0.0f;
float globalMs[kNumInputs] = {0.0f, 0.0f, 0.0f};
float privateMs[kNumInputs] = {0.0f, 0.0f, 0.0f};
double maxShare[kNumInputs] = {0.0, 0.0, 0.0};
int binsUsed[kNumInputs] = {0, 0, 0};
for (int which = 0; which < kNumInputs; ++which) {
fillInput(which, h_in.data(), kBytes);
histogramCpu(h_in.data(), h_want.data(), kBytes);
inputSkew(h_want.data(), kBytes, &maxShare[which], &binsUsed[which]);
CUDA_CHECK(
cudaMemcpy(d_in, h_in.data(), kBytes, cudaMemcpyHostToDevice));
if (checkHistograms(kInputNames[which], d_in, d_bins, h_bins.data(),
h_want.data(), binBytes) != EXIT_SUCCESS) {
status = EXIT_FAILURE;
break;
}
// The floor is measured once, on the first input, because sumBytes
// does the same work whatever the bytes say. It runs from the same
// grid over the same buffer as the two rows it is a floor for.
if (which == 0) {
sumBytes<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials, kBytes);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_partials.data(), d_partials, partialBytes,
cudaMemcpyDeviceToHost));
unsigned long long got = 0ull;
for (size_t t = 0; t < threads; ++t) {
got += h_partials[t];
}
const unsigned long long want = sumBytesCpu(h_in.data(), kBytes);
if (got != want) {
std::fprintf(stderr,
"sumBytes wrong: partials sum to %llu, want "
"%llu\n",
got, want);
status = EXIT_FAILURE;
break;
}
floorMs = timeKernel([&] {
sumBytes<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials,
kBytes);
});
}
// Zeroed before each timing, never inside it. Thirteen runs then pile
// up in the bins and nothing reads them, which the fourth
// static_assert at the top of this file is what makes safe to do.
CUDA_CHECK(cudaMemset(d_bins, 0, binBytes));
globalMs[which] = timeKernel([&] {
histogramGlobal<<<kBlocks, kThreadsPerBlock>>>(d_in, d_bins,
kBytes);
});
CUDA_CHECK(cudaMemset(d_bins, 0, binBytes));
privateMs[which] = timeKernel([&] {
histogramPrivate<<<kBlocks, kThreadsPerBlock>>>(d_in, d_bins,
kBytes);
});
}
// A real branch returning EXIT_FAILURE, not an assert(), because CI
// builds Release and NDEBUG deletes assert(). This is arithmetic and not
// a claim about the card: a zero or negative elapsed time means the event
// pair never separated, and every figure derived from it would be
// infinite or nonsense.
if (status == EXIT_SUCCESS) {
for (int which = 0; which < kNumInputs; ++which) {
if (globalMs[which] <= 0.0f || privateMs[which] <= 0.0f) {
std::fprintf(stderr,
"timing broken on %s: %.6f and %.6f ms, and "
"neither can be zero\n",
kInputNames[which],
static_cast<double>(globalMs[which]),
static_cast<double>(privateMs[which]));
status = EXIT_FAILURE;
}
}
if (floorMs <= 0.0f) {
std::fprintf(stderr, "timing broken on sumBytes: %.6f ms\n",
static_cast<double>(floorMs));
status = EXIT_FAILURE;
}
}
if (status == EXIT_SUCCESS) {
const double mib = static_cast<double>(kBytes) / (1024.0 * 1024.0);
const double mergeCut =
static_cast<double>(kBytes) /
(static_cast<double>(kBlocks) * static_cast<double>(kBins));
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name,
prop.major, prop.minor);
std::printf(
"Shared memory: %zu bytes per block by default, %zu with the "
"opt-in,\n %zu per SM\n",
prop.sharedMemPerBlock, prop.sharedMemPerBlockOptin,
prop.sharedMemPerMultiprocessor);
std::printf("Input: %zu bytes (%.2f MiB), %d bins of %zu bytes\n",
kBytes, mib, kBins, sizeof(unsigned int));
std::printf(
"Grid: %d blocks x %d threads. Mean of %d runs after %d "
"warm-ups.\n\n",
kBlocks, kThreadsPerBlock, kTimedRuns, kWarmupRuns);
std::printf("Atomic updates that reach global memory, per run\n");
std::printf(" histogramGlobal %12zu\n", kBytes);
std::printf(" histogramPrivate %12zu (a factor of %.1f)\n\n",
static_cast<size_t>(kBlocks) * static_cast<size_t>(kBins),
mergeCut);
std::printf("Reading every byte and doing nothing else with it\n");
std::printf(" sumBytes %10.3f ms %8.1f GB/s\n\n",
static_cast<double>(floorMs), bandwidthGBs(floorMs));
std::printf("%-15s %8s %10s %12s %13s %8s\n", "input", "max bin",
"bins used", "global (ms)", "private (ms)", "speedup");
std::printf("%-15s %8s %10s %12s %13s %8s\n", "--------------",
"-------", "---------", "-----------", "------------",
"-------");
for (int which = 0; which < kNumInputs; ++which) {
std::printf(
"%-15s %7.2f%% %10d %12.3f %13.3f %8.1f\n", kInputNames[which],
maxShare[which], binsUsed[which],
static_cast<double>(globalMs[which]),
static_cast<double>(privateMs[which]),
static_cast<double>(globalMs[which] / privateMs[which]));
}
std::printf(
"\nEvery row read the same %zu bytes from the same grid and the "
"same\nblock, through the same two kernels. Only the contents of "
"the buffer\nchanged between rows, so the spread down each "
"column is contention.\n",
kBytes);
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_bins));
CUDA_CHECK(cudaFree(d_partials));
return status;
}