Case Study 2: A Reproducible Benchmark for the Heat Solver

"If you can not measure it, you can not improve it." — commonly attributed to Lord Kelvin (William Thomson)

Executive Summary

Chapter 29 is about to make the heat solver's stencil faster, and Chapter 38 will present the whole project as a paper with a performance section. Both need a single, defensible number: the solver's serial baseline per-step time, measured so carefully that when you later claim "the optimized version is $2.4\times$ faster," a skeptical reviewer cannot wave it away. A lone toc() reading cannot carry that weight — it is one sample from a noisy distribution, taken on a machine that was doing other things.

This case study builds that measurement instrument. You will design a benchmark harness around the project's timers module — warm-up, repetition, a working-copy reset between runs, a robust statistic with a reported spread, and a guard against the compiler deleting the very loop you are timing — and then use cachegrind to diagnose why the stencil is slow, handing Chapter 29 a precise target: a memory-bound loop with a bad-loop-order cache signature. The deliverable is not a faster solver (that is next chapter); it is a trustworthy, reproducible baseline and a diagnosis, which is the only honest foundation for optimizing.

Skills applied

• Building a warm-up / repeat / robust-statistic harness on the timers module (§28.1, §28.5).
• Choosing and reporting the minimum, median, and spread of repeated timings (§28.5).
• Defeating dead-code elimination so the benchmark measures real work (§28.5).
• Diagnosing a memory-bound kernel with valgrind --tool=cachegrind and cg_annotate (§28.3, §28.4).
• Recording the environment (compiler, flags, machine) for reproducibility (§28.5; forward to Ch. 30, 37).

Background

The solver's time-stepping loop calls step (the frozen Chapter 6 signature) n_steps times; step's cost is dominated by the five-point stencil inside laplacian (Case Study 1's discipline would have you confirm that with gprof first — and §28.2's representative profile does: ~78% in laplacian). The figure of merit is seconds per step. The naïve approach — tic; loop; toc; divide — produces a number, but not a defensible one. We will fix that in five phases.

Phase 1 — Warm up, repeat, and reset the state

A benchmark that runs the timed code once measures start-up as much as steady state. The harness warms up, then repeats, and — because each solver run mutates the field — restores a clean working copy before every repetition. That reset is the subtlety a first attempt always misses: time the same work each rep, from the same initial state.

subroutine bench_solver(u0, alpha, dt, n_steps, n_rep, times)
  use timers, only: tic, toc
  real(dp), intent(in)  :: u0(:,:), alpha, dt
  integer,  intent(in)  :: n_steps, n_rep
  real(dp), intent(out) :: times(n_rep)
  real(dp), allocatable :: u(:,:)
  integer :: rep, n
  ! warm-up: one full untimed run to prime caches, pages, and CPU frequency
  u = u0
  do n = 1, n_steps
    call step_raw(u, alpha, dt)
  end do
  ! timed repetitions, each from a fresh copy of the initial field
  do rep = 1, n_rep
    u = u0                       ! RESET: identical work every repetition
    call tic()
    do n = 1, n_steps
      call step_raw(u, alpha, dt)
    end do
    times(rep) = toc()
  end do
end subroutine bench_solver

Here step_raw is the interior update on a bare array (the field_t wrapper adds nothing to time); times comes back with n_rep samples ready for statistics. The per-step time is times(rep)/n_steps, but we do not reduce to a single number yet — that is Phase 2.

Correctness anchor. While benchmarking, keep a cheap invariant check so you never optimize a broken run: the held hot edge means maxval(u) must stay exactly 100.0 (the discrete maximum principle, r ≤ 1/4). If a "faster" variant ever changes that number, it changed the physics — stop and diagnose. Correctness is not negotiable for speed.

Phase 2 — Reduce to a statistic, with its spread

Given the sample vector, report the minimum (the least-disturbed run — the truest cost of the code itself), the median (typical behaviour, robust to a single hiccup), and the spread (so the reader knows how noisy the machine was). On a known sample set the arithmetic is fully determinate:

real(dp) :: t(5) = [0.100_dp, 0.102_dp, 0.098_dp, 0.101_dp, 0.099_dp]
! sorted: 0.098, 0.099, 0.100, 0.101, 0.102
print '(a, f6.3)', 'min    = ', minval(t)        ! 0.098
print '(a, f6.3)', 'median = ', t_median(t)      ! 0.100  (middle of five)
print '(a, f6.3)', 'range  = ', maxval(t)-minval(t)  ! 0.004

min    =  0.098
median =  0.100
range  =  0.004

A range that is a few percent of the median (here $0.004/0.100 = 4\%$) is a healthy, quiet measurement; a range comparable to the median means the machine was too busy to trust — go quiet it and re-run. Reporting the min and the median and the spread is three numbers instead of one, and it is the difference between a measurement and a guess with good posture.

Phase 3 — Defeat dead-code elimination and pin the environment

Two ways a benchmark lies about being fast. First, if the timed loop's result is never used, an optimizing compiler may delete it — you would time nothing and celebrate. The solver is safe here because u is read back by the maxval check, but make the dependence explicit and unpredictable so no future refactor can strand it:

real(dp) :: sink
sink = 0.0_dp
do rep = 1, n_rep
  u = u0
  call tic()
  do n = 1, n_steps
    call step_raw(u, alpha, dt)
  end do
  times(rep) = toc()
  sink = sink + u(2,2)          ! consume a live result each rep
end do
print *, times, sink           ! printing sink keeps every run's work alive

Second, the environment must be controlled and recorded. Benchmark the build you will ship (-O2 or -O3, never -O0), on a quiet machine, ideally with a fixed CPU frequency, and change one variable at a time. Then write down what you did, because a number without its conditions is unreproducible:

That table is the reproducibility habit of Chapter 37, applied to performance: the "works on my machine" problem is a measurement problem too.

Phase 4 — Diagnose why: cachegrind on the stencil

The baseline tells you how fast; to tell Chapter 29 what to fix, diagnose the bottleneck's cause. §28.3's arithmetic-intensity estimate already flags the stencil as memory-bound ($I \approx 0.4$ flop/byte). Cachegrind confirms it and localises it. Run it on a small grid (cachegrind simulates, so it is slow):

$ gfortran -std=f2018 -Wall -O2 -g stencil_demo.f90 -o stencil_demo
$ valgrind --tool=cachegrind ./stencil_demo
$ cg_annotate cachegrind.out.<pid>

Two representative summaries — the solver's good loop order (inner over the first, contiguous index) versus a bad order (inner over the second index) on the identical computation:

                        D1 miss rate     LLd miss rate
  good order (i inner)      6.3%              0.9%
  bad  order (j inner)     18.8%              9.8%

Same D refs in both — identical work — but the bad order roughly triples the D1 misses and multiplies the last-level misses tenfold, each miss a stall of hundreds of cycles. cg_annotate pins those misses to the stencil's load line. This is the mechanism behind Chapter 27's "10× from loop order," now visible as a number. The diagnosis Chapter 29 receives is precise: a memory-bound stencil whose cost is cache misses, to be fixed by loop order and blocking — not by touching the arithmetic.

cachegrind vs perf. Cachegrind's miss rates are a reproducible guide to locality, but its absolute counts model a generic cache, not your CPU. For the real time cost on this machine, cross-check with perf stat -e cache-misses,cache-references ./solver, which reads the hardware counters. Use cachegrind for why and where; use perf for how much.

Phase 5 — The deliverable

You now hand Chapter 29 two things a single toc() could never provide: a reproducible baseline (min 0.098 s/run, median 0.100, range 0.004, with the flags and machine recorded) and a diagnosis (memory-bound stencil, bad-order cache signature). Every optimization in the next chapter will be judged against this baseline, re-run through this exact harness, and reported with the same three statistics. That is what makes a performance claim survive review — in Chapter 38's paper and anywhere else.

Discussion Questions

1. The harness resets u = u0 before every timed repetition. What exactly would go wrong — in both the timing and the physics — if you forgot the reset and let each run continue from the previous one's final state?
2. We report the minimum, the median, and the range. Construct a scenario where reporting only the minimum would mislead a reader, and one where reporting only the mean would.
3. The dead-code-elimination guard consumes u(2,2) into a sink. Why is consuming a single element enough to keep the whole time-stepping loop alive, and when might it not be enough?
4. Cachegrind is deterministic but models a generic cache; perf is exact but noisy. For the claim "loop reordering cut cache misses in half," which tool belongs in the paper, and why might you cite both?

Your Turn: Extensions

• Option A (build). Complete bench_solver into a compilable program: bundle the kinds, timers, heat_types, and heat_solver modules from project-checkpoint.f90, add the statistics of Phase 2, and print a reproducibility table like Phase 3's. Verify the maxval(u) = 100.00 invariant holds every run.
• Option B (design). Add a sweep mode that benchmarks the solver at grid sizes $128, 256, 512, 1024$ and prints per-step time versus grid size. Predict the scaling (per-step cost $\sim n^2$ interior points) and check whether your measured curve matches — and where cache effects bend it.
• Option C (diagnose). Deliberately write the stencil in the bad loop order, benchmark both orders through your harness, and run cachegrind on each. Report the wall-time ratio and the D1-miss-rate ratio, and confirm the slow one is slow for the reason cachegrind says.

Key Takeaways

• A defensible baseline is warm-up + repetition + a robust statistic + a reported spread + a recorded environment — not a single toc(). Three numbers (min, median, range) and a flags/machine table.
• Reset mutated state between repetitions, or you benchmark different work each time.
• Consume the result so the optimizer cannot delete the loop you are timing.
• Separate the two questions: the harness answers how fast (the baseline); cachegrind answers why (memory-bound, bad-order cache misses). Chapter 29 needs both to optimize the right thing the right way.