Why doesn't lufact return Triangular Matrices for its data?

[Chris Rackauckas]

Why doesn't lufact return Triangular Matrices?

A = rand(1000,1000)

B = lufact(A,Val{false})

C = B[:U]
D = Base.LinAlg.UpperTriangular(B[:U])

The returned C is just a normal 1000x1000 array, not the same as D. I know there isn't a storage advantage to using Triangular matrices, but it seems there's a speed benefit to using Triangular matrices (for example I checked D*D instead of C*C), and maybe further improvements to Triangular matrices could be made after the array buffer implementation. 

[Andreas Noack]

The problem is that our Triangular types are square and U might not be square, e.g.

julia> lufact(randn(3,4))[:U]

3×4 Array{Float64,2}:

 -2.98058  -0.234937   1.49172   1.7675

  0.0      -1.15066   -1.72942   1.26475

  0.0       0.0       -1.58043  -0.0808058

We could loosen then type to allow for trapezoid matrices but it would make them more complicated and require more checks.