getindex execution time is linear in the number of blocks

[mzgubic]

using LinearAlgebra
using BlockDiagonals

blocksizes = [2, 4, 8, 16, 32, 64, 128, 256, 512, 1024]
btimes = []
mtimes = []

for nblocks in blocksizes
    b = BlockDiagonal([rand(1, 1) for _ in 1:nblocks])
    m = rand(nblocks, nblocks)

    t = @benchmark $b[$nblocks, $nblocks]
    push!(btimes, mean(t).time)

    t = @benchmark $m[$nblocks, $nblocks]
    push!(mtimes, mean(t).time)
end

plot(blocksizes, btimes; label="BlockDiagonal", xlabel="matrix size", ylabel="time (ns)", title="execution time for getindex()")
plot!(blocksizes, mtimes; label= "Matrix")

This scaling is a problem when BlockDiagonals hit fallback implementations for matrices which often use getindex. Example in the wild #116.

The reason for this is the implementation of getindex which iterates over the blocks to find the right location - the sizes of blocks are not known in advance.

One way to improve this relationship would be to store the cumulative block sizes in both dimensions, and then do a bisection search, which would bring the relationship to O(LogN)

[jishnub]

BlockArrays does something like this, where the array axis stores the block sizes. I wonder if it makes sense for BlockDiagonal to subtype an AbstractBlockArray?