code/day18-constant-memory/constant_memory.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 18: a convolution filter in constant memory against the same filter in
// global memory, swept across five problem sizes.
//
// The point is not that constant memory is fast. It is that the constant
// cache charges by the number of distinct addresses a warp asks for, so the
// same 32 floats are nearly free when every lane wants the same tap and
// expensive when the lanes want 32 different taps. One kernel pair covers
// both: rotMask is a runtime argument, 0 for the broadcast pattern and
// kTapMask for the per-lane pattern, so the two rows execute the same
// instructions and only the address pattern changes.
//
// Both kernels in a row read the same inputs, apply the same taps and write
// the same outputs, and both are checked against the same CPU reference
// before either time is printed. Neither can win by doing less work, which is
// the discipline day 11 exists to teach.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o constant_memory constant_memory.cu
// Run: ./constant_memory human readable
// ./constant_memory --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 <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <vector>
#include <cuda_runtime.h>
// The one error macro. This file is standalone, so it carries its own
// verbatim copy. `err_` has a trailing underscore so it cannot collide with a
// variable at the call site, and the do/while makes the macro one statement.
#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 taps is 128 bytes, which is one 128-byte cache line and a rounding error
// against the 64 KiB of constant memory. That is the honest case: a real
// convolution filter is tiny, and the question is whether where you put
// something this small can be measured at all.
constexpr int kFilterTaps = 32;
constexpr int kTapMask = kFilterTaps - 1;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kRelTolerance = 1e-5f;
// Table 31 of the compute-capability appendix gives 64 KB of constant memory
// for every compute capability from 7.5 to 12.x, so this bound is the same on
// every card the course targets.
constexpr size_t kConstantBytes = 64ull * 1024ull;
// The sweep. 4 Ki elements is 32 KiB of traffic and sits inside L2 on any
// card in the matrix; 8 Mi elements is 64 MiB and cannot. The middle sizes
// are there so the crossover has somewhere to happen.
constexpr int kNumSizes = 5;
constexpr size_t kSizes[kNumSizes] = {1ull << 12, 1ull << 15, 1ull << 18,
1ull << 21, 1ull << 23};
constexpr int kPatterns = 2;
constexpr int kExpectedRows = kNumSizes * kPatterns;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert((kFilterTaps & kTapMask) == 0,
"the rotation wraps with & kTapMask, so the tap count has to be "
"a power of two");
static_assert(kFilterTaps >= 32,
"the per-lane row claims a warp asks for 32 distinct taps, "
"which needs at least 32 of them to exist");
static_assert(kFilterTaps * sizeof(float) <= kConstantBytes,
"the filter has to fit in the 64 KiB of constant memory");
static_assert(kSizes[0] >= kThreadsPerBlock,
"the smallest sweep size must fill at least one block");
// The filter, in constant memory. CUDA-CODE-STYLE.md fixes h_, d_ and u_ for
// pointers and says nothing about a __constant__ symbol, so this file names
// it beside them: c_ is the constant-memory copy of the same 32 floats that
// d_filter holds in global memory.
__constant__ float c_filter[kFilterTaps];
// out[i] = sum over k of filter[(rot + k) & kTapMask] * in[i + k], with the
// filter read from constant memory.
//
// One thread owns one output element. At every k the 32 lanes of a warp read
// 32 consecutive `in` elements, which is 128 contiguous bytes and four
// sectors, so the input side is coalesced and identical in both kernels. The
// only thing that differs between them is where `filter` lives.
//
// rot is threadIdx.x & rotMask, and rotMask arrives at runtime so the
// compiler cannot fold it away. With rotMask 0 every lane reads the same tap
// at each k, which is the broadcast the constant cache is built for. With
// rotMask kTapMask the 32 lanes read 32 different taps, which the Best
// Practices Guide says the constant path serializes.
//
// Launch assumption: `in` holds n + kFilterTaps - 1 elements, so no lane
// reads past the end and the tap loop needs no bounds check of its own.
// snippet: constant-kernel
__global__ void conv1dConstant(const float* __restrict__ in,
float* __restrict__ out, size_t n, int rotMask) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
const int rot = static_cast<int>(threadIdx.x) & rotMask;
float acc = 0.0f;
for (int k = 0; k < kFilterTaps; ++k) {
acc += c_filter[(rot + k) & kTapMask] * in[i + k];
}
out[i] = acc;
}
}
// end snippet
// The same arithmetic with the filter in global memory. `filter` is a const
// pointer marked __restrict__, which the programming guide says compiles to a
// read-only cache load, the PTX ld.global.nc that __ldg() emits. On Turing
// that read-only path and L1 are the same 96 KiB of unified cache, so this is
// not a second cache next to L1, it is L1 told the data will not change.
//
// Warp addresses and the launch assumption are the same as above.
// snippet: global-kernel
__global__ void conv1dGlobal(const float* __restrict__ in,
const float* __restrict__ filter,
float* __restrict__ out, size_t n, int rotMask) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
const int rot = static_cast<int>(threadIdx.x) & rotMask;
float acc = 0.0f;
for (int k = 0; k < kFilterTaps; ++k) {
acc += filter[(rot + k) & kTapMask] * in[i + k];
}
out[i] = acc;
}
}
// end snippet
// CPU reference. Written for obvious correctness, not speed: plain loops, no
// OpenMP, no intrinsics. It never allocates; the caller owns every buffer.
//
// The rotation is a property of the lane, so the reference has to know the
// launch shape to reproduce it. Element i is owned by thread i % blockDim.x
// because every block is full and blockDim.x is kThreadsPerBlock on every
// launch in this program. Change the block size and this line changes too.
//
// Accumulates in double where the kernel accumulates in float: the
// reference's job is to be right, not to match bit for bit.
static void conv1dCpu(const float* in, const float* filter, float* out,
size_t n, int rotMask) {
for (size_t i = 0; i < n; ++i) {
const int rot = static_cast<int>(i % kThreadsPerBlock) & rotMask;
double acc = 0.0;
for (int k = 0; k < kFilterTaps; ++k) {
acc += static_cast<double>(filter[(rot + k) & kTapMask]) *
static_cast<double>(in[i + static_cast<size_t>(k)]);
}
out[i] = static_cast<float>(acc);
}
}
// Returns the first index where got and want differ by more than the relative
// tolerance, or n if they agree everywhere. Returning the index rather than a
// bool is the point: "wrong at 512" names the block, "wrong" does not.
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;
}
// 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.
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;
}
// The compulsory traffic for one launch: every input element read once, every
// output element written once. It counts elements, not load instructions,
// because each input element is read by up to kFilterTaps threads and those
// repeats are meant to be served by the caches. Both kernels in a row move
// this same figure, so it is computed here from n alone and never per kernel,
// where the two could drift apart.
static size_t movedBytes(size_t n) {
return (2 * n + kFilterTaps - 1) * sizeof(float);
}
static double bandwidthGBs(size_t n, float ms) {
const double bytes = static_cast<double>(movedBytes(n));
return bytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// argc and argv are here for --json and for nothing else. Every size in the
// sweep is a constant at the top of this file, so the page and the program
// cannot disagree about what was measured.
int main(int argc, char** argv) {
const bool json = (argc > 1 && std::strcmp(argv[1], "--json") == 0);
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
const size_t maxElems = kSizes[kNumSizes - 1];
const size_t inElems = maxElems + kFilterTaps - 1;
const size_t filterBytes = kFilterTaps * sizeof(float);
// Small whole numbers and quarters, so every product is a multiple of
// 0.25 and every 32-term sum is exact in float. A mismatch below can only
// be an index bug, never a rounding difference.
std::vector<float> h_filter(kFilterTaps);
std::vector<float> h_in(inElems);
std::vector<float> h_out(maxElems);
std::vector<float> h_want(maxElems);
std::vector<float> h_broadcast(maxElems);
for (int k = 0; k < kFilterTaps; ++k) {
h_filter[static_cast<size_t>(k)] = static_cast<float>(k % 5) - 2.0f;
}
for (size_t i = 0; i < inElems; ++i) {
h_in[i] = static_cast<float>(i % 17) * 0.25f;
}
float* d_in = nullptr;
float* d_out = nullptr;
float* d_filter = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, inElems * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_out, maxElems * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_filter, filterBytes));
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), inElems * sizeof(float),
cudaMemcpyHostToDevice));
// The same 32 floats land in both places. cudaMemcpyToSymbol takes the
// symbol itself, not its address: writing &c_filter compiles and then
// fails at runtime with `invalid device symbol`.
// snippet: fill-filter
CUDA_CHECK(cudaMemcpyToSymbol(c_filter, h_filter.data(), filterBytes));
CUDA_CHECK(cudaMemcpy(d_filter, h_filter.data(), filterBytes,
cudaMemcpyHostToDevice));
// end snippet
if (!json) {
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name,
prop.major, prop.minor);
std::printf(
"L2 %d KiB, persisting L2 max %d B, policy window max %d B\n",
prop.l2CacheSize / 1024, prop.persistingL2CacheMaxSize,
prop.accessPolicyMaxWindowSize);
if (prop.persistingL2CacheMaxSize == 0) {
std::printf(
"This device sets aside no L2 for persisting accesses, so "
"cudaAccessPolicyWindow has nothing to reserve.\n");
}
std::printf("Filter: %d taps, %zu bytes, of %zu bytes of constant\n",
kFilterTaps, filterBytes, kConstantBytes);
std::printf(
"%d threads per block, %d warm-ups, %d timed runs, mean "
"reported, no copies inside the measurement\n\n",
kThreadsPerBlock, kWarmupRuns, kTimedRuns);
std::printf("%9s %-9s %9s %10s %7s %12s\n", "n", "lanes", "const ms",
"global ms", "ratio", "global GB/s");
std::printf("%9s %-9s %9s %10s %7s %12s\n", "--------", "---------",
"--------", "---------", "-----", "-----------");
}
const int rotMasks[kPatterns] = {0, kTapMask};
const char* patternNames[kPatterns] = {"broadcast", "per-lane"};
int rows = 0;
size_t badIndex = 0;
size_t badElems = 0;
float badGot = 0.0f;
float badWant = 0.0f;
const char* badVariant = nullptr;
// The second experiment gate. rotOf stays 0 while the per-lane row really
// is a different access pattern from the broadcast row; a non-zero value
// means rotMask changed nothing and the two rows measured the same thing
// twice, which is the failure mode day 11 exists to warn about.
size_t rotSame = 0;
size_t rotOf = 0;
for (int s = 0; s < kNumSizes && badVariant == nullptr && rotOf == 0; ++s) {
const size_t n = kSizes[s];
const size_t outBytes = n * sizeof(float);
const int blocks =
static_cast<int>((n + kThreadsPerBlock - 1) / kThreadsPerBlock);
for (int p = 0; p < kPatterns && badVariant == nullptr; ++p) {
const int rotMask = rotMasks[p];
conv1dCpu(h_in.data(), h_filter.data(), h_want.data(), n, rotMask);
if (p == 0) {
std::memcpy(h_broadcast.data(), h_want.data(), outBytes);
}
const float msConst = timeKernel([&] {
conv1dConstant<<<blocks, kThreadsPerBlock>>>(d_in, d_out, n,
rotMask);
});
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, outBytes,
cudaMemcpyDeviceToHost));
const size_t badConst =
firstMismatch(h_out.data(), h_want.data(), n, kRelTolerance);
// Read the bad value out now. The next launch overwrites d_out
// and the copy after it overwrites h_out, so by the time the
// report is printed this element is gone.
const float gotConst = (badConst == n) ? 0.0f : h_out[badConst];
const float msGlobal = timeKernel([&] {
conv1dGlobal<<<blocks, kThreadsPerBlock>>>(d_in, d_filter,
d_out, n, rotMask);
});
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, outBytes,
cudaMemcpyDeviceToHost));
const size_t badGlobal =
firstMismatch(h_out.data(), h_want.data(), n, kRelTolerance);
// A real branch, not an assert. CI builds Release, Release
// defines NDEBUG, and NDEBUG deletes assert(), so the check would
// be missing from exactly the build that matters.
if (badConst != n) {
badVariant = "conv1dConstant";
badIndex = badConst;
badElems = n;
badGot = gotConst;
badWant = h_want[badConst];
} else if (badGlobal != n) {
badVariant = "conv1dGlobal";
badIndex = badGlobal;
badElems = n;
badGot = h_out[badGlobal];
badWant = h_want[badGlobal];
} else {
if (p == kPatterns - 1) {
size_t same = 0;
for (size_t e = 0; e < n; ++e) {
if (h_broadcast[e] == h_want[e]) {
++same;
}
}
if (same * 2 > n) {
rotSame = same;
rotOf = n;
}
}
++rows;
if (json) {
std::printf(
"{\"pattern\":\"%s\",\"n\":%zu,\"const_ms\":%.6f,"
"\"global_ms\":%.6f,\"ratio\":%.4f,"
"\"const_gbps\":%.3f,\"global_gbps\":%.3f,"
"\"bytes\":%zu}\n",
patternNames[p], n, static_cast<double>(msConst),
static_cast<double>(msGlobal),
static_cast<double>(msConst / msGlobal),
bandwidthGBs(n, msConst), bandwidthGBs(n, msGlobal),
movedBytes(n));
} else {
std::printf("%9zu %-9s %9.4f %10.4f %7.2f %12.1f\n", n,
patternNames[p], static_cast<double>(msConst),
static_cast<double>(msGlobal),
static_cast<double>(msConst / msGlobal),
bandwidthGBs(n, msGlobal));
}
}
}
}
if (badVariant == nullptr && rotOf == 0 && !json) {
std::printf(
"\nBoth kernels in a row read the same %d taps and the same\n"
"n + %d inputs, and write the same n outputs, so neither can win\n"
"by doing less work. Only where the filter lives changes.\n",
kFilterTaps, kFilterTaps - 1);
}
// One cleanup block. Every path below it, including both failure paths,
// has already freed every allocation.
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_out));
CUDA_CHECK(cudaFree(d_filter));
if (badVariant != nullptr) {
std::fprintf(stderr, "%s wrong at %zu of %zu: got %.9g, want %.9g\n",
badVariant, badIndex, badElems,
static_cast<double>(badGot), static_cast<double>(badWant));
return EXIT_FAILURE;
}
if (rotOf != 0) {
std::fprintf(stderr,
"the per-lane row is not a different access pattern: "
"%zu of %zu outputs at n = %zu match the broadcast row, "
"so rotMask changed nothing and both rows measured the "
"same access twice\n",
rotSame, rotOf, rotOf);
return EXIT_FAILURE;
}
if (rows != kExpectedRows) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's table and "
"this program disagree\n",
rows, kExpectedRows);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}