Chapter 21 — Key Takeaways (Linear Algebra and LAPACK)

A one-page reference to the chapter where you stop writing linear algebra and start calling it. Keep it beside you the first dozen times you invoke LAPACK.

The vocabulary

Term Meaning
LAPACK Linear Algebra PACKage — the Fortran library of dense/banded solvers, eigensolvers, and factorizations
BLAS Basic Linear Algebra Subprograms — the vector/matrix kernels LAPACK is built on
BLAS levels 1 = vector–vector ($O(n)$, memory-bound); 2 = matrix–vector ($O(n^2)$, memory-bound); 3 = matrix–matrix ($O(n^3)$ on $O(n^2)$ data, compute-bound, near-peak)
leading dimension (lda) the declared first dimension of a 2-D array — the column stride in memory, not the size you use
info LAPACK status: 0 success, < 0 the (-info)-th argument was illegal (your bug), > 0 a numerical failure
ipiv integer pivot array recording partial-pivoting row interchanges
sparse format store only nonzeros + indices (COO / CSR / CSC): $O(n)$ memory for an $O(n)$-nonzero matrix

The LAPACK naming scheme (read any routine)

  d   ge   sv        precision + matrix type + computation
  |   |    |
  |   |    +--- sv=solve  trf=LU/Cholesky  trs=solve-w/-factors  ev=eigen  svd=SVD  gels=least-sq
  |   +-------- ge=general  sy=symmetric  po=SPD  gt=tridiagonal  tr=triangular  he=Hermitian
  +------------ s=real32  d=real64  c=complex32  z=complex64
Name Decodes to Does
dgesv double · general · solve solve $A\mathbf{x}=\mathbf{b}$
dsyev double · symmetric · eigen eigenvalues/vectors of a symmetric matrix
dgesvd double · general · SVD singular value decomposition
dgtsv double · tridiagonal · solve solve a tridiagonal system in $O(n)$
dgels double · general · least-squares least-squares (QR); better than normal equations
dgetrf/dgetrs LU factorize / solve-with-factors reuse one factorization for many right-hand sides

The two calls to memorize

! Solve A x = b.  a -> LU factors, b -> solution x (BOTH overwritten). Check info!
call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)   ! nrhs=1, lda=ldb=n for a full n x n

! Eigenvalues of a symmetric matrix, with the WORKSPACE QUERY:
allocate(work(1))
call dsyev('N', 'U', n, a, n, w, work, -1, info)  ! query: optimal size -> work(1)
lwork = int(work(1)); deallocate(work); allocate(work(lwork))
call dsyev('N', 'U', n, a, n, w, work, lwork, info)   ! w = eigenvalues, ascending

Reading a LAPACK reference page (the drill)

When you meet a routine you have not used, do not memorize it — read its page in five passes:

  1. Decode the name (precision · matrix type · computation) to confirm it does what you want.
  2. Find the argument list and mark each input / output / input-output — especially which arrays get overwritten.
  3. Spot the leading dimensions (lda, ldb, ldu, …) and pass the declared first dimensions.
  4. Check for a workspace query (work, lwork): if present, call once with lwork = -1, then allocate.
  5. Read the info codes for this routine's > 0 meaning (singular? not converged?), and always test it.

Which routine / when

You want to… Reach for
solve one (or a few) dense $A\mathbf{x}=\mathbf{b}$ dgesv
solve a tridiagonal system dgtsv ($O(n)$, not dgesv's $O(n^3)$)
eigenvalues of a symmetric matrix dsyev
the SVD of any matrix dgesvd
least squares (overdetermined) dgels (QR) or SVD — not normal equations for high degree
many right-hand sides, same $A$ dgetrf once, then dgetrs per RHS
a big sparse system a sparse solver (SuperLU/MUMPS/PETSc, or CG/GMRES) — never dgesv
multiply matrices, large the BLAS dgemm (or matmul); never a hand loop for production

Pitfalls

  • dgesv overwrites BOTH a and b. a becomes the LU factors, b becomes the solution. Copy them first if you need the originals (e.g. to compute a residual).
  • Always check info. A failed, unchecked LAPACK call returns a buffer of garbage that looks like an answer.
  • Column-major, twice. reshape([...], [n,n]) fills column-major (the transpose of what you typed) — use order=[2,1]. Printing a whole matrix streams it column-major — print row by row.
  • lda is the declared first dimension, not the used size. For a submatrix of a big array, lda stays the big array's first dimension.
  • undefined reference to 'dgesv_' is a LINK error, not a compile error. Add -llapack -lblas.
  • No interface = silent argument bugs. External LAPACK calls are not argument-checked by -Wall; a swapped argument compiles and then corrupts memory. Check info, or use a checked wrapper (stdlib solve).
  • Normal equations square the condition number ($\kappa(A^{\mathsf T}A)\approx\kappa(A)^2$). Prefer dgels/SVD.

Compile flags introduced

$ gfortran -std=f2018 -Wall solve.f90 -o solve -llapack -lblas    # link LAPACK, then BLAS
$ gfortran ... -lopenblas                                          # tuned BLAS+LAPACK, no code change

Numbers worth carrying

  • Dense LU (dgesv) costs about $\tfrac{2}{3}n^3$ flops; a matmul, about $2n^3$.
  • Only Level-3 BLAS ($O(n^3)$ work on $O(n^2)$ data) reaches near-peak speed; Levels 1–2 are memory-bound.
  • A condition number $\kappa(A)\approx10^{k}$ costs about $k$ of your ~16 double-precision digits.
  • A tuned BLAS beats a hand-written matrix multiply by ~10× and up — measured in Ch. 29.

Project piece added this chapter

The solver gains an optional implicit stepper: one backward-Euler step is the tridiagonal system $A\mathbf{u}^{n+1}=\mathbf{u}^{n}$ (diagonal $1+2r$, off-diagonals $-r$), solved with dgesv:

call dgesv(n, 1, a, n, ipiv, rhs, n, info)   ! rhs (u^n) -> u^{n+1}

This is the optional path; the default explicit stepping (and the CFL limit implicit stepping escapes) is built in Chapter 24. For a real 1-D run, swap the dense dgesv for the tridiagonal dgtsv. Saved as heat-solver/heat_implicit.f90.