11# This file is a part of Julia. License is MIT: https://julialang.org/license
22
33# Compressed sparse columns data structure
4- # Assumes that no zeros are stored in the data structure
4+ # No assumptions about stored zeros in the data structure
55# Assumes that row values in rowval for each column are sorted
66# issorted(rowval[colptr[i]:(colptr[i+1]-1)]) == true
7+ # Assumes that 1 <= colptr[i] <= colptr[i+1] for i in 1..n
8+ # Assumes that nnz <= length(rowval) < typemax(Ti)
9+ # Assumes that 0 <= length(nzval) < typemax(Ti)
710
811"""
912 SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
@@ -20,8 +23,8 @@ struct SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
2023 rowval:: Vector{Ti} # Row indices of stored values
2124 nzval:: Vector{Tv} # Stored values, typically nonzeros
2225
23- function SparseMatrixCSC {Tv,Ti} (m:: Integer , n:: Integer , colptr:: Vector{Ti} , rowval :: Vector{Ti} ,
24- nzval:: Vector{Tv} ) where {Tv,Ti<: Integer }
26+ function SparseMatrixCSC {Tv,Ti} (m:: Integer , n:: Integer , colptr:: Vector{Ti} ,
27+ rowval :: Vector{Ti} , nzval:: Vector{Tv} ) where {Tv,Ti<: Integer }
2528 @noinline throwsz (str, lbl, k) =
2629 throw (ArgumentError (" number of $str ($lbl ) must be ≥ 0, got $k " ))
2730 m < 0 && throwsz (" rows" , ' m' , m)
3235function SparseMatrixCSC (m:: Integer , n:: Integer , colptr:: Vector , rowval:: Vector , nzval:: Vector )
3336 Tv = eltype (nzval)
3437 Ti = promote_type (eltype (colptr), eltype (rowval))
38+ sparse_check_Ti (m, n, Ti)
39+ sparse_check (n, colptr, rowval, nzval)
40+ # silently shorten rowval and nzval to usable index positions.
41+ maxlen = abs (widemul (m, n))
42+ isbitstype (Ti) && (maxlen = min (maxlen, typemax (Ti) - 1 ))
43+ length (rowval) > maxlen && resize! (rowval, maxlen)
44+ length (nzval) > maxlen && resize! (nzval, maxlen)
3545 SparseMatrixCSC {Tv,Ti} (m, n, colptr, rowval, nzval)
3646end
3747
48+ function sparse_check_Ti (m:: Integer , n:: Integer , Ti:: Type )
49+ @noinline throwTi (str, lbl, k) =
50+ throw (ArgumentError (" $str ($lbl = $k ) does not fit in Ti = $(Ti) " ))
51+ 0 ≤ m && (! isbitstype (Ti) || m ≤ typemax (Ti)) || throwTi (" number of rows" , " m" , m)
52+ 0 ≤ n && (! isbitstype (Ti) || n ≤ typemax (Ti)) || throwTi (" number of columns" , " n" , n)
53+ end
54+
55+ function sparse_check (n:: Integer , colptr:: Vector{Ti} , rowval, nzval) where Ti
56+ sparse_check_length (" colptr" , colptr, n+ 1 , String) # don't check upper bound
57+ ckp = Ti (1 )
58+ ckp == colptr[1 ] || throw (ArgumentError (" $ckp == colptr[1] != 1" ))
59+ @inbounds for k = 2 : n+ 1
60+ ck = colptr[k]
61+ ckp <= ck || throw (ArgumentError (" $ckp == colptr[$(k- 1 ) ] > colptr[$k ] == $ck " ))
62+ ckp = ck
63+ end
64+ sparse_check_length (" rowval" , rowval, ckp- 1 , Ti)
65+ sparse_check_length (" nzval" , nzval, 0 , Ti) # we allow empty nzval !!!
66+ end
67+ function sparse_check_length (rowstr, rowval, minlen, Ti)
68+ len = length (rowval)
69+ len >= minlen || throw (ArgumentError (" $len == length($rowstr ) < $minlen " ))
70+ ! isbitstype (Ti) || len < typemax (Ti) ||
71+ throw (ArgumentError (" $len == length($rowstr ) >= $(typemax (Ti)) " ))
72+ end
73+
3874size (S:: SparseMatrixCSC ) = (S. m, S. n)
3975
4076# Define an alias for views of a SparseMatrixCSC which include all rows and a unit range of the columns.
@@ -497,7 +533,10 @@ function sparse(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{
497533 " length(I) (=$(length (I)) ) == length(J) (= $(length (J)) ) == length(V) (= " ,
498534 " $(length (V)) )" )))
499535 end
500-
536+ Tj = Ti
537+ while isbitstype (Tj) && coolen >= typemax (Tj)
538+ Tj = widen (Tj)
539+ end
501540 if m == 0 || n == 0 || coolen == 0
502541 if coolen != 0
503542 if n == 0
@@ -509,13 +548,13 @@ function sparse(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{
509548 SparseMatrixCSC (m, n, fill (one (Ti), n+ 1 ), Vector {Ti} (), Vector {Tv} ())
510549 else
511550 # Allocate storage for CSR form
512- csrrowptr = Vector {Ti } (undef, m+ 1 )
551+ csrrowptr = Vector {Tj } (undef, m+ 1 )
513552 csrcolval = Vector {Ti} (undef, coolen)
514553 csrnzval = Vector {Tv} (undef, coolen)
515554
516555 # Allocate storage for the CSC form's column pointers and a necessary workspace
517556 csccolptr = Vector {Ti} (undef, n+ 1 )
518- klasttouch = Vector {Ti } (undef, n)
557+ klasttouch = Vector {Tj } (undef, n)
519558
520559 # Allocate empty arrays for the CSC form's row and nonzero value arrays
521560 # The parent method called below automagically resizes these arrays
@@ -580,25 +619,28 @@ transposition," ACM TOMS 4(3), 250-269 (1978) inspired this method's use of a pa
580619counting sorts.
581620"""
582621function sparse! (I:: AbstractVector{Ti} , J:: AbstractVector{Ti} ,
583- V:: AbstractVector{Tv} , m:: Integer , n:: Integer , combine, klasttouch:: Vector{Ti } ,
584- csrrowptr:: Vector{Ti } , csrcolval:: Vector{Ti} , csrnzval:: Vector{Tv} ,
585- csccolptr:: Vector{Ti} , cscrowval:: Vector{Ti} , cscnzval:: Vector{Tv} ) where {Tv,Ti<: Integer }
622+ V:: AbstractVector{Tv} , m:: Integer , n:: Integer , combine, klasttouch:: Vector{Tj } ,
623+ csrrowptr:: Vector{Tj } , csrcolval:: Vector{Ti} , csrnzval:: Vector{Tv} ,
624+ csccolptr:: Vector{Ti} , cscrowval:: Vector{Ti} , cscnzval:: Vector{Tv} ) where {Tv,Ti<: Integer ,Tj <: Integer }
586625
587626 require_one_based_indexing (I, J, V)
627+ sparse_check_Ti (m, n, Ti)
628+ sparse_check_length (" I" , I, 0 , Tj)
588629 # Compute the CSR form's row counts and store them shifted forward by one in csrrowptr
589- fill! (csrrowptr, Ti (0 ))
630+ fill! (csrrowptr, Tj (0 ))
590631 coolen = length (I)
632+ min (length (J), length (V)) >= coolen || throw (ArgumentError (" J and V need length >= length(I) = $coolen " ))
591633 @inbounds for k in 1 : coolen
592634 Ik = I[k]
593635 if 1 > Ik || m < Ik
594636 throw (ArgumentError (" row indices I[k] must satisfy 1 <= I[k] <= m" ))
595637 end
596- csrrowptr[Ik+ 1 ] += Ti (1 )
638+ csrrowptr[Ik+ 1 ] += Tj (1 )
597639 end
598640
599641 # Compute the CSR form's rowptrs and store them shifted forward by one in csrrowptr
600- countsum = Ti (1 )
601- csrrowptr[1 ] = Ti (1 )
642+ countsum = Tj (1 )
643+ csrrowptr[1 ] = Tj (1 )
602644 @inbounds for i in 2 : (m+ 1 )
603645 overwritten = csrrowptr[i]
604646 csrrowptr[i] = countsum
@@ -613,7 +655,8 @@ function sparse!(I::AbstractVector{Ti}, J::AbstractVector{Ti},
613655 throw (ArgumentError (" column indices J[k] must satisfy 1 <= J[k] <= n" ))
614656 end
615657 csrk = csrrowptr[Ik+ 1 ]
616- csrrowptr[Ik+ 1 ] = csrk + Ti (1 )
658+ @assert csrk >= Tj (1 ) " index into csrcolval exceeds typemax(Ti)"
659+ csrrowptr[Ik+ 1 ] = csrk + Tj (1 )
617660 csrcolval[csrk] = Jk
618661 csrnzval[csrk] = V[k]
619662 end
@@ -626,21 +669,21 @@ function sparse!(I::AbstractVector{Ti}, J::AbstractVector{Ti},
626669 # Minimizing extraneous communication and nonlocality of reference, primarily by using
627670 # only a single auxiliary array in this step, is the key to this method's performance.
628671 fill! (csccolptr, Ti (0 ))
629- fill! (klasttouch, Ti (0 ))
630- writek = Ti (1 )
672+ fill! (klasttouch, Tj (0 ))
673+ writek = Tj (1 )
631674 newcsrrowptri = Ti (1 )
632- origcsrrowptri = Ti (1 )
675+ origcsrrowptri = Tj (1 )
633676 origcsrrowptrip1 = csrrowptr[2 ]
634677 @inbounds for i in 1 : m
635- for readk in origcsrrowptri: (origcsrrowptrip1- Ti (1 ))
678+ for readk in origcsrrowptri: (origcsrrowptrip1- Tj (1 ))
636679 j = csrcolval[readk]
637680 if klasttouch[j] < newcsrrowptri
638681 klasttouch[j] = writek
639682 if writek != readk
640683 csrcolval[writek] = j
641684 csrnzval[writek] = csrnzval[readk]
642685 end
643- writek += Ti (1 )
686+ writek += Tj (1 )
644687 csccolptr[j+ 1 ] += Ti (1 )
645688 else
646689 klt = klasttouch[j]
@@ -654,23 +697,24 @@ function sparse!(I::AbstractVector{Ti}, J::AbstractVector{Ti},
654697 end
655698
656699 # Compute the CSC form's colptrs and store them shifted forward by one in csccolptr
657- countsum = Ti (1 )
700+ countsum = Tj (1 )
658701 csccolptr[1 ] = Ti (1 )
659702 @inbounds for j in 2 : (n+ 1 )
660703 overwritten = csccolptr[j]
661704 csccolptr[j] = countsum
662705 countsum += overwritten
706+ countsum <= typemax (Ti) || throw (ArgumentError (" more than typemax(Ti)-1 == $(typemax (Ti)- 1 ) entries" ))
663707 end
664708
665709 # Now knowing the CSC form's entry count, resize cscrowval and cscnzval if necessary
666- cscnnz = countsum - Ti (1 )
710+ cscnnz = countsum - Tj (1 )
667711 length (cscrowval) < cscnnz && resize! (cscrowval, cscnnz)
668712 length (cscnzval) < cscnnz && resize! (cscnzval, cscnnz)
669713
670714 # Finally counting-sort the row and nonzero values from the CSR form into cscrowval and
671715 # cscnzval. Tracking write positions in csccolptr corrects the column pointers.
672716 @inbounds for i in 1 : m
673- for csrk in csrrowptr[i]: (csrrowptr[i+ 1 ]- Ti (1 ))
717+ for csrk in csrrowptr[i]: (csrrowptr[i+ 1 ]- Tj (1 ))
674718 j = csrcolval[csrk]
675719 x = csrnzval[csrk]
676720 csck = csccolptr[j+ 1 ]
@@ -683,16 +727,16 @@ function sparse!(I::AbstractVector{Ti}, J::AbstractVector{Ti},
683727 SparseMatrixCSC (m, n, csccolptr, cscrowval, cscnzval)
684728end
685729function sparse! (I:: AbstractVector{Ti} , J:: AbstractVector{Ti} ,
686- V:: AbstractVector{Tv} , m:: Integer , n:: Integer , combine, klasttouch:: Vector{Ti } ,
687- csrrowptr:: Vector{Ti } , csrcolval:: Vector{Ti} , csrnzval:: Vector{Tv} ,
688- csccolptr:: Vector{Ti} ) where {Tv,Ti<: Integer }
730+ V:: AbstractVector{Tv} , m:: Integer , n:: Integer , combine, klasttouch:: Vector{Tj } ,
731+ csrrowptr:: Vector{Tj } , csrcolval:: Vector{Ti} , csrnzval:: Vector{Tv} ,
732+ csccolptr:: Vector{Ti} ) where {Tv,Ti<: Integer ,Tj <: Integer }
689733 sparse! (I, J, V, m, n, combine, klasttouch,
690734 csrrowptr, csrcolval, csrnzval,
691735 csccolptr, Vector {Ti} (), Vector {Tv} ())
692736end
693737function sparse! (I:: AbstractVector{Ti} , J:: AbstractVector{Ti} ,
694- V:: AbstractVector{Tv} , m:: Integer , n:: Integer , combine, klasttouch:: Vector{Ti } ,
695- csrrowptr:: Vector{Ti } , csrcolval:: Vector{Ti} , csrnzval:: Vector{Tv} ) where {Tv,Ti<: Integer }
738+ V:: AbstractVector{Tv} , m:: Integer , n:: Integer , combine, klasttouch:: Vector{Tj } ,
739+ csrrowptr:: Vector{Tj } , csrcolval:: Vector{Ti} , csrnzval:: Vector{Tv} ) where {Tv,Ti<: Integer ,Tj <: Integer }
696740 sparse! (I, J, V, m, n, combine, klasttouch,
697741 csrrowptr, csrcolval, csrnzval,
698742 Vector {Ti} (undef, n+ 1 ), Vector {Ti} (), Vector {Tv} ())
@@ -826,7 +870,7 @@ adjoint!(X::SparseMatrixCSC{Tv,Ti}, A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} = f
826870
827871function ftranspose (A:: SparseMatrixCSC{Tv,Ti} , f:: Function ) where {Tv,Ti}
828872 X = SparseMatrixCSC (A. n, A. m,
829- Vector {Ti} (undef , A. m+ 1 ),
873+ ones (Ti , A. m+ 1 ),
830874 Vector {Ti} (undef, nnz (A)),
831875 Vector {Tv} (undef, nnz (A)))
832876 halfperm! (X, A, 1 : A. n, f)
@@ -1045,7 +1089,7 @@ function permute!(X::SparseMatrixCSC{Tv,Ti}, A::SparseMatrixCSC{Tv,Ti},
10451089 _checkargs_sourcecompatdest_permute! (A, X)
10461090 _checkargs_sourcecompatperms_permute! (A, p, q)
10471091 C = SparseMatrixCSC (A. n, A. m,
1048- Vector {Ti} (undef , A. m + 1 ),
1092+ ones (Ti , A. m + 1 ),
10491093 Vector {Ti} (undef, nnz (A)),
10501094 Vector {Tv} (undef, nnz (A)))
10511095 _checkargs_permutationsvalid_permute! (p, C. colptr, q, X. colptr)
@@ -1064,7 +1108,7 @@ function permute!(A::SparseMatrixCSC{Tv,Ti}, p::AbstractVector{<:Integer},
10641108 q:: AbstractVector{<:Integer} ) where {Tv,Ti}
10651109 _checkargs_sourcecompatperms_permute! (A, p, q)
10661110 C = SparseMatrixCSC (A. n, A. m,
1067- Vector {Ti} (undef , A. m + 1 ),
1111+ ones (Ti , A. m + 1 ),
10681112 Vector {Ti} (undef, nnz (A)),
10691113 Vector {Tv} (undef, nnz (A)))
10701114 workcolptr = Vector {Ti} (undef, A. n + 1 )
@@ -1135,11 +1179,11 @@ function permute(A::SparseMatrixCSC{Tv,Ti}, p::AbstractVector{<:Integer},
11351179 q:: AbstractVector{<:Integer} ) where {Tv,Ti}
11361180 _checkargs_sourcecompatperms_permute! (A, p, q)
11371181 X = SparseMatrixCSC (A. m, A. n,
1138- Vector {Ti} (undef , A. n + 1 ),
1182+ ones (Ti , A. n + 1 ),
11391183 Vector {Ti} (undef, nnz (A)),
11401184 Vector {Tv} (undef, nnz (A)))
11411185 C = SparseMatrixCSC (A. n, A. m,
1142- Vector {Ti} (undef , A. m + 1 ),
1186+ ones (Ti , A. m + 1 ),
11431187 Vector {Ti} (undef, nnz (A)),
11441188 Vector {Tv} (undef, nnz (A)))
11451189 _checkargs_permutationsvalid_permute! (p, C. colptr, q, X. colptr)
@@ -2331,15 +2375,32 @@ function _setindex_scalar!(A::SparseMatrixCSC{Tv,Ti}, _v, _i::Integer, _j::Integ
23312375 # Column j does not contain entry A[i,j]. If v is nonzero, insert entry A[i,j] = v
23322376 # and return. If to the contrary v is zero, then simply return.
23332377 if ! iszero (v)
2334- insert! (A. rowval, searchk, i)
2335- insert! (A. nzval, searchk, v)
2378+ nz = A. colptr[A. n+ 1 ]
2379+ # throw exception before state is partially modified
2380+ ! isbitstype (Ti) || nz < typemax (Ti) ||
2381+ throw (ArgumentError (" nnz(A) going to exceed typemax(Ti) = $(typemax (Ti)) " ))
2382+
2383+ # if nnz(A) < length(rowval/nzval): no need to grow rowval and preserve values
2384+ _insert! (A. rowval, searchk, i, nz)
2385+ _insert! (A. nzval, searchk, v, nz)
23362386 @simd for m in (j + 1 ): (A. n + 1 )
2337- @inbounds A. colptr[m] += 1
2387+ @inbounds A. colptr[m] += Ti ( 1 )
23382388 end
23392389 end
23402390 return A
23412391end
23422392
2393+ # insert item at position pos, shifting only from pos+1 to nz
2394+ function _insert! (v:: Vector , pos:: Integer , item, nz:: Integer )
2395+ if nz > length (v)
2396+ insert! (v, pos, item)
2397+ else # nz < length(v)
2398+ Base. unsafe_copyto! (v, pos+ 1 , v, pos, nz - pos)
2399+ v[pos] = item
2400+ v
2401+ end
2402+ end
2403+
23432404function Base. fill! (V:: SubArray {Tv, <: Any , <: SparseMatrixCSC , Tuple{Vararg{Union{Integer, AbstractVector{<: Integer }},2 }}}, x) where Tv
23442405 A = V. parent
23452406 I, J = V. indices
0 commit comments