Gauss-Jordan

Why Gauss-Jordan Matters

Solving a linear system becomes especially transparent when row operations turn its coefficient matrix into the identity. The remaining right-hand-side entries then display the unknowns directly. Gauss–Jordan elimination pursues exactly that goal: normalize each pivot and remove the other entries in its column, both above and below it.

The same row operations can act on several right-hand sides at once. If an identity matrix is included among them, its transformed columns become the inverse of the original coefficient matrix. This makes Gauss–Jordan useful for understanding reduced row-echelon form, solving small systems, and demonstrating matrix inverses.

This lesson uses partial pivoting, a scale-relative floating-point tolerance, and plain nested arrays. The worked example carries one system vector and an identity block together, producing a solution and an inverse in one pass. The implemented interface targets nonsingular square coefficient matrices; the variants explain how general row reduction extends to rank-deficient or rectangular problems.

Gauss–Jordan elimination uses reversible row operations to make every accepted pivot equal to one and every other entry in its column equal to zero. For a nonsingular square matrix, reducing $[A\mid B]$ gives $[I\mid X]$ with $AX=B$; choosing $B=I$ gives the inverse.

Core Concepts and Terminology

A matrix is in reduced row-echelon form, or RREF, when zero rows are at the bottom, each nonzero row's leading entry is one, these leading ones move rightward down the rows, and each leading one is the only nonzero entry in its column. Gaussian elimination commonly stops at upper-triangular or echelon form. Gauss–Jordan continues until the pivot columns are fully reduced.

For the square nonsingular case, use a coefficient matrix $A\in\mathbb R^{n\times n}$ and a right-hand-side block $B\in\mathbb R^{n\times q}$, with $q\ge1$. Each column of $B$ represents a separate system with the same coefficients. The algorithm transforms the entire augmented matrix of width $n+q$:

$$ [A\mid B]\longrightarrow[I_n\mid X],\qquad AX=B. $$

Row exchanges, multiplication of a row by a nonzero scalar, and adding a multiple of another row preserve all represented systems. In particular, reducing $[A\mid I_n]$ transforms the identity block into $A^{-1}$. Auburn University's elementary-matrix notes explain this connection between row operations and inverse matrices.

The worked data are

$$ A=\begin{pmatrix}0&2\\1&3\end{pmatrix},\qquad b=\begin{pmatrix}4\\7\end{pmatrix},\qquad B=[b\mid I_2]=\begin{pmatrix}4&1&0\\7&0&1\end{pmatrix}. $$

The augmented matrix therefore has two coefficient columns followed by three right-hand-side columns. Only coefficient columns are pivot candidates in this implementation. The first result column solves $Ax=b$; the next two result columns form $A^{-1}$.

Term	Meaning	Worked example or role
Pivot column	Coefficient column being reduced	Columns $0$ and $1$, never the appended data columns
Partial pivoting	Largest absolute candidate in the active column selects the row	Column $0$ chooses row $1$ because its entry is $1$
Pivot normalization	Divide the entire row by its saved pivot	Divide the second pivot row by $2$
Elimination factor	Entry removed from another row after normalization	Factor $3$ removes the entry above the second pivot
Reduced pivot column	Unit vector with one leading one	Final coefficient columns are the two columns of $I_2$
Multiple right-hand sides	Several systems carried through the same operations	$B$ contains $b$ and both identity columns
Inverse block	Transformed identity columns	$A^{-1}=\begin{pmatrix}-3/2&1\\1/2&0\end{pmatrix}$
Residual	Error measured against the original data	Check $Ax-b$ and $AA^{-1}-I$
Relative tolerance $\tau$	Numerical pivot acceptance rule	Default $\tau=10^{-12}$ relative to the original coefficient scale

How Gauss-Jordan Works, Step by Step

The function accepts a nonempty square $A$, a matching finite right-hand-side matrix $B$ with at least one column, and a finite tolerance $0\lt\tau\lt1$. It returns a fresh full augmented matrix, approximately $[I\mid X]$, while preserving both input matrices.

Copy $A$ and $B$ side by side. Compute coefficient scale $\alpha=\max_{ij}\lvert A_{ij}\rvert$ and threshold $\tau\alpha$.
For pivot column $k$, search rows $k,\ldots,n-1$ for the largest absolute coefficient. Use the first maximum to make ties deterministic.
Reject if its magnitude is at most the threshold. Otherwise swap the complete augmented rows.
Save the pivot value before changing any entry. Divide every entry in that row, including all right-hand sides, by the saved value. The pivot becomes one.
For every other row $i\ne k$, save its current coefficient $f$ in column $k$. Subtract $f$ times the normalized pivot row from the entire row, then explicitly set its pivot-column entry to zero.
Continue to the next coefficient column. At completion, extract the right-hand-side blocks and verify them using the original matrices.

The coefficient scale excludes $B$: the pivot test concerns the coefficient system. This simple global threshold may reject an invertible but poorly scaled matrix. It is a numerical policy, not an exact singularity or rank certificate.

The two saved scalars matter. If the pivot is reread after its own entry has been divided, it may already equal one. If an elimination factor is reread after its column has been zeroed, it no longer represents the row operation. Saving both values makes every update use the intended row operation.

Gauss-Jordan control flow

flowchart TD
    A["Read coefficient matrix A, block B, and tolerance"] --> C{"Valid finite shapes and tolerance?"}
    C -- No --> D["Reject input"]
    C -- Yes --> E["Copy augmented matrix; compute scale; k = 0"]
    E --> F{"k less than n?"}
    F -- Yes --> G["Select largest absolute active-column pivot"]
    G --> H{"Pivot exceeds threshold?"}
    H -- No --> I["Report singular or numerically rejected system"]
    H -- Yes --> J["Swap complete augmented rows"]
    J --> K["Save pivot; normalize entire pivot row"]
    K --> L["Eliminate column from every other row"]
    L --> M["Increment k"]
    M --> F
    F -- No --> N["Return identity and transformed right-hand-side blocks"]

Worked Example with a Full Trace

Append $b$ and $I_2$ to the coefficient matrix. The two separators distinguish the system vector from the identity block:

$$ [A\mid b\mid I]=\left[\begin{array}{cc|c|cc}0&2&4&1&0\\1&3&7&0&1\end{array}\right]. $$

The scale is $\alpha=3$, giving threshold $3\times10^{-12}$.

Pivot column $0$. The candidate magnitudes are $0$ and $1$. Swap rows $0$ and $1$, carrying all five entries:

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

The saved pivot is $1$. Dividing the first row by $1$ leaves its entries $1,3,7,0,1$. The other row's elimination factor is $0$, so subtracting its zero multiple changes nothing. Column $0$ is now the first identity column.

Pivot column $1$. Only the second row remains active, with pivot $2$. Save $2$ and divide the full row by it. Its entries become $0/2=0$, $2/2=1$, $4/2=2$, $1/2$, and $0/2=0$:

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

Eliminate the entry above this pivot with factor $3$. The first row updates entry by entry to $1-3\cdot0=1$, $3-3\cdot1=0$, $7-3\cdot2=1$, $0-3(1/2)=-3/2$, and $1-3\cdot0=1$:

$$ [I\mid x\mid A^{-1}]=\left[\begin{array}{cc|c|cc}1&0&1&-3/2&1\\0&1&2&1/2&0\end{array}\right]. $$

There is no back substitution: the coefficient block is already the identity. Read $x=(1,2)^{\mathsf T}$ and

$$ A^{-1}=\begin{pmatrix}-3/2&1\\1/2&0\end{pmatrix}. $$

Direct substitution gives $Ax=(4,7)^{\mathsf T}$. The inverse check is

$$ AA^{-1}=\begin{pmatrix}0(-3/2)+2(1/2)&0\cdot1+2\cdot0\\1(-3/2)+3(1/2)&1\cdot1+3\cdot0\end{pmatrix} =\begin{pmatrix}1&0\\0&1\end{pmatrix}. $$

All these sample values are exactly representable in binary floating point, so the printed residuals are zero. General input matrices can produce nonzero rounding error.

Worked Example

Trace It by Hand

The second pivot row before normalization is $[0,2\mid4\mid1,0]$. Predict every entry after normalization, then the first row after eliminating its coefficient $3$ above the pivot.

Divide all five entries by $2$, giving $[0,1\mid2\mid1/2,0]$. Subtract three times this row from $[1,3\mid7\mid0,1]$, obtaining $[1,0\mid1\mid-3/2,1]$. Explain why scaling only the coefficient entries would produce a different solution even if the final coefficient block looked like an identity.

Implementation in Python, C++, and Java

Each panel is a complete program with the same gauss_jordan operation. It supports any positive number of right-hand-side columns for a square coefficient matrix that passes the pivot tolerance. It returns the full reduced augmented workspace, not just its right block. Singular, rank-deficient, or tolerance-rejected coefficient matrices cause an exception; the routine does not classify consistency or choose free variables.

The multiplication and error-check helpers assume compatible nonempty matrix dimensions. The sample extracts one solution column and the inverse block. It reports the maximum absolute entry of $Ax-b$ and of $AA^{-1}-I$ using the original inputs. The programs use double precision and tolerance $10^{-12}$, with residual threshold $10^{-10}$ for this small example.

All three print five identical lines: x = [1.000000, 2.000000]; inverse = [[-1.500000, 1.000000], [0.500000, 0.000000]]; solve residual = 0.000e+00; inverse residual = 0.000e+00; check = PASS. Save the Java panel as Main.java.

Finite input checks cannot prevent arithmetic overflow for extreme magnitudes. The routines target moderate, sensibly scaled inputs and do not estimate conditioning. For broader numerical work, use scaled residuals and condition estimates rather than interpreting an absolute sample threshold as a universal accuracy guarantee.

from math import isfinite


def gauss_jordan(a: list[list[float]], b: list[list[float]],
                 tol: float = 1e-12) -> list[list[float]]:
    n = len(a)
    if n == 0 or len(b) != n or any(len(row) != n for row in a):
        raise ValueError("require square A and matching B rows")
    q = len(b[0])
    if q == 0 or any(len(row) != q for row in b):
        raise ValueError("B must have a positive consistent column count")
    if not isfinite(tol) or not 0 < tol < 1:
        raise ValueError("require 0 < tol < 1")
    if any(not isfinite(v) for matrix in (a, b) for row in matrix for v in row):
        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)
    width = n + q
    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]
        divisor = aug[k][k]
        for j in range(width):
            aug[k][j] /= divisor
        for i in range(n):
            if i != k:
                factor = aug[i][k]
                for j in range(width):
                    aug[i][j] -= factor * aug[k][j]
                aug[i][k] = 0.0
    return aug


def multiply(a: list[list[float]], b: list[list[float]]) -> list[list[float]]:
    out = [[0.0] * len(b[0]) for _ in a]
    for i in range(len(a)):
        for k in range(len(b)):
            for j in range(len(b[0])):
                out[i][j] += a[i][k] * b[k][j]
    return out


def max_error(a: list[list[float]], b: list[list[float]]) -> float:
    return max(abs(a[i][j] - b[i][j]) for i in range(len(a)) for j in range(len(a[0])))


def show(a: list[list[float]]) -> str:
    return "[" + ", ".join("[" + ", ".join(f"{v:.6f}" for v in row) + "]" for row in a) + "]"


def main() -> None:
    a = [[0.0, 2.0], [1.0, 3.0]]
    b = [[4.0, 1.0, 0.0], [7.0, 0.0, 1.0]]
    reduced = gauss_jordan(a, b)
    x = [[row[2]] for row in reduced]
    inv = [row[3:] for row in reduced]
    solve_error = max_error(multiply(a, x), [[4.0], [7.0]])
    inverse_error = max_error(multiply(a, inv), [[1.0, 0.0], [0.0, 1.0]])
    print("x = [" + ", ".join(f"{row[0]:.6f}" for row in x) + "]")
    print(f"inverse = {show(inv)}")
    print(f"solve residual = {solve_error:.3e}")
    print(f"inverse residual = {inverse_error:.3e}")
    print("check = " + ("PASS" if max(solve_error, inverse_error) <= 1e-10 else "FAIL"))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>

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

Matrix gauss_jordan(const Matrix& a, const Matrix& b, double tol = 1e-12) {
    const std::size_t n = a.size();
    if (n == 0 || b.size() != n) throw std::invalid_argument("invalid row counts");
    const std::size_t q = b[0].size();
    if (q == 0) throw std::invalid_argument("B must have columns");
    if (!std::isfinite(tol) || tol <= 0 || tol >= 1) throw std::invalid_argument("invalid tolerance");
    Matrix aug(n, std::vector<double>(n + q));
    double scale = 0;
    for (std::size_t i = 0; i < n; ++i) {
        if (a[i].size() != n || b[i].size() != q) throw std::invalid_argument("invalid matrix shapes");
        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]));
        }
        for (std::size_t j = 0; j < q; ++j) {
            if (!std::isfinite(b[i][j])) throw std::invalid_argument("inputs must be finite");
            aug[i][n + j] = b[i][j];
        }
    }
    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]);
        const double divisor = aug[k][k];
        for (std::size_t j = 0; j < n + q; ++j) aug[k][j] /= divisor;
        for (std::size_t i = 0; i < n; ++i) if (i != k) {
            const double factor = aug[i][k];
            for (std::size_t j = 0; j < n + q; ++j) aug[i][j] -= factor * aug[k][j];
            aug[i][k] = 0;
        }
    }
    return aug;
}

Matrix multiply(const Matrix& a, const Matrix& b) {
    Matrix out(a.size(), std::vector<double>(b[0].size(), 0));
    for (std::size_t i = 0; i < a.size(); ++i)
        for (std::size_t k = 0; k < b.size(); ++k)
            for (std::size_t j = 0; j < b[0].size(); ++j) out[i][j] += a[i][k] * b[k][j];
    return out;
}

double max_error(const Matrix& a, const Matrix& b) {
    double error = 0;
    for (std::size_t i = 0; i < a.size(); ++i)
        for (std::size_t j = 0; j < a[0].size(); ++j) error = std::max(error, std::abs(a[i][j] - b[i][j]));
    return error;
}

void show(const Matrix& a) {
    std::cout << '[' << std::fixed << std::setprecision(6);
    for (std::size_t i = 0; i < a.size(); ++i) {
        if (i != 0) std::cout << ", ";
        std::cout << '[';
        for (std::size_t j = 0; j < a[0].size(); ++j) {
            if (j != 0) std::cout << ", ";
            std::cout << a[i][j];
        }
        std::cout << ']';
    }
    std::cout << ']';
}

int main() {
    const Matrix a = { {0, 2}, {1, 3}};
    const Matrix b = { {4, 1, 0}, {7, 0, 1}};
    const Matrix reduced = gauss_jordan(a, b);
    Matrix x(2, std::vector<double>(1)), inv(2, std::vector<double>(2));
    for (std::size_t i = 0; i < 2; ++i) {
        x[i][0] = reduced[i][2];
        for (std::size_t j = 0; j < 2; ++j) inv[i][j] = reduced[i][3 + j];
    }
    const double solve_error = max_error(multiply(a, x), Matrix{ {4}, {7}});
    const double inverse_error = max_error(multiply(a, inv), Matrix{ {1, 0}, {0, 1}});
    std::cout << std::fixed << std::setprecision(6) << "x = [" << x[0][0] << ", " << x[1][0] << "]\n";
    std::cout << "inverse = "; show(inv); std::cout << '\n';
    std::cout << std::scientific << std::setprecision(3)
              << "solve residual = " << solve_error << '\n' << "inverse residual = " << inverse_error << '\n';
    std::cout << "check = " << (std::max(solve_error, inverse_error) <= 1e-10 ? "PASS" : "FAIL") << '\n';
    return 0;
}

import java.util.Locale;

public class Main {
    public static double[][] gaussJordan(double[][] a, double[][] b, double tol) {
        int n = a.length;
        if (n == 0 || b.length != n) throw new IllegalArgumentException("invalid row counts");
        int q = b[0].length;
        if (q == 0) throw new IllegalArgumentException("B must have columns");
        if (!Double.isFinite(tol) || tol <= 0 || tol >= 1) throw new IllegalArgumentException("invalid tolerance");
        double[][] aug = new double[n][n + q];
        double scale = 0;
        for (int i = 0; i < n; ++i) {
            if (a[i].length != n || b[i].length != q) throw new IllegalArgumentException("invalid matrix shapes");
            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]));
            }
            for (int j = 0; j < q; ++j) {
                if (!Double.isFinite(b[i][j])) throw new IllegalArgumentException("inputs must be finite");
                aug[i][n + j] = b[i][j];
            }
        }
        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;
            double divisor = aug[k][k];
            for (int j = 0; j < n + q; ++j) aug[k][j] /= divisor;
            for (int i = 0; i < n; ++i) if (i != k) {
                double factor = aug[i][k];
                for (int j = 0; j < n + q; ++j) aug[i][j] -= factor * aug[k][j];
                aug[i][k] = 0;
            }
        }
        return aug;
    }

    public static double[][] multiply(double[][] a, double[][] b) {
        double[][] out = new double[a.length][b[0].length];
        for (int i = 0; i < a.length; ++i)
            for (int k = 0; k < b.length; ++k)
                for (int j = 0; j < b[0].length; ++j) out[i][j] += a[i][k] * b[k][j];
        return out;
    }

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

    private static void show(double[][] a) {
        System.out.print("[");
        for (int i = 0; i < a.length; ++i) {
            if (i != 0) System.out.print(", ");
            System.out.print("[");
            for (int j = 0; j < a[0].length; ++j) {
                if (j != 0) System.out.print(", ");
                System.out.printf(Locale.ROOT, "%.6f", a[i][j]);
            }
            System.out.print("]");
        }
        System.out.print("]");
    }

    public static void main(String[] args) {
        double[][] a = { {0, 2}, {1, 3}};
        double[][] b = { {4, 1, 0}, {7, 0, 1}};
        double[][] reduced = gaussJordan(a, b, 1e-12);
        double[][] x = new double[2][1], inv = new double[2][2];
        for (int i = 0; i < 2; ++i) {
            x[i][0] = reduced[i][2];
            for (int j = 0; j < 2; ++j) inv[i][j] = reduced[i][3 + j];
        }
        double solveError = maxError(multiply(a, x), new double[][]{ {4}, {7}});
        double inverseError = maxError(multiply(a, inv), new double[][]{ {1, 0}, {0, 1}});
        System.out.printf(Locale.ROOT, "x = [%.6f, %.6f]%n", x[0][0], x[1][0]);
        System.out.print("inverse = "); show(inv); System.out.println();
        System.out.printf(Locale.ROOT, "solve residual = %.3e%ninverse residual = %.3e%n", solveError, inverseError);
        System.out.println("check = " + (Math.max(solveError, inverseError) <= 1e-10 ? "PASS" : "FAIL"));
    }
}

Complexity, Correctness, and Invariants

For $n$ coefficient rows and $q$ right-hand sides, the augmented width is $n+q$. Each of the $n$ pivot stages normalizes one full row and updates $n-1$ other full rows. Consequently, successful elimination takes $\Theta(n^2(n+q))$ time, including $O(n^2)$ pivot searching. The copied augmented workspace uses $\Theta(n(n+q))$ storage. With a single right-hand side or $q=O(n)$, these become cubic time and quadratic storage.

The shown code updates every other row even when its factor is zero. For a fixed accepted shape, best, average, and worst operation counts therefore have the same asymptotic bounds. Rejection at the first pivot follows validation and copying and can finish in $\Theta(n(n+q))$ time. An average including rejected inputs requires an input distribution.

For algebraic correctness, let $E$ represent the cumulative product of elementary row operations. At every stage, the workspace equals $[EA\mid EB]$. Each accepted pivot permits a reversible normalization, and row subtraction is reversible as well. After all coefficient columns become identity columns, $EA=I$, so $E=A^{-1}$ and the right block is $EB=A^{-1}B$. This directly proves both the solve and inverse interpretations in exact arithmetic.

A second invariant explains why later stages preserve earlier columns. Before stage $k$, the first $k$ coefficient columns are the corresponding identity columns. Every active pivot row below those completed pivots has zeros in them, so swapping active rows, scaling the chosen row, and subtracting its multiples leave those columns intact. The new stage establishes one additional identity column.

Floating-point arithmetic perturbs these identities. Partial pivoting, residual checks, scaling, and conditioning analysis all have distinct roles. RREF is unique in exact arithmetic, but floating-point results and rank decisions depend on tolerance and scaling.

Method	Time	Workspace or storage	Best use and assumptions
Gauss–Jordan with $q$ right-hand sides	$\Theta(n^2(n+q))$	$\Theta(n(n+q))$ here	Nonsingular square coefficient block accepted by tolerance
Gaussian elimination plus back substitution	$\Theta(n^3)$ for one solve	$\Theta(n^2)$	Usually less arithmetic when only one solution is needed
Retained pivoted LU plus $q$ triangular solves	$\Theta(n^3+n^2q)$	$\Theta(n^2+nq)$	Factor once and reuse for right-hand sides
Augment with identity for inverse	$\Theta(n^3)$	$\Theta(n^2)$	Computes all inverse columns together
QR or SVD methods	Dimension-dependent	Dimension-dependent	Appropriate for least squares or numerical rank analysis

After $k$ completed stages, the first $k$ coefficient columns are identity columns and the entire augmented workspace records the same cumulative row operations on $A$ and every right-hand side. Each accepted normalized pivot adds one more identity column; at completion, $[EA\mid EB]=[I\mid A^{-1}B]$ in exact arithmetic.

Variants, Pitfalls, and Optimizations

Normalize the entire row. Dividing only the coefficient entries leaves the right-hand sides inconsistent with the transformed equations. In the example, $[0,2\mid4\mid1,0]$ must become $[0,1\mid2\mid1/2,0]$. Keeping the appended values unchanged would yield the wrong solution $(-5,4)^{\mathsf T}$, whose original-system residual has infinity norm $4$.

Eliminate above and below. Stopping after entries below each pivot become zero gives echelon or triangular form, not necessarily RREF. Gauss–Jordan removes every other entry in a pivot column. For a single ordinary solve, Gaussian elimination and back substitution often provide the needed answer with fewer operations.

General RREF needs column skipping. A rectangular or rank-deficient reducer keeps separate pivot-row and pivot-column indices. If a coefficient column has no usable pivot, it skips that column and may identify a free variable. The supplied function instead requires a pivot in every coefficient column and rejects otherwise. A row of zero coefficients with a nonzero right-hand side signals inconsistency in exact arithmetic; zero rows can also encode redundant equations. Numerical classification requires appropriate scaling and tolerance policies.

Keep the bookkeeping aligned. Swap the whole augmented row and restrict pivot searches to coefficient columns. Moving only part of a row changes the represented systems, while choosing pivots from the appended identity or system columns misidentifies variables. Save normalization divisors and elimination factors before overwriting their entries.

Do not confuse inversion with solving. The inverse block is useful for teaching the row-operation identity and for applications that explicitly need inverse entries. When the objective is merely solving one or many systems, a suitable retained factorization generally avoids forming the inverse. The same nonsingular square $A$ can be used with many $B$ columns in the present function, but repeated calls repeat elimination.

Arithmetic and structure matter. Exact rational RREF is useful for symbolic calculations but may involve large fractions. Finite-field elimination uses nonzero pivots and field inverses instead of floating tolerances. Sparse Gauss–Jordan can generate fill-in and lose sparsity; choose a specialized sparse solver when storage or scale requires one. A global relative threshold can reject badly scaled invertible matrices and cannot replace a condition estimate.

Normalize and eliminate across every augmented column, preserve whole-row swaps, and save the old divisor and factor before updating. The coefficient block must be reduced above and below each pivot; a tolerance-based rejection alone does not classify rank, consistency, or numerical accuracy.

Practice Problems and Real-World Use

Specify which columns are coefficients and which are right-hand sides before reducing. Verify each extracted result against the original data. The tasks below connect RREF, systems, and inverses while keeping small answers independently checkable.

Exercise or application	Main idea or check	Difficulty
Solve the worked system	Reduce $[A\mid b]$ to read $x=(1,2)^{\mathsf T}$ and check $Ax=b$	Easy
Compute and verify the inverse	Reduce $[A\mid I]$ and verify $AA^{-1}=I$ for the sample	Easy
Process multiple right-hand sides	For the sample, add $b_2=(2,4)^{\mathsf T}$ and recover its solution $(1,1)^{\mathsf T}$	Medium
Compare echelon and reduced forms	Explain why $\begin{pmatrix}1&3\\0&2\end{pmatrix}$ is triangular but not RREF	Easy
Extend to a redundant system	Design a general reducer with skipped columns and identify a free variable in a consistent rank-deficient example	Hard
Symbolic parameter relationships	Use exact rational RREF to express pivot variables in terms of free variables	Medium
Small calibration or circuit systems	Solve several measurement vectors with the same coefficients; compare Gauss–Jordan with factorization reuse	Medium

For independent verification, construct small well-conditioned matrices with known solution blocks and compare the entire returned workspace with $[I\mid X]$. Include row permutations, negative pivots, multiple right-hand-side widths, malformed blocks, nonfinite values, and singular inputs. Exact determinant-based solutions for tiny systems help catch errors that matching forward and inverse code could otherwise share.

Cookie Consent