Exercises: Why Fortran Is Fast

This is the chapter where performance stops being a slogan and becomes something you can predict, explain, and measure. The exercises reflect that: some ask you to reason about memory and the compiler, some ask you to port or fix code, and a few ask you to estimate — the back-of-the-envelope arithmetic that separates people who guess about performance from people who know. Do the estimates even when they feel imprecise; an order-of-magnitude answer you can defend beats a precise number you cannot.

Difficulty: ⭐ warm-up · ⭐⭐ standard · ⭐⭐⭐ deeper. Solutions: worked solutions to the daggered (†) and odd-numbered problems are in appendices/answers-to-selected.md; the compilable ones are in code/exercise-solutions.f90. Never run a "predict the output" program before you have predicted it — that is the whole exercise. And remember the chapter's honesty rule: any timing you produce is yours to measure; the outputs you can compute by hand.

Part A — Warm-ups ⭐

27.1 † In one line each, say what -O0, -O2, -O3, and -Ofast do, and name the single most important thing that changes between -O2 and -O3, and between -O3 and -Ofast.

27.2 For a Fortran 2D array a(:,:), which loop index belongs on the inner loop, and in one sentence, why?

27.3 † State the no-aliasing advantage in one sentence, from the compiler's point of view.

27.4 What does -Ofast add on top of -O3, and give one reason to be cautious about it in numerical code. (Link your answer to Chapter 20.)

27.5 † Name the gfortran flag that lists the loops the compiler successfully vectorized, and the one that lists loops it failed to vectorize with a reason.

27.6 True or false, with a one-sentence justification: "At -O3, gfortran unrolls all your loops automatically."

Part B — Memory and Loop Order ⭐⭐

27.7 † A cache line is 64 bytes and a real(dp) is 8 bytes. For an array a(512, 512) of real(dp): (a) how many bytes is one column, and how many cache lines does it span? (b) When you loop across a row (inner loop over the last index), roughly how many cache lines do you touch to read one full row of 512 elements, and why is that the slow pattern?

27.8 Write the cache-friendly (with-the-grain) triple-nested loop that initializes a 3D array u(nx, ny, nz) of real(dp) to zero using explicit indices. Which index is innermost, and which outermost?

27.9 † Explain what "memory-bound" means, and why, for a memory-bound loop, reducing the number of arithmetic operations may not make it any faster. What would make it faster?

27.10 Predict which of these two nests is faster on a large a(n,n), and state the direction of the factor (you need not give an exact number):

! Nest A
do i = 1, n
  do j = 1, n
    a(i,j) = a(i,j) * 2.0_dp
  end do
end do

! Nest B
do j = 1, n
  do i = 1, n
    a(i,j) = a(i,j) * 2.0_dp
  end do
end do

Part C — The No-Aliasing Advantage and pure/elemental ⭐⭐

27.11 † Why can a Fortran compiler vectorize c = a + b (on distinct assumed-shape arrays) with no runtime check, while a C compiler compiling the same loop over pointers may have to insert one — or serialize entirely? What C keyword narrows the gap, and why is it weaker than Fortran's guarantee?

27.12 Which of the following procedures can legally be declared pure? For each, say yes or no and why. (a) a function returning x*x + 1; (b) a subroutine that writes a value into an intent(out) argument; (c) a function that prints a debug line and returns a sum; (d) a function that reads and updates a module counter each call.

27.13 † Take this scalar function and make it elemental, then show the one-line call that applies it to a whole array t(:):

function to_kelvin(celsius) result(k)
  real(dp), intent(in) :: celsius
  real(dp)             :: k
  k = celsius + 273.15_dp
end function to_kelvin

27.14 In your own words, list two distinct optimizations a compiler may perform on a call to a pure function that it could not perform on an otherwise-identical impure function, and give a one-line reason for each.

Part D — Find the Bug ⭐⭐

27.15 † This program compiles and may even print a "reasonable" answer at -O0, yet it is broken. Find the defect, name the rule it violates, and explain why the printed result may change when you recompile at -O3.

pure subroutine shift_avg(a, b)          ! b(i) := average of a's two neighbors
  real(dp), intent(in)  :: a(:)
  real(dp), intent(out) :: b(:)
  integer :: i
  b = a
  do i = 2, size(a) - 1
    b(i) = 0.5_dp * (a(i-1) + a(i+1))
  end do
end subroutine shift_avg
! ...called as:
call shift_avg(v, v)

27.16 This loop will not vectorize no matter how you compile it. Identify the reason (name the kind of dependence), and say whether it is a fixable coding issue or an inherent property of the algorithm:

do i = 2, n
  x(i) = x(i-1) + h * f(i)          ! a running (prefix) sum
end do

27.17 † A colleague reports that a loop they expected to vectorize did not, and shows you the report line missed: not vectorized: possible dependence between data-refs. Give two distinct plausible causes and, for each, one thing they could check or change.

Part E — Port It ⭐⭐

27.18 Here is a pure-Python daxpy loop. Port it to a pure Fortran subroutine with intent on every argument, and add a one-sentence note on why the Fortran version vectorizes where the Python loop cannot.

def axpy(a, x, y):          # y <- a*x + y, elementwise
    for i in range(len(y)):
        y[i] = a * x[i] + y[i]
    return y

27.19 † Here is a NumPy five-point smoother of a 2D array. Port the interior update to Fortran using array sections (a whole-array statement), and separately as an explicit loop nest in the cache-friendly order. Then estimate the class of speedup you would expect over a pure-Python double loop doing the same thing (order of magnitude, framed as "typically…").

new[1:-1, 1:-1] = 0.25 * (u[:-2, 1:-1] + u[2:, 1:-1] + u[1:-1, :-2] + u[1:-1, 2:])

Part F — Back of the Envelope ⭐⭐⭐

Order-of-magnitude estimates. Show your reasoning; the exact number is not the point, the reasoning is.

27.20 A real(dp) array a(4000, 4000). (a) What is its total size in bytes, and in mebibytes? (b) A typical L1 data cache is 32 KB. Does one row (4000 elements) fit in L1? Does one column? (c) What does your answer imply about which loop order will keep the most reuse in cache?

27.21 † Suppose vectorization performs 4 real(dp) operations per instruction, but your loop is memory-bound — limited by how fast the array streams from memory, not by arithmetic. Will turning on vectorization give you a 4× speedup? Explain what actually bounds the loop, and what would have to be true for the full 4× to appear.

27.22 Estimate the arithmetic intensity (floating-point operations performed per byte of data moved) of one five-point-stencil update u(i,j) = u(i,j) + r*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)-4*u(i,j)). Count the flops and the real(dp) values read/written per interior point (assume no caching for the estimate), then argue whether the stencil is likely memory-bound or compute-bound.

27.23 † Your solver runs in 5 hours. You discover its hot loop is in the wrong (against-the-grain) order and fixing it — a two-line change taking 30 minutes — gives a 10× speedup on that loop, which is 80% of the runtime. Estimate the new total runtime (use Amdahl-style reasoning: only 80% speeds up 10×). After how many runs has the 30-minute fix paid for itself?

Part G — Design It and Read the Report ⭐⭐ / ⭐⭐⭐

27.24 Extend the Project Checkpoint: add a third variant of the stencil sweep written as a single whole-array statement using array sections (unew(2:nx-1,2:ny-1) = u(2:nx-1,2:ny-1) + r*( … )). Predict how its speed should compare to the explicit i-inner loop, and explain your prediction in terms of what the compiler sees.

27.25 † Design a minimal experiment, using system_clock, that demonstrates the -O0-versus--O3 difference on a single vectorizable loop. Say (a) what your program computes, (b) what it prints, (c) why you must print a checksum (or otherwise consume the result), and (d) how you would compile it twice to make the comparison.

27.26 On Compiler Explorer (godbolt.org) you compile a c = a + b kernel at -O3 -march=native and see instructions named vaddpd. At -O2 you see addsd. In one or two sentences, explain what each name tells you about whether the loop vectorized.

Part H — Interleaved ⭐⭐

Reaching back to earlier chapters.

27.27 † (Chapter 5 + 27.) The Chapter 5 whole-array Laplacian is written with array sections: lap(2:n-1,2:n-1) = u(1:n-2,2:n-1) + u(3:n,2:n-1) + u(2:n-1,1:n-2) + u(2:n-1,3:n) - 4*u(2:n-1,2:n-1). In what memory order does evaluating this statement traverse u, and is that traversal cache-friendly? How does writing it as sections (versus a hand loop) affect what the compiler may assume about aliasing?

27.28 (Chapter 20 + 27.) You rebuild your validated solver with -Ofast and the final printed temperature changes in its sixth decimal place. Is this a bug? Explain what -Ofast did, connect it to Chapter 20's point that floating-point addition is not associative, and say how you would decide whether the change is acceptable.

Solutions to the daggered and odd-numbered problems are in appendices/answers-to-selected.md; the compilable ones (27.8, 27.13, 27.18, 27.19, 27.24) are in code/exercise-solutions.f90. Every "expected output" in those solutions was computed by hand — compile them and check, and if your machine disagrees, you have found something worth reporting.