sage: A = matrix(CDF, [[1, 2, 3], [4, 5, 6]])
sage: C = transpose(A)*A
sage: P, L, U = C.LU()
sage: P
[0.0 0.0 1.0]
[1.0 0.0 0.0]
[0.0 1.0 0.0]
sage: L
[ 1.0 0.0 0.0]
[ 0.62962962963 1.0 0.0]
[0.814814814815 0.5 1.0]
sage: U
[ 27.0 36.0 45.0]
[ 0.0 -0.666666666667 -1.33333333333]
[ 0.0 0.0 -0.0]
sage: L*U
[27.0 36.0 45.0]
[17.0 22.0 27.0]
[22.0 29.0 36.0]
sage: P*C
[27.0 36.0 45.0]
[17.0 22.0 27.0]
[22.0 29.0 36.0]
sage: transpose(P)*L*U
[17.0 22.0 27.0]
[22.0 29.0 36.0]
[27.0 36.0 45.0]
sage: C
[17.0 22.0 27.0]
[22.0 29.0 36.0]
[27.0 36.0 45.0]
sage: transpose(P)*L*U == C
False
sage: transpose(P)*L*U - C
[-3.5527136788e-15 0.0 0.0]
[ 0.0 0.0 0.0]
[ 0.0 0.0 0.0]
sage: O = matrix(CDF, 3, 3, 0)
sage: transpose(P)*L*U - C == O
False
sage: transpose(P)*L*U - C - O
[-3.5527136788e-15 0.0 0.0]
[ 0.0 0.0 0.0]
[ 0.0 0.0 0.0]