code/day11-coalescing/coalescing.cuThis is the source used by the lesson and its recorded evidence. Compile commands and expected output live in the directory README.
// Day 11: memory coalescing, measured honestly.
//
// The point of this program is NOT that strided access is slow. It is that
// measuring it wrong is easy. Every variant here touches every element of the
// buffer exactly once, so the byte count is identical and the only thing that
// changes is the order the addresses arrive in.
//
// Getting that right took two attempts. The first version indexed with
// (i * stride) % n. With n a power of two that map is not a bijection: at
// stride 16 it touches half the 32-byte sectors and at stride 32 a quarter,
// so those rows reported two and four times the bandwidth the hardware
// actually moved. They were also the two rows the lesson called "the floor".
// The map below is a permutation for every stride, verified before use.
//
// Build: nvcc -O3 -arch=sm_75 -o coalescing coalescing.cu
// Run: ./coalescing human readable
// ./coalescing --json one JSON object per row, for CI
//
// 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 <cstring>
#include <cuda_runtime.h>
#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 Mi floats, 256 MiB per buffer. A power of two so that n / stride is
// exact for every stride we test.
static const size_t kElems = 64ull * 1024ull * 1024ull;
static const int kThreadsPerBlock = 256;
static const int kWarmupRuns = 3;
static const int kTimedRuns = 10;
// snippet: permuted-kernel
// One kernel serves every stride, so the instruction mix is identical across
// rows and only the address pattern differs.
//
// The launch is a 2D grid: x walks `m = n / stride` columns, y walks `stride`
// rows. Element index j = col * stride + row. Sweeping col over [0, m) and
// row over [0, stride) hits every j in [0, n) exactly once, so this is a
// permutation, not a sampling.
//
// Consecutive threads in a warp have consecutive `col`, so their addresses
// are `stride` elements apart. At stride 1 that is 4 bytes apart, which is
// perfectly coalesced. At stride 32 it is 128 bytes apart, so every lane
// lands in its own 32-byte sector.
//
// There is no modulo and no division here. An earlier version used `% n` on a
// runtime size_t, which nvcc cannot strength-reduce to a mask, so the strided
// rows paid for a 64-bit division the coalesced row did not. That broke the
// claim that the access pattern is the only variable.
__global__ void copyPermuted(const float* __restrict__ in,
float* __restrict__ out, size_t m, int stride) {
size_t col = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (col < m) {
size_t j = col * static_cast<size_t>(stride) + blockIdx.y;
out[j] = in[j];
}
}
// end snippet: permuted-kernel
// The CPU-correct, GPU-wrong pattern, kept as its own kernel because it is a
// different teaching point. Each thread takes a contiguous run of elements.
// Still a permutation: thread t owns [t * chunk, (t + 1) * chunk).
__global__ void copyChunked(const float* __restrict__ in,
float* __restrict__ out, size_t n, int chunk) {
size_t t = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
size_t base = t * static_cast<size_t>(chunk);
for (int k = 0; k < chunk; ++k) {
size_t idx = base + static_cast<size_t>(k);
if (idx < n) {
out[idx] = in[idx];
}
}
}
// Times one launch with CUDA events. A host-side clock around a launch
// measures the launch, not the kernel, because launches are asynchronous.
// See day 9.
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
// Warm up. The first launch of each kernel pays a module load cost, so
// every kernel that gets timed also gets warmed.
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 variant reads kElems floats and writes kElems floats. That equality is
// the whole point of the program, so the figure is computed once from one
// constant rather than per kernel, where the two could drift apart.
static double bandwidthGBs(float ms) {
const double bytes = 2.0 * static_cast<double>(kElems) * sizeof(float);
return bytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
int main(int argc, char** argv) {
bool json = (argc > 1 && std::strcmp(argv[1], "--json") == 0);
int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
float* d_in = nullptr;
float* d_out = nullptr;
const size_t bytes = kElems * sizeof(float);
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
CUDA_CHECK(cudaMemset(d_in, 1, bytes));
const int strides[] = {1, 2, 4, 8, 16, 32};
const int kNumStrides = sizeof(strides) / sizeof(strides[0]);
const int kChunk = 32;
const int kNumRows = kNumStrides + 1;
if (!json) {
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name,
prop.major, prop.minor);
std::printf("Elements: %zu floats, %.0f MiB per buffer\n\n", kElems,
static_cast<double>(bytes) / (1024.0 * 1024.0));
std::printf("%-26s %11s %11s\n", "pattern", "time (ms)", "GB/s");
std::printf("%-26s %11s %11s\n", "--------------------------",
"----------", "----------");
}
int rows = 0;
for (int s = 0; s < kNumStrides; ++s) {
const int stride = strides[s];
const size_t m = kElems / static_cast<size_t>(stride);
const int gridX =
static_cast<int>((m + kThreadsPerBlock - 1) / kThreadsPerBlock);
dim3 grid(gridX, stride);
float ms = timeKernel([&] {
copyPermuted<<<grid, kThreadsPerBlock>>>(d_in, d_out, m, stride);
});
++rows;
char label[64];
std::snprintf(label, sizeof(label), "permuted (stride %d)", stride);
if (json) {
std::printf(
"{\"pattern\":\"stride\",\"stride\":%d,\"ms\":%.6f,"
"\"gbps\":%.3f,\"bytes\":%zu}\n",
stride, ms, bandwidthGBs(ms), 2 * bytes);
} else {
std::printf("%-26s %11.3f %11.1f\n", label, ms, bandwidthGBs(ms));
}
}
const size_t kChunkSz = static_cast<size_t>(kChunk);
const size_t chunkThreads = (kElems + kChunkSz - 1) / kChunkSz;
const int chunkBlocks = static_cast<int>(
(chunkThreads + kThreadsPerBlock - 1) / kThreadsPerBlock);
float msChunk = timeKernel([&] {
copyChunked<<<chunkBlocks, kThreadsPerBlock>>>(d_in, d_out, kElems,
kChunk);
});
++rows;
if (json) {
std::printf(
"{\"pattern\":\"chunked\",\"chunk\":%d,\"ms\":%.6f,"
"\"gbps\":%.3f,\"bytes\":%zu}\n",
kChunk, msChunk, bandwidthGBs(msChunk), 2 * bytes);
} else {
std::printf("%-26s %11.3f %11.1f\n", "chunked (32 per thread)", msChunk,
bandwidthGBs(msChunk));
std::printf(
"\nEvery row above moved exactly %.0f MiB: each variant is a\n"
"permutation of the same buffer, so no row can win by doing less\n"
"work. If a coalescing benchmark shows strided access winning,\n"
"check that invariant first.\n",
2.0 * static_cast<double>(bytes) / (1024.0 * 1024.0));
}
// The row count is checked so the lesson's prose cannot drift from what
// the program actually prints. A real branch, not an assert: CI builds
// Release, Release defines NDEBUG, and NDEBUG deletes assert(), so the
// check would be absent from exactly the build that matters.
if (rows != kNumRows) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's table and "
"this program disagree\n",
rows, kNumRows);
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_out));
return EXIT_FAILURE;
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_out));
return 0;
}