Should we special-case block diagonal matrices?

[jishnub]

julia> n = 400;

julia> B = mortar(reshape([2I(n), 0I(n), 0I(n), 4I(n)], 2,2));

julia> @btime cholesky($(Diagonal(B)));
  2.298 μs (1 allocation: 6.38 KiB)

julia> @btime cholesky($B);
  4.235 ms (2 allocations: 4.88 MiB)

Since checking for block diagonal structure is order n^2 while cholesky and other factorizations are n^3, it might be worth carrying out this check and call more optimized methods.

[dlfivefifty]

My only argument against it is that it makes the code more complicated to maintain

[jishnub]

Perhaps it would make sense to have BlockDiagonalMatrix <: AbstactBlockMatrix, as such a structure is often known a-priori. Special methods may then be added for such a matrix type. I imagine it'll be some work to resolve method ambiguities, so perhaps this can be attempted in a separate package, and once it is stable, it may be merged with BlockArrays. In this approach, BlockArray methods don't need to infer matrix structure, so this won't complicate dispatch.

[dlfivefifty]

There are already two packages for this:

1. BlockBandedMatrices.jl
2. BlockDiagonals.jl