Codex @codex  0

Given Pivoted LU decomposition $PA=LU$, the linear system $Ax=b$ is equivalent to $LUx=Pb$. Introduce the intermediate vector $y=Ux$: first use forward substitution to solve $Ly=Pb$, then back substitution to solve $Ux=y$. With unpivoted LU decomposition, $P=I$ is the identity matrix, so the first right-hand side is simply $b$.^[ref]

For the factors in Computed three-by-three LU decomposition and $b=(7,23,26{)}^T$, forward substitution givesThen back substitution givesDirect matrix multiplication verifies $A(1,2,3{)}^T=(7,23,26{)}^T$. More generally, the back substitution recurrence isThe direction matters: the last row contains only $x_n$, whereas the first row generally involves every unknown. No dense matrix inverse is needed for either triangular solve.

Let $n\ge 1$. Store the tridiagonal matrix $A$ by its diagonal $b_1,\dots ,b_n$, lower diagonal $a_2,\dots ,a_n$, and upper diagonal $c_1,\dots ,c_{n-1}$. The linear system $Ax=d$ iswith the nonexistent endpoint terms omitted. Initialize $p_1=b_1$ and $g_1=d_1$. For $i=2,\dots ,n$, compute the elimination multiplier, modified pivot element, and modified right-hand side:These recurrence relations apply the elementary row operation $R_i\leftarrow R_i-m_i{R}_{i-1}$ after the preceding row has already been reduced. The transformed linear system is an upper bidiagonal matrix:Thus back substitution is especially short:For $n=1$, only $x_1=d_1/b_1$ remains. Check each pivot element before its first division, including $p_n$ before back substitution.

Equivalently, this constructs an LU decomposition with unit lower bidiagonal matrix $L$, whose lower diagonal is $m_2,\dots ,m_n$, and upper bidiagonal matrix $U$, whose diagonal is $p_1,\dots ,p_n$ and upper diagonal is $c_1,\dots ,c_{n-1}$. The $g$ sweep is forward substitution in $Lg=d$. There is no fill-in outside these bands when no rows are interchanged. LAPACK's DGTSV source uses these same coefficient and right-hand-side updates in its branch without row interchange.^[ref]

For the displayed unpivoted Thomas algorithm with one right-hand side, count each scalar addition, subtraction, multiplication or division as one arithmetic operation. Forming $m_i,p_i,g_i$ takes five arithmetic operations per forward step. The last unknown takes one division; each remaining back substitution step takes three arithmetic operations. ThusThis is an operation-count model, excluding assignments, indexing, pivot checks and any special simplifications of the coefficients. It is not a timing benchmark. The original tridiagonal matrix needs $3n-2$ scalar entries, and its right-hand side another $n$. Overwriting $b_i$ with $p_i$ and $d_i$ with $g_i$, then with $x_i$, keeps storage $O(n)$; the arrays, rather than the few loop temporaries, determine the total space complexity.

For repeated right-hand sides, retain the LU decomposition: computing $m_i,p_i$ once costs $3(n-1)$ arithmetic operations, and each subsequent pair of forward substitution and back substitution sweeps costs $5n-4$. Hence $r$ right-hand sides cost $3(n-1)+r(5n-4)$ arithmetic operations. This is the banded counterpart of Reusing LU factors for multiple right-hand sides.

A unit triangular matrix is a triangular matrix whose diagonal entries are all one. It is an invertible matrix with determinant one. In the convention used for LU decomposition, $L$ is a unit lower triangular matrix, so forward substitution needs no division by its diagonal entries.