Case Study 2: Designing a Kernel the Compiler Will Vectorize

"The fastest code is the code you never had to hand-tune, because you handed the compiler a problem it could solve."

Executive Summary

Where the first case study diagnosed someone else's slow loop, this one asks you to build a fast one from the start — and, just as importantly, to build the small apparatus that proves it is fast. You will design a numerical kernel and a minimal "does-it-vectorize?" harness so that every principle of this chapter is working together on purpose: distinct arrays so the no-aliasing advantage applies, pure/elemental helpers so the optimizer is free, the first-index-inner loop order so memory streams with the grain, and the optimization report as the acceptance test. Then you will reason about the ceiling — the arithmetic intensity that decides whether vectorization can help at all — so you know when this effort pays and when the loop is memory-bound no matter what. This is the design mindset Chapter 29 builds into a full discipline.

Skills applied: structuring code for auto-vectorization (§27.1, §27.4); guaranteeing non-aliasing by construction (§27.3); cache-friendly loop order (§27.2); reading -fopt-info as an acceptance test (§27.5); arithmetic-intensity / roofline reasoning to predict the ceiling; honest timing (measured, not asserted).

Background

Your task is a kernel used in every simulation and benchmark suite: the vector triad, $c_i = a_i + s\, b_i$, the STREAM benchmark's namesake and the simplest kernel that stresses both arithmetic and memory. It is deliberately humble, because the goal is not the algorithm — it is the design method. If you can make the compiler vectorize a triad and can prove that it did, you can do it for the stencils and reductions that matter, because the discipline is identical.

We set four design goals up front, each traceable to a section of this chapter:

Phase 1 — Write the Kernel to Vectorize

The kernel, written to hit all four goals at once. Note what is absent: no pointers, no aliasing, no side effects, no exotic control flow — nothing for the compiler to be afraid of.

module triad_mod
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
contains
  ! c = a + s*b.  Distinct arrays (no aliasing), pure (no side effects),
  ! whole-array (structure handed to the compiler). Practically begs to vectorize.
  pure subroutine triad(a, b, s, c)
    real(dp), intent(in)  :: a(:), b(:), s
    real(dp), intent(out) :: c(:)
    c = a + s * b
  end subroutine triad
end module triad_mod

Verify correctness on values you can check in your head before you care about speed:

program check_triad
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use triad_mod, only: triad
  implicit none
  real(dp) :: a(4) = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp]
  real(dp) :: b(4) = [1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp]
  real(dp) :: c(4)
  call triad(a, b, 3.0_dp, c)        ! c = a + 3*b = [4, 5, 6, 7]
  print '(a, 4f6.1)', 'triad c = a + 3*b : ', c
end program check_triad

$ gfortran -std=f2018 -Wall -O3 -march=native check_triad.f90 -o check && ./check
triad c = a + 3*b :    4.0   5.0   6.0   7.0

The hand check: c = [1+3·1, 2+3·1, 3+3·1, 4+3·1] = [4, 5, 6, 7]. Correctness first; now we make the same kernel prove its speed.

Phase 2 — Build the "Does It Vectorize?" Harness

A benchmark you cannot trust is worse than none. The minimal honest harness times the kernel over a large array, runs it several times, and — crucially — consumes the result so the compiler cannot delete the work as dead code. (The rigorous version — warm-up passes, discarding the first run, reporting variance — is Chapter 28; here we build the honest skeleton.)

program bench_triad
  use, intrinsic :: iso_fortran_env, only: dp => real64, int64
  use triad_mod, only: triad
  implicit none
  integer,  parameter :: n = 2**20              ! 1048576 (~1e6), fits a default integer
  integer,  parameter :: reps = 50
  real(dp), allocatable :: a(:), b(:), c(:)
  integer(int64) :: t0, t1, rate
  real(dp) :: checksum
  integer :: k

  allocate(a(n), b(n), c(n))
  a = 1.0_dp;  b = 2.0_dp                        ! c should be 1 + 3*2 = 7 everywhere
  call system_clock(count_rate=rate)

  call system_clock(t0)
  do k = 1, reps
    call triad(a, b, 3.0_dp, c)
  end do
  call system_clock(t1)

  checksum = c(1) + c(n)                          ! consume the result -> no dead-code elimination
  print '(a, f0.1)',    'checksum (should be 14.0) = ', checksum
  print '(a, f10.4, a)', 'elapsed (illustrative)    = ', real(t1 - t0, dp)/real(rate, dp), ' s'
end program bench_triad

The checksum is exact and hand-computable — every element of c is 1 + 3·2 = 7, so c(1) + c(n) = 14.0 — which both guards against dead-code elimination and gives us a correctness anchor. The elapsed time is illustrative and is not shown as a number: this book runs nothing.

$ gfortran -std=f2018 -Wall -O3 -march=native bench_triad.f90 -o bench && ./bench
checksum (should be 14.0) = 14.0
elapsed (illustrative)    =     <hardware-dependent>  s

Phase 3 — Make the Compiler Report Its Verdict

Design without verification is hope. Turn the optimization report into an acceptance test: the kernel is "done" only when the compiler confirms the loop vectorized.

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

"32 byte vectors" means four real(dp) per instruction — an AVX register's worth. If instead you see nothing (and -fopt-info-vec-missed names a reason), the kernel failed the acceptance test and you go back to the four goals to find which one you broke. Cross-check on Compiler Explorer (godbolt.org): the inner loop should be built from packed instructions such as vaddpd/vfmadd…pd, not the scalar addsd/mulsd. Seeing the packed instructions is the design succeeding, made visible.

The reasoning that matters: the four design goals are not superstition — each removes a specific reason the compiler would otherwise refuse to vectorize. Distinct arrays remove the aliasing doubt; purity removes the side-effect doubt; the counted whole-array form removes the structure doubt; and correct loop order makes the vectorized loads worth doing. Break any one and the report tells you exactly which, because the report is downstream of the same four conditions.

Phase 4 — Know the Ceiling: Arithmetic Intensity

Here is the sober part, and the mark of an engineer rather than a cargo-culter: vectorization has a ceiling, and for the triad the ceiling is low. Vectorization speeds up arithmetic; if a loop spends its time waiting on memory, four-wide arithmetic finishes its work only to wait four times as eagerly for the next data. Whether that happens is set by the arithmetic intensity — flops performed per byte moved.

Count it for the triad c = a + s*b, per element: 2 flops (one multiply, one add) and, in real(dp), 3 array touches — read a, read b, write c — at 8 bytes each, so 24 bytes. That is $2 / 24 \approx 0.083$ flops per byte: a tiny intensity. The triad is emphatically memory-bound. Its speed is set by memory bandwidth, and vectorizing the arithmetic, while correct and worth doing, cannot lift it past what the memory system can deliver.

The lesson is strategic. For low-intensity kernels (triad, stencil), the returns come from moving less memory — better loop order, cache blocking, fusing passes — the material of Chapter 29, not from squeezing more flops. For high-intensity kernels (matmul), vectorization and register tiling pay off enormously, which is why a tuned BLAS dgemm is so fast and why you should call it rather than reproduce it (Chapter 21). Knowing which regime you are in before you optimize is the roofline intuition that Chapter 38 formalizes.

Phase 5 — Generalize the Method

Lift the triad's lessons into a checklist you can apply to any kernel you write:

1. Design for non-aliasing. Separate output arrays from input arrays; never pass one array as both a written and a read argument. The no-aliasing advantage is free only if you honor the contract.
2. Keep the hot loop pure. No I/O, no global-state updates, no impure calls in the inner loop. Mark helpers pure/elemental so the compiler enforces it and is freed by it.
3. Loop in the memory grain. Inner loop over the first index, always, for column-major arrays.
4. Prefer whole-array or simple counted loops. Hand the compiler structure, not a pointer puzzle.
5. Make the compiler confirm it. -fopt-info-vec is the acceptance test; packed assembly on Godbolt is the receipt.
6. Then check the ceiling. Estimate arithmetic intensity to know whether you are memory- or compute-bound, so you optimize the thing that actually limits you.

Run that checklist and you will write kernels that are fast for the right reason — because they gave a very good compiler a problem it was designed to solve — and you will be able to prove it, which is the whole difference between engineering and hoping.

Discussion Questions

1. The triad has arithmetic intensity ≈ 0.08 flop/byte and vectorizes cleanly. Reconcile these: how can a loop be "fully vectorized" and yet see little speedup from vectorization? What is actually limiting it?
2. Goal 1 (distinct arrays) and Goal 2 (purity) both remove a "doubt" the compiler would otherwise have. Describe the specific doubt each removes, and what the compiler is forced to do while the doubt remains.
3. Your harness prints a checksum "to prevent dead-code elimination." Explain the failure mode if you omit it: what might the compiler legally do to a benchmark whose result is never used, and how would that corrupt your timing?
4. Given the intensity table, which of the three kernels most rewards -march=native (wider vectors), and which would barely notice it? Defend your ranking.

Your Turn: Extensions

• Option A. Replace the triad with the five-point stencil kernel and repeat the whole method: write it to the four goals, build the harness with an exact checksum (use the maximum-principle bound from the checkpoint — maxval stays on the hot edge), and confirm vectorization with -fopt-info-vec.
• Option B. Deliberately break one design goal at a time (alias the arrays; add a print to a helper; reverse the loop order) and record, for each, what -fopt-info-vec-missed reports. Build a small table mapping each broken goal to the compiler's complaint — a diagnostic Rosetta Stone.
• Option C. Compute the arithmetic intensity of your own solver's hot loop, place it on the intensity table above, and write a one-paragraph prediction of whether loop-order/vectorization or memory-reduction/blocking will help it more — then, after Chapter 29, check your prediction against a measurement.

Key Takeaways

• Fast kernels are designed, not discovered: distinct arrays, pure hot loops, first-index-inner order, and simple structure give the compiler everything it needs to vectorize.
• A benchmark must consume its result (a checksum) or the optimizer may delete the very work you are timing.
• Vectorization has a ceiling set by arithmetic intensity. Low-intensity kernels (triad, stencil) are memory-bound — optimize memory traffic, not flops; high-intensity kernels (matmul) reward vectorization, and the tuned library rewards it most.
• The optimization report is the acceptance test that turns a design intention into a verified fact.