I have an algorithm where I have to perform A * x and A \ x a lot of times for the same A::AbstractMatrix. Importantly, A can be dense, triangular, sparse, etc.
I would like to save a generic decomposition object for A which makes this fast and numerically stable.
First I tried lu(A), but that does not work for eg LowerTriangular. So now I am special-casing for these things. Is there a generic solution?
Maybe add a specialization for lu(::LowerTriangular) in the stdlib? It’s hard to imagine a better generic solution than lu.
Except maybe QR?
For unequal sizes, sure.
A is square (sorry, forgot to say).
In fact for A::Diagonal or A::AbstractTriangular, I can just leave it alone. Would an lu method that does this make sense? Or complement with an L = I. I am not sure.
Isn’t that the point of factorize?
That’s a good point, unfortunately it is not type stable. A version that would be less eager to find the best solution could work better for me.
If you’re factoring a matrix I’m guessing the dynamic dispatch overhead is the least of your problems. Can’t you get by with a function barrier or a parametric type?
Thanks, I will benchmarks and see how it goes. This is all for CI code by the way, where I am also checking inference, but I think putting the function barrier in the right place could work around this.