A naive implementation of Gaussian elimination in Python

gaussian.py

def display(matrix):
    """Simple output for matrix."""
    
    for row in matrix:
        print(repr(row))

def leading(row):
    """Return the leading index of a row."""

    for i, element in enumerate(row):
        if element != 0:
            return i

def augment(matrix, solution):
    """Add a solution element to each row."""

    for row, scalar in zip(matrix, solution):
        row.append(scalar)

def annihilate(matrix, row, column):
    """Annihilate the column of each row below that specified."""

    for i in range(row + 1, len(matrix)):
        factor = matrix[i][column] / matrix[row][column]
        for j in range(column, len(matrix[i])):
            matrix[i][j] -= factor * matrix[row][j]

def reduce(matrix):
    """Fill backwards to convert to RREF."""

    for i in reversed(range(len(matrix))):  # Row we are using to reduce
        index = leading(matrix[i])
        for j in reversed(range(i)):  # Row we are reducing
            factor = matrix[j][index] / matrix[i][index]
            for k in range(index, len(matrix[j])):
                matrix[j][k] -= factor * matrix[i][k]

def build(augmented):
    """Factor out the solutions from the augmented matrix."""

    final = []
    for row in augmented:
        final.append([row[-1] / row[leading(row)]])
    return final

def solve(matrix, solution):
    """Solve a matrix equation Ax = b for x."""

    matrix = [row[:] for row in matrix]
    print("Starting matrix:")
    display(matrix)
    print()

    augment(matrix, solution)
    print("Augmented:")
    display(matrix)
    print()

    for i in range(len(matrix)):
        annihilate(matrix, i, i)
    print("Row echelon:")
    display(matrix)
    print()

    reduce(matrix)
    print("Reduced row echelon:")
    display(matrix)
    print()

    result = build(matrix)
    print("Solution:")
    display(result)
    print()


A = [[ 1.0347,  0.8884,  1.4384, -0.1022, -0.0301],
     [ 0.7269, -1.1471,  0.3252, -0.2414, -0.1649],
     [-0.3034, -1.0689, -0.7549,  0.3192,  0.6277],
     [ 0.2939, -0.8095,  1.3703,  0.3129,  1.0933],
     [-0.7873, -2.9443, -1.7115, -0.8649,  1.1093]]

b = [-0.8637, 0.0774, -1.2141, -1.1135, -0.0068]

solve(A, b)

# Test
# import numpy
# AA = numpy.matrix(A)
# xx = numpy.matrix(x)
# print(AA * xx)