Back to Data Structures

Gaussian Elimination

October 16, 2026 Wasil Zafar 29 min read

Concepts, invariants, and Python, C++, and Java implementations of Gaussian Elimination

Table of Contents

  1. Why Gaussian Elimination Matters
  2. Core Concepts and Terminology
  3. How Gaussian Elimination Works, Step by Step
  4. Worked Example with a Full Trace
  5. Implementation in Python, C++, and Java
  6. Complexity, Correctness, and Invariants
  7. Variants, Pitfalls, and Optimizations
  8. Practice Problems and Real-World Use

Why Gaussian Elimination Matters

Three linked measurements can produce three equations in the same three unknowns. Solving one equation at a time is awkward because changing an unknown affects the others. Gaussian elimination systematically removes variables until the last equation contains only one unknown, then reconstructs the rest by substitution.

The method is a foundation of numerical linear algebra. Circuit models, structural calculations, interpolation systems, and many other applications eventually require solving $Ax=b$. Recording elimination's multipliers and row exchanges also leads to a reusable LU factorization.

This lesson develops a square-system solver with partial pivoting. The worked example starts with a zero diagonal and requires another row swap in the next column, so pivot selection remains visible throughout. You will implement the same algorithm in Python, C++, and Java, preserve the original equations, and inspect both an absolute residual and a scaled error diagnostic. The discussion distinguishes exact algebraic correctness from floating-point accuracy.

Gaussian elimination solves a linear system by applying reversible row operations to produce an equivalent upper-triangular system, then using back substitution. Partial pivoting chooses a usable pivot at every column and helps control the size of elimination multipliers.

Core Concepts and Terminology

Let $A$ be an $n\times n$ matrix, $b$ a vector of length $n$, and $x$ the unknown vector. The equations are $\sum_j A_{ij}x_j=b_i$. Store a working augmented matrix $[A\mid b]$, so swapping a whole row also moves its right-hand side.

Allowed elementary row operations exchange two equations, multiply an equation by a nonzero scalar, or add a multiple of another equation. These preserve the solution set in exact arithmetic. Our implementation uses exchanges and row subtraction during forward elimination; it does not normalize each pivot to one.

At pivot column $k$, partial pivoting selects a row $p\ge k$ that maximizes $\lvert U_{pk}\rvert$, where $U$ is the current coefficient workspace. After exchanging rows $p$ and $k$, eliminate a lower entry with multiplier

$$ f_{ik}=\frac{U_{ik}}{U_{kk}},\qquad R_i\longleftarrow R_i-f_{ik}R_k. $$

Here $R_i$ includes the right-hand side. The sign of a pivot does not matter; its magnitude must be usable. In exact arithmetic, largest-magnitude column selection makes $\lvert f_{ik}\rvert\le1$.

Use zero-based indices. The worked system is

$$ \begin{aligned} 2x_1+x_2&=7,\\ x_0+x_1&=3,\\ 2x_0+x_2&=5, \end{aligned}\qquad A=\begin{pmatrix}0&2&1\\1&1&0\\2&0&1\end{pmatrix},\quad b=\begin{pmatrix}7\\3\\5\end{pmatrix}. $$

A unique solution exists exactly when the square coefficient matrix is nonsingular. A zero pivot before row selection does not establish singularity: the current column may contain a suitable entry lower down. If no usable pivot remains, the numerical routine rejects the solve; further analysis is needed to classify rank or consistency.

From original equations to upper-triangular equationsThe original augmented matrix 0,2,1,7;1,1,0,3;2,0,1,5 is converted to 2,0,1,5;0,2,1,7;0,0,-1,-3. The pivots are 2,2,-1.Eliminate variables while preserving the equationsThe last column is the right-hand side and moves with every row.Original augmented system021711032015Whole-row swaps + row subtraction → equivalent systemUpper-triangular system2015021700-1-3Zeros below the diagonal allow back substitution.Last row gives x2; then solve x1; finally solve x0.
TermDefinitionWorked example or role
Augmented matrixCoefficients with right-hand side appended$[A\mid b]$ has three rows and four columns
Active columnPivot column restricted to rows not yet finalizedAt $k=1$, compare only rows $1$ and $2$
Partial pivotingLargest absolute active-column entry determines the row swapFirst select magnitude $2$ in row $2$
MultiplierFactor used to cancel a lower pivot-column entry$1/2$ in both nonzero elimination updates
Upper-triangular formEntries below the diagonal are zeroFinal diagonal entries are $2,2,-1$
Back substitutionSolve unknowns from the last row upwardRecover $x_2$, then $x_1$, then $x_0$
Fill-inAn initial zero becomes nonzero during eliminationA working coefficient becomes $-1/2$
Residual $r$Error in the original equations$r=Ax-b$
Backward errorRelative data perturbation explaining the computed solutionA scaled residual, distinct from solution error

How Gaussian Elimination Works, Step by Step

The code accepts a nonempty square system with finite double-precision inputs and a finite relative pivot tolerance $0\lt\tau\lt1$. Its default is $\tau=10^{-12}$. It copies the system and computes the original coefficient scale

$$ \alpha=\max_{i,j}\lvert A_{ij}\rvert. $$

The pivot threshold is $\tau\alpha$. This is a simple scale-relative heuristic. Uniformly rescaling a moderate system rescales the threshold, but disparate row or column scales can still cause an invertible matrix to be rejected.

  1. Build a fresh augmented workspace with the coefficients followed by the right-hand side.
  2. For each column $k=0,\ldots,n-1$, search active rows for the largest absolute pivot candidate. Break ties by keeping the first maximum.
  3. Reject if the selected magnitude is at most $\tau\alpha$. The failure message covers singular or numerically rejected systems without claiming an exact rank classification.
  4. Swap the entire selected row with row $k$.
  5. For each row $i\gt k$, compute the multiplier from the old entry, set that pivot-column entry to zero, and update every later column, including the appended right-hand side.
  6. After the forward phase, solve upward from the last row:

$$ x_i=\frac{g_i-\sum_{j=i+1}^{n-1}U_{ij}x_j}{U_{ii}}, $$

where $g$ is the transformed right-hand side.

  1. Use the original $A$ and $b$ to calculate diagnostics. Reusing overwritten coefficients would check a different representation and could conceal mistakes.

Partial pivoting is the principle used in general LU-based linear-system solvers such as LAPACK's DGESV. The programs here expose the operations in plain nested arrays rather than calling a numerical package.

One full Gaussian elimination stepSelect the magnitude-two first-column pivot, swap full augmented rows zero and two, then eliminate the magnitude-one entry using factor one half. Row two has zero multiplier.First pivot column: select, swap, eliminateCompare absolute values 0, 1, 2; choose row 2.Before column 0021711032015After row 0 ↔ row 2201511030217After lower-row updates201501-1/21/20217Swap the whole row, including the right-hand side.Row 1 multiplier = 1/2; row 2 multiplier = 0.Row 1: [1, 1, 0 | 3] − (1/2)[2, 0, 1 | 5] = [0, 1, −1/2 | 1/2].
Gaussian Elimination control flow
flowchart TD
    A["Read A, b, and tolerance"] --> B{"Valid finite square system?"}
    B -- No --> C["Reject input"]
    B -- Yes --> D["Copy augmented system; compute scale; k = 0"]
    D --> E{"k less than n?"}
    E -- Yes --> F["Choose largest absolute active-column pivot"]
    F --> G{"Pivot magnitude exceeds threshold?"}
    G -- No --> H["Report singular or numerically rejected solve"]
    G -- Yes --> I["Swap whole augmented rows"]
    I --> J["Eliminate all rows below the pivot"]
    J --> K["Increment k"]
    K --> E
    E -- No --> L["Back substitute from last row upward"]
    L --> M["Evaluate original-system residual and backward error"]
    M --> N["Return solution and diagnostics"]

Worked Example with a Full Trace

The initial workspace is

$$ \left[\begin{array}{ccc|c}0&2&1&7\\1&1&0&3\\2&0&1&5\end{array}\right]. $$

Here $\alpha=2$, so the threshold is $2\times10^{-12}$.

Column $k=0$. The active pivot magnitudes are $0,1,2$. Swap rows $0$ and $2$:

$$ \left[\begin{array}{ccc|c}2&0&1&5\\1&1&0&3\\0&2&1&7\end{array}\right]. $$

For row $1$, the multiplier is $1/2$. After explicitly writing the leading zero, update the remaining entries: $1-(1/2)0=1$, $0-(1/2)1=-1/2$, and $3-(1/2)5=1/2$. The workspace becomes

$$ \left[\begin{array}{ccc|c}2&0&1&5\\0&1&-1/2&1/2\\0&2&1&7\end{array}\right]. $$

Row $2$ has multiplier $0$, so its later entries remain $2,1,7$. This completes the first column.

Column $k=1$. Compare the active magnitudes $1$ and $2$. A second pivot swap is necessary: exchange rows $1$ and $2$.

$$ \left[\begin{array}{ccc|c}2&0&1&5\\0&2&1&7\\0&1&-1/2&1/2\end{array}\right]. $$

The remaining lower multiplier is $1/2$. Write its column-one entry as zero, then update $-1/2-(1/2)1=-1$ and $1/2-(1/2)7=-3$:

$$ [U\mid g]=\left[\begin{array}{ccc|c}2&0&1&5\\0&2&1&7\\0&0&-1&-3\end{array}\right]. $$

Column $k=2$. The final pivot is $-1$. Its magnitude is $1$, so it is accepted; no lower row remains. The two swaps and all lower-row updates are now complete.

Back substitution. The last row gives $x_2=(-3)/(-1)=3$. The middle row gives $x_1=(7-1\cdot3)/2=2$. For the first row, the running numerator starts at $5$, stays $5$ after subtracting $0\cdot2$, then becomes $2$ after subtracting $1\cdot3$. Thus $x_0=2/2=1$.

Original-system check. The row products are $0\cdot1+2\cdot2+1\cdot3=7$, $1\cdot1+1\cdot2+0\cdot3=3$, and $2\cdot1+0\cdot2+1\cdot3=5$. The residual is zero in this exactly representable sample. The final pivot product is $-4$; two row exchanges preserve its sign, giving $\det(A)=-4$ as another algebraic check.

Second pivot and complete back substitutionSwap the remaining active rows to put two in the second pivot. A half-multiple subtraction leaves final row 0,0,-1,-3. Back substitution gives 3,2,1 in reverse variable order.A second pivot swap is requiredAt column 1, active entries have magnitudes 1 and 2.Swap rows 1 and 22015021701-1/21/2Eliminate with factor 1/22015021700-1-3Last pivot is −1: magnitude 1 is usable.Last rowx2 = (−3) / (−1) = 3Middle rowx1 = (7 − 3) / 2 = 2First rowx0 = (5 − 0 × 2 − 3) / 2 = 1Original-system check: A × [1, 2, 3] = [7, 3, 5] = b.
Worked Example

Trace It by Hand

After the first column, the lower augmented rows are $[0,1,-1/2\mid1/2]$ and $[0,2,1\mid7]$. Predict the next pivot row, multiplier, and final lower row before reading on.

Select the row with column-one magnitude $2$ and swap the entire augmented rows. The lower multiplier is $1/2$, and its new row is $[0,0,-1\mid-3]$. The negative pivot is valid because its magnitude exceeds the threshold. Finish by recovering $x_2=3$, $x_1=2$, and $x_0=1$.

Implementation in Python, C++, and Java

Every panel implements the same pivot searches, full augmented-row swaps, elimination order, and back substitution. All inputs remain unchanged. The solver rejects malformed shapes, nonfinite solve inputs, invalid tolerances, and pivots below its numerical threshold. It is a unique-solution square-system solver; it does not return a rank, a nullspace basis, or a least-squares fit.

The diagnostics routine assumes compatible original-system dimensions. It computes $r=Ax-b$ and the normwise backward-error diagnostic

$$ \eta_\infty=\frac{\lVert r\rVert_\infty}{\lVert A\rVert_\infty\lVert x\rVert_\infty+\lVert b\rVert_\infty}, $$

where the matrix infinity norm is the maximum absolute row sum. If the denominator and residual are both zero, the diagnostic is defined as zero. This interpretation is discussed in LAPACK's error analysis for linear systems.

All programs print the same four lines: x = [1.000000, 2.000000, 3.000000]; residual_inf = 0.000e+00; backward_error = 0.000e+00; check = PASS. The sample check requires an absolute residual at most $10^{-10}$ and backward error at most $10^{-12}$. Save the Java panel as Main.java.

The routines target moderate, sensibly scaled values. Finite input validation does not prevent intermediate overflow on extreme magnitudes. For general problems, scale the diagnostics appropriately and consider a condition estimate; a small backward error alone does not guarantee small error in the unknown vector.

from math import isfinite


def solve(a: list[list[float]], b: list[float], tol: float = 1e-12) -> list[float]:
    n = len(a)
    if n == 0 or len(b) != n or any(len(row) != n for row in a):
        raise ValueError("require a nonempty square system")
    if not isfinite(tol) or not 0 < tol < 1:
        raise ValueError("require 0 < tol < 1")
    if any(not isfinite(v) for row in a for v in row) or any(not isfinite(v) for v in b):
        raise ValueError("inputs must be finite")
    aug = [a[i][:] + [b[i]] for i in range(n)]
    scale = max(abs(v) for row in a for v in row)
    for k in range(n):
        pivot = max(range(k, n), key=lambda i: abs(aug[i][k]))
        if abs(aug[pivot][k]) <= tol * scale:
            raise ValueError("singular or numerically rejected system")
        aug[k], aug[pivot] = aug[pivot], aug[k]
        for i in range(k + 1, n):
            factor = aug[i][k] / aug[k][k]
            aug[i][k] = 0.0
            for j in range(k + 1, n + 1):
                aug[i][j] -= factor * aug[k][j]
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        value = aug[i][n]
        for j in range(i + 1, n):
            value -= aug[i][j] * x[j]
        x[i] = value / aug[i][i]
    return x


def diagnostics(a: list[list[float]], x: list[float],
                b: list[float]) -> tuple[float, float]:
    absolute = norm_a = norm_b = 0.0
    for i, row in enumerate(a):
        value = 0.0
        row_sum = 0.0
        for j, entry in enumerate(row):
            value += entry * x[j]
            row_sum += abs(entry)
        absolute = max(absolute, abs(value - b[i]))
        norm_a = max(norm_a, row_sum)
        norm_b = max(norm_b, abs(b[i]))
    norm_x = max(abs(v) for v in x)
    denominator = norm_a * norm_x + norm_b
    return absolute, absolute / denominator if denominator != 0 else 0.0


def main() -> None:
    a = [[0.0, 2.0, 1.0], [1.0, 1.0, 0.0], [2.0, 0.0, 1.0]]
    b = [7.0, 3.0, 5.0]
    x = solve(a, b)
    absolute, backward = diagnostics(a, x, b)
    print("x = [" + ", ".join(f"{v:.6f}" for v in x) + "]")
    print(f"residual_inf = {absolute:.3e}")
    print(f"backward_error = {backward:.3e}")
    print("check = " + ("PASS" if absolute <= 1e-10 and backward <= 1e-12 else "FAIL"))


if __name__ == "__main__":
    main()
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <utility>
#include <vector>

using Matrix = std::vector<std::vector<double>>;
using Vector = std::vector<double>;

Vector solve(const Matrix& a, const Vector& b, double tol = 1e-12) {
    const std::size_t n = a.size();
    if (n == 0 || b.size() != n) throw std::invalid_argument("invalid system shape");
    if (!std::isfinite(tol) || tol <= 0 || tol >= 1) throw std::invalid_argument("invalid tolerance");
    Matrix aug(n, Vector(n + 1));
    double scale = 0;
    for (std::size_t i = 0; i < n; ++i) {
        if (a[i].size() != n) throw std::invalid_argument("matrix must be square");
        for (std::size_t j = 0; j < n; ++j) {
            if (!std::isfinite(a[i][j])) throw std::invalid_argument("inputs must be finite");
            aug[i][j] = a[i][j];
            scale = std::max(scale, std::abs(a[i][j]));
        }
        if (!std::isfinite(b[i])) throw std::invalid_argument("inputs must be finite");
        aug[i][n] = b[i];
    }
    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(aug[i][k]) > std::abs(aug[pivot][k])) pivot = i;
        if (std::abs(aug[pivot][k]) <= tol * scale) throw std::runtime_error("singular or numerically rejected system");
        std::swap(aug[k], aug[pivot]);
        for (std::size_t i = k + 1; i < n; ++i) {
            const double factor = aug[i][k] / aug[k][k];
            aug[i][k] = 0;
            for (std::size_t j = k + 1; j <= n; ++j) aug[i][j] -= factor * aug[k][j];
        }
    }
    Vector x(n, 0);
    for (std::size_t index = n; index > 0; --index) {
        const std::size_t i = index - 1;
        double value = aug[i][n];
        for (std::size_t j = i + 1; j < n; ++j) value -= aug[i][j] * x[j];
        x[i] = value / aug[i][i];
    }
    return x;
}

std::pair<double, double> diagnostics(const Matrix& a, const Vector& x, const Vector& b) {
    double absolute = 0, norm_a = 0, norm_b = 0, norm_x = 0;
    for (std::size_t i = 0; i < a.size(); ++i) {
        double value = 0, row_sum = 0;
        for (std::size_t j = 0; j < a[i].size(); ++j) {
            value += a[i][j] * x[j];
            row_sum += std::abs(a[i][j]);
        }
        absolute = std::max(absolute, std::abs(value - b[i]));
        norm_a = std::max(norm_a, row_sum);
        norm_b = std::max(norm_b, std::abs(b[i]));
    }
    for (double v : x) norm_x = std::max(norm_x, std::abs(v));
    const double denominator = norm_a * norm_x + norm_b;
    return {absolute, denominator != 0 ? absolute / denominator : 0};
}

int main() {
    const Matrix a = { {0, 2, 1}, {1, 1, 0}, {2, 0, 1}};
    const Vector b = {7, 3, 5};
    const Vector x = solve(a, b);
    const auto errors = diagnostics(a, x, b);
    std::cout << "x = [" << std::fixed << std::setprecision(6);
    for (std::size_t i = 0; i < x.size(); ++i) {
        if (i != 0) std::cout << ", ";
        std::cout << x[i];
    }
    std::cout << "]\n" << std::scientific << std::setprecision(3)
              << "residual_inf = " << errors.first << '\n'
              << "backward_error = " << errors.second << '\n';
    std::cout << "check = " << (errors.first <= 1e-10 && errors.second <= 1e-12 ? "PASS" : "FAIL") << '\n';
    return 0;
}
import java.util.Locale;

public class Main {
    public static double[] solve(double[][] a, double[] b, double tol) {
        int n = a.length;
        if (n == 0 || b.length != n) throw new IllegalArgumentException("invalid system shape");
        if (!Double.isFinite(tol) || tol <= 0 || tol >= 1) throw new IllegalArgumentException("invalid tolerance");
        double[][] aug = new double[n][n + 1];
        double scale = 0;
        for (int i = 0; i < n; ++i) {
            if (a[i].length != n) throw new IllegalArgumentException("matrix must be square");
            for (int j = 0; j < n; ++j) {
                if (!Double.isFinite(a[i][j])) throw new IllegalArgumentException("inputs must be finite");
                aug[i][j] = a[i][j];
                scale = Math.max(scale, Math.abs(a[i][j]));
            }
            if (!Double.isFinite(b[i])) throw new IllegalArgumentException("inputs must be finite");
            aug[i][n] = b[i];
        }
        for (int k = 0; k < n; ++k) {
            int pivot = k;
            for (int i = k + 1; i < n; ++i)
                if (Math.abs(aug[i][k]) > Math.abs(aug[pivot][k])) pivot = i;
            if (Math.abs(aug[pivot][k]) <= tol * scale) throw new IllegalArgumentException("singular or numerically rejected system");
            double[] row = aug[k]; aug[k] = aug[pivot]; aug[pivot] = row;
            for (int i = k + 1; i < n; ++i) {
                double factor = aug[i][k] / aug[k][k];
                aug[i][k] = 0;
                for (int j = k + 1; j <= n; ++j) aug[i][j] -= factor * aug[k][j];
            }
        }
        double[] x = new double[n];
        for (int i = n - 1; i >= 0; --i) {
            double value = aug[i][n];
            for (int j = i + 1; j < n; ++j) value -= aug[i][j] * x[j];
            x[i] = value / aug[i][i];
        }
        return x;
    }

    public static double[] diagnostics(double[][] a, double[] x, double[] b) {
        double absolute = 0, normA = 0, normB = 0, normX = 0;
        for (int i = 0; i < a.length; ++i) {
            double value = 0, rowSum = 0;
            for (int j = 0; j < a[i].length; ++j) {
                value += a[i][j] * x[j];
                rowSum += Math.abs(a[i][j]);
            }
            absolute = Math.max(absolute, Math.abs(value - b[i]));
            normA = Math.max(normA, rowSum);
            normB = Math.max(normB, Math.abs(b[i]));
        }
        for (double v : x) normX = Math.max(normX, Math.abs(v));
        double denominator = normA * normX + normB;
        return new double[]{absolute, denominator != 0 ? absolute / denominator : 0};
    }

    public static void main(String[] args) {
        double[][] a = { {0, 2, 1}, {1, 1, 0}, {2, 0, 1}};
        double[] b = {7, 3, 5};
        double[] x = solve(a, b, 1e-12);
        double[] errors = diagnostics(a, x, b);
        System.out.print("x = [");
        for (int i = 0; i < x.length; ++i) {
            if (i != 0) System.out.print(", ");
            System.out.printf(Locale.ROOT, "%.6f", x[i]);
        }
        System.out.println("]");
        System.out.printf(Locale.ROOT, "residual_inf = %.3e%nbackward_error = %.3e%n", errors[0], errors[1]);
        System.out.println("check = " + (errors[0] <= 1e-10 && errors[1] <= 1e-12 ? "PASS" : "FAIL"));
    }
}

Complexity, Correctness, and Invariants

For an accepted $n\times n$ system, pivot searches cost $O(n^2)$ overall. Elimination visits $n(n-1)/2$ lower-row targets, with successively shorter trailing updates. The resulting work is $\Theta(n^3)$; back substitution and dense residual evaluation each add $\Theta(n^2)$. Copying the augmented matrix requires $\Theta(n^2)$ storage, plus $\Theta(n)$ for the solution.

The shown loops do not skip zero multipliers. Consequently, among successful solves of a fixed dimension, best, average, and worst asymptotic times are all $\Theta(n^3)$, including diagonal inputs. A rejection at the first pivot can finish after $\Theta(n^2)$ validation and construction. An average over successful and rejected inputs needs a distribution; it is not assumed here.

In exact arithmetic, a row swap simply reorders equations, and subtracting a multiple of a pivot equation has a reversible inverse operation. The transformed system therefore has the same solutions. After a pivot column is eliminated, earlier zero columns stay zero because all active rows already have zeros there. Repeating this establishes an equivalent upper-triangular system. Nonzero pivots then make each back-substitution division valid and determine a unique solution.

Floating-point error must be assessed separately. Partial pivoting controls the current multipliers but does not guarantee small element growth for every matrix. A small backward error means that small perturbations of the original data can explain the computed solution; the actual solution error also depends on conditioning. The pivot threshold is a heuristic acceptance rule, not a proof of singularity or an accuracy bound.

MethodTimeWorkspace or storageKey assumptions and purpose
Dense elimination with partial pivoting$\Theta(n^3)$ on successful solves$\Theta(n^2)$ hereSquare system accepted by the tolerance
Back substitution after triangularization$\Theta(n^2)$$\Theta(n)$ outputUpper-triangular matrix with usable diagonal entries
Reuse retained pivoted LU factors$\Theta(n^2)$ per additional right-hand side$\Theta(n^2)$ factorsFactorization and permutation already available
Gauss–Jordan elimination$\Theta(n^3)$$\Theta(n^2)$Eliminates above pivots too; useful for reduced form or inverse demonstrations
Cholesky factorization and solve$\Theta(n^3)$ factorization$\Theta(n^2)$Symmetric positive-definite matrix; exploits structure
Sparse direct eliminationPattern-dependentDepends on factor fill-inOriginal nonzero count alone does not determine cost

After column $k$ is processed, the augmented equations remain equivalent to the original system in exact arithmetic, and every entry below the first $k+1$ pivot positions is zero. Whole-row swaps and elimination preserve this invariant; accepted nonzero pivots then justify solving the triangular equations from bottom to top.

Variants, Pitfalls, and Optimizations

Choose pivots repeatedly and by magnitude. The worked example needs a swap in both columns $0$ and $1$. Looking only for an initial row permutation is insufficient. A negative pivot is acceptable; compare its absolute value with the threshold. Dividing by a tiny pivot without pivoting can create huge multipliers and destructive cancellation, as the next diagram illustrates.

Preserve the right-hand side and fill-in. A row exchange moves the complete augmented row, not just the coefficients. Every subtraction updates the right-hand side as well. Newly created nonzero entries are part of the transformed equations. Forcing the original sparse zero pattern to remain unchanged alters the system unless an approximation is intentional.

Distinguish failure categories. In exact arithmetic, an echelon row $[0,\ldots,0\mid c]$ with $c\ne0$ establishes inconsistency. A redundant all-zero row can occur in a consistent system with free variables. The supplied solver rejects when its next pivot is unusable and does not continue with column skipping to distinguish these cases. A floating-point threshold rejection can also occur for an invertible matrix.

Factor once for repeated solves. Store the multipliers, row permutations, and upper-triangular factor to form pivoted LU. Then each additional right-hand side requires triangular solves instead of repeating the cubic factorization. Avoid forming $A^{-1}$ merely to solve $Ax=b$. The current routine returns only the solution, so repeated calls repeat elimination.

Use more suitable structure when available. Symmetric positive-definite matrices support Cholesky. Large sparse matrices need sparse ordering and factor storage, or an iterative method with appropriate preconditioning. Sparse elimination can introduce substantial fill-in. Rectangular least-squares problems and numerical rank questions are usually better handled with QR or SVD than by extending this unique-square-solve interface blindly.

Choose the arithmetic deliberately. Exact rational elimination gives exact answers but can grow large numerators and denominators; fraction-free methods can help. Over a finite field, test pivots for nonzero field values and use modular inverses rather than floating tolerances. In floating point, row or column scaling may help, but partial pivoting and a fixed tolerance are not substitutes for conditioning analysis. Iterative refinement can reuse factors while recomputing residuals from the original matrix.

Partial pivoting avoids a tiny leading pivotFor rows 1e-20,1 and 1,1 with right-hand side 1,2, pivoting chooses the magnitude-one entry. Unsafe elimination with the tiny pivot creates multiplier 1e20 and may produce solution 0,1 with residual norm one.A tiny pivot can destroy an unpivoted calculationSystem: 10⁻²⁰ x0 + x1 = 1; x0 + x1 = 2.Correct: swap to use the magnitude-one pivotMultiplier ≈ 10⁻²⁰; computed solution ≈ [1, 1].The original-system residual is tiny.Row selection avoids dividing by the tiny first entry.Unsafe: keep the tiny pivot and ignore toleranceMultiplier = 10²⁰; subtraction loses small terms.Typical double result: [0, 1], with residual norm 1.A nonzero pivot can still be a poor numerical choice.Partial pivoting controls multipliers; conditioning still matters.

Search the active column at every step, move the whole augmented row, and test pivot magnitudes rather than signed values. Retain fill-in and evaluate errors against the original system; do not interpret a tolerance rejection as an exact rank classification or a small residual as a universal accuracy guarantee.

Practice Problems and Real-World Use

Start with the equations, their units, and the expected matrix structure. Verify small systems by direct substitution, then test the numerical diagnostics independently. These exercises cover both the algorithm's control flow and the limits of its square-system interface.

Exercise or applicationMain idea or checkDifficulty
Work the three-equation exampleTwo row swaps give pivots $2,2,-1$ and solution $[1,2,3]^{\mathsf T}$Easy
Triangular-system exerciseVerify that back substitution alone solves a supplied upper-triangular systemEasy
Compare two singular systemsUse proportional coefficient rows with consistent and inconsistent right-hand sides; explain why the routine rejects bothMedium
Tiny leading pivotTest $\begin{pmatrix}10^{-20}&1\\1&1\end{pmatrix}$ with $b=(1,2)^{\mathsf T}$; compare pivoted and unsafe unpivoted arithmeticMedium
Repeated circuit-model solvesRetain factors for several right-hand sides and check residuals for every solveMedium
Polynomial interpolation at distinct nodesForm a small Vandermonde system and recover coefficients; inspect conditioning as nodes become closeMedium
Sparse finite-difference equationsTrack fill-in and compare dense reference answers with an appropriate structured or sparse solverHard

Construct independent tests by selecting a known vector, multiplying by a well-conditioned matrix to obtain $b$, and checking recovery. Include row permutations, negative pivots, zero right-hand sides, malformed shapes, nonfinite inputs, and scale changes. For tiny dimensions, exact rational or determinant-based solutions provide an oracle independent of the elimination implementation.