code/day13-shared-memory/shared_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 13: shared memory and tiling, measured on a transpose.
//
// Four kernels over one matrix, all launched with the same grid, the same
// block and the same four elements per thread, so the only thing that
// changes from row to row is the order the addresses arrive in:
//
// copyRowMajor the ceiling. No transpose, both ends coalesced.
// transposeNaive day 12's kernel. Reads rows, writes columns.
// transposeTiledStatic a 32x32 tile in __shared__, both ends coalesced.
// transposeTiledDynamic the same tile, sized at launch instead.
//
// Every row moves exactly 2 * kRows * kCols * sizeof(float) bytes, and that
// figure is computed once from the constants rather than per kernel, so no
// row can win by doing less work. Day 11 exists because a published
// benchmark got that wrong.
//
// The program then probes the per-block shared memory cap: it launches a
// kernel asking for more than the default, prints whatever the runtime says,
// opts in with cudaFuncSetAttribute, and launches the same kernel again.
// Day 74 needs that call and the failure it prevents has no obvious cause.
//
// What this program does not measure: bank conflicts. The tile below is
// [32][32] with no padding, so the second loop puts all 32 lanes of a warp
// on one bank. That is deliberate and day 15 is where it gets measured and
// fixed.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o shared_memory shared_memory.cu
// Run: ./shared_memory
//
// 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 <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)
// 2003 x 3001, both prime, so no tile size divides either side and the grid
// rounds up on both axes. Every guard in every kernel therefore runs on
// every launch. The matrix is also not square, so a kernel that swaps row
// and column produces the wrong shape rather than a plausible wrong answer.
constexpr size_t kRows = 2003;
constexpr size_t kCols = 3001;
// A 32 x 32 tile worked by 32 x 8 threads, so each thread moves four
// elements four rows apart. 32 columns is one warp's worth of consecutive
// floats, which is the entire point of the tile, and 32 * 8 is 256 threads,
// the course default block size.
constexpr int kTileDim = 32;
constexpr int kBlockRows = 8;
constexpr int kThreadsPerBlock = kTileDim * kBlockRows;
// The dynamic kernel asks for exactly the bytes the static one declares, so
// the two differ only in where the tile's size comes from.
constexpr size_t kTileBytes = sizeof(float) * kTileDim * kTileDim;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert(kThreadsPerBlock <= 1024,
"a thread block holds at most 1024 threads on every compute "
"capability this course targets");
static_assert(kTileDim % kBlockRows == 0,
"the load and store loops step by kBlockRows and have to land "
"on kTileDim exactly, or part of the tile is never touched");
static_assert(kRows % kTileDim != 0 && kCols % kTileDim != 0,
"both sides must divide unevenly by the tile, or the grid does "
"not round up and the bounds checks never run");
static_assert(kTileBytes <= 48 * 1024,
"a static __shared__ array has to fit the 48 KiB per-block "
"default; more than that needs the dynamic form plus "
"cudaFuncSetAttribute, which is what probeSharedMemoryCap() "
"below is about");
// Integer ceiling division. constexpr so one function can size the grid at
// run time and appear in a static_assert. A static_assert rather than an
// assert because CI builds Release, Release defines NDEBUG, and NDEBUG
// deletes assert(), so a check that can be made at compile time is made
// there.
static constexpr size_t ceilDiv(size_t a, size_t b) {
return (a + b - 1) / b;
}
static_assert(ceilDiv(kRows, static_cast<size_t>(kTileDim)) <= 65535,
"gridDim.y stops at 65535, unlike gridDim.x which goes to "
"2^31 - 1. A taller matrix than that needs its rows folded "
"into x");
// out[row][col] = in[row][col]. No transpose at all.
//
// This row is the ceiling. It moves the same bytes as the three transposes
// below with both ends coalesced, so it says what "as fast as this can go"
// means on the card in front of you rather than on a spec sheet.
//
// Memory: one warp is 32 consecutive values of col inside one row, so its 32
// addresses cover 128 contiguous bytes on the load and 128 on the store.
// Four 32-byte sectors each way.
//
// Launch assumption: the grid covers ceilDiv(cols, 32) by ceilDiv(rows, 32)
// tiles and the block is 32 by 8.
__global__ void copyRowMajor(const float* __restrict__ in,
float* __restrict__ out, size_t rows,
size_t cols) {
const size_t col = blockIdx.x * static_cast<size_t>(kTileDim) + threadIdx.x;
const size_t rowBase =
blockIdx.y * static_cast<size_t>(kTileDim) + threadIdx.y;
for (int j = 0; j < kTileDim; j += kBlockRows) {
const size_t row = rowBase + static_cast<size_t>(j);
if (row < rows && col < cols) {
out[row * cols + col] = in[row * cols + col];
}
}
}
// out[col][row] = in[row][col]. Day 12's kernel, reproduced here so all four
// rows of the table come out of one run on one card.
//
// Memory: the load is coalesced, 32 consecutive floats of one input row. The
// store is not. A warp holds one row index and 32 consecutive column
// indices, so its 32 output addresses are `rows` elements apart, which at
// 2003 floats is 8,012 bytes. Every lane lands in its own 32-byte sector and
// its own 128-byte line: 32 transactions where the load needed 4.
//
// There is no index order that fixes this. Move which subscript comes from
// threadIdx.x and the strided end moves from the store to the load. It does
// not go away, because out[c][r] = in[r][c] disagrees with itself.
//
// Launch assumption: same grid and block as copyRowMajor.
__global__ void transposeNaive(const float* __restrict__ in,
float* __restrict__ out, size_t rows,
size_t cols) {
const size_t col = blockIdx.x * static_cast<size_t>(kTileDim) + threadIdx.x;
const size_t rowBase =
blockIdx.y * static_cast<size_t>(kTileDim) + threadIdx.y;
for (int j = 0; j < kTileDim; j += kBlockRows) {
const size_t row = rowBase + static_cast<size_t>(j);
if (row < rows && col < cols) {
out[col * rows + row] = in[row * cols + col];
}
}
}
// out[col][row] = in[row][col], staged through a 32 x 32 tile in shared
// memory so that both global accesses are coalesced.
//
// One block owns one 32 x 32 square of the input. The first loop reads that
// square along its rows, 32 consecutive floats per warp, and writes it into
// the tile. The barrier makes every one of those writes visible to every
// thread in the block. The second loop reads the tile down its columns,
// which shared memory serves without a memory transaction, and writes 32
// consecutive floats of one output row.
//
// Memory: 4 sectors per warp on the global load and 4 on the global store,
// the same as copyRowMajor. The strided access still happens, it happens in
// shared memory, and shared memory has no sectors to waste.
//
// The output block sits at the transposed position, so col and rowBase are
// recomputed from the other blockIdx before the store. Reusing the input's
// col and row there gives a program that compiles, runs, and transposes
// nothing.
//
// Barrier rule: every thread reaches the __syncthreads(), including threads
// whose element is out of range. The guard covers the load, never the
// barrier, and nothing returns above it.
//
// Bank conflicts: the tile is [32][32] with no padding, so tile[x][y] for a
// fixed y and 32 consecutive x lands on one bank 32 times over. Day 15
// measures that and adds the column that fixes it.
//
// Launch assumption: same grid and block as copyRowMajor.
// snippet: tiled-kernel
__global__ void transposeTiledStatic(const float* __restrict__ in,
float* __restrict__ out, size_t rows,
size_t cols) {
__shared__ float tile[kTileDim][kTileDim];
size_t col = blockIdx.x * static_cast<size_t>(kTileDim) + threadIdx.x;
size_t rowBase = blockIdx.y * static_cast<size_t>(kTileDim) + threadIdx.y;
for (int j = 0; j < kTileDim; j += kBlockRows) {
const size_t row = rowBase + static_cast<size_t>(j);
if (row < rows && col < cols) {
tile[threadIdx.y + j][threadIdx.x] = in[row * cols + col];
}
}
__syncthreads();
col = blockIdx.y * static_cast<size_t>(kTileDim) + threadIdx.x;
rowBase = blockIdx.x * static_cast<size_t>(kTileDim) + threadIdx.y;
for (int j = 0; j < kTileDim; j += kBlockRows) {
const size_t row = rowBase + static_cast<size_t>(j);
if (row < cols && col < rows) {
out[row * rows + col] = tile[threadIdx.x][threadIdx.y + j];
}
}
}
// end snippet: tiled-kernel
// The same kernel with the tile sized at launch instead of at compile time.
//
// `extern __shared__` declares one unsized array per block. Its size is the
// third argument of the execution configuration, in bytes, which this course
// has left off every launch until now. The array is one dimensional because
// the compiler is not told the shape, so the [row][col] indexing above
// becomes an explicit multiply here. That is the whole cost of the dynamic
// form, and the whole benefit is that one kernel can serve tile sizes chosen
// by the caller.
//
// Launch assumption: same grid and block as copyRowMajor, plus at least
// kTileBytes of dynamic shared memory.
__global__ void transposeTiledDynamic(const float* __restrict__ in,
float* __restrict__ out, size_t rows,
size_t cols) {
extern __shared__ float tile[];
size_t col = blockIdx.x * static_cast<size_t>(kTileDim) + threadIdx.x;
size_t rowBase = blockIdx.y * static_cast<size_t>(kTileDim) + threadIdx.y;
for (int j = 0; j < kTileDim; j += kBlockRows) {
const size_t row = rowBase + static_cast<size_t>(j);
if (row < rows && col < cols) {
tile[(threadIdx.y + j) * kTileDim + threadIdx.x] =
in[row * cols + col];
}
}
__syncthreads();
col = blockIdx.y * static_cast<size_t>(kTileDim) + threadIdx.x;
rowBase = blockIdx.x * static_cast<size_t>(kTileDim) + threadIdx.y;
for (int j = 0; j < kTileDim; j += kBlockRows) {
const size_t row = rowBase + static_cast<size_t>(j);
if (row < cols && col < rows) {
out[row * rows + col] =
tile[threadIdx.x * kTileDim + threadIdx.y + j];
}
}
}
// Writes then reads back one float per requested element of a dynamically
// sized shared array, so nothing in it can be optimised away. Its only job
// is to be launched with more shared memory than the per-block default
// allows, once without permission and once with it.
//
// Memory: irrelevant. This kernel is a launch-time probe, not a benchmark,
// and its timing is never reported.
//
// Launch assumption: the third execution-configuration argument is at least
// floats * sizeof(float) bytes, and floats is at least one.
__global__ void touchDynamicShared(float* __restrict__ out, int floats) {
extern __shared__ float scratch[];
for (int k = static_cast<int>(threadIdx.x); k < floats;
k += static_cast<int>(blockDim.x)) {
scratch[k] = static_cast<float>(k);
}
__syncthreads();
if (threadIdx.x == 0) {
out[0] = scratch[floats - 1];
}
}
// CPU reference. Two nested loops in row-major order, written for obvious
// correctness rather than speed: no blocking, no OpenMP, no intrinsics. It
// never allocates; the caller owns every buffer.
//
// The loop nest and the kernels' grids are two independent spellings of one
// mapping, which is the only reason comparing them proves anything.
static void transposeCpu(const float* in, float* out, size_t rows,
size_t cols) {
for (size_t row = 0; row < rows; ++row) {
for (size_t col = 0; col < cols; ++col) {
out[col * rows + row] = in[row * cols + col];
}
}
}
// Returns the first index where got and want differ, or n if they agree
// everywhere. The comparison is exact rather than a tolerance: a transpose
// moves bits and computes nothing, so there is no multiply for the GPU to
// fuse and nothing for the two sides to round differently. A mismatch is an
// indexing bug and can be nothing else.
static size_t firstMismatch(const float* got, const float* want, size_t n) {
for (size_t i = 0; i < n; ++i) {
if (got[i] != want[i]) {
return i;
}
}
return n;
}
// Prints the first disagreement as a row and a column of the output, because
// "wrong at 6,005" names nothing a reader can look at and "wrong at row 2,
// column 0" names the first off-diagonal element. Returns 1 on a mismatch so
// main() can total the failures up and still take a single exit.
static int report(const char* name, const float* got, const float* want,
size_t outCols) {
const size_t n = kRows * kCols;
const size_t bad = firstMismatch(got, want, n);
if (bad == n) {
return 0;
}
std::fprintf(stderr,
"%s wrong at output row %zu, column %zu: got %.9g, want "
"%.9g\n",
name, bad / outCols, bad % outCols, got[bad], want[bad]);
return 1;
}
// 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; the alternative is four copies of the event boilerplate,
// which is how a warm-up goes missing from one of them.
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 kRows * kCols floats and writes kRows * kCols
// floats. That equality is the point of the program, so the figure is
// computed once from the constants rather than per kernel, where the two
// could drift apart.
static double bandwidthGBs(float ms) {
const double bytes =
2.0 * static_cast<double>(kRows * kCols) * sizeof(float);
return bytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// The per-block shared memory cap, and the call that lifts it.
//
// A kernel gets the default allowance without asking. Anything above it is
// opt-in, per kernel, and the refusal arrives at launch rather than at
// compile time, carrying an error string that says nothing about shared
// memory. Day 74's multi-stage pipeline needs more than the default, which is
// why the probe sits on the page that introduces __shared__.
//
// Returns 1 if the opt-in launch was also refused, because then the probe
// proved nothing and the run is not evidence for anything.
static int probeSharedMemoryCap(size_t defaultBytes, size_t optInBytes,
float* d_probe) {
const size_t requestBytes = defaultBytes + 4096;
if (requestBytes > optInBytes) {
std::printf(
"Skipping the shared memory probe: this device offers no opt-in\n"
"headroom above its %zu byte default.\n",
defaultBytes);
return 0;
}
const int floats = static_cast<int>(requestBytes / sizeof(float));
// Neither launch is wrapped in CUDA_CHECK. The first is expected to be
// refused and this function's job is to print the refusal, not to exit on
// it. cudaGetLastError() also clears the error, so the second check
// answers for its own launch and not for the first one.
// snippet: shared-optin
touchDynamicShared<<<1, kThreadsPerBlock, requestBytes>>>(d_probe, floats);
const cudaError_t beforeOptIn = cudaGetLastError();
std::printf("%zu bytes of dynamic shared memory, no opt-in: %s\n",
requestBytes, cudaGetErrorString(beforeOptIn));
// cudaFuncSetAttribute's third parameter is an int, so the narrowing is
// written out rather than left to the call. requestBytes is bounded by
// optInBytes a few lines up, so it fits.
CUDA_CHECK(cudaFuncSetAttribute(touchDynamicShared,
cudaFuncAttributeMaxDynamicSharedMemorySize,
static_cast<int>(requestBytes)));
touchDynamicShared<<<1, kThreadsPerBlock, requestBytes>>>(d_probe, floats);
const cudaError_t afterOptIn = cudaGetLastError();
std::printf("%zu bytes after cudaFuncSetAttribute: %s\n",
requestBytes, cudaGetErrorString(afterOptIn));
// end snippet: shared-optin
if (beforeOptIn == cudaSuccess || afterOptIn == cudaSuccess) {
CUDA_CHECK(cudaDeviceSynchronize());
}
if (beforeOptIn == cudaSuccess) {
std::printf(
"This device accepted %zu bytes with no opt-in, so the default\n"
"cap is not where this lesson says it is. Say so on the page.\n",
requestBytes);
}
if (afterOptIn != cudaSuccess) {
std::fprintf(stderr,
"the opt-in launch was refused as well, so the probe "
"proved nothing\n");
return 1;
}
return 0;
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
const size_t n = kRows * kCols;
const size_t bytes = n * sizeof(float);
const double bufferMiB = static_cast<double>(bytes) / (1024.0 * 1024.0);
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name, prop.major,
prop.minor);
std::printf("Matrix: %zu x %zu floats, %.1f MiB per buffer\n", kRows, kCols,
bufferMiB);
// Both members are size_t, so both take %zu. Passing a size_t to %d is
// undefined behaviour in a varargs call and only prints the right number
// by little-endian accident.
std::printf(
"Shared memory per block: %zu bytes default, %zu bytes opt-in\n",
prop.sharedMemPerBlock, prop.sharedMemPerBlockOptin);
std::printf("Shared memory per SM: %zu bytes\n",
prop.sharedMemPerMultiprocessor);
std::printf("Tile: %d x %d floats, %zu bytes per block\n\n", kTileDim,
kTileDim, kTileBytes);
// Values are row * 1000 + col, so every element names its own position
// and a swapped subscript produces a number that could not have come
// from that cell. The largest value is 2,005,000, well inside the
// 16,777,216 integers a float represents exactly, so the comparison
// below can be exact.
std::vector<float> h_in(n);
std::vector<float> h_want(n);
std::vector<float> h_got(n);
for (size_t row = 0; row < kRows; ++row) {
for (size_t col = 0; col < kCols; ++col) {
h_in[row * kCols + col] = static_cast<float>(row * 1000 + col);
}
}
transposeCpu(h_in.data(), h_want.data(), kRows, kCols);
float* d_in = nullptr;
float* d_out = nullptr;
float* d_probe = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
CUDA_CHECK(cudaMalloc(&d_probe, sizeof(float)));
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), bytes, cudaMemcpyHostToDevice));
const dim3 block(kTileDim, kBlockRows);
const dim3 grid(static_cast<unsigned int>(
ceilDiv(kCols, static_cast<size_t>(kTileDim))),
static_cast<unsigned int>(
ceilDiv(kRows, static_cast<size_t>(kTileDim))));
// Correctness first, one launch each, every result compared against a
// reference computed on the host. Failures are counted rather than
// returned from here, so that one cleanup block below runs on every
// path.
int failures = 0;
copyRowMajor<<<grid, block>>>(d_in, d_out, kRows, kCols);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_got.data(), d_out, bytes, cudaMemcpyDeviceToHost));
failures += report("copyRowMajor", h_got.data(), h_in.data(), kCols);
transposeNaive<<<grid, block>>>(d_in, d_out, kRows, kCols);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_got.data(), d_out, bytes, cudaMemcpyDeviceToHost));
failures += report("transposeNaive", h_got.data(), h_want.data(), kRows);
transposeTiledStatic<<<grid, block>>>(d_in, d_out, kRows, kCols);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_got.data(), d_out, bytes, cudaMemcpyDeviceToHost));
failures +=
report("transposeTiledStatic", h_got.data(), h_want.data(), kRows);
// The third execution-configuration argument, for the first time in this
// course: bytes of dynamic shared memory per block. Leaving it off gives
// a kernel whose extern __shared__ array has length zero, and the first
// write runs off the end of it.
transposeTiledDynamic<<<grid, block, kTileBytes>>>(d_in, d_out, kRows,
kCols);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_got.data(), d_out, bytes, cudaMemcpyDeviceToHost));
failures +=
report("transposeTiledDynamic", h_got.data(), h_want.data(), kRows);
// Only launches go inside the timed regions. The allocations and the
// copies are above, so these numbers are kernel time and nothing else.
const float msCopy = timeKernel(
[&] { copyRowMajor<<<grid, block>>>(d_in, d_out, kRows, kCols); });
const float msNaive = timeKernel(
[&] { transposeNaive<<<grid, block>>>(d_in, d_out, kRows, kCols); });
const float msStatic = timeKernel([&] {
transposeTiledStatic<<<grid, block>>>(d_in, d_out, kRows, kCols);
});
const float msDynamic = timeKernel([&] {
transposeTiledDynamic<<<grid, block, kTileBytes>>>(d_in, d_out, kRows,
kCols);
});
std::printf("%-24s %11s %11s %10s\n", "kernel", "time (ms)", "GB/s",
"vs copy");
std::printf("%-24s %11s %11s %10s\n", "------------------------",
"----------", "----------", "---------");
std::printf("%-24s %11.3f %11.1f %9.2fx\n", "copyRowMajor", msCopy,
bandwidthGBs(msCopy), 1.0);
std::printf("%-24s %11.3f %11.1f %9.2fx\n", "transposeNaive", msNaive,
bandwidthGBs(msNaive), msNaive / msCopy);
std::printf("%-24s %11.3f %11.1f %9.2fx\n", "transposeTiledStatic",
msStatic, bandwidthGBs(msStatic), msStatic / msCopy);
std::printf("%-24s %11.3f %11.1f %9.2fx\n", "transposeTiledDynamic",
msDynamic, bandwidthGBs(msDynamic), msDynamic / msCopy);
std::printf(
"\nEvery row above moved exactly %.1f MiB. All four launch the same\n"
"grid and the same block and read and write every element once, so\n"
"no row can win by doing less work.\n\n",
2.0 * bufferMiB);
failures += probeSharedMemoryCap(prop.sharedMemPerBlock,
prop.sharedMemPerBlockOptin, d_probe);
// One cleanup block, reached on every path, before the single exit. The
// failure count is a real branch and not an assert(): CI builds Release,
// Release defines NDEBUG, and NDEBUG deletes assert(), so a check
// written that way is absent from exactly the build that matters.
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_out));
CUDA_CHECK(cudaFree(d_probe));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed, listed above\n", failures);
return EXIT_FAILURE;
}
std::printf("all four kernels match the CPU reference at %zu elements\n",
n);
return EXIT_SUCCESS;
}