Don't access uninitialized indices in tril!/triu! for numeric arrays

[jishnub]

This specializes the generic triu! and tril! methods for arrays of numbers, where a zero is known to exist. In such cases, we don't need to read the possibly uninitialized values of the elements at the indices that are to be assigned to. After this, the following would be possible:

julia> M = Matrix{BigFloat}(undef, 3, 3)
3×3 Matrix{BigFloat}:
 #undef  #undef  #undef
 #undef  #undef  #undef
 #undef  #undef  #undef

julia> triu!(M)
3×3 Matrix{BigFloat}:
 #undef  #undef  #undef
   0.0   #undef  #undef
   0.0     0.0   #undef

[mikmoore]

It seems like the current implementation is insufficiently general. I expect it to fail sometimes on matrices with Union types. A Union can be <:Number but not have an applicable zero. Observe the following:

julia> x = fill!(Matrix{Union{Float64,ComplexF64}}(undef,3,3), 0.0)
3×3 Matrix{Union{Float64, ComplexF64}}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> x isa AbstractMatrix{<:Number}
true

julia> x[1,1] = zero(eltype(x))
ERROR: MethodError: no method matching Union{Float64, ComplexF64}(::Int64)

The above matrix would fit the signature of the provided specialization but would then throw an error because zero(Union{Float64,ComplexF64}) is not defined.

Would it work to assign each entry by checking whether isassigned, then use either zero on that entry (if assigned) or zero on the eltype if not? I assume the isassigned would be compiled away for arrays that cannot hold #undef elements? Or you could try the reverse direction of using zero(eltype()) when valid and only checking the specific indices when not?

Not sure what the best resolution is.

[jishnub]

Interesting observation! This seems to also impact other structured matrix types:

julia> Diagonal(Union{Float64,ComplexF64}[1,2])
ERROR: MethodError: no method matching Union{Float64, ComplexF64}(::Int64)

where the off-diagonal zeros are not defined. In that case, there isn't a parent to probe, either.

Still, this PR appears to break this case, as the following used to work until now:

julia> UpperTriangular(Matrix{Union{Float64,ComplexF64}}(undef, 2, 2))
2×2 UpperTriangular{Union{Float64, ComplexF64}, Matrix{Union{Float64, ComplexF64}}}:
 0.0+0.0im  0.0+0.0im
     ⋅      0.0+0.0im

[mikmoore]

The Diagonal issue isn't what you show (your error was in constructing the underlying Vector, since it doesn't know whether 1 should be converted to Float64 or ComplexF64). There is no issue with construction with that resolved. However, the same error arises from a subsequent show (like in the REPL's end-of-input display):

julia> Diagonal(Union{Float64,ComplexF64}[1.0,2.0])
2×2 Diagonal{Union{Float64, ComplexF64}, Vector{Union{Float64, ComplexF64}}}:
Error showing value of type Diagonal{Union{Float64, ComplexF64}, Vector{Union{Float64, ComplexF64}}}:
ERROR: MethodError: no method matching Union{Float64, ComplexF64}(::Int64)

julia> Diagonal([1.0,2.0]) # why did it even need to know what `zero(eltype(...))` is? it's not shown
2×2 Diagonal{Float64, Vector{Float64}}:
 1.0   ⋅
  ⋅   2.0

Of course, a matrix with structural zeros requires a valid zero to sample them, so these are fragile matrices to work with in any case. I'm sure lots of other operations break too. So perhaps this shouldn't be held up on that account? Hard to say.

From the stacktrace, it appears that show is sampling the off-diagonals. It's unclear why this is necessary, since the structural zeros are printed as dots.

[jishnub]

Perhaps we should use zero on the eltype for concrete number types? This should work around the issue with Unions.

[jishnub]

Gentle bump

[dkarrasch]

I wonder if the haszero trait should be declared public, so that other packages such as StaticArrays could indicate their types as suitable. But we could think about that later, I guess.

Otherwise LGTM.

[jishnub]

@mikmoore I wonder if you foresee any issues with the current implementation?

[mikmoore]

Looks good to me. The last "arguably useful" holdout I foresee is union types.

One might attempt Base.promote_union on a type in haszero and _zero. This would allow this to resolve types like Union{Float64, ComplexF64} even if an entry is #undef. But it actually appears that concrete Unions can't be #undef (?), so this may be irrelevant...

One could go further and remove Missing and Nothing from union types too, in performing the assessment. Or, as a more generic alternative, one could sequentially check every type in the a union and simply use the first one that haszero (erroring if none do).

I could see the Missing and maybe Nothing cases making a difference, but none of these seem likely to make a big difference in practice. The PR looks reasonable as it is and exploring my ideas here can always be done in the future if the desire arises.