Block diagonal and sparse matrices

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.)