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

…ays (#52528)

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
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
```

stdlib/LinearAlgebra/src/dense.jl

@@ -106,6 +106,11 @@ norm1(x::Union{Array{T},StridedVector{T}}) where {T<:BlasReal} =
 norm2(x::Union{Array{T},StridedVector{T}}) where {T<:BlasFloat} =
     length(x) < NRM2_CUTOFF ? generic_norm2(x) : BLAS.nrm2(x)
 
+# Conservative assessment of types that have zero(T) defined for themselves
+haszero(::Type) = false
+haszero(::Type{T}) where {T<:Number} = isconcretetype(T)
+@propagate_inbounds _zero(M::AbstractArray{T}, i, j) where {T} = haszero(T) ? zero(T) : zero(M[i,j])
+
 """
     triu!(M, k::Integer)
 
@@ -136,7 +141,7 @@ function triu!(M::AbstractMatrix, k::Integer)
     m, n = size(M)
     for j in 1:min(n, m + k)
         for i in max(1, j - k + 1):m
-            @inbounds M[i,j] = zero(M[i,j])
+            @inbounds M[i,j] = _zero(M, i,j)
         end
     end
     M
@@ -173,12 +178,13 @@ function tril!(M::AbstractMatrix, k::Integer)
     require_one_based_indexing(M)
     m, n = size(M)
     for j in max(1, k + 1):n
-        @inbounds for i in 1:min(j - k - 1, m)
-            M[i,j] = zero(M[i,j])
+        for i in 1:min(j - k - 1, m)
+            @inbounds M[i,j] = _zero(M, i,j)
         end
     end
     M
 end
+
 tril(M::Matrix, k::Integer) = tril!(copy(M), k)
 
 """

stdlib/LinearAlgebra/test/triangular.jl

@@ -897,6 +897,35 @@ end
     end
 end
 
+@testset "tril!/triu! for non-bitstype matrices" begin
+    @testset "numeric" begin
+        M = Matrix{BigFloat}(undef, 3, 3)
+        tril!(M)
+        L = LowerTriangular(ones(3,3))
+        copytrito!(M, L, 'L')
+        @test M == L
+
+        M = Matrix{BigFloat}(undef, 3, 3)
+        triu!(M)
+        U = UpperTriangular(ones(3,3))
+        copytrito!(M, U, 'U')
+        @test M == U
+    end
+    @testset "array elements" begin
+        M = fill(ones(2,2), 4, 4)
+        tril!(M)
+        L = LowerTriangular(fill(fill(2,2,2),4,4))
+        copytrito!(M, L, 'L')
+        @test M == L
+
+        M = fill(ones(2,2), 4, 4)
+        triu!(M)
+        U = UpperTriangular(fill(fill(2,2,2),4,4))
+        copytrito!(M, U, 'U')
+        @test M == U
+    end
+end
+
 @testset "avoid matmul ambiguities with ::MyMatrix * ::AbstractMatrix" begin
     A = [i+j for i in 1:2, j in 1:2]
     S = SizedArrays.SizedArray{(2,2)}(A)