Chapter 29: Optimization Techniques — Loop Optimization, SIMD, and Cache-Friendly Code

"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 E. Knuth, Structured Programming with go to Statements (1974)

Overview

You have a profiler's verdict in hand. Chapter 28 taught you to measure before you touch anything — to find the one loop where your program actually spends its time, and to know whether that loop is starved for data or starved for arithmetic. Now, and only now, do you get to optimize it. This is the chapter that makes code fast on purpose.

The temptation, once you have permission to optimize, is to reach for the clever tricks first — hand-written SIMD intrinsics, four-deep manual loop unrolling, a bit-twiddling reciprocal. Resist it. The largest, most reliable speedups in numerical Fortran come from a short list of unglamorous transformations that give the compiler what it needs: walk memory in the right order, touch each byte as few times as possible, and promise the compiler the independence and contiguity it cannot prove on its own. Do those, compile with optimization on, and you will capture most of the available speed. The exotic techniques exist, and we will meet them, but they belong to Knuth's "critical 3%" — reached for rarely, only when the profiler proves you need them, and always paid for in readability.

This chapter is also, deliberately, an honest one about limits. Some optimizations that sound powerful — cache blocking a single stencil sweep, hand-vectorizing a memory-bound loop — buy almost nothing, and a good engineer knows why before wasting a day on them. And the chapter ends where every performance story honestly ends: with the recognition that for a whole class of problems, the fastest thing you can do is stop writing loops and call a library that a team of specialists spent decades tuning. We measured that library sitting under NumPy back in Chapter 1; here we admit it beats us, and learn when to let it.

In this chapter, you will learn to:

• Reorder, fuse, and fission loops so they respect Fortran's column-major memory and move as little data as possible.
• Apply cache blocking (tiling) to a kernel with data reuse — and recognize the kernels where it does nothing.
• Structure code so the compiler auto-vectorizes it into SIMD instructions, and use do concurrent to say "these iterations are independent" out loud.
• Use the contiguous attribute and understand Fortran's aliasing guarantees well enough to help the optimizer — and to know when, rarely, you must tune by hand.
• Judge when further optimization is not worth it, and why a tuned BLAS still beats your best loop.

Learning Paths

How to read this chapter by track. - 🔬 Scientist ("my solver is too slow") — §29.1 (loop order and fusion) and §29.5 (when to stop, and call BLAS) are the highest-value pages you will read all book. Skim §29.4's aliasing theory. - 📖 Standard — read straight through; this is the payoff of Part VII, turning the why of Chapter 27 and the measure of Chapter 28 into do this. - 🔧 Legacy — the Jacobi and matrix-multiply kernels here are exactly the shapes buried in the old codes of Part IV; §29.1 is how you speed them up without rewriting the physics. - ⚡ HPC — this is your chapter. Read every section closely, especially §29.3 (do concurrent is your bridge to Chapter 33) and §29.2's honest accounting of when tiling pays.

⚠️ A standing note on the numbers in this chapter. Every speedup figure here is illustrative — an order of magnitude typical of the transformation, not a measurement from a specific run. We have not executed any of this code; the correctness values (what the programs print) are hand-computed and exact, but the timings depend on your CPU, its caches, your compiler, and its flags. The one rule that outranks every technique in this chapter is Chapter 28's: measure it yourself. When we write "roughly 5× faster," read it as "expect a large factor here, and go find out what it is on your machine."

29.1 Loops: Reordering, Fusion, and Fission

Almost all the time in a numerical program is spent in a handful of loops, and almost all of their time is spent waiting for data to arrive from memory. So the first and most productive thing you can do to a hot loop is not reduce its arithmetic — it is arrange for the data it needs to already be close by when it asks. Three loop transformations do exactly that, and none of them changes a single computed value.

Loop reordering for column-major

You met this idea as theory in Chapter 5 and as a measured effect in Chapter 27; here it becomes a habit. Fortran stores a two-dimensional array column by column — a(1,1), a(2,1), a(3,1), … march down the first column before the second begins — so two elements that are adjacent in memory differ in their first index. A CPU never fetches one number from memory; it fetches a whole cache line (typically 64 bytes, eight real(dp) values) at a time, betting that you will want the neighbors next. Walk memory in order and you use all eight; jump around and you use one and throw the line away.

The consequence is a rule with no exceptions in dense Fortran: the innermost loop should run over the first array index. Consider filling or sweeping a matrix:

! FAST — inner loop over the first index i: consecutive iterations touch
!        consecutive memory (down a column). The cache line is fully used.
do j = 1, n
  do i = 1, n
    a(i,j) = f(i,j)
  end do
end do

! SLOW — inner loop over the second index j: each iteration jumps a whole
!        column (n elements) forward in memory. On a large array this can be
!        several times slower, for an identical result.
do i = 1, n
  do j = 1, n
    a(i,j) = f(i,j)
  end do
end do

The two nests compute exactly the same array. On a small grid you will never notice the difference; on a large one the slow order can cost you a factor of several, because it turns one cache line's worth of useful work into eight cache misses.

⚡ Performance Note: the size of the penalty scales with how badly you overrun the cache. When n is small enough that a whole column fits in the L1 or L2 cache, the wrong order is nearly free — the strided data is still resident. As n grows past the cache, every jump becomes a fresh trip to main memory, and the gap widens to the illustrative "several ×" you will read about. This is why you measure at the problem's real size (Chapter 28): a benchmark on a grid that fits in cache will lie to you about the grid that does not.

🐍 Python Comparison: NumPy has the identical hazard wearing the opposite default. NumPy arrays are row-major (C order), so its fast inner axis is the last one — a[i, :] is contiguous, a[:, j] is strided. When you move an array across the Chapter 15 boundary with f2py, this is the whole ballgame: pass a C-ordered NumPy array to Fortran expecting column-major and either a silent transpose or an expensive copy happens. The lesson is the same in both languages — know which index is contiguous, and loop over it innermost — but the answer is mirrored.

Loop fusion

Definition (loop fusion). Loop fusion (also called loop jamming) merges two adjacent loops that run over the same index range into a single loop, so that each data element is touched once instead of once per loop. The payoff is memory traffic: a value produced in the first loop can be consumed in the second while it is still in a register or cache line, rather than written out to memory and read back. Fusion also removes the loop-overhead of a second sweep and, often, an entire temporary array.

Your own solver has a textbook opportunity for this. The Chapter 24 step computes the whole Laplacian into a temporary array and then sweeps the field again to apply it:

lap = laplacian(field%u, field%dx, field%dy)          ! sweep 1: write the whole lap(:,:)
field%u(2:nx-1,2:ny-1) = field%u(2:nx-1,2:ny-1)  &    ! sweep 2: read u and lap, write u
                       + alpha*dt * lap(2:nx-1,2:ny-1)

That is two passes over an $n \times n$ field and a full temporary array lap that is written once and read once — pure overhead. Fuse the two passes into one loop nest, computing each point's Laplacian and applying it in the same breath, and the temporary array vanishes:

do j = 2, ny-1
  do i = 2, nx-1                                       ! one sweep: compute AND apply
    u_new(i,j) = u(i,j) + rx*(u(i-1,j) - 2.0_dp*u(i,j) + u(i+1,j))  &
                        + ry*(u(i,j-1) - 2.0_dp*u(i,j) + u(i,j+1))
  end do
end do

Two clean wins hide in that rewrite. First, fusion: one pass over memory instead of two, and no lap array allocated, filled, and re-read. Second — a smaller but real classic — the divisions 1/dx**2 and 1/dy**2 have been hoisted out of the loop into the constants rx = alpha*dt/dx**2 and ry = alpha*dt/dy**2, computed once. A floating-point divide costs many times what a multiply does; doing it once instead of $n^2$ times per step is free speed. (A good compiler at -O2 may hoist a loop-invariant divide for you, but writing it explicitly costs nothing and guarantees it.)

💡 Intuition: think of memory bandwidth as a narrow pipe and your arrays as water you must pull through it. Every optimization in this section is really the same move — pull less water through the pipe. Loop fusion pulls the field through once instead of twice. Correct loop order pulls it through in full buckets instead of teaspoons. Neither does less arithmetic; both do less waiting, and waiting is where the time goes.

Loop fission

Fusion's mirror image is occasionally what you want instead.

Definition (loop fission). Loop fission (also loop distribution) splits a single loop into two or more loops over the same range, each doing part of the original body. You fission a loop to isolate a troublesome part — most often to separate a portion the compiler can vectorize from a portion it cannot (a call to a non-inlinable function, an if with a data-dependent branch), so the clean part can run at full SIMD speed. Fission can also relieve register pressure or reduce the number of distinct memory streams a single loop juggles, when too many arrays in one body cause the hardware to thrash.

Fusion and fission are not rivals; they are tools for opposite symptoms. Fuse when two loops share data and you are paying to move it twice. Fission when one loop is doing two jobs and the messy job is holding the tidy one back. The compiler applies both automatically in many cases, but it is conservative — it will not fuse loops if it cannot prove the merge is safe, and your knowledge of the code often exceeds its proof.

🔄 Check Your Understanding. 1. In a real(dp) array a(n,n), which of a(i,j) and a(i+1,j) is adjacent in memory to a(i,j), and what does that imply about which loop index belongs innermost? 2. The fused step above eliminated the lap temporary. Roughly how much memory traffic per step did that save, in units of "one pass over the field"? 3. You have a loop whose body is y(i) = sqrt(x(i)); if (flag(i)) call log_it(y(i)). Why might loop fission make the sqrt part faster?

Answers

1. a(i+1,j) is adjacent (the first index strides by one element; the second strides by a whole column). So the loop over i must be innermost for unit-stride, cache-friendly access. 2. About two passes: it removed one full write of lap and one full read of lap (the temporary), leaving essentially one read of the old field and one write of the new — roughly halving the traffic of the two-sweep version. 3. The sqrt line is a clean, vectorizable, unit-stride kernel; the call log_it is an unpredictable branch to a procedure the compiler cannot inline, which blocks vectorization of the whole loop. Fissioning them lets the sqrt loop run at full SIMD width while the rare logging runs in its own loop.

29.2 Cache Blocking (Tiling)

Loop order fixes how you walk memory. Cache blocking fixes how much of it you try to keep close at once.

Definition (cache blocking; tiling). Cache blocking, or tiling, restructures a loop nest so that it works on the data in small blocks (tiles) that fit inside a fast cache level, finishing all the work that reuses a block before moving on. Instead of sweeping an entire large array end to end — evicting the front of it from cache long before a later pass needs it again — you sweep a sub-block small enough to stay resident, exhaust its reuse, and only then advance. The transformation adds outer loops over the blocks and shrinks the inner loops to a block's width.

The reason blocking matters — and the reason it sometimes does not — comes down to a single ratio: how many arithmetic operations a kernel does for each byte it loads. This is its arithmetic intensity, and it decides everything.

🚪 Threshold Concept — most numerical kernels are memory-bound, so optimization is mostly about moving less data, not doing less arithmetic. A modern CPU core can perform arithmetic far faster than main memory can feed it operands — the gap is often an order of magnitude or more. Draw the two ceilings on one graph (arithmetic rate on one axis, memory bandwidth on the other) and you have the roofline model: every kernel runs at whichever ceiling it hits first. A kernel with low arithmetic intensity — like the five-point stencil, which does a handful of adds per value fetched — hits the memory ceiling and is called memory-bound; no amount of cleverer arithmetic speeds it up, because the arithmetic units are already idling, waiting for data. A kernel with high arithmetic intensity — like matrix multiply, which reuses each loaded value $O(n)$ times — can hit the arithmetic ceiling and is compute-bound. Once you internalize this, optimization stops being a grab-bag of tricks and becomes one question: which ceiling am I under, and what moves me toward it? For memory-bound code, the answer is almost always "touch less memory." For compute-bound code, "keep the arithmetic units fed" — which is what blocking is for.

Matrix multiply is blocking's home turf, because it has reuse to capture. Computing $C = AB$ for $n \times n$ matrices does $O(n^3)$ arithmetic over only $O(n^2)$ data, so each element of $A$ and $B$ is used $n$ times. The naive triple loop uses an element, moves on, and by the time it needs that element again the cache has long since evicted it — so it is re-fetched from main memory, over and over. Blocking the three loops so they multiply small sub-tiles keeps each tile of $A$ and $B$ in cache for all $n_b$ of its reuses:

! Blocked matrix multiply C = A*B. Outer loops step by the block size nb;
! inner loops sweep within one tile, which stays resident in cache.
c = 0.0_dp
do jj = 1, n, nb
  do kk = 1, n, nb
    do ii = 1, n, nb
      do j = jj, min(jj+nb-1, n)
        do k = kk, min(kk+nb-1, n)
          do i = ii, min(ii+nb-1, n)          ! i innermost: column-major for C and A
            c(i,j) = c(i,j) + a(i,k)*b(k,j)
          end do
        end do
      end do
    end do
  end do
end do

The computed result is identical to the naive nest and to the intrinsic matmul — blocking only reorders when each multiply-add happens, not which ones. On a large matrix the payoff is real and often large, because you have converted a flood of cache misses into a trickle. code/example-03-blocked-matmul.f90 runs this on a hand-checkable $4 \times 4$ case (block size 2) and confirms it against matmul.

Now the honesty. Applied to your heat solver's single stencil sweep, cache blocking buys you almost nothing, and it is important to understand why so you do not waste an afternoon on it. The stencil is memory-bound: it streams through the field roughly once per step, doing about five flops per value, with each value reused only by its four immediate neighbors. A good loop order already delivers those neighbors in the same cache lines; there is no large-scale reuse for a tile to capture, because the kernel has no $O(n)$ reuse to begin with. Blocking helps when a kernel revisits data many times; a Jacobi sweep visits each point a small constant number of times and then leaves. The place blocking does return to the stencil is temporal blocking — fusing several time steps over a spatial tile so the field is reused across steps before leaving cache — but that is a substantially more complex, correctness-fraught transformation (the tiles must overlap to carry the halo forward), and it belongs to the specialist codes, not to a first optimization pass.

⚡ Performance Note: the practical decision rule falls straight out of the roofline. Before you block a loop, ask Chapter 28's question: is this kernel memory-bound or compute-bound? Tile the compute-bound, high-reuse kernels (dense linear algebra, and in 3D some stencils whose working set no longer fits in cache). Leave the memory-bound, low-reuse kernels alone — for them, loop order and fusion (§29.1) are the whole game, and you are already near the memory ceiling. And note the deeper lesson that §29.5 will collect on: the moment you find yourself hand-blocking a matrix multiply, you should be asking whether to call a library that has already done it far better than you will.

🔗 Connection: real BLAS and LAPACK implementations — OpenBLAS, Intel MKL, the reference kernels behind Chapter 21's dgesv — are blocked not for one cache level but for all of them at once (L1, L2, L3, registers), with a separately tuned tile size per microarchitecture. The blocked multiply above is a toy model of what those libraries do in hand-written assembly. Seeing the toy is worth it; shipping the toy is not.

29.3 SIMD Vectorization and do concurrent

The transformations so far reduce waiting. Vectorization reduces the number of instructions the CPU must issue to do the arithmetic, by doing several elements' worth at once.

Definition (SIMD). SIMD — Single Instruction, Multiple Data — is a hardware capability in which one machine instruction applies the same operation to a whole small vector of operands packed into a wide register. A 256-bit AVX register holds four real(dp) values; a 512-bit AVX-512 register holds eight. A single SIMD add can therefore add four or eight pairs of doubles in the time a scalar add handles one. Chapter 27 introduced vectorization as one of the things the compiler does; here we define the mechanism and, more importantly, how to make it happen.

Definition (auto-vectorization). Auto-vectorization is the compiler's automatic transformation of an ordinary scalar loop into one that uses SIMD instructions, processing several iterations per instruction — without you writing any special syntax, intrinsics, or assembly. You write a plain do loop; the compiler, at -O2/-O3, recognizes that the iterations are independent and unit-stride and emits vector instructions. This is the good path: portable, readable source that runs at vector speed. The bad path — writing the SIMD intrinsics by hand — is a last resort we return to in §29.4.

The compiler will auto-vectorize a loop only when it can prove the transformation is safe and profitable. Your job is to write loops that make the proof easy. The conditions are:

• No loop-carried dependency. Iteration $i$ must not need a value that iteration $i-1$ wrote. y(i) = a*x(i) + y(i) vectorizes; y(i) = y(i-1) + x(i) (a running sum) does not, because each step needs the previous result.
• Unit-stride, contiguous access. The loop should walk memory in steps of one element — which, in Fortran, means the inner loop over the first index (§29.1). Strided or gathered access defeats or slows vectorization.
• A countable trip count. The compiler must know, or compute at entry, how many iterations there are — a plain do i = 1, n qualifies; a do while with a data-dependent exit usually does not.
• No aliasing between what you read and what you write. Here Fortran hands you a gift C cannot: procedure arguments are assumed not to alias (Chapter 27), so the compiler may assume the input and output arrays are distinct and vectorize freely.
• A simple body. No un-inlinable calls, minimal branching. elemental and pure procedures (Chapter 6) inline and vectorize; a call to a side-effecting subroutine blocks the whole loop (and is exactly what §29.1's loop fission isolates).

Written this way, the stencil's inner loop meets every condition, and the array-section form of the Laplacian you have used since Chapter 5 is if anything easier for the compiler, because whole-array syntax states independence outright. You do not usually need to do anything beyond writing clean, unit-stride, dependency-free loops and turning optimization on.

⚡ Performance Note — confirm it, do not assume it. Whether a loop actually vectorized is not a matter of faith. Ask the compiler for its optimization report (Chapter 27): gfortran -O3 -fopt-info-vec prints which loops it vectorized, and -fopt-info-vec-missed prints the ones it declined and the reason ("not enough data-parallelism", "possible aliasing", "unsupported use in outer loop"). Reading that report is how you turn "I hope this is fast" into "the compiler vectorized the inner loop; the outer one is memory-bound, as expected." The full menu of flags — -O3, -march=native to unlock AVX-512 on your specific chip — is Chapter 30's subject; for now, know that vectorization needs -O2 or higher and the right -march to use the widest registers your CPU has.

do concurrent for performance

Sometimes you know two iterations are independent but the compiler cannot prove it, and it conservatively refuses to vectorize. Fortran gives you a way to say so in the language itself. You met do concurrent as a control construct back in Chapter 4; here is its performance purpose.

Definition (do concurrent, performance use). do concurrent is a loop form in which you assert to the compiler that the iterations may be executed in any order, or all at once, with no iteration depending on another. It is a promise you make and the compiler trusts — it does not re-prove independence, so it is free to vectorize, unroll, or (with the right flags or compiler) run the iterations across multiple threads or a GPU. Crucially, do concurrent does not by itself guarantee parallel or vector execution; it removes the compiler's need to prove independence and thereby permits optimizations it would otherwise forgo. The independence promise is yours to keep — break it (let one iteration read another's write) and the result is undefined.

The stencil update expressed with do concurrent says exactly what is true — every interior point's new value depends only on old neighbors, so the points may be computed in any order:

do concurrent (j = 2:ny-1, i = 2:nx-1)               ! "these are all independent"
  u_new(i,j) = u(i,j) + rx*(u(i-1,j) - 2.0_dp*u(i,j) + u(i+1,j))  &
                      + ry*(u(i,j-1) - 2.0_dp*u(i,j) + u(i,j+1))
end do

Note that it writes into a separate u_new: the promise of independence is true only because no iteration reads a value another iteration wrote. Overwrite u in place and you would break the promise (and, as Chapter 24 warned, silently change the physics from Jacobi to Gauss–Seidel). Kept honestly, this one construct is your on-ramp to parallelism: the same loop that a compiler may vectorize today is the loop that Chapter 33 will hand to a team of OpenMP threads and Chapter 35 may offload to a GPU, with little or no change to the body.

⚠️ Common Pitfall — do concurrent is a promise, not a parallel directive. Two mistakes recur. First, people expect do concurrent to automatically run in parallel; on most compilers today, plain do concurrent compiles to a serial loop that is merely easier to vectorize — you need -ftree-parallelize-loops, or a compiler like nvfortran targeting a GPU, to actually spread it across hardware. Treat the speed it buys as "better vectorization," and use the explicit models of Part VIII when you want guaranteed parallelism. Second, and worse, people put a hidden dependency inside it — a running total, an in-place update — and the compiler, trusting the promise it does not check, produces wrong answers only sometimes, only at high optimization, only on some machines. If the iterations are not truly independent, do not use do concurrent.

🔄 Check Your Understanding. 1. Name two properties an inner loop must have for the compiler to auto-vectorize it. 2. What does do concurrent promise the compiler, and what does the compiler give back in return? 3. Why is y(i) = y(i-1) + x(i) not vectorizable as written, while y(i) = y(i) + a*x(i) is?

Answers

1. Any two of: iterations are independent (no loop-carried dependency); access is unit-stride/contiguous (inner loop over the first index); the trip count is countable; no aliasing between read and write arrays; the body has no un-inlinable calls or unpredictable branches. 2. You promise that the iterations are independent and may run in any order; in return the compiler is freed from proving independence itself and may vectorize, reorder, unroll, or (with flags/compiler) parallelize the loop. It does not automatically make it parallel. 3. The first has a loop-carried dependency: iteration $i$ needs y(i-1), the result of iteration $i-1$, so the iterations must run in order — a scalar chain. The second's iterations are independent (each touches only its own y(i)), so several can be done at once in a SIMD register.

29.4 Contiguous Pointers, Aliasing, and When Hand-Tuning Wins

Vectorization rests on two assumptions the compiler must be able to make: that memory is walked in unit strides, and that the arrays you read do not overlap the arrays you write. Modern Fortran gives you two ways to guarantee those assumptions when the compiler cannot see them for itself.

The contiguous attribute

An assumed-shape dummy argument (real(dp), intent(in) :: u(:,:)) is wonderfully general — it accepts an array section like field%u(2:n-1, ::2), which may be strided or non-contiguous in memory. That generality has a cost: at the top of the procedure the compiler cannot assume unit stride, so it may insert a runtime check, or make a contiguous copy-in/copy-out of the argument, or simply decline to vectorize. When you know the array is contiguous, you can say so with the contiguous attribute from Chapter 11:

subroutine stencil(u, u_new, rx, ry)
  real(dp), intent(in),    contiguous :: u(:,:)       ! promise: unit stride, no copy needed
  real(dp), intent(inout), contiguous :: u_new(:,:)
  real(dp), intent(in) :: rx, ry
  ! ... unit-stride inner loop, now free to vectorize without a contiguity check ...
end subroutine stencil

contiguous is most valuable on pointers. A pointer can be associated with a strided slice of a larger target, so by default the compiler must assume a pointer array is non-contiguous and load through it carefully. Declaring the pointer contiguous promises it always points at a unit-stride block, restoring the unit-stride assumption vectorization needs. This is precisely why the Project Checkpoint uses contiguous pointers for its two field buffers: it wants both the O(1) pointer swap and the unit-stride guarantee.

Aliasing hints

Fortran's default that procedure arguments do not alias is the no-aliasing advantage Chapter 27 built its case on — it is why the compiler may assume the stencil's input and output are distinct and vectorize the loop, with no restrict keyword needed as in C. But you can forfeit that advantage. Pointers and targets can alias: if you pass the same target to two pointer dummies, or overlap two sections, the compiler must assume the worst and serializes. The practical rule follows from the whole shape of the book: prefer allocatable to pointer (Chapter 11), keep your reads and writes in genuinely distinct arrays, and you keep the no-aliasing gift that makes Fortran fast. Reach for pointers only when you need them — and when you do, contiguous and honest non-overlap are how you pay the performance back.

When hand-tuning beats the compiler (rarely)

Everything so far has been about helping the compiler — better loop order, fewer passes, honest promises of independence and contiguity. The question inevitably arises: should you ever go over the compiler's head and tune by hand — write the SIMD intrinsics yourself, unroll the loop manually, insert prefetch instructions?

The honest answer is rarely, and only after the profiler proves it. A modern compiler at -O3 -march=native is a formidable optimizer that has internalized decades of these tricks; for the overwhelming majority of loops, hand-tuning produces code that is longer, less portable, harder to maintain, and no faster — or slower, because you defeated an optimization the compiler would have done better. The cases where hand-tuning genuinely wins are narrow and specialist:

• The compiler provably fails to vectorize a hot, proven-critical loop, and the optimization report (-fopt-info-vec-missed) shows a reason you cannot fix through source structure.
• A kernel needs a SIMD operation with no straightforward scalar expression (a shuffle, a horizontal reduction, a specific rounding mode).
• You are writing the innermost kernel of a library that thousands of programs will call — a BLAS dgemm, an FFT butterfly — where a 20% gain, amortized over the world, is worth a month of one expert's assembly.

That last case is the tell. If your loop is important enough to hand-tune, it is usually important enough that someone has already done it and put it in a library. Which is exactly where §29.5 is headed.

🐛 Find the Bug. A well-meaning colleague "optimized" the stencil by making its buffers pointers so they could swap them cheaply, and now the vectorized -O3 build is slower than the plain -O2 one, and the optimization report says "loop not vectorized: possible aliasing." Here is the declaration:

fortran real(dp), pointer :: u(:,:), u_new(:,:) ! ... u and u_new each point into one big allocated slab, via sections ... call stencil(u, u_new, rx, ry) ! stencil takes intent(in) u(:,:), intent(inout) u_new(:,:)

What did the pointers cost, and what one attribute fixes it?

Answer

Plain pointer arrays are assumed possibly non-contiguous and possibly aliasing — so the compiler cannot assume unit stride or that u and u_new are distinct, and it refuses to vectorize the inner loop (and may insert copy-in/copy-out). The fix is to declare them contiguous: real(dp), contiguous, pointer :: u(:,:), u_new(:,:). That restores the unit-stride guarantee; and because they point at separate slabs, non-overlap holds. (If they genuinely could overlap, contiguous would not be enough — you would have to redesign so the read and write arrays are distinct, which they must be for a correct Jacobi step anyway.) One attribute, and the vectorizer comes back.

29.5 When to Stop: Diminishing Returns and Why BLAS Wins

Optimization has no natural end — there is always one more percent to chase — so knowing when to stop is as much a professional skill as knowing how to start. Three honest ideas govern the decision.

Diminishing returns are real and steep. The speedups in this book's techniques are not equal. Turning on -O2/-O3 and fixing loop order routinely recovers a large factor for near-zero effort and zero readability cost. Fusion and hoisting add a solid further chunk for a small, local rewrite. Blocking a compute-bound kernel can be large but costs real complexity. Hand-written SIMD intrinsics might add a final tens-of-percent for a large cost in readability and portability. Plotted as effort against speedup, the curve bends hard: the first hour buys most of the win, and each hour after buys less. Knuth's "critical 3%" is precisely the small region where the curve is still worth climbing — and the corollary is that for the other 97% of your code, the fastest and cheapest thing is clean, correct, readable code that you never optimized at all.

Readability is a cost you are spending. A 2× speedup that turns a ten-line kernel anyone can read into a sixty-line kernel only you understand is often a bad trade, because the sixty-line version will rot: the next person (possibly you, in a year) will be afraid to touch it, bugs will hide in it, and its "optimization" will outlive the hardware it was tuned for. The discipline of Part IX — that scientific software must be maintained, sometimes for decades — means an optimization must justify its complexity, not just its speedup. Optimize the proven hot loop; leave the other 97% legible.

⚡ Performance Note — you cannot beat the memory ceiling. For a memory-bound kernel there is a hard upper bound on how fast it can possibly go: the time to move its data through the memory system, once. If your stencil already streams the field the minimum number of times, in the right order, you are at the roofline's memory ceiling, and further arithmetic cleverness is physically incapable of helping — the arithmetic units are already idle. Recognizing that you have hit the ceiling is what tells you to stop (or to change the algorithm, or to add hardware — more memory bandwidth, or Part VIII's parallelism). This is why Chapter 28's memory-bound / compute-bound diagnosis is the prerequisite to this whole chapter: it tells you which ceiling you are under, and therefore whether there is any room left above you at all.

And the deepest lesson: a tuned BLAS still beats you. We close the performance part where we opened the book. Back in Chapter 21 you wrote your own matrix multiply and then called LAPACK's, and the promise was made that the library would win. Here is why it wins, in the language of this chapter. A production dgemm from OpenBLAS or Intel MKL is: blocked for every level of the cache hierarchy at once (registers, L1, L2, L3), with tile sizes tuned per microarchitecture; hand-vectorized in assembly to saturate the SIMD units and hide latency by software pipelining; and structured to run at a large fraction of the machine's theoretical peak arithmetic rate — a number a naive triple loop reaches perhaps a few percent of. The gap between your best-blocked Fortran matrix multiply and a tuned dgemm is routinely an order of magnitude, and it represents person-decades of specialist effort you cannot and should not reproduce.

🔗 Connection: this is the same theme from a new angle. In Chapter 1 we noted that NumPy is fast because it calls down into compiled Fortran and C — BLAS and LAPACK. Now you can say precisely what that compiled code is doing that your loop is not: multi-level cache blocking and hand-vectorization to hit the arithmetic roofline. The professional move is not to compete with it but to recognize its shape — "this is a dense linear-algebra kernel" — and call it. Hand-roll the physics that is uniquely yours (your stencil, your boundary conditions); never hand-roll the dgemm underneath it.

📜 From History: the BLAS (Basic Linear Algebra Subprograms) were specified in the late 1970s precisely so that the machine-specific optimization of a handful of core operations — a dot product, a matrix multiply — could be done once, by a vendor's experts, behind a fixed interface, and every scientific program above them would inherit the speed for free. That 1979 design decision is why, fifty years later, the right answer to "how do I make my matrix multiply fast" is still "don't; call dgemm." The interface outlived a dozen generations of hardware, and the tuning moved underneath it without a single caller changing a line. It is the same standardization instinct that Chapter 1 credited for Fortran's own longevity.

🔄 Check Your Understanding. 1. You have optimized a kernel and the profiler shows it is memory-bound and already streaming its data the minimum number of times. Should you keep optimizing the arithmetic? Why or why not? 2. Give two reasons a tuned dgemm beats a well-blocked matrix multiply you wrote yourself. 3. State Knuth's "critical 3%" idea in your own words as a rule for where to spend optimization effort.

Answers

1. No. If it is memory-bound and at the minimum traffic, it is at the roofline's memory ceiling — the arithmetic units are already idle waiting for data, so making the arithmetic cleverer cannot help. To go faster you must move less data (change the algorithm), add memory bandwidth, or parallelize. 2. Any two: it is blocked for all cache levels simultaneously (not just one); it is hand-vectorized in assembly to saturate the SIMD units; it is tuned per microarchitecture; it hides memory latency by software pipelining — together reaching a large fraction of peak versus a naive loop's few percent. 3. Profile first; spend optimization effort only on the small fraction of code that the measurement proves is hot (the "critical 3%"), and leave the other 97% clean and readable — optimizing it would cost complexity for a speedup too small to matter.

Project Checkpoint

Your solver computes correct physics (Chapter 24) and you have profiled it (Chapter 28) and found — as every explicit stencil code does — that essentially all the time is inside step. Now we optimize step, applying this chapter's techniques in the order of their payoff, and we do it under one iron constraint: the optimized step must produce numerically identical results to the Chapter 24 version. An optimization that changes the answer is not an optimization; it is a bug.

Four changes, from largest payoff to smallest:

1. Correct loop order (§29.1). The inner loop runs over the first index i, so the sweep walks each column contiguously. This is the single biggest, most reliable win — do it first, always.
2. Loop fusion (§29.1). Compute the Laplacian and apply the update in one sweep, with the divisions hoisted into rx = alpha*dt/dx**2 and ry = alpha*dt/dy**2. The whole lap(:,:) temporary disappears, halving memory traffic.
3. Two buffers, swapped with contiguous pointers (§29.4). Instead of allocating a work array inside step every call (an $O(n^2)$ allocation and copy per step), the driver owns two persistent buffers and swaps two contiguous pointers each step — an $O(1)$ swap, no data copied, unit stride guaranteed.
4. do concurrent (§29.3). The interior update is expressed as do concurrent, asserting independence so the compiler may vectorize now and parallelize in Chapter 33 behind the same code.

Here is the optimized kernel and its buffer swap (the full self-contained program, hand-verified, is code/project-checkpoint.f90):

subroutine stencil(u, u_new, rx, ry)                  ! fused, column-major, contiguous
  real(dp), intent(in),    contiguous :: u(:,:)
  real(dp), intent(inout), contiguous :: u_new(:,:)   ! interior written; Dirichlet edges preserved
  real(dp), intent(in) :: rx, ry
  integer :: i, j, nx, ny
  nx = size(u,1);  ny = size(u,2)
  do concurrent (j = 2:ny-1, i = 2:nx-1)              ! independent: reads OLD u, writes NEW u_new
    u_new(i,j) = u(i,j) + rx*(u(i-1,j) - 2.0_dp*u(i,j) + u(i+1,j))  &
                        + ry*(u(i,j-1) - 2.0_dp*u(i,j) + u(i,j+1))
  end do
end subroutine stencil

! ... in the driver, per step: swap the two contiguous field pointers (O(1), no copy) ...
call stencil(u, u_new, rx, ry)
tmp => u;  u => u_new;  u_new => tmp

Both buffers have their Dirichlet edges set once before the time loop, so swapping preserves them forever and u_new's intent(inout) (not intent(out), which would wipe the edges) leaves them untouched. On the Chapter 24 test — a $5\times5$ plate, hot top edge at $100$, $\alpha = 1$, $\Delta t = 0.2$ (so $r = 0.2 \le \tfrac14$) — two steps produce, hand-computed, exactly the Chapter 24 field:

after step 1:                          after step 2:
  100.00 100.00 100.00 100.00 100.00     100.00 100.00 100.00 100.00 100.00
    0.00  20.00  20.00  20.00   0.00       0.00  28.00  32.00  28.00   0.00
    0.00   0.00   0.00   0.00   0.00       0.00   4.00   4.00   4.00   0.00
    0.00   0.00   0.00   0.00   0.00       0.00   0.00   0.00   0.00   0.00
    0.00   0.00   0.00   0.00   0.00       0.00   0.00   0.00   0.00   0.00

That the numbers match Chapter 24 to the digit is the whole point: every transformation here — reorder, fuse, hoist, swap, do concurrent — merely changes when and in what order the same arithmetic happens, never which arithmetic. Because each output cell is an independent expression in old-field values (not a reduction), reordering the iterations is bit-for-bit identical, not merely close. That is how you verify an optimization: run the naive and optimized versions on the same input and assert maxval(abs(u_opt - u_naive)) == 0 (for this dependency-free kernel) — and had the kernel been a sum, you would assert a small tolerance instead, because reassociating a floating-point reduction can change the last bit.

Now re-measure — and read the result honestly. Wrap the step in the tic/toc timer from Chapter 28's timers module and time a long run on a large grid. The illustrative shape of what you will see (Tier 2 — your numbers will differ; measure them):

The precise figures are meaningless out of context; the pattern is the lesson. Loop order alone recovers the largest factor. Fusion adds a solid further slice. Vectorization adds a final one — and then the curve flattens hard, because the stencil is memory-bound: once you are streaming the field once, in order, vectorized, you are near the memory roofline and there is little left to win from this kernel on this hardware. That flattening is not failure; it is the diminishing-returns signal of §29.5 telling you to stop here and, when you need more, reach for Part VIII's parallelism — which the do concurrent you just wrote has already prepared you for. Notice, too, what we did not do: no cache blocking (§29.2 explained it buys a memory-bound single sweep almost nothing) and no hand-written intrinsics (§29.4: not worth it here). Choosing not to apply a technique, for a stated reason, is as much a part of optimization as applying one.

Summary

This chapter turned Chapter 28's measurement into action: the concrete transformations that make numerical Fortran fast, in order of payoff, and the judgment to know when to stop.

Key rules and numbers to carry forward:

• Inner loop over the first index — column-major, unit stride. The one rule with no exceptions.
• Optimize in payoff order: -O2/-O3 + loop order → fusion + hoisting → vectorization → (rarely) anything else. The first hour buys most of the speedup.
• Roofline / memory-bound: most stencils are memory-bound; you cannot out-compute the memory ceiling, so optimize by moving less data, and stop when you reach it.
• do concurrent is a promise, not a parallel directive — keep the promise (true independence, no in-place update) or get undefined results.
• Verify every optimization: identical result required. Dependency-free kernels → bit-identical (== 0); reductions → within tolerance (FP reassociation).
• Don't hand-roll dgemm. For dense linear algebra, call the tuned BLAS (Chapter 21); it beats you by an order of magnitude.

The heat solver's step now runs in the correct loop order, in a single fused sweep, over swapped contiguous buffers, expressed as do concurrent — fast, still exactly correct, and ready to go parallel.

Spaced Review

Retrieval practice on the chapters this one is built from: arrays and memory layout (Chapter 5), why Fortran is fast (Chapter 27), and profiling (Chapter 28). Answer before peeking.

1. (Ch. 5) Fortran stores a(n,n) column-major. Which two elements are adjacent in memory, a(i,j) and a(i+1,j), or a(i,j) and a(i,j+1) — and how does that answer dictate every loop-order decision in this chapter?
   

   Answer

   `a(i,j)` and `a(i+1,j)` are adjacent — the *first* index strides by one element; the second strides by a whole column of `n` elements. Therefore the inner loop must run over the first index `i` for unit-stride access, which is the foundation of §29.1 (reorder), §29.3 (vectorization needs unit stride), and §29.4 (`contiguous`).

2. (Ch. 5) Why can a whole-array or array-section statement like u_new(2:n-1,2:n-1) = … be easier for the compiler to vectorize than the equivalent hand-written do loop?
   

   Answer

   Array syntax states the operation on whole (sub)arrays at once, which declares that the elements are independent and the operation is elementwise — exactly the information the compiler needs to vectorize. A `do` loop conveys the same intent only if the compiler can *prove* independence and unit stride; the array form hands it that proof for free (and is why array sections remain a first-class habit even in optimized code).

3. (Ch. 27) Fortran assumes procedure arguments do not alias — the "no-aliasing advantage." How does that assumption directly enable the vectorization of §29.3, and what Fortran feature can forfeit it?
   

   Answer

   Because the compiler may assume the input array and the output array of the stencil do not overlap, it can load, compute, and store several elements at once (SIMD) without worrying that a store changed a value it is about to load. That is exactly what auto-vectorization needs. **Pointers/targets** can alias (two pointers may reference the same storage), which forfeits the guarantee — hence the §29.4 advice to prefer `allocatable`, keep read/write arrays distinct, and mark pointers `contiguous`.

4. (Ch. 28) Before optimizing a loop, Chapter 28 said to classify it as memory-bound or compute-bound. Why is that classification the prerequisite to this chapter — what different advice does each answer lead to?
   

   Answer

   It tells you which roofline ceiling you are under, and therefore which optimizations can possibly help. **Memory-bound** (low arithmetic intensity, e.g., a stencil) → reduce memory traffic: fix loop order, fuse, minimize passes; cleverer arithmetic and cache blocking of a single sweep do nothing. **Compute-bound** (high reuse, e.g., matmul) → keep the arithmetic units fed: cache blocking and vectorization pay, and calling a tuned BLAS pays most. Optimizing without this classification wastes effort on the wrong ceiling.

5. (Ch. 28) You apply an optimization and the wall-clock time drops, but you changed loop order, ran once, and did not warm the cache. Name two methodology errors and how Chapter 28 said to avoid them.
   

   Answer

   Any two of: (i) *single run* — one measurement has no notion of variance; repeat and report a best/median over several runs. (ii) *cold cache / no warm-up* — the first run pays one-time costs (page faults, cache fills); discard a warm-up iteration. (iii) *unrepresentative size* — if the benchmark grid fits in cache but the real one does not, the speedup will not transfer; measure at the real problem size. (iv) *no baseline control* — confirm both versions compute the same result before trusting the timing.

What's Next

You have made the source fast: the right loop order, one fused sweep, vectorizable and independent. But the same source can run at wildly different speeds depending on what you tell the compiler to do with it — and so far we have leaned on -O3 and -march=native as incantations without unpacking them. Chapter 30 is that unpacking: the optimization levels -O2, -O3, and the sharp-edged -Ofast; -march=native and what it trades for portability; -flto for whole-program optimization; the differences between gfortran, Intel's ifx, and NVIDIA's nvfortran; and the professional discipline of recording the exact flags so a result is reproducible. Fast code and the right flags are two halves of one craft; you have the first half. Let's get the compiler to finish the job.