from numpy import array
from scipy.linalg import lu
# define a square matrix
A = array([[3,2,1], [9,3,4], [4,2,5]])
print(A)
# LU decomposition
P, L, U = lu(A)
print(P)
print(L)
print(U)
# reconstruct
B = P.dot(L).dot(U)
print(B)