Skip to content

Commit

Permalink
multiplication of sparse triangular matrices #35609 #35610 #35642 (#3…
Browse files Browse the repository at this point in the history 
…5659)

* multiplication of sparse triangular matrices
* removed disambiguity and improved test cases
* test cases for `nnz, nzrange, rowvals, nonzeros`
  • Loading branch information
@KlausC
KlausC committed Jun 1, 2020
1 parent 932a1ec commit cfb9b55
Show file tree
Hide file tree
Showing 5 changed files with 223 additions and 62 deletions.
137 changes: 85 additions & 52 deletions stdlib/SparseArrays/src/linalg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,9 +44,9 @@ function mul!(C::StridedVecOrMat, A::AbstractSparseMatrixCSC, B::Union{StridedVe
end
C
end
*(A::AbstractSparseMatrixCSC{TA}, x::StridedVector{Tx}) where {TA,Tx} =
*(A::SparseMatrixCSCUnion{TA}, x::StridedVector{Tx}) where {TA,Tx} =
(T = promote_op(matprod, TA, Tx); mul!(similar(x, T, size(A, 1)), A, x, true, false))
*(A::AbstractSparseMatrixCSC{TA}, B::AdjOrTransStridedOrTriangularMatrix{Tx}) where {TA,Tx} =
*(A::SparseMatrixCSCUnion{TA}, B::AdjOrTransStridedOrTriangularMatrix{Tx}) where {TA,Tx} =
(T = promote_op(matprod, TA, Tx); mul!(similar(B, T, (size(A, 1), size(B, 2))), A, B, true, false))

function mul!(C::StridedVecOrMat, adjA::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, B::Union{StridedVector,AdjOrTransStridedOrTriangularMatrix}, α::Number, β::Number)
Expand Down Expand Up @@ -125,7 +125,7 @@ function mul!(C::StridedVecOrMat, X::AdjOrTransStridedOrTriangularMatrix, A::Abs
end
C
end
*(X::AdjOrTransStridedOrTriangularMatrix, A::AbstractSparseMatrixCSC{TvA}) where {TvA} =
*(X::AdjOrTransStridedOrTriangularMatrix, A::SparseMatrixCSCUnion{TvA}) where {TvA} =
(T = promote_op(matprod, eltype(X), TvA); mul!(similar(X, T, (size(X, 1), size(A, 2))), X, A, true, false))

function mul!(C::StridedVecOrMat, X::AdjOrTransStridedOrTriangularMatrix, adjA::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, α::Number, β::Number)
Expand Down Expand Up @@ -182,11 +182,20 @@ end
# Sparse matrix multiplication as described in [Gustavson, 1978]:
# http://dl.acm.org/citation.cfm?id=355796

*(A::AbstractSparseMatrixCSC, B::AbstractSparseMatrixCSC) = spmatmul(A,B)
*(A::AbstractSparseMatrixCSC, B::Adjoint{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(A, copy(B))
*(A::AbstractSparseMatrixCSC, B::Transpose{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(A, copy(B))
*(A::Transpose{<:Any,<:AbstractSparseMatrixCSC}, B::AbstractSparseMatrixCSC) = spmatmul(copy(A), B)
*(A::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, B::AbstractSparseMatrixCSC) = spmatmul(copy(A), B)
const SparseTriangular{Tv,Ti} = Union{UpperTriangular{Tv,<:SparseMatrixCSCUnion{Tv,Ti}},LowerTriangular{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}}
const SparseOrTri{Tv,Ti} = Union{SparseMatrixCSCUnion{Tv,Ti},SparseTriangular{Tv,Ti}}

*(A::SparseOrTri, B::AbstractSparseVector) = spmatmulv(A, B)
*(A::SparseOrTri, B::SparseColumnView) = spmatmulv(A, B)
*(A::SparseOrTri, B::SparseVectorView) = spmatmulv(A, B)
*(A::SparseMatrixCSCUnion, B::SparseMatrixCSCUnion) = spmatmul(A,B)
*(A::SparseTriangular, B::SparseMatrixCSCUnion) = spmatmul(A,B)
*(A::SparseMatrixCSCUnion, B::SparseTriangular) = spmatmul(A,B)
*(A::SparseTriangular, B::SparseTriangular) = spmatmul1(A,B)
*(A::SparseOrTri, B::Adjoint{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(A, copy(B))
*(A::SparseOrTri, B::Transpose{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(A, copy(B))
*(A::Transpose{<:Any,<:AbstractSparseMatrixCSC}, B::SparseOrTri) = spmatmul(copy(A), B)
*(A::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, B::SparseOrTri) = spmatmul(copy(A), B)
*(A::Adjoint{<:Any,<:AbstractSparseMatrixCSC}, B::Adjoint{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(copy(A), copy(B))
*(A::Transpose{<:Any,<:AbstractSparseMatrixCSC}, B::Transpose{<:Any,<:AbstractSparseMatrixCSC}) = spmatmul(copy(A), copy(B))

Expand All @@ -198,16 +207,15 @@ end
# done by a quicksort of the row indices or by a full scan of the dense result vector.
# The last is faster, if more than ≈ 1/32 of the result column is nonzero.
# TODO: extend to SparseMatrixCSCUnion to allow for SubArrays (view(X, :, r)).
function spmatmul(A::AbstractSparseMatrixCSC{TvA,TiA}, B::AbstractSparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB}
function spmatmul(A::SparseOrTri{TvA,TiA}, B::Union{SparseOrTri{TvB,TiB},SparseVectorUnion{TvB,TiB},SubArray{TvB,<:Any,<:AbstractSparseArray{TvB,TiB}}}) where {TvA,TiA,TvB,TiB}

Tv = promote_op(matprod, TvA, TvB)
Ti = promote_type(TiA, TiB)
mA, nA = size(A)
nB = size(B, 2)
nA == size(B, 1) || throw(DimensionMismatch())

rowvalA = rowvals(A); nzvalA = nonzeros(A)
rowvalB = rowvals(B); nzvalB = nonzeros(B)
nnzC = max(estimate_mulsize(mA, nnz(A), nA, nnz(B), nB) * 11 ÷ 10, mA)
nnzC = min(estimate_mulsize(mA, nnz(A), nA, nnz(B), nB) * 11 ÷ 10 + mA, mA*nB)
colptrC = Vector{Ti}(undef, nB+1)
rowvalC = Vector{Ti}(undef, nnzC)
nzvalC = Vector{Tv}(undef, nnzC)
Expand All @@ -221,45 +229,8 @@ function spmatmul(A::AbstractSparseMatrixCSC{TvA,TiA}, B::AbstractSparseMatrixCS
resize!(rowvalC, nnzC)
resize!(nzvalC, nnzC)
end
colptrC[i] = ip0 = ip
k0 = ip - 1
for jp in nzrange(B, i)
nzB = nzvalB[jp]
j = rowvalB[jp]
for kp in nzrange(A, j)
nzC = nzvalA[kp] * nzB
k = rowvalA[kp]
if xb[k]
nzvalC[k+k0] += nzC
else
nzvalC[k+k0] = nzC
xb[k] = true
rowvalC[ip] = k
ip += 1
end
end
end
if ip > ip0
if prefer_sort(ip-k0, mA)
# in-place sort of indices. Effort: O(nnz*ln(nnz)).
sort!(rowvalC, ip0, ip-1, QuickSort, Base.Order.Forward)
for vp = ip0:ip-1
k = rowvalC[vp]
xb[k] = false
nzvalC[vp] = nzvalC[k+k0]
end
else
# scan result vector (effort O(mA))
for k = 1:mA
if xb[k]
xb[k] = false
rowvalC[ip0] = k
nzvalC[ip0] = nzvalC[k+k0]
ip0 += 1
end
end
end
end
colptrC[i] = ip
ip = spcolmul!(rowvalC, nzvalC, xb, i, ip, A, B)
end
colptrC[nB+1] = ip
end
Expand All @@ -272,6 +243,68 @@ function spmatmul(A::AbstractSparseMatrixCSC{TvA,TiA}, B::AbstractSparseMatrixCS
return C
end

# process single rhs column
function spcolmul!(rowvalC, nzvalC, xb, i, ip, A, B)
rowvalA = rowvals(A); nzvalA = nonzeros(A)
rowvalB = rowvals(B); nzvalB = nonzeros(B)
mA = size(A, 1)
ip0 = ip
k0 = ip - 1
@inbounds begin
for jp in nzrange(B, i)
nzB = nzvalB[jp]
j = rowvalB[jp]
for kp in nzrange(A, j)
nzC = nzvalA[kp] * nzB
k = rowvalA[kp]
if xb[k]
nzvalC[k+k0] += nzC
else
nzvalC[k+k0] = nzC
xb[k] = true
rowvalC[ip] = k
ip += 1
end
end
end
if ip > ip0
if prefer_sort(ip-k0, mA)
# in-place sort of indices. Effort: O(nnz*ln(nnz)).
sort!(rowvalC, ip0, ip-1, QuickSort, Base.Order.Forward)
for vp = ip0:ip-1
k = rowvalC[vp]
xb[k] = false
nzvalC[vp] = nzvalC[k+k0]
end
else
# scan result vector (effort O(mA))
for k = 1:mA
if xb[k]
xb[k] = false
rowvalC[ip0] = k
nzvalC[ip0] = nzvalC[k+k0]
ip0 += 1
end
end
end
end
end
return ip
end

# special cases of same twin Upper/LowerTriangular
spmatmul1(A, B) = spmatmul(A, B)
function spmatmul1(A::UpperTriangular, B::UpperTriangular)
UpperTriangular(spmatmul(A, B))
end
function spmatmul1(A::LowerTriangular, B::LowerTriangular)
LowerTriangular(spmatmul(A, B))
end
# exploit spmatmul for sparse vectors and column views
function spmatmulv(A, B)
spmatmul(A, B)[:,1]
end

# estimated number of non-zeros in matrix product
# it is assumed, that the non-zero indices are distributed independently and uniformly
# in both matrices. Over-estimation is possible if that is not the case.
Expand Down Expand Up @@ -765,7 +798,7 @@ function _ldiv!(U::UpperTriangularWrapped, B::StridedVecOrMat)
end

(\)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = ldiv!(L, Array(B))
(*)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B))
#(*)(L::TriangularSparse, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B))

## end of triangular

Expand Down
17 changes: 15 additions & 2 deletions stdlib/SparseArrays/src/sparsematrix.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,8 +109,15 @@ julia> nnz(A)
```
"""
nnz(S::AbstractSparseMatrixCSC) = Int(getcolptr(S)[size(S, 2) + 1] - 1)
nnz(S::ReshapedArray{T,1,<:AbstractSparseMatrixCSC}) where T = nnz(parent(S))
count(pred, S::AbstractSparseMatrixCSC) = count(pred, nzvalview(S)) + pred(zero(eltype(S)))*(prod(size(S)) - nnz(S))
nnz(S::ReshapedArray{<:Any,1,<:AbstractSparseMatrixCSC}) = nnz(parent(S))
nnz(S::UpperTriangular{<:Any,<:AbstractSparseMatrixCSC}) = nnz1(S)
nnz(S::LowerTriangular{<:Any,<:AbstractSparseMatrixCSC}) = nnz1(S)
nnz(S::SparseMatrixCSCView) = nnz1(S)
nnz1(S) = sum(length.(nzrange.(Ref(S), axes(S, 2))))

function count(pred, S::AbstractSparseMatrixCSC)
count(pred, nzvalview(S)) + pred(zero(eltype(S)))*(prod(size(S)) - nnz(S))
end

"""
nonzeros(A)
Expand Down Expand Up @@ -138,6 +145,8 @@ julia> nonzeros(A)
"""
nonzeros(S::SparseMatrixCSC) = getfield(S, :nzval)
nonzeros(S::SparseMatrixCSCView) = nonzeros(S.parent)
nonzeros(S::UpperTriangular{<:Any,<:SparseMatrixCSCUnion}) = nonzeros(S.data)
nonzeros(S::LowerTriangular{<:Any,<:SparseMatrixCSCUnion}) = nonzeros(S.data)

"""
rowvals(A::AbstractSparseMatrixCSC)
Expand All @@ -164,6 +173,8 @@ julia> rowvals(A)
"""
rowvals(S::SparseMatrixCSC) = getfield(S, :rowval)
rowvals(S::SparseMatrixCSCView) = rowvals(S.parent)
rowvals(S::UpperTriangular{<:Any,<:SparseMatrixCSCUnion}) = rowvals(S.data)
rowvals(S::LowerTriangular{<:Any,<:SparseMatrixCSCUnion}) = rowvals(S.data)

"""
nzrange(A::AbstractSparseMatrixCSC, col::Integer)
Expand All @@ -186,6 +197,8 @@ column. In conjunction with [`nonzeros`](@ref) and
"""
nzrange(S::AbstractSparseMatrixCSC, col::Integer) = getcolptr(S)[col]:(getcolptr(S)[col+1]-1)
nzrange(S::SparseMatrixCSCView, col::Integer) = nzrange(S.parent, S.indices[2][col])
nzrange(S::UpperTriangular{<:Any,<:SparseMatrixCSCUnion}, i::Integer) = nzrangeup(S.data, i)
nzrange(S::LowerTriangular{<:Any,<:SparseMatrixCSCUnion}, i::Integer) = nzrangelo(S.data, i)

function Base.isstored(A::AbstractSparseMatrixCSC, i::Integer, j::Integer)
@boundscheck checkbounds(A, i, j)
Expand Down
36 changes: 28 additions & 8 deletions stdlib/SparseArrays/src/sparsevector.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,23 +32,43 @@ SparseVector(n::Integer, nzind::Vector{Ti}, nzval::Vector{Tv}) where {Tv,Ti} =

# Define an alias for a view of a whole column of a SparseMatrixCSC. Many methods can be written for the
# union of such a view and a SparseVector so we define an alias for such a union as well
const SparseColumnView{T} = SubArray{T,1,<:AbstractSparseMatrixCSC,Tuple{Base.Slice{Base.OneTo{Int}},Int},false}
const SparseVectorUnion{T} = Union{SparseVector{T}, SparseColumnView{T}}
const AdjOrTransSparseVectorUnion{T} = LinearAlgebra.AdjOrTrans{T, <:SparseVectorUnion{T}}
const SparseColumnView{Tv,Ti} = SubArray{Tv,1,<:AbstractSparseMatrixCSC{Tv,Ti},Tuple{Base.Slice{Base.OneTo{Int}},Int},false}
const SparseVectorView{Tv,Ti} = SubArray{Tv,1,<:AbstractSparseVector{Tv,Ti},Tuple{Base.Slice{Base.OneTo{Int}}},false}
const SparseVectorUnion{Tv,Ti} = Union{SparseVector{Tv,Ti}, SparseColumnView{Tv,Ti}, SparseVectorView{Tv,Ti}}
const AdjOrTransSparseVectorUnion{Tv,Ti} = LinearAlgebra.AdjOrTrans{Tv, <:SparseVectorUnion{Tv,Ti}}

### Basic properties

size(x::SparseVector) = (getfield(x, :n),)
nnz(x::SparseVector) = length(nonzeros(x))
count(f, x::SparseVector) = count(f, nonzeros(x)) + f(zero(eltype(x)))*(length(x) - nnz(x))

# implement the nnz - nzrange - nonzeros - rowvals interface for sparse vectors

nnz(x::SparseVector) = length(nonzeros(x))
function nnz(x::SparseColumnView)
rowidx, colidx = parentindices(x)
return length(nzrange(parent(x), colidx))
end
nnz(x::SparseVectorView) = nnz(x.parent)

"""
nzrange(x::SparseVectorUnion, col)
Give the range of indices to the structural nonzero values of a sparse vector.
The column index `col` is ignored (assumed to be `1`).
"""
function nzrange(x::SparseVectorUnion, j::Integer)
j == 1 ? (1:nnz(x)) : throw(BoundsError(x, (":", j)))
end

nonzeros(x::SparseVector) = getfield(x, :nzval)
function nonzeros(x::SparseColumnView)
rowidx, colidx = parentindices(x)
A = parent(x)
@inbounds y = view(nonzeros(A), nzrange(A, colidx))
return y
end
nonzeros(x::SparseVectorView) = nonzeros(parent(x))

nonzeroinds(x::SparseVector) = getfield(x, :nzind)
function nonzeroinds(x::SparseColumnView)
Expand All @@ -57,12 +77,12 @@ function nonzeroinds(x::SparseColumnView)
@inbounds y = view(rowvals(A), nzrange(A, colidx))
return y
end
nonzeroinds(x::SparseVectorView) = nonzeroinds(parent(x))

rowvals(x::SparseVectorUnion) = nonzeroinds(x)

indtype(x::SparseColumnView) = indtype(parent(x))
function nnz(x::SparseColumnView)
rowidx, colidx = parentindices(x)
return length(nzrange(parent(x), colidx))
end
indtype(x::SparseVectorView) = indtype(parent(x))

## similar
#
Expand Down
80 changes: 80 additions & 0 deletions stdlib/SparseArrays/test/sparse.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2996,4 +2996,84 @@ end
end
end

@testset "Multiplying with triangular sparse matrices #35609 #35610" begin
n = 10
A = sprand(n, n, 5/n)
U = UpperTriangular(A)
L = LowerTriangular(A)
AM = Matrix(A)
UM = Matrix(U)
LM = Matrix(L)
Y = A * U
@test Y ≈ AM * UM
@test typeof(Y) == typeof(A)
Y = A * L
@test Y ≈ AM * LM
@test typeof(Y) == typeof(A)
Y = U * A
@test Y ≈ UM * AM
@test typeof(Y) == typeof(A)
Y = L * A
@test Y ≈ LM * AM
@test typeof(Y) == typeof(A)
Y = U * U
@test Y ≈ UM * UM
@test typeof(Y) == typeof(U)
Y = L * L
@test Y ≈ LM * LM
@test typeof(Y) == typeof(L)
Y = L * U
@test Y ≈ LM * UM
@test typeof(Y) == typeof(A)
Y = U * L
@test Y ≈ UM * LM
@test typeof(Y) == typeof(A)
end

#testing the sparse matrix/vector access functions nnz, nzrange, rowvals, nonzeros
@testset "generic sparse matrix access functions" begin
I = [1,3,4,5, 1,3,4,5, 1,3,4,5];
J = [4,4,4,4, 5,5,5,5, 6,6,6,6];
V = [14,34,44,54, 15,35,45,55, 16,36,46,56];
A = sparse(I, J, V, 9, 9);
AU = UpperTriangular(A)
AL = LowerTriangular(A)
b = SparseVector(9, I[1:4], V[1:4])
c = view(A, :, 5)
d = view(b, :)

@testset "nnz $n" for (n, M, nz) in (("A", A, 12), ("AU", AU, 11), ("AL", AL, 3),
("b", b, 4), ("c", c, 4), ("d", d, 4))
@test nnz(M) == nz
@test_throws BoundsError nzrange(M, 0)
@test_throws BoundsError nzrange(M, size(M, 2) + 1)
end
@testset "nzrange(A, $i)" for (i, nzr) in ((1,1:0),(4,1:4),(5,5:8),(6,9:12),(9,13:12))
@test nzrange(A, i) == nzr
end
@testset "nzrange(AU, $i)" for (i, nzr) in ((2,1:0),(4,1:3),(5,5:8),(6,9:12),(8,13:12))
@test nzrange(AU, i) == nzr
end
@testset "nzrange(AL, $i)" for (i, nzr) in ((3,1:0),(4,3:4),(5,8:8),(6,13:12),(7,13:12))
@test nzrange(AL, i) == nzr
end
@test nzrange(b, 1) == 1:4
@test nzrange(c, 1) == 1:4
@test nzrange(d, 1) == 1:4

@test rowvals(A) == I
@test rowvals(AL) == I
@test rowvals(AL) == I
@test rowvals(b) == I[1:4]
@test rowvals(c) == I[5:8]
@test rowvals(d) == I[1:4]

@test nonzeros(A) == V
@test nonzeros(AU) == V
@test nonzeros(AL) == V
@test nonzeros(b) == V[1:4]
@test nonzeros(c) == V[5:8]
@test nonzeros(d) == V[1:4]
end

end # module
Loading

2 comments on commit cfb9b55

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Executing the daily package evaluation, I will reply here when finished:

@nanosoldier runtests(ALL, isdaily = true)

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Your package evaluation job has completed - possible new issues were detected. A full report can be found here. cc @maleadt

Please sign in to comment.