Case Study 1: Chasing the Wrong Loop

"Bottlenecks occur in surprising places, so don't try to second guess and put in a speed hack until you've proven that's where the bottleneck is." — Rob Pike, Rule 1 of programming

Executive Summary

A research team proudly announces that they have "optimized the particle simulator." They spent a week rewriting the time integrator — the routine that advances each particle's position and velocity — replacing a clean loop with a hand-unrolled, hand-vectorized version they are visibly pleased with. The wall-clock time of the whole simulation improved by about three percent. They are baffled: the integrator itself is now three times faster, so where did the speedup go?

It went nowhere, because the integrator was never the problem. This case study is a guided post-mortem. You will confirm the physics is correct on a hand-checkable case, fix the team's broken homemade timer, profile the code the way they should have from the start, and read the flat profile and call graph to discover that a single $O(n^2)$ pairwise loop — untouched by their week of work — owns ninety percent of the runtime. Then you will compute, in advance, the largest speedup their integrator work could ever have delivered, and see that three percent was roughly the ceiling. The lesson is the chapter's thesis: measure first, and optimize the loop the profiler points at — not the one you find most interesting.

Skills applied

• Timing a region correctly with system_clock, and diagnosing a broken timer (§28.1).
• Distinguishing wall-clock from CPU time to rule out I/O (§28.1).
• Running gprof and reading a flat profile and a call graph to locate the hot loop (§28.2).
• Classifying the hot loop and recognising an $O(n^2)$ algorithm as the real lever (§28.3).
• Computing the size of the prize — the speedup ceiling — before spending effort (§28.5, the 80/20 rule).

Background

The simulator integrates $n$ gravitating particles. Each step does two things: compute the pairwise forces (every particle feels every other — an $O(n^2)$ double loop), then advance each particle (an $O(n)$ single loop). Here is the shape of the two routines, simplified to their cost structure. The expensive one:

pure function total_pe(x, y, z) result(pe)      ! the O(n^2) pairwise loop — the real hot spot
  use, intrinsic :: iso_fortran_env, only: dp => real64
  real(dp), intent(in) :: x(:), y(:), z(:)
  real(dp) :: pe, r
  integer  :: i, j, n
  n  = size(x)
  pe = 0.0_dp
  do i = 1, n - 1
    do j = i + 1, n                              ! every pair once -> n(n-1)/2 iterations
      r  = sqrt((x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2)
      pe = pe + 1.0_dp / r
    end do
  end do
end function total_pe

The cheap one — the integrator the team lovingly rewrote — is a single pass over the particles, a few multiply-adds each: $O(n)$ work, a rounding error next to the $O(n^2)$ force loop for any interesting $n$. That asymmetry is the whole story, and it was visible in the algorithm before anyone ran anything.

Phase 1 — Confirm the physics on a hand-checkable case

Never profile a wrong program; you will optimize a bug. First, pin the computation to a case you can check by hand. Three particles at $(0,0,0)$, $(3,0,0)$, and $(0,4,0)$ form a 3–4–5 triangle, so the three pairwise distances are exactly $3$, $4$, and $5$, and the total inverse-distance potential is

$$ \text{PE} = \frac{1}{3} + \frac{1}{4} + \frac{1}{5} = \frac{20 + 15 + 12}{60} = \frac{47}{60} = 0.78333\ldots $$

real(dp) :: x(3) = [0.0_dp, 3.0_dp, 0.0_dp]
real(dp) :: y(3) = [0.0_dp, 0.0_dp, 4.0_dp]
real(dp) :: z(3) = [0.0_dp, 0.0_dp, 0.0_dp]
print '(a, f10.6)', 'total PE = ', total_pe(x, y, z)

total PE =   0.783333

The physics is right. Now, and only now, is it worth asking where the time goes.

Phase 2 — Fix the broken timer

The team "measured" with a homemade timer that "kept giving weird, sometimes negative numbers," so they stopped trusting it and switched to eyeballing a stopwatch. Here is what they wrote:

integer :: c0, c1                       ! BUG 1: default integer -> low-res, wraps in seconds/minutes
real    :: secs
call system_clock(c0)
call run_simulation()
call system_clock(c1)
secs = (c1 - c0) / 1000.0               ! BUG 2: count_rate GUESSED as 1000, not read
print *, 'seconds:', secs               ! BUG 3: never divided by the *actual* rate

Three defects, all from §28.1. The default-integer clock is coarse and wraps quickly, so a long run can end with c1 < c0 and a negative elapsed time — the "weird numbers." The rate is hard-coded to 1000 instead of read from the call, so even when it doesn't wrap the seconds are wrong by whatever factor the true rate differs. The fix is the canonical idiom:

use, intrinsic :: iso_fortran_env, only: dp => real64, int64
integer(int64) :: c0, c1, rate
real(dp)       :: secs
call system_clock(count_rate=rate)      ! READ the rate
call system_clock(count=c0)
call run_simulation()
call system_clock(count=c1)
secs = real(c1 - c0, dp) / real(rate, dp)

Port it. The team's data scientist timed the Python prototype with timeit, which auto-repeats and reports the best of several runs. The Fortran equivalent is exactly the warm-up-and-repeat harness of §28.5 wrapped around this idiom — same discipline, four more lines. A correct timer is the price of admission; it does not yet tell you where the time goes, which is the next phase.

Sanity check. With the fixed timer, also print cpu_time. Wall and CPU come back nearly equal — the run is compute-only, not waiting on disk — so we have a genuine computation to profile, not an I/O stall.

Phase 3 — Profile it, and read the evidence

Compile the release build with profiling and run gprof:

$ gfortran -std=f2018 -Wall -O2 -pg -g nbody.f90 -o nbody
$ ./nbody && gprof ./nbody gmon.out > profile.txt

The flat profile (representative — your numbers will differ) is unambiguous:

  %   cumulative   self              self     total
 time   seconds   seconds    calls  ms/call  ms/call  name
 90.2      7.22     7.22     1000     7.22     7.22  total_pe_
  4.1      7.55     0.33     1000     0.33     0.33  integrate_
  3.5      7.83     0.28     1000     0.28     0.28  apply_bc_
  2.2      8.01     0.18        1   180.00  8010.00  MAIN__

The call graph confirms total_pe_ is a leaf called once per step from MAIN__, with no hidden children inflating it. The verdict: total_pe is 90% of the run; integrate is 4%. The team spent a week on the 4% line.

Phase 4 — Compute the size of the prize before optimizing

This is the step that would have saved the week. If a routine is a fraction $p$ of the runtime and you speed it by a factor $K$, the whole-program speedup is at most

$$ S = \frac{1}{(1 - p) + p/K}. $$

Plug in the integrator's $p = 0.041$. Even making it infinitely fast ($K \to \infty$) gives $S = 1/(1 - 0.041) = 1.043$ — a 4.3% ceiling. Their actual $3\times$ ($K = 3$) gives $S = 1/(0.959 + 0.041/3) = 1/0.973 = 1.028$, about 2.8% — exactly the "three percent" they observed. The measurement was never mysterious; it was the arithmetic of the 80/20 rule.

Now do the same for the real prize. total_pe is $p = 0.902$. A mere $2\times$ on it yields $S = 1/(0.098 + 0.902/2) = 1/0.549 = 1.82$ — an 82% whole-program speedup, from half the effort. And because the loop is $O(n^2)$, the largest lever is not micro-optimization at all but algorithm: a tree-based $O(n\log n)$ method (Barnes–Hut) changes the exponent, not the constant. The profiler pointed at the loop; the loop's complexity pointed at the biggest win of all.

Phase 5 — The right move

The disciplined sequence, in order: the physics is verified (Phase 1); the timer is trustworthy (Phase 2); the profiler names total_pe as 90% of the cost (Phase 3); the prize arithmetic says that is where 80%+ of the achievable speedup lives (Phase 4). Only now do you optimize — and you optimize total_pe, not the integrator. How to make that specific memory-and-compute loop faster — better memory access, vectorization, and the algorithmic change — is the subject of Chapter 29. The contribution of this chapter is everything that comes before the optimizing: knowing, with evidence, that total_pe is the loop to chase and the integrator was a beautiful waste of a week.

Discussion Questions

1. The team's integrator rewrite was genuinely $3\times$ faster and genuinely correct. In what sense was it still a mistake? What should have happened before the week of work began?
2. The $O(n^2)$ pairwise loop is the hot spot at every $n$ large enough to matter. Argue why algorithmic improvement (changing the exponent) is a different and often larger lever than the micro-optimizations of Chapter 29 — and when micro-optimization is nonetheless the right choice.
3. The broken timer sometimes printed negative times. Explain the exact mechanism, and why using int64 counters would have prevented the symptom even without fixing the guessed rate.
4. Suppose the flat profile had instead shown four routines at ~22% each. How would that change your optimization strategy, and what does it say about the "80/20" assumption?

Your Turn: Extensions

• Option A (analyse). Take the representative flat profile above and compute, for each of the four routines, the maximum whole-program speedup obtainable by making that one infinitely fast. Rank them. Which single routine should a second week of effort target, and what is its ceiling?
• Option B (instrument). Add a system_clock-based timer around each of the three phase routines separately (force, integrate, boundary) so the code prints its own per-phase breakdown every 100 steps — a poor man's profiler that needs no -pg. When is this inline instrumentation preferable to gprof?
• Option C (port). Reproduce the team's mistake in miniature: write both loops, use the correct timer, and confirm that speeding the $O(n)$ integrator barely moves the total while speeding the $O(n^2)$ loop moves it a lot. Report the two whole-program speedups and compare them to the Phase 4 ceilings.

Key Takeaways

• Intuition about where time goes is reliably wrong. The team optimized the routine they found interesting; the profiler named a different one. Measure first, every time.
• A broken timer is worse than no timer — it produces confident nonsense. Read count_rate, use int64 counters, subtract, divide.
• Compute the size of the prize before paying for it. A routine that is 4% of the runtime caps the whole program's speedup at ~4%, no matter how clever the rewrite; the flat profile plus the formula $S = 1/((1-p)+p/K)$ tells you this in ten seconds.
• The profiler finds the hot loop; the loop's complexity finds the biggest win. For the $O(n^2)$ pairwise force, the largest lever is algorithmic, not a speed hack — a distinction Chapter 29 builds on.