import numpy as np
from lu import LU

def solve(A,b):
    '''
    Solve linear system using LU w/o pivoting
    '''
    A1=np.copy(A)
    LU(A1)
    y=np.zeros_like(b)
    y[0]=b[0]
    for i in range(1,len(b)):
        y[i]=b[i]-np.inner(A1[i,:i],y[:i])
    x=np.zeros_like(b)
    x[-1]=y[-1]/A1[-1,-1]
    for i in range(len(b)-2,-1,-1):
        x[i]=(y[i]-np.inner(A1[i,i+1:],x[i+1:]))/A1[i,i]
    return x

if __name__=='__main__':
    N = 100
    A = np.random.randn(N,N)
    x = np.random.randn(N)
    b = np.matmul(A,x)
    x1= solve(A,b)
    print('Relative forward  error: %e' % (np.linalg.norm(x-x1,np.Inf)/np.linalg.norm(x,np.Inf)))
    print('Relative backward error: %e' % (np.linalg.norm(b-np.matmul(A,x1),np.Inf)/np.linalg.norm(b,np.Inf)))

    # example 2.13 in the book
    epsilon = 1.0E-20
    A=np.array([[epsilon,1],[1,2]])
    b=np.array([1.,4.])
    x=solve(A,b)
    print(x)