Should we special-case block diagonal matrices?

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.

Start with the benchmark in the issue, comparing cholesky($(Diagonal(B))) with cholesky($B) for the block-diagonal mortar matrix. Investigate the relevant cholesky and block-structure entry points; done would require a decided approach for detecting block-diagonal matrices and evidence that the optimized factorization path improves the shown case without regressing other factorizations.