Case Study 1: Reading the Compiler's Mind

"The first step in optimization is not to make the code faster. It is to find out why it is slow."

Executive Summary

You have inherited a working, correct, and painfully slow Fortran routine: a 2D image/field smoother that replaces each interior cell with the average of its four neighbors — the nearest-neighbor loop at the heart of countless simulations and image filters. It gives the right answer and takes far too long. This case study is a diagnosis, not a rewrite. We will read the code the way the compiler reads it, find the three distinct reasons it cannot go fast — a loop order against the memory grain, a helper that blocks vectorization, and an aliasing hazard inherited from a C translation — and fix each while proving the output never changes. By the end you will be able to look at a slow array loop and predict, before touching a profiler, most of what is wrong with it.

Skills applied: column-major loop order and the cache line (§27.2); the no-aliasing advantage and how a C port can violate it (§27.3); pure/elemental as vectorization enablers (§27.4); reading the optimization report to confirm a diagnosis (§27.5); the honesty that timings are measured, not asserted.

Background

The inherited routine came, as these things do, from a well-meaning translation of a C reference implementation. Here it is, cleaned up just enough to compile with modern flags but otherwise faithful to what was handed over:

module smoother_slow
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
contains
  ! A helper that (unnecessarily) touches module state, and a smoother that
  ! loops in the wrong order and is called with the SAME array in and out.
  integer :: call_count = 0                     ! module variable -> impurity

  function avg4(a, i, j) result(m)
    real(dp), intent(in) :: a(:,:)
    integer,  intent(in) :: i, j
    real(dp) :: m
    call_count = call_count + 1                 ! side effect! blocks optimization
    m = 0.25_dp * (a(i-1,j) + a(i+1,j) + a(i,j-1) + a(i,j+1))
  end function avg4

  subroutine smooth(u, unew)
    real(dp), intent(in)  :: u(:,:)
    real(dp), intent(out) :: unew(:,:)
    integer :: i, j, nx, ny
    nx = size(u,1);  ny = size(u,2)
    unew = u
    do i = 2, nx-1                              ! outer over rows
      do j = 2, ny-1                            ! INNER over the LAST index -> against the grain
        unew(i,j) = avg4(u, i, j)
      end do
    end do
  end subroutine smooth
end module smoother_slow

It is called, in the driver, as call smooth(field, field) — the same array for input and output, exactly as the C original did with a single buffer. Three separate problems are hiding in plain sight. We take them in the order a performance engineer would.

Phase 1 — Characterize: Is It Even the Loop?

The discipline of Chapter 28 is measure before you touch — do not optimize a routine until a profiler has told you it is the bottleneck. Assume that step is done and smooth is genuinely where the time goes. Now ask the sharper question this chapter equips you to answer: why is it slow?

Begin with the memory access pattern, because for a nearest-neighbor loop the arithmetic is trivial (a few adds and a multiply per cell) and the cost is almost entirely memory traffic — the loop is memory-bound (§27.2). So the first thing to check is whether it walks memory with the grain. It does not. The inner loop runs over j, the last index. For a Fortran column-major array, consecutive inner iterations unew(i,2), unew(i,3), unew(i,4), … are a full column apart in memory — a different cache line almost every write, and the neighbor reads u(i,j-1), u(i,j+1) stride the same way. The routine fetches 64-byte cache lines and uses one value from each. Diagnosis 1: the loop nest is inside-out.

Phase 2 — Fix the Loop Order

The fix is mechanical and is the highest-return change available: swap the loops so the inner one runs over the first index.

do j = 2, ny-1                                  ! outer over columns
  do i = 2, nx-1                                ! INNER over the first index -> with the grain
    unew(i,j) = 0.25_dp * (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1))
  end do
end do

Now consecutive inner iterations unew(2,j), unew(3,j), unew(4,j) are adjacent in memory, each fetched cache line is used fully, and the u(i-1,j)/u(i+1,j) neighbors are one element away. Same answer, with the grain. On a grid larger than the last-level cache you would typically see this single change buy a several-fold speedup — the exact factor is yours to measure, and you would measure it the way Chapter 28 prescribes.

To be sure the fix is only a fix, verify the output is unchanged on a tiny grid you can check by hand. Take a 4×4 plate, hot top row u(1,:) = 100, everything else 0, one Jacobi sweep:

$ # one smoothing sweep, interior of a 4x4 hot-top plate
interior after one sweep (rows 2..3, cols 2..3):
   25.00   25.00
    0.00    0.00

The hand check: unew(2,2) = 0.25*(u(1,2)+u(3,2)+u(2,1)+u(2,3)) = 0.25*(100+0+0+0) = 25, and unew(2,3) likewise = 25; the row-3 interior reads only zeros, so it stays 0. Both loop orders produce this identical field, because each output cell reads only old neighbors — visitation order cannot change the result. The speed changed; the science did not.

Phase 3 — The Aliasing Hazard

The second defect is more dangerous than slowness, because it threatens correctness under optimization. The routine is called call smooth(field, field) — the same array as both the intent(in) argument u and the intent(out) argument unew. That aliases a written argument with a read one, which the Fortran standard forbids (§27.3): it is undefined behavior.

Why did it "work" in C with one buffer? Because the C original updated in place deliberately, accepting that each cell reads some already-updated neighbors — a Gauss–Seidel-flavored sweep, a different (though often still valid) algorithm. Translated literally to Fortran and then compiled at -O3, the compiler trusts the no-alias guarantee and may load blocks of the original u before any writes land, computing the Jacobi result instead. The two optimization levels can now disagree, and the standard says the program was never valid. This is Diagnosis 2, and it is exactly the bug from §27.3's "Find the Bug."

The fix is to make the intent explicit and honest: use two distinct buffers and ping-pong between them (call smooth(a, b) then call smooth(b, a)), which is what the corrected checkpoint does. If in-place Gauss–Seidel was genuinely intended, it must be written as a single intent(inout) argument, so the in-place reads are in the source where the compiler can see them — never smuggled in through aliasing.

The reasoning that matters: an aliasing bug is worse than a slow loop because a slow loop is visibly wrong and an aliasing bug is invisibly wrong — it can pass every test at -O0 and silently produce a different answer in the optimized production build. The no-aliasing advantage is a gift with a contract: the compiler optimizes freely because you promised distinctness, so breaking the promise breaks the program in the hardest-to-find way.

Phase 4 — The Helper That Blocks Vectorization

The third defect is the avg4 helper. It increments a module counter, call_count, on every call — a side effect. That single line makes avg4 impure, and an impure function called in the inner loop forbids vectorization: the compiler must call it, in order, once per cell, because each call might change observable state (§27.4). The loop cannot become packed SIMD instructions while a side-effecting call sits in its body.

Two fixes, both correct:

• Inline the arithmetic into the loop (as Phase 2 already did), removing the call entirely — simplest, and it is what you want for a body this small.
• Mark the helper pure (pure function avg4(...)) and delete the counter. The compiler enforces purity: leaving the call_count line in a pure function is a compile error, which is the language catching the exact impurity that was blocking you. Purity is not a hint; it is a checked promise.

Confirm the diagnosis rather than trusting it. Ask the compiler what it did:

$ gfortran -std=f2018 -O3 -march=native -fopt-info-vec smooth.f90 -o smooth
smooth.f90:NN:NN: optimized: loop vectorized using 32 byte vectors

Before the fixes, that line is absent (and -fopt-info-vec-missed names the reason — a call it could not analyze, or a dependence). After the fixes, the inner loop reports as vectorized. The report turns "I think it should be faster now" into "the compiler confirms the loop is vectorized." (The exact line/column and wording vary by gfortran version; read it for the yes/no, not the precise text.)

Phase 5 — Put It Together and Re-measure

The corrected smoother is short, with the grain, alias-free, and vectorizable:

pure subroutine smooth(u, unew)                 ! distinct in/out; never call smooth(x, x)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp), intent(in)  :: u(:,:)
  real(dp), intent(out) :: unew(:,:)
  integer :: i, j, nx, ny
  nx = size(u,1);  ny = size(u,2)
  unew = u
  do j = 2, ny-1
    do i = 2, nx-1
      unew(i,j) = 0.25_dp * (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1))
    end do
  end do
end subroutine smooth

Three defects, three fixes, and a table that is the whole case study in miniature:

The final speedup is a number you measure, not one we assert — and you would measure it with the methodology of Chapter 28, then push further (blocking, do concurrent) with Chapter 29. But the diagnosis — the reading of the compiler's mind — is done, and it took no profiler at all: just the three ideas of this chapter, applied to the code as written.

Discussion Questions

1. Of the three defects, which is the most dangerous to leave in a production scientific code, and why? (Consider: which one can pass all your tests and still be wrong?)
2. The C original updated in place on purpose. Was the Fortran translation's smooth(field, field) a bug in the translation or a bug in the original intent? How would you find out, and how would you write the in-place version legally in Fortran?
3. The avg4 helper's side effect was a counter — arguably harmless. Construct an argument for why the compiler is nonetheless right to refuse to vectorize a loop containing it. What could a side-effecting call legitimately do that would make reordering unsafe?
4. You fixed the loop order and the report now says "loop vectorized." Does that guarantee the routine is now fast? What else could still bottleneck it?

Your Turn: Extensions

• Option A. Take the original slow module and, without changing the algorithm's result, apply the three fixes one at a time, confirming after each with -fopt-info-vec and the tiny 4×4 hand check. Keep a log of what each fix changed in the report.
• Option B. Write the legal in-place (Gauss–Seidel) version with a single intent(inout) argument. Verify by hand on the 4×4 plate that it gives a different interior than the Jacobi version (because it reads updated neighbors), and explain why that is a different algorithm, not a bug.
• Option C. Instrument both the slow and fixed versions with system_clock on a large grid, run each a few times, and report your own measured factor with a one-paragraph note on the conditions (grid size, compiler, flags) — the reproducibility habit of Chapter 37.

Key Takeaways

• Most "slow array loop" diagnoses are one of three things, and you can spot all three by reading the code: wrong loop order (memory), aliasing (correctness under optimization), and a side-effecting call in the hot loop (vectorization).
• The aliasing defect is the one to fear, because it is invisible: correct at -O0, wrong at -O3, and never diagnosed. The no-aliasing advantage is a contract — keep written arguments distinct.
• The optimization report converts belief into evidence. Fix, then ask the compiler whether it worked, then measure — never optimize by faith.