Pivoted LU decomposition (source code)

= Pivoted LU decomposition
{c}
{title2=$PA=LU$}

We fix the convention $PA=LU$: the <permutation matrix> $P$ acts on the rows of $A$, and the same row ordering must be applied to $b$ before the <Triangular solves for LU decomposition>. Equivalently, $A=P^TLU$, since $P^{-1}=P^T$. LAPACK documents its factors as $A=PLU$ instead, so check the defining equation rather than just the letter $P$.<DGETRS (LAPACK 3.3.1) by LAPACK project>[LAPACK DGETRS, A argument: factors from DGETRF satisfy $A=PLU$]{ref} <Partial pivoting> selects the largest absolute candidate in the current column, among rows $k,\ldots,n$, before forming <elimination multipliers>.<Fundamentals of Numerical Computation by Tobin A. Driscoll and others>[§2.6.1, Definition 2.6.1; §2.6.2, Definition 2.6.2 and equation (2.6.1): pivot selection and reordered-row factorization]{ref}

A zero first <pivot element> does not imply a <singular matrix>. For example,
$$
A_0=\begin{pmatrix}0&1\\2&3\end{pmatrix},\qquad \det A_0=-2,\qquad P=\begin{pmatrix}0&1\\1&0\end{pmatrix},\qquad PA_0=\underbrace{\begin{pmatrix}1&0\\0&1\end{pmatrix}}_{L}\underbrace{\begin{pmatrix}2&3\\0&1\end{pmatrix}}_{U}.
$$
Unpivoted <Gaussian elimination> would divide $2$ by zero. Indeed, no $A_0=LU$ with a unit <lower triangular matrix> $L$ exists: $u_{11}=0$ would force $(LU)_{21}=l_{21}u_{11}=0$, contradicting $(A_0)_{21}=2$. One row interchange resolves the obstruction. For $b_0=(1,5)^T$, use $Pb_0=(5,1)^T$; <back substitution> then gives $x=(1,1)^T$. Using the original $b_0$ with these factors would solve a different <linear system>.

<Partial pivoting> bounds each <elimination multiplier> by one in absolute value in exact arithmetic. This helps control <roundoff error>, but intermediate entries can still grow; it is not an unconditional guarantee of <numerical stability> for every <matrix>.<Fundamentals of Numerical Computation by Tobin A. Driscoll and others>[§2.6.4, Example 2.6.6 and the concluding discussion: stability of unpivoted and pivoted factorization]{ref} It also does not remove the original problem's <condition number>. For the <diagonal matrix> $D=\operatorname{diag}(1,\varepsilon)$ with $0<\varepsilon<1$, the <condition number> in the <Euclidean norm> is $\kappa_2(D)=1/\varepsilon$. Row permutations preserve this <condition number>, because they preserve <Euclidean norms>. Even an exact solve changes from $(1,1)^T$ to $(1,1+\eta/\varepsilon)^T$ when the right-hand side changes from $(1,\varepsilon)^T$ to $(1,\varepsilon+\eta)^T$. This sensitivity is a property of the <linear system>, not a zero-<pivot element> failure.<Fundamentals of Numerical Computation by Tobin A. Driscoll and others>[§2.8, Definition 2.8.1, and §2.8.1, Theorem 2.8.1, equations (2.8.5)–(2.8.6): condition number and right-hand-side sensitivity]{ref}