Chapter 28: Profiling and Benchmarking — Finding the Bottleneck Before Optimizing

"Rule 2. Measure. Don't tune for speed until you've measured, and even then don't unless one part of the code overwhelms the rest." — Rob Pike, Notes on Programming in C (1989)

Overview

In Chapter 27 you learned why Fortran is fast — first-class arrays, column-major layout, the no-aliasing advantage — and you saw, in principle, how a change of loop order could turn a routine ten times faster. "In principle" is the operative phrase. This chapter is where principle becomes measurement. Before you touch a single line to make it faster, you are going to learn to answer three questions with numbers instead of guesses: How long does my code actually take? Where, exactly, does that time go? And is the bottleneck the arithmetic or the memory?

That order — measure, locate, classify, and only then optimize — is not a suggestion. It is the professional discipline that separates engineers who make code fast from hobbyists who make code different. Every experienced performance programmer has the same scar: a week spent lovingly hand-tuning a routine that, it turned out, accounted for two percent of the runtime, while the real culprit sat untouched a few functions away. The human intuition for where a program spends its time is famously, reliably wrong. Rob Pike's second rule of programming, in the epigraph above, is the whole chapter in one line: measure first. A profiler is not a luxury you reach for when a code is mysteriously slow; it is the instrument you reach for before you decide anything is slow at all.

This is a performance chapter, so we will be precise about a promise made throughout the book. Every piece of timing code here is exact, standard Fortran that compiles with gfortran -std=f2018 -Wall, and you should type it and run it. But every timing number — every "0.8 seconds," every profiler percentage, every cache-miss rate — is illustrative. We cannot run code as we write, and even if we could, the numbers would be ours, not yours. Performance is a property of a specific machine on a specific day, so treat all figures here as representative shapes, clearly labelled, and trust only what your own measurements tell you.

In this chapter, you will learn to:

• Time a region of code correctly with system_clock and cpu_time, and know precisely what each measures.
• Tell wall-clock time from CPU time, and read the gap between them as a clue about I/O and parallelism.
• Run gprof on a real program and read its flat profile and call graph to find the hot loop.
• Decide whether a kernel is memory-bound or compute-bound, and pick the right kind of fix.
• Use valgrind --tool=cachegrind to see the cache misses behind a slow loop order.
• Build a benchmark you can defend — warmed up, repeated, statistically honest, reproducibly configured.

Learning Paths

How to read this chapter by track. - 🔬 Scientist ("my code is slow and I don't know why") — this is your chapter. Read §28.1 and §28.2 closely; they are the two tools you will use every week. §28.5's methodology will save you from fooling yourself. - 📖 Standard — read straight through; profiling is the bridge from "why Fortran is fast" (Ch. 27) to "how to make your Fortran fast" (Ch. 29). - 🔧 Legacy — profiling is exactly how you find the load-bearing 3% of an inherited 100,000-line code before you dare change anything. §28.2 and the 80/20 rule (§28.5) are your orientation. - ⚡ HPC — you know the tools; skim §28.1. Spend your time on §28.3 (memory- vs compute-bound) and §28.4 (cachegrind), the analysis that tells you whether to chase bandwidth or flops.

28.1 Measure First: cpu_time, system_clock, and Wall vs CPU Time

Optimization without measurement is superstition. So the first tool in the box is the humblest: a way to ask the program how long something took. Fortran gives you two standard intrinsics for this, and the difference between them is not pedantry — it is the difference between two genuinely different quantities.

Definition (profiling). Profiling is measuring where a running program spends its resources — time, memory traffic, cache misses — so you can find the small part of the code that dominates the cost. A profiler answers "where does the time go?"; a timer (this section) answers the simpler "how long did this take?". You almost always start with a timer to confirm there is a problem, then reach for a profiler to locate it.

cpu_time — processor seconds consumed

The intrinsic subroutine cpu_time returns the amount of processor time your program has used so far, in seconds, as a real number:

real(dp) :: t0, t1
call cpu_time(t0)
! ... work to be timed ...
call cpu_time(t1)
print '(a, f8.3, a)', 'CPU time: ', t1 - t0, ' s'

Two rules govern it. First, only differences are meaningful — the absolute value of t0 is an arbitrary origin, so you always subtract two readings. Second, cpu_time measures time the CPU spent on your program, not time elapsed on the clock on the wall. If your program sleeps, waits for a disk read, or blocks on the network, that idle time is not charged to CPU time. And on most implementations, if your program runs on several cores at once, the CPU time is summed across them — so a program that runs for 2 seconds of real time on four busy cores may report roughly 8 seconds of CPU time. Hold that fact; in Chapter 33 it becomes the way you prove your OpenMP parallelism is actually using the cores you asked for.

system_clock — the wall-clock idiom

For elapsed real time — the seconds a stopwatch on the wall would measure — you use system_clock. This is the timing idiom you will write hundreds of times, so learn it exactly once, correctly.

Definition (system_clock). The intrinsic subroutine system_clock(count, count_rate, count_max) reports an integer tick count from a processor clock. count is the current tick; count_rate is the number of ticks per second; count_max is the largest value count reaches before wrapping back to zero. Elapsed wall-clock time in seconds is real(count_end - count_start) / real(count_rate). All three arguments are intent(out) integers (you may request any subset by keyword), and you must read count_rate from the call — never assume its value.

Here is the canonical form, and the one used everywhere in this book. Note the deliberate use of 64-bit integers:

use, intrinsic :: iso_fortran_env, only: dp => real64, int64
integer(int64) :: c_start, c_end, c_rate
real(dp)       :: wall

call system_clock(count_rate=c_rate)     ! ticks per second — read it, don't guess
call system_clock(count=c_start)
! ... work to be timed ...
call system_clock(count=c_end)
wall = real(c_end - c_start, dp) / real(c_rate, dp)
print '(a, f10.6, a)', 'wall time: ', wall, ' s'

Work through the arithmetic once by hand so the idiom is never mysterious. Suppose count_rate comes back as $1{,}000{,}000$ (microsecond ticks), and your two counts are $3{,}200{,}000$ and $3{,}700{,}000$. Then the elapsed time is $(3{,}700{,}000 - 3{,}200{,}000)/1{,}000{,}000 = 0.5$ seconds. That is all system_clock ever does: subtract two integers and divide by the rate. The only ways to get it wrong are to forget the division, to divide by a guessed rate, or to mishandle the wrap-around — which we turn to now.

⚠️ Common Pitfall — the two integer-kind traps. Why 64-bit integers, and why read the rate? The resolution and wrap point of the clock depend on the integer kind of the arguments. Request the clock through 64-bit (int64) integers and gfortran gives you a high-resolution monotonic clock — typically nanosecond ticks, so count_rate on the order of $10^9$, with a count_max so large the counter would take on the order of centuries to wrap. Request it through default 32-bit integers and you may get a coarse millisecond clock that also wraps in a matter of weeks or less — and if it wraps mid-measurement, count_end < count_start and your elapsed time goes negative. Two habits make this a non-issue: always use int64 counters, and always read count_rate from the call rather than hard-coding a number. If you must support a clock that can wrap, guard it explicitly:

fortran integer(int64) :: c_start, c_end, c_max, ticks call system_clock(count_rate=c_rate, count_max=c_max) ! ... take c_start, do work, take c_end ... ticks = c_end - c_start if (ticks < 0_int64) ticks = ticks + c_max + 1_int64 ! the clock wrapped once wall = real(ticks, dp) / real(c_rate, dp)

With int64 the guard essentially never fires (that is why we use int64), but writing it down is how you show you understand what the counter is.

Wall time vs CPU time: reading the gap

The two intrinsics measure two quantities, and the relationship between them is diagnostic. On a quiet machine running a single-threaded, compute-only program, they roughly agree. When they disagree, the disagreement tells you something:

⚡ Performance Note: This little table is a free diagnostic you get before any profiler. Time a slow program both ways. If wall $\gg$ CPU, no amount of loop tuning will help — your program is waiting, not computing, and the fix is in the I/O (Part VI) or the scheduling, not the arithmetic. If wall $\approx$ CPU and both are large, now you have a computation to profile. Measuring both costs four lines and rules out a whole category of wasted effort.

A reusable timers module

Writing the six-line idiom inline every time is error-prone, so we wrap it once — and this is the first brick of this chapter's Project Checkpoint, the canonical timers module the solver will use from here on. It offers the tic()/toc() pair that MATLAB and NumPy users will recognise:

module timers
  use, intrinsic :: iso_fortran_env, only: dp => real64, int64
  implicit none
  private
  public :: tic, toc

  integer(int64) :: start_count = 0_int64      ! module state: persists between calls (saved)
  integer(int64) :: rate        = 1_int64
  logical        :: have_rate   = .false.

contains

  subroutine tic()                             ! start (or restart) the stopwatch
    if (.not. have_rate) then
      call system_clock(count_rate=rate)       ! read the rate once
      have_rate = .true.
    end if
    call system_clock(count=start_count)
  end subroutine tic

  function toc() result(seconds)               ! elapsed wall seconds since the last tic()
    real(dp) :: seconds
    integer(int64) :: now
    call system_clock(count=now)
    seconds = real(now - start_count, dp) / real(rate, dp)
  end function toc

end module timers

The module-level variables carry state between calls — the one legitimate use of the saved-variable behaviour we warned about as a trap in Chapter 6. Here it is deliberate and encapsulated: tic stashes the start count, toc reads the clock again and returns the difference in seconds. Note that neither procedure can be pure — they read a global clock and touch module state, which are exactly the side effects pure forbids.

🔄 Check Your Understanding. 1. Your system_clock reading returns count_rate = 1000 and elapsed counts of 47. How many seconds is that, and what does the coarse rate warn you about your measurement? 2. A program reports wall time 12 s and CPU time 3 s. Is it compute-bound? What is it doing?

Answers

(1) $47/1000 = 0.047$ s. A rate of only 1000 ticks/second means millisecond resolution — timing anything shorter than a few milliseconds is mostly quantisation noise; either time a longer region or use an int64 clock. (2) No — wall $\gg$ CPU means the CPU was idle three-quarters of the time. It is I/O-bound (or waiting on something), so optimizing arithmetic is pointless; look at the reads/writes.

28.2 Profiling with gprof

A timer tells you the whole took 8 seconds. It cannot tell you that 6.2 of those seconds were inside one innermost loop. For that you need a profiler, and the classic, portable, no-install one for compiled code is gprof, which ships with the GNU toolchain alongside gfortran.

gprof works in two complementary ways at once. It instruments your code — the compiler inserts a tiny bookkeeping call at the entry of every procedure, so gprof can count how many times each was called and who called whom. And it samples — while the program runs, a timer interrupts it a hundred or so times a second and records which procedure the program counter is in, building up a statistical picture of where the time goes. Call counts are exact; times are statistical, which is why a profiled run must be long enough to gather many samples before its percentages mean anything.

The workflow

Three steps: compile-and-link with -pg, run the program (which drops a file named gmon.out in the working directory), then run gprof on the pair.

$ gfortran -std=f2018 -Wall -O2 -pg -g heat_profiled.f90 -o heat_profiled
$ ./heat_profiled            # runs normally, and writes gmon.out on exit
$ gprof ./heat_profiled gmon.out > profile.txt

The -pg flag must appear on both the compile and the link step. We add -g for symbol names and, crucially, keep the optimization flag we intend to ship with (-O2 here) — profiling a slow debug build is one of the great time-wasters, addressed in the pitfall below.

Reading a flat profile

The flat profile is the first thing gprof prints, and often the only thing you need. It ranks procedures by how much time was spent inside each one, ignoring their callees. Here is a representative flat profile for an instrumented run of the heat solver — constructed for illustration; your numbers will differ:

Flat profile:

Each sample counts as 0.01 seconds.
  %   cumulative   self              self     total
 time   seconds   seconds    calls  ms/call  ms/call  name
 78.0      3.90     3.90     2000     1.95     1.95  laplacian_
 17.0      4.75     0.85     2000     0.43     2.38  step_
  5.0      5.00     0.25        1   250.00  5000.00  MAIN__

Read it column by column, because every column answers a question:

• % time — the fraction of total runtime spent inside this procedure. laplacian_ is 78% of the program. This column alone usually tells you what to optimize.
• self seconds — the raw seconds spent inside the procedure itself, not its callees. The flat profile is sorted by this.
• cumulative seconds — a running total of self down the list; the last row equals the whole runtime.
• calls — how many times it was invoked (exact, from instrumentation). laplacian_ ran 2000 times, once per step.
• self ms/call and total ms/call — milliseconds per call, excluding and including callees. step_'s total per-call (2.38 ms) is much larger than its self (0.43 ms) — because most of each step call is spent down inside laplacian.

Notice the trailing underscores: gfortran decorates external procedure names, so laplacian appears as laplacian_ and the main program as MAIN__. That is normal; it is how you know you are reading Fortran symbols. The verdict here is unambiguous and typical of numerical code: nearly four-fifths of the time is in one function. That function is where the entire optimization budget should go, and everything else is noise.

Reading a call graph

Sometimes a function is hot not because it is slow but because it is called from an unexpected place, or a million times. The call graph shows the caller/callee structure with time attributed along the edges. A representative excerpt — again illustrative:

index  % time    self  children    called     name
                                                 <spontaneous>
[1]    100.0     0.25      4.75        1         MAIN__ [1]
                 0.85      3.90     2000/2000    step_ [2]
-----------------------------------------------
                 0.85      3.90     2000/2000    MAIN__ [1]
[2]     95.0     0.85      3.90     2000         step_ [2]
                 3.90      0.00     2000/2000    laplacian_ [3]
-----------------------------------------------
                 3.90      0.00     2000/2000    laplacian_ [3]
[2]     ...
[3]     78.0     3.90      0.00     2000         laplacian_ [3]

Each block describes one procedure (the line flush with its [index]). The lines above that line are its callers; the lines below are its callees. The self/children columns split its time between the procedure itself and everything it calls. Reading the middle block: step_ was called 2000 times, all from MAIN__; it spent 0.85 s in itself and 3.90 s in its children; and its single child is laplacian_, which accounts for all 3.90 s. The call graph confirms the flat profile's story and adds the why: step is expensive only because it calls laplacian, so the hot loop lives inside laplacian.

⚠️ Common Pitfall — never profile the -O0 build. The single most common profiling mistake is compiling with -pg but without your optimization flags — profiling the unoptimized debug build. The results lie in two directions at once. A -O0 build leaves in function calls that -O2 would inline away, so cheap helper functions show up as hot when in the real build they vanish entirely. And the relative balance between routines shifts, because the optimizer speeds up some far more than others. Rule: profile the build you will ship. Compile -O2 -pg -g, not -O0 -pg. The one twist is that heavy inlining can merge small functions into their callers, so a hot inlined kernel may appear as time in the caller — a reason to keep the kernel a pure function during profiling and read the call graph, not just the flat list.

🐍 Python Comparison: Python's profiler is cProfile, and its output is the direct analogue of gprof's flat profile — functions ranked by cumulative and total time. The workflow is the same (python -m cProfile myscript.py), the discipline is identical (measure, find the hot spot), and the lesson transfers completely. What does not transfer is the scale: because the pure-Python hot loop you find is often 50–100× slower than the compiled equivalent, the fix frequently is not "tune the Python" but "move this loop into Fortran and call it from Python" — the f2py workflow of Chapter 15. Profiling is where you decide which loop earns that treatment.

gprof is the portable classic, but it is not the only game. On Linux, perf (perf stat, perf record) reads the CPU's hardware performance counters directly — cache misses, branch mispredictions, instructions per cycle — with almost no overhead and no recompile; gprofng, Intel VTune, and Linaro MAP go further still. We stay with gprof here because it is everywhere and its concepts — flat profile, call graph, self vs total — are the vocabulary all the others share.

28.3 Identifying Hot Loops; Memory-Bound vs Compute-Bound

The flat profile pointed at laplacian. But a function is not yet a target; the target is the loop inside it. Zooming from function to loop, and then asking what is holding that loop back, is the analytical heart of this chapter.

Definition (hot loop). A hot loop — or hot spot — is the small region of code, very often a single innermost loop, where a program spends the large majority of its time. Numerical codes are extraordinarily concentrated this way: it is common for one stencil sweep, one matrix kernel, or one particle-interaction loop to own 80–99% of the runtime. Find it, and you have found essentially all of the available speedup; optimize anywhere else and you are polishing brass on a sinking ship.

Inside laplacian, the hot loop is the sweep over the interior grid points — the five-point stencil applied to every cell. Once you have located it, the decisive question is not "how do I make this loop shorter?" but "what is this loop waiting on?" There are two answers, and they demand opposite fixes.

Definition (compute-bound vs memory-bound). A loop is compute-bound when its speed is limited by how fast the processor can perform arithmetic — the floating-point units are the bottleneck, and data arrives faster than it can be crunched. A loop is memory-bound when its speed is limited by how fast data can be moved from memory into the processor — the arithmetic units sit idle, starved, waiting for the next cache line. The same hardware runs both; which one you are depends on the ratio of arithmetic to data movement in your loop.

That ratio has a name.

Definition (arithmetic intensity). Arithmetic intensity is the number of floating-point operations a loop performs per byte of data it moves from memory, $I = \text{flops} / \text{bytes}$. Low intensity (few flops per byte) means the loop is starved for data and is memory-bound; high intensity (many flops per byte) means the loop keeps the arithmetic units busy and is compute-bound. It is the single most useful number for predicting which class a kernel falls into.

Apply it to the stencil. Each interior cell reads its own value and four neighbours — but on a good sweep, neighbours are shared between adjacent cells and already sitting in cache, so the fresh traffic is roughly one new real(dp) (8 bytes) loaded and one stored per cell: on the order of 16 bytes. The arithmetic is a handful of additions and a couple of multiplies — call it ~6 flops. That is an arithmetic intensity of roughly $6/16 \approx 0.4$ flops per byte: very low. The stencil is emphatically memory-bound. The processor could do the six flops in a heartbeat; it spends its life waiting for the next row of the grid to arrive from memory.

🚪 Threshold Concept — in scientific computing, the bottleneck is usually the memory, not the math. Newcomers optimize the arithmetic: fewer multiplies, a clever algebraic simplification, a lookup table. For the stencil — and for most of the array-sweeping kernels at the core of simulation — this is almost useless, because the arithmetic was never the bottleneck. The processor is waiting on data. The wins come from moving data less and reusing it more: the right loop order so you stream through memory contiguously (Ch. 27), cache blocking so data is reused while it is still in cache (Ch. 29), and shrinking the data itself. Once you internalise that a modern CPU can do dozens of flops in the time it takes to fetch one number from main memory, the entire craft of optimization reorganises around feeding the processor rather than unburdening it.

The contrast makes it concrete. A naïve dense matrix multiply, $C = AB$, does $O(n^3)$ flops on $O(n^2)$ data, so each element loaded is reused $O(n)$ times — high arithmetic intensity, and it can be made strongly compute-bound. That is exactly why a tuned BLAS dgemm (Chapter 21) blocks the computation to keep sub-tiles resident in cache and then saturates the floating-point units: it has arranged for the math, not the memory, to be the limit. Your hand-rolled triple loop moves the same data far more times and stays memory-bound — which is the deep reason, promised back in Chapter 21, that you cannot beat the library by trying harder at the arithmetic.

How do you tell which class you are in, without a hardware lab? Three tests, in increasing rigour:

1. The back-of-the-envelope intensity, as above. Count fresh bytes and flops per iteration. Below ~1 flop/byte, assume memory-bound until proven otherwise.
2. The "add arithmetic" experiment. Add some cheap, non-eliminable arithmetic to the loop body (say, another multiply-add on data already in registers) and re-time. If the loop barely slows down, the FP units had spare capacity — you were memory-bound. If it slows in proportion, you were compute-bound.
3. The cache-miss measurement — §28.4. A high last-level-cache miss rate is the direct fingerprint of a memory-bound loop.

🔗 Connection: The classification is not academic — it routes you to the right chapter. Memory-bound (the stencil, most array kernels)? Your fixes are loop order and cache blocking, in Chapter 29. Compute-bound (dense linear algebra, tight transcendental-function loops)? Your fixes are SIMD vectorization and better algorithms — or, best of all, calling a tuned library. The "roofline" picture that unifies both — a ceiling set by memory bandwidth rising to a plateau set by peak flops — you will meet in the capstone's performance analysis, Chapter 38.

🔄 Check Your Understanding. 1. A loop reads two real(dp) arrays and writes one, doing a single multiply-add per element (this is the classic y = y + a*x, "axpy"). Estimate its arithmetic intensity and classify it. 2. You add an extra x = x*1.0000001_dp inside a hot loop and the runtime is unchanged. What have you learned?

Answers

(1) Per element: ~24 bytes moved (two loaded, one stored, 8 bytes each) for 2 flops (one multiply, one add) → $I \approx 2/24 \approx 0.08$ flop/byte. Deeply memory-bound — this is why BLAS Level-1 operations like axpy are bandwidth-limited and cannot be sped up by cleverer arithmetic. (2) The floating-point units had idle capacity that the extra flop soaked up "for free," so the loop is memory-bound — it is waiting on data, not on the arithmetic.

28.4 Cache Behaviour with valgrind --tool=cachegrind

The intensity estimate says the stencil is memory-bound. cachegrind lets you see it — to watch a bad loop order generate a storm of cache misses, and a good one avoid them. It is the measurement that makes Chapter 27's "loop order can be 10×" claim tangible.

You met valgrind in Chapter 13 as the tool that catches memory errors — leaks, invalid reads, uninitialised values. Valgrind is really a family of tools sharing one engine, and cachegrind is the one that models the cache. Instead of running your instructions on the real CPU, cachegrind runs them on a simulated machine with a modelled cache hierarchy (a first-level instruction cache I1, a first-level data cache D1, and a last-level cache LL) and counts every memory reference and every miss. Because it is a simulation, its counts are deterministic and reproducible — the same run gives the same numbers every time, unlike wall-clock timing — at the price of running your program perhaps 20–100× slower. So you run it on a small problem.

$ gfortran -std=f2018 -Wall -O2 -g stencil_demo.f90 -o stencil_demo
$ valgrind --tool=cachegrind ./stencil_demo
$ cg_annotate cachegrind.out.12345          # 12345 = the PID it printed

Cachegrind prints a summary to the terminal and writes a detailed cachegrind.out.<pid> file that cg_annotate turns into a per-line, per-function breakdown. Here is a representative summary for the good loop order (inner loop over the first array index, striding contiguously through column-major memory) — illustrative numbers; yours will differ:

==12345== D   refs:      420,000,000  (300,000,000 rd + 120,000,000 wr)
==12345== D1  misses:     26,300,000  ( 25,000,000 rd +   1,300,000 wr)
==12345== LLd misses:      3,900,000  (  3,600,000 rd +     300,000 wr)
==12345== D1  miss rate:         6.3% (        8.3%   +         1.1%  )
==12345== LLd miss rate:         0.9% (        1.2%   +         0.3%  )

Now the bad loop order (inner loop over the second index, striding across rows — jumping a full column's worth of memory every iteration), on the identical computation:

==12345== D   refs:      420,000,000  (300,000,000 rd + 120,000,000 wr)
==12345== D1  misses:     78,900,000  ( 76,000,000 rd +   2,900,000 wr)
==12345== LLd misses:     41,000,000  ( 40,000,000 rd +   1,000,000 wr)
==12345== D1  miss rate:        18.8% (       25.3%   +         2.4%  )
==12345== LLd miss rate:         9.8% (       13.3%   +         0.8%  )

Same number of data references — the arithmetic and the array are identical — but the bad order roughly triples the D1 misses and multiplies the last-level-cache misses tenfold. Every one of those extra LL misses is a stall of hundreds of cycles while the CPU waits for main memory. That is the entire mechanism behind the "10× from loop order" headline: not more work, just work that constantly falls out of cache.

cg_annotate localises it to the exact source line. A representative annotated excerpt shows the per-line data reads (Dr) and D1 read misses (D1mr) against the stencil's load:

--------------------------------------------------------------------------------
         Dr        D1mr          line
--------------------------------------------------------------------------------
100,000,000  19,000,000          lap(i,j) = ( u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) &   ! bad order
          0           0                     - 4.0_dp*u(i,j) ) / h2

Nineteen million of a hundred million reads on that one line missed the first-level cache — a smoking gun. Fix the loop order and the same line's D1mr collapses. This is the loop you will actually reorder in Chapter 29; cachegrind is how you prove, before and after, that the reorder did what you claimed.

⚠️ Common Pitfall — cachegrind models a cache, not your cache. Because cachegrind simulates a generic cache hierarchy (it auto-detects sizes but idealises the behaviour), its miss rates are a faithful guide to your algorithm's locality but its absolute counts are not a promise about your specific CPU. Use it to compare two versions of your code (the ratio is rock-solid and reproducible), not to predict wall-clock time. For the real hardware truth — actual cycles lost to actual misses on your processor — use perf stat -e cache-misses,cache-references ./prog, which reads the on-chip counters. The two are complementary: cachegrind tells you why and where with perfect reproducibility; perf tells you how much it cost on this machine.

🔄 Check Your Understanding. Two versions of a kernel show identical D refs but version B has a D1 miss rate of 22% versus version A's 6%. Which version has the better loop order for a Fortran array, and why are the reference counts identical?

Answer

Version A. Identical D refs means both do the same arithmetic on the same data — the amount of work is unchanged. The only difference is locality: A's lower miss rate means it accesses memory contiguously (inner loop over the first, fastest-varying index — column-major order, Ch. 5/27), so each cache line is fully used before eviction. B strides across memory and evicts lines before reusing them. Same work, different memory order, wildly different speed.

28.5 Benchmarking Methodology: Warm-up, Repetition, Variance, and the 80/20 Rule

You can now time a region, profile a program, find the hot loop, and classify it. The last skill is the one that keeps you honest: producing a benchmark number you would stake a claim on. A single toc() reading is not a measurement — it is an anecdote.

Definition (benchmark methodology). Benchmark methodology is the disciplined practice of measuring performance so the numbers are meaningful and reproducible: warm up before timing, repeat the measurement many times, report a robust statistic (minimum or median, not a lone sample) together with its spread, control the environment (release build, quiet machine, fixed CPU frequency, same input), and change one thing at a time. Skip these and your "benchmark" measures the operating system's mood, not your code.

Consider the forces that corrupt a naïve single measurement:

• Cold caches and cold pages. The very first pass over a fresh array pays for cache misses and first-touch page faults that later passes do not. Timing the first iteration measures start-up, not steady state.
• CPU frequency scaling. Modern processors idle at a low clock and ramp up ("turbo") under load, taking milliseconds to reach full speed — so an early measurement catches the CPU half-asleep. Conversely, under sustained load they may throttle down for heat.
• Interference. Another process, a background indexer, an email client — anything that steals the core for a moment inflates a measurement.

Definition (warm-up). A warm-up is one or more untimed executions of the code before you start timing, run to pull the data into cache, fault in the memory pages, and bring the CPU to its steady operating frequency — so the timed runs measure the code's steady-state performance, not its start-up transient.

The pattern that neutralises all of this is warm up, then repeat, then take a robust statistic:

use timers, only: tic, toc
real(dp) :: best, t
integer  :: rep
integer, parameter :: n_warm = 3, n_rep = 20

do rep = 1, n_warm            ! warm-up: run but do not time
  call kernel(a, b, c)
end do

best = huge(1.0_dp)           ! we will keep the fastest (see below)
do rep = 1, n_rep
  call tic()
  call kernel(a, b, c)
  t = toc()
  best = min(best, t)         ! minimum = the least-disturbed run
end do
print '(a, es12.4, a)', 'best of 20: ', best, ' s'

Which statistic? For "how fast can this code go," the minimum is often the right choice: every disturbance (a stolen core, a cache eviction) can only make a run slower, so the fastest run is the one least contaminated by noise, the closest to the true cost. When you care about typical behaviour under real conditions, the median is more honest, because it is robust to outliers in a way the mean is not — a single 10× hiccup barely moves the median but can wreck the mean. Report the spread too (the range, or the median with an interquartile range); a benchmark without a notion of its own variance is a point with no error bar. What you should almost never report is a single run, or a mean dominated by one outlier.

⚠️ Common Pitfall — the compiler deleted your benchmark. The nastiest benchmarking bug is a fast result that is a lie. If you compute a value in your timed loop and never use it, an optimizing compiler is entitled to notice the result is dead and delete the entire computation — leaving you timing an empty loop and celebrating an infinite speedup. The cure is to use the result in a way the compiler cannot predict away: accumulate it and print it after the timer stops.

fortran real(dp) :: sink sink = 0.0_dp do rep = 1, n_rep call tic() call kernel(a, b, c) t = toc() sink = sink + c(1) ! consume a result so the loop can't be optimized away end do print *, best, sink ! printing sink keeps the whole computation alive

The same trap catches the warm-up loop and, subtly, any benchmark that recomputes an identical result each iteration (the compiler may hoist it out and run it once). Make each iteration's work observably matter.

The 80/20 rule — why this all pays off

Underneath the whole chapter sits an empirical regularity strong enough to plan around.

Definition (the 80/20 rule). The 80/20 rule (the Pareto principle applied to performance): the large majority of a program's runtime is spent in a small minority of its code — very often roughly 80% of the time in 20% of the code, and in numerical software the split is frequently far more extreme, 95/5 or 99/1. The practical consequence: profile to find that small hot fraction, pour your effort there, and deliberately ignore the rest.

Donald Knuth put the discipline most memorably, and his phrasing even carries the numbers:

📜 From History. "We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil. Yet we should not pass up our opportunities in that critical 3%." — Donald Knuth, "Structured Programming with go to Statements" (1974). The famous half-sentence is almost always quoted alone, but the whole thought is the better lesson: most code should be written for clarity and left alone; the profiler's job is to find the 3% that is worth making ugly for speed. Optimize before you profile and you will spend your 3% of effort on the wrong 3% of code.

The rule also sets a ceiling you should compute before optimizing. If the profiler says a routine is 20% of your runtime, then making it infinitely fast — removing it entirely — speeds the whole program by at most 20% (a fivefold-shorter part leaves the other 80% untouched). That ceiling on the payoff from speeding one part is formalised as Amdahl's Law, which governs parallel speedup and which you will meet properly in Chapter 31; for now, the sequential version of its lesson is enough — always check the size of the prize before paying for it. Optimizing the 78%-of-runtime laplacian can, at best, roughly quintuple the whole solver; optimizing the 5% MAIN__ can, at best, speed it up by 5%. The profile told you which prize is worth chasing.

🔄 Check Your Understanding. A profile shows function assemble at 8% of runtime. Your colleague spent a week making it 4× faster. What is the most the whole program could have sped up, and what should they have profiled for first?

Answer

Making 8% of the runtime 4× faster removes three-quarters of that 8%, i.e. 6 percentage points — the whole program gets at most ~6% faster (new time $\approx 0.92 + 0.08/4 = 0.94$, a 1.06× speedup). A week for 6%. They should have optimized the largest entry in the flat profile first; the 80/20 rule and the Amdahl ceiling both say the prize was elsewhere.

Project Checkpoint

This chapter's increment: instrument the solver and find its hot loop. In Chapter 24 the solver became real; in Chapter 27 you reasoned about why loop order matters. Now you measure it — the checkpoint that turns the solver from "works" into "works, and here is exactly where its time goes."

Two additions. First, adopt the timers module from §28.1 as the project's canonical timers.f90, and use it to time the time-stepping loop:

use timers, only: tic, toc
real(dp) :: elapsed
integer  :: n

call tic()
do n = 1, n_steps
  call step(field, alpha, dt)          ! the frozen Ch. 6 signature, unchanged
end do
elapsed = toc()
print '(a, i0, a, f10.4, a)', 'ran ', n_steps, ' steps in ', elapsed, ' s'
print '(a, es12.4, a)', 'per step: ', elapsed / real(n_steps, dp), ' s'

The per-step time is the number you will track for the rest of Part VII — the single figure of merit the optimization chapters will drive down. (Its value is machine-dependent, so the harness above prints it rather than asserting it; the compilable version in code/project-checkpoint.f90 also runs a tiny hand-checkable $5\times5$ case whose output is exact, so you can confirm the instrumented solver still computes the right physics — $28, 32, 28$ across the second row after two steps, precisely as in Chapter 24.)

Second, profile it. Compile the driver with -O2 -pg -g, run it, and read the flat profile. You will find — as the representative profile in §28.2 showed — that the sweep inside laplacian dominates: the hot loop is the five-point stencil. That is not a surprise; it is a confirmation, and confirming your intuition with a measurement (rather than trusting it) is the entire point. You now have, in hand, the one loop worth optimizing and the evidence that it is worth optimizing.

That evidence is precisely what Chapter 29 consumes: it reorders and blocks this loop, and re-runs this harness to prove the speedup. And the before/after numbers you record here become the "serial baseline" of the capstone's performance analysis in Chapter 38. Measure now; optimize next; and never again optimize a loop you have not first proven is hot.

Summary

This chapter replaced guessing with measuring: time it, profile it, classify it, and only then optimize it.

The system_clock idiom, memorised: read count_rate from the call, use int64 counters, subtract two counts, divide by the rate. Never guess the rate; never assume the clock cannot wrap.

The methodology checklist: warm up, repeat ($\ge$ ~20), take the minimum (best case) or median (typical), report the spread, build with -O2/-O3, quiet the machine, change one thing at a time, consume the result so the compiler can't delete the benchmark.

The two things to remember: first, the profiler almost always contradicts your intuition about where the time goes — so measure first, optimize second, never the reverse. Second, the 80/20 rule means the prize lives in one hot loop; compute the size of the prize (its share of the runtime) before you spend a day chasing it.

Spaced Review

Retrieval practice on two earlier chapters this one leans on. Answer before expanding.

1. (Ch. 6.) The solver's laplacian is a pure function. What two things is a pure procedure forbidden to do, and why does that both help the optimizer and mean tic/toc cannot be pure?
   

   Answer

   A `pure` procedure has no side effects: it may not modify anything outside itself (no changing global/module state, no `intent(out)`/`inout` on non-result data beyond its arguments' declared intents) and may not perform I/O. That guarantee lets the compiler reorder, cache, or parallelize calls freely. `tic`/`toc` both read the system clock and `tic` writes module state — external side effects — so they are inherently impure.

2. (Ch. 6.) State the intent of each argument of the frozen step(field, alpha, dt), and why field must be intent(inout).
   

   Answer

   `field` is `intent(inout)` — it is read (the old temperatures) and written (the new ones), updated in place. `alpha` and `dt` are `intent(in)` — read-only parameters. The compiler rejects any attempt to write `alpha`/`dt`, catching a class of bugs before the code runs.

3. (Ch. 27.) In cache terms, why does looping with the first array index innermost beat looping with the second index innermost, on a Fortran 2D array?
   

   Answer

   Fortran is column-major: consecutive first-index elements (`u(i,j)`, `u(i+1,j)`) are adjacent in memory. Making the first index innermost walks memory contiguously, so each cache line loaded is fully used before eviction (few misses). Second-index-innermost strides by a whole column each step, touching one element per cache line and evicting lines before reuse — the high-miss pattern §28.4 measures.

4. (Ch. 27.) What is the no-aliasing advantage, and how does it let a Fortran compiler produce faster code than an equivalent C loop?
   

   Answer

   Fortran assumes distinct procedure arguments do not overlap in memory (do not alias). So when a routine reads one array and writes another, the compiler may safely reorder, vectorize, and keep values in registers, knowing a write can't secretly change a value it's about to read. C must assume any two pointers might alias (unless told `restrict`), forcing conservative reloads. Same algorithm, more compiler freedom, faster code.

5. (Ch. 27.) You read a compiler optimization report (-fopt-info) and it says the stencil loop was not vectorized. How does that finding combine with a profiler and cachegrind to guide your next move?
   

   Answer

   The three are complementary. The profiler says the stencil is the hot loop (worth effort); cachegrind says it is memory-bound (a high miss rate); the opt-report says it did not vectorize (a missed compute optimization). Together they say: the ceiling is memory bandwidth, so fix locality first (loop order/blocking, Ch. 29); vectorization will help only once the data is arriving fast enough to feed it. Measurement, not the report alone, sets the priority.

What's Next

You have the bottleneck in your sights: the stencil loop inside laplacian, proven hot by the profiler, proven memory-bound by its arithmetic intensity and its cache-miss rate. Chapter 29 is where you finally get to make it fast — reordering the loop for column-major memory, blocking it so data is reused while still in cache, and coaxing the compiler to vectorize it — and then re-running the very harness you built here to prove each change earned its keep. You measured; now you optimize. In that order, always.