LAPACK
/ ELL-a-pack /
Almost every serious piece of scientific software, whether you wrote it in MATLAB, Python, R, or Julia, eventually leans on the same small set of battle-tested routines for the heavy lifting of linear algebra. LAPACK — short for Linear Algebra PACKage — is that set. It is a free, open, Fortran library, decades in the making, that solves the everyday workhorses of numerical computing: systems of linear equations, least-squares problems, eigenvalues and eigenvectors, and singular value decompositions. When you call solve, eig, or svd in a high-level language, a LAPACK routine is usually doing the work.
What makes LAPACK matter, beyond being free, is HOW it is built. It expresses its algorithms in terms of BLAS calls — and crucially, in terms of the Level-3 BLAS, the matrix-matrix operations that move large blocks of numbers at once. This 'blocked' design means LAPACK does its arithmetic on chunks that fit in cache, so it gets most of the speed of a tuned BLAS for free. Its routines follow a tidy naming scheme: a name like dgesv reads as d (double precision), ge (general matrix), sv (solve a linear system); dpotrf is d, po (symmetric positive-definite), trf (triangular factorization, i.e. Cholesky). Once you learn the code, the library's catalogue becomes readable. LAPACK is the successor to the older LINPACK and EISPACK, redesigned in the late 1980s specifically so that the inner loops would be cache-friendly on modern machines.
The deeper value is correctness and care. LAPACK's routines were written by leading numerical analysts and embody backward-error analysis, careful pivoting, and edge-case handling (overflow, underflow, scaling) that a hand-rolled version almost never gets right. Many routines return condition-number estimates or error bounds alongside the answer, so you can tell how much to trust it — which is exactly the trustworthiness discipline this field cares about. The honest caveats: LAPACK targets DENSE matrices on a single shared-memory machine (sparse and massively-parallel problems need other libraries like SuiteSparse, PETSc, or ScaLAPACK), and it is, after all, written in Fortran's column-major convention, which row-major callers (C, Python) must respect or silently transpose.
The routine dgesv solves A x = b for a general double-precision matrix. Internally it calls dgetrf (LU factorization with partial pivoting, expressed as blocked Level-3 BLAS) then dgetrs (forward and back substitution). MATLAB's backslash, SciPy's linalg.solve, and R's solve all ultimately route here — you almost never write the LU yourself.
One LAPACK call (dgesv) chains a blocked LU factorization with triangular solves, all on top of BLAS.
LAPACK is not magic against ill-conditioning: it gives a backward-stable answer, the exact solution to a nearby problem, but if your matrix is ill-conditioned the forward error can still be large. That is the problem's fault, not LAPACK's — and several routines will tell you so via a returned condition estimate.