Exercises: Linear Algebra and LAPACK

These exercises live or die on getting the calling convention exactly right, so most of them want you at a keyboard with the LAPACK reference page open beside you. The habit to build is the one this chapter is really about: state the problem in standard form, hand it to the library, and check the answer — a residual you computed and confirmed is a result you own. Several problems ask you to link -llapack -lblas; if the link fails, reread §21.5 before you touch the code.

Difficulty: ⭐ warm-up · ⭐⭐ standard · ⭐⭐⭐ deeper. Solutions: worked solutions to the daggered (†) and odd-numbered problems are in appendices/answers-to-selected.md; the computational ones also appear as compilable code in code/exercise-solutions.f90. Try every problem before you look. Compile with gfortran -std=f2018 -Wall (add -llapack -lblas for anything that calls LAPACK).

Part A — Matrix Storage and the Intrinsics ⭐

21.1 † Predict the exact printout, then compile to confirm. Pay attention to the fill order:

integer :: m(2, 3)
m = reshape([1, 2, 3, 4, 5, 6], [2, 3])     ! no order= argument
print '(3i4)', m(1, :)                        ! row 1
print '(3i4)', m(2, :)                        ! row 2

Write the 2×3 matrix m represents, and explain why it is not the matrix whose rows read 1 2 3 and 4 5 6.

21.2 By hand, compute matmul(a, b) for a = [[2, 1], [0, 3]] and b = [[1, 4], [2, 1]] (each written row by row). Show all four entries. Then compute the elementwise a * b and confirm it differs.

21.3 † For a square matrix a, the product matmul(transpose(a), a) is always symmetric. Explain why in one sentence, and verify it for a = [[1, 2], [3, 4]] by computing both matmul(transpose(a), a) and its transpose and showing they match.

21.4 dot_product(x, y) returns a scalar; matmul(m, x) returns a vector when m is a matrix and x a vector. For m of shape (4, 4) and x of length 4, what is the shape of matmul(m, x)? And what is dot_product([1.0_dp, -1.0_dp, 2.0_dp], [3.0_dp, 3.0_dp, 1.0_dp])?

Part B — LAPACK Naming and the Calling Convention ⭐⭐

21.5 † Decode each LAPACK name into precision, matrix type, and computation: (a) dgesv; (b) sgetrf; (c) dsyev; (d) zheev; (e) dgels. (One computes an LU factorization; one finds eigenvalues of a complex Hermitian matrix; one solves a least-squares problem.)

21.6 Write out the full dgesv argument list in order, and for each argument state whether it is input, output, or input/output. Which two arrays does dgesv overwrite, and with what?

21.7 † You declare real(dp) :: a(100, 100) but only fill and solve the top-left 5×5 block. What value must you pass as lda to dgesv, and why is it not 5? What goes wrong if you pass 5?

21.8 A dgesv call returns info = 3. Then a second call, on a different system, returns info = -7. Interpret each value: what happened, and in which case is the bug likely in your code rather than in the data?

Part C — Port It ⭐⭐

Translate the snippet to modern Fortran that calls LAPACK, and note what the library makes you do by hand that the high-level tool did for you.

21.9 † Port this NumPy to Fortran with dgesv, and print the solution:

import numpy as np
A = np.array([[3.0, 2.0], [1.0, 2.0]])
b = np.array([7.0, 5.0])
x = np.linalg.solve(A, b)

What is x? Name three things NumPy did that your Fortran must now do explicitly.

21.10 Port w = np.linalg.eigvalsh(A) (eigenvalues of a real symmetric matrix, ascending) to a dsyev call, including the workspace query. What does eigvalsh assume about A that tells you to use dsyev rather than the general dgeev?

21.11 † Port the MATLAB one-liner x = A \ b (solve the system) to Fortran. Beyond the syntax, name the two behavioral differences you must account for: what happens to A, and what memory layout does your matrix already have that MATLAB's does not.

Part D — Find the Bug ⭐⭐

Each snippet is wrong. Say what happens — a link error, a silent wrong answer, or a crash — and fix it.

21.12 †

real(dp) :: a(3,3), b(3)
integer  :: ipiv(3), info
! ... fill a and b ...
call dgesv(3, 1, a, 3, b, ipiv, 3, info)     ! solve A x = b

21.13 A programmer solves, then checks the residual with the same a:

call dgesv(n, 1, a, n, ipiv, b, n, info)
resid = matmul(a, b) - rhs_saved             ! expecting ~0

The residual comes back large and nonsensical even though info == 0. Why?

21.14 † This is meant to solve 2x + 3y = 8, x + 4y = 9 (intended matrix [[2,3],[1,4]]), whose solution is (1, 2) — but it prints the wrong answer:

a = reshape([2.0_dp, 3.0_dp, 1.0_dp, 4.0_dp], [2, 2])   ! meant to be [[2,3],[1,4]]
b = [8.0_dp, 9.0_dp]
call dgesv(2, 1, a, 2, ipiv, b, 2, info)

What system did it actually solve (write the matrix), and what is the one-token fix?

21.15 This links and runs but sometimes prints garbage:

call dgesv(n, 1, a, n, ipiv, b, n, info)
print '(f10.4)', b        ! the solution

What critical step is missing, and what should the code do when it is present?

Part E — Design It (Extend the Solver) ⭐⭐⭐

21.16 † (Design it — solver.) Extend the Project Checkpoint from two interior nodes to n. Assemble the tridiagonal matrix a(n,n) with 1+2r on the diagonal and -r off it, build the right-hand side for a hot left edge (u_L = 1) and cold right edge (u_R = 0) starting from u^n = 0, and solve with dgesv for n = 3, r = 1. What are the three interior temperatures, and is the profile physically sensible?

21.17 (Design it — solver.) Rewrite 21.16 to call the tridiagonal solver dgtsv(n, nrhs, dl, d, du, b, ldb, info) instead of dgesv. What three short arrays replace the full a(n,n), and why is this the honest choice for a 1-D implicit step (think about the operation count as n grows)?

21.18 † (Design it.) The two-dimensional implicit heat step couples each grid point to its four neighbors. Explain why, for an nx × ny grid, the matrix A becomes (nx·ny) × (nx·ny), and why it is sparse (mostly zeros). For nx = ny = 1000, how many rows does A have, and why does this rule out dgesv?

Part F — Back of the Envelope ⭐⭐⭐

21.19 † (Cost of a solve.) LU factorization of a dense n × n matrix costs about $\tfrac{2}{3}n^3$ floating-point operations. For n = 1000, estimate the operation count. At a sustained $10^{10}$ flop/s, how long does the dgesv take? Compare to the $2n^3$ of a single matmul at the same n.

21.20 (Memory.) A dense n × n matrix of real(dp) needs $8n^2$ bytes. For n = 10{,}000, how many GiB is that ($1\ \text{GiB} = 1024^3$ bytes)? Now suppose the matrix has only 5 nonzeros per row; roughly how much memory would a sparse (CSR) representation need instead, and what factor have you saved?

21.21 † (Arithmetic intensity.) A Level-1 BLAS daxpy on length-n vectors does about 2n flops touching about 3n numbers; a Level-3 dgemm on n × n matrices does about 2n³ flops touching about 3n² numbers. Compute the flops-per-number ratio for each as a function of n, and explain in one sentence why only Level 3 can approach the processor's peak speed.

Part G — Interleaved (Chapters 5, 16, and 20) ⭐⭐

21.22 † (Ch. 5.) Fortran stores a(n, n) column-major. Write the doubly-nested loop that fills it with the memory grain, and state the one-sentence reason this same layout is why a Fortran matrix needs no transposition before you hand it to dgesv.

21.23 (Ch. 16.) In the ecosystem chapter, LAPACK and BLAS were named but not called. What does LAPACK depend on, in what order must you list them to the linker, and how would you declare the dependency in an fpm project instead of typing flags?

21.24 † (Ch. 20.) A matrix A has condition number $\kappa(A) \approx 10^{8}$. You solve A x = b in double precision, whose machine epsilon is about $10^{-16}$. Roughly how many correct significant decimal digits can you expect in x, and why does a small residual not by itself guarantee an accurate solution?

21.25 (Ch. 20.) Distinguish two failure modes of a linear solve: dgesv returning info > 0 (an exactly singular matrix) versus a matrix that is merely ill-conditioned. Which one does info catch, which one does it not, and what LAPACK routine estimates the second?

21.26 † (Synthesis.) Write a complete program that solves $A = \begin{bmatrix} 2 & 1 & 1 \\ 1 & 3 & 2 \\ 1 & 0 & 4 \end{bmatrix}$, $\mathbf{b} = (3, 2, 9)$ with dgesv, prints the solution, and confirms maxval(abs(matmul(a_orig, x) - b_orig)) < 10^{-10}. What is x, and why must you save a_orig and b_orig before the call?

21.27 (Ch. 5.) Explain why you should almost never form the inverse A⁻¹ explicitly (say, to compute x = matmul(A_inv, b)) when all you want is x. What does dgesv do instead, and which LAPACK routine would you use if you genuinely needed the inverse?

21.28 † (Design + Ch. 20.) Fit a straight line $y = c_1 + c_2 x$ to the points $(0,1), (1,3), (2,5), (3,8)$ by forming the normal equations $A^{\mathsf{T}}A\,\mathbf{c} = A^{\mathsf{T}}\mathbf{y}$ (with transpose and matmul) and solving the 2×2 system with dgesv. What are c_1 and c_2? Then say, citing conditioning (Ch. 20), why the normal equations are numerically inferior to solving the least-squares problem directly with dgels or the SVD.

Solutions to the daggered and odd-numbered problems are in appendices/answers-to-selected.md; the four computational ones (21.9, 21.16, 21.26, 21.28) are worked as runnable, LAPACK-linked code in code/exercise-solutions.f90. Design problems 21.16–21.18 admit variations — as long as your matrix is assembled correctly, your info is checked, and your residual is small, you are right.