A = matrix([[2,1,0],[1,2,1],[0,1,2]])
L = matrix([[1,0,0],[1/2,1,0],[0,2/3,1]])
U = matrix([[2,1,0],[0,3/2,1],[0,0,4/3]])
print "L*U is",L*U