# Insert your code for LU_unb_var5 from
# 6.3 Solving A x b via LU factorization and triangular solves

import numpy as np
import laff
import flame

# Form the matrices
m = 200

L = np.matrix( np.tril(np.random.rand(m,m), -1 ) + np.eye(m) )
U = np.matrix( np.triu(np.random.rand(m,m) ) + np.eye(m) * m )

A = L * U

# Create a large, random solution vector x
x = np.matrix( np.random.rand(m,1) )

# 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


# Later, we are also going to solve A x = b2.  Here we create that b2:
x2 = np.matrix( np.random.rand(m,1) )
b2 = A * x2

# Set the block size, we'll use it later. 
# Since we're just messing around with blocked algorithms,
# we set this totally arbitrarily
nb = 20

# Insert code here

# recreate matrix A
A = L * U

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

# apply blocked LU to matrix A
# remember nb holds our block size
LU_blk_var5( A, nb )

# Compare A to the original L and U matrices
print( 'Maximum value of (A - L - U) after factorization' )
print( np.max( np.abs( A - np.tril(L,-1) - U ) ) ) #The "-1" ignores the diagonal

def Solve( A, b ):
    
    # insert appropriate calls to routines you have written here!
    # remember the variable nb holds our block size

# 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( '2-Norm of Updated b - original x' )
print( np.linalg.norm(b - x) )


# insert appropriate calls here

# Test your calls
print( '2-Norm of updated b2 - original x2' )
print( np.linalg.norm(b2 - x2) )