Case Study 2: Optimizing a Matrix Multiply — and Learning to Call BLAS Instead

"The competent programmer is fully aware of the strictly limited size of his own skull; therefore he approaches the programming task in full humility." — Edsger W. Dijkstra

Executive Summary

Case Study 1 repaired someone else's loop. This one builds — you will write a dense matrix multiply, then optimize it up the full ladder of this chapter: correct loop order, then cache blocking, then vectorization, measuring the gain at each rung. It is the chapter's showcase because matrix multiply is the one kernel where blocking genuinely shines (unlike the memory-bound stencil), so you get to see tiling earn its complexity. And then comes the lesson the whole performance part has been building toward since Chapter 1: you will compare your best, hard-won, blocked, vectorized Fortran against a tuned BLAS dgemm (Chapter 21) — and it will beat you by an order of magnitude. The point of the exercise is not to win; it is to understand exactly why you lose, and therefore when to stop writing loops and call the library. That is the most valuable thing this chapter can teach.

Skills applied

  • Building and correctness-checking a dense matrix multiply against the matmul intrinsic (§29.1).
  • Reordering the triple loop for column-major, unit-stride access (§29.1).
  • Implementing cache blocking on a genuinely compute-bound, high-reuse kernel (§29.2).
  • Reasoning about arithmetic intensity and the roofline to know which optimizations can pay (§29.2).
  • Deciding when to stop and call a tuned BLAS instead — the §29.5 climax and the Chapter 21 anchor.

Background

Matrix multiply, $C = AB$, is the beating heart of dense numerical computing — it underlies linear solves, least squares, neural-network layers, and much of the physics in Part V. It is also the textbook example of a compute-bound, high-reuse kernel: for $n\times n$ matrices it does $O(n^3)$ arithmetic over only $O(n^2)$ data, so each element of $A$ and $B$ is reused $n$ times. That reuse is exactly what cache blocking exists to capture, and exactly what a naive triple loop throws away. We will build it up in stages, verifying at every step that the answer never changes, and use a hand-checkable test throughout: with $A(i,k)=i$ and $B(k,j)=j$, the product is $C(i,j) = \sum_k i\,j = n\,i\,j$, which for $n=4$ is $C(i,j)=4ij$.

Phase 1 — A correct baseline

Start with the most literal translation of the mathematical definition — the order you would write on a blackboard, i outermost:

subroutine matmul_naive(a, b, c, n)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  integer,  intent(in)  :: n
  real(dp), intent(in)  :: a(n,n), b(n,n)
  real(dp), intent(out) :: c(n,n)
  integer :: i, j, k
  c = 0.0_dp
  do i = 1, n                              ! i OUTERMOST: the "math" order
    do j = 1, n
      do k = 1, n
        c(i,j) = c(i,j) + a(i,k)*b(k,j)
      end do
    end do
  end do
end subroutine matmul_naive

It is correct — verify it against matmul on the $4\times4$ test — but it is slow, and for a column-major reason. With k innermost, a(i,k) strides across a whole row per step (column-major: a(i,k) and a(i,k+1) are n elements apart), so the access to A misses cache relentlessly. On a large matrix this naive order reaches only a few percent of the machine's peak arithmetic rate (illustratively; measure yours). That gap between a few percent and 100% is the entire story of this case.

Phase 2 — Fix the loop order

Reorder so the innermost loop indexes the first array index. A standard cache-friendly order for column-major matmul puts i innermost (the jki order), so both c(i,j) and a(i,k) are walked with unit stride down their columns:

subroutine matmul_ordered(a, b, c, n)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  integer,  intent(in)  :: n
  real(dp), intent(in)  :: a(n,n), b(n,n)
  real(dp), intent(out) :: c(n,n)
  integer  :: i, j, k
  real(dp) :: bkj
  c = 0.0_dp
  do j = 1, n
    do k = 1, n
      bkj = b(k,j)                         ! b(k,j) is loop-invariant in i: hoist it
      do i = 1, n                          ! i INNERMOST -> unit stride down columns
        c(i,j) = c(i,j) + a(i,k)*bkj
      end do
    end do
  end do
end subroutine matmul_ordered

Same product, delivered in unit-stride sweeps. On a large matrix this reorder alone is a large multiplier over the naive order, and the inner loop is now a clean, vectorizable axpy (c(:,j) += a(:,k)*bkj) the compiler will happily SIMD-ize. But we are still leaving reuse on the table: each column of A is re-read from memory for every column of C, and once n is large enough that a column no longer stays in cache between uses, those re-reads go all the way to main memory.

Phase 3 — Cache blocking, where it finally pays

Here is where tiling earns its keep — because, unlike the stencil, this kernel has $O(n)$ reuse to capture. Block all three loops so the multiply works on sub-tiles small enough to stay resident in cache for all of a tile's reuses:

subroutine matmul_blocked(a, b, c, n, nb)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  integer,  intent(in)  :: n, nb
  real(dp), intent(in)  :: a(n,n), b(n,n)
  real(dp), intent(out) :: c(n,n)
  integer :: i, j, k, ii, jj, kk
  c = 0.0_dp
  do jj = 1, n, nb
    do kk = 1, n, nb
      do ii = 1, n, nb                     ! outer loops step by the block size
        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 within the tile
              c(i,j) = c(i,j) + a(i,k)*b(k,j)
            end do
          end do
        end do
      end do
    end do
  end do
end subroutine matmul_blocked

Now a tile of A and a tile of B are loaded once and reused nb times before being evicted, converting main-memory re-reads into cache hits. On a large matrix, with nb tuned to the L1/L2 size (§29.2's back-of-envelope, then measured), this lifts the kernel to a substantial fraction of peak — illustratively into the tens of percent, well above the ordered version. This is the payoff §29.2 promised for a compute-bound, high-reuse kernel: blocking is worth its complexity here precisely because there is reuse to capture, which the memory-bound stencil never had.

Verify, as always, that none of this changed the answer:

program verify_matmul
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter :: n = 4
  real(dp) :: a(n,n), b(n,n), c_naive(n,n), c_block(n,n), c_ref(n,n)
  integer  :: i, j
  do j = 1, n
    do i = 1, n
      a(i,j) = real(i,dp);  b(i,j) = real(j,dp)    ! -> C(i,j) = 4ij
    end do
  end do
  call matmul_naive  (a, b, c_naive, n)
  call matmul_blocked(a, b, c_block, n, 2)
  c_ref = matmul(a, b)
  print '(a)', 'C = A*B  (blocked):'
  do i = 1, n
    print '(4f8.1)', c_block(i,:)
  end do
  print '(a, f6.1)', 'max |naive   - matmul| = ', maxval(abs(c_naive - c_ref))
  print '(a, f6.1)', 'max |blocked - matmul| = ', maxval(abs(c_block - c_ref))
  ! (the three matmul_* subroutines as above, in the same file)
end program verify_matmul

$ gfortran -std=f2018 -Wall -O3 verify_matmul.f90 -o vm && ./vm
C = A*B  (blocked):
     4.0     8.0    12.0    16.0
     8.0    16.0    24.0    32.0
    12.0    24.0    36.0    48.0
    16.0    32.0    48.0    64.0
max |naive   - matmul| =    0.0
max |blocked - matmul| =    0.0

Hand-check a cell: $C(2,3) = \sum_{k=1}^4 A(2,k)B(k,3) = \sum_{k=1}^4 2\cdot3 = 4\cdot6 = 24$, matching the printed row 2, column 3. Every order — naive, blocked, and the intrinsic — computes the identical product, so both differences are exactly $0$. (With these integer values the sums are exact in real(dp); note that for general reals, reordering the $k$-accumulation can change the last bit, since floating-point addition is not associative — a <= tol check would be the honest assertion there, per §29.3.)

Phase 4 — The reckoning: dgemm beats you

You now have a blocked, vectorized matrix multiply that hits a respectable fraction of peak, and you are proud of it. Link against a tuned BLAS and call dgemm on the same matrices (Chapter 21):

! C := 1.0*A*B + 0.0*C  via the BLAS Level-3 routine
call dgemm('N', 'N', n, n, n, 1.0_dp, a, n, b, n, 0.0_dp, c, n)

$ gfortran -std=f2018 -Wall -O3 mm_bench.f90 -o mmb -lopenblas && ./mmb

On a large matrix, OpenBLAS or Intel MKL dgemm reaches roughly 80% of the machine's peak — and your best blocked loop, perhaps 10–20%. That is an order-of-magnitude gap, and no reasonable amount of further hand-tuning closes it. Why does the library win so decisively? Three reasons, each a technique from this chapter taken to a level you should not attempt by hand:

What dgemm does This chapter's version
Blocks for every cache level at once (registers, L1, L2, L3), each tile size tuned per CPU model We blocked for one level, with one hand-guessed nb
Inner microkernel hand-written in assembly with explicit SIMD and software pipelining to hide latency We relied on the compiler's auto-vectorizer
Packs tiles into contiguous buffers to guarantee unit stride and TLB friendliness We passed the arrays as-is

These are person-decades of specialist effort behind a fixed interface, re-tuned for each generation of hardware. Your job is not to reproduce them.

⚡ Performance Note — recognize the shape, call the library. The professional skill on display is pattern recognition: the moment you see a dense matrix multiply (or a solve, an eigenproblem, an FFT), the right move is not to optimize a loop but to recognize the shape and call the tuned routine. Hand-write the physics that is uniquely yours — your stencil, your boundary conditions, your model — and never hand-write the dgemm underneath it. This is the same division of labor Chapter 1 pointed at when it noted NumPy is fast because it calls this exact compiled code.

Phase 5 — The professional conclusion

So what was the point of building the blocked multiply if the answer is "call dgemm"? Two things. First, you now understand what dgemm is doing — multi-level blocking and vectorization to hit the arithmetic roofline — so it is no longer magic, and you can reason about its performance, its memory needs, and when even it is not the right tool (very small matrices, very sparse ones). Second, and more general: you learned to recognize a library-shaped problem. The optimization ladder in this chapter is for the kernels that are uniquely yours; for the standard ones, the fastest code you can write is the call that hands the work to someone who already climbed the ladder further than you ever will. Knowing which is which — that is the judgment §29.5 was really teaching.

Discussion Questions

  1. Cache blocking transformed the matrix multiply but did almost nothing for the heat stencil in Case Study 1. Explain the difference in one sentence using arithmetic intensity, and give the general rule for when blocking is worth trying.
  2. Your blocked loop reached ~15% of peak; dgemm reached ~80%. Is the remaining gap worth your time to close by hand? Under what rare circumstance (§29.4) would writing assembly intrinsics be the right call — and who, realistically, is the person who should write them?
  3. dgemm's reassociation of the sum can make its result differ from your naive loop's in the last bit. For a physics simulation, when does that matter and when is it noise? How would you decide the tolerance in a regression test (Chapter 20)?
  4. The chapter's thesis is "help the compiler, don't outsmart it — except in the critical 3%." Where does calling BLAS fit that thesis: is it helping the compiler, replacing it, or something else?

Your Turn: Extensions

  • Option A (build + measure). Implement all three matmul_* routines, verify them against matmul, then time each on a $1000\times1000$ matrix with the Chapter 28 timers. Compute the fraction of your CPU's theoretical peak each reaches (peak $\approx$ cores $\times$ clock $\times$ flops/cycle). Do you see the naive-few-percent → ordered → blocked climb?
  • Option B (tune). Sweep the block size nb over $\{16, 32, 64, 128, 256\}$ for matmul_blocked on a large matrix and plot time vs nb. Where is the sweet spot, and how does it relate to your L1/L2 cache size (§29.2's estimate)? Note how the best nb is found, not derived.
  • Option C (call the library). Link OpenBLAS and call dgemm on the same matrix; measure the speedup over your best blocked version and the fraction of peak it reaches. Then write the one-sentence commit message you would use when you delete your blocked loop in favor of the dgemm call.

Key Takeaways

  • Matrix multiply is blocking's home turf — high arithmetic intensity, $O(n)$ reuse — so tiling earns its complexity here, exactly where it did nothing for the memory-bound stencil.
  • The optimization ladder is real: naive (a few % of peak) → correct loop order → cache blocking (tens of %) → and still the tuned library beats you.
  • A tuned dgemm reaches ~80% of peak by blocking for every cache level, hand-vectorized assembly, and per-architecture tuning — an order of magnitude over your best loop, and person-decades you should not reproduce.
  • Recognize the shape, call the library. For dense linear algebra, the fastest code you can write is the BLAS/LAPACK call (Chapter 21). Hand-write only the physics that is uniquely yours.
  • Verify every step against matmul: bit-exact for integer test data; within tolerance for general reals, because reordering a floating-point reduction can move the last bit.