Tridiagonal elimination recurrences (source code)

= Tridiagonal elimination recurrences
{title2=$m_i=a_i/p_{i-1},\quad p_i=b_i-m_ic_{i-1}$}

Let $n\geq1$. Store the <tridiagonal matrix> $A$ by its diagonal $b_1,\ldots,b_n$, lower diagonal $a_2,\ldots,a_n$, and upper diagonal $c_1,\ldots,c_{n-1}$. The <linear system> $Ax=d$ is
$$
a_ix_{i-1}+b_ix_i+c_ix_{i+1}=d_i,\qquad 1\leq i\leq n,
$$
with the nonexistent endpoint terms omitted. Initialize $p_1=b_1$ and $g_1=d_1$. For $i=2,\ldots,n$, compute the <elimination multiplier>, modified <pivot element>, and modified right-hand side:
$$
m_i=\frac{a_i}{p_{i-1}},\qquad p_i=b_i-m_ic_{i-1},\qquad g_i=d_i-m_ig_{i-1}.
$$
These <recurrence relations> apply the <elementary row operation> $R_i\leftarrow R_i-m_iR_{i-1}$ after the preceding row has already been reduced. The transformed <linear system> is an upper <bidiagonal matrix>:
$$
p_ix_i+c_ix_{i+1}=g_i\quad(1\leq i<n),\qquad p_nx_n=g_n.
$$
Thus <back substitution> is especially short:
$$
x_n=\frac{g_n}{p_n},\qquad x_i=\frac{g_i-c_ix_{i+1}}{p_i}\quad(i=n-1,n-2,\ldots,1).
$$
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>.

\Image[/tridiagonal-elimination-band-structure.png]
{title=Tridiagonal elimination preserves the narrow band: only one row and one diagonal coefficient change per step}

Equivalently, this constructs an <LU decomposition> with unit lower <bidiagonal matrix> $L$, whose lower diagonal is $m_2,\ldots,m_n$, and upper <bidiagonal matrix> $U$, whose diagonal is $p_1,\ldots,p_n$ and upper diagonal is $c_1,\ldots,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.<DGTSV (LAPACK 3.12.1) by LAPACK project>[DGTSV, LAPACK 3.12.1, Fortran lines 175–180: multiplier and two updates without row interchange]{ref}