GPU correctness and honest benchmarks

gpu-kernel-basics got a kernel to run. This skill is about the two questions everyone asks next — in code review, in an interview, in the first week of a rendering job:

1. Is it right? Does the GPU output match a trusted CPU reference, within a tolerance you can justify?
2. Is it faster, and by how much? Measured how, on what, including what?

The order matters. A fast wrong answer isn't a speedup, and a speedup nobody can reproduce isn't evidence. The habits here — compare with documented tolerances, time with synchronisation and warm-up, summarise with medians, write down hardware, toolkit, sizes and method, and report "unavailable" or "slower" plainly — are what make your numbers trusted. They build on profiling-benchmarking (warm-up, min/median, perf), which this lesson doesn't repeat.

1. Why GPU and CPU results differ

Bit-identical output between a CPU loop and a GPU kernel is the exception. These are the usual, legitimate reasons:

• Reassociation. A CPU loop adds left to right; a GPU reduction adds in a tree (pairs, then pairs of pairs, then per-block partial sums). Float addition isn't associative, so the two give different roundings.
• FMA contraction. nvcc fuses a * b + c into one fused multiply-add (-fmad=true is the default): one rounding instead of two. A CPU build contracts only if the compiler is allowed to and the target has FMA (-march=x86-64-v3 with GCC's default -ffp-contract=fast in GNU mode).
• Fast maths. --use_fast_math swaps in intrinsics (__expf, __sinf, __fdividef), flushes denormals to zero (-ftz=true) and uses approximate division and square root. Even without it, CUDA's expf, sinf… have documented ULP error bounds (the CUDA Programming Guide lists them) that differ from glibc's.
• Atomics. atomicAdd on floats accumulates in whatever order threads arrive — the result changes from run to run in the last bits.

And these are bugs that a tolerance must not hide: a missing bounds check, a shared-memory race without __syncthreads(), uninitialised device memory. compute-sanitizer (--tool memcheck, racecheck, initcheck, synccheck) finds them; run it on the correctness test, never on the timed run.

Here is the effect on real numbers. The pairwise loop is what a block reduction does; the double sum is a better reference than either float version:

#include <print>
#include <vector>

float pairwise(std::vector<float> level) {
    while (level.size() > 1) {
        std::vector<float> next((level.size() + 1) / 2);
        for (std::size_t i = 0; i < next.size(); ++i)
            next[i] = level[2 * i] + (2 * i + 1 < level.size() ? level[2 * i + 1] : 0.0f);
        level = std::move(next);
    }
    return level[0];
}

int main() {
    std::vector<float> v(1'000'000, 0.1f);

    float forward = 0.0f;
    for (float x : v) forward += x;

    float kahan = 0.0f, carry = 0.0f;  // compensated summation
    for (float x : v) {
        float y = x - carry;
        float t = kahan + y;
        carry = (t - kahan) - y;
        kahan = t;
    }

    double reference = 0.0;
    for (float x : v) reference += x;

    std::println("forward  {:.4f}", forward);
    std::println("pairwise {:.4f}", pairwise(v));
    std::println("kahan    {:.4f}", kahan);
    std::println("double   {:.4f}", reference);
}

The naive CPU loop is the worst of the four (off by almost 1%). When a GPU reduction disagrees with a float CPU loop, don't assume the GPU is wrong: compute the reference in double (or with Kahan / pairwise summation), and compare both float implementations against it.

2. Tolerances: absolute, relative, ULP

Three ways to say "close enough", each right in a different regime:

Production code combines them: NumPy's assert_allclose and PyTorch's torch.testing.assert_close accept an element when it passes the absolute or the relative test; Google Test's EXPECT_FLOAT_EQ allows 4 ULPs.

ULP distance from bit patterns

For floats of the same sign, IEEE 754 ordering matches the ordering of the bit patterns read as integers — adjacent floats have adjacent patterns. So the number of floats between a and b is a subtraction, once negative floats (sign bit + magnitude) are folded onto the same line:

The special values need explicit rules, and they're where comparison code goes wrong:

• NaN never matches, not even NaN. Two identical NaN patterns have key distance 0 — check for NaN before any bit arithmetic.
• ±0 are equal: a == b handles it, and so does the key.
• Infinities match only the same infinity. inf - inf is NaN, so test a == b first, then reject any remaining infinity.

#include <algorithm>
#include <bit>
#include <cmath>
#include <cstdint>
#include <limits>
#include <print>

std::int64_t key(float x) {
    std::uint32_t u = std::bit_cast<std::uint32_t>(x);
    std::int64_t m = u & 0x7fffffffu;
    return (u >> 31) ? -m : m;
}

std::uint64_t ulps(float a, float b) {
    if (std::isnan(a) || std::isnan(b)) return std::numeric_limits<std::uint64_t>::max();
    std::int64_t d = key(a) - key(b);
    return static_cast<std::uint64_t>(d < 0 ? -d : d);
}

bool close(float a, float b, double atol, double rtol, std::uint64_t max_ulps) {
    if (a == b) return true;                           // ±0, equal infinities
    if (std::isinf(a) || std::isinf(b)) return false;  // any other infinity
    double d = std::abs(double(a) - double(b));        // NaN fails every <=
    return d <= atol || d <= rtol * std::max(std::abs(double(a)), std::abs(double(b))) ||
           ulps(a, b) <= max_ulps;
}

int main() {
    float tiny = std::numeric_limits<float>::denorm_min();
    float nan = std::numeric_limits<float>::quiet_NaN();
    float inf = std::numeric_limits<float>::infinity();
    std::println("ulps(1, next 1)   = {}", ulps(1.0f, std::nextafter(1.0f, 2.0f)));
    std::println("ulps(-tiny, tiny) = {}", ulps(-tiny, tiny));
    std::println("ulps(1, 2)        = {}  (a whole binade: 2^23)", ulps(1.0f, 2.0f));
    std::println("close(1e-9, 0)    rel only: {}, abs 1e-6: {}", close(1e-9f, 0.0f, 0, 1e-5, 0),
                 close(1e-9f, 0.0f, 1e-6, 1e-5, 0));
    std::println("close(nan, nan) = {}, close(inf, inf) = {}, close(inf, max) = {}", close(nan, nan, 1, 1, 1000),
                 close(inf, inf, 0, 0, 0), close(inf, std::numeric_limits<float>::max(), 1e38, 1, 1000));
}

Choosing and documenting a tolerance

A tolerance is a claim, so derive it rather than widening it until the test goes green:

1. Compute a high-precision reference (double, or the CPU code with compensated sums) for a few representative inputs.
2. Measure the error distribution of both the float CPU path and the GPU path against it: worst absolute error, worst ULPs, and where they occur.
3. Set the tolerance a small factor above what you observed, with a reason tied to the algorithm ("sum of up to 4096 splat contributions per pixel; reassociation bound ≈ n·ε·Σ|x|"), and write it next to the test.
4. Report the worst element — its index, both values and its error — on failure, not just "1342 mismatches". The worst element is where you start debugging.

// tests/composite_gpu_test.cpp — the tolerance is part of the test's documentation.
// Per-pixel alpha compositing of ≤ 4096 splats: reassociation and FMA give
// ≤ 64 ULP on the fixtures (measured: worst 19 ULP); |x| < 1e-6 is black.
constexpr Tolerance kCompositeTol{.abs = 1e-6, .rel = 0.0, .ulps = 64};

When results must be bitwise reproducible (golden images in CI, a determinism requirement), don't widen tolerances — make the GPU deterministic: fixed-order reductions instead of float atomics, -fmad=false if needed, no fast maths, and the same build on both sides. That costs speed; say so.

3. Measuring GPU time honestly

Launches are asynchronous

A kernel launch returns as soon as the work is queued. Host timers around a launch measure the launch (microseconds), not the kernel:

Two correct ways to time:

// (a) GPU-side: events in the same stream as the kernel. Measures device time only.
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, stream);
composite<<<grid, block, 0, stream>>>(d_splats, d_image, n);
cudaEventRecord(stop, stream);
cudaEventSynchronize(stop);                  // wait until stop has happened
float ms = 0.0f;
cudaEventElapsedTime(&ms, start, stop);      // ~0.5 µs resolution
if (cudaError_t e = cudaGetLastError(); e != cudaSuccess) report_launch_error(e);

// (b) End-to-end, what a user waits for: upload + kernel + download, wall clock.
auto t0 = std::chrono::steady_clock::now();
cudaMemcpyAsync(d_splats, h_splats, bytes_in, cudaMemcpyHostToDevice, stream);
composite<<<grid, block, 0, stream>>>(d_splats, d_image, n);
cudaMemcpyAsync(h_image, d_image, bytes_out, cudaMemcpyDeviceToHost, stream);
cudaStreamSynchronize(stream);
auto t1 = std::chrono::steady_clock::now();

Report both when they differ much — "kernel 4.1 ms, end-to-end 11.8 ms" tells the reader where the time goes. A CPU-vs-GPU comparison that a user cares about is end-to-end; leaving transfers out flatters the GPU. (Use pinned host memory, cudaMallocHost, or the async copies are not truly asynchronous — and say which you used.)

Warm-up is bigger on a GPU

The first iterations pay for things that never recur: creating the CUDA context (hundreds of ms), loading the module and JIT-compiling PTX if the binary lacks SASS for this GPU, first-touch allocation, cold L2, and GPU clocks ramping up from idle. Run several untimed iterations, then time many. Drop warm-up by position (the first runs), never by value — dropping the largest samples after sorting is a different, biased statistic.

Summarising: median, MAD, a high percentile

profiling-benchmarking used the minimum for small deterministic CPU functions. For GPU work the run-to-run spread is part of the story (clocks, contention, other processes on the GPU), so report a robust centre and a spread:

• median — the typical run; insensitive to a few hiccups.
• MAD — median absolute deviation from the median; 1.4826 × MAD estimates the standard deviation for normal data, so "more than 3 scaled MADs away" is a robust outlier rule.
• p90 / p99 — the tail, which is what a frame-time budget cares about. Percentiles have several definitions; say which (NumPy's default "linear" interpolates between closest ranks).

#include <algorithm>
#include <cmath>
#include <print>
#include <vector>

double median_sorted(const std::vector<double>& v) {
    std::size_t m = v.size() / 2;
    return v.size() % 2 ? v[m] : (v[m - 1] + v[m]) / 2.0;
}

int main() {
    std::vector<double> run_order{180, 9.5, 5.1, 4.1, 4.0, 4.2, 4.1, 4.0, 6.3, 4.1, 4.0, 4.2};
    const std::size_t warmup = 3;

    std::vector<double> x(run_order.begin() + warmup, run_order.end());  // by position
    std::ranges::sort(x);
    double med = median_sorted(x);

    std::vector<double> dev;
    for (double v : x) dev.push_back(std::abs(v - med));
    std::ranges::sort(dev);
    double mad = median_sorted(dev);

    double h = 0.9 * static_cast<double>(x.size() - 1);  // NumPy 'linear'
    std::size_t lo = static_cast<std::size_t>(h), hi = std::min(lo + 1, x.size() - 1);
    double p90 = x[lo] + (h - static_cast<double>(lo)) * (x[hi] - x[lo]);

    double mean = 0;
    for (double v : x) mean += v;
    mean /= static_cast<double>(x.size());

    int outliers = 0;
    for (double v : x) outliers += std::abs(v - med) > 3 * 1.4826 * mad;

    std::println("n={} median={:.2f} MAD={:.2f} p90={:.2f} mean={:.2f} outliers={}", x.size(), med, mad, p90, mean,
                 outliers);
}

The mean (4.33) is pulled up by the one 6.3 ms hiccup; the median (4.10) isn't. The outlier rule flags that hiccup.

Memory is a result too

Peak device memory often decides whether a renderer fits a scene. Measure it: cudaMemGetInfo(&free, &total) before and after (coarse: other processes count), or allocate with cudaMallocAsync from a memory pool and read cudaMemPoolGetAttribute(pool, cudaMemPoolAttrUsedMemHigh, &bytes) — the pool's high-water mark. Report it with the sizes it was measured at.

4. Profilers: Nsight Systems first, Nsight Compute second

Start with the timeline. If a frame takes 30 ms and the kernel 4 ms, no amount of kernel tuning helps; the timeline shows the synchronous copies, per-frame cudaMalloc, or idle gaps while the CPU prepares the next launch. Mark your own phases with NVTX ranges so they appear on the timeline:

#include <nvtx3/nvtx3.hpp>
void render_frame(const Scene& s) {
    nvtx3::scoped_range r{"render_frame"};
    { nvtx3::scoped_range u{"upload"};   upload(s); }
    { nvtx3::scoped_range k{"composite"}; launch_composite(s); }
    { nvtx3::scoped_range d{"download"}; download(s); }
}

--stats=true prints summary tables after the run (illustrative excerpt):

 ** CUDA GPU Kernel Summary (cuda_gpu_kern_sum):
 Time (%)  Total Time (ns)  Instances  Avg (ns)    Med (ns)    Name
 --------  ---------------  ---------  ----------  ----------  ----------------
     91.2      123,004,112         30   4,100,137   4,098,560  composite(...)
      8.8       11,870,240         30     395,675     394,880  project_points(...)

 ** CUDA GPU MemOps Summary (by Time) (cuda_gpu_mem_time_sum):
 Time (%)  Total Time (ns)  Count  Operation
 --------  ---------------  -----  ------------------
     71.4      214,551,808     30  [CUDA memcpy HtoD]
     28.6       85,904,384     30  [CUDA memcpy DtoH]

Here copies take more than twice the kernel time — the next optimisation is keeping data resident on the GPU, not the kernel.

Nsight Compute then explains one kernel: the SpeedOfLight section says how close it is to peak memory bandwidth or compute, MemoryWorkloadAnalysis shows whether accesses are coalesced, Occupancy and the warp-state statistics show why warps stall. It replays the kernel several times, by default locks clocks to base (--clock-control base) and flushes caches between passes, so its durations are for analysis — don't put them in a timing table next to your event timings.

5. Reporting

A benchmark report is short and complete. The four facts a reader needs are hardware, toolkit (driver, CUDA, compilers and flags), sizes, and method (what was timed, warm-up, number of runs, statistic, spread) — plus the correctness result and the limitations you know about:

# Splat compositing, CPU vs GPU
- Hardware: Ryzen 9 7950X (16 cores); RTX 4070 12 GB
- Toolkit: CUDA 13.0, driver 580; GCC 16 -O2 -march=native; nvcc -O3 -arch=sm_89
- Sizes: 1,000,000 splats, 1920x1080, fixture garden-1m
- Method: 5 warm-up + 30 timed runs; median and MAD; GPU end-to-end incl. transfers (pinned)
- Correctness: PASS, all 2,073,600 pixels within abs 1e-6 / 64 ULP; worst 19 ULP at (811, 402)
- Result: GPU is 8.00x faster than CPU (median 48.0 ms vs 6.0 ms)
- Limitations: one scene; laptop on AC power; CPU path single-threaded

Rules that keep a report honest:

• No correctness, no speedup. If the comparison failed, say so and don't state a speedup.
• Slower is a result. "GPU is 2.4x slower at 10k splats (transfers dominate), 6x faster at 2M" is more useful than either number alone.
• Unavailable is not passed. No GPU, no driver, nvcc but no device: the GPU section says unavailable with the reason, the tests are skipped (CTest's SKIP_RETURN_CODE, GoogleTest's GTEST_SKIP()), and no speedup is stated — even if old numbers are lying around.
• Ratios go one way. Speedup = old time / new time. Write the two medians next to the ratio so the reader can check the direction.

This app's sandbox has no GPU and no nvcc, and your machine may have nvcc 13 without a working driver right now. That is exactly the "unavailable" case: build the CUDA code, run the CPU reference and the report tooling, and record the GPU section as unavailable until hardware is there. The build task implements this report library.

6. What interviews probe

• "Your GPU result differs from the CPU in the 6th digit. Is it a bug?" — reassociation/FMA/fast-math vs races and bounds; how you'd decide (double reference, worst element, ULP distribution, compute-sanitizer).
• "How did you time it?" — events vs synchronised wall clock, warm-up, transfers in or out, how many runs, which statistic.
• "Why is it slower than you expected?" — the timeline first (nsys), then the kernel (ncu), and a bandwidth estimate: bytes moved / time against the GPU's peak.

In the project

Ungraded practice in projects/scene-viewer/gpu/ (stage 5); the evidence is your own build output and, when you have a GPU, your report.

1. Read gpu/bench.cpp: the per-field tolerance, counting visibility mismatches on culling boundaries separately, CUDA-event phases plus end-to-end wall time, and the "unavailable" path that exits 77.
2. Replace its median with your build task's statistics (median, MAD, p90, warm-up dropped) and add a ULP-based comparison for depth.
3. On a machine with an NVIDIA GPU, save the Markdown report with the hardware, driver and toolkit it prints; explain why the GPU can be slower end to end at 10 000 splats. Capture nsys profile --stats=true and note where the time goes.