I have a code that does something like this:

Z = inv(ABCDE)

Here, A is a diagonal matrix. B is block diagonal, its diagonal blocks are 2x2. C is block diagonal, its diagonal blocks are 4x4. Subsequent matrices have larger blocks, but the whole thing is “compatible” so that ABCDE is indeed block diagonal, and the inverse Z is therefore also block diagonal.

In MATLAB, I’d use sparse matrices throughout. In my experience, if a sparse matrix is block diagonal, MATLAB seems to do the right thing.

In Julia, the BlockDiagonals package doesn’t cope too well with these increasing block sizes. The sparse matrix package does, but I have yet to find a way to convince it to invert the resulting block diagonal matrix.

My favorite solution would be to somehow get Julia to invert these sparse matrices that happen to be block diagonal, but if I cannot do it, I’ll have to implement my own BlockDiag class that better understands blocks. Does anyone know how to efficiently invert a block diagonal matrix encoded as a sparse matrix in Julia?

Thanks,

S