Chapter 29 — Key Takeaways (Optimization Techniques)

A one-page reference for loop reordering, fusion/fission, cache blocking, SIMD/auto-vectorization, do concurrent, contiguous, and the judgment to stop — the chapter that turns measurement into speed.

The optimization ladder (largest payoff first)

Loop reordering — the one rule with no exceptions

$$ \text{inner loop over the FIRST index} \;\Longrightarrow\; \text{unit stride, full cache lines} $$

• Fortran is column-major: a(i,j) and a(i+1,j) are adjacent in memory; a(i,j) and a(i,j+1) are a whole column apart.
• Inner loop over i (first index) = fast; over j (second index) = several × slower on a large array — for the identical result.
• NumPy is row-major (opposite): its fast inner axis is the last one. Mind this at the f2py boundary.

do j = 1, n
  do i = 1, n          ! i inner: consecutive iterations touch consecutive memory
    a(i,j) = ...
  end do
end do

Fusion, fission, hoisting

• Hoist divides (many × a multiply's cost); doing rx = alpha*dt/dx**2 once beats $n^2$ divides per step.

Cache blocking (tiling) — and when NOT to

• Blocking works a small tile that stays cache-resident through all its reuses, then advances.
• Helps compute-bound, high-reuse kernels (matrix multiply: $O(n^3)$ work / $O(n^2)$ data → $O(n)$ reuse).
• Does almost nothing for a memory-bound single stencil sweep (low reuse). Loop order + fusion are its whole story.
• Decide by arithmetic intensity (flops per byte) and the roofline:

SIMD and auto-vectorization

• SIMD = one instruction, many operands: AVX register = 4 real(dp); AVX-512 = 8.
• Auto-vectorization = compiler emits SIMD from a plain loop (no intrinsics), at -O2+.
• A loop vectorizes when it is: independent (no loop-carried dependency), unit-stride, countable, non-aliasing, simple body (no un-inlinable calls).

y(i) = a*x(i) + y(i)      ! vectorizes (independent)
y(i) = y(i-1) + x(i)      ! does NOT (loop-carried dependency)

• Confirm it: gfortran -O3 -fopt-info-vec (what vectorized) / -fopt-info-vec-missed (what didn't + why).

do concurrent (performance use)

do concurrent (j = 2:ny-1, i = 2:nx-1)     ! "these iterations are independent"
  u_new(i,j) = ...                          ! reads OLD u, writes NEW u_new
end do

• A promise the compiler trusts, not a parallel directive. It permits vectorization/parallelization; it does not guarantee parallel execution (needs -ftree-parallelize-loops, or nvfortran for GPU).
• Break the promise (in-place update, a running total) → undefined result. Keep read/write arrays distinct.
• Your on-ramp to Chapter 33 (OpenMP) and Chapter 35 (GPU).

contiguous and aliasing (§29.4)

• real(dp), intent(in), contiguous :: u(:,:) promises unit stride — the compiler skips a contiguity check / copy-in-out and may vectorize. Most valuable on pointers (which may point at strided slices).
• Fortran assumes procedure arguments do not alias (the Ch. 27 advantage) — that is why the stencil vectorizes. Pointers/targets can alias and forfeit it; prefer allocatable, keep read/write arrays distinct, mark pointers contiguous.

Verifying an optimization (correctness is paramount)

Compile flags introduced

Pitfalls

• Inner loop over the second index — strided, column-major-hostile; several × slow. The #1 numerical-Fortran bug.
• A wasted pass / temporary — a full array written and re-read is pure traffic for a memory-bound kernel.
• do concurrent with a hidden dependency — undefined result; the compiler trusts, does not check.
• Plain pointer in a hot loop — assumed non-contiguous and possibly aliasing → no vectorization; add contiguous.
• Blocking a memory-bound sweep / hand-writing intrinsics — effort for near-zero gain; wrong rung of the ladder.
• Optimizing before measuring — the profiled 3% is where effort belongs; the other 97% stays readable.

Numbers / rules worth memorizing

• Inner loop over the first index. The single rule with no exceptions.
• Most stencils are memory-bound: optimize by moving less data; you cannot out-compute the memory ceiling.
• A tuned dgemm beats a hand-rolled matrix multiply by ~an order of magnitude — call the library.
• Optimize in payoff order; the first hour buys most of the speedup (Knuth's "critical 3%").

The heat-solver piece added

The stencil step is now optimized: correct loop order, one fused sweep (no lap temporary, divisions hoisted), two buffers swapped via contiguous pointers (O(1), no copy), expressed as do concurrent — fast, bit-for-bit identical to Chapter 24, and ready to go parallel behind the same interface. Next, Chapter 30 picks the flags that finish the job.