code/day85-python/vector_add.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 85: the C++ half of "one kernel, four Python libraries".
//
// This file is the reference implementation and the Compiler Explorer
// artifact. four_ways.py reads the kernel out of this file by its snippet
// markers and hands the same characters to cuda.core, CuPy and PyCUDA, so
// there is exactly one copy of the source in the day and the four launchers
// cannot silently drift from it. Numba is the exception: its kernel is
// written in Python and is a reimplementation, which is why the lesson keeps
// it in its own column.
//
// The inputs come from an index hash rather than an RNG so that this program
// and the NumPy reference in four_ways.py produce bit-identical bytes without
// sharing a seed or a file. Values land in [1, 2), so every sum lands in
// [2, 4) and the addition rounds off a real bit instead of being exact.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o vector_add vector_add.cu
// Run: ./vector_add
//
// VERIFIED: Tesla T4, driver 580.173.02, CUDA 12.6, 2026-09-02. See
// evidence/cpp-2026-09-02.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)
// 4 Mi elements plus 611, so the bounds check runs on every launch and each
// buffer is 16 MiB. Three buffers at that size fit inside what Compiler
// Explorer will run in 20 seconds, and four_ways.py uses the same count so
// the two programs' GB/s columns mean the same thing.
constexpr size_t kElems = (1ull << 22) + 611ull;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kRelTolerance = 1e-5f;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
// out[i] = a[i] + b[i]. One thread owns one element.
//
// Memory: consecutive threads take consecutive elements, so one warp's 32
// addresses cover 128 contiguous bytes on each of the three buffers. The
// kernel moves 3n floats and does n flops, which is why the number worth
// printing is GB/s and not GFLOP/s.
//
// Launch assumption: gridDim.x * blockDim.x >= n. The guard is a branch, not
// an early return, because an early return above a barrier hangs a kernel and
// the habit is easier to never form than to unlearn.
// snippet: kernel
__global__ void vectorAdd(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
// Knuth's multiplicative hash, one xor-shift, then the low 23 bits read as a
// mantissa. Returns a float in [1, 2) with all 23 mantissa bits set from the
// index. The scale is 2^-23, a power of two, so the multiply is exact and
// NumPy's uint32 arithmetic reproduces this bit for bit.
static float sample(size_t i) {
unsigned int h = static_cast<unsigned int>(i) * 2654435761u;
h ^= h >> 15;
return 1.0f + static_cast<float>(h & 0x7FFFFFu) * (1.0f / 8388608.0f);
}
// 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. 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.
// Copied verbatim from the course's style reference; the warm-up lives inside
// it so no caller can forget one.
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;
}
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)\n", prop.name, prop.major,
prop.minor);
const size_t bytes = kElems * sizeof(float);
const int blocks =
static_cast<int>((kElems + kThreadsPerBlock - 1) / kThreadsPerBlock);
std::printf(
"n = %zu, %d threads per block, %d blocks, %zu MiB per "
"buffer\n",
kElems, kThreadsPerBlock, blocks, bytes >> 20);
std::vector<float> h_a(kElems);
std::vector<float> h_b(kElems);
std::vector<float> h_out(kElems);
std::vector<float> h_want(kElems);
for (size_t i = 0; i < kElems; ++i) {
h_a[i] = sample(2 * i);
h_b[i] = sample(2 * i + 1);
}
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));
vectorAdd<<<blocks, kThreadsPerBlock>>>(d_a, d_b, d_out, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
vectorAddCpu(h_a.data(), h_b.data(), h_want.data(), kElems);
const size_t bad =
firstMismatch(h_out.data(), h_want.data(), kElems, kRelTolerance);
// Only the launches sit inside the timed region. Allocation and both
// copies are above it, so this number is kernel time and nothing else.
// 3n floats move: two reads and one write.
float ms = 0.0f;
if (bad == kElems) {
ms = timeKernel([&] {
vectorAdd<<<blocks, kThreadsPerBlock>>>(d_a, d_b, d_out, kElems);
});
}
// Free before reporting, so the failure path frees too. Every cudaMalloc
// has a matching cudaFree before every return, including the one taken
// when the answer is wrong.
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_out));
if (bad != kElems) {
std::fprintf(stderr, "vectorAdd wrong at %zu: got %.9g, want %.9g\n",
bad, h_out[bad], h_want[bad]);
return EXIT_FAILURE;
}
const double movedBytes = 3.0 * static_cast<double>(kElems) * sizeof(float);
std::printf("all %zu elements match the CPU reference\n", kElems);
std::printf(
"vectorAdd: %.4f ms, %.1f GB/s (mean of %d runs, copies not "
"included)\n",
ms, movedBytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9,
kTimedRuns);
return EXIT_SUCCESS;
}