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`

.

1 Like

`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`

?

1 Like

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.