Triangular solves for LU decomposition (source code)

= Triangular solves for LU decomposition
{c}
{title2=$Ly=Pb,\quad Ux=y$}

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$.<Fundamentals of Numerical Computation by Tobin A. Driscoll and others>[§2.4.4, equation (2.4.8) and Algorithm 2.4.2, and §2.6.3: two triangular solves and permutation of b]{ref}

For the factors in <Computed three-by-three LU decomposition> and $b=(7,23,26)^T$, <forward substitution> gives
$$
y_1=7,\qquad y_2=23-2(7)=9,\qquad y_3=26+7-3(9)=6.
$$
Then <back substitution> gives
$$
x_3=6/2=3,\qquad x_2=(9-3)/3=2,\qquad x_1=(7-2-3)/2=1.
$$
Direct <matrix multiplication> verifies $A(1,2,3)^T=(7,23,26)^T$. More generally, the <back substitution> recurrence is
$$
x_i=\frac{y_i-\sum_{j=i+1}^{n}u_{ij}x_j}{u_{ii}},\qquad i=n,n-1,\ldots,1.
$$
The 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.