Adding a BlockDiagonal and Diagonal matrix fails for non-square matrices

[mzgubic]

Similar to #95 and #67 in the sense that the code is written for square blocks only.

julia> b = BlockDiagonal([ones(2, 3), ones(3, 2)])
5×5 BlockDiagonal{Float64, Matrix{Float64}}:
 1.0  1.0  1.0  0.0  0.0
 1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  1.0
 0.0  0.0  0.0  1.0  1.0
 0.0  0.0  0.0  1.0  1.0

julia> b + Diagonal(rand(5))
ERROR: BoundsError: attempt to access 3×2 Matrix{Float64} at index [3, 3]
Stacktrace:
 [1] getindex
   @ ./array.jl:862 [inlined]
 [2] +(B::BlockDiagonal{Float64, Matrix{Float64}}, M::Diagonal{Float64, Vector{Float64}})
   @ BlockDiagonals ~/JuliaEnvs/BlockDiagonals.jl/src/base_maths.jl:54
 [3] top-level scope
   @ REPL[8]:1

This would require us to decide what to return (BlockDiagonal vs Matrix) based on the shapes of the arrays, which would be a type instability.