Gaussian Elimination with Partial Pivoting

Gaussian elimination is a direct method for solving a linear system of equations.

A linear system is a set of simultaneous equations (linear) in several variables. In theory, solving such a system algebraically is straightforward. First, we eliminate the first variable either by substitution or appropriate row operations. We then continue to eliminate the second variable and so on. Finally, we’re left with one equation for the final variable which we can easily solve. We then work backwards, finding the remaining unknowns one by one.

When using floating point arithmetic, there can be pitfalls to this method. Each floating point operation introduces a small roundoff error. These errors may seem negligable when considered in isolation, but they can (depending on the characteristics of the matrix) build up during the course of the calculation and lead to significant computational errors.

The formal analysis of error growth is a key topic in Numerical Analysis, and one of the pioneers of such analysis was J.H. Wilkinson.

One of the refinements to the Gaussian Elimination method which Wilkinson studied was Partial Pivoting. This modification costs very little in terms of computational overhead but is extremely effective in controlling error growth. For example, see How Accurate is Gaussian Elimination, by N Higham.

The Problem

For the case where we have $n$ equations in $n$ unknowns, the problem is written mathematically as follows.

Given matrix $A \in R^{n \times n}$ and vector $b \in R^n$, find $x \in R^n$ such that:
$Ax = b$

Implementation

We shall now conclude with an example C++ implementation of GEPP (Gaussian Elimination with Partial Pivoting).

#include <iostream>
#include <utility>

template <int N>
void MatrixSolve(double* x, double A[][N], double* b)
{
    int i, j, k;
    double U[N][N+1];

    // Copy A to U and augment with vector b.
    for (i = 0; i < N; i++)
    {
        U[i][N] = b[i];
        for (j = 0; j < N; j++)
            U[i][j] = A[i][j];
    }
    
    // Factorisation stage
    for (k = 0; k < N; k++)
    {
        // Find the best pivot
        int p = k;
        double maxPivot = 0.0;
        for (i = k; i < N; i++)
        {
            if (fabs(U[i][k] > maxPivot))
            {
                maxPivot = fabs(U[i][k]);
                p = i;
            }
        }

        // Swap rows k and p
        if (p != k)
        {
            for (i = k; i < N + 1; i++)
                std::swap(U[p][i], U[k][i]);
        }

        // Elimination of variables
        for (i = k + 1; i < N; i++)
        {
            double m = U[i][k] / U[k][k];
            for (j = k; j < N + 1; j++)
                U[i][j] -= m * U[k][j];
        }
    }

    // Back substitution
    for (int k = N - 1; k >= 0; k--)
    {
        double sum = U[k][N];
        for (int j = k + 1; j < N; j++)
            sum -= U[k][j] * x[j];
        x[k] = sum / U[k][k];
    }
}