from __future__ import division
from pylab import *

def lu(A):
    """LU decomposition of A without pivoting (Doolittle's method)."""
    n = shape(A)[0]
    L = eye(n)
    U = zeros([n,n])
    for k in range(n):
        U[k,k:] = A[k,k:] - dot(L[k,:k], U[:k,k:])        
        L[k+1:,k] = (A[k+1:,k] - dot(L[k+1:,:k], U[:k,k]))/U[k,k]
    return L, U

def uSolve(U,b):
    """Solve an upper triangular system Ux = b by back substitution."""
    n = size(b)
    x = zeros(n)
    # the reversed() function lets you iterate backwards over a list
    for i in reversed(range(n)):
        x[i] = (b[i] - dot(U[i,i+1:], x[i+1:]))/U[i,i]
    return x
    
def lSolve(L,b):
    """Solve a lower triangular system Lx = b by forward substitution."""
    n = size(b)
    x = zeros(n)
    for i in range(n):
        x[i] = 0 # replace this line
    return x