Chapter 21: Linear Algebra: Solving Ax=b, Matrix Operations, and Using LAPACK

"The purpose of computing is insight, not numbers." — Richard W. Hamming, Numerical Methods for Scientists and Engineers (1962)

Overview

Back in Chapter 1 we made a claim that probably sounded like advocacy: that when your Python calls numpy.linalg.solve, or your MATLAB session inverts a matrix, the actual arithmetic very often happens inside compiled Fortran that has been tuned for four decades. In Chapter 16 we named that Fortran: LAPACK, the Linear Algebra PACKage, and the BLAS beneath it. This is the chapter where the claim stops being trivia and becomes a skill. You are going to call LAPACK directly, from your own Fortran program, to solve a system of equations, to find the eigenvalues of a matrix, and to decompose a matrix into its singular values — the three operations that sit underneath an enormous fraction of scientific computing. This is the LAPACK we met in Chapter 16, now doing real work.

The lesson underneath the mechanics is the more important one, and it runs against a beginner's instinct. Faced with "solve $A\mathbf{x} = \mathbf{b}$," the natural impulse of someone who has just learned about arrays and loops is to write the solver — to code up Gaussian elimination, because you know how it works and it feels like the honest thing to do. This chapter will teach you to resist that impulse for any serious computation. Not because you cannot write Gaussian elimination — you can, and we will — but because a correct, naive implementation will be numerically fragile and, for a matrix of any size, some tens of times slower than the routine a specialist wrote, tested against decades of pathological cases, and tuned to the exact shape of your processor's cache. The mark of an expert here is knowing what not to build. Your job is to state the problem in the standard form, hand it to the library, and check the answer.

In this chapter, you will learn to:

• Lay matrices out the way Fortran and LAPACK both expect them — column-major — and use the matmul, dot_product, and transpose intrinsics for the small cases where they are the right tool.
• Write a matrix multiply from scratch, then measure it against the intrinsic and the BLAS in your head, and understand the three levels of the BLAS that explain the gap.
• Call dgesv to solve a dense linear system, with every argument — the leading dimension, the pivot array ipiv, and the status flag info — understood, not copied on faith.
• Read a LAPACK reference page and decode any routine name from its letters, then drive the workspace-query pattern used by dsyev for eigenvalues and dgesvd for the singular value decomposition.
• Link your program against LAPACK and BLAS (-llapack -lblas), and know when to reach for a tuned implementation (OpenBLAS, MKL) or, when your matrix is mostly zeros, a sparse solver instead.

Learning Paths

How to read this chapter by track. - 🔬 Scientist — this is one of the highest-value chapters in the book for you. Read §21.1 and §21.3 closely; §21.3 is the one you will use tomorrow. Skim §21.2, and treat §21.6 as a map for when your matrices get large. - 📖 Standard — read straight through. §21.3–21.4 are the payoff of the whole "call the library" arc. - 🔧 Legacy — the LAPACK you will inherit is called exactly this way; the F77-style external interface in §21.3 is why argument-order mistakes are silent, which is your daily hazard. - ⚡ HPC — §21.2 (BLAS levels), §21.5 (OpenBLAS/MKL), and the sparse pointer in §21.6 are your map to where the FLOPS actually are. The tuned-BLAS-beats-your-loop point returns, measured, in Chapter 29.

21.1 Matrices in Memory, and the Array Intrinsics You Already Have

A matrix, to Fortran, is just a rank-2 array. You met it in Chapter 5: real(dp) :: a(3, 3) is a 3×3 matrix, a(i, j) is the element in row i, column j, and the whole apparatus of sections and whole-array operations applies. There is nothing new to learn about declaring a matrix. What you must hold firmly in mind — because it governs every library call in this chapter — is how the matrix sits in memory.

Fortran stores arrays in column-major order: the first index varies fastest, so the elements of a(:,:) are laid down a(1,1), a(2,1), a(3,1), a(1,2), … — one whole column contiguously, then the next. (C, C++, and NumPy default to the opposite, row-major.) You learned this in §5.6 as a performance idea, the reason the inner loop belongs on the first index. Here it acquires a second life: LAPACK and the BLAS were written in Fortran and expect column-major matrices, so the layout you already respect for speed is also the layout the library demands. This is a happy coincidence you should not have to think about — until the day you share a matrix with C or NumPy, when it becomes the bug that eats an afternoon.

The immediate consequence shows up the moment you try to type a matrix into your source code. Suppose you want the matrix

$$ B = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}. $$

The tempting thing is reshape([1, 2, 3, 4], [2, 2]). But reshape pours the list into the array in storage order — column by column — so it fills b(1,1)=1, b(2,1)=2, b(1,2)=3, b(2,2)=4, giving you

$$ \begin{bmatrix} 1 & 3 \\ 2 & 4 \end{bmatrix}, $$

the transpose of what you wrote. To make the source read like the mathematics, pass order=[2, 1], which tells reshape to fill the second dimension (the columns) fastest — that is, row by row:

b = reshape([ 1.0_dp, 2.0_dp,   &
              3.0_dp, 4.0_dp ], [2, 2], order=[2, 1])   ! really [[1,2],[3,4]]

⚠️ Common Pitfall: reshape([1,2,3,4], [2,2]) does not give you the matrix whose rows are 1 2 and 3 4. Column-major fill makes it the transpose. Either write the elements in column order, or use order=[2,1] and write them row by row. We use order=[2,1] throughout this book so the source matches the maths — but the machine still stores it column-major, and that is what LAPACK reads.

For the small dense operations, the intrinsics from Chapter 5 are exactly what you want, and you should prefer them over any hand-written loop. matmul(a, b) is the matrix product $C_{ij} = \sum_k a_{ik} b_{kj}$; `dot_product(x, y)` is the scalar $\sum_i x_i y_i$; transpose(a) swaps rows and columns. Here they are on hand-checkable data:

program matrix_storage
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: b(2,2), c(2,2), prod(2,2)
  real(dp) :: u(3) = [1.0_dp, 2.0_dp, 3.0_dp]
  real(dp) :: v(3) = [4.0_dp, 5.0_dp, 6.0_dp]
  integer  :: i

  b = reshape([ 1.0_dp, 2.0_dp,  3.0_dp, 4.0_dp ], [2, 2], order=[2, 1])   ! [[1,2],[3,4]]
  c = reshape([ 5.0_dp, 6.0_dp,  7.0_dp, 8.0_dp ], [2, 2], order=[2, 1])   ! [[5,6],[7,8]]
  prod = matmul(b, c)

  print '(a)', 'matmul(b, c):'
  do i = 1, 2
    print '(2f7.1)', prod(i, :)                 ! print row i, so it reads normally
  end do
  print '(a, f7.1)', 'dot_product(u, v) = ', dot_product(u, v)
end program matrix_storage

$ gfortran -std=f2018 -Wall -O2 matrix_storage.f90 -o matstore && ./matstore
matmul(b, c):
   19.0   22.0
   43.0   50.0
dot_product(u, v) =    32.0

Work the product by hand to trust it: the top-left entry is (row 1 of $B$)·(column 1 of $C$) $= 1\cdot5 + 2\cdot7 = 19$; the top-right is $1\cdot6 + 2\cdot8 = 22$; the bottom row is $3\cdot5+4\cdot7 = 43$ and $3\cdot6+4\cdot8 = 50$. And $\mathbf{u}\cdot\mathbf{v} = 1\cdot4 + 2\cdot5 + 3\cdot6 = 4 + 10 + 18 = 32$. Notice the printing detail: because the array is stored column-major, printing prod whole would stream it in memory order and show you the transpose; printing prod(i, :) one row at a time is how you get output that reads like a matrix. The same column-major fact, biting twice in one short program — once on input via reshape, once on output via print — is a fair warning of how insistently it will matter.

🐍 Python Comparison: In NumPy you would write B @ C for the matrix product and np.dot(u, v) for the inner product, and NumPy stores B row-major by default. When you eventually hand a NumPy array to a Fortran routine (via f2py, Chapter 15), you pass order='F' to get column-major, or you accept a silent transpose. The intrinsics themselves are a wash for correctness; the memory order is where the two languages actually disagree.

🔄 Check Your Understanding. 1. You write a = reshape([1.0_dp, 0.0_dp, 0.0_dp, 1.0_dp], [2,2]) intending the identity. Do you get it? 2. What is the shape of matmul(m, x) when m is (3,3) and x is a length-3 vector?

Answers

1. Yes, by luck — the identity is symmetric, so its column-major and row-major fills are identical. Try it with a non-symmetric matrix and the trap reappears. 2. A length-3 rank-1 array; matmul accepts a matrix times a vector and returns a vector.

21.2 Writing Your Own Matrix Multiply — and Why You Should Not Ship It

You should write a matrix multiply once in your life, to know that you can, and to feel in your hands why the library exists. Here it is, in the loop order that respects column-major layout (inner loop over the first index i, so each innermost pass walks down a column with the memory grain):

pure subroutine my_matmul(a, b, c)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  real(dp), intent(in)  :: a(:,:), b(:,:)
  real(dp), intent(out) :: c(:,:)
  integer :: i, j, k
  c = 0.0_dp
  do j = 1, size(b, 2)              ! each column of the result
    do k = 1, size(a, 2)            ! accumulate rank-1 updates
      do i = 1, size(a, 1)          ! inner loop over the FIRST index (column-major)
        c(i, j) = c(i, j) + a(i, k) * b(k, j)
      end do
    end do
  end do
end subroutine my_matmul

This is correct. It even has the cache-friendly loop order that a naive i, j, k version would get wrong. Compile it, test it against matmul, and it will agree to the last bit for small matrices. And you should still never ship it for large ones, because a tuned BLAS will beat it by a large factor — often more than ten times, sometimes far more — and the reason is instructive. Your triple loop touches each element of a and b many times but does only a couple of arithmetic operations each time it fetches them from memory; for large matrices it spends most of its life waiting for memory, not computing. A tuned matrix multiply restructures the work into small blocks that fit in cache, so that once a block of a is loaded it is reused for many multiplications before being evicted — turning a memory-bound computation into a compute-bound one that can approach the processor's peak arithmetic rate. It also uses the CPU's vector (SIMD) instructions deliberately, prefetches the next block while computing the current one, and was tuned to the exact cache sizes of real processors. That is a specialist's month of work, and it is sitting in a library you can link in one flag.

The library in question is the BLAS, and it is organized into three levels by exactly the ratio that governs the performance above — how much arithmetic you do per element of memory you touch.

Definition (BLAS levels). The BLAS (Basic Linear Algebra Subprograms) are grouped into three levels by the shape of their operands. Level 1 is vector–vector work — $O(n)$ operations on $O(n)$ data (e.g. daxpy, $\mathbf{y} \leftarrow a\mathbf{x} + \mathbf{y}$; ddot, an inner product). Level 2 is matrix–vector — $O(n^2)$ operations on $O(n^2)$ data (e.g. dgemv, $\mathbf{y} \leftarrow \alpha A\mathbf{x} + \beta\mathbf{y}$). **Level 3** is matrix–matrix — $O(n^3)$ operations on only $O(n^2)$ data (e.g. dgemm, $C \leftarrow \alpha AB + \beta C$). Level 3 is special: because it does $n$ operations for every element it loads, it can reuse data in cache and reach near-peak speed. Level 1 and 2 are memory-bound and cannot. This ratio — arithmetic per byte, called arithmetic intensity — is why "cast your computation as Level-3 BLAS calls" is the single most repeated piece of HPC advice.

So the honest replacement for my_matmul is not a cleverer loop; it is dgemm, the Level-3 BLAS matrix–matrix routine, or the matmul intrinsic (which a good compiler may implement by calling the BLAS, or may inline and vectorize itself). For small matrices — say, a 3×3 rotation, or the little systems in this chapter's examples — matmul is perfectly fine and clearer than a BLAS call. The crossover where hand-rolled loops become a real mistake arrives sooner than beginners expect; by the time a matrix is a few hundred on a side, the tuned routine is winning comfortably.

⚡ Performance Note: The gap between a naive triple loop and a tuned dgemm is not a constant factor you can ignore — it grows with the matrix, because the naive version becomes ever more memory-bound while dgemm stays compute-bound. We measure exactly this, and see a hand-written loop lose to the library, in Chapter 29. The lesson lands harder when it is your own loop on the losing side.

🚪 Threshold Concept. The expert instinct in numerical computing is not "I can implement this." It is "someone has implemented this better than I can afford to, and my job is to call it correctly and check the result." Writing your own linear-algebra kernels for production is, with rare exceptions, a way to be slower and more wrong at the same time. This inverts the pride of a new programmer, and internalizing it is a large part of becoming a computational scientist rather than a hobbyist.

21.3 Solving $A\mathbf{x} = \mathbf{b}$ with dgesv

Here is the heart of the chapter. You have a system of linear equations — $n$ equations in $n$ unknowns — written in matrix form as $A\mathbf{x} = \mathbf{b}$, with $A$ a known $n \times n$ matrix, $\mathbf{b}$ a known right-hand-side vector, and $\mathbf{x}$ the unknown you want. The professional way to find $\mathbf{x}$ in Fortran is one call to LAPACK's dgesv. Before we make the call, decode its name, because LAPACK's naming scheme is a language you can learn to read in five minutes and then never need a manual for again.

Definition (LAPACK routine naming). A LAPACK routine name is built from letters that encode, in order: the data type, the matrix type, and the computation. The first letter is the precision — s single real, d double real, c single complex, z double complex. The next two letters name the matrix structure — ge general, sy symmetric, po symmetric positive-definite, gt general tridiagonal, tr triangular, and so on. The final letters name the job — sv solve a linear system, trf compute a triangular (LU/Cholesky) factorization, trs solve using that factorization, ev eigenvalues and eigenvectors, svd singular value decomposition. So dgesv parses as double + general + solve: "solve a general linear system in double precision." Once you can read the name, you can guess the routine you need and confirm it in the reference — dsyev is double/symmetric/eigenvalues, dgesvd is double/general/SVD, sgetrf is single/general/LU-factorization.

The interface, which you will see in every LAPACK-using program you ever read, is:

call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)

Each argument earns its place. Read them slowly, because the whole skill is here:

• n — the order of the system: $A$ is n×n. (Input.)
• nrhs — the number of right-hand sides. You usually have one vector $\mathbf{b}$, so nrhs = 1; but LAPACK will solve for many right-hand sides at once if you pass $B$ as an n×nrhs matrix. (Input.)
• a — the coefficient matrix $A$, an lda×n array. On exit it is overwritten with the $L$ and $U$ factors of its LU decomposition. Your original $A$ is gone unless you saved a copy. (Input/output.)
• lda — the leading dimension of a (defined below).
• ipiv — an integer array of length n, the pivot indices. LAPACK does Gaussian elimination with partial pivoting — it reorders rows for numerical stability — and records the swaps here. On exit it tells you which rows were interchanged; you rarely read it, but you must pass it. (Output.)
• b — on entry the right-hand side(s), an ldb×nrhs array; on exit it is overwritten with the solution $\mathbf{x}$. This is the answer. (Input/output.)
• ldb — the leading dimension of b.
• info — the status flag (defined below). (Output.)

Two of those arguments — lda/ldb and info — are the ones newcomers get wrong, so they get their own definitions.

Definition (leading dimension, lda). The leading dimension of a two-dimensional array is the number of rows in the array as it is declared in memory — not necessarily the number of rows you are using. Because Fortran stores column-major, consecutive columns of a are lda elements apart in memory; the leading dimension is the stride from one column to the next. If you declare real(dp) :: a(n, n) and use all of it, then lda = n. But if you declare real(dp) :: a(100, 100) and solve only the top-left n×n block, you must pass lda = 100 — the declared first dimension — so LAPACK can step correctly from column to column through the larger array. Getting lda wrong does not usually crash; it silently reads the wrong elements. When in doubt, lda is the first number in the array's declaration.

💡 Intuition: Think of the matrix as living in a filing cabinet whose drawers (columns) are a fixed width apart. lda is that fixed width — the distance you jump to reach the next column — and it is a property of the cabinet (the declared array), not of how many folders you happen to be using in each drawer. That is why a submatrix of a big array keeps the big array's lda.

Definition (info). Every LAPACK computational routine returns an integer info as its last argument, reporting what happened. info = 0 means success. info < 0 means the (-info)-th argument had an illegal value — for info = -4, your fourth argument (lda) was invalid; this is almost always a bug in your call. info > 0 means a numerical problem specific to the routine: for dgesv, info = i means $U_{ii}$ came out exactly zero, so $A$ is exactly singular and no solution was computed. Always check info. A routine that fails and is not checked hands you a buffer of garbage that looks exactly like an answer.

Now the worked example. We solve

$$ \begin{aligned} x + y + z &= 6 \\ 2y + 5z &= -4 \\ 2x + 5y - z &= 27 \end{aligned} \qquad\Longleftrightarrow\qquad A = \begin{bmatrix} 1 & 1 & 1 \\ 0 & 2 & 5 \\ 2 & 5 & -1 \end{bmatrix},\quad \mathbf{b} = \begin{bmatrix} 6 \\ -4 \\ 27 \end{bmatrix}. $$

The solution, which you can verify by substitution, is $\mathbf{x} = (5, 3, -2)$: check the first equation, $5 + 3 - 2 = 6$; the second, $2\cdot3 + 5\cdot(-2) = 6 - 10 = -4$; the third, $2\cdot5 + 5\cdot3 - (-2) = 10 + 15 + 2 = 27$. Here is the program that lets LAPACK find it:

program solve_linear_system
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer, parameter :: n = 3
  real(dp) :: a(n,n), a_orig(n,n)      ! keep a copy: dgesv destroys a
  real(dp) :: b(n),   b_orig(n)        ! keep a copy: dgesv overwrites b with x
  integer  :: ipiv(n), info

  a = reshape([ 1.0_dp,  1.0_dp,  1.0_dp,   &
                0.0_dp,  2.0_dp,  5.0_dp,   &
                2.0_dp,  5.0_dp, -1.0_dp ], [n, n], order=[2, 1])
  b = [ 6.0_dp, -4.0_dp, 27.0_dp ]

  a_orig = a
  b_orig = b

  call dgesv(n, 1, a, n, ipiv, b, n, info)     ! nrhs=1, lda=ldb=n

  if (info /= 0) then
    print '(a, i0)', 'dgesv failed, info = ', info
    error stop 1
  end if

  print '(a)', 'Solution vector x:'
  print '(f10.4)', b                            ! b now holds x
  print '(a, es9.2)', 'Max residual: ', maxval(abs(matmul(a_orig, b) - b_orig))
end program solve_linear_system

$ gfortran -std=f2018 -Wall solve_linear_system.f90 -o solve -llapack -lblas && ./solve
Solution vector x:
    5.0000
    3.0000
   -2.0000
Max residual:  0.00E+00

The residual line is the discipline that separates a working scientist from a lucky one. After the solve, b holds $\mathbf{x}$, but a has been overwritten by its LU factors — which is exactly why we saved a_orig and b_orig before the call. The largest-magnitude entry of $A\mathbf{x} - \mathbf{b}$, computed from the originals, tells us how well the reported $\mathbf{x}$ actually satisfies the equations. Here every number in the elimination happens to be an exact binary fraction, so the residual is exactly zero; for a general system it will be a tiny number near the rounding floor, and if it ever comes back large, the solve — or your setup — is wrong, and you have caught it.

⚠️ Common Pitfall — the silent argument. LAPACK routines are plain external procedures with no Fortran interface visible to your compiler, so -Wall cannot check the arguments of a dgesv call. Swap two arguments — pass ldb where nrhs belongs, or a rank-1 b of the wrong length — and the program compiles without a murmur and then reads or writes the wrong memory at run time. This is the number-one LAPACK bug. Two defenses: check info every single time, and prefer an interface that does get checked — the fortran-lang stdlib ships a stdlib_linalg module whose solve(A, b) wrapper is type-checked and returns the solution as a function result. Use the raw dgesv when you must (legacy code, exotic options); reach for a checked wrapper when you can.

🐍 Python Comparison: numpy.linalg.solve(A, b) is this exact dgesv, wrapped. NumPy makes a copy so your A survives, allocates the pivot array for you, raises LinAlgError instead of returning info > 0, and hides the leading dimension entirely. Everything NumPy does for you, you are now doing by hand — which is why the Fortran is longer, and also why the Fortran has no per-call allocation or copy you did not ask for. When the solve is the inner loop of a simulation run a million times, that control is the point.

🔗 Connection — conditioning. dgesv returning info = 0 means it finished, not that the answer is accurate. If $A$ is ill-conditioned — nearly singular, so the columns almost line up — then small rounding errors in the data or the arithmetic are amplified into large errors in $\mathbf{x}$, and a small residual can still hide a wrong answer. This is the conditioning idea from Chapter 20: a condition number $\kappa(A) \approx 10^{k}$ costs you about $k$ of your ~16 significant decimal digits in double precision. LAPACK gives you dgecon to estimate $\kappa(A)$ from the factorization; when a solution matters, estimate the condition number, do not just trust info = 0.

🔄 Check Your Understanding. 1. After call dgesv(...), where is the solution, and what happened to your original matrix? 2. dgesv returns info = -6. What does that tell you, and whose fault is it likely to be? 3. Why must you pass ipiv even though you never read it?

Answers

1. The solution overwrites b; a now holds the $L$ and $U$ factors, so the original $A$ is gone unless you copied it. 2. The sixth argument (b) had an illegal value — an info < 0 is a bad argument in your call (probably a wrong shape or leading dimension), not a numerical failure. 3. dgesv uses partial pivoting for stability and needs somewhere to record the row interchanges; ipiv is required workspace-plus-output even when you do not inspect it.

21.4 Eigenvalues, the SVD, and Reading a Reference Page

Solving $A\mathbf{x} = \mathbf{b}$ is one of three workhorses. The other two are the eigenvalue problem — find the $\lambda$ and $\mathbf{v}$ with $A\mathbf{v} = \lambda\mathbf{v}$ — and the singular value decomposition, $A = U\Sigma V^{\mathsf{T}}$, which factors any matrix into a rotation, a scaling by the singular values, and another rotation. Eigenvalues tell you a system's natural frequencies, a matrix's stability, the principal axes of a dataset; the SVD underlies least-squares fitting, data compression, principal component analysis, and the honest way to handle a rank-deficient matrix. LAPACK does both, and calling them teaches the one new mechanic these routines add: the workspace query.

The routine for a real symmetric matrix's eigenvalues is dsyev (double / symmetric / eigenvalues):

call dsyev(jobz, uplo, n, a, lda, w, work, lwork, info)

• jobz — 'N' to compute eigenvalues only, 'V' to compute eigenvectors as well.
• uplo — 'U' or 'L': which triangle of the symmetric a you filled (LAPACK reads only one).
• n, a, lda — the order, the matrix, and its leading dimension, as before. On exit, if jobz='V', a holds the orthonormal eigenvectors in its columns; otherwise a is destroyed.
• w — a real array of length n; on exit, the eigenvalues in ascending order.
• work, lwork — a real work array and its length. This is the new part.
• info — status, as always: 0 success, < 0 bad argument, > 0 the algorithm failed to converge.

Many LAPACK routines need scratch space whose optimal size depends on the matrix and the machine, so rather than guess, you ask. Call the routine once with lwork = -1 — the workspace query — and it does no real work but writes the optimal lwork into work(1). You read that, allocate a work array of that size, and call again for real:

program eigenvalues
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter   :: n = 3
  real(dp)              :: a(n,n), w(n)
  real(dp), allocatable :: work(:)
  integer               :: lwork, info

  a = reshape([ 2.0_dp, 1.0_dp, 0.0_dp,   &
                1.0_dp, 2.0_dp, 1.0_dp,   &
                0.0_dp, 1.0_dp, 2.0_dp ], [n, n], order=[2, 1])

  allocate(work(1))
  call dsyev('N', 'U', n, a, n, w, work, -1, info)   ! query: how much workspace?
  lwork = int(work(1))
  deallocate(work); allocate(work(lwork))            ! now provide it

  call dsyev('N', 'U', n, a, n, w, work, lwork, info)
  if (info /= 0) then
    print '(a, i0)', 'dsyev failed, info = ', info
    error stop 1
  end if

  print '(a)', 'eigenvalues (ascending):'
  print '(f10.4)', w
end program eigenvalues

$ gfortran -std=f2018 -Wall eigenvalues.f90 -o eig -llapack -lblas && ./eig
eigenvalues (ascending):
    0.5858
    2.0000
    3.4142

The matrix is the symmetric tridiagonal $\begin{bmatrix} 2 & 1 & 0\\ 1 & 2 & 1\\ 0 & 1 & 2\end{bmatrix}$, whose eigenvalues have the closed form $2 + \sqrt{2}\cos(k\pi/4)$ for $k = 1, 2, 3$ — that is, $2 - \sqrt{2} \approx 0.5858$, exactly $2$, and $2 + \sqrt{2} \approx 3.4142$. LAPACK returns them ascending in w, matching to four decimals. (This tridiagonal is not a random choice: it is precisely the shape of the discrete Laplacian you have been building since Chapter 5, and its eigenvalues are the mathematical reason the explicit heat scheme in Chapter 24 has a stability limit at all.)

The SVD routine, dgesvd, is a longer call in the same spirit, and it is the perfect example for the real skill of this section: reading a reference page rather than memorizing an interface. Its signature is

call dgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)

and you do not need it in your head. You need to know how to read it: m and n are the dimensions of the possibly-rectangular a; s receives the singular values (descending); u and vt receive the left and right singular vectors $U$ and $V^{\mathsf{T}}$; jobu and jobvt (with values 'A', 'S', 'O', 'N') choose how much of $U$ and $V^{\mathsf{T}}$ to compute; and work/lwork drive the same query pattern as dsyev. When you meet a routine you have not used, you open its reference page — from the man pages, the Netlib documentation, or the Intel MKL reference — and read the argument list exactly as we just read dgesv's: which are input, which are overwritten on output, which are the leading dimensions, and which is info. That is a transferable skill; memorizing dgesvd's thirteen arguments is not.

🧩 Try It Yourself: Before reading on, decode these three routine names from their letters, then check yourself against §21.3's scheme: dpotrf, zgeev, sgels. What precision, what matrix type, and what computation is each? (One is a Cholesky factorization; one finds the eigenvalues of a general complex matrix; one solves a least-squares problem.) Being able to read a name you have never seen is most of LAPACK fluency.

📜 From History: LAPACK did not appear from nowhere. Its ancestors were LINPACK (linear systems) and EISPACK (eigenproblems), Fortran libraries from the 1970s built on the newly standardized BLAS — whose Level-1 vector routines were published in 1979. Those libraries were superb for their era's machines but were built around Level-1 and Level-2 operations, and as processors grew fast caches in the 1980s, that made them memory-bound. LAPACK, first released in 1992 by a team including Jack Dongarra and Jim Demmel, was a ground-up redesign to express its algorithms in terms of Level-3 BLAS block operations so they would run near peak speed on cache-based and shared-memory machines — the very arithmetic-intensity argument from §21.2, made into a library. It has been maintained on Netlib for over thirty years, and it is the reason that "call LAPACK" is a complete answer to "how do I do linear algebra fast." When Fortran's critics call the language a relic, the software running underneath their own numerical tools is a thirty-year-old Fortran library that nobody has managed to beat.

21.5 Linking BLAS and LAPACK

A program that calls dgesv will compile without LAPACK, because the compiler treats dgesv as an external procedure to be resolved later. It will fail to link, with an error like undefined reference to 'dgesv_' (note the trailing underscore the Fortran compiler adds to the symbol). The fix is to tell the linker where the routines live:

$ gfortran -std=f2018 -Wall solve.f90 -o solve -llapack -lblas

-llapack links the reference LAPACK library; -lblas links the BLAS it depends on. Order matters on most linkers: a library must appear after the object that uses it, and since LAPACK calls the BLAS, put -llapack before -lblas. On a Linux box you typically install them first — for example sudo apt install liblapack-dev libblas-dev on Debian or Ubuntu — which places liblapack.so and libblas.so where -l can find them.

Which BLAS you link is a performance decision that does not change a line of your code. The reference BLAS (-lblas) is the plain, correct, unoptimized implementation from Netlib — fine for correctness and small problems, but it leaves most of your processor idle. A tuned BLAS implements the same interface far faster:

• OpenBLAS — a free, open-source, aggressively tuned BLAS/LAPACK; link -lopenblas (it usually provides the LAPACK symbols too, so it can replace both -llapack -lblas). The common default for good performance on commodity hardware.
• Intel MKL (Math Kernel Library) — Intel's heavily optimized BLAS/LAPACK, typically the fastest on Intel CPUs; linked through a documented set of -l flags or Intel's link-line advisor.
• Vendor libraries — Apple's Accelerate, AMD's AOCL, NVIDIA's libraries for GPUs — same interface, tuned for their silicon.

Because they all present the same BLAS/LAPACK interface, you write your dgesv call once and choose your speed at link time. Swapping reference BLAS for OpenBLAS on a large dgemm can be a tenfold difference or more, with your source untouched — the cleanest performance win in this book.

🔗 Connection: This is the ecosystem from Chapter 16 paying off. There we noted LAPACK/BLAS as the libraries "everything depends on"; here you link them. And if you build with fpm, you declare the dependency once — a link = ["lapack", "blas"] line in fpm.toml — and stop typing the flags. The Chapter 29 punchline is the same story from the performance side: the tuned library beats the loop you would write, so link the library.

🐛 Find the Bug. A colleague reports that their program "won't compile" and pastes: text /usr/bin/ld: /tmp/cc8kQ.o: in function 'MAIN__': solve.f90:(.text+0x2a1): undefined reference to 'dgesv_' collect2: error: ld returned 1 exit status They insist the code is right. Are they compiling wrong?

Answer

The code compiled fine — this is a link error, not a compile error (note ld, the linker, and undefined reference). They forgot -llapack -lblas on the command line, so the linker cannot find dgesv. The trailing underscore in dgesv_ is normal Fortran name-mangling. Add the libraries and it links.

21.6 Sparse Matrices, Briefly

Everything so far assumes a dense matrix — every one of the $n^2$ entries stored and, in general, nonzero. That assumption breaks for the largest problems. The linear system for a two-dimensional heat step on a $1000 \times 1000$ grid has a million unknowns; the dense matrix would be a million-by-million array of real(dp), which is eight terabytes — absurd — and yet almost every entry is zero, because each grid point couples only to its four neighbors. A matrix that is overwhelmingly zeros is sparse, and it demands a different representation and different solvers.

Definition (sparse format). A sparse storage format records only the nonzero entries of a matrix, plus enough index information to know where they belong, so an $n \times n$ matrix with only $O(n)$ nonzeros costs $O(n)$ memory instead of $O(n^2)$. The common formats are COO (coordinate: three arrays of row index, column index, value — simplest to build), CSR (compressed sparse row: values and column indices ordered row by row, with a pointer to where each row starts — the workhorse for matrix–vector products), and CSC (the column-major analogue, natural in Fortran). The trade is always the same: a huge saving in memory and in arithmetic, paid for with indirect indexing that is harder to write and less cache-friendly per element than a dense array.

You do not store a sparse matrix in a Fortran 2D array, and you do not solve a sparse system with dgesv — dgesv would treat all those zeros as real work and cost $O(n^3)$. Instead you reach for a library built for sparsity. Direct sparse solvers factor the matrix while trying to keep it sparse: SuperLU, UMFPACK (part of SuiteSparse), MUMPS, and Intel MKL's PARDISO are the standard names. Iterative solvers never factor at all; they only need to multiply the matrix by vectors, and reach the solution by successive approximation — conjugate gradients (CG) for symmetric positive-definite systems, GMRES for general ones, usually with a preconditioner. PETSc, a large Fortran-and-C toolkit from Argonne, wraps a great many of these behind one interface and is the common choice for serious sparse work on clusters.

For the special case that matters most to our project — a tridiagonal matrix, nonzero only on the diagonal and its two neighbors — LAPACK itself has a dedicated dense-storage routine, dgtsv (double / general tridiagonal / solve). It takes three short arrays for the sub-, main, and super-diagonals and solves in $O(n)$ time instead of $O(n^3)$. That is the honest tool for the one-dimensional implicit heat step of this chapter's Project Checkpoint — and knowing it exists, rather than throwing the full dense dgesv at a tridiagonal system, is exactly the "call the right specialist" judgment this chapter is teaching.

🔄 Check Your Understanding. 1. Why is dgesv the wrong tool for a sparse matrix with a million rows? 2. What does CSR store, and roughly how much memory does it need for an $n \times n$ matrix with five nonzeros per row?

Answers

1. dgesv is a dense solver: it stores all $n^2$ entries and does $O(n^3)$ work, treating the zeros as real. For a million rows that is terabytes of storage and an impossible operation count. 2. CSR stores the nonzero values, their column indices, and a per-row start pointer — here about $5n$ values plus $5n$ column indices plus $n{+}1$ pointers, i.e. $O(n)$ memory, not $O(n^2)$.

Project Checkpoint

This checkpoint is an optional, advanced side path — and the payoff of the whole LAPACK arc. The default way your solver marches heat forward in time is explicit stepping, built in Chapter 24: each new temperature is a simple weighted average of the old neighbors, with no linear algebra at all. Its price is the CFL limit — the timestep must stay below a stability bound, or the simulation blows up. Implicit stepping removes that limit: you may take arbitrarily large timesteps stably, but each step now requires solving a linear system. That is where LAPACK enters the solver. Here we take one implicit (backward-Euler) step of the one-dimensional heat equation and solve it with dgesv.

Discretize $\frac{\partial u}{\partial t} = \alpha\frac{\partial^2 u}{\partial x^2}$ with backward Euler in time and the three-point second difference in space. Writing $r = \alpha\,\Delta t/\Delta x^2$, the update for each interior node becomes an equation in the new temperatures:

$$ -r\,u_{i-1}^{n+1} + (1 + 2r)\,u_i^{n+1} - r\,u_{i+1}^{n+1} = u_i^{n}. $$

Collect those equations and you have a tridiagonal system $A\,\mathbf{u}^{n+1} = \mathbf{u}^{n}$, with $(1+2r)$ on the diagonal and $-r$ on the off-diagonals; fixed-temperature (Dirichlet) boundaries move onto the right-hand side. Take the smallest illustration — two interior nodes, a hot left edge $u_L = 1$, a cold right edge $u_R = 0$, a cold start $\mathbf{u}^n = (0, 0)$, and $r = 1$ (deliberately above the explicit CFL bound of $r \le \tfrac12$, to show implicit stepping is stable anyway):

program implicit_heat_step
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter :: n = 2           ! two interior nodes
  real(dp), parameter :: r = 1.0_dp      ! r = alpha*dt/dx^2  (above the explicit CFL limit)
  real(dp) :: a(n,n), rhs(n)
  integer  :: ipiv(n), info, i

  a = 0.0_dp                             ! assemble the tridiagonal A
  do i = 1, n
    a(i, i) = 1.0_dp + 2.0_dp*r          ! diagonal 1+2r
    if (i > 1) a(i, i-1) = -r            ! sub-diagonal -r
    if (i < n) a(i, i+1) = -r            ! super-diagonal -r
  end do
  rhs = [ r*1.0_dp, r*0.0_dp ]           ! u^n=0; hot-left/cold-right BCs -> [1.0, 0.0]

  call dgesv(n, 1, a, n, ipiv, rhs, n, info)   ! rhs <- u^{n+1}
  if (info /= 0) error stop 'implicit step: singular system'

  print '(a)', 'u^{n+1} after one implicit step:'
  print '(f9.4)', rhs
end program implicit_heat_step

With $r = 1$ the matrix is $\begin{bmatrix} 3 & -1 \\ -1 & 3 \end{bmatrix}$ and the right-hand side is $(1, 0)$, so (determinant $8$) the solution is $\mathbf{u}^{n+1} = (\tfrac{3}{8}, \tfrac{1}{8}) = (0.3750, 0.1250)$ — warm near the hot edge, cooler toward the cold one, exactly the physics you expect, and stable even though $r$ is twice the explicit limit. dgesv overwrites rhs with that solution:

u^{n+1} after one implicit step:
   0.3750
   0.1250

For two nodes this is overkill, and for a real one-dimensional run you would use the tridiagonal dgtsv (§21.6) rather than the dense dgesv; the point of the checkpoint is that the implicit step is a linear solve, and LAPACK does the solving. This is the LAPACK anchor at its climax: the library we named in Chapter 16, called with correct arguments, doing load-bearing work inside your own simulation. It feeds forward to Chapter 24, where the explicit path becomes the default core and this implicit path stands ready for the stiff, large-timestep runs where it earns its extra cost.

Summary

This chapter turned "the numerical libraries are Fortran" from a fact you were told into a routine you can call.

The two things to memorize: first, the dgesv calling sequence and that it overwrites both a and b — so copy them if you need them and always check info. Second, the naming scheme, because once you can read dgesv, dsyev, and dgesvd from their letters, the whole 1,700-routine library is a reference page away, not a memorization task.

Spaced Review

Three chapters ago you learned floating point; sixteen chapters ago you met the ecosystem; and the array foundations are further back still. Retrieve them.

1. (Ch. 5) Fortran stores a(3,3) column-major. Which loop nest — inner over i or inner over j — fills it with the memory grain, and why does LAPACK's expecting column-major storage follow from the same fact?
   Answer

   Inner loop over the first index i walks down a column through adjacent memory (cache-friendly); do j …; do i … is the fast order. LAPACK is Fortran, so it assumes the same first-index-fastest layout — which is why a matrix you build in Fortran needs no rearrangement to hand to dgesv, but a NumPy array does.
2. (Ch. 5) What is the difference between a * b and matmul(a, b) for two rank-2 arrays, and which one is a matrix product?
   Answer

   a * b is the elementwise product, multiplying matching positions; matmul(a, b) is the matrix product $\sum_k a_{ik}b_{kj}$. For linear algebra you want matmul (or a BLAS call), never *.
3. (Ch. 16) In the ecosystem chapter, what were LAPACK and BLAS described as, and what tool would you use to declare them as dependencies of a project?
   Answer

   The foundational numerical libraries "everything depends on," written in Fortran. fpm, the Fortran Package Manager, lets you declare the build and its library links (e.g. a link entry for lapack/blas) in fpm.toml.
4. (Ch. 16 + 20) Reference BLAS and OpenBLAS present the same interface; why can swapping them change your runtime tenfold without changing a digit of your answer — and what would change your answer's accuracy instead?
   Answer

   They implement the same mathematical operations, so results agree (to rounding); OpenBLAS is just tuned (blocking, SIMD, cache use) to run near peak. Accuracy is governed not by which BLAS you pick but by the conditioning of the problem and the precision you chose (Ch. 20) — an ill-conditioned matrix loses digits no library can restore.

What's Next

You have the workhorse of numerical linear algebra in hand: state a system in standard form, call LAPACK, check the answer. Chapter 22 turns from algebra to calculus done numerically — approximating derivatives with finite differences and integrals with the trapezoidal, Simpson, and Gaussian-quadrature rules, and measuring how fast each converges as you refine the step. It is the other half of the numerical toolkit, and it feeds directly into the differential-equation solvers of Chapters 23 and 24 — where the tridiagonal systems you just learned to solve become the beating heart of the implicit methods. Let's compute a derivative without ever writing down its formula.