code/day08-grid-stride/grid_stride.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 8: bounds checks and grid-stride loops.
//
// Adds two vectors three ways at ten sizes, from 0 up to 2^22, and prints
// which of the three is right at each size:
//
// A one thread per element, grid sized to n (day 5's program)
// B one thread per element, grid sized to the GPU
// C a grid-stride loop, grid sized to the GPU
//
// B is the column to watch. It is right while the grid it was given covers n
// and silently wrong after, and no CUDA call reports anything either way. C is
// right at every size and under every launch configuration, which is the whole
// argument for the loop.
//
// Nothing here is timed. Day 9 is where timing arrives, for the reason day 5
// gives: a clock around a first CUDA program measures context creation and two
// PCIe copies rather than the kernel.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o grid_stride grid_stride.cu
// Run: ./grid_stride
//
// 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 <cmath>
#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)
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarpSize = 32; // 32 on every GPU this course targets
// 611 is 13 x 47, so no block size divides it. Day 5 used it to make one
// bounds check run; here it is one row of ten.
constexpr size_t kOddElems = 611;
constexpr size_t kMaxElems = 4ull * 1024ull * 1024ull; // 2^22
constexpr float kRelTolerance = 1e-5f;
// Table 30 of the programming guide's compute-capability appendix caps a
// grid's x dimension at 2^31 - 1 blocks and its y and z dimensions at 65,535,
// on every architecture this course targets. README.md carries the URL.
constexpr size_t kMaxGridX = 2147483647ull;
// Returned by runOne when every element it compared was right. Distinct from
// any index, so "right at n = 0" and "wrong at index 0" cannot collide.
constexpr size_t kNoMismatch = static_cast<size_t>(-1);
// Blocks a grid needs to give one thread to each of n elements. constexpr so
// the static_asserts below can call it, and used at run time by the launches.
constexpr size_t blocksFor(size_t n) {
return (n + kThreadsPerBlock - 1) / kThreadsPerBlock;
}
static_assert(kThreadsPerBlock % kWarpSize == 0,
"block size must be a whole number of warps");
static_assert(kOddElems % kThreadsPerBlock != 0,
"the odd size must not divide by the block size, or the ragged "
"last block never happens");
static_assert(blocksFor(kOddElems) == 3,
"611 at 256 threads per block is 3 blocks");
static_assert(blocksFor(0) == 0,
"a grid sized to n is empty at n = 0, which is not a launch");
static_assert(blocksFor(kMaxElems) < kMaxGridX,
"a grid sized to kMaxElems must fit the grid x-dimension cap");
// out[i] = a[i] + b[i]. One thread owns one element. This is day 5's kernel,
// unchanged.
//
// Memory: consecutive threads take consecutive elements, so one warp's 32
// addresses cover 128 contiguous bytes.
//
// Launch assumption, and it is the subject of this lesson:
// gridDim.x * blockDim.x >= n. Give it a smaller grid and it writes the first
// gridDim.x * blockDim.x elements, leaves the rest untouched, and returns
// without reporting anything.
// snippet: one-per-element
__global__ void addOnePerElement(const float* a, const float* b, float* out,
size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
out[i] = a[i] + b[i];
}
}
// end snippet
// out[i] = a[i] + b[i] for every i, whatever grid you launch.
//
// Memory: the step is the whole grid, so inside one iteration consecutive
// threads still hold consecutive addresses and one warp still covers 128
// contiguous bytes. The pattern the loop repeats is the coalesced one, which
// is what day 11 measures.
//
// Launch assumption: none. Any grid with at least one thread is correct,
// including <<<1, 1>>>, which is how you serialise the kernel to debug it.
//
// The loop condition is the bounds check. There is no second guard, and the
// cast on blockDim.x matters for the same reason it does in the index: both
// blockDim.x and gridDim.x are unsigned int, so the product wraps at 2^32
// threads without it.
// snippet: grid-stride
__global__ void addGridStride(const float* a, const float* b, float* out,
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) {
out[i] = a[i] + b[i];
}
}
// end snippet
// CPU reference. Written for obvious correctness, not speed: plain loop, no
// OpenMP, no intrinsics. It never allocates; the caller owns every buffer.
static void vectorAddCpu(const float* a, const float* b, float* out, size_t n) {
for (size_t i = 0; i < n; ++i) {
out[i] = a[i] + b[i];
}
}
// Returns the first index where got and want differ by more than the relative
// tolerance, or n if they agree everywhere.
//
// The comparison is written as !(diff <= tol) and not as (diff > tol), because
// runOne fills the output with a NaN sentinel first and every comparison
// against a NaN is false. Under (diff > tol) an element nobody wrote would
// pass; under this form it fails, which is the whole point of the sentinel.
static size_t firstMismatch(const float* got, const float* want, size_t n,
float relTolerance) {
for (size_t i = 0; i < n; ++i) {
const float scale = (want[i] == 0.0f) ? 1.0f : std::fabs(want[i]);
if (!(std::fabs(got[i] - want[i]) <= relTolerance * scale)) {
return i;
}
}
return n;
}
// Runs one kernel under one launch configuration over the first n elements and
// returns the index of the first wrong element, or kNoMismatch when all n are
// right.
//
// The output is filled with a NaN sentinel first, so an element no thread
// wrote fails the comparison instead of matching a zero by luck. At n = 0
// there is nothing to compare, so the check becomes: did the sentinel at index
// 0 survive? A kernel that writes out[0] when n is 0 has an off-by-one.
//
// A launch the driver refuses is reported through *launchErr and the
// comparison is skipped, because an empty grid is one of the rows in the table
// rather than a bug in this program. The caller reads *launchErr first.
static size_t runOne(bool useGridStride, int blocks, int threads,
const float* d_a, const float* d_b, float* d_out,
float* h_out, const float* h_want, size_t n,
cudaError_t* launchErr) {
const size_t sentinelElems = (n == 0) ? 1 : n;
const size_t sentinelBytes = sentinelElems * sizeof(float);
CUDA_CHECK(cudaMemset(d_out, 0xFF, sentinelBytes));
if (useGridStride) {
addGridStride<<<blocks, threads>>>(d_a, d_b, d_out, n);
} else {
addOnePerElement<<<blocks, threads>>>(d_a, d_b, d_out, n);
}
// Read rather than checked, for the reason above. Reading it also clears
// the error state, so a refused launch here cannot poison the next row.
*launchErr = cudaGetLastError();
if (*launchErr != cudaSuccess) {
return kNoMismatch;
}
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out, d_out, sentinelBytes, cudaMemcpyDeviceToHost));
if (n == 0) {
return std::isnan(h_out[0]) ? kNoMismatch : 0;
}
const size_t bad = firstMismatch(h_out, h_want, n, kRelTolerance);
return (bad == n) ? kNoMismatch : bad;
}
// One cell of the table: what happened, in at most thirteen characters.
static void renderCell(char* buf, size_t bufLen, size_t bad,
cudaError_t launchErr) {
if (launchErr != cudaSuccess) {
std::snprintf(buf, bufLen, "refused");
} else if (bad == kNoMismatch) {
std::snprintf(buf, bufLen, "ok");
} else {
std::snprintf(buf, bufLen, "wrong@%zu", bad);
}
}
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);
// The grid for columns B and C is sized to the GPU, not to the data. Ask
// the driver how many blocks of this kernel fit on one SM, multiply by the
// SM count, and that grid covers the machine once. n appears nowhere in
// it, which is the point: the same three numbers serve every row below.
// snippet: device-grid
int blocksPerSm = 0;
CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
&blocksPerSm, addGridStride, kThreadsPerBlock, 0));
const int deviceBlocks = blocksPerSm * prop.multiProcessorCount;
const size_t totalThreads =
static_cast<size_t>(deviceBlocks) * kThreadsPerBlock;
// end snippet
std::printf(
"device grid: %d blocks x %d threads = %zu threads "
"(%d blocks per SM x %d SMs)\n",
deviceBlocks, kThreadsPerBlock, totalThreads, blocksPerSm,
prop.multiProcessorCount);
// The two sizes that matter most are read off the hardware rather than
// written down: exactly one thread each, and one more element than that.
// The list is in run order, not sorted, because where the device grid
// lands among the constants depends on the card.
const size_t sizes[] = {0,
1,
31,
32,
kOddElems,
1024,
1025,
totalThreads,
totalThreads + 1,
kMaxElems};
constexpr int kNumSizes =
static_cast<int>(sizeof(sizes) / sizeof(sizes[0]));
size_t maxElems = 0;
for (int c = 0; c < kNumSizes; ++c) {
if (sizes[c] > maxElems) {
maxElems = sizes[c];
}
}
const size_t bytes = maxElems * sizeof(float);
// Small whole numbers, so every sum is exact on both sides and a mismatch
// below can only be an index bug. h_want is computed once at the largest
// size; every smaller case compares against its prefix.
std::vector<float> h_a(maxElems);
std::vector<float> h_b(maxElems);
std::vector<float> h_out(maxElems);
std::vector<float> h_want(maxElems);
for (size_t i = 0; i < maxElems; ++i) {
h_a[i] = static_cast<float>(i);
h_b[i] = static_cast<float>(2 * i);
}
vectorAddCpu(h_a.data(), h_b.data(), h_want.data(), maxElems);
float* d_a = nullptr;
float* d_b = nullptr;
float* d_out = nullptr;
CUDA_CHECK(cudaMalloc(&d_a, bytes));
CUDA_CHECK(cudaMalloc(&d_b, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), bytes, cudaMemcpyHostToDevice));
std::printf("\nA = one thread per element, grid sized to n (day 5)\n");
std::printf("B = one thread per element, grid sized to the GPU\n");
std::printf("C = grid-stride loop, grid sized to the GPU\n\n");
std::printf("%12s %11s %13s %13s %13s\n", "n", "blocks(A)", "A", "B", "C");
int failures = 0;
cudaError_t zeroGridErr = cudaSuccess;
for (int c = 0; c < kNumSizes; ++c) {
const size_t n = sizes[c];
const int dataBlocks = static_cast<int>(blocksFor(n));
cudaError_t errA = cudaSuccess;
cudaError_t errB = cudaSuccess;
cudaError_t errC = cudaSuccess;
const size_t badA =
runOne(false, dataBlocks, kThreadsPerBlock, d_a, d_b, d_out,
h_out.data(), h_want.data(), n, &errA);
const size_t badB =
runOne(false, deviceBlocks, kThreadsPerBlock, d_a, d_b, d_out,
h_out.data(), h_want.data(), n, &errB);
const size_t badC =
runOne(true, deviceBlocks, kThreadsPerBlock, d_a, d_b, d_out,
h_out.data(), h_want.data(), n, &errC);
if (n == 0) {
zeroGridErr = errA;
}
char cellA[24];
char cellB[24];
char cellC[24];
renderCell(cellA, sizeof(cellA), badA, errA);
renderCell(cellB, sizeof(cellB), badB, errB);
renderCell(cellC, sizeof(cellC), badC, errC);
std::printf("%12zu %11d %13s %13s %13s\n", n, dataBlocks, cellA, cellB,
cellC);
// Gate 1. The grid-stride kernel is right at every size and its launch
// is never refused, because its grid never depends on n.
if (errC != cudaSuccess || badC != kNoMismatch) {
std::fprintf(stderr,
"C wrong at n = %zu: first bad index %zu, launch %s\n",
n, badC, cudaGetErrorString(errC));
++failures;
}
// Gate 2. One thread per element is right whenever the grid was sized
// to n and the driver accepted the launch. Day 5's program, restated
// as a check.
if (errA == cudaSuccess && badA != kNoMismatch) {
std::fprintf(stderr, "A wrong at n = %zu: first bad index %zu\n", n,
badA);
++failures;
}
// Gate 3, the prediction, written before the first run. One thread per
// element on a grid sized to the GPU is right exactly while
// n <= totalThreads, and when it is wrong the first wrong index is
// totalThreads itself, because that is the first element no thread
// owns. A run that disagrees in either direction means the model on
// the page is wrong, not that the machine is unlucky.
const size_t wantB = (n <= totalThreads) ? kNoMismatch : totalThreads;
if (badB != wantB) {
char wantCell[24];
renderCell(wantCell, sizeof(wantCell), wantB, cudaSuccess);
std::fprintf(stderr, "B at n = %zu: predicted %s, measured %s\n", n,
wantCell, cellB);
++failures;
}
}
if (zeroGridErr != cudaSuccess) {
std::printf(
"\nA at n = 0 asks for a grid of 0 blocks. This driver refused "
"it: %s\n",
cudaGetErrorString(zeroGridErr));
} else {
std::printf(
"\nA at n = 0 asks for a grid of 0 blocks. This driver accepted "
"it and nothing was written.\n");
}
// The same grid-stride kernel, one size, four launch configurations. This
// is the property the rest of the course leans on: day 10 sweeps block
// sizes without touching a kernel, because a grid-stride kernel does not
// care what grid it was given. <<<1, 1>>> is the debugging case, one
// thread walking all 611 elements in order.
const int probeBlocks[] = {1, 1, deviceBlocks, 10000};
const int probeThreads[] = {1, 256, kThreadsPerBlock, 32};
constexpr int kNumProbes =
static_cast<int>(sizeof(probeBlocks) / sizeof(probeBlocks[0]));
std::printf("\nC alone, n = %zu, four launch configurations:\n", kOddElems);
for (int p = 0; p < kNumProbes; ++p) {
cudaError_t err = cudaSuccess;
const size_t bad =
runOne(true, probeBlocks[p], probeThreads[p], d_a, d_b, d_out,
h_out.data(), h_want.data(), kOddElems, &err);
char cellText[24];
renderCell(cellText, sizeof(cellText), bad, err);
std::printf(" <<<%5d, %4d>>> %s\n", probeBlocks[p], probeThreads[p],
cellText);
if (err != cudaSuccess || bad != kNoMismatch) {
++failures;
}
}
// Free before reporting, so the failure path frees too.
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_out));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
std::printf(
"\nall %d sizes and %d launch configurations agree with the "
"CPU reference where the model says they should\n",
kNumSizes, kNumProbes);
return EXIT_SUCCESS;
}