If you already have you matrix in component form, i.e. you have it broken into submatrices, then it is straightforward to construct a block diagonal matrix of any size.
Clear[BlockDiagonal]
BlockDiagonal[a : {_?MatrixQ ..}] :=
ArrayFlatten[
DiagonalMatrix[Range[Length@a]] /. i_Integer?Positive :> a[[i]]]
(* Convenience notation *)
Unprotect[CirclePlus];
a_?MatrixQ ⊕ b__?MatrixQ := BlockDiagonal[{a, b}]
Protect[CirclePlus];
For convenience,
Clear[ZeroMatrix]
ZeroMatrix[dim_Integer?Positive] := ConstantArray[0, {dim, dim}]
ZeroMatrix[dims : {_Integer?Positive, _Integer?Positive}] :=
ConstantArray[0, dims]
Which is used as follows:
A = {{a, b}, {c, d}};
B = {{e, f}, {g, h}};
A ⊕ ZeroMatrix[3]⊕B

Edit:
Per the comments, if you want to interlace the blocks from two different matrices, you can either specify the blocks individually, as above, or you can make use of a feature of SparseArray
to capture the blocks automatically.
Clear[BlockMatrixRiffle];
BlockMatrixRiffle[A_?MatrixQ, B_?MatrixQ]:=
Module[{nonzA, nonzB},
nonzA = FindClusters[SparseArray[A]["NonzeroPositions"], Method -> "Agglomerate"];
nonzB = FindClusters[SparseArray[B]["NonzeroPositions"], Method -> "Agglomerate"];
BlockDiagonal@Riffle[
A[[##]]&@@@Map[Span[Min[#], Max[#]]&, Transpose/@nonzA, {2}],
B[[##]]&@@@Map[Span[Min[#], Max[#]]&, Transpose/@nonzB, {2}]
]
]
(Edit 2: changed clustering method to Agglomerate which is hierarchical, and more appropriate for this case.)
SparseArray
only stores elements that are different from a common base which is usually set to zero. So, when you turn a block diagonal matrix into a SparseArray
you can retrieve the non-zero positions by
SparseArray[A ⊕ B]["NonzeroPositions"]
(* {{1, 1}, {1, 2}, {2, 1}, {2, 2}, {3, 3}, {3, 4}, {4, 3}, {4, 4}} *)
(Note: this will not find zero blocks, so if those exist, then some other method must be used.) Now, we use FindClusters
to group them into blocks:
FindClusters[%, Method -> "Agglomerate"]
(* {{{1, 1}, {1, 2}, {2, 1}, {2, 2}}, {{3, 3}, {3, 4}, {4, 3}, {4, 4}}} *)
From there, I turn them into min/max ranges for use in Part
:
Map[{Min[#], Max[#]}&, Transpose /@ %]
(* {{{1, 2}, {1, 2}}, {{3, 4}, {3, 4}}} *)
So, when I Apply
A[[##]]
to the above, I get the equivalent of
{A[[{1,2}, {1,2}]], A[[{3,4}, {3, 4}]]}
for both the matrices supplied to BlockMatrixRiffle
which then interlaces the two matrices.
mat1 = Table[a[i, j], {i, 2}, {j, 2}] ⊕ Table[b[i, j], {i, 2}, {j, 2}];
mat2 = Table[c[i, j], {i, 3}, {j, 3}] ⊕ Table[d[i, j], {i, 2}, {j, 2}];
BlockMatrixRiffle[mat1, mat2]

Note, if the number of blocks in B
exceeds the number of blocks in A
, those extra blocks will be ignored. But, if the reverse is true, multiple copies of the blocks from B
will be included in the result. This can be prevented by a little extra processing, but I leave that as an exercise.