Operator fusion/matrix chain multiplication at constant `(u, p, t)`-slices

[danielwe]

Say you have C = ComposedOperator(A::MatrixOperator, B::MatrixOperator). You're passing it to an iterative solver. Depending on the shapes of A and B, precomputing the product Cfused = A.A * B.A and passing that to the solver instead could be faster in some cases. You can't define C this way from the outset because A and B depend on (u, p, t), but you can always fuse ahead of any sequence of operations that lives on a constant (u, p, t)-slice.

Opt-in functionality to support this might be useful. For more complicated structures it may try to approximate the optimal solution or it may rely on a flag on each composite node to decide whether to fuse.

We touched on this in #12 (comment).

Low prioirity.

[ChrisRackauckas]

https://en.wikipedia.org/wiki/Matrix_chain_multiplication#:~:text=Matrix%20chain%20multiplication%20(or%20the,of%20the%20matrix%20multiplications%20involved.

Eagerly doing the matrix*vector product first is always the right thing to do when it's dense. When it's sparse, you may want to do a fused operation C first. And example of this is with PDE stencils, doing a 2-dimensional convolution with proper blocking is faster than doing two separate 1-dimensional stencils.