'''
This script solves a random nxn-system AX=B, first by a standard algorithm,
and then by using an LU-factorization A=LU and solving the two triangular
systems LY=B and then UX=Y.

Note that LU is a lot faster, but that finding the LU-factorization takes
a long time. For this reason, LU is more efficient if we have to solve systems
AX=B many times with the same A but for different right sides B
'''
import numpy as np
from scipy.linalg import lu, solve_triangular
import time

# Generate a random nxn-matrix and a nx1 vector
n = 5_000
np.random.seed(0)
A = np.random.rand(n, n)
B = np.random.rand(n)

# Perform LU decomposition
start = time.time()
P, L, U = lu(A) #gives A=PLU
new_B = np.dot(np.linalg.inv(P), B) #rewrite as LU=P^{-1}B
decomposition_time = time.time() - start
print(f"Time for performing LU: {decomposition_time} seconds")
        
#Solve with standard numpy.linalg.solve
start = time.time()
X_standard = np.linalg.solve(A, B)
standard_time = time.time() - start
print(f"Time for standard method: {standard_time} seconds")

#Solve using LU
start = time.time()
Y = solve_triangular(L, new_B, lower=True)
X_lu = solve_triangular(U, Y)
lu_time = time.time() - start
print(f"Time for LU method: {lu_time} seconds")

#Uncomment to print solutions
#print("Solutions using standard algorithm: \n", X_standard)
#print("Solutions using LU-factorization: \n",X_lu)