using LinearAlgebra, BlockDiagonals
A = rand(3, 2)
B = rand(2, 1)
C = rand(3, 1)
mul!(C, A, B) # works
A = BlockDiagonal([A, A])
B = BlockDiagonal([B, B])
C = BlockDiagonal([C, C])
mul!(C, A, B) # errors
# what it should be:
for i in eachindex(C.blocks)
mul!(C.blocks[i], A.blocks[i], B.blocks[i])
end
C == A*B
To fix this, we should change the isequal_blocksizes check with an actual proper check that the sizes work together.
Alternatively a more lazy version would just check that length(C.blocks) == length(A.blocks) == length(B.blocks), and then just call mul! which performs a size check anyways.
The
isequal_blocksizes
check here https://github.com/JuliaArrays/BlockDiagonals.jl/blob/9b07a077d37510184857659d77666832db19b2eb/src/linalg.jl#L141-L142 is much too strict. See for example this code:To fix this, we should change the
isequal_blocksizes
check with an actual proper check that the sizes work together.Alternatively a more lazy version would just check that
length(C.blocks) == length(A.blocks) == length(B.blocks)
, and then just callmul!
which performs a size check anyways.