BLAS and LAPACK
/ BLASS / LAY-pack /
If you write your own matrix multiply or your own Gaussian elimination as a naive triple loop, it will be correct but often ten or a hundred times slower than it should be, and possibly less accurate. BLAS and LAPACK are the standard, battle-tested libraries that solve this once and for all: BLAS provides the low-level building blocks and LAPACK provides the high-level solvers built on top of them. They are the engine inside MATLAB, NumPy, R, Julia and essentially every scientific computing environment.
BLAS, the Basic Linear Algebra Subprograms, is a standardized set of operations organized in three levels: Level 1 is vector-vector work (like a dot product, O(n) flops), Level 2 is matrix-vector work (like A x, O(n^2) flops), and Level 3 is matrix-matrix work (like A B, O(n^3) flops). The crucial insight is that Level 3 operations do O(n^3) arithmetic on only O(n^2) data, so they reuse each number many times and can be made cache-efficient — the heavily tuned implementations (OpenBLAS, Intel MKL, BLIS) run near the hardware's peak speed. LAPACK, the Linear Algebra PACKage, builds the algorithms you actually call — LU, Cholesky, QR, eigenvalues, SVD, least squares — and crucially it reorganizes them as BLOCKED algorithms that spend most of their time in fast Level-3 BLAS rather than slow element-by-element loops.
The practical lesson is the mantra 'don't roll your own': for standard dense linear algebra, call LAPACK (which calls a tuned BLAS) rather than writing the loops yourself, because the library is faster, more accurate, more thoroughly tested, and already handles pivoting and edge cases. The blocking is the key idea worth understanding even if you never touch the code: it exists because on modern hardware the bottleneck is data movement through the memory hierarchy, not flops, so algorithms are restructured to maximize the arithmetic done per byte fetched from memory. The one honest caveat: for sparse or highly structured matrices the dense BLAS/LAPACK is the wrong tool — you reach for sparse libraries instead — and results can differ in the last bits across BLAS implementations because floating-point addition is not associative.
NumPy's numpy.linalg.solve(A, b) calls LAPACK's dgesv, which factors A as P A = L U with partial pivoting via blocked Level-3 BLAS and then back-substitutes — far faster and safer than a hand-written loop.
A one-line high-level call dispatches to decades of tuned, blocked, pivoting-aware library code.
Blocking exists because modern performance is limited by memory movement, not flop count; for sparse or structured matrices, though, dense BLAS/LAPACK is the wrong tool. Results can vary in the last bits across implementations since floating-point addition is not associative.