From 06be3bd8bd7c7b1df2bd313c62fef422b1e67f91 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 28 Apr 2025 12:43:55 +0530 Subject: [PATCH 1/4] Copy matrices in `triu`/`tril` if no zero exists --- src/dense.jl | 11 +++++++++++ src/diagonal.jl | 4 ++-- src/generic.jl | 10 ++++++++++ test/diagonal.jl | 7 +++++++ 4 files changed, 30 insertions(+), 2 deletions(-) diff --git a/src/dense.jl b/src/dense.jl index 9e7eff4f..49533151 100644 --- a/src/dense.jl +++ b/src/dense.jl @@ -124,6 +124,17 @@ haszero(::Type) = false haszero(::Type{T}) where {T<:Number} = isconcretetype(T) haszero(::Type{Union{Missing,T}}) where {T<:Number} = haszero(T) @propagate_inbounds _zero(M::AbstractArray{T}, inds...) where {T} = haszero(T) ? zero(T) : zero(M[inds...]) +function zero!(M::AbstractArray{T}) where {T} + if haszero(T) + fill!(M, zero(T)) + else + for i in eachindex(M) + v = @inbounds M[i] + z = zero(v) + @inbounds M[i] = z + end + end +end """ triu!(M, k::Integer) diff --git a/src/diagonal.jl b/src/diagonal.jl index 3f72c16c..0d9d469f 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -279,7 +279,7 @@ function triu!(D::Diagonal{T}, k::Integer=0) where T throw(ArgumentError(string("the requested diagonal, $k, must be at least ", "$(-n + 1) and at most $(n + 1) in an $n-by-$n matrix"))) elseif k > 0 - fill!(D.diag, zero(T)) + zero!(D.diag) end return D end @@ -290,7 +290,7 @@ function tril!(D::Diagonal{T}, k::Integer=0) where T throw(ArgumentError(LazyString(lazy"the requested diagonal, $k, must be at least ", lazy"$(-n - 1) and at most $(n - 1) in an $n-by-$n matrix"))) elseif k < 0 - fill!(D.diag, zero(T)) + zero!(D.diag) end return D end diff --git a/src/generic.jl b/src/generic.jl index 21bcd4cb..cde15157 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -466,6 +466,11 @@ julia> triu(a,-3) """ function triu(M::AbstractMatrix, k::Integer = 0) d = similar(M) + if !haszero(eltype(M)) + # since the zero would need to be evaluated from the elements, + # we copy the array to avoid undefined references in triu! + copy!(d, M) + end A = triu!(d,k) if iszero(k) copytrito!(A, M, 'U') @@ -509,6 +514,11 @@ julia> tril(a,-3) """ function tril(M::AbstractMatrix,k::Integer=0) d = similar(M) + if !haszero(eltype(M)) + # since the zero would need to be evaluated from the elements, + # we copy the array to avoid undefined references in tril! + copy!(d, M) + end A = tril!(d,k) if iszero(k) copytrito!(A, M, 'L') diff --git a/test/diagonal.jl b/test/diagonal.jl index 2b7cfd4d..ca9a2141 100644 --- a/test/diagonal.jl +++ b/test/diagonal.jl @@ -819,6 +819,13 @@ end @test all(x -> size(x) == (2,2), D) @test D == D1 * D2 end + + @testset "triu/tril" begin + D = Diagonal(fill(ones(2,2), 3)) + M = Matrix{eltype(D)}(D) + @test triu(D,1) == triu(M,1) + @test tril(D,-1) == tril(M,-1) + end end @testset "Eigensystem for block diagonal (issue #30681)" begin From 2fd0a12edf09445c1d3a64036e76b661f3385328 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 28 Apr 2025 13:08:47 +0530 Subject: [PATCH 2/4] Remove diagonal changes --- src/dense.jl | 11 ----------- src/diagonal.jl | 4 ++-- test/dense.jl | 9 +++++++++ test/diagonal.jl | 7 ------- 4 files changed, 11 insertions(+), 20 deletions(-) diff --git a/src/dense.jl b/src/dense.jl index 49533151..9e7eff4f 100644 --- a/src/dense.jl +++ b/src/dense.jl @@ -124,17 +124,6 @@ haszero(::Type) = false haszero(::Type{T}) where {T<:Number} = isconcretetype(T) haszero(::Type{Union{Missing,T}}) where {T<:Number} = haszero(T) @propagate_inbounds _zero(M::AbstractArray{T}, inds...) where {T} = haszero(T) ? zero(T) : zero(M[inds...]) -function zero!(M::AbstractArray{T}) where {T} - if haszero(T) - fill!(M, zero(T)) - else - for i in eachindex(M) - v = @inbounds M[i] - z = zero(v) - @inbounds M[i] = z - end - end -end """ triu!(M, k::Integer) diff --git a/src/diagonal.jl b/src/diagonal.jl index 0d9d469f..3f72c16c 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -279,7 +279,7 @@ function triu!(D::Diagonal{T}, k::Integer=0) where T throw(ArgumentError(string("the requested diagonal, $k, must be at least ", "$(-n + 1) and at most $(n + 1) in an $n-by-$n matrix"))) elseif k > 0 - zero!(D.diag) + fill!(D.diag, zero(T)) end return D end @@ -290,7 +290,7 @@ function tril!(D::Diagonal{T}, k::Integer=0) where T throw(ArgumentError(LazyString(lazy"the requested diagonal, $k, must be at least ", lazy"$(-n - 1) and at most $(n - 1) in an $n-by-$n matrix"))) elseif k < 0 - zero!(D.diag) + fill!(D.diag, zero(T)) end return D end diff --git a/test/dense.jl b/test/dense.jl index ee659748..ef41980c 100644 --- a/test/dense.jl +++ b/test/dense.jl @@ -1408,4 +1408,13 @@ end @test 2^A == 2^Matrix(A) end +@testset "triu/tril for block matrices" begin + O = ones(2,2) + Z = zero(O) + M = fill(O, 3, 3) + res = fill(Z, size(M)) + res[1,2] = res[1,3] = res[2,3] = O + @test triu(GenericArray(M),1) == res +end + end # module TestDense diff --git a/test/diagonal.jl b/test/diagonal.jl index ca9a2141..2b7cfd4d 100644 --- a/test/diagonal.jl +++ b/test/diagonal.jl @@ -819,13 +819,6 @@ end @test all(x -> size(x) == (2,2), D) @test D == D1 * D2 end - - @testset "triu/tril" begin - D = Diagonal(fill(ones(2,2), 3)) - M = Matrix{eltype(D)}(D) - @test triu(D,1) == triu(M,1) - @test tril(D,-1) == tril(M,-1) - end end @testset "Eigensystem for block diagonal (issue #30681)" begin From b9a0dfb010a6338ebc0080d02652564cb4709c5d Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 28 Apr 2025 15:55:59 +0530 Subject: [PATCH 3/4] Add tril test --- test/dense.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/dense.jl b/test/dense.jl index ef41980c..99e0b429 100644 --- a/test/dense.jl +++ b/test/dense.jl @@ -1415,6 +1415,7 @@ end res = fill(Z, size(M)) res[1,2] = res[1,3] = res[2,3] = O @test triu(GenericArray(M),1) == res + @test tril(GenericArray(M),-1) == res' end end # module TestDense From 8154f8e03c6c6629d9dee361f6a6fe3c2d427dab Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Mon, 28 Apr 2025 16:09:10 +0530 Subject: [PATCH 4/4] dispatch on haszero --- src/generic.jl | 32 ++++++++++++++++++++------------ 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/src/generic.jl b/src/generic.jl index cde15157..afbd7b9f 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -464,13 +464,9 @@ julia> triu(a,-3) 1.0 1.0 1.0 1.0 ``` """ -function triu(M::AbstractMatrix, k::Integer = 0) +triu(M::AbstractMatrix, k::Integer = 0) = _triu(M, Val(haszero(eltype(M))), k) +function _triu(M::AbstractMatrix, ::Val{true}, k::Integer) d = similar(M) - if !haszero(eltype(M)) - # since the zero would need to be evaluated from the elements, - # we copy the array to avoid undefined references in triu! - copy!(d, M) - end A = triu!(d,k) if iszero(k) copytrito!(A, M, 'U') @@ -482,6 +478,14 @@ function triu(M::AbstractMatrix, k::Integer = 0) end return A end +function _triu(M::AbstractMatrix, ::Val{false}, k::Integer) + d = similar(M) + # since the zero would need to be evaluated from the elements, + # we copy the array to avoid undefined references in triu! + copy!(d, M) + A = triu!(d,k) + return A +end """ tril(M, k::Integer = 0) @@ -512,13 +516,9 @@ julia> tril(a,-3) 1.0 0.0 0.0 0.0 ``` """ -function tril(M::AbstractMatrix,k::Integer=0) +tril(M::AbstractMatrix,k::Integer=0) = _tril(M, Val(haszero(eltype(M))), k) +function _tril(M::AbstractMatrix, ::Val{true}, k::Integer) d = similar(M) - if !haszero(eltype(M)) - # since the zero would need to be evaluated from the elements, - # we copy the array to avoid undefined references in tril! - copy!(d, M) - end A = tril!(d,k) if iszero(k) copytrito!(A, M, 'L') @@ -530,6 +530,14 @@ function tril(M::AbstractMatrix,k::Integer=0) end return A end +function _tril(M::AbstractMatrix, ::Val{false}, k::Integer) + d = similar(M) + # since the zero would need to be evaluated from the elements, + # we copy the array to avoid undefined references in tril! + copy!(d, M) + A = tril!(d,k) + return A +end """ triu!(M)