the BLAS levels
/ BLAS rhymes with 'glass' /
Numerical linear algebra is built from a small set of recurring building blocks — adding vectors, multiplying a matrix by a vector, multiplying two matrices. Rather than have every program reinvent (and badly optimise) these, the community agreed on a standard catalogue of them called the BLAS, the Basic Linear Algebra Subprograms, and hardware vendors ship hand-tuned versions for each chip. The BLAS are organised into three levels, and the level number turns out to predict exactly how fast each operation can run on real hardware.
Level 1 is vector-vector operations: things like y = alpha*x + y (a scaled add, called axpy) or a dot product. For vectors of length n they do O(n) flops and touch O(n) data, so their arithmetic intensity is about 1 flop per element loaded — hopelessly memory-bound. Level 2 is matrix-vector operations, like y = A*x. They do O(n^2) flops and move O(n^2) data, again about 1 flop per element — still memory-bound. Level 3 is matrix-matrix operations, the crown being the general matrix multiply C = alpha*A*B + beta*C (called gemm). Here is the magic: O(n^3) flops over only O(n^2) data, so the arithmetic intensity grows like n. Each value, once loaded into cache, is reused about n times. That is the operation that can be blocked, vectorized, and parallelised until it runs at near the machine's peak.
The practical lesson is profound and a little counterintuitive: to make a numerical algorithm fast, you want to express as much of its work as possible as level-3 BLAS, even if that means doing the same flops in a different order. This is exactly why LAPACK (the higher-level library of solvers and factorizations) is rebuilt around block algorithms that call gemm in their inner loop — a blocked LU or Cholesky spends most of its flops in level-3 calls and so inherits their near-peak speed, while a naive textbook version drowns in memory-bound level-1 and level-2 work. 'Cast it as matrix-matrix multiply' is one of the most powerful heuristics in high-performance numerics.
Three ways to add a rank-one update to a matrix illustrate the levels. Looping element by element is level-1-like and crawls at memory speed. The level-2 call ger (A = A + x*y^T) is better but still ~1 flop/byte. Accumulating many such updates as one gemm (A = A + X*Y^T with X, Y tall-skinny) is level-3 and runs near peak — identical flops, but reuse per loaded value goes from O(1) to O(n).
Level 1, 2, 3 = vector, matrix-vector, matrix-matrix; only level 3 has the reuse to hit peak.
BLAS is an interface, not one implementation — the same gemm call may run 50x faster from a vendor-tuned library (MKL, OpenBLAS) than from a reference version, and the speed comes from blocking and vectorization, not from doing fewer flops. Always link a tuned BLAS; never hand-roll your own matrix multiply.