Case Study 1: Following the Thread

DataField.Dev

Case Study 1: Following the Thread

"Every high-level convenience is a low-level library wearing a nice coat."

Executive Summary

You have written numpy.linalg.solve(A, b) a hundred times without asking what happens underneath. In this case study we pull that thread all the way down — from a one-line Python call, through LAPACK, through BLAS, to the hardware-tuned kernel that does the actual arithmetic — and back up into the Fortran you will write in Chapter 21. The goal is not to call anything new (we cannot call LAPACK until Chapter 21); it is to read the stack: to learn to look at any numerical convenience and see the Fortran floor holding it up. By the end you will be able to name every layer between a scientist's Python and a processor's fused multiply-add, and you will have verified, by hand and in compilable Fortran, the operation that sits at the very bottom.

Skills applied: naming LAPACK and BLAS and their layering (§16.1); reading LAPACK routine names (§16.1); placing the tuned-BLAS decision (§16.1, and forward to §29); connecting the Python interop of Chapter 15 to the library floor of §16.1; a sanity check by "multiplying back" (a Level-2 BLAS idea).

Background

A data scientist on your team fits models in Python. Her code is unremarkable:

import numpy as np
x = np.linalg.solve(A, b)     # solve the linear system A x = b

She believes, reasonably, that this is "Python." It is not — or rather, the np.linalg.solve part is a thin Python wrapper over a tower of compiled libraries, most of it Fortran, some of it older than she is. When her solve is slow, or gives a surprising result on an ill-conditioned matrix, or when she asks why installing a different math library made her code twice as fast overnight, the answers all live below the Python. Our job is to descend the tower one floor at a time.

Phase 1 — The Top Floor: the Python Convenience

Start where she starts. np.linalg.solve accepts NumPy arrays and returns a solution vector; it validates shapes, checks for squareness, and — crucially — hands the real work to a compiled routine. NumPy does not implement Gaussian elimination in Python; that would be hopelessly slow, exactly the "pure-Python loop falls off a cliff" effect you measured in Chapter 15. Instead it marshals the arrays into the memory layout the compiled library expects and calls down.

Note one detail that will matter two floors down: NumPy arrays are C-contiguous (row-major) by default, while the library it is about to call is Fortran, which is column-major (Chapter 5). Somewhere in the descent, someone accounts for that transpose — the same order='F' issue you met wrapping your own kernel in Chapter 15. The convenience layer's real job is bookkeeping: shapes, types, contiguity, error checks. The arithmetic is elsewhere.

Phase 2 — The Solver Floor: LAPACK

The "elsewhere" is LAPACK. For a general dense system in double precision, np.linalg.solve ultimately drives LAPACK's gesv family — the routine you can now read on sight:

   D  GE  SV
   │  │   └─ solve  (LU-factor the matrix, then solve for x)
   │  └───── general matrix
   └──────── double precision, real

dgesv does two things: it factors $A$ into a product of triangular matrices (an LU factorization), then uses that factorization to solve for $\mathbf{x}$. This is the routine — precisely this name — that you will call yourself in Chapter 21. Seeing it here, holding up a NumPy call your colleague runs daily, is the whole point of naming LAPACK in this chapter rather than waiting: it is not exotic; it is already load-bearing in code you use.

LAPACK also explains one of her mysteries. On an ill-conditioned matrix, dgesv can return a solution with large error, and LAPACK exposes this through an info return code and companion routines that estimate the conditioning — machinery the Python layer often hides. The instinct "the library will just handle it" is exactly the instinct Chapter 20 will teach you to distrust. The floor is solid, but you have to read its warnings.

Phase 3 — The Engine Room: BLAS

LAPACK is not, however, where the arithmetic actually runs. LAPACK's algorithms are deliberately written to express their heavy work as calls to BLAS — and specifically to Level-3 BLAS, the matrix-matrix kernels, because that is where hardware runs fastest (§16.1). The LU factorization inside dgesv is blocked: it peels off a panel of columns, factors it, and then updates the rest of the matrix with a big matrix-matrix multiply (dgemm) and triangular solve (dtrsm) — Level-3 calls that do $O(n^3)$ arithmetic on $O(n^2)$ data.

This is the layer where the flops live. If you counted the floating-point operations your colleague's solve performs, essentially all of them happen inside a handful of BLAS calls, not in LAPACK's own code and certainly not in NumPy's Python. LAPACK is the orchestration; BLAS is the engine.

Phase 4 — The Foundation: the Tuned Kernel

One floor lower still is the answer to her overnight-speedup mystery. The BLAS interface is fixed, but the implementation is swappable (§16.1's threshold concept). Her original NumPy may have shipped with the reference BLAS or a modest kernel; installing a build linked against OpenBLAS or Intel MKL replaced the dgemm in the engine room with one blocked for her CPU's cache and issuing vector instructions by hand — same interface, same answer, several times the speed. She changed nothing in her Python and nothing in LAPACK; she swapped the foundation, and the whole tower got faster.

This is worth stating plainly to her, because it reframes performance work: you rarely make numerical code fast by rewriting your own loops; you make it fast by ensuring the tuned library is the one doing the work. That lesson returns, from the other direction, in Chapter 29, where your hand-rolled matrix multiply loses to exactly this tuned dgemm.

Phase 5 — Back Up to Fortran: the Sanity Check

Descend once more and you reach the arithmetic itself — a multiply and an add, repeated. You cannot call dgesv yet, but you can do the operation at the very bottom of the stack, and use it to check a solution the way a careful numericist always does: multiply back. If $\mathbf{x}$ really solves $A\mathbf{x}=\mathbf{b}$, then computing $A\mathbf{x}$ must reproduce $\mathbf{b}$. That check is a matrix-vector product — a Level-2 BLAS gemv — and the intrinsic matmul does it:

! cs01-verify-solution.f90 — check a hand-solved 2x2 system by multiplying back.
!   gfortran -std=f2018 -Wall cs01-verify-solution.f90 -o verify && ./verify
program verify_solution
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: a(2,2), x(2)

  ! System:  2x + 1y = 3
  !          1x + 3y = 5     ->  hand solution x = 0.8, y = 1.4
  a = reshape([2.0_dp, 1.0_dp, 1.0_dp, 3.0_dp], [2,2])   ! column-major
  x = [0.8_dp, 1.4_dp]

  print '(a)', 'A x  (must reproduce b = [3.00, 5.00]):'
  print '(2f6.2)', matmul(a, x)
end program verify_solution

! Hand computation:
!   (A x)(1) = 2*0.8 + 1*1.4 = 1.6 + 1.4 = 3.0
!   (A x)(2) = 1*0.8 + 3*1.4 = 0.8 + 4.2 = 5.0
!
! Expected output:
! A x  (must reproduce b = [3.00, 5.00]):
!   3.00  5.00

The output reproduces $\mathbf{b}$, so the solution checks out. In Chapter 21 you will let dgesv find $\mathbf{x}$ instead of solving by hand — but you will verify it exactly like this, by multiplying back, all the way down at the bottom of the tower your colleague never has to think about.

Discussion Questions

Your colleague says "I write Python, not Fortran." List, from top to bottom, every layer between her np.linalg.solve call and the processor's arithmetic, and say which layers are Fortran.
She installs a new math library and her code speeds up 2× with no code change. Using the layers above, explain exactly what changed and what did not.
Why does verifying a solution by computing $A\mathbf{x}$ and comparing to $\mathbf{b}$ catch a broad class of errors — and what kind of error (hint: conditioning) might it fail to reveal? (Foreshadows Chapter 20.)

Your Turn: Extensions

Option A. Pick a numerical convenience in a language you use — R's solve, MATLAB's backslash \, SciPy's linalg.eig, Julia's \ — and trace its documented backend. How many of them bottom out in LAPACK/BLAS? Write the tower for one of them.
Option B. Extend the sanity-check program to a $3 \times 3$ system: choose $A$ and $\mathbf{x}$, compute $\mathbf{b} = A\mathbf{x}$ by hand, then verify with matmul. Confirm your hand arithmetic matches the program's output (predict before you run).
Option C. Read the first paragraph of the reference documentation for dgesv (search "LAPACK dgesv"). List its arguments and match each to a concept you already know: which is $A$, which is $\mathbf{b}$, which returns the answer, and what does info report? This is your on-ramp to Chapter 21.

Key Takeaways

A high-level numerical call is a tower: Python convenience → LAPACK (the solver, e.g. dgesv) → BLAS (the engine, Level-3 kernels) → a hardware-tuned implementation. Most of the tower, and essentially all the flops, are Fortran and the libraries built on it.
You can read a LAPACK name on sight now, which means the "magic" under NumPy is legible, not mysterious.
Performance in this world usually comes from which library does the work, not from rewriting your own loops — swap in a tuned BLAS and the whole tower accelerates.
Always sanity-check a solution by multiplying back; it is a cheap Level-2 operation and it catches real mistakes. Reading a stack — seeing the Fortran floor under the Python — is the skill this chapter buys you.