import numpy as np
import laff
import flame

# Applying the LU Factorization to a matrix is the process of
# computing a unit lower triangular matrix, L, and upper 
# triangular matrix, U, such that A = L U.  To avoid nasty fractions
# in these caclulations, we create our matrix A from two matrices, L and U, 
# whose elements we know to be integer valued.

L = np.matrix( ' 1, 0, 0. 0;\
                -2, 1, 0, 0;\
                 1,-3, 1, 0;\
                 2, 3,-1, 1' )

U = np.matrix( ' 2,-1, 3,-2;\
                 0,-2, 1,-1;\
                 0, 0, 1, 2;\
                 0, 0, 0, 3' )

A = L * U

print( 'A = ' )
print( A )

# create a solution vector x

x = np.matrix( '-1;\
                 2;\
                 1;\
                -2' )

# store the original value of x

xold = np.matrix( np.copy( x ) )

# create a solution vector b so that A x = b

b = A * x
print( 'b = ' )
print( b )

# Much later, we are also going to solve A x = y.  Here we create that y:
x2 = np.matrix( '1;\
                 2;\
                 -2;\
                 1' )
y = A * x2
print( 'y = ' )
print( y )

# Insert code here

# recreate matrix A
A = L * U

# recreate the right-hand side
b = A * xold

# apply Gaussian elimination to matrix A
LU_unb_var5( A )


print( 'A after factorization' )
print( A )
print( 'Original L' )
print( L )
print( 'Original U' )
print( U )

# Insert code here

print( A )
print( b )
Ltrsv_unb_var1( A, b )

print( 'updated b' )
print( b )

x[ 3 ] = b[ 3 ] / A[ 3,3 ]
x[ 2 ] = ( b[ 2 ] - A[ 2,3 ] * x[ 3 ] ) / A[ 2,2 ]
x[ 1 ] = ( b[ 1 ] - A[ 1,2 ] * x[ 2 ] - A[ 1,3 ] * x[ 3 ] ) / A[ 1,1 ]
x[ 0 ] = ( b[ 0 ] - A[ 0,1 ] * x[ 1 ] - A[ 0,2 ] * x[ 2 ]- A[ 0,3 ] * x[ 3 ] ) / A[ 0,0 ]

print( 'x = ' )
print( x )

print( 'x - xold' )
print( x - xold )

# Insert code here

# just to be sure, let's start over.  We'll recreate A, x, and b, run all the routines, and
# then compare the updated b to the original vector x.

A = L * U
b = A * x

LU_unb_var5( A )
Ltrsv_unb_var1( A, b )
Utrsv_unb_var1( A, b )

print( 'updated b' )
print( b )
print( 'original x' )
print( x )
print( 'b - x' )
print( b - x )


def Solve( A, b ):
    
    # insert appropriate calls to routines you have written here!
    
    

# just to be sure, let's start over.  We'll recreate A, x, and b, run all the routines, and
# then compare the updated b to the original vector x.

A = L * U
b = A * x

Solve( A, b )

print( 'updated b' )
print( b )
print( 'original x' )
print( x )
print( 'b - x' )
print( b - x )


# insert appropriate calls here.


print( 'y = ' )
print( y )

print( 'x2 - y =' )
print( x2 - y )