Specialise `tr` when matrix is square, but blocks not all square

[mjp98]

Currently tr is computed block by block, which works when all blocks are square. For consistency with LinearAlgebra this should probably work as long as the total matrix is square.

julia> using BlockDiagonals

julia> B = BlockDiagonal([randn(1,2),rand(2,1)])
3×3 BlockDiagonal{Float64, Matrix{Float64}}:
 -1.2767  -1.60808  0.0
  0.0      0.0      0.368744
  0.0      0.0      0.381816

julia> tr(Matrix(B))
-0.8948857867073675

julia> tr(B)
ERROR: DimensionMismatch: matrix is not square: dimensions are (1, 2)
Stacktrace:
  [1] checksquare
    @ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/LinearAlgebra.jl:234 [inlined]
  [2] tr(A::Matrix{Float64})
    @ LinearAlgebra /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/dense.jl:345
  [3] _mapreduce(f::typeof(tr), op::typeof(Base.add_sum), #unused#::IndexLinear, A::Vector{Matrix{Float64}})
    @ Base ./reduce.jl:435
  [4] _mapreduce_dim
    @ ./reducedim.jl:365 [inlined]
  [5] #mapreduce#765
    @ ./reducedim.jl:357 [inlined]
  [6] mapreduce
    @ ./reducedim.jl:357 [inlined]
  [7] #_sum#775
    @ ./reducedim.jl:999 [inlined]
  [8] _sum
    @ ./reducedim.jl:999 [inlined]
  [9] #sum#773
    @ ./reducedim.jl:995 [inlined]
 [10] sum
    @ ./reducedim.jl:995 [inlined]
 [11] tr(B::BlockDiagonal{Float64, Matrix{Float64}})
    @ BlockDiagonals ~/.julia/packages/BlockDiagonals/gqKhZ/src/linalg.jl:10
 [12] top-level scope
    @ REPL[180]:1