import numpy

def backsubstitute(A, b):
    n,i = b.shape
    if i != 1:
        raise ValueError("b should be shape (n, 1), but is %s" % (b.shape,))
    i,j = A.shape
    if i != n or j != n:
        raise ValueError("A should be shape (%d, %d), but is %s" % (n,n,A.shape))

    x = b.copy()

    for i in range(n-1,-1,-1):
        x[i] = b[i]
        for j in range(n-1, i, -1):
            x[i] -= A[i,j] * x[j]
        x[i] /= A[i,i]

    return x

def gauss_eliminate(A_in, b_in):
    # Kopien der Matrix und des Vektors anlegen
    A = A_in.copy()
    b = b_in.copy()

    # Gausselimination durchfuehren
    # ....

    return A, b

def solve(A, b):
    Ares, bres = gauss_eliminate(A, b)
    x = backsubstitute(Ares, bres)
    return x

N = 5
A1 = numpy.random.random((N,N))
A2 = A1.copy()
A2[0,0] = 0.0

b2 = b1 = numpy.random.random((N,1))

import scipy.linalg

print "----------"
print "Solve eq 1"
print "----------"
print "A1 * x1 = b1"
print "A1 = %s" % A1
print "b1 = %s" % b1
x1_scipy = scipy.linalg.solve(A1, b1)
print "x1_scipy = %s" % x1_scipy
x1_own = solve(A1, b1)
print "x1_own = %s" % x1_own

diff = abs(x1_own - x1_scipy)
print "Maximal error = %s" % diff.max()
print

print "----------"
print "Solve eq 2"
print "----------"
print "A2 * x2 = b2"
print "A2 = %s" % A2
print "b2 = %s" % b2
x2_scipy = scipy.linalg.solve(A2,b2)
print "x1_scipy = %s" % x2_scipy
x2_own = solve(A2,b2)
print "x2_own = %s" % x2_own

diff = abs(x1_own - x1_scipy)
print "Maximal error = %s" % diff.max()