Block diagonal and sparse matrices

[stevengj]

@loisel

A while back, I was writing a preconditioner for PetSc and I was using the various matrix primitives built into PetSc. I think I multiplied some sparse matrices together. I knew that theoretically this could be done efficiently, and I could have written special-purpose code to make sure it was done efficiently, but I naively assumed PetSc’s primitive would be clever enough and unfortunately it didn’t work out. That preconditioner ended up being slow.

Side note: If you have a preconditioner ABC that is a product of sparse matrices, normally I would express this as a linear operator x \mapsto A(B(Cx)) to avoid the risk of losing sparsity in the product. (I think in Petsc you can do this as a “composite” preconditioner, or by defining a custom “shell” preconditioner with your own matvec. In Julia, you can implement such “lazy” products easily with LinearMaps.jl or LinearOperators.jl.)