Skip to content

Commit 3d80e68

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

File tree

4 files changed

+234
-129
lines changed

4 files changed

+234
-129
lines changed

‎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

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

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

‎doc/stdlib/arrays.rst

+3-1
Original file line numberDiff line numberDiff line change
@@ -855,12 +855,14 @@ 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
862862
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``\ .
863863

864+
For additional documentation and an expert driver, see ``Base.SparseArrays.sparse!``\ .
865+
864866
.. function:: sparsevec(I, V, [m, combine])
865867

866868
.. Docstring generated from Julia source

0 commit comments

Comments
 (0)