Pasted

import numpy
 
eps = 2**-52
 
def LU(A):
    U = numpy.matrix(A)
    n = U.shape[0]
    assert len(U.shape)==2
    assert n==U.shape[1]
    L = numpy.matrix(numpy.identity(n))
    for j in range(n-1):
        assert abs(U[j,j])>=eps
        for i in range(j+1, n):
            mult=U[i,j]/U[j,j]
            for k in range(n):
                U[i,k] -= mult*U[j,k]
            L[i,j]=mult
    return L, U
    
def backsolve(A, b):
    b = numpy.matrix(b).transpose()
    L, U = LU(A)
    n = L.shape[0]
    assert b.shape[0]==n
    assert b.shape[1]==1
    c = numpy.matrix(numpy.zeros((n, 1)))
    x = c.copy()
    for i in range(n):
        mult = c[i] = b[i]
        for k in range(i,n):
            b[k]-= mult*L[k,i]
    for i in range(n-1, -1, -1):
        mult = x[i] = c[i]/U[i,i]
        for k in range(i):
            c[k]-=mult*U[k,i]
    return x