Case Study 1: The Loop That Strode the Wrong Way

"The real problem is that programmers have spent far too much time worrying about efficiency in the wrong places and at the wrong times." — Donald E. Knuth

Executive Summary

A colleague's image-smoothing routine — a Jacobi averaging sweep, structurally identical to your heat solver's stencil — is "unbearably slow" on the large images it was written for, and they are convinced they need to rewrite it in C. They do not. This case study is a guided diagnosis and repair of existing code: you will profile it to confirm where the time goes (Chapter 28), recognize from a cache profile that it is memory-bound and access-hostile, find the two problems — a loop nest that strides across memory the wrong way, and a two-pass structure that hauls a whole temporary array through the memory system for nothing — and fix both without changing a single output pixel. Then, crucially, you will stop, because a roofline estimate proves there is nothing left to win from this kernel on this hardware. The whole repair is two edits; the skill is knowing which two.

Skills applied

• Profiling to localize the hot loop and classify it memory- vs compute-bound (§29.2, and Chapter 28).
• Recognizing a column-major loop-order defect and correcting it (§29.1).
• Applying loop fusion to eliminate a temporary array and its memory traffic (§29.1).
• Verifying that an optimization is bit-for-bit correctness-preserving (§29.3, the checkpoint's rule).
• Using the roofline / arithmetic-intensity argument to decide when to stop (§29.2, §29.5).

Background

The routine smooths a 2D field by replacing each interior cell with the average of its four neighbors — one Jacobi sweep — and it is written like this:

subroutine smooth_slow(a, b, n)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  integer,  intent(in)  :: n
  real(dp), intent(in)  :: a(n,n)
  real(dp), intent(out) :: b(n,n)
  real(dp) :: nbr(n,n)                     ! a whole temporary array
  integer  :: i, j
  b = a
  do i = 2, n-1                            ! <-- outer over i
    do j = 2, n-1                          ! <-- inner over j  (STRIDES columns!)
      nbr(i,j) = a(i-1,j) + a(i+1,j) + a(i,j-1) + a(i,j+1)
    end do
  end do
  do i = 2, n-1                            ! second pass over the whole interior
    do j = 2, n-1
      b(i,j) = 0.25_dp * nbr(i,j)
    end do
  end do
end subroutine smooth_slow

The arithmetic is right — but two performance sins hide in plain sight, and on a $4000\times4000$ image they cost dearly.

Phase 1 — Measure first: where does the time go?

Before changing anything, obey Chapter 28's first law and measure. A gprof flat profile of a run that smooths a large image many times is unambiguous:

So the entire optimization target is smooth_slow; nothing else is worth a minute of your time (that is the 80/20 rule doing its job). Next, classify the loop. Run it under valgrind --tool=cachegrind and the last-level-cache miss rate is high — far higher than the handful of misses the arithmetic would need. The kernel does about four adds and one multiply per cell, over an array far larger than cache: it is memory-bound, and it is missing cache far more than a well-ordered sweep should. That single fact — hot, memory-bound, cache-hostile — points straight at the access pattern.

Sanity check. Before hunting further, confirm the kernel is memory-bound by estimate, not just by profile. It does ~5 flops per interior cell and moves ~5 real(dp) reads plus a write — call it ~48 bytes — per cell. That is roughly $0.1$ flop/byte: deeply memory-bound (§29.2). A memory-bound kernel's speed is set by how it moves data, so the fix will be about access order and traffic, never about the arithmetic.

Phase 2 — Diagnose and fix the loop order

Look at the inner loop: it runs over j, the second index. Fortran stores a(n,n) column-major, so a(i,j) and a(i,j+1) are a whole column — n elements — apart in memory. Each inner iteration therefore jumps n elements forward, using one value from each fetched cache line and discarding the other seven. That is the high miss rate cachegrind saw. The fix is to swap the loops so the inner one runs over i, the first index, walking each column contiguously:

do j = 2, n-1                              ! outer over j
  do i = 2, n-1                            ! inner over i  -> unit stride, cache-friendly
    nbr(i,j) = a(i-1,j) + a(i+1,j) + a(i,j-1) + a(i,j+1)
  end do
end do

Not one computed value changes — the assignment is the same for every (i,j); only the order of visiting them does. On a large image this reorder alone typically recovers a large factor (illustratively several ×; measure yours), because it converts a stream of cache misses into full cache-line reuse.

⚠️ Common Pitfall — the wrong order is invisible on a small test. Had the colleague benchmarked on a $200\times200$ image that fits in cache, both loop orders would have run at the same speed and hidden the bug entirely. The defect only bites when the array overflows cache — which is exactly the size the routine was written for. Measure at the real problem size, or the benchmark lies to you.

Phase 3 — Diagnose and fix the wasted pass

The reorder fixed how the code walks memory; the two-pass structure is about how much. The routine builds a whole nbr(:,:) array in the first sweep and reads it back in the second — an extra full write and full read of an $n\times n$ array, pure traffic for a memory-bound kernel. Fuse the two sweeps into one, and the temporary vanishes:

subroutine smooth_fast(a, b, n)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  integer,  intent(in)  :: n
  real(dp), intent(in)  :: a(n,n)
  real(dp), intent(out) :: b(n,n)
  integer  :: i, j
  b = a                                    ! carry the (unchanged) edges
  do j = 2, n-1
    do i = 2, n-1                          ! fused: average AND scale in one sweep
      b(i,j) = 0.25_dp * (a(i-1,j) + a(i+1,j) + a(i,j-1) + a(i,j+1))
    end do
  end do
end subroutine smooth_fast

One sweep, no nbr array, unit stride, and the scaling folded in. For a memory-bound kernel the speedup is roughly the traffic saved (§29.2): removing a full write and read of the temporary is a further solid gain on top of the reorder.

Now prove it is still correct. Both versions applied to a $5\times5$ field with a single hot cell a(3,3)=8 must produce the identical plus-sign of smoothed values:

program verify_smooth
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter :: n = 5
  real(dp) :: a(n,n), bs(n,n), bf(n,n)
  integer  :: i
  a = 0.0_dp;  a(3,3) = 8.0_dp
  call smooth_slow(a, bs, n)
  call smooth_fast(a, bf, n)
  print '(a)', 'smoothed field (fast):'
  do i = 1, n
    print '(5f6.1)', bf(i,:)
  end do
  print '(a, f6.1)', 'max |slow - fast| = ', maxval(abs(bs - bf))
  ! (smooth_slow / smooth_fast as above, in the same file)
end program verify_smooth

$ gfortran -std=f2018 -Wall -O3 verify_smooth.f90 -o vs && ./vs
smoothed field (fast):
   0.0   0.0   0.0   0.0   0.0
   0.0   0.0   2.0   0.0   0.0
   0.0   2.0   0.0   2.0   0.0
   0.0   0.0   2.0   0.0   0.0
   0.0   0.0   0.0   0.0   0.0
max |slow - fast| =    0.0

Hand-check the hot cell's neighbors: each of b(2,3), b(3,2), b(3,4), b(4,3) averages the single value $8$ with three zeros, $0.25\times8 = 2$; the center b(3,3) averages its four (still-zero) neighbors to $0$. The heat has spread into a plus sign, and max |slow - fast| is exactly $0$ — because both kernels compute the identical expression for every cell, the optimization is bit-for-bit correctness-preserving, not merely close.

Phase 4 — Know when to stop

The colleague, delighted, now wants to cache-block the sweep and hand-write AVX intrinsics. Talk them out of it. The roofline argument from §29.2 settles it: this kernel is memory-bound at ~$0.1$ flop/byte, and after the reorder and fusion it now streams the field essentially the minimum number of times, in order. It is sitting at the memory ceiling. Cache blocking a single sweep captures reuse that a five-point average does not have; hand-vectorizing arithmetic that is already idling behind memory latency changes nothing. There is no more speed to extract from this kernel on this hardware — the next honest step, if it is still too slow, is Part VIII's parallelism (more memory bandwidth across more cores), not more single-core cleverness.

⚡ Performance Note. The most valuable optimization skill on display here is not the two fixes — it is the decision to stop after them. Two edits, both improving readability (a fused loop is shorter than a two-pass one), captured essentially all the available speed; a week of intrinsics would have added complexity and near-zero performance. That is Knuth's epigraph in action: the effort was spent in the right place (the profiled 97.8%) at the right time (after measuring), and then it stopped.

Phase 5 — Prevent the recurrence

Fixing this routine is easy; teaching the codebase not to regrow the bug is the professional move. Two cheap habits:

• A loop-order lint rule. In review, flag any nested array loop whose inner index is not the first. It is the single most common performance defect in numerical Fortran, and it is mechanical to spot.
• A size-honest benchmark. Keep a benchmark that runs at the production image size, not a toy that fits in cache, so a future reorder regression shows up as a real slowdown instead of hiding.

Discussion Questions

1. The two fixes each preserved the output exactly, yet one (loop order) is invisible in the source's meaning while the other (fusion) actually shortens it. Why is it a happy accident here that the fastest version is also the most readable — and when is that not the case?
2. The kernel was memory-bound at ~$0.1$ flop/byte. Construct a plausible variant of a smoothing routine that would be compute-bound instead, and say how your optimization strategy would change.
3. Cachegrind showed the miss rate but not the cause. What is the danger of "optimizing to the profiler" — chasing a metric like cache misses — versus reasoning about the access pattern? When could lowering the miss rate fail to lower the runtime?
4. The colleague wanted to rewrite in C. Given Fortran's no-aliasing advantage (Chapter 27), argue why the C rewrite would most likely have been slower, not faster, unless they added restrict everywhere.

Your Turn: Extensions

• Option A (analyze). Instrument both smooth_slow and smooth_fast with the tic/toc timer from the Chapter 28 timers module, run each on a $2000\times2000$ field, and report the measured speedup of (i) the reorder and (ii) the fusion, separately. Do your numbers land in the illustrative "several ×, then a solid further slice" shape? (They are your numbers — the point is to get them.)
• Option B (profile). Reproduce the diagnosis: build smooth_slow, confirm with gprof that it is the hot routine, and with cachegrind that the slow order has a markedly higher miss rate than the fast order. Seeing the miss rate fall when you swap the loops is the whole lesson made concrete.
• Option C (extend). Add a -fopt-info-vec build of smooth_fast and read the report: did the inner loop vectorize? Then make the buffers pointer (without contiguous) and watch the report change to "possible aliasing." Restore vectorization with the contiguous attribute (§29.4).

Key Takeaways

• Measure before you touch anything. The profile said 97.8% of time was in one routine; everything else was noise. Optimize the profiled hot loop, nothing else.
• The wrong loop order is the most common performance bug in numerical Fortran — inner loop over the second index strides across column-major memory. Swap it; the answer is unchanged, the speed is transformed.
• Fuse away temporaries for a memory-bound kernel: a whole extra array written and read is pure traffic. The fused loop is usually shorter as well as faster.
• Verify bit-for-bit. A dependency-free reorder or fusion changes which order the same arithmetic runs, not which arithmetic — so max |before - after| must be exactly $0$.
• Stop at the roofline. Once a memory-bound kernel streams its data the minimum number of times, in order, it is at the memory ceiling; blocking and hand-vectorization buy nothing. Knowing to stop is the skill.