I have an algorithm where I have to perform
A * x and
A \ x a lot of times for the same
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
A is square (sorry, forgot to say).
In fact for
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
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.