When solving Ax=b with a sparse A, the output (along with the inverse A^{-1}) is generically dense even if b is sparse. (See also the thread Sparse solve with sparse rhs - #3 by stevengj) Julia’s generic sparse-matrix type assumes that this is the case and therefore always computes a dense output from A \ b or lu(A) \ b.
Of course, for the extremely special sparsity pattern of block-diagonal A this is not the case, and you indeed have a sparse (block-diagonal) inverse and sparse solutions for sparse right-hand-sides. But Julia’s SparseMatrixCSC will not exploit such “structured sparsity” — it is designed exclusively for unstructured sparsity.
Offhand, I can see at least three options for you:
- Write your own block-diagonal arithmetic specialized for your case.
- Submit a PR to BlockDiagonals.jl to support multiplying block-diagonal matrices with unequal — but compatible — block sizes, e.g. m \times m blocks with km \times km blocks. This is more work in some ways because the BlockDiagonals package probably supports more generic block sizes than you need, but may have the greatest long-term benefit.
- Multiply the matrices as unstructured sparse matrices, but write a constructor that then converts the result to block-diagonal. (This is probably the least-efficient option.)