From edd935f0bd9eaec163547baabb155c981dc5eab3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Thu, 18 Mar 2021 23:11:16 +0100 Subject: [PATCH 01/10] first commit --- base/abstractarray.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/base/abstractarray.jl b/base/abstractarray.jl index 50333a7144871..79a807da7eef7 100644 --- a/base/abstractarray.jl +++ b/base/abstractarray.jl @@ -323,6 +323,16 @@ function _all_match_first(f::F, inds, A, B...) where F<:Function end _all_match_first(f::F, inds) where F<:Function = true +""" + eachstoredindex(A) + +Returns an iterable over the indices of `A` where the values are structurally non-zero. +It falls back to `eachindex(A)` and can be redefined by array types. +`eachstoredindex(A)` is not guaranteed to return the same shape as `eachindex(A)`. +""" +eachstoredindex(A) = eachindex(A) + + # keys with an IndexStyle keys(s::IndexStyle, A::AbstractArray, B::AbstractArray...) = eachindex(s, A, B...) From dbea613bea3e2f4c787e1261ae3569a6b8c7bd35 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 19 Mar 2021 00:07:24 +0100 Subject: [PATCH 02/10] diagonal --- stdlib/LinearAlgebra/src/diagonal.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 10e8e64e01ae7..52a71dc683269 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -111,6 +111,11 @@ end parent(D::Diagonal) = D.diag +function Base.eachstoredindex(D::Diagonal) + n = length(D.diag) + return 1:n+1:(n*n) +end + ishermitian(D::Diagonal{<:Real}) = true ishermitian(D::Diagonal{<:Number}) = isreal(D.diag) ishermitian(D::Diagonal) = all(ishermitian, D.diag) From 0d2d90d7d2edb40305bd9b96f3b494337179ce16 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Sun, 21 Mar 2021 15:26:14 +0100 Subject: [PATCH 03/10] Update stdlib/LinearAlgebra/src/diagonal.jl Co-authored-by: Daniel Karrasch --- stdlib/LinearAlgebra/src/diagonal.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 52a71dc683269..6fd1fb0e9820e 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -111,10 +111,7 @@ end parent(D::Diagonal) = D.diag -function Base.eachstoredindex(D::Diagonal) - n = length(D.diag) - return 1:n+1:(n*n) -end +Base.eachstoredindex(D::Diagonal) = diagind(D) ishermitian(D::Diagonal{<:Real}) = true ishermitian(D::Diagonal{<:Number}) = isreal(D.diag) From deaab22b44681f99f5dc6fb624174f839184b454 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Sun, 21 Mar 2021 16:18:45 +0100 Subject: [PATCH 04/10] implement eachstoredindex for sparse vecs mats --- stdlib/SparseArrays/Project.toml | 4 +++- stdlib/SparseArrays/src/sparsematrix.jl | 11 +++++++++++ stdlib/SparseArrays/src/sparsevector.jl | 3 +++ stdlib/SparseArrays/test/sparse.jl | 2 ++ stdlib/SparseArrays/test/sparsevector.jl | 2 ++ 5 files changed, 21 insertions(+), 1 deletion(-) diff --git a/stdlib/SparseArrays/Project.toml b/stdlib/SparseArrays/Project.toml index 53d4a9f064ad3..31cf5b2156899 100644 --- a/stdlib/SparseArrays/Project.toml +++ b/stdlib/SparseArrays/Project.toml @@ -3,12 +3,14 @@ uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [extras] Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Dates", "Test", "InteractiveUtils"] +test = ["Dates", "Test", "InteractiveUtils", "Printf"] diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index 7bd988d881152..5bac41d1b7d3e 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -1525,6 +1525,17 @@ function findnz(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} return (I, J, V) end +function Base.eachstoredindex(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} + numnz = nnz(S) + indices = Vector{CartesianIndex{2}}(undef, numnz) + count = 1 + @inbounds for col = 1 : size(S, 2), k = getcolptr(S)[col] : (getcolptr(S)[col+1]-1) + indices[count] = CartesianIndex(rowvals(S)[k], col) + count += 1 + end + return indices +end + function _sparse_findnextnz(m::AbstractSparseMatrixCSC, ij::CartesianIndex{2}) row, col = Tuple(ij) col > size(m, 2) && return nothing diff --git a/stdlib/SparseArrays/src/sparsevector.jl b/stdlib/SparseArrays/src/sparsevector.jl index 777be897ea7df..b3f33089b1fda 100644 --- a/stdlib/SparseArrays/src/sparsevector.jl +++ b/stdlib/SparseArrays/src/sparsevector.jl @@ -77,6 +77,9 @@ function nonzeroinds(x::SparseColumnView) @inbounds y = view(rowvals(A), nzrange(A, colidx)) return y end + +Base.eachstoredindex(x::SparseVector) = getfield(x, :nzind) + nonzeroinds(x::SparseVectorView) = nonzeroinds(parent(x)) rowvals(x::SparseVectorUnion) = nonzeroinds(x) diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index ef5159d41031d..b5565b532eb3d 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -1833,7 +1833,9 @@ end # test that stored zeros are still stored zeros in the diagonal S = sparse([1,3],[1,3],[0.0,0.0]); V = diag(S) @test nonzeroinds(V) == [1,3] + @test Base.eachstoredindex(V) == [1,3] @test nonzeros(V) == [0.0,0.0] + @test V[Base.eachstoredindex(V)] == nonzeros(V) end @testset "expandptr" begin diff --git a/stdlib/SparseArrays/test/sparsevector.jl b/stdlib/SparseArrays/test/sparsevector.jl index 602818678908b..87d28a1bec2ee 100644 --- a/stdlib/SparseArrays/test/sparsevector.jl +++ b/stdlib/SparseArrays/test/sparsevector.jl @@ -31,7 +31,9 @@ x1_full[SparseArrays.nonzeroinds(spv_x1)] = nonzeros(spv_x1) @test count(!iszero, x) == 3 @test nnz(x) == 3 @test SparseArrays.nonzeroinds(x) == [2, 5, 6] + @test Base.eachstoredindex(x) == [2, 5, 6] @test nonzeros(x) == [1.25, -0.75, 3.5] + @test nonzeros(x) == x[Base.eachstoredindex(x)] @test count(SparseVector(8, [2, 5, 6], [true,false,true])) == 2 end From 5b3609a4febcf85397bd5228ec0869cb366b234d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Sun, 21 Mar 2021 16:56:34 +0100 Subject: [PATCH 05/10] implemented for wrappers on sparse matrices --- stdlib/SparseArrays/src/sparsematrix.jl | 15 +++++++++++++++ stdlib/SparseArrays/test/sparse.jl | 8 ++++++++ 2 files changed, 23 insertions(+) diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index 5bac41d1b7d3e..ce17a0b50713a 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -1536,6 +1536,21 @@ function Base.eachstoredindex(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} return indices end +function Base.eachstoredindex(S::Union{Adjoint{T, ST}, Transpose{T, ST}}) where {T, ST <: AbstractSparseMatrixCSC{T}} + numnz = nnz(S) + indices = Vector{CartesianIndex{2}}(undef, numnz) + count = 1 + @inbounds for col = 1 : size(S, 2), k = getcolptr(S)[col] : (getcolptr(S)[col+1]-1) + indices[count] = CartesianIndex(col, rowvals(S)[k]) + count += 1 + end + return indices +end + +function Base.eachstoredindex(S::Union{Symmetric{T, ST}, LinearAlgebra.AbstractTriangular{T, ST}, UpperHessenberg{T, ST}, Hermitian{T, ST}}) where {T, ST <: AbstractSparseMatrixCSC{T}} + return Base.eachstoredindex(parent(S)) +end + function _sparse_findnextnz(m::AbstractSparseMatrixCSC, ij::CartesianIndex{2}) row, col = Tuple(ij) col > size(m, 2) && return nothing diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index b5565b532eb3d..e23f4c40e3b76 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -2877,6 +2877,14 @@ end @test SparseMatrixCSC(at(wr(A))) == Matrix(at(wr(B))) end + @testset "eachstoredindex($(wr))" for wr in (UpperTriangular, LowerTriangular, + UnitUpperTriangular, UnitLowerTriangular, + Hermitian, (Hermitian, :L), Symmetric, (Symmetric, :L), Transpose, Adjoint) + S = dowrap(wr, A) + sum(S[Base.eachstoredindex(S)]) == sum(S) + sum(S[Base.eachstoredindex(S)]) == sum(Matrix(S)) + end + @test sparse([1,2,3,4,5]') == SparseMatrixCSC([1 2 3 4 5]) @test sparse(UpperTriangular(A')) == UpperTriangular(B') @test sparse(Adjoint(UpperTriangular(A'))) == Adjoint(UpperTriangular(B')) From da28eb7392e05794103c0c0b4d4117eb7c2da02e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Sun, 21 Mar 2021 17:55:46 +0100 Subject: [PATCH 06/10] call parent array in adjoint/transpose --- stdlib/SparseArrays/src/sparsematrix.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index ce17a0b50713a..60cb0e56b1f7e 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -1537,11 +1537,12 @@ function Base.eachstoredindex(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} end function Base.eachstoredindex(S::Union{Adjoint{T, ST}, Transpose{T, ST}}) where {T, ST <: AbstractSparseMatrixCSC{T}} - numnz = nnz(S) + P = parent(S) + numnz = nnz(P) indices = Vector{CartesianIndex{2}}(undef, numnz) count = 1 - @inbounds for col = 1 : size(S, 2), k = getcolptr(S)[col] : (getcolptr(S)[col+1]-1) - indices[count] = CartesianIndex(col, rowvals(S)[k]) + @inbounds for col = 1 : size(P, 2), k = getcolptr(P)[col] : (getcolptr(P)[col+1]-1) + indices[count] = CartesianIndex(col, rowvals(P)[k]) count += 1 end return indices From 76b4f2577d74041d228f16914fad5af56dd83145 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Sun, 21 Mar 2021 22:03:13 +0100 Subject: [PATCH 07/10] remove deps on Printf in SparseArrays --- stdlib/SparseArrays/Project.toml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/stdlib/SparseArrays/Project.toml b/stdlib/SparseArrays/Project.toml index 31cf5b2156899..53d4a9f064ad3 100644 --- a/stdlib/SparseArrays/Project.toml +++ b/stdlib/SparseArrays/Project.toml @@ -3,14 +3,12 @@ uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [extras] Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Dates", "Test", "InteractiveUtils", "Printf"] +test = ["Dates", "Test", "InteractiveUtils"] From 3ec4c865c99d0ff4f777eed22261431357bdead2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Wed, 24 Mar 2021 15:40:11 +0100 Subject: [PATCH 08/10] Update stdlib/SparseArrays/src/sparsematrix.jl Co-authored-by: Daniel Karrasch --- stdlib/SparseArrays/src/sparsematrix.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index 60cb0e56b1f7e..f408d07e2f76b 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -1529,7 +1529,7 @@ function Base.eachstoredindex(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} numnz = nnz(S) indices = Vector{CartesianIndex{2}}(undef, numnz) count = 1 - @inbounds for col = 1 : size(S, 2), k = getcolptr(S)[col] : (getcolptr(S)[col+1]-1) + @inbounds for col = 1 : size(S, 2), k in nzrange(S, col) indices[count] = CartesianIndex(rowvals(S)[k], col) count += 1 end From bc6d944eff9e415c01944a355489ab1864808c73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Wed, 24 Mar 2021 15:42:48 +0100 Subject: [PATCH 09/10] homogeneize col access --- stdlib/SparseArrays/src/sparsematrix.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index f408d07e2f76b..02ac5b50c878a 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -1515,7 +1515,7 @@ function findnz(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} V = Vector{Tv}(undef, numnz) count = 1 - @inbounds for col = 1 : size(S, 2), k = getcolptr(S)[col] : (getcolptr(S)[col+1]-1) + @inbounds for col in 1:size(S, 2), k in nzrange(S, col) I[count] = rowvals(S)[k] J[count] = col V[count] = nonzeros(S)[k] @@ -1529,7 +1529,7 @@ function Base.eachstoredindex(S::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti} numnz = nnz(S) indices = Vector{CartesianIndex{2}}(undef, numnz) count = 1 - @inbounds for col = 1 : size(S, 2), k in nzrange(S, col) + @inbounds for col in 1:size(S, 2), k in nzrange(S, col) indices[count] = CartesianIndex(rowvals(S)[k], col) count += 1 end From df097859603b276e3278c1a46f1afb9fb3f4f672 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Wed, 24 Mar 2021 15:56:04 +0100 Subject: [PATCH 10/10] add NEWS entry --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index e0885b713f7ae..d832c81e2cef2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -34,6 +34,7 @@ New library functions * `isunordered(x)` returns true if `x` is value that is normally unordered, such as `NaN` or `missing`. * New macro `Base.@invokelatest f(args...; kwargs...)` provides a convenient way to call `Base.invokelatest(f, args...; kwargs...)` ([#37971]) * New macro `Base.@invoke f(arg1::T1, arg2::T2; kwargs...)` provides an easier syntax to call `invoke(f, Tuple{T1,T2}, arg1, arg2; kwargs...)` ([#38438]) +* New function `Base.eachstoredindex` returns an index iterator over the structural non-zero indices of a container ([#40103]) New library features --------------------