From 0c7863bf26a4600ff460afbfc113f1a0c8b2aab5 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Mon, 14 Aug 2017 22:47:56 +0200 Subject: [PATCH 1/6] diag of SparseMatrixCSC should always return SparseVector --- base/sparse/linalg.jl | 2 +- base/sparse/sparsematrix.jl | 2 -- test/sparse/sparse.jl | 10 ++++++++++ 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/base/sparse/linalg.jl b/base/sparse/linalg.jl index e1403d2967cb8..ad56b6ed057cc 100644 --- a/base/sparse/linalg.jl +++ b/base/sparse/linalg.jl @@ -879,7 +879,7 @@ for f in (:\, :Ac_ldiv_B, :At_ldiv_B) if m == n if istril(A) if istriu(A) - return ($f)(Diagonal(A), B) + return ($f)(Diagonal(Vector(diag(A))), B) else return ($f)(LowerTriangular(A), B) end diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index 212565d4cddcc..9195d5c837ef3 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -3412,8 +3412,6 @@ function trace(A::SparseMatrixCSC{Tv}) where Tv s end -diag(A::SparseMatrixCSC{Tv}) where {Tv} = Tv[d for d in SpDiagIterator(A)] - function diagm(v::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} if size(v,1) != 1 && size(v,2) != 1 throw(DimensionMismatch("input should be nx1 or 1xn")) diff --git a/test/sparse/sparse.jl b/test/sparse/sparse.jl index 714c8fea0f582..52eb705dfc2a9 100644 --- a/test/sparse/sparse.jl +++ b/test/sparse/sparse.jl @@ -1323,6 +1323,16 @@ end @test diagm(sparse(ones(5,1))) == speye(5) end +@testset "diag" begin + for T in (Float64, Complex128) + S = sprand(T, 5, 5, 0.5) + A = Matrix(S) + @test diag(S)::SparseVector{T,Int} == diag(A) + @test diag(S, 1)::SparseVector{T,Int} == diag(A, 1) + @test diag(S, -1)::SparseVector{T,Int} == diag(A, -1) + end +end + @testset "expandptr" begin A = speye(5) @test Base.SparseArrays.expandptr(A.colptr) == collect(1:5) From 7af9ab2289096fe9427a9ebfaef616e36c0bde88 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Wed, 16 Aug 2017 16:57:38 +0200 Subject: [PATCH 2/6] rewrite SpDiagIterator to accept any diagonal --- base/sparse/sparsematrix.jl | 43 +++++++++++++++++++++++++------------ test/sparse/sparse.jl | 17 ++++++++++----- 2 files changed, 41 insertions(+), 19 deletions(-) diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index 9195d5c837ef3..be4a55f5df93b 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -3382,23 +3382,38 @@ end ## diag and related using an iterator -mutable struct SpDiagIterator{Tv,Ti} +struct SpDiagIterator{Tv,Ti} A::SparseMatrixCSC{Tv,Ti} - n::Int + d::Int # diagonal to iterate over end -SpDiagIterator(A::SparseMatrixCSC) = SpDiagIterator(A,minimum(size(A))) -length(d::SpDiagIterator) = d.n -start(d::SpDiagIterator) = 1 -done(d::SpDiagIterator, j) = j > d.n +start(it::SpDiagIterator) = it.d < 0 ? (1-it.d, 1) : (1, it.d+1) +done(it::SpDiagIterator, rc) = rc[1] > size(it.A,1) || rc[2] > size(it.A,2) +function next(it::SpDiagIterator{Tv}, rc) where Tv + r, c = rc + A = it.A + r1 = Int(A.colptr[c]) + r2 = Int(A.colptr[c+1]-1) + (r1 > r2) && (return (zero(Tv), (r+1, c+1))) + r1 = searchsortedfirst(A.rowval, r, r1, r2, Forward) + (((r1 > r2) || (A.rowval[r1] != r)) ? zero(Tv) : A.nzval[r1], (r+1, c+1)) +end -function next(d::SpDiagIterator{Tv}, j) where Tv - A = d.A - r1 = Int(A.colptr[j]) - r2 = Int(A.colptr[j+1]-1) - (r1 > r2) && (return (zero(Tv), j+1)) - r1 = searchsortedfirst(A.rowval, j, r1, r2, Forward) - (((r1 > r2) || (A.rowval[r1] != j)) ? zero(Tv) : A.nzval[r1], j+1) +function diag(A::SparseMatrixCSC{Tv,Ti}, d::Int=0) where {Tv,Ti} + m, n = size(A) + if !(-m <= d <= n) + throw(ArgumentError("requested diagonal, $d, out of bounds in matrix of size ($m, $n)")) + end + l = d < 0 ? min(m+d,n) : min(n-d,m) + ind = Vector{Ti}(); sizehint!(ind, min(l, nnz(A))) + val = Vector{Tv}(); sizehint!(val, min(l, nnz(A))) + for (i, v) in enumerate(SpDiagIterator(A, d)) + if !iszero(v) + push!(ind, i) + push!(val, v) + end + end + return SparseVector{Tv,Ti}(l, ind, val) end function trace(A::SparseMatrixCSC{Tv}) where Tv @@ -3406,7 +3421,7 @@ function trace(A::SparseMatrixCSC{Tv}) where Tv throw(DimensionMismatch("expected square matrix")) end s = zero(Tv) - for d in SpDiagIterator(A) + for d in SpDiagIterator(A,0) s += d end s diff --git a/test/sparse/sparse.jl b/test/sparse/sparse.jl index 52eb705dfc2a9..f47427e0e0bb6 100644 --- a/test/sparse/sparse.jl +++ b/test/sparse/sparse.jl @@ -1325,11 +1325,18 @@ end @testset "diag" begin for T in (Float64, Complex128) - S = sprand(T, 5, 5, 0.5) - A = Matrix(S) - @test diag(S)::SparseVector{T,Int} == diag(A) - @test diag(S, 1)::SparseVector{T,Int} == diag(A, 1) - @test diag(S, -1)::SparseVector{T,Int} == diag(A, -1) + S1 = sprand(T, 5, 5, 0.5) + S2 = sprand(T, 10, 5, 0.5) + S3 = sprand(T, 5, 10, 0.5) + for S in (S1, S2, S3) + A = Matrix(S) + @test diag(S)::SparseVector{T,Int} == diag(A) + for k in -5:5 + @test diag(S, k)::SparseVector{T,Int} == diag(A, k) + end + @test_throws ArgumentError diag(S, -size(S,1)-1) + @test_throws ArgumentError diag(S, size(S,2)+1) + end end end From 00905db51cd679a724f18b0fb894c5e5d7a5a4da Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Thu, 17 Aug 2017 02:23:06 +0200 Subject: [PATCH 3/6] remove SpDiagIterator --- NEWS.md | 2 ++ base/sparse/sparsematrix.jl | 55 +++++++++++++++---------------------- test/sparse/sparse.jl | 2 +- 3 files changed, 25 insertions(+), 34 deletions(-) diff --git a/NEWS.md b/NEWS.md index 4f001ae743fec..0be83e8075b2d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -312,6 +312,8 @@ Deprecated or removed * `Base.cpad` has been removed; use an appropriate combination of `rpad` and `lpad` instead ([#23187]). + * `Base.SparseArrays.SpDiagIterator` has been removed ([#23261]). + Command-line option changes --------------------------- diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index be4a55f5df93b..43c711669e216 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -3380,51 +3380,40 @@ function expandptr(V::Vector{<:Integer}) res end -## diag and related using an iterator - -struct SpDiagIterator{Tv,Ti} - A::SparseMatrixCSC{Tv,Ti} - d::Int # diagonal to iterate over -end - -start(it::SpDiagIterator) = it.d < 0 ? (1-it.d, 1) : (1, it.d+1) -done(it::SpDiagIterator, rc) = rc[1] > size(it.A,1) || rc[2] > size(it.A,2) -function next(it::SpDiagIterator{Tv}, rc) where Tv - r, c = rc - A = it.A - r1 = Int(A.colptr[c]) - r2 = Int(A.colptr[c+1]-1) - (r1 > r2) && (return (zero(Tv), (r+1, c+1))) - r1 = searchsortedfirst(A.rowval, r, r1, r2, Forward) - (((r1 > r2) || (A.rowval[r1] != r)) ? zero(Tv) : A.nzval[r1], (r+1, c+1)) -end function diag(A::SparseMatrixCSC{Tv,Ti}, d::Int=0) where {Tv,Ti} m, n = size(A) if !(-m <= d <= n) throw(ArgumentError("requested diagonal, $d, out of bounds in matrix of size ($m, $n)")) end - l = d < 0 ? min(m+d,n) : min(n-d,m) - ind = Vector{Ti}(); sizehint!(ind, min(l, nnz(A))) - val = Vector{Tv}(); sizehint!(val, min(l, nnz(A))) - for (i, v) in enumerate(SpDiagIterator(A, d)) - if !iszero(v) - push!(ind, i) - push!(val, v) - end + if d <= 0 + rrange = (1-d):min(m, min(m,n)-d) + crange = 1:min(n, m+d) + else # d > 0 + rrange = 1:min(m, n-d) + crange = (1+d):min(n, min(m,n)+d) + end + ind = Vector{Ti}() + val = Vector{Tv}() + for (i, (r, c)) in enumerate(zip(rrange, crange)) + r1 = Int(A.colptr[c]) + r2 = Int(A.colptr[c+1]-1) + r1 > r2 && continue + r1 = searchsortedfirst(A.rowval, r, r1, r2, Forward) + ((r1 > r2) || (A.rowval[r1] != r)) && continue + push!(ind, i) + push!(val, A.nzval[r1]) end - return SparseVector{Tv,Ti}(l, ind, val) + return SparseVector{Tv,Ti}(length(rrange), ind, val) end function trace(A::SparseMatrixCSC{Tv}) where Tv - if size(A,1) != size(A,2) - throw(DimensionMismatch("expected square matrix")) - end + n = checksquare(A) s = zero(Tv) - for d in SpDiagIterator(A,0) - s += d + for i in 1:n + s += A[i,i] end - s + return s end function diagm(v::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} diff --git a/test/sparse/sparse.jl b/test/sparse/sparse.jl index f47427e0e0bb6..79fd719a194d4 100644 --- a/test/sparse/sparse.jl +++ b/test/sparse/sparse.jl @@ -1331,7 +1331,7 @@ end for S in (S1, S2, S3) A = Matrix(S) @test diag(S)::SparseVector{T,Int} == diag(A) - for k in -5:5 + for k in -size(S,1):size(S,2) @test diag(S, k)::SparseVector{T,Int} == diag(A, k) end @test_throws ArgumentError diag(S, -size(S,1)-1) From a77a2066e19a535f7983a97087db401023926bba Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Thu, 17 Aug 2017 10:11:23 +0200 Subject: [PATCH 4/6] Int -> Integer, simplifications --- base/sparse/sparsematrix.jl | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index 43c711669e216..c2a8b3a0c165b 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -3381,21 +3381,17 @@ function expandptr(V::Vector{<:Integer}) end -function diag(A::SparseMatrixCSC{Tv,Ti}, d::Int=0) where {Tv,Ti} +function diag(A::SparseMatrixCSC{Tv,Ti}, d::Integer=0) where {Tv,Ti} m, n = size(A) if !(-m <= d <= n) throw(ArgumentError("requested diagonal, $d, out of bounds in matrix of size ($m, $n)")) end - if d <= 0 - rrange = (1-d):min(m, min(m,n)-d) - crange = 1:min(n, m+d) - else # d > 0 - rrange = 1:min(m, n-d) - crange = (1+d):min(n, min(m,n)+d) - end + l = d < 0 ? min(m+d,n) : min(n-d,m) + r, c = d <= 0 ? (-d, 0) : (0, d) # start row/col -1 ind = Vector{Ti}() val = Vector{Tv}() - for (i, (r, c)) in enumerate(zip(rrange, crange)) + for i in 1:l + r += 1; c += 1 r1 = Int(A.colptr[c]) r2 = Int(A.colptr[c+1]-1) r1 > r2 && continue @@ -3404,7 +3400,7 @@ function diag(A::SparseMatrixCSC{Tv,Ti}, d::Int=0) where {Tv,Ti} push!(ind, i) push!(val, A.nzval[r1]) end - return SparseVector{Tv,Ti}(length(rrange), ind, val) + return SparseVector{Tv,Ti}(l, ind, val) end function trace(A::SparseMatrixCSC{Tv}) where Tv From 17eef91775db066f46a9329d6462f3f21203cfd1 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Thu, 17 Aug 2017 11:28:50 +0200 Subject: [PATCH 5/6] test that stored zeros are still stored zeros in the diag --- test/sparse/sparse.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/sparse/sparse.jl b/test/sparse/sparse.jl index 79fd719a194d4..3753f49b8a9d0 100644 --- a/test/sparse/sparse.jl +++ b/test/sparse/sparse.jl @@ -1338,6 +1338,10 @@ end @test_throws ArgumentError diag(S, size(S,2)+1) end 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 V.nzind == [1,3] + @test V.nzval == [0.0,0.0] end @testset "expandptr" begin From 5530e4d03b9edb46bd38bf233a40137dde754b2c Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Thu, 17 Aug 2017 22:14:08 +0200 Subject: [PATCH 6/6] remove type instability for non Int d --- base/sparse/sparsematrix.jl | 9 +++++---- test/sparse/sparse.jl | 4 ++-- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index c2a8b3a0c165b..0c501c8c72bfa 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -3383,11 +3383,12 @@ end function diag(A::SparseMatrixCSC{Tv,Ti}, d::Integer=0) where {Tv,Ti} m, n = size(A) - if !(-m <= d <= n) - throw(ArgumentError("requested diagonal, $d, out of bounds in matrix of size ($m, $n)")) + k = Int(d) + if !(-m <= k <= n) + throw(ArgumentError("requested diagonal, $k, out of bounds in matrix of size ($m, $n)")) end - l = d < 0 ? min(m+d,n) : min(n-d,m) - r, c = d <= 0 ? (-d, 0) : (0, d) # start row/col -1 + l = k < 0 ? min(m+k,n) : min(n-k,m) + r, c = k <= 0 ? (-k, 0) : (0, k) # start row/col -1 ind = Vector{Ti}() val = Vector{Tv}() for i in 1:l diff --git a/test/sparse/sparse.jl b/test/sparse/sparse.jl index 3753f49b8a9d0..38e25d0a3a5a8 100644 --- a/test/sparse/sparse.jl +++ b/test/sparse/sparse.jl @@ -1330,9 +1330,9 @@ end S3 = sprand(T, 5, 10, 0.5) for S in (S1, S2, S3) A = Matrix(S) - @test diag(S)::SparseVector{T,Int} == diag(A) + @test diag(S)::SparseVector{T,Int} == diag(A) for k in -size(S,1):size(S,2) - @test diag(S, k)::SparseVector{T,Int} == diag(A, k) + @test diag(S, k)::SparseVector{T,Int} == diag(A, k) end @test_throws ArgumentError diag(S, -size(S,1)-1) @test_throws ArgumentError diag(S, size(S,2)+1)