one-file-projects/gauss.py

def gauss(A,b,offset=(0,0)):
    ox,oy = offset
    for i in range(oy,len(A)):
        if A[i][ox] != 0:
            A[i], A[oy] = A[oy], A[i]
            b[i], b[oy] = b[oy], b[i]
            break

    pprint(A)

    for i in range(oy+1,len(A)):
        if A[i][ox] != 0:
            c = A[i][0]/A[oy][0]
            A[i] = [ (A[i][x] - A[oy][x]*c) for x in range(ox,len(A[i])) ]
            b[i] = b[i] - b[oy]*c

    if oy < len(A)-1:
        A,b = gauss(A,b,offset=(ox+1,oy+1))

    b[oy] /= A[oy][ox]
    x = b[oy]
    for y in range(oy):
        b[y] -= A[y][ox]*x

    return A,b


def vprint(A):
    print("Vector ----")
    for i in range(len(A)):
        print("%10.4f" % A[i])
    print("-----")

def pprint(A):
    print("Matrix ----")
    for i in range(len(A)):
        print( " ".join(["%10.4f" % (A[i][j]) for j in range(len(A[i])) ]) )
    print("-----")

A = [ [    0,     0, 0,-1,-1, 1 ],
      [   -1,     1, 0, 0, 0, 0 ],
      [    0,    -1, 1, 0, 0, 0 ],
      [-1000,     0, 0, 1, 0, 0 ],
      [    0, -1000, 0, 0, 1, 0 ],
      [    0,     0, 0, 0, 0, 1 ] ]

b = [ 0,
      0,
      0,
      0,
      0,
      5 ]

pprint(A)
vprint(b)


# _, x = gauss(A,b)

from numpy.linalg import solve
x = solve(A,b)

vprint(x)