Skip to content

Commit ee9ae74

Browse files
Sacha0Sacha0
committed
Replace the LGPL-licensed sparse() parent method with an MIT-licensed version. See #13001 and #14631.
1 parent 2f5975c commit ee9ae74

File tree

5 files changed

+240
-130
lines changed

5 files changed

+240
-130
lines changed

‎NEWS.md

+2
Original file line numberDiff line numberDiff line change
@@ -91,6 +91,8 @@ Library improvements
9191
* Rank one update and downdate functions, `lowrankupdate`, `lowrankupdate!`, `lowrankdowndate`,
9292
and `lowrankdowndate!`, for dense Cholesky factorizations ([#14243],[#14424])
9393

94+
* All `sparse` methods now retain provided numerical zeros as structural nonzeros.
95+
9496
* New `foreach` function for calling a function on every element of a collection when
9597
the results are not needed.
9698

‎base/sparse/csparse.jl

-125
Original file line numberDiff line numberDiff line change
@@ -13,131 +13,6 @@
1313
# Section 2.4: Triplet form
1414
# http://www.cise.ufl.edu/research/sparse/CSparse/
1515

16-
"""
17-
sparse(I,J,V,[m,n,combine])
18-
19-
Create a sparse matrix `S` of dimensions `m x n` such that `S[I[k], J[k]] = V[k]`.
20-
The `combine` function is used to combine duplicates. If `m` and `n` are not
21-
specified, they are set to `maximum(I)` and `maximum(J)` respectively. If the
22-
`combine` function is not supplied, `combine` defaults to `+` unless the elements
23-
of `V` are Booleans in which case `combine` defaults to `|`. All elements of `I` must
24-
satisfy `1 <= I[k] <= m`, and all elements of `J` must satisfy `1 <= J[k] <= n`.
25-
"""
26-
function sparse{Tv,Ti<:Integer}(I::AbstractVector{Ti},
27-
J::AbstractVector{Ti},
28-
V::AbstractVector{Tv},
29-
nrow::Integer, ncol::Integer,
30-
combine::Union{Function,Base.Func})
31-
32-
N = length(I)
33-
if N != length(J) || N != length(V)
34-
throw(ArgumentError("triplet I,J,V vectors must be the same length"))
35-
end
36-
if N == 0
37-
return spzeros(eltype(V), Ti, nrow, ncol)
38-
end
39-
40-
# Work array
41-
Wj = Array(Ti, max(nrow,ncol)+1)
42-
43-
# Allocate sparse matrix data structure
44-
# Count entries in each row
45-
Rnz = zeros(Ti, nrow+1)
46-
Rnz[1] = 1
47-
nz = 0
48-
for k=1:N
49-
iind = I[k]
50-
iind > 0 || throw(ArgumentError("all I index values must be > 0"))
51-
iind <= nrow || throw(ArgumentError("all I index values must be ≤ the number of rows"))
52-
if V[k] != 0
53-
Rnz[iind+1] += 1
54-
nz += 1
55-
end
56-
end
57-
Rp = cumsum(Rnz)
58-
Ri = Array(Ti, nz)
59-
Rx = Array(Tv, nz)
60-
61-
# Construct row form
62-
# place triplet (i,j,x) in column i of R
63-
# Use work array for temporary row pointers
64-
@simd for i=1:nrow; @inbounds Wj[i] = Rp[i]; end
65-
66-
@inbounds for k=1:N
67-
iind = I[k]
68-
jind = J[k]
69-
jind > 0 || throw(ArgumentError("all J index values must be > 0"))
70-
jind <= ncol || throw(ArgumentError("all J index values must be ≤ the number of columns"))
71-
p = Wj[iind]
72-
Vk = V[k]
73-
if Vk != 0
74-
Wj[iind] += 1
75-
Rx[p] = Vk
76-
Ri[p] = jind
77-
end
78-
end
79-
80-
# Reset work array for use in counting duplicates
81-
@simd for j=1:ncol; @inbounds Wj[j] = 0; end
82-
83-
# Sum up duplicates and squeeze
84-
anz = 0
85-
@inbounds for i=1:nrow
86-
p1 = Rp[i]
87-
p2 = Rp[i+1] - 1
88-
pdest = p1
89-
90-
for p = p1:p2
91-
j = Ri[p]
92-
pj = Wj[j]
93-
if pj >= p1
94-
Rx[pj] = combine(Rx[pj], Rx[p])
95-
else
96-
Wj[j] = pdest
97-
if pdest != p
98-
Ri[pdest] = j
99-
Rx[pdest] = Rx[p]
100-
end
101-
pdest += one(Ti)
102-
end
103-
end
104-
105-
Rnz[i] = pdest - p1
106-
anz += (pdest - p1)
107-
end
108-
109-
# Transpose from row format to get the CSC format
110-
RiT = Array(Ti, anz)
111-
RxT = Array(Tv, anz)
112-
113-
# Reset work array to build the final colptr
114-
Wj[1] = 1
115-
@simd for i=2:(ncol+1); @inbounds Wj[i] = 0; end
116-
@inbounds for j = 1:nrow
117-
p1 = Rp[j]
118-
p2 = p1 + Rnz[j] - 1
119-
for p = p1:p2
120-
Wj[Ri[p]+1] += 1
121-
end
122-
end
123-
RpT = cumsum(Wj[1:(ncol+1)])
124-
125-
# Transpose
126-
@simd for i=1:length(RpT); @inbounds Wj[i] = RpT[i]; end
127-
@inbounds for j = 1:nrow
128-
p1 = Rp[j]
129-
p2 = p1 + Rnz[j] - 1
130-
for p = p1:p2
131-
ind = Ri[p]
132-
q = Wj[ind]
133-
Wj[ind] += 1
134-
RiT[q] = j
135-
RxT[q] = Rx[p]
136-
end
137-
end
138-
139-
return SparseMatrixCSC(nrow, ncol, RpT, RiT, RxT)
140-
end
14116

14217
# Compute the elimination tree of A using triu(A) returning the parent vector.
14318
# A root node is indicated by 0. This tree may actually be a forest in that

‎base/sparse/sparsematrix.jl

+219-1
Original file line numberDiff line numberDiff line change
@@ -335,7 +335,225 @@ function sparse_IJ_sorted!{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector
335335
return SparseMatrixCSC(m, n, colptr, I, V)
336336
end
337337

338-
## sparse() can take its inputs in unsorted order (the parent method is now in csparse.jl)
338+
"""
339+
sparse(I, J, V,[ m, n, combine])
340+
341+
Create a sparse matrix `S` of dimensions `m x n` such that `S[I[k], J[k]] = V[k]`. The
342+
`combine` function is used to combine duplicates. If `m` and `n` are not specified, they
343+
are set to `maximum(I)` and `maximum(J)` respectively. If the `combine` function is not
344+
supplied, `combine` defaults to `+` unless the elements of `V` are Booleans in which case
345+
`combine` defaults to `|`. All elements of `I` must satisfy `1 <= I[k] <= m`, and all
346+
elements of `J` must satisfy `1 <= J[k] <= n`. Numerical zeros in (`I`, `J`, `V`) are
347+
retained as structural nonzeros.
348+
349+
For additional documentation and an expert driver, see `Base.SparseArrays.sparse!`.
350+
"""
351+
function sparse{Tv,Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{Tv}, m::Integer, n::Integer, combine)
352+
coolen = length(I)
353+
if length(J) != coolen || length(V) != coolen
354+
throw(ArgumentError(string("the first three arguments' lengths must match, ",
355+
"length(I) (=$(length(I))) == length(J) (= $(length(J))) == length(V) (= ",
356+
"$(length(V)))")))
357+
end
358+
359+
if m == 0 || n == 0 || coolen == 0
360+
if coolen != 0
361+
if n == 0
362+
throw(ArgumentError("column indices J[k] must satisfy 1 <= J[k] <= n"))
363+
elseif m == 0
364+
throw(ArgumentError("row indices I[k] must satisfy 1 <= I[k] <= m"))
365+
end
366+
end
367+
SparseMatrixCSC(m, n, ones(Ti, n+1), Vector{Ti}(), Vector{Tv}())
368+
else
369+
# Allocate storage for CSR form
370+
csrrowptr = Vector{Ti}(m+1)
371+
csrcolval = Vector{Ti}(coolen)
372+
csrnzval = Vector{Tv}(coolen)
373+
374+
# Allocate storage for the CSC form's column pointers and a necessary workspace
375+
csccolptr = Vector{Ti}(n+1)
376+
klasttouch = Vector{Ti}(n)
377+
378+
# Allocate empty arrays for the CSC form's row and nonzero value arrays
379+
# The parent method called below automagically resizes these arrays
380+
cscrowval = Vector{Ti}()
381+
cscnzval = Vector{Tv}()
382+
383+
sparse!(I, J, V, m, n, combine, klasttouch,
384+
csrrowptr, csrcolval, csrnzval,
385+
csccolptr, cscrowval, cscnzval )
386+
end
387+
end
388+
389+
"""
390+
sparse!{Tv,Ti<:Integer}(
391+
I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{Tv},
392+
m::Integer, n::Integer, combine, klasttouch::Vector{Ti},
393+
csrrowptr::Vector{Ti}, csrcolval::Vector{Ti}, csrnzval::Vector{Tv},
394+
[csccolptr::Vector{Ti}], [cscrowval::Vector{Ti}, cscnzval::Vector{Tv}] )
395+
396+
Parent of and expert driver for `sparse`; see `sparse` for basic usage. This method
397+
allows the user to provide preallocated storage for `sparse`'s intermediate objects and
398+
result as described below. This capability enables more efficient successive construction
399+
of `SparseMatrixCSC`s from coordinate representations, and also enables extraction of an
400+
unsorted-column representation of the result's transpose at no additional cost.
401+
402+
This method consists of three major steps: (1) Counting-sort the provided coordinate
403+
representation into an unsorted-row CSR form including repeated entries. (2) Sweep through
404+
the CSR form, simultaneously calculating the desired CSC form's column-pointer array,
405+
detecting repeated entries, and repacking the CSR form with repeated entries combined;
406+
this stage yields an unsorted-row CSR form with no repeated entries. (3) Counting-sort the
407+
preceding CSR form into a fully-sorted CSC form with no repeated entries.
408+
409+
Input arrays `csrrowptr`, `csrcolval`, and `csrnzval` constitute storage for the
410+
intermediate CSR forms and require `length(csrrowptr) >= m + 1`,
411+
`length(csrcolval) >= length(I)`, and `length(csrnzval >= length(I))`. Input
412+
array `klasttouch`, workspace for the second stage, requires `length(klasttouch) >= n`.
413+
Optional input arrays `csccolptr`, `cscrowval`, and `cscnzval` constitute storage for the
414+
returned CSC form `S`. `csccolptr` requires `length(csccolptr) >= n + 1`. If necessary,
415+
`cscrowval` and `cscnzval` are automatically resized to satisfy
416+
`length(cscrowval) >= nnz(S)` and `length(cscnzval) >= nnz(S)`; hence, if `nnz(S)` is
417+
unknown at the outset, passing in empty vectors of the appropriate type (`Vector{Ti}()`
418+
and `Vector{Tv}()` respectively) suffices, or calling the `sparse!` method
419+
neglecting `cscrowval` and `cscnzval`.
420+
421+
On return, `csrrowptr`, `csrcolval`, and `csrnzval` contain an unsorted-column
422+
representation of the result's transpose.
423+
424+
You may reuse the input arrays' storage (`I`, `J`, `V`) for the output arrays
425+
(`csccolptr`, `cscrowval`, `cscnzval`). For example, you may call
426+
`sparse!(I, J, V, csrrowptr, csrcolval, csrnzval, I, J, V)`.
427+
428+
For the sake of efficiency, this method performs no argument checking beyond
429+
`1 <= I[k] <= m` and `1 <= J[k] <= n`. Use with care. Testing with `--check-bounds=yes`
430+
is wise.
431+
432+
This method runs in `O(m, n, length(I))` time. The HALFPERM algorithm described in
433+
F. Gustavson, "Two fast algorithms for sparse matrices: multiplication and permuted
434+
transposition," ACM TOMS 4(3), 250-269 (1978) inspired this method's use of a pair of
435+
counting sorts.
436+
437+
Performance note: As of January 2016, `combine` should be a functor for this method to
438+
perform well. This caveat may disappear when the work in `jb/functions` lands.
439+
"""
440+
function sparse!{Tv,Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti},
441+
V::AbstractVector{Tv}, m::Integer, n::Integer, combine, klasttouch::Vector{Ti},
442+
csrrowptr::Vector{Ti}, csrcolval::Vector{Ti}, csrnzval::Vector{Tv},
443+
csccolptr::Vector{Ti}, cscrowval::Vector{Ti}, cscnzval::Vector{Tv} )
444+
445+
# Compute the CSR form's row counts and store them shifted forward by one in csrrowptr
446+
fill!(csrrowptr, 0)
447+
coolen = length(I)
448+
@inbounds for k in 1:coolen
449+
Ik = I[k]
450+
if 1 > Ik || m < Ik
451+
throw(ArgumentError("row indices I[k] must satisfy 1 <= I[k] <= m"))
452+
end
453+
csrrowptr[Ik+1] += 1
454+
end
455+
456+
# Compute the CSR form's rowptrs and store them shifted forward by one in csrrowptr
457+
countsum = 1
458+
csrrowptr[1] = 1
459+
@inbounds for i in 2:(m+1)
460+
overwritten = csrrowptr[i]
461+
csrrowptr[i] = countsum
462+
countsum += overwritten
463+
end
464+
465+
# Counting-sort the column and nonzero values from J and V into csrcolval and csrnzval
466+
# Tracking write positions in csrrowptr corrects the row pointers
467+
@inbounds for k in 1:coolen
468+
Ik, Jk = I[k], J[k]
469+
if 1 > Jk || n < Jk
470+
throw(ArgumentError("column indices J[k] must satisfy 1 <= J[k] <= n"))
471+
end
472+
csrk = csrrowptr[Ik+1]
473+
csrrowptr[Ik+1] = csrk+1
474+
csrcolval[csrk] = Jk
475+
csrnzval[csrk] = V[k]
476+
end
477+
# This completes the unsorted-row, has-repeats CSR form's construction
478+
479+
# Sweep through the CSR form, simultaneously (1) caculating the CSC form's column
480+
# counts and storing them shifted forward by one in csccolptr; (2) detecting repeated
481+
# entries; and (3) repacking the CSR form with the repeated entries combined.
482+
#
483+
# Minimizing extraneous communication and nonlocality of reference, primarily by using
484+
# only a single auxiliary array in this step, is the key to this method's performance.
485+
fill!(csccolptr, 0)
486+
fill!(klasttouch, 0)
487+
writek = 1
488+
newcsrrowptri = 1
489+
origcsrrowptri = 1
490+
origcsrrowptrip1 = csrrowptr[2]
491+
@inbounds for i in 1:m
492+
for readk in origcsrrowptri:(origcsrrowptrip1-1)
493+
j = csrcolval[readk]
494+
if klasttouch[j] < newcsrrowptri
495+
klasttouch[j] = writek
496+
if writek != readk
497+
csrcolval[writek] = j
498+
csrnzval[writek] = csrnzval[readk]
499+
end
500+
writek += 1
501+
csccolptr[j+1] += 1
502+
else
503+
klt = klasttouch[j]
504+
csrnzval[klt] = combine(csrnzval[klt], csrnzval[readk])
505+
end
506+
end
507+
newcsrrowptri = writek
508+
origcsrrowptri = origcsrrowptrip1
509+
origcsrrowptrip1 != writek && (csrrowptr[i+1] = writek)
510+
i < m && (origcsrrowptrip1 = csrrowptr[i+2])
511+
end
512+
513+
# Compute the CSC form's colptrs and store them shifted forward by one in csccolptr
514+
countsum = 1
515+
csccolptr[1] = 1
516+
@inbounds for j in 2:(n+1)
517+
overwritten = csccolptr[j]
518+
csccolptr[j] = countsum
519+
countsum += overwritten
520+
end
521+
522+
# Now knowing the CSC form's entry count, resize cscrowval and cscnzval if necessary
523+
cscnnz = countsum - 1
524+
length(cscrowval) < cscnnz && resize!(cscrowval, cscnnz)
525+
length(cscnzval) < cscnnz && resize!(cscnzval, cscnnz)
526+
527+
# Finally counting-sort the row and nonzero values from the CSR form into cscrowval and
528+
# cscnzval. Tracking write positions in csccolptr corrects the column pointers.
529+
@inbounds for i in 1:m
530+
for csrk in csrrowptr[i]:(csrrowptr[i+1]-1)
531+
j = csrcolval[csrk]
532+
x = csrnzval[csrk]
533+
csck = csccolptr[j+1]
534+
csccolptr[j+1] = csck+1
535+
cscrowval[csck] = i
536+
cscnzval[csck] = x
537+
end
538+
end
539+
540+
SparseMatrixCSC(m, n, csccolptr, cscrowval, cscnzval)
541+
end
542+
function sparse!{Tv,Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti},
543+
V::AbstractVector{Tv}, m::Integer, n::Integer, combine, klasttouch::Vector{Ti},
544+
csrrowptr::Vector{Ti}, csrcolval::Vector{Ti}, csrnzval::Vector{Tv},
545+
csccolptr::Vector{Ti} )
546+
sparse!(I, J, V, m, n, combine, klasttouch,
547+
csrrowptr, csrcolval, csrnzval,
548+
csccolptr, Vector{Ti}(), Vector{Tv}() )
549+
end
550+
function sparse!{Tv,Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti},
551+
V::AbstractVector{Tv}, m::Integer, n::Integer, combine, klasttouch::Vector{Ti},
552+
csrrowptr::Vector{Ti}, csrcolval::Vector{Ti}, csrnzval::Vector{Tv} )
553+
sparse!(I, J, V, m, n, combine, klasttouch,
554+
csrrowptr, csrcolval, csrnzval,
555+
Vector{Ti}(n+1), Vector{Ti}(), Vector{Tv}() )
556+
end
339557

340558
dimlub(I) = isempty(I) ? 0 : Int(maximum(I)) #least upper bound on required sparse matrix dimension
341559

‎doc/stdlib/arrays.rst

+4-2
Original file line numberDiff line numberDiff line change
@@ -855,11 +855,13 @@ Sparse Vectors and Matrices
855855
Sparse vectors and matrices largely support the same set of operations as their
856856
dense counterparts. The following functions are specific to sparse arrays.
857857

858-
.. function:: sparse(I,J,V,[m,n,combine])
858+
.. function:: sparse(I, J, V,[ m, n, combine])
859859

860860
.. Docstring generated from Julia source
861861
862-
Create a sparse matrix ``S`` of dimensions ``m x n`` such that ``S[I[k], J[k]] = V[k]``\ . The ``combine`` function is used to combine duplicates. If ``m`` and ``n`` are not specified, they are set to ``maximum(I)`` and ``maximum(J)`` respectively. If the ``combine`` function is not supplied, ``combine`` defaults to ``+`` unless the elements of ``V`` are Booleans in which case ``combine`` defaults to ``|``\ . All elements of ``I`` must satisfy ``1 <= I[k] <= m``\ , and all elements of ``J`` must satisfy ``1 <= J[k] <= n``\ .
862+
Create a sparse matrix ``S`` of dimensions ``m x n`` such that ``S[I[k], J[k]] = V[k]``\ . The ``combine`` function is used to combine duplicates. If ``m`` and ``n`` are not specified, they are set to ``maximum(I)`` and ``maximum(J)`` respectively. If the ``combine`` function is not supplied, ``combine`` defaults to ``+`` unless the elements of ``V`` are Booleans in which case ``combine`` defaults to ``|``\ . All elements of ``I`` must satisfy ``1 <= I[k] <= m``\ , and all elements of ``J`` must satisfy ``1 <= J[k] <= n``\ . Numerical zeros in (``I``\ , ``J``\ , ``V``\ ) are retained as structural nonzeros.
863+
864+
For additional documentation and an expert driver, see ``Base.SparseArrays.sparse!``\ .
863865

864866
.. function:: sparsevec(I, V, [m, combine])
865867

0 commit comments

Comments
 (0)