Simple gauss-jordan matrix solver. Somewhat useful; might expand in the future or something.

linalg.py

def Matrix (T):
    class M:
        def __init__ (self, *args):
            if (type(args[0]) == type(self)):
                self.rows = [ [ x for x in row ] for row in args[0].rows ]
            elif type(args[0] == str):
                self.rows = [ map(T, row.strip().split()) for row in args[0].strip().split('\n') ]
            else:
                self.rows = [ map(T, row) for row in args[0] ]
        def __repr__ (self):
            return "\n\t%s\n"%("\n\t".join([ "\t".join([ "%4s"%x for x in row ]) for row in self.rows ]))
        def add_row (self, i, j, k):
            self.rows[i] = [ a + b * k for a, b in zip(self.rows[i], self.rows[j]) ]
        def scale_row (self, i, k):
            self.rows[i] = [ a * k for a in self.rows[i] ]
        def swap_rows (self, i, j):
            self.rows[i], self.rows[j] = self.rows[j], self.rows[i]

        def reduce_gauss_jordan (self, operator):
            N = len(self.rows)
            for i in range(N):
                if self.rows[i][i] == 0:
                    j = i
                    while self.rows[j][i] == 0 and j < N:
                        j += 1
                    if n < N:
                        operator.swap(i, j)
                    else:
                        continue
                if self.rows[i][i] != 1:
                    operator.scale(i, 1 / self.rows[i][i]) 
                for j in range(N):
                    if i != j:
                        operator.add(j, i, -self.rows[j][i])
        def clone (self):
            return M(self)
        def rref (self):
            class Operator:
                def __init__(self, matrix):
                    self.matrix = matrix
                def scale (self, i, k):
                    self.matrix.scale_row(i, k)
                def swap (self, i, j):
                    self.matrix.swap_rows(i, j)
                def add(self, i, j, k):
                    self.matrix.add_row(i, j, k)
                def apply (self):
                    self.matrix.reduce_gauss_jordan(self)
                    return self.matrix
            return Operator(self.clone()).apply()

        def generate_rref(self):
            class Operator:
                def __init__(self,matrix):
                    self.matrix = matrix
                def scale (self, i, k):
                    print("Scale row %s by %s"%(i, k))
                    self.matrix.scale_row(i, k)
                    print(self.matrix)
                def swap (self, i, j):
                    print("Swap rows %s, %s"%(i, j))
                    self.matrix.swap_rows(i, j)
                    print(self.matrix)
                def add (self, i, j, k):
                    print("Add row %s += row %s * %s"%(i, j, k))
                    self.matrix.add_row(i, j, k)
                    print(self.matrix)
                def apply (self):
                    print("Initial: (generating RREF via gauss-jordan)")
                    print(self.matrix)
                    self.matrix.reduce_gauss_jordan(self)
                    return self.matrix
            return Operator(self.clone()).apply()
    return M

# Can subsitute any type in here - float, decimal, Fraction, some custom type that builds an AST, etc...
from fractions import Fraction
M = Matrix(Fraction)('''
    1 -1 -1 0 0
    -8 -6 0 1 0
    -8 0 -3 1 -1
''')

if __name__ == '__main__':
    print("Initial matrix: %s"%M)
    print("Solved matrix in RREF: %s"%M.rref())
    print("Result: %s"%(M.generate_rref()))

output.txt

Initial matrix: 
       1      -1      -1       0       0
      -8      -6       0       1       0
      -8       0      -3       1      -1

Solved matrix in RREF: 
       1       0       0    -1/10   1/15
       0       1       0    -1/30   -4/45
       0       0       1    -1/15   7/45

Initial: (generating RREF via gauss-jordan)

       1      -1      -1       0       0
      -8      -6       0       1       0
      -8       0      -3       1      -1

Add row 1 += row 0 * 8

       1      -1      -1       0       0
       0     -14      -8       1       0
      -8       0      -3       1      -1

Add row 2 += row 0 * 8

       1      -1      -1       0       0
       0     -14      -8       1       0
       0      -8     -11       1      -1

Scale row 1 by -1/14

       1      -1      -1       0       0
       0       1     4/7    -1/14      0
       0      -8     -11       1      -1

Add row 0 += row 1 * 1

       1       0    -3/7    -1/14      0
       0       1     4/7    -1/14      0
       0      -8     -11       1      -1

Add row 2 += row 1 * 8

       1       0    -3/7    -1/14      0
       0       1     4/7    -1/14      0
       0       0    -45/7    3/7      -1

Scale row 2 by -7/45

       1       0    -3/7    -1/14      0
       0       1     4/7    -1/14      0
       0       0       1    -1/15   7/45

Add row 0 += row 2 * 3/7

       1       0       0    -1/10   1/15
       0       1     4/7    -1/14      0
       0       0       1    -1/15   7/45

Add row 1 += row 2 * -4/7

       1       0       0    -1/10   1/15
       0       1       0    -1/30   -4/45
       0       0       1    -1/15   7/45

Result: 
       1       0       0    -1/10   1/15
       0       1       0    -1/30   -4/45
       0       0       1    -1/15   7/45

[Finished in 0.2s]