[julia-dev] 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.

Post by 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.