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?