Avoiding temporaries in hvcat with scalar functions

Many times I need to build matrices by hvcating combinations of other matrices, for example

M = [ω*I-A -D; -D' B-C]

for some general matrices A, B, C, D. Let’s assume they are dense (i.e. of type Matrix), but quite large. Building M like above creates many temporaries, because each block is materialized before being passed to hvcat. But note that the operations that define the blocks are all scalar functions. That is, they can be written as

M = [f₁₁.(M₁₁s...) f₁₂.(M₁₂s...); f₂₁.(M₂₁s...) f₂₂.(M₂₂s...)]

for some matrices Mᵢⱼs... and functions fᵢⱼ. This allows in principle to build M without materializing any intermediate matrices. We can do that e.g. by preallocating a matriz M = Matrix{promote_type(...)}(undef, computed_dimensions...) and doing a loop for each block, filling in element by element. That approach is fast, but of course quite horrible to write and read. Ideally I would like to be able to write

M = @lazycat [ω*I-A -D; -D' B-C]
M = @lazycat [f₁₁(M₁₁s...) f₁₂(M₁₂s...); f₂₁(M₂₁s...) f₂₂(M₂₂s...)]

My question: anybody knows of a package with a macro similar to @lazycat?

The lovely package LazyArrays seems closest, but cannot currently do the above hvcating, as far a I can tell. I also tried BlockArrays, but that seems to be the complementary to LazyArrays: it can do the hvcating, but not the broadcasting…

1 Like