gauss.py

gauss.py

from numpy import *

def triangularize(A):
    while A.size > 0:
        i = abs(A[:, 0]).argmax() # find pivot row
        if A[i, 0]: # false only if first column is all 0s
            A[[0, i]] = A[[i, 0]] # swap first row and pivot row
            A[0, :] /= float(A[0, 0]) # normalize first row
            A[1:, :] -= A[0, :] * A[1:, :] # eliminate first column
        A = A[1:, 1:] # advance to lower-right submatrix

def gauss(A):
    triangularize(A) # lower triangularize
    triangularize(A.transpose()) # upper triangularize

# examples

A = array([[1, 2, 0, 0], [0, 0, 1, -1], [0, 0, 0, 0]])
print A
gauss(A)
print A