LinearAlgebra: diagzero for non-OneTo axes

[jishnub]

Currently, the off-diagonal zeros for a block-Diagonal matrix is computed using diagzero, which calls zeros for the sizes of the elements. This returns an Array, unless one specializes diagzero for the custom Diagonal matrix type.

This PR defines a zeroslike function that dispatches on the axes of the elements, which lets packages specialize on the axes to return custom AbstractArrays. Choosing to specialize on the eltype avoids the need to specialize on the container, and allows packages to return appropriate types for custom axis types.

With this,

julia> LinearAlgebra.zeroslike(::Type{S}, ax::Tuple{SOneTo, Vararg{SOneTo}}) where {S<:SMatrix} = SMatrix{map(length, ax)...}(ntuple(_->zero(eltype(S)), prod(length, ax)))

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

julia> D[1,2] # now an SMatrix
2×3 SMatrix{2, 3, Int64, 6} with indices SOneTo(2)×SOneTo(3):
 0  0  0
 0  0  0

julia> LinearAlgebra.zeroslike(::Type{S}, ax::Tuple{SOneTo, Vararg{SOneTo}}) where {S<:MMatrix} = MMatrix{map(length, ax)...}(ntuple(_->zero(eltype(S)), prod(length, ax)))

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

julia> D[1,2] # now an MMatrix
2×3 MMatrix{2, 3, Int64, 6} with indices SOneTo(2)×SOneTo(3):
 0  0  0
 0  0  0

The reason this can't be the default behavior is that we are not guaranteed that there exists a similar method that accepts the combination of axes. This is why we have to fall back to using the sizes, unless a specialized method is provided by a package.

One positive outcome of this is that indexing into such a block-diagonal matrix will now usually be type-stable, which mitigates JuliaLang/LinearAlgebra.jl#927 to some extent (although it doesn't resolve the issue).

I've also updated the getindex for Bidiagonal to use diagzero, instead of the similarly defined bidiagzero function that it was using. Structured block matrices may now use diagzero uniformly to generate the zero elements.

[mikmoore]

Can the LinearAlgebra.AbstractTriangular family also benefit from this PR? In fact, they appear to not even be correct in some cases:

julia> UpperTriangular([ones(i,j) for i in 1:2, j in 1:2])
2×2 UpperTriangular{Matrix{Float64}, Matrix{Matrix{Float64}}}:
 [1.0;;]    [1.0 1.0]
     ⋅      [1.0 1.0; 1.0 1.0]

julia> UpperTriangular([ones(i,j) for i in 1:2, j in 1:2])[1,2]
1×2 Matrix{Float64}:
 1.0  1.0

julia> UpperTriangular([ones(i,j) for i in 1:2, j in 1:2])[2,1] # output has wrong size!
1×2 Matrix{Float64}:
 0.0  0.0

[jishnub]

Sounds reasonable to me, but that may be a separate PR

[dkarrasch]

LGTM. Should this mechanism be communicated to package maintainers in NEWS.md?