How to make rotationally invariant solves fast?

Rotationally invariant operators are typically represented by ModalInterlace which combines matrices acting only on Fourier modes into a diagonal-block-banded matrix (that is, block banded matrix with diagonal blocks, i.e. subblockbandwidths are (0,0)).

At the moment we aren't taking advantage of this structure. There are two approaches:

1. Do it at the level of factorizing a BandedBlockBandedMatrix by checking whether the subblockbandwidths are (0,0). The easiest way would be to copy to a banded matrix, do a QR, and then copy the data back.
2. Introduce a DiagonalBlockBandedMatrix though this doesn't seem to have any obvious benefits over (1).
3. Do it directly on a ModalInterlace. This has the benefit of being able to call qr on each of the banded operators separately (and potentially taking advantage of @threads). Though probably (1) could also be parallelised.

@TSGut @ioannisPApapadopoulos any thoughts? Do you think this is important right now? It won't help with the variable coefficients problem but could be important for fractional DEs and time-stepping.

Start by reading ModalInterlace, BandedBlockBandedMatrix factorization, and the qr paths mentioned in the issue. Compare the three proposed approaches, including whether copying to a banded matrix or parallelizing separate banded operators is viable. Done means selecting and implementing a fast path for rotationally invariant solves, with evidence that it improves the relevant fractional differential equation or time-stepping workloads.