Chapter 27: Why Fortran Is Fast — Compiler Optimization, Memory Layout, and the No-Aliasing Advantage

"Anyone can build a fast CPU. The trick is to build a fast system." — Seymour Cray

Overview

Everything in this book so far has quietly promised that Fortran is fast. In Chapter 1 we said the reason traces to two design decisions from the 1950s — first-class arrays and a rule against pointer aliasing. In Chapter 5 we drew the memory layout and said, without proof, that getting a loop's order wrong could cost you a factor of ten. In Chapter 6 we marked procedures pure and elemental and hinted that this "helps the optimizer." We have been writing IOUs. This chapter pays them all.

Here is the thing to understand before we begin: you are not the one who makes Fortran fast — the compiler is. You write c = a + b, readable and close to the mathematics, and somewhere between your source and the running program a remarkable piece of software rewrites that line into instructions that add eight numbers at once, streaming through memory at the speed of the hardware. Your job is not to outsmart the compiler. Your job is to understand what it is trying to do, and then to hand it code it can optimize — the right loop order, the right promises, the right structure — and to read its reports when it tells you where it succeeded and where it gave up. That partnership, between a programmer who understands the machine and a compiler that has spent sixty years learning to exploit Fortran's guarantees, is the whole of high-performance Fortran. Performance is not accidental, and after this chapter it will not be mysterious either.

In this chapter, you will learn to:

• Explain what an optimizing compiler actually does — the optimization levels -O0 through -O3, inlining, loop unrolling, and vectorization — and say which of these gfortran does for you by default and which you must ask for.
• Trace the column-major loop-order effect from Chapter 5 to its cause in the CPU's cache line, and predict which of two loop nests will fly and which will crawl.
• State the no-aliasing advantage precisely: what the Fortran standard promises the compiler about procedure arguments, why C could not make the same promise, and what C's restrict keyword is trying to recover.
• Say exactly why pure and elemental are optimization licenses — promises that unlock reorderings the compiler could not otherwise risk.
• Read a compiler optimization report and look at the generated assembly, so you can see whether a loop vectorized instead of guessing.

A promise about honesty, made once and kept throughout. gfortran is not installed on the machine this book was written on, and no timing in this chapter was measured — every speedup figure is an illustrative order of magnitude, framed as "typical; measure it yourself," never a benchmark we ran. Program outputs — the actual numbers a correct program prints — are still computed by hand and exact. When you see "10×," read it as "several-fold, and you will confirm the real factor on your own hardware." The direction of every effect is certain; the magnitude is yours to measure, and Chapter 28 will teach you to measure it properly.

Learning Paths

How to read this chapter by track. - 🔬 Scientist ("my code is finally correct — now make it fast") — §27.2 (loop order) is the single highest-return page in the book for you, and §27.4 (pure/elemental) is nearly free speed. Skim §27.5 until you need it. - 📖 Standard — read straight through. §27.3 (the no-aliasing advantage) is the language-design heart of the chapter and the deepest "why Fortran" answer in the whole book. - 🔧 Legacy ("why was the old code already fast?") — §27.3 explains it: FORTRAN's no-alias rule predates C's restrict by four decades. §27.1 gives you the flags to rebuild an inherited code with optimization on. - ⚡ HPC — all of it, and it is foundational. §27.2 and §27.5 are daily tools; everything here is the groundwork for the optimization discipline of Chapter 29 and the flags of Chapter 30.

27.1 What the Compiler Actually Does

Start by dismantling a misconception. When people say "Fortran is fast," they sometimes imagine Fortran is somehow closer to the hardware than other languages — a low-level language, like assembly with better manners. This is exactly backwards. Fortran is a high-level language, further from the machine in many ways than C. Its speed does not come from the programmer touching the hardware. It comes from the programmer not touching the hardware, and instead handing a very good compiler a very clear description of the computation — one the compiler is then free to reshape aggressively, because Fortran's rules guarantee the reshaping is safe.

So the first thing to understand about performance is the optimizing compiler: the program that translates your source into machine code and, along the way, rewrites it — sometimes drastically — into a faster form that computes the same result. That last clause is the contract. Every optimization the compiler performs must be meaning-preserving: the fast program must produce the answer the slow program would have. The compiler's freedom is exactly the set of rewrites it can prove are safe, and — this is the theme of the whole chapter — Fortran's design hands it more provably-safe rewrites than most languages can.

💡 Intuition: think of the compiler as a tireless, literal-minded assistant who will rewrite your arithmetic into any equivalent form that runs faster, but who will never make a change unless it can prove the answer stays the same. Your job is to write code whose safety is easy to prove. A clear loop over a distinct array is easy; a tangle of pointers that might overlap is hard. Fortran, by design, makes your code easy to prove things about.

Optimization levels: the dial from -O0 to -O3

You control how hard the compiler works with a single flag, first met in Chapter 2:

The jump that matters most for numerical work is from -O2 to -O3, because that is where vectorization switches on. The jump you must respect is to -Ofast, because it quietly breaks a promise this book has been keeping since Chapter 20: that floating-point arithmetic follows IEEE 754 exactly. -ffast-math allows the compiler to pretend addition is associative — to compute (a + b) + c as a + (b + c) — which is not true of finite-precision floats and can move your answer in the last digits or, in a badly conditioned computation, much further. We meet that trade honestly in Chapter 30; for now, know that -O2 and -O3 preserve your arithmetic and -Ofast does not.

Compile and run. The commands you will use in this chapter, in order of aggressiveness:

console $ gfortran -std=f2018 -Wall -O2 kernel.f90 -o kernel # safe, fast, IEEE-clean $ gfortran -std=f2018 -Wall -O3 -march=native kernel.f90 -o kernel # vectorized for THIS cpu

-march=native tells the compiler to use every instruction your specific processor supports — wider vectors, newer instructions — at the cost of portability of the binary. The flags are the subject of Chapter 30; we use -O2/-O3 here and explain the rest there.

Three transformations worth naming

Under those flags, the compiler performs hundreds of distinct optimizations. Three are worth naming now, because they recur through the rest of Part VII.

Inlining. A procedure call is not free: arguments are arranged, control jumps elsewhere, a stack frame is built, and on return it is all torn down. For a small, hot procedure called millions of times, that overhead can dominate the actual work.

Definition (inlining). Inlining is the optimization that replaces a call to a procedure with a copy of the procedure's body, pasted directly into the caller. It removes the call overhead entirely, and — far more importantly — it lets the compiler optimize the pasted-in code together with the surrounding code: constants flow in, dead branches vanish, and neighboring operations combine. gfortran inlines small procedures at -O2 and is more aggressive at -O3. Inlining is often the enabling optimization: it is what turns a pure function call inside a loop (§27.4) into code the compiler can vectorize.

Loop unrolling. A loop has bookkeeping — increment the counter, test the bound, branch back — paid once per iteration. Loop unrolling reduces that tax by doing several iterations' worth of work per pass, cutting the number of tests and branches and exposing more independent operations the processor can run at once. Here is a subtlety worth knowing, and a good example of the difference between what a compiler can do and what it does by default: gfortran does not unroll loops at any -O level unless you ask. You request it explicitly with -funroll-loops (or it is enabled by profile-guided optimization). Some other compilers unroll more eagerly. This is exactly the kind of fact the optimization reports of §27.5 let you stop guessing about.

Vectorization. This is the big one, and the reason -O3 matters. A modern CPU has instructions that operate on several numbers simultaneously — one instruction that adds four (or eight) pairs of doubles in the time a scalar instruction adds one pair. These are called SIMD instructions, for Single Instruction, Multiple Data.

Definition (vectorization). Vectorization is the optimization that transforms a loop performing one operation per iteration into a loop performing that operation on several array elements at once, using the CPU's SIMD (Single Instruction, Multiple Data) instructions. A loop that adds a(i) + b(i) one element at a time becomes a loop that adds four or eight elements per instruction. The width depends on the hardware (and the -march flag); the speedup is bounded by that width, so vectorization alone can make a memory-friendly numerical loop several times faster. This is an introduction — the discipline of writing code the compiler will actually vectorize (alignment, loop structure, cache blocking, do concurrent) is the subject of Chapter 29. Here we establish what vectorization is and the two things Fortran does to make it possible: whole-array structure (§27.2) and the no-aliasing guarantee (§27.3).

Vectorization is not magic and it is not guaranteed. The compiler will vectorize a loop only when it can prove the iterations are independent — that computing element i does not depend on the result of element i-1, and that the arrays being read and written do not overlap. Those two conditions are precisely where Fortran's design pays off, and they organize the next two sections: §27.2 is about the memory access pattern that makes vectorized loads worthwhile, and §27.3 is about the aliasing guarantee that makes vectorization legal in the first place.

🧩 Try It Yourself. Open Compiler Explorer at godbolt.org, choose a recent gfortran from the compiler dropdown, and paste a tiny subroutine that does c = a + b on assumed-shape arrays. Compile first with -O2, then with -O3 -march=native, and watch the generated assembly on the right change. You are looking for instruction names ending in pd (packed double — vectorized, several numbers at once) replacing ones ending in sd (scalar double — one at a time). You do not need to read assembly fluently; you need only to see the packed instructions appear. That visual is worth more than any paragraph here.

🔄 Check Your Understanding. 1. Which optimization level first turns on auto-vectorization in gfortran, and what does -Ofast add on top of -O3 that you must be careful about? 2. What is inlining, and why is it often the optimization that enables others? 3. True or false: at -O3, gfortran unrolls your loops by default.

Answers

1. -O3 turns on auto-vectorization. -Ofast adds -ffast-math, which lets the compiler reorder floating-point arithmetic (treating it as associative) and assume no NaN/Inf — it is faster but can change your numerical results, breaking the IEEE-754 guarantees of Chapter 20. 2. Inlining pastes a called procedure's body into the caller, removing call overhead and letting the compiler optimize the body together with the surrounding code (constant propagation, vectorization across the former call boundary). It frequently unlocks optimizations that the call would otherwise have blocked. 3. False. gfortran unrolls only when you ask with -funroll-loops (or under profile-guided optimization); no plain -O level unrolls loops for you.

27.2 Column-Major Access Patterns: Loop Order as a 10× Difference

This is the section Chapter 5 has been pointing at since we first drew the memory diagram. You already know the rule — Fortran stores a 2D array one column at a time; the first index varies fastest; put the inner loop over the first index — and you have applied it in every checkpoint since. Now we explain why it is worth a factor of several to ten, and we make the explanation precise enough that you can predict the effect for any loop nest you write.

The cache line is the unit of memory traffic

The whole story is one hardware fact. When your program reads a single number from main memory, the processor does not fetch that one number. It fetches a whole contiguous block — a cache line — and holds it in a small, fast on-chip memory (the cache) on the bet that you will soon want the neighbors.

Definition (cache line). A cache line is the fixed-size contiguous block of memory — on typical hardware 64 bytes, which is eight real(dp) values — that the processor transfers between main memory and its cache as a single unit. You never load one real(dp); you load the 64-byte line containing it. Code that then uses the other seven values before the line is evicted gets them essentially for free; code that uses one value and moves on has paid for eight and used one. The gap between main-memory latency (hundreds of CPU cycles) and cache latency (a few cycles) is so large that, for most array code, the bottleneck is not arithmetic — it is memory traffic. Such code is called memory-bound.

That last word is the key to all of Part VII. A dense numerical kernel typically does a few arithmetic operations per array element — an add, a multiply — each taking roughly one cycle, while waiting for that element to arrive from memory can cost hundreds of cycles if it is not already in cache. The processor spends most of its time not computing but waiting. So the way to go fast is not to do less arithmetic; it is to waste less memory bandwidth — to make sure every cache line you pay to fetch is used to the last byte before it is thrown away. And whether that happens is decided entirely by the order in which your loops walk the array.

The two loop nests, and why one wins

Recall the layout from Chapter 5, now annotated with cache lines:

Fortran column-major layout of a(i,j), and how it fills cache lines
(8 real(dp) per 64-byte line shown as [......])

  memory:  a(1,1) a(2,1) a(3,1) ... a(8,1) | a(1,2) a(2,2) ...
           \___________ one cache line ___________/ \__ next line __
                  (all of column 1, contiguous)        (column 2 ...)

  Going DOWN a column  (inner loop over i):   ....->  walks WITHIN a line
  Going ACROSS a row   (inner loop over j):   jumps one whole column per step
                                              -> a different cache line every access

Consider summing an n × n array. The two loop nests below compute the identical number; the compiler is free to vectorize either; and one can be many times faster than the other:

! WITH the grain: inner loop over the FIRST index. Consecutive iterations
! touch consecutive memory -> each cache line is used fully before the next.
total = 0.0_dp
do j = 1, n            ! outer: over columns
  do i = 1, n          ! inner: over rows  <-- first index varies fastest
    total = total + a(i, j)
  end do
end do

! AGAINST the grain: inner loop over the LAST index. Consecutive iterations
! jump a whole column (n elements) apart -> a new cache line almost every time.
total = 0.0_dp
do i = 1, n            ! outer: over rows
  do j = 1, n          ! inner: over columns  <-- last index varies fastest
    total = total + a(i, j)
  end do
end do

Trace the memory traffic. In the first nest, the inner loop reads a(1,j), a(2,j), a(3,j), … — adjacent addresses. The first read of a column fetches a cache line of eight values; the next seven reads are already in cache, free. The processor streams down the column at full speed. In the second nest, the inner loop reads a(i,1), a(i,2), a(i,3), … — each a full column's length (n values, 8n bytes) apart. For a large n, each read lands on a different cache line, so the processor fetches 64 bytes, uses 8, and throws the rest away — then does it again. It is doing the same additions but moving eight times as much memory, and when the array is larger than the cache, it stalls on main memory almost every iteration.

Same array. Same arithmetic. Same answer. The only difference is which do line is on the inside — and for a large array the against-the-grain version can run several times, and in the worst case around ten times, slower.

⚡ Performance Note. The "10×" is an illustrative order of magnitude, not a promise, and not a number we measured. The true factor depends on n, the cache sizes, the element type, and the compiler — for a small array that fits in cache it may be a modest difference; for one much larger than the last-level cache it can exceed 10×. What is not negotiable is the direction: first-index-inner is with the grain, last-index-inner is against it, always. The Project Checkpoint has you set up exactly this measurement; Chapter 28 makes the timing rigorous.

🚪 Threshold Concept — arithmetic is cheap; memory is expensive. The instinct from a first programming course is that "fast code does less arithmetic." For dense array computation on modern hardware, that instinct is wrong, and unlearning it is a threshold you cross once. The floating-point units are so fast, and main memory so comparatively slow, that most numerical loops spend their time waiting for data, not computing. Once you see a loop nest as a pattern of memory accesses first and a pattern of arithmetic second — once "is this cache-friendly?" becomes the question you ask before "how many flops is this?" — you have started to think the way high-performance code demands. Arrays are Fortran's superpower precisely because the language's column-major layout is simple enough to hold in your head and reason about, every time.

The compiler cannot save you here

You might hope the compiler would notice a bad loop order and swap the loops for you — the optimization is called loop interchange, and compilers do sometimes perform it. But you cannot rely on it. The compiler will interchange loops only when it can prove the swap is safe and profitable, and many real loop bodies — with function calls, data dependencies, or complex indexing — defeat that analysis. Loop order is one of the few performance decisions that is squarely your responsibility, and it is the highest-leverage one you make. This is why we drilled it in Chapter 5 and drill it again here: get the loop order right by reflex, and you will never leave this particular factor of ten on the table.

🐍 Python Comparison. NumPy hides column-major-versus-row-major behind its array object, and defaults to row-major (C order) — the opposite of Fortran. A NumPy sum like A.sum() is fast because it walks memory in whatever order A is actually stored, inside a compiled loop. But the moment you pass arrays between Python and Fortran, the layout mismatch becomes yours to manage: a NumPy array handed to a Fortran routine must be in Fortran (column-major) order, or you get a silent transpose or a slow copy. This is the order='F' issue you will handle directly in Chapter 15. And when you drop out of NumPy's vectorized calls into an explicit Python for loop over elements — as a stencil sweep forces you to — you lose not only vectorization but the memory-layout control entirely, and the pure-Python loop runs commonly 50–100× slower than the Fortran. Fortran and Python are better together: keep the memory-order-sensitive inner loop in Fortran.

27.3 The No-Aliasing Advantage

We come to the deepest answer this book gives to the question "why Fortran?" — the one we sketched informally in Chapter 1 and promised to make rigorous here. It is not about arithmetic and not about memory layout. It is about a promise the language makes to the compiler that most languages cannot make, and it is the reason a Fortran loop and a C loop that look identical can compile to code of very different speed.

First, recall the hazard, defined back in Chapter 1.

🔗 Connection — aliasing, from Chapter 1. Two names alias when they refer to the same storage. If a compiler cannot rule out that x and y might be the same memory, it must assume that writing through x could change y — which forbids a great many reorderings and optimizations. We met this informally in Chapter 1 §1.3; now we make it the center of the performance story.

What the Fortran standard promises

Here is the rule, in plain words first. When you pass arguments into a Fortran procedure, and the procedure writes to one of them, the language forbids that written argument from secretly being the same storage as another argument the procedure reads. The programmer is responsible for never creating such an overlap; in exchange, the compiler is entitled to assume the overlap never happens and to optimize accordingly.

Definition (the no-aliasing advantage). The no-aliasing advantage is the optimization freedom a Fortran compiler gains from the standard's prohibition on argument aliasing: if a dummy argument is defined (written) within a procedure, the program is not permitted to have it associated with — aliased to — any other dummy argument or accessible entity that the procedure also references. Because the standard makes this the programmer's responsibility, the compiler may assume distinctness for free, and is therefore free to load, compute, reorder, and store in whatever sequence is fastest — including vectorizing. It is a language-level guarantee, enforced by the contract of the standard rather than checked at run time. Violating it does not produce a diagnostic; it produces undefined behavior, so the discipline is: never alias a written argument.

The consequence is concrete. Consider the humble operation from Chapter 1 — scale one array and add it to another, the kernel the numerical-linear-algebra world calls axpy ($\mathbf{y} \leftarrow a\mathbf{x} + \mathbf{y}$). In Fortran:

pure subroutine axpy(a, x, y)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp), intent(in)    :: a, x(:)
  real(dp), intent(inout) :: y(:)
  y = y + a * x            ! whole-array; compiler KNOWS x and y are distinct
end subroutine axpy

The standard guarantees x and y do not overlap. The compiler therefore knows that reading x(i) can never be affected by writing y(i-1), so it can load a whole SIMD register's worth of x, another of y, fuse-multiply-add them, and store the result — no iteration waiting on the one before. The loop vectorizes, and it vectorizes without a single runtime check.

Why C could not, and what restrict recovers

Now the same routine in C:

void axpy(int n, double a, const double *x, double *y) {
    for (int i = 0; i < n; i++)
        y[i] = y[i] + a * x[i];
}

x and y are pointers, and in C, by default, any two pointers might point into the same array. The compiler must assume the worst: that y and x could overlap, so that writing y[i] might change a value of x the loop is about to read. Historically — in C89, before any fix existed — this simply forbade vectorization: the compiler had to execute the loop in strict order, one element at a time, serialized. This is the origin of the decades-old folklore that "Fortran is faster than C for numerics." It was not folklore; it was this.

Modern C compilers are cleverer than that: rather than give up, they often emit a runtime aliasing check — a few instructions at the top of the loop that test whether the arrays actually overlap — and generate two versions of the loop, a fast vectorized one taken when they do not overlap and a safe scalar one when they do. That recovers much of the speed, but at a cost Fortran never pays: the check itself, the branch, the duplicated code (which bloats the instruction cache), and — crucially — the fact that the analysis often cannot be done at all across a function-call boundary or through several layers of pointers, where the compiler falls back to the conservative, serialized code.

So C99 added a keyword to let the programmer opt back into what Fortran guarantees by default:

void axpy(int n, double a, const double *restrict x, double *restrict y) {
    for (int i = 0; i < n; i++)
        y[i] = y[i] + a * x[i];   /* restrict PROMISES x and y do not overlap */
}

Definition (restrict). restrict is a C99 type qualifier by which the C programmer promises the compiler that, for the lifetime of the pointer, the object it points to is accessed only through that pointer — i.e., it does not alias. It is the manual, opt-in, per-pointer, easy-to-get-wrong equivalent of the guarantee Fortran makes automatically, by default, for every procedure argument in the language. If the C programmer lies — passes overlapping pointers marked restrict — the behavior is undefined, exactly as in Fortran, but with no language-level habit of correctness to lean on.

There is the difference, stated exactly. In C, non-aliasing is a promise the programmer must remember to make, pointer by pointer, and must not get wrong. In Fortran, non-aliasing is the default state of the world, guaranteed by the standard for every argument, so the compiler optimizes aggressively from the first line without being asked. Fortran programmers get, for free and everywhere, the optimization C programmers must request explicitly and can silently break.

🚪 Threshold Concept — the compiler optimizes because Fortran forbids aliasing. This is the sentence to carry out of the chapter. Fortran's celebrated speed for numerical work is not because it is low-level, not because its arithmetic is special, and not only because of column-major arrays. It is, more than anything else, because the language makes a promise about memory — that the things a procedure writes are distinct from the things it reads — and that promise is exactly the piece of information a compiler needs to reorder and vectorize freely. The language gives up a little flexibility (you may not alias arguments) and buys, in return, an enormous amount of optimization headroom. Once you see that the speed is the guarantee, you understand why sixty years of "surely we can replace Fortran now" has kept arriving at the same answer. This is why Fortran is not dead: the design decision that makes it fast is baked so deep that matching it means rebuilding the same promise, which is what C had to do with restrict and what every serious numerical language wrestles with.

🐛 Find the Bug. This subroutine averages each interior element's two neighbors — a tiny 1D smoother. It is called two ways. One call is a lurking disaster. Which, and why does it matter more at -O3 than at -O0?

fortran pure subroutine smooth(a, b) ! b(i) := average of a's neighbors use, intrinsic :: iso_fortran_env, only: dp => real64 implicit none 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 smooth ! ... call smooth(x, y) ! call 1 call smooth(x, x) ! call 2 <-- ?

Diagnosis

Call 2 is nonconforming. b is written (intent(out)) and a is read (intent(in)); passing the same array x for both aliases a written argument with a read one, which the Fortran standard forbids. The result is undefined behavior. Why it bites harder at -O3: at -O0 the compiler runs the loop literally, element by element, so b(i) reads whatever a(i-1) currently holds — and since b is a, a(i-1) was just overwritten in the previous iteration, giving a "cascading" result that at least is deterministic. At -O3, trusting the no-alias guarantee, the compiler may vectorize — loading a block of the original a values at once, before any writes — producing a different answer. Neither answer is "the bug"; the bug is that the program was never valid, and the two optimization levels disagreeing is the symptom. Fix: never alias a written argument — use distinct arrays (call smooth(x, y)), and if you genuinely need in-place behavior, write it explicitly with a temporary. This is the no-aliasing rule working exactly as designed: your side of the contract is to keep written arguments distinct.

📜 From History. Fortran's anti-aliasing rule is as old as procedure arguments themselves; it was in the language long before anyone spoke of "optimizing compilers" as a discipline. It was not adopted for speed — it came from Fortran's argument-passing model — but it turned out to be one of the most valuable optimization guarantees a language ever accidentally made. C, designed for systems programming where pointers must be free to alias, made the opposite choice for good reasons of its own, and spent a quarter century (until restrict in C99) building back a way to say what Fortran had always said. The lesson is one the whole book keeps teaching: Fortran's "old" design decisions were, for numerical computing, quietly and durably right.

🔄 Check Your Understanding. 1. State the no-aliasing advantage in one sentence, from the compiler's point of view. 2. What does C's restrict keyword do, and why is it needed when Fortran needs no equivalent? 3. Is aliasing a written Fortran dummy argument a compile-time error? What happens if you do it?

Answers

1. The Fortran standard guarantees that a procedure's written arguments do not overlap its read arguments, so the compiler may assume distinctness and reorder/vectorize loads and stores freely, with no runtime check. 2. restrict lets a C programmer promise, per pointer, that it does not alias — recovering the optimization Fortran gets by default. It is needed because C pointers may alias by default; Fortran arguments may not, so Fortran needs no keyword. 3. No — it is generally not diagnosed at compile time (the compiler cannot see how you will call the procedure). It is undefined behavior: the program may give different answers at different optimization levels, or appear to work and fail later. The programmer must keep written arguments distinct.

27.4 pure and elemental: Optimization Licenses

Back in Chapter 6 §6.4 you marked procedures pure and elemental and we promised the payoff would come "in Chapter 27." Here it is. These attributes are not decoration and not merely documentation: they are licenses you grant the compiler — promises about your code's behavior that unlock optimizations the compiler could not otherwise justify.

Recall the meaning (defined in Chapter 6; we only use it here):

🔗 Connection — pure and elemental, from Chapter 6. A pure procedure promises it has no side effects: it does not modify global or host state, perform I/O, or do anything observable but return its result through its arguments. An elemental procedure is written for scalar arguments but may be called on whole arrays, applying elementwise, and is automatically pure. The compiler enforces both — try to print inside a pure procedure and the build fails. See Chapter 6 §6.4 for the full rules; this section is about what the compiler does with those promises.

Why purity is speed

A pure function's promise — "calling me changes nothing except my result, and my result depends only on my inputs" — is exactly what an optimizer needs to hear, for four concrete reasons:

• Hoisting out of loops. If a pure function is called inside a loop with arguments that do not change across iterations, the compiler may call it once, before the loop, and reuse the result. It could never risk this for an impure function that might, say, advance a random-number generator or write a log line on each call.
• Common-subexpression elimination. Two calls to a pure function with the same arguments must return the same value, so the compiler may compute it once and share it. Impure calls must each be executed, in case they differ.
• Reordering and removal. A pure call whose result is unused can be deleted entirely; pure calls with no dependency between them can be reordered into whatever schedule runs fastest.
• Parallel and vectorized application. Because a pure call cannot interfere with any other, the compiler (and you) may run many of them at once. This is why the standard requires the body of a do concurrent loop, and any procedure it calls, to be free of side effects — the very promise pure encodes. We use do concurrent for real in Chapter 29.

Here is the mechanism in miniature. Suppose a loop transforms each element through a small function:

do i = 1, n
  b(i) = scale * transform(x(i)) + offset
end do

If transform is pure, the compiler may inline it (§27.1), see that scale and offset are loop-invariant, hoist them, and vectorize the whole loop — the pure promise plus inlining is what makes the vectorization legal, because the compiler now knows the inlined body has no hidden side effect that must happen once per iteration in order. If transform were impure — if it might touch a module variable — the compiler would have to call it, in order, once per element, and the loop would not vectorize. Same loop; the attribute is the difference.

elemental is vectorization-shaped

elemental goes one step further: it tells the compiler the function is meant to be applied independently to every element of an array. That is the exact shape of a vectorizable loop. Write your scalar-to-scalar transformation once as elemental, apply it to a whole array, and you have handed the compiler a loop it is practically begging to vectorize — with the independence guaranteed by the elemental/pure contract rather than something it must prove.

elemental function activate(x) result(y)   ! automatically pure
  use, intrinsic :: iso_fortran_env, only: dp => real64
  real(dp), intent(in) :: x
  real(dp)             :: y
  y = max(0.0_dp, x)                        ! a "ReLU"; scalar body, array-ready
end function activate
! ...
field = activate(field)                     ! whole-array; vectorizable elementwise

⚡ Performance Note. Marking your small mathematical helpers pure and your scalar-to-scalar transformations elemental — whenever they honestly qualify — is one of the cheapest speedups available: zero code change to the algorithm, a real gift to the optimizer, and better documentation besides. It is not a guarantee of vectorization (the loop body must still be vectorizable and the access pattern cache-friendly), but it removes an obstacle the compiler would otherwise trip over. Treat the attributes as habit, not afterthought.

🔗 Connection — and yet a tuned library still wins. Purity and good loop order will make your loops fast, but for the classic dense-linear-algebra kernels — matrix multiply, solving $A\mathbf{x} = \mathbf{b}$ — a hand-tuned BLAS/LAPACK will still beat your best pure loop, because it adds cache blocking, register tiling, and hand-written vector code on top of everything here. That is not a failure of your Fortran; it is decades of specialist tuning you should reuse, not reproduce. We met LAPACK in Chapter 21 and return to "why the library wins" in Chapter 29. Know which loops to optimize yourself and which to hand to dgemm.

27.5 Reading a Compiler Optimization Report

Everything so far has been about what the compiler can do. This section is about finding out what it actually did — because the alternative, guessing, is how people waste afternoons "optimizing" a loop the compiler already vectorized while ignoring the one it silently gave up on. Stop guessing. Ask.

Definition (optimization report). An optimization report is diagnostic output the compiler emits, on request, describing which optimizations it applied to which lines of your source and — often more usefully — which it declined to apply and why. It turns the optimizer from a black box into a conversation: you can see that the loop on line 40 vectorized, that the one on line 55 did not "because of a possible data dependence," and act on the specific reason.

gfortran: the -fopt-info family

gfortran reports through the -fopt-info family of flags. The workhorses:

The reports go to standard error, or to a file if you append =filename. Because vectorization only happens with optimization on, you pair these with -O3:

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

That line says the loop starting at line 14 was vectorized with 32-byte vectors — four real(dp) per instruction (an AVX register). When a loop you expected to vectorize does not, ask why:

$ gfortran -std=f2018 -O3 -fopt-info-vec-missed loop.f90 -o loop
loop.f90:22:15: missed: couldn't vectorize loop
loop.f90:23:10: missed: not vectorized: possible dependence between data-refs

That second message — "possible dependence between data-refs" — is the compiler telling you it could not prove the loop's iterations are independent. In Fortran that is often a real dependency in your algorithm (a recurrence, a(i) = a(i-1) + …); in code translated from C it can be an aliasing worry that Fortran's rules should have removed, a sign you have accidentally reintroduced a dependence. Either way, the report points you at the exact line.

⚠️ A caution about exact wording. The precise text of these messages varies between gfortran versions — "loop vectorized using N byte vectors," the phrasing of a missed reason, the line/column format — so the lines above are representative, not a fixed contract. Treat the report as a signal to read (did this loop vectorize, yes or no, and if not what reason did it name?), not a string to match. Your gfortran will phrase it its own way.

Seeing it yourself: assembly and Compiler Explorer

The report tells you what the compiler decided; the generated assembly shows you the result. You do not need to read assembly fluently — you need to recognize two families of instruction, and Compiler Explorer (godbolt.org) makes them easy to see side by side with your source:

• Scalar floating-point instructions end in sd — addsd, mulsd, vaddsd — for "scalar double," one number per instruction.
• Packed (vectorized) instructions end in pd — addpd, vaddpd, vfmadd…pd — for "packed double," several numbers per instruction.

A loop full of …pd instructions vectorized; a loop of …sd did not. That is the whole diagnostic, and it takes ten seconds once you know where to look.

🔗 Connection — other compilers, and Chapter 30. The Intel compilers (ifx/ifort) produce a report with -qopt-report (writing .optrpt files), and NVIDIA's nvfortran uses -Minfo. Same idea, different flags. Chapter 30 surveys the compiler-specific flags in full; the point common to all of them is the discipline: make the compiler tell you what it did, then optimize what it actually missed — which is exactly the measure-first ethic Chapter 28 builds into a method.

🔄 Check Your Understanding. 1. Which gfortran flag tells you why a loop failed to vectorize, and what optimization level must accompany it? 2. On Compiler Explorer, what instruction-name feature distinguishes a vectorized floating-point loop from a scalar one? 3. Why is reading the optimization report a better use of time than "optimizing by intuition"?

Answers

1. -fopt-info-vec-missed, paired with -O3 (vectorization is off at lower levels, so there would be nothing to report). 2. Vectorized loops use packed instructions ending in pd (e.g. vaddpd); scalar loops use instructions ending in sd (e.g. addsd). 3. Because the compiler often already did what you were about to do by hand (so you would gain nothing), or gave up on a loop for a specific, fixable reason (which the report names) — either way, the report tells you where the real opportunity is, so you spend effort where it pays.

Project Checkpoint

Every chapter so far has built the heat solver; this one measures it. Your checkpoint is to take the solver's stencil sweep — the interior update from Chapter 24 — write it in the two loop orders (inner loop over i, then inner loop over j), time both, and confirm with your own eyes that column-major layout costs what §27.2 says it costs.

The program below is self-contained and compilable. It does two things. First, a correctness check on a tiny 5×5 plate — the exact setup from the Chapter 24 checkpoint (hot top edge at 100, everything else 0, CFL-safe ratio $r = \alpha\,\Delta t/\Delta x^2 = 0.2$) — running one explicit step in both loop orders and proving they give the identical field. Second, a timing run on a larger grid, sweeping many times in each order. The correctness numbers are exact and hand-computed; the timings are illustrative — you will run this yourself and record your own factor.

module kinds
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)
end module kinds

module stencil_bench
  use kinds, only: dp
  implicit none
contains
  ! One explicit (FTCS) interior sweep, inner loop over the FIRST index i
  ! (down columns) -> with the grain of column-major memory.
  subroutine step_i_inner(u, unew, r)
    real(dp), intent(in)  :: u(:,:), r
    real(dp), intent(out) :: unew(:,:)
    integer :: i, j, nx, ny
    nx = size(u,1);  ny = size(u,2)
    unew = u                                   ! copy Dirichlet boundary; interior overwritten
    do j = 2, ny-1                             ! outer: columns
      do i = 2, nx-1                           ! inner: first index  <-- fast
        unew(i,j) = u(i,j) + r*( u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) - 4.0_dp*u(i,j) )
      end do
    end do
  end subroutine step_i_inner

  ! The SAME arithmetic, inner loop over the LAST index j (across rows)
  ! -> against the grain: a new cache line almost every access.
  subroutine step_j_inner(u, unew, r)
    real(dp), intent(in)  :: u(:,:), r
    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: rows
      do j = 2, ny-1                           ! inner: last index  <-- slow
        unew(i,j) = u(i,j) + r*( u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) - 4.0_dp*u(i,j) )
      end do
    end do
  end subroutine step_j_inner
end module stencil_bench

program checkpoint
  use, intrinsic :: iso_fortran_env, only: int64
  use kinds,         only: dp
  use stencil_bench, only: step_i_inner, step_j_inner
  implicit none
  real(dp), allocatable :: a(:,:), b(:,:)
  real(dp) :: fa(5,5), fi(5,5), fj(5,5)             ! distinct in/out arrays: never alias
  real(dp), parameter :: r = 0.2_dp                 ! CFL-safe (Chapter 24): r <= 1/4
  integer,  parameter :: n = 1000, nsweeps = 100    ! nsweeps even
  integer(int64) :: t0, t1, rate
  integer :: s

  ! --- correctness on a 5x5 plate: one step, both orders, must be identical ---
  fa = 0.0_dp;  fa(1,:) = 100.0_dp                  ! hot top edge (row i=1)
  call step_i_inner(fa, fi, r)                      ! fi = one step, i-inner
  call step_j_inner(fa, fj, r)                      ! fj = one step, j-inner (fa read-only both times)
  print '(a, 3f8.2)', 'i-inner, interior row 2:  ', fi(2,2:4)
  print '(a, 3f8.2)', 'j-inner, interior row 2:  ', fj(2,2:4)
  print '(a, es9.1)', 'max |difference|        = ', maxval(abs(fi - fj))

  ! --- timing on a large grid: many sweeps each order (ILLUSTRATIVE) ---
  allocate(a(n,n), b(n,n))
  call system_clock(count_rate=rate)

  a = 0.0_dp;  a(1,:) = 100.0_dp;  b = a
  call system_clock(t0)
  do s = 1, nsweeps, 2
    call step_i_inner(a, b, r);  call step_i_inner(b, a, r)   ! ping-pong buffers, no copy
  end do
  call system_clock(t1)
  print '(a, f9.2, a)',  'maxval(u) after i-inner = ', maxval(a), '   (bounded by the hot edge)'
  print '(a, f10.4, a)', 'i-inner elapsed (with the grain)    = ', real(t1-t0, dp)/real(rate, dp), ' s'

  a = 0.0_dp;  a(1,:) = 100.0_dp;  b = a
  call system_clock(t0)
  do s = 1, nsweeps, 2
    call step_j_inner(a, b, r);  call step_j_inner(b, a, r)
  end do
  call system_clock(t1)
  print '(a, f9.2, a)',  'maxval(u) after j-inner = ', maxval(a), '   (same field, same max)'
  print '(a, f10.4, a)', 'j-inner elapsed (against the grain) = ', real(t1-t0, dp)/real(rate, dp), ' s'
end program checkpoint

$ gfortran -std=f2018 -Wall -O3 -march=native project-checkpoint.f90 -o checkpoint && ./checkpoint

The exact, hand-computed part of the output — the correctness check — is:

i-inner, interior row 2:     20.00   20.00   20.00
j-inner, interior row 2:     20.00   20.00   20.00
max |difference|        =   0.0E+00
maxval(u) after i-inner =    100.00   (bounded by the hot edge)
maxval(u) after j-inner =    100.00   (same field, same max)

Two things are provably exact here, and worth the hand check. The interior values are 20.00. With $r=0.2$, the hot top row at 100, and everything else 0, one step gives each interior point of row 2 the value $0 + 0.2\,(100 + 0 + 0 + 0 - 0) = 20$ — identical to the Chapter 24 trace, and identical in both loop orders because each output element reads only old neighbors, so visitation order cannot change the result. The maximum stays 100.00 for any number of sweeps: with $r \le 1/4$ the update is a weighted average with non-negative weights summing to one (a convex combination), so no interior value can exceed the maximum of the fixed boundary — the discrete maximum principle. That maxval(u) = 100.0 is what keeps the timed loop's work from being optimized away and gives us an exact value to check.

The timing lines are deliberately not shown as numbers. They are hardware-dependent, and this book runs no code. What you will typically see when you run it: the j-inner (against-the-grain) sweep takes noticeably longer than the i-inner one — commonly on the order of a few times slower for a $1000\times1000$ grid, and the gap widens as you grow n past your cache. Run it, read the two elapsed lines, and compute the ratio yourself; that number, on your machine, is the payoff of Chapter 5 made concrete. If the two are nearly equal, your grid may be fitting in cache — increase n and try again.

This is a rough measurement, not a rigorous one: no warm-up, one trial, wall-clock via system_clock. That is exactly the gap Chapter 28 closes, where you will learn to warm the cache, repeat, and separate wall time from CPU time. And the fix — making the slow order fast, and pushing further with blocking and do concurrent — is Chapter 29. For now you have done the essential first thing: you have made the invisible visible, and confirmed that loop order alone moves the clock.

Summary

This chapter explained the mechanism behind a claim the book has made since page one: that Fortran is fast, and that the speed is a consequence of design, not luck.

The two things to memorize. First: the inner loop runs over the first index, because Fortran is column-major and the cache line is the unit of memory traffic — this single habit is the largest performance factor you personally control. Second: the compiler optimizes because Fortran forbids aliasing — the language's promise that written arguments are distinct is the deepest reason its numerical code is fast, and it is why C needed restrict to catch up. Everything in Chapter 29 is built on these two facts.

Spaced Review

Retrieval practice on the two chapters this one rests upon. Answer before peeking.

1. (Ch. 5) In memory, which element immediately follows a(1000, 1) in a real(dp) :: a(1000, 1000), and what does that imply for which loop index belongs on the inside?
   

   Answer

   `a(1, 2)` — after the last element of column 1 comes the first element of column 2, because Fortran stores the first index fastest (column-major, [Chapter 5](../../part-01-foundations/chapter-05-arrays/index.md)). So the **inner loop must run over the first index `i`**, sweeping down a column, so consecutive iterations touch adjacent memory and fill each cache line before moving on.

2. (Ch. 5) Why is the whole-array statement c = a + b not merely shorter than the equivalent loop, but better information for the compiler?

   Answer

   It states the entire operation at once over known-shape, known-distinct arrays, with no possibility of aliasing among a, b, c (Fortran's rule) — exactly the structure the compiler needs to vectorize. A hand-written loop hides that structure and forces the compiler to reconstruct it. The readable form and the fast form are the same form; that is the point of Chapter 5.

3. (Ch. 20) The checkpoint compares two arrays with maxval(abs(fb - fj)) and gets exactly zero, even though Chapter 20 warns against comparing floats for exact equality. Why is exact zero legitimate here?
   

   Answer

   Because the two loop orders perform the *identical arithmetic operations* on each element — each output reads only old neighbors, so the same additions and multiplications produce bit-for-bit identical results regardless of visitation order. Chapter 20's caution is about comparing floats produced by *different* computations (where rounding differs); here the computations are the same, so the difference is genuinely, exactly `0.0`.

4. (Ch. 20) Why can compiling with -Ofast (which turns on -ffast-math) change the result of a reduction like sum(x), and how does that connect to catastrophic cancellation?
   

   Answer

   `-ffast-math` lets the compiler treat floating-point addition as associative and reorder it — e.g., summing in a different order or in parallel partial sums. Because finite-precision addition is *not* associative ([Chapter 20](../../part-05-numerical-methods/chapter-20-floating-point/index.md)), a different order yields a different rounding error, so the sum can change in its last digits — and in an ill-conditioned sum with heavy cancellation, much more than the last digits. That is why `-Ofast` is opt-in and must be validated against your numerics.

5. (Ch. 5 + 20) You mark a small function pure to help the optimizer, and it does floating-point work. Does pure change the numerical result the function computes?

   Answer

   No. pure restricts side effects (no I/O, no global-state changes), not arithmetic; the function computes exactly the same value. What it changes is the compiler's freedom — to hoist, share, reorder, and vectorize calls — which affects speed, not the result. (Aggressive flag choices like -ffast-math can change results; the pure attribute alone does not.)

What's Next

You have seen why Fortran is fast and confirmed that loop order moves the clock — but you measured it crudely, with a single untimed guess of a run. Before you optimize anything for real, you need to measure properly: to find, with evidence rather than intuition, exactly which loop in a large program is eating the time. That is the discipline of Chapter 28: timing with cpu_time and system_clock, profiling with gprof to read a flat profile and a call graph, distinguishing memory-bound from compute-bound loops, and the benchmarking methodology — warm-up, repetition, variance — that turns a number you ran into a number you can trust. Measure first, optimize second, and never the other way around. Let's learn to measure.