◄ WORLD V · SONNY 5DART 111 · a helldive at the net

THE THOMAS ALGORITHM solve a tridiagonal system in one pass

When a linear system is tridiagonal — only the diagonal and its two neighbours are nonzero, as in every 1-D heat, spline, or finite-difference problem — you do not need full Gaussian elimination. One forward sweep and one back-substitution solve it in O(n). It is just Gauss specialised, and it is named for a man who never published it.

THE TECHNIQUE forward eliminate, then back-substitute

Sweep top to bottom, eliminating each sub-diagonal entry by combining rows (carrying modified super-diagonal and right-hand-side values); then sweep bottom to top, reading off each unknown. No pivoting, no fill-in — the three diagonals stay three diagonals. Below: a small tridiagonal system solved. live demo


  

HISTORY & CREDIT named for a report no one read

“Thomas published the tridiagonal algorithm” — he published wrote it in an unpublished internal report; the name came from someone else. cited

1949 · Llewellyn H. Thomas (Columbia’s Watson Lab) writes the method in an internal report, Elliptic Problems in Linear Difference Equations over a Network (often mis-titled “Differential”) — never formally published.
the name · David M. Young coined “Thomas algorithm” later; per W. F. Ames the method was independently discovered several times — it is simply Gaussian elimination on a tridiagonal matrix.
the caution · without pivoting it can be unstable — but for diagonally dominant systems (the usual case) it is exact and stable.
Thomas himself · is far more famous for the Thomas-Fermi atom model and for computing the fine structure of hydrogen — the tridiagonal solver was a footnote.

One of the most-used numerical kernels in the world, from a report its author never bothered to publish. Thomas 1949 (unpublished)

RECOMMEND FOR I-13 three diagonals, three 1-D arrays — sidesteps PS-004

A tridiagonal matrix is only three diagonals, so it never needs a 2-D array — and the small system solves exactly:

$ i13 run thomas.i13 # 2x0-x1=1, -x0+2x1-x2=3, -x1+2x2=3 x = [3, 5, 4] # forward sweep + back substitution, all f64
Recommend: nothing new — and it is a poster child for sidestepping PS-004: hold the sub-, main-, and super-diagonals as three separate 1-D f64 arrays (plus the right-hand side), so the missing 2-D array never bites; every step is f64 multiply / subtract / divide (verified x = [3, 5, 4]). O(n) time, O(n) space.
Note: it is Gauss / Householder (089) narrowed to a band — the same row-elimination machinery, but where a general solve wants a matrix (PS-004), the tridiagonal case runs on the arrays the corpus already has.