# Block diagonal and sparse matrices

**URL:** <https://discourse.julialang.org/t/block-diagonal-and-sparse-matrices/105551>\
**Category:** Numerics\
**Tags:** linearalgebra\
**Created:** [October 29, 2023, 11:13pm UTC](https://discourse.julialang.org/t/block-diagonal-and-sparse-matrices/105551 "2023-10-29T23:13:12Z")\
**Posts on this page:** 1\
**Showing post:** 5

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 30, 2023, 3:14pm UTC](https://discourse.julialang.org/t/block-diagonal-and-sparse-matrices/105551/5 "2023-10-30T15:14:26Z")

</div>

> [@loisel](#):
>
> I obtain a sparse LU decomposition, but then L\F is dense even if F is sparse.

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](https://discourse.julialang.org/t/sparse-solve-with-sparse-rhs/81595/3)) 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](https://github.com/JuliaArrays/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.)

---

_[View the full topic](https://discourse.julialang.org/t/block-diagonal-and-sparse-matrices/105551)._