vector_add.cu

This 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 5: vector addition, end to end, with the error checking already in.
//
// Nothing here is timed. A host clock around a first CUDA program measures
// context creation and two PCIe copies, not the kernel, so there is no
// cudaEvent and no clock below. Day 9 takes that apart.
//
// README.md has the three edits the lesson asks you to make, and says where
// each resulting failure surfaces.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o vector_add vector_add.cu
// Run:   ./vector_add
//
// 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`.
// snippet: check-macro
#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)
// end snippet

// 611 is 13 x 47, so no block size divides it. At 256 threads per block that
// is 3 blocks, 768 threads and 157 threads with no element, which is what
// makes the kernel's bounds check run on every launch instead of never.
// kRelTolerance is not zero because a GPU may fuse a multiply and an add
// where the host compiler does not; on these inputs both sides are exact.
constexpr size_t kElems = 611;
constexpr int kThreadsPerBlock = 256;  // 8 warps
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. Day 11 measures what that is worth;
// nothing here depends on it yet.
//
// Launch assumption: gridDim.x * blockDim.x >= n. No early return above the
// guard: that habit hangs a kernel once a __syncthreads() sits below it,
// which happens on day 13.
// 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

// 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 whole 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;
}

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);
    const size_t launched = static_cast<size_t>(blocks) * kThreadsPerBlock;

    std::printf("n = %zu, %d threads per block, %d blocks, %zu threads\n",
                kElems, kThreadsPerBlock, blocks, launched);
    std::printf("%zu threads have no element to add\n", launched - kElems);

    // Every host pointer carries h_ and every device pointer d_, because the
    // prefix is the only thing standing between you and passing one where the
    // other belongs. The values are small whole numbers, so every sum is
    // exact and a mismatch below can only be an index bug.
    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] = static_cast<float>(i);
        h_b[i] = static_cast<float>(2 * i);
    }

    // cudaMalloc takes the address of the pointer, not the pointer, because
    // it writes the new device address back through it. The address it writes
    // means something to the GPU and nothing to this process: reading d_a[0]
    // here is not slow, it is a segmentation fault.
    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));

    // Destination, source, bytes, direction. The direction has to agree with
    // the pointers; the runtime documents a mismatch as undefined behaviour.
    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);

    // The launch returns void, so these two lines are the only way to hear
    // about it. cudaGetLastError() reports a launch the driver refused, such
    // as a block size above 1024; cudaDeviceSynchronize() waits and reports
    // what the kernel did, such as an illegal address. They catch different
    // things and skipping either is how a wrong number reaches your screen.
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    // A device-to-host cudaMemcpy is itself a sync point, so the synchronize
    // above is not needed for correctness. It is there to attribute a kernel
    // failure to the kernel's line instead of to this one.
    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);

    // Free before reporting, so the failure path frees too. Every cudaMalloc
    // has a matching cudaFree before every return, including the one you take
    // when the answer is wrong. This is the pairing the lesson is about, and
    // an early return that skips it would teach the opposite.
    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;
    }
    std::printf("all %zu elements match the CPU reference\n", kElems);
    return EXIT_SUCCESS;
}