LU Decomposition

Why LU Decomposition Matters

Imagine a circuit model whose connections and component values stay fixed while the input voltages change. Each experiment produces a new right-hand side in $Ax=b$, but the coefficient matrix $A$ stays the same. Repeating all of Gaussian elimination for every experiment wastes work. LU decomposition saves the elimination work in two triangular matrices, so each later solve needs only two substitution passes.

The letters describe the shapes: $L$ is lower triangular, and $U$ is upper triangular. Row swaps make the method work when a diagonal candidate is zero or inconveniently small. We will record those swaps in a permutation $P$ and use the convention $PA=LU$. The same factors will solve two different systems, and independent residuals will check both the factorization and the answers.

By the end, you should be able to explain what a stored elimination multiplier means, trace a pivot that moves earlier multipliers, implement a reusable factor-and-solve interface in three languages, and distinguish numerical rejection from an exact proof of singularity.

Core Concepts and Terminology

The implementation accepts a nonempty finite square matrix $A\in\mathbb{R}^{n\times n}$ and a relative pivot tolerance $0\lt\tau\lt1$. On success it returns a unit lower triangular matrix $L$, an upper triangular matrix $U$, and a row-order array $p$. It targets matrices whose pivots pass the chosen tolerance; it rejects singular matrices and may also reject nonsingular matrices with very small pivots relative to the input scale.

Use zero-based indices throughout. Define $P$ by $P_{i,p_i}=1$, with all other entries zero. Consequently, row $i$ of $PA$ is original row $p_i$ of $A$, and $(Pb)_i=b_{p_i}$. This array is the final row ordering, rather than a log of individual swaps. The code applies it directly without constructing a dense permutation matrix.

Triangular shape makes substitution possible. In $Ly=c$, row $i$ contains only the already-known values $y_0,\ldots,y_{i-1}$ and the new value $y_i$, whose coefficient is $1$. In $Ux=y$, start at the last row and work upward, dividing by each nonzero diagonal entry of $U$.

The LAPACK DGETRF reference documents partial pivoting, unit lower triangular factors, and a production implementation based on matrix blocks. Its displayed convention is $A=P_{\mathrm{LAPACK}}LU$; our convention uses $P=P_{\mathrm{LAPACK}}^{-1}$, so check the direction of a library's permutation before copying its solve formula. LAPACK's pivot indices describe sequential interchanges, unlike the final ordering returned here.

How LU Decomposition Works, Step by Step

Start with $U=A$, $L=I$, and $p=(0,1,\ldots,n-1)$. Keep the original scale $\alpha=\max_{i,j}\lvert A_{ij}\rvert$ fixed throughout factorization. Reject $\alpha=0$ immediately. At the beginning of stage $k$, the columns before $k$ already have zeros below their diagonal in $U$.

1. Select a pivot. Search rows $k,\ldots,n-1$ for a maximum of $\lvert U_{ik}\rvert$. Keep the first row on a tie, so the implementations agree. If the chosen pivot satisfies $\lvert U_{rk}\rvert/\alpha\le\tau$, stop with an error before dividing.
2. Move the equation and its history. If $r\ne k$, swap entire rows $k$ and $r$ of $U$, swap $p_k$ and $p_r$, and swap only entries in columns $0,\ldots,k-1$ between rows $k$ and $r$ of $L$. Those entries contain previously recorded multipliers. Leave the diagonal and not-yet-filled columns of $L$ in place.
3. Store each multiplier. For every $i\gt k$, set $L_{ik}=U_{ik}/U_{kk}$. Save this ratio before setting $U_{ik}=0$.
4. Update the trailing row. For every $j\gt k$, replace $U_{ij}$ by $U_{ij}-L_{ik}U_{kj}$. The pivot row is unchanged. In this Doolittle-style arrangement, we do not normalize its diagonal to $1$.
5. Advance. Repeat for $k+1$. The final stage still checks its pivot even though no rows remain below it. Return the completed factors and permutation.

Once the factors exist, a solve uses

$$ y_i=b_{p_i}-\sum_{j=0}^{i-1}L_{ij}y_j, \qquad x_i=\frac{y_i-\sum_{j=i+1}^{n-1}U_{ij}x_j}{U_{ii}}. $$

Evaluate the first recurrence from top to bottom and the second from bottom to top. Empty sums are zero. A later right-hand side reuses the same $L$, $U$, and $p$; only $y$ and $x$ are recomputed.

flowchart TD
    A["Validate A and tolerance; initialize L, U, p"] --> B{"More pivot columns?"}
    B -->|Yes| C["Choose largest active column entry"]
    C --> D{"Pivot passes relative test?"}
    D -->|No| E["Report rejected pivot"]
    D -->|Yes| F["Swap U rows, p entries, and earlier L entries"]
    F --> G["Store multipliers; update trailing submatrix"]
    G --> H["Advance pivot column"]
    H --> B
    B -->|No| I["Keep factors for reuse"]
    I --> J["Permute b; forward solve L; back solve U"]
    J --> K{"Another right-hand side?"}
    K -->|Yes| J
    K -->|No| M["Check original-system residuals"]

Worked Example with a Full Trace

Use

$$ A=\begin{pmatrix}0&2&1\\1&1&0\\2&1&1\end{pmatrix}, \qquad b=\begin{pmatrix}7\\3\\7\end{pmatrix}. $$

Initialization. Set $U=A$, $L=I_3$, $p=(0,1,2)$, and $\alpha=2$. The right-hand side remains untouched during factorization.

Stage $k=0$: first row swap and two eliminations. The candidate magnitudes are $(0,1,2)$, so choose row $2$. Swap rows $0$ and $2$ of $U$, and update $p=(2,1,0)$. There are no earlier columns of $L$ to swap. At this point,

$$ U=\begin{pmatrix}2&1&1\\1&1&0\\0&2&1\end{pmatrix}. $$

The pivot ratio is $2/2=1$, which passes the tolerance. For row $1$, store $m_{10}=1/2$ in $L_{10}$ and subtract half of row $0$: $(1,1,0)-(1/2)(2,1,1)=(0,1/2,-1/2)$. For row $2$, store $m_{20}=0$; its entries stay $(0,2,1)$. The completed first stage is

$$ L=\begin{pmatrix}1&0&0\\1/2&1&0\\0&0&1\end{pmatrix}, \qquad U=\begin{pmatrix}2&1&1\\0&1/2&-1/2\\0&2&1\end{pmatrix}, \qquad p=(2,1,0). $$

Stage $k=1$: move an earlier multiplier. The active column entries are $1/2$ and $2$, so swap rows $1$ and $2$ of $U$. Update $p=(2,0,1)$. Crucially, also swap $L_{10}=1/2$ with $L_{20}=0$; leave both diagonal ones alone. Immediately after this swap,

$$ L=\begin{pmatrix}1&0&0\\0&1&0\\1/2&0&1\end{pmatrix}, \qquad U=\begin{pmatrix}2&1&1\\0&2&1\\0&1/2&-1/2\end{pmatrix}. $$

The pivot ratio is again $1$. Store $m_{21}=(1/2)/2=1/4$. Row $2$ becomes $(0,1/2,-1/2)-(1/4)(0,2,1)=(0,0,-3/4)$.

Stage $k=2$: final pivot check. The remaining pivot is $-3/4$. Its relative magnitude is $(3/4)/2=3/8\gt10^{-12}$, so it passes. No row remains to eliminate. The final factors are

$$ L=\begin{pmatrix}1&0&0\\0&1&0\\1/2&1/4&1\end{pmatrix}, \qquad U=\begin{pmatrix}2&1&1\\0&2&1\\0&0&-3/4\end{pmatrix}. $$

Check the factorization directly:

$$ LU=\begin{pmatrix}2&1&1\\0&2&1\\1&1&0\end{pmatrix} =PA, \qquad P=\begin{pmatrix}0&0&1\\1&0&0\\0&1&0\end{pmatrix}. $$

Forward substitution. Permute the original vector once: $Pb=(7,7,3)^T$. Then $y_0=7$, $y_1=7$, and $y_2=3-(1/2)7-(1/4)7=-9/4$.

Back substitution. Compute $x_2=(-9/4)/(-3/4)=3$, then $x_1=(7-3)/2=2$, and finally $x_0=(7-2-3)/2=1$. Hence $x=(1,2,3)^T$, and direct multiplication gives $Ax=(7,3,7)^T=b$.

Reuse the factors. For $b^{(2)}=(0,1,2)^T$, the permuted vector is $(2,0,1)^T$, the forward result is $y^{(2)}=(2,0,0)^T$, and the back solve gives $x^{(2)}=(1,0,0)^T$. No elimination is repeated. Two row swaps also give $\det(A)=(-1)^2(2)(2)(-3/4)=-3$.

Implementation in Python, C++, and Java

Each panel is a complete program using ordinary nested arrays, explicit partial pivoting, and the same sample data. Run Python 3 directly; compile the C++ program as C++14 or newer; save the Java program as Main.java. The factorization copies $A$ and returns separate $L$, $U$, and $p$ arrays. The solve function accepts an unchanged factorization result and preserves its factors and the input vector.

The default tolerance is $\tau=10^{-12}$. The pivot test compares against the original maximum entry magnitude, so uniformly scaling $A$ leaves the test unchanged in exact arithmetic. This simple threshold is a policy for rejecting small pivots, not a condition-number estimate or a general numerical-rank algorithm. Equilibration or a different tolerance may be appropriate when row and column scales differ greatly. The programs reject invalid dimensions, nonfinite data, rejected pivots, and detected nonfinite arithmetic; they assume real values representable as ordinary floating-point numbers.

The factor error is $\max_{i,j}\lvert(LU)_{ij}-A_{p_i,j}\rvert$. The solve error is the maximum absolute component of $Ax-b$, computed against the original input. These are absolute residuals. The sample uses an absolute pass threshold of $10^{-10}$ because its data have modest scale; the threshold is not a universal accuracy criterion. Printing converts signed zero to ordinary zero for consistent output without changing computations.

import math


def lu_factor(a, tol=1e-12):
    if not math.isfinite(tol) or not 0.0 < tol < 1.0:
        raise ValueError("tol must be finite and between 0 and 1")
    n = len(a)
    if n == 0 or any(len(row) != n for row in a):
        raise ValueError("A must be a nonempty square matrix")
    if any(not math.isfinite(v) for row in a for v in row):
        raise ValueError("A must contain finite values")
    u = [list(map(float, row)) for row in a]
    l = [[float(i == j) for j in range(n)] for i in range(n)]
    p = list(range(n))
    scale = max(abs(v) for row in u for v in row)
    if scale == 0.0:
        raise ValueError("zero matrix is singular")
    for k in range(n):
        pivot = k
        for i in range(k + 1, n):
            if abs(u[i][k]) > abs(u[pivot][k]):
                pivot = i
        if abs(u[pivot][k]) / scale <= tol:
            raise ValueError("pivot rejected at the chosen tolerance")
        if pivot != k:
            u[k], u[pivot] = u[pivot], u[k]
            p[k], p[pivot] = p[pivot], p[k]
            for j in range(k):
                l[k][j], l[pivot][j] = l[pivot][j], l[k][j]
        for i in range(k + 1, n):
            l[i][k] = u[i][k] / u[k][k]
            u[i][k] = 0.0
            for j in range(k + 1, n):
                u[i][j] -= l[i][k] * u[k][j]
                if not math.isfinite(u[i][j]):
                    raise ValueError("nonfinite factorization arithmetic")
    return l, u, p


def lu_solve(factors, b):
    # factors must be an unchanged result of lu_factor.
    l, u, p = factors
    n = len(p)
    if len(b) != n or any(not math.isfinite(v) for v in b):
        raise ValueError("b must have n finite entries")
    y = [0.0] * n
    for i in range(n):
        value = b[p[i]]
        for j in range(i):
            value -= l[i][j] * y[j]
        y[i] = value  # L has diagonal entries 1.
        if not math.isfinite(y[i]):
            raise ValueError("nonfinite forward substitution")
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        value = y[i]
        for j in range(i + 1, n):
            value -= u[i][j] * x[j]
        x[i] = value / u[i][i]
        if not math.isfinite(x[i]):
            raise ValueError("nonfinite back substitution")
    return x


def factor_error(a, factors):
    l, u, p = factors
    n = len(p)
    return max(abs(sum(l[i][k] * u[k][j] for k in range(n))
                   - a[p[i]][j]) for i in range(n) for j in range(n))


def solve_error(a, x, b):
    return max(abs(sum(row[j] * x[j] for j in range(len(x))) - b[i])
               for i, row in enumerate(a))


def vector_text(v):
    return "[" + ", ".join(f"{(0.0 if value == 0.0 else value):.6f}" for value in v) + "]"


def matrix_text(a):
    return "[" + ", ".join(vector_text(row) for row in a) + "]"


def main():
    a = [[0.0, 2.0, 1.0], [1.0, 1.0, 0.0], [2.0, 1.0, 1.0]]
    b = [7.0, 3.0, 7.0]
    b2 = [0.0, 1.0, 2.0]
    factors = lu_factor(a)
    l, u, p = factors
    x = lu_solve(factors, b)
    x2 = lu_solve(factors, b2)
    fe = factor_error(a, factors)
    se = solve_error(a, x, b)
    re = solve_error(a, x2, b2)
    print("permutation = " + str(p))
    print("L = " + matrix_text(l))
    print("U = " + matrix_text(u))
    print("x = " + vector_text(x))
    print("reused x = " + vector_text(x2))
    print(f"factor residual = {fe:.3e}")
    print(f"solve residual = {se:.3e}")
    print(f"reuse residual = {re:.3e}")
    print("check = " + ("PASS" if max(fe, se, re) < 1e-10 else "FAIL"))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>

using Matrix = std::vector<std::vector<double>>;
struct Factors {
    Matrix l, u;
    std::vector<std::size_t> p;
};

Factors lu_factor(const Matrix& a, double tol = 1e-12) {
    if (!std::isfinite(tol) || !(tol > 0.0 && tol < 1.0))
        throw std::invalid_argument("tol must be finite and between 0 and 1");
    const std::size_t n = a.size();
    if (n == 0) throw std::invalid_argument("A must be nonempty");
    double scale = 0.0;
    for (const auto& row : a) {
        if (row.size() != n) throw std::invalid_argument("A must be square");
        for (double v : row) {
            if (!std::isfinite(v)) throw std::invalid_argument("A must be finite");
            scale = std::max(scale, std::abs(v));
        }
    }
    if (scale == 0.0) throw std::runtime_error("zero matrix is singular");
    Factors f{Matrix(n, std::vector<double>(n, 0.0)), a,
              std::vector<std::size_t>(n)};
    std::iota(f.p.begin(), f.p.end(), 0);
    for (std::size_t i = 0; i < n; ++i) f.l[i][i] = 1.0;
    for (std::size_t k = 0; k < n; ++k) {
        std::size_t pivot = k;
        for (std::size_t i = k + 1; i < n; ++i)
            if (std::abs(f.u[i][k]) > std::abs(f.u[pivot][k])) pivot = i;
        if (std::abs(f.u[pivot][k]) / scale <= tol)
            throw std::runtime_error("pivot rejected at the chosen tolerance");
        if (pivot != k) {
            std::swap(f.u[k], f.u[pivot]);
            std::swap(f.p[k], f.p[pivot]);
            for (std::size_t j = 0; j < k; ++j)
                std::swap(f.l[k][j], f.l[pivot][j]);
        }
        for (std::size_t i = k + 1; i < n; ++i) {
            f.l[i][k] = f.u[i][k] / f.u[k][k];
            f.u[i][k] = 0.0;
            for (std::size_t j = k + 1; j < n; ++j) {
                f.u[i][j] -= f.l[i][k] * f.u[k][j];
                if (!std::isfinite(f.u[i][j]))
                    throw std::runtime_error("nonfinite factorization arithmetic");
            }
        }
    }
    return f;
}

std::vector<double> lu_solve(const Factors& f, const std::vector<double>& b) {
    // f must be an unchanged result of lu_factor.
    const std::size_t n = f.p.size();
    if (b.size() != n) throw std::invalid_argument("b must have n entries");
    for (double v : b)
        if (!std::isfinite(v)) throw std::invalid_argument("b must be finite");
    std::vector<double> y(n), x(n);
    for (std::size_t i = 0; i < n; ++i) {
        double value = b[f.p[i]];
        for (std::size_t j = 0; j < i; ++j) value -= f.l[i][j] * y[j];
        y[i] = value;
        if (!std::isfinite(y[i])) throw std::runtime_error("nonfinite forward substitution");
    }
    for (std::size_t i = n; i-- > 0;) {
        double value = y[i];
        for (std::size_t j = i + 1; j < n; ++j) value -= f.u[i][j] * x[j];
        x[i] = value / f.u[i][i];
        if (!std::isfinite(x[i])) throw std::runtime_error("nonfinite back substitution");
    }
    return x;
}

double factor_error(const Matrix& a, const Factors& f) {
    double error = 0.0;
    const std::size_t n = a.size();
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) {
            double value = 0.0;
            for (std::size_t k = 0; k < n; ++k) value += f.l[i][k] * f.u[k][j];
            error = std::max(error, std::abs(value - a[f.p[i]][j]));
        }
    return error;
}

double solve_error(const Matrix& a, const std::vector<double>& x,
                   const std::vector<double>& b) {
    double error = 0.0;
    for (std::size_t i = 0; i < a.size(); ++i) {
        double value = 0.0;
        for (std::size_t j = 0; j < x.size(); ++j) value += a[i][j] * x[j];
        error = std::max(error, std::abs(value - b[i]));
    }
    return error;
}

std::string vector_text(const std::vector<double>& v) {
    std::ostringstream s;
    s << std::fixed << std::setprecision(6) << '[';
    for (std::size_t i = 0; i < v.size(); ++i) {
        if (i) s << ", ";
        s << (v[i] == 0.0 ? 0.0 : v[i]);
    }
    s << ']';
    return s.str();
}

std::string matrix_text(const Matrix& a) {
    std::string s = "[";
    for (std::size_t i = 0; i < a.size(); ++i) {
        if (i) s += ", ";
        s += vector_text(a[i]);
    }
    return s + "]";
}

int main() {
    const Matrix a = { {0.0, 2.0, 1.0}, {1.0, 1.0, 0.0}, {2.0, 1.0, 1.0} };
    const std::vector<double> b = {7.0, 3.0, 7.0}, b2 = {0.0, 1.0, 2.0};
    const auto f = lu_factor(a);
    const auto x = lu_solve(f, b), x2 = lu_solve(f, b2);
    const double fe = factor_error(a, f), se = solve_error(a, x, b),
                 re = solve_error(a, x2, b2);
    std::cout << "permutation = [";
    for (std::size_t i = 0; i < f.p.size(); ++i) {
        if (i) std::cout << ", ";
        std::cout << f.p[i];
    }
    std::cout << "]\nL = " << matrix_text(f.l) << "\nU = " << matrix_text(f.u)
              << "\nx = " << vector_text(x) << "\nreused x = " << vector_text(x2)
              << std::scientific << std::setprecision(3)
              << "\nfactor residual = " << fe << "\nsolve residual = " << se
              << "\nreuse residual = " << re
              << "\ncheck = " << (std::max({fe, se, re}) < 1e-10 ? "PASS" : "FAIL") << '\n';
}

import java.util.Arrays;
import java.util.Locale;

public class Main {
    static class Factors {
        final double[][] l, u;
        final int[] p;
        Factors(double[][] l, double[][] u, int[] p) {
            this.l = l; this.u = u; this.p = p;
        }
    }

    static Factors luFactor(double[][] a, double tol) {
        if (!Double.isFinite(tol) || !(tol > 0.0 && tol < 1.0))
            throw new IllegalArgumentException("tol must be finite and between 0 and 1");
        int n = a.length;
        if (n == 0) throw new IllegalArgumentException("A must be nonempty");
        double scale = 0.0;
        for (double[] row : a) {
            if (row == null || row.length != n)
                throw new IllegalArgumentException("A must be square");
            for (double v : row) {
                if (!Double.isFinite(v)) throw new IllegalArgumentException("A must be finite");
                scale = Math.max(scale, Math.abs(v));
            }
        }
        if (scale == 0.0) throw new IllegalArgumentException("zero matrix is singular");
        double[][] l = new double[n][n], u = new double[n][];
        int[] p = new int[n];
        for (int i = 0; i < n; ++i) {
            l[i][i] = 1.0; u[i] = a[i].clone(); p[i] = i;
        }
        for (int k = 0; k < n; ++k) {
            int pivot = k;
            for (int i = k + 1; i < n; ++i)
                if (Math.abs(u[i][k]) > Math.abs(u[pivot][k])) pivot = i;
            if (Math.abs(u[pivot][k]) / scale <= tol)
                throw new IllegalArgumentException("pivot rejected at the chosen tolerance");
            if (pivot != k) {
                double[] row = u[k]; u[k] = u[pivot]; u[pivot] = row;
                int index = p[k]; p[k] = p[pivot]; p[pivot] = index;
                for (int j = 0; j < k; ++j) {
                    double value = l[k][j]; l[k][j] = l[pivot][j]; l[pivot][j] = value;
                }
            }
            for (int i = k + 1; i < n; ++i) {
                l[i][k] = u[i][k] / u[k][k];
                u[i][k] = 0.0;
                for (int j = k + 1; j < n; ++j) {
                    u[i][j] -= l[i][k] * u[k][j];
                    if (!Double.isFinite(u[i][j]))
                        throw new IllegalArgumentException("nonfinite factorization arithmetic");
                }
            }
        }
        return new Factors(l, u, p);
    }

    static double[] luSolve(Factors f, double[] b) {
        // f must be an unchanged result of luFactor.
        int n = f.p.length;
        if (b.length != n) throw new IllegalArgumentException("b must have n entries");
        for (double v : b)
            if (!Double.isFinite(v)) throw new IllegalArgumentException("b must be finite");
        double[] y = new double[n], x = new double[n];
        for (int i = 0; i < n; ++i) {
            double value = b[f.p[i]];
            for (int j = 0; j < i; ++j) value -= f.l[i][j] * y[j];
            y[i] = value;
            if (!Double.isFinite(y[i])) throw new IllegalArgumentException("nonfinite forward substitution");
        }
        for (int i = n - 1; i >= 0; --i) {
            double value = y[i];
            for (int j = i + 1; j < n; ++j) value -= f.u[i][j] * x[j];
            x[i] = value / f.u[i][i];
            if (!Double.isFinite(x[i])) throw new IllegalArgumentException("nonfinite back substitution");
        }
        return x;
    }

    static double factorError(double[][] a, Factors f) {
        double error = 0.0;
        int n = a.length;
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j) {
                double value = 0.0;
                for (int k = 0; k < n; ++k) value += f.l[i][k] * f.u[k][j];
                error = Math.max(error, Math.abs(value - a[f.p[i]][j]));
            }
        return error;
    }

    static double solveError(double[][] a, double[] x, double[] b) {
        double error = 0.0;
        for (int i = 0; i < a.length; ++i) {
            double value = 0.0;
            for (int j = 0; j < x.length; ++j) value += a[i][j] * x[j];
            error = Math.max(error, Math.abs(value - b[i]));
        }
        return error;
    }

    static String vectorText(double[] v) {
        StringBuilder s = new StringBuilder("[");
        for (int i = 0; i < v.length; ++i) {
            if (i > 0) s.append(", ");
            s.append(String.format(Locale.ROOT, "%.6f", v[i] == 0.0 ? 0.0 : v[i]));
        }
        return s.append(']').toString();
    }

    static String matrixText(double[][] a) {
        StringBuilder s = new StringBuilder("[");
        for (int i = 0; i < a.length; ++i) {
            if (i > 0) s.append(", ");
            s.append(vectorText(a[i]));
        }
        return s.append(']').toString();
    }

    public static void main(String[] args) {
        double[][] a = { {0.0, 2.0, 1.0}, {1.0, 1.0, 0.0}, {2.0, 1.0, 1.0} };
        double[] b = {7.0, 3.0, 7.0}, b2 = {0.0, 1.0, 2.0};
        Factors f = luFactor(a, 1e-12);
        double[] x = luSolve(f, b), x2 = luSolve(f, b2);
        double fe = factorError(a, f), se = solveError(a, x, b), re = solveError(a, x2, b2);
        System.out.println("permutation = " + Arrays.toString(f.p));
        System.out.println("L = " + matrixText(f.l));
        System.out.println("U = " + matrixText(f.u));
        System.out.println("x = " + vectorText(x));
        System.out.println("reused x = " + vectorText(x2));
        System.out.printf(Locale.ROOT, "factor residual = %.3e%nsolve residual = %.3e%nreuse residual = %.3e%n",
                          fe, se, re);
        System.out.println("check = " + (Math.max(fe, Math.max(se, re)) < 1e-10 ? "PASS" : "FAIL"));
    }
}

All three programs print the same nine lines:

• permutation = [2, 0, 1]
• L = [[1.000000, 0.000000, 0.000000], [0.000000, 1.000000, 0.000000], [0.500000, 0.250000, 1.000000]]
• U = [[2.000000, 1.000000, 1.000000], [0.000000, 2.000000, 1.000000], [0.000000, 0.000000, -0.750000]]
• x = [1.000000, 2.000000, 3.000000]
• reused x = [1.000000, 0.000000, 0.000000]
• factor residual = 0.000e+00
• solve residual = 0.000e+00
• reuse residual = 0.000e+00
• check = PASS

Complexity, Correctness, and Invariants

For a successful dense $n\times n$ input, these implementations take $\Theta(n^3)$ factorization time in the best, average, and worst cases. They execute the trailing updates even when a multiplier is zero, so an identity matrix still follows the same cubic loop structure. The leading arithmetic count is about $2n^3/3$ floating-point operations, counting a multiplication and subtraction separately. Pivot searches and accumulated row/prefix swaps contribute $O(n^2)$ work.

The returned factors and working copy require $\Theta(n^2)$ space in total, with $O(n)$ permutation storage. Each solve takes $\Theta(n^2)$ time and $\Theta(n)$ additional space. Factor once and solve $q$ right-hand sides for total time $\Theta(n^3+qn^2)$; solving them sequentially needs only one solve workspace at a time. Refactoring for every right-hand side instead costs $\Theta(qn^3)$. An early pivot rejection can stop sooner, although input validation and copying already take quadratic work on valid square inputs.

For correctness, reason first in exact arithmetic with nonzero accepted pivots. At each stage, keep the evolving identity $PA=LU$, with identity entries in the unfilled part of $L$ and the current partially reduced matrix in $U$. Initially $P=I$, $L=I$, and $U=A$, so the identity holds.

A swap matrix $S$ exchanging active rows $k$ and $r$ changes $P$ to $SP$ and $U$ to $SU$. Swapping the earlier entries of those rows in $L$ is equivalent, in its current partially filled shape, to replacing $L$ by $SLS$. Therefore $(SLS)(SU)=SLU=SPA$. This is why the stored multipliers must move and why swapping whole rows of $L$ would be wrong.

Next, subtracting $m_{ik}$ times pivot row $k$ from row $i$ of $U$ is multiplication by an elementary elimination matrix $E$. Its inverse adds that multiple back. Recording $m_{ik}$ in the new column of $L$ effects $L\leftarrow LE^{-1}$ while $U\leftarrow EU$, preserving the product. Already completed columns of $U$ remain zero below their diagonal because both active rows have zeros there. At termination, $U$ is upper triangular and $L$ is unit lower triangular.

Finally, $Ly=Pb$ and $Ux=y$ imply $LUx=Pb$, hence $PAx=Pb$, and invertibility of $P$ gives $Ax=b$. Nonzero pivots make the triangular solution unique. Floating-point arithmetic yields approximate identities instead of exact ones; partial pivoting limits multipliers to magnitude at most $1$ but does not remove all possible element growth or sensitivity to ill conditioning.

Variants, Pitfalls, and Optimizations

Permutation direction is part of the interface. With $PA=LU$, solve using $Pb$, meaning b[p[i]]. For the final ordering $(2,0,1)$, this gives $(b_2,b_0,b_1)$; applying the inverse order would generally give a different vector. A library's sequential pivot log must be applied as swaps, not treated as this final indexing array.

A later pivot must move earlier multipliers. In the worked example, omitting the prefix swap leaves the wrong matrix

$$ L_{\mathrm{bad}}=\begin{pmatrix}1&0&0\\1/2&1&0\\0&1/4&1\end{pmatrix}. $$

Using it with the otherwise correct final $U$ and $p$ gives a factor residual of $1$. For the original $b$, forward substitution gives $y_{\mathrm{bad}}=(7,7/2,17/8)^T$, and the resulting solution is $x_{\mathrm{bad}}=(10/3,19/6,-17/6)^T$. Its original-system residual is $\lVert Ax_{\mathrm{bad}}-b\rVert_\infty=7/2$. A triangular-looking result alone cannot establish correctness.

Choice of factorization and storage. Doolittle fixes the diagonal of $L$ to ones; Crout instead fixes the diagonal of $U$. Packed LU stores strict-lower multipliers and upper-triangular entries in a single matrix, with the unit diagonal of $L$ implicit. In packed storage a whole row swap moves both types of stored data, while our separate arrays require a whole-row swap in $U$ and a prefix swap in $L$.

Dense, sparse, and structured cases. Blocking groups trailing updates into matrix products for better cache use and optimized BLAS execution. Sparse LU adds ordering and symbolic analysis to control fill-in, since zeros can become nonzeros. Cholesky is attractive for symmetric positive definite systems; QR or SVD better fits many least-squares and rank-sensitive tasks. Unpivoted LU needs additional guarantees; a zero leading entry does not by itself make a matrix singular.

Numerical checks and reuse. A small pivot under this threshold may reflect poor scaling rather than exact singularity: $\operatorname{diag}(1,10^{-14})$ is nonsingular but rejected with $\tau=10^{-12}$. Conversely, passing the test does not certify good conditioning. A useful scale-aware solve diagnostic is

$$ \eta=\frac{\lVert Ax-b\rVert_\infty} {\lVert A\rVert_\infty\lVert x\rVert_\infty+\lVert b\rVert_\infty}, \qquad \lVert A\rVert_\infty=\max_i\sum_j\lvert A_{ij}\rvert. $$

If the denominator is zero, treat an exactly zero residual separately. A small backward residual can coexist with a large solution error for an ill-conditioned matrix. Iterative refinement reuses the factors to solve a correction equation from a newly evaluated residual. Reuse applies while $A$ stays fixed; changing its entries generally requires new factors. To solve $Ax=b$, apply the factors directly instead of first forming an explicit inverse.

Practice Problems and Real-World Use

Start with hand traces where every multiplier is a simple fraction. Then test a factorization independently from its solve routine: verify triangular shape, the permutation, $LU\approx PA$, and $Ax\approx b$ for known solutions. Include a negative one-by-one matrix, a swap at a later pivot, a nonsingular matrix rejected at the default tolerance, and a singular matrix whose final pivot is zero.

For a reusable solver, choose several known vectors $x^{(t)}$, compute $b^{(t)}=Ax^{(t)}$, factor $A$ once, and recover every vector from the same factors. Check residuals with an appropriate scale normalization as the matrix size and magnitude change. This separates factorization bugs from permutation or substitution bugs.