Chapter 29 Exercises — Optimization Techniques

These exercises turn the chapter's techniques into muscle memory: reorder a loop for column-major access, fuse two passes into one, block a matrix multiply, coax the compiler into vectorizing, and — as often as not — decide not to optimize because the profiler says there is nothing to win. The through-line is the chapter's iron rule: an optimization must never change the answer. Several problems ask you to prove that it did not.

The difficulty tiers are:

  • ⭐ Foundational — one idea, short code, predict-then-run.
  • ⭐⭐ Applied — combine loop transformations, blocking, or vectorization into working code and measure.
  • ⭐⭐⭐ Challenge — a roofline estimate, a subtle dependency, or a design decision; expect to think.

Problems marked † have full worked solutions in appendices/answers-to-selected.md (and, for code, in this chapter's code/exercise-solutions.f90). Odd-numbered problems are also solved there. Everything compiles with gfortran -std=f2018 -Wall. Two standing rules apply with special force in a performance chapter: you may not trust a speedup you have not measured (Chapter 28's methodology — warm up, repeat, measure at the real size), and you must verify correctness before you trust timing — a fast wrong answer is worthless.

Part A — Type, Compile, and Run (predict first)

A1 ⭐† Type in example-01-loop-order.f90. Before compiling, predict the interior Laplacian and the value of max |fast - slow|. Then compile and run. In one sentence, say why the two loop orders print the identical field, and in a second sentence say what would differ if the grid were $4000\times4000$ instead of $5\times5$.

A2 ⭐ Predict the output of example-02-fusion-doconcurrent.f90 — the field after one step and both max differences — before running it. Which of the three methods (unfused, fused, do concurrent) computes a different answer, and why is the correct answer "none of them"?

A3 ⭐⭐† Run example-03-blocked-matmul.f90. Predict the $4\times4$ result and max |blocked - matmul| first. Then change the block size nb from 2 to 3, and to 4 (one whole tile). Does the printed matrix change? Explain what nb does and does not affect.

A4 ⭐⭐ Take the Project Checkpoint's optimized stencil. Without running it, hand-compute the field after three steps (continue the $28,32,28$ / $4,4,4$ pattern one more step), then confirm by adding a third call stencil + swap to project-checkpoint.f90. Which interior cells are non-zero after step 3?

Part B — Loop Transformations

B5 ⭐† The fragment below is a two-pass update: it builds a whole grad2(:) array, then applies it. Fuse it into a single loop with no temporary array, hoisting the loop-invariant c = k/dx**2 out of the loop. Confirm your fused version prints the identical y.

do i = 2, n-1
  grad2(i) = ( x(i-1) - 2.0_dp*x(i) + x(i+1) ) / dx**2
end do
do i = 2, n-1
  y(i) = x(i) + k*grad2(i)
end do

B6 ⭐⭐ Count the divisions. In the two-pass fragment of B5, how many floating-point divisions execute per call for n = 1000? In your fused, hoisted version, how many? Explain why a divide is worth hoisting when a multiply usually is not.

B7 ⭐⭐† Loop fission. The loop below mixes a clean, vectorizable computation with a rare, branch-and-call side effect that blocks vectorization of the whole loop:

do i = 1, n
  s(i) = sqrt(x(i)*x(i) + y(i)*y(i))
  if (s(i) > threshold) call record_outlier(i, s(i))
end do

Split it into two loops so the sqrt computation can vectorize, and explain what fission bought and what it cost (there is a cost — name it).

B8 ⭐⭐⭐ Fusion is not always legal. Give a concrete pair of adjacent loops over the same range that you cannot fuse without changing the result, and state the dependency that forbids it. (Hint: consider a second loop that reads an element the first loop has not written yet at that iteration.)

Part C — Cache Blocking and the Roofline

C9 ⭐† Compute the arithmetic intensity (flops per real(dp) byte loaded, order of magnitude) of (a) the five-point stencil update and (b) an $n\times n$ matrix multiply. Which is memory-bound, which is compute-bound, and which one is worth cache blocking? One or two sentences each.

C10 ⭐⭐† Using example-03-blocked-matmul.f90 as a starting point, write a blocked(a, b, c, n, nb) subroutine and verify it against matmul for nb = 2, 3, 4 on a $6\times6$ matrix of your choice. Report max |blocked - matmul| for each nb. Then state: on a large matrix, why does the right nb matter for speed even though it never matters for the answer?

C11 ⭐⭐ Back-of-envelope block size. Your L1 data cache is $32\ \mathrm{KiB}$. For a blocked matrix multiply that must keep three tiles resident (one each of $A$, $B$, $C$), roughly how large can a square real(dp) tile be before it overflows L1? Show the arithmetic. Why is the real best block size usually found by measurement, not this estimate?

C12 ⭐⭐⭐† Roofline crossover. A CPU core sustains $50\ \mathrm{Gflop/s}$ of real(dp) arithmetic and $20\ \mathrm{GB/s}$ of memory bandwidth. Compute the machine balance (flops per byte at which it switches from memory-bound to compute-bound). A kernel does $0.25$ flops per byte loaded — memory- or compute-bound? What does that tell you about whether to optimize its arithmetic or its memory traffic?

Part D — SIMD and do concurrent

D13 ⭐† For each loop body, say whether it can auto-vectorize as written, and if not, name the obstacle: (a) y(i) = 2.0_dp*x(i) + z(i); (b) y(i) = y(i-1) + x(i); (c) y(i) = x(i); if (mask(i)) call f(y(i)); (d) total = total + x(i).

D14 ⭐⭐ The running sum s(i) = s(i-1) + x(i) has a loop-carried dependency. It nonetheless has a fast parallel form (a prefix sum / scan). Without writing it, explain in two sentences why the naive loop cannot vectorize but the total sum(x) (a plain reduction) can be computed with SIMD by partial sums.

D15 ⭐⭐† Rewrite the stencil sweep as a do concurrent loop. Then state precisely the promise you are making to the compiler, and exhibit a one-line change to the body that would break that promise (make the iterations no longer independent) and thereby make the do concurrent produce undefined results.

D16 ⭐⭐ You compile with gfortran -O3 -fopt-info-vec-missed and see: "loop not vectorized: possible aliasing." The loop is subroutine axpy(y, x, a, n) with y and x both intent(inout) pointer arrays. Explain the message and give a source change (an attribute or a redesign) that lets it vectorize.

Part E — Design It (extend the solver)

E17 ⭐⭐† This is the Project Checkpoint as an exercise. Starting from Chapter 24's step, produce an optimized version that (1) uses the correct loop order, (2) fuses the Laplacian and update, (3) hoists the divisions, and (4) is expressed as do concurrent. Then write a verification harness: run the Chapter 24 step and your optimized version on the same random field and error stop unless maxval(abs(u_opt - u_ref)) == 0.0_dp. Why is exactly zero (not a tolerance) the right assertion here?

E18 ⭐⭐ Give the solver's stencil kernel contiguous dummy arguments and explain, in terms of the optimization report, what that attribute lets the compiler stop worrying about. On what kind of actual argument would passing a non-contiguous array to a contiguous dummy be a bug — and how would -fcheck=all help?

E19 ⭐⭐⭐ Two verifications. The stencil update is dependency-free, so the optimized version is bit-identical to the naive one. A reduction — say, computing the total heat sum(u) — is not: SIMD partial sums reassociate the additions and can change the last bit. Write both assertions: a bit-exact == 0 check for the stencil, and a tolerance check <= tol for the reduction, and justify the tolerance you chose in terms of machine epsilon (Chapter 20).

Part F — Port It / Modernize It

F20 ⭐⭐† Port it. Here is a NumPy kernel a colleague wrote (row-major, C order):

for j in range(1, m-1):          # NumPy: last axis is contiguous
    row_new[j] = row[j] + c*(row[j-1] - 2*row[j] + row[j+1])

Port the whole 2D version to Fortran on a field%u(:,:), and — this is the point — say which Fortran loop index must be innermost for cache efficiency, and how that differs from the innermost axis you would choose in NumPy. One kernel, mirrored memory order.

F21 ⭐⭐ Modernize it. The FORTRAN 77 matrix multiply below has the wrong loop order for Fortran's column-major storage (the innermost loop strides across rows). Modernize it to the book's style and reorder the loops for unit-stride access without changing the result. Which loop becomes innermost, and why?

      DO 30 I = 1, N
      DO 20 J = 1, N
      DO 10 K = 1, N
      C(I,J) = C(I,J) + A(I,K)*B(K,J)
   10 CONTINUE
   20 CONTINUE
   30 CONTINUE

Part G — Back of the Envelope

G22 ⭐† Memory traffic. For an $n\times n$ real(dp) field, the Chapter 24 two-pass step reads and writes the field and a full lap temporary; the fused step does neither for lap. Estimate the bytes of memory traffic per step for each, for $n = 1000$, and the fractional saving. Why does that fraction directly predict the speedup for a memory-bound kernel?

G23 ⭐⭐ SIMD ceiling. An AVX register is $256$ bits; AVX-512 is $512$. How many real(dp) values fit in each? What is the ideal vectorization speedup for a perfectly SIMD-friendly loop on each, and name two concrete reasons you will not actually reach it.

G24 ⭐⭐† Divisions saved. Your solver runs $n = 500$, so the interior is $498\times498$ points, for $10{,}000$ time steps. How many floating-point divisions does hoisting the two stencil divisions (1/dx**2, 1/dy**2) out of the inner loop save over the whole run? Express it as a count and comment on whether it is worth the one-line change.

Part H — Interleaved (earlier chapters)

H25 ⭐† (Ch. 5) A colleague's kernel is do i = 1, n; do j = 1, n; a(i,j) = b(i,j) + c(i,j); end do; end do and it is several times slower than it should be on a large array. Name the problem in one phrase, give the corrected loop nesting, and justify it from Fortran's column-major storage.

H26 ⭐⭐ (Ch. 27) The no-aliasing advantage is why the stencil vectorizes. Explain, in three or four sentences, how passing two distinct intent array arguments lets the compiler vectorize, and how switching those arguments to pointer (that might alias) would force it to serialize — the exact hazard §29.4 guards against with contiguous.

H27 ⭐⭐† (Ch. 28) You profile the solver and gprof reports 96% of time in step, 3% in I/O, 1% elsewhere. Cachegrind shows step has a high last-level-cache miss rate. Which optimization from this chapter do you apply first, and why does the profile tell you not to bother optimizing the I/O?

H28 ⭐⭐⭐ (Ch. 21) You have hand-written and cache-blocked a matrix multiply, and it reaches perhaps 8% of your CPU's peak real(dp) rate. A tuned dgemm from OpenBLAS reaches ~80%. Estimate the speedup of switching to dgemm, list three things it does that your blocked loop does not, and state the one-line professional conclusion of §29.5.

Solutions to the starred, odd-numbered, and †-marked problems are in appendices/answers-to-selected.md; the compilable ones are in code/exercise-solutions.f90. If your output disagrees with a hand-computed "Expected output," verify correctness before you touch timing — a fast wrong answer has taught you nothing.