Case Study 1: The Solver That Divides by Zero

"The textbook algorithm and the reliable algorithm are not the same algorithm."

Executive Summary

You have inherited a small linear-system solver — twenty-odd lines of clean, readable Gaussian elimination that a predecessor wrote rather than link a library. It works. It has worked for years, on the systems it has been fed. This study reads that code, understands exactly what it does, and then hands it a system that a first-year linear-algebra student would solve in their head — and watches it produce NaN. The bug is not a typo; it is the algorithm, which omits the one step (partial pivoting) that turns textbook Gaussian elimination into a reliable one. We then replace the whole routine with a single dgesv call, confirm it solves the case that broke the original, and draw the lesson the chapter has been building toward: the reason you call LAPACK is not laziness or speed alone — it is that a specialist has already handled the failure modes you have not thought of yet.

Skills applied: reading and analyzing an existing numerical routine; the dgesv calling convention and partial pivoting (§21.3); interpreting info (§21.3); the leading dimension (§21.3); the connection between pivoting and conditioning (§21.3 and Chapter 20); the "don't ship your own" judgment (§21.2).

Background

The routine solves $A\mathbf{x} = \mathbf{b}$ by the method every numerical-methods course teaches: forward elimination to reduce $A$ to upper-triangular form, then back substitution to read off the unknowns. It is correct mathematics. Here it is, as you found it — read it and satisfy yourself that it does what it claims before we break it:

! Inherited: textbook Gaussian elimination, NO pivoting. Reads clearly; fails quietly.
subroutine naive_solve(a, b, x, n)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  integer,  intent(in)    :: n
  real(dp), intent(inout) :: a(n, n), b(n)
  real(dp), intent(out)   :: x(n)
  integer  :: k, i
  real(dp) :: m
  do k = 1, n - 1                       ! forward elimination
    do i = k + 1, n
      m = a(i, k) / a(k, k)             ! multiplier — divides by the diagonal PIVOT
      a(i, k:n) = a(i, k:n) - m * a(k, k:n)
      b(i)      = b(i)      - m * b(k)
    end do
  end do
  do i = n, 1, -1                       ! back substitution
    x(i) = (b(i) - dot_product(a(i, i+1:n), x(i+1:n))) / a(i, i)
  end do
end subroutine naive_solve

The whole story is in one line: m = a(i, k) / a(k, k). The algorithm divides by the diagonal element a(k, k) — the pivot — and never asks whether that is a safe thing to do.

Phase 1 — Confirm It Works on a Benign System

Never study a bug without first confirming the code works when it should. Feed it a well-behaved 2×2 system,

$$ \begin{bmatrix} 2 & 1 \\ 1 & 3 \end{bmatrix}\mathbf{x} = \begin{bmatrix} 4 \\ 7 \end{bmatrix}, \qquad \text{solution } \mathbf{x} = (1, 2). $$

Trace it by hand: the first pivot is a(1,1) = 2, comfortably nonzero. The multiplier is $m = 1/2$; row 2 becomes $[0,\ 2.5]$ with right-hand side $7 - 0.5\cdot4 = 5$; back substitution gives $x_2 = 5/2.5 = 2$ and $x_1 = (4 - 1\cdot2)/2 = 1$. The routine returns $(1, 2)$, exactly right. So far the inheritance looks sound, and you can see why nobody has touched it: on ordinary systems, it is correct.

Phase 2 — The System That Breaks It

Now a system that is, if anything, easier — you can solve it by inspection:

$$ \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}\mathbf{x} = \begin{bmatrix} 1 \\ 2 \end{bmatrix}. $$

The first equation says $y = 1$; the second then says $x + 1 = 2$, so $x = 1$. The answer is $(1, 1)$, and there is nothing pathological about the system — its determinant is $-1$, it is perfectly well-conditioned, it has a unique solution. Hand it to naive_solve and watch the first line detonate:

program break_it
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: a(2,2), b(2), x(2)
  a = reshape([ 0.0_dp, 1.0_dp, 1.0_dp, 1.0_dp ], [2, 2], order=[2, 1])
  b = [ 1.0_dp, 2.0_dp ]
  call naive_solve(a, b, x, 2)
  print *, 'naive_solve gives x =', x
end program break_it

The very first pivot is a(1,1) = 0. The multiplier m = a(2,1) / a(1,1) is $1 / 0$, which in IEEE arithmetic is $+\infty$ (gfortran does not trap it by default). From there the poison spreads: $\infty \cdot 0$ is `NaN`, the eliminated row fills with `NaN` and $-\infty$, and back substitution divides one indefinite quantity by another. The program prints not the answer but its wreckage:

 naive_solve gives x =                       NaN                       NaN

(The exact spacing of NaN is compiler-specific; the point is that it is not a number.) The routine did not warn you, did not stop, did not return an error code. It computed NaN with total confidence and handed it back as if it were an answer. In a larger program that value would flow downstream and quietly corrupt everything it touched. This is the failure mode the library exists to prevent.

Phase 3 — The One Missing Idea: Partial Pivoting

The system is fine; the algorithm is incomplete. Textbook Gaussian elimination assumes the pivot is nonzero; reliable Gaussian elimination guarantees it by reordering. Before eliminating with column $k$, partial pivoting searches that column for its largest-magnitude entry and swaps that row up into the pivot position. For our system, column 1 is $(0, 1)$; the larger magnitude is the $1$ in row 2, so pivoting swaps the rows:

$$ \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix}\mathbf{x} = \begin{bmatrix} 2 \\ 1 \end{bmatrix}, $$

and now the pivot is $1$, the elimination is trivial, and back substitution gives $(1, 1)$ — the right answer. Partial pivoting costs almost nothing (a search and a row swap per column) and it is the difference between a solver that works on the systems it happens to be given and one that works on every nonsingular system. It is precisely the step naive_solve omits, and precisely the step LAPACK's dgesv performs automatically, recording the swaps in its ipiv array.

Phase 4 — Replace the Routine with dgesv

You could add pivoting to naive_solve — it is a good exercise. You should not ship it. The robust, tested, pivoting-by-default solver already exists and is one call away. Here is the replacement, on the exact system that produced NaN:

program robust_solve
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer, parameter :: n = 2
  real(dp) :: a(n,n), b(n)
  integer  :: ipiv(n), info

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

  call dgesv(n, 1, a, n, ipiv, b, n, info)      ! partial pivoting handled internally
  if (info /= 0) then
    print '(a, i0)', 'dgesv reported info = ', info
    error stop 1
  end if
  print '(a)', 'x ='
  print '(f8.4)', b
end program robust_solve

$ gfortran -std=f2018 -Wall robust_solve.f90 -o robust -llapack -lblas && ./robust
x =
  1.0000
  1.0000

The system that reduced the inherited routine to NaN, dgesv solves without comment, because pivoting is not an optional extra to LAPACK — it is built in. And note the discipline the library enforces: it hands you an info you are expected to check. Had the matrix been genuinely singular, info would have come back positive and you would have known, rather than discovering it three modules downstream when a plot came out wrong.

Phase 5 — The Deeper Danger: Pivots That Are Small, Not Zero

A zero pivot at least fails loudly, in NaN. The subtler menace is a pivot that is merely tiny. Consider

$$ \begin{bmatrix} 10^{-20} & 1 \\ 1 & 1 \end{bmatrix}\mathbf{x} = \begin{bmatrix} 1 \\ 2 \end{bmatrix}, $$

whose true solution is very close to $(1, 1)$. naive_solve does not divide by zero here — the pivot is $10^{-20}$, not $0$ — so it runs to completion and returns a number. That is worse. The multiplier becomes $1/10^{-20} = 10^{20}$, an enormous factor that, when it multiplies the pivot row and subtracts, utterly swamps the modest entries carrying the actual answer: computing $1 - 10^{20}$ in double precision simply loses the $1$, because it falls below the rounding resolution of $10^{20}$. This is catastrophic cancellation, the hazard of Chapter 20, and it leaves the routine confidently reporting a badly wrong $x$ — often something near $0$ instead of $1$. Partial pivoting avoids the whole disaster by refusing to use a tiny pivot when a larger one is available in the column: dgesv swaps the rows, uses the $1$ as the pivot, keeps every multiplier at or below $1$ in magnitude, and returns the correct $(1, 1)$.

The pattern is the argument of the whole chapter, made concrete. The inherited code is not wrong in any line; it faithfully implements the algorithm as it appears in a textbook. It is incomplete in the way that real numerical software cannot afford to be, and the missing piece — pivoting for stability — is exactly the kind of hard-won robustness that decades of LAPACK development have baked in. You call the library not because you could not write elimination, but because the version worth trusting is not the one in the textbook.

Discussion Questions

1. naive_solve divides by a(k, k) with no check. Add a guard that calls error stop if abs(a(k,k)) is below a small threshold. Does that make the routine correct, or merely loud? What can pivoting do that a guard cannot?
2. The zero-pivot system is perfectly well-conditioned — its answer is not sensitive to small perturbations — yet the naive routine fails on it. Explain, in one sentence, why conditioning (a property of the matrix) and stability (a property of the algorithm) are different things.
3. dgesv records its row swaps in ipiv but you never read it. Sketch what the routine would have to return to you instead if it did not keep that array — and why keeping it is cheaper than the alternative.

Your Turn: Extensions

• Option A. Add partial pivoting to naive_solve yourself: before each elimination step, find maxloc(abs(a(k:n, k))), swap that row up in both a and b, then proceed. Test it on all three systems in the table and confirm it now matches dgesv. (This is the best possible way to appreciate what the library does for free.)
• Option B. Instrument the tiny-pivot case: solve the $10^{-20}$ system with both routines and print the residual $\lVert A\mathbf{x} - \mathbf{b}\rVert$ from each (using a saved copy of $A$). Which routine's residual betrays the wrong answer, and which hides it?
• Option C. Replace the dgesv call with the fortran-lang stdlib solve wrapper (stdlib_linalg), which type-checks its arguments. Note how much of the boilerplate — the ipiv array, the info check, the copy of a — disappears, and decide when you would still want the raw dgesv.

Key Takeaways

• The reason to call LAPACK is not only speed; it is robustness you did not write — partial pivoting, the step that separates textbook Gaussian elimination from a solver you can trust on every nonsingular system.
• A zero pivot fails loudly (NaN); a tiny pivot fails silently (catastrophic cancellation, Chapter 20), which is worse. Pivoting prevents both by never using a small pivot when a larger one is available.
• dgesv performs pivoting automatically, records it in ipiv, and forces you to receive an info status — turning a silent failure into a checkable one.
• Inherited numerical code can be perfectly readable and still be subtly incomplete. When you find a hand-rolled solver, the right move is usually not to fix it but to replace it with the library and delete it.