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

Skip to first unread message

Chris Rackauckas

Sep 27, 2016, 12:08:15 PM9/27/16
to julia-dev
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

Sep 27, 2016, 7:33:19 PM9/27/16
to juli...@googlegroups.com
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.
Reply all
Reply to author
0 new messages