Support broadcasting over structured block matrices

[jishnub]

Fix JuliaLang/LinearAlgebra.jl#985

After this, broadcasting over structured block matrices with matrix-valued elements works:

julia> D = Diagonal([[1 2; 3 4], [5 6; 7 8]])
2×2 Diagonal{Matrix{Int64}, Vector{Matrix{Int64}}}:
 [1 2; 3 4]      ⋅     
     ⋅       [5 6; 7 8]

julia> D .+ D
2×2 Diagonal{Matrix{Int64}, Vector{Matrix{Int64}}}:
 [2 4; 6 8]      ⋅     
     ⋅       [10 12; 14 16]

julia> cos.(D)
2×2 Matrix{Matrix{Float64}}:
 [0.855423 -0.110876; -0.166315 0.689109]  [1.0 0.0; 0.0 1.0]
 [1.0 0.0; 0.0 1.0]                        [0.928384 -0.069963; -0.0816235 0.893403]

Such operations show up when using BlockArrays.

The implementation is a bit hacky as it uses 0I as the zero element in fzero, which isn't really the correct zero if the blocks are rectangular. Nonetheless, this works, as fzero is only used to determine if the structure is preserved.

[dkarrasch]

Awesome! Should we perhaps have some failing tests, like matrices of matrices whose elements don't match in size, or something?

[jishnub]

I'm actually struggling to come up with an example that specifically fails despite this, and isn't related to the properties of the elements (or missing UniformScaling methods). Perhaps I'll add some tests that are broken due to missing methods? These would hopefully be fixed in the future.

[dkarrasch]

I was thinking of something like

A = Diagonal([zeros(2,3) for _ in 1:3]) # or zeros(3,3)
B = Diagonal([zeros(2,2) for _ in 1:3])
@testthrows SomeError A .+ B

[jishnub]

I've added a test with incompatible element sizes, and also one that is directly related to this PR's changes:

julia> D = Diagonal([ones(1,1), ones(1,1)])
2×2 Diagonal{Matrix{Float64}, Vector{Matrix{Float64}}}:
 [1.0;;]     ⋅   
    ⋅     [1.0;;]

julia> size.(D)
ERROR: MethodError: no method matching size(::UniformScaling{Float64})
The function `size` exists, but no method is defined for this combination of argument types.
[...]

Thankfully, this isn't a regression, as the operation didn't work before either. In fact, the following works now:

julia> ndims.(D)
2×2 Matrix{Int64}:
 2  2
 2  2