diff --git a/base/constants.jl b/base/constants.jl index 77310ba5c4e0d..3ee09c807a3dc 100644 --- a/base/constants.jl +++ b/base/constants.jl @@ -73,7 +73,7 @@ const golden = φ for T in (MathConst, Rational, Integer, Number) ^(::MathConst{:e}, x::T) = exp(x) end -for T in (Range, BitArray, SparseMatrixCSC, StridedArray, AbstractArray) +for T in (Range, BitArray, CompressedSparseMatrix, StridedArray, AbstractArray) .^(::MathConst{:e}, x::T) = exp(x) end ^(::MathConst{:e}, x::AbstractMatrix) = expm(x) diff --git a/base/deprecated.jl b/base/deprecated.jl index c756eae219f12..3d57d9438b8cf 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -413,8 +413,8 @@ eval(Sys, :(@deprecate shlib_list dllist)) @deprecate degrees2radians deg2rad @deprecate radians2degrees rad2deg -@deprecate spzeros(m::Integer) spzeros(m, m) -@deprecate spzeros(Tv::Type, m::Integer) spzeros(Tv, m, m) +@deprecate spzeros(m::Integer; format=CSC) spzeros(m, m; format=format) +@deprecate spzeros(Tv::Type, m::Integer; format=CSC) spzeros(Tv, m, m; format=format) @deprecate myindexes localindexes diff --git a/base/exports.jl b/base/exports.jl index b972eb93dc407..00b9e529cfc5e 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -85,6 +85,9 @@ export SharedMatrix, SharedVector, SparseMatrixCSC, + SparseMatrixCSR, + CSR, + CSC, StatStruct, StepRange, StridedArray, diff --git a/base/sparse.jl b/base/sparse.jl index 24d46e9b0a72f..f27740b6c3814 100644 --- a/base/sparse.jl +++ b/base/sparse.jl @@ -3,14 +3,14 @@ include("sparse/abstractsparse.jl") module SparseMatrix importall Base -import Base.NonTupleType, Base.float, Base.Order, Base.Sort.Forward +import Base.NonTupleType, Base.float, Base.Order, Base.Sort.Forward, Base.showarray, Base.dims2string, Base.with_output_limit -export SparseMatrixCSC, +export SparseMatrixCSC, SparseMatrixCSR, CompressedSparseMatrix, blkdiag, dense, diag, diagm, droptol!, dropzeros!, etree, full, getindex, ishermitian, issparse, issym, istril, istriu, nnz, setindex!, sparse, sparsevec, spdiagm, speye, spones, sprand, sprandbool, sprandn, spzeros, symperm, trace, tril, tril!, - triu, triu!, nonzeros + triu, triu!, nonzeros, writemime, summary, showarray include("sparse/sparsematrix.jl") include("sparse/csparse.jl") diff --git a/base/sparse/abstractsparse.jl b/base/sparse/abstractsparse.jl index 6c6726d678b4e..2b5014b09c728 100644 --- a/base/sparse/abstractsparse.jl +++ b/base/sparse/abstractsparse.jl @@ -7,3 +7,7 @@ issparse(A::AbstractArray) = false issparse(S::AbstractSparseArray) = true indtype{Tv,Ti}(S::AbstractSparseArray{Tv,Ti}) = Ti + +const CSR = :csr +const CSC = :csc + diff --git a/base/sparse/csparse.jl b/base/sparse/csparse.jl index e1410689debc5..11e4a4c6fbe22 100644 --- a/base/sparse/csparse.jl +++ b/base/sparse/csparse.jl @@ -12,11 +12,28 @@ # Based on Direct Methods for Sparse Linear Systems, T. A. Davis, SIAM, Philadelphia, Sept. 2006. # Section 2.4: Triplet form # http://www.cise.ufl.edu/research/sparse/CSparse/ -function sparse{Tv,Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, - V::AbstractVector{Tv}, - nrow::Integer, ncol::Integer, combine::Function) - if length(I) == 0; return spzeros(eltype(V),nrow,ncol); end +function sparse{Tv,Ti<:Integer}(nrow::Integer, ncol::Integer, data::TripletSparseStore{Tv,Ti}, combine::Function; format=CSC) + valsptype(format) + I = data.rows + J = data.cols + V = data.vals + + (length(I) == 0) && (return spzeros(eltype(V), nrow, ncol, format=format)) + + if format === CSC + (ncdim, cdim, RpT, RiT, RxT) = sparse_compress(I, J, V, nrow, ncol, combine) + SparseMatrixCSC(ncdim, cdim, RpT, RiT, RxT) + else + (ncdim, cdim, RpT, RiT, RxT) = sparse_compress(J, I, V, ncol, nrow, combine) + SparseMatrixCSR(cdim, ncdim, RpT, RiT, RxT) + end +end +sparse{Tv,Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{Tv}, nrow::Integer, ncol::Integer, combine::Function) = + sparse(nrow, ncol, TripletSparseStore(I, J, V), combine, format=CSC) + +function sparse_compress{Tv,Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{Tv}, + nrow::Integer, ncol::Integer, combine::Function) # Work array Wj = Array(Ti, max(nrow,ncol)+1) @@ -85,7 +102,7 @@ function sparse{Tv,Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, RiT = Array(Ti, anz) RxT = Array(Tv, anz) - # Reset work array to build the final colptr + # Reset work array to build the final indptr Wj[1] = 1 for i=2:(ncol+1); Wj[i] = 0; end for j = 1:nrow @@ -111,7 +128,7 @@ function sparse{Tv,Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, end end - return SparseMatrixCSC(nrow, ncol, RpT, RiT, RxT) + (nrow, ncol, RpT, RiT, RxT) end ## Transpose @@ -120,91 +137,86 @@ end # Section 2.5: Transpose # http://www.cise.ufl.edu/research/sparse/CSparse/ -# NOTE: When calling transpose!(S,T), the colptr in the result matrix T must be set up +# NOTE: When calling transpose!(S,T), the indptr in the result matrix T must be set up # by counting the nonzeros in every row of S. -function transpose!{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, T::SparseMatrixCSC{Tv,Ti}) - (mT, nT) = size(T) - colptr_T = T.colptr - rowval_T = T.rowval +function transpose!{Tc<:CompressedSparseMatrix}(S::Tc, T::Tc) + indptr_T = getindptr(T) + indval_T = getindval(T) nzval_T = T.nzval nnzS = nnz(S) - colptr_S = S.colptr - rowval_S = S.rowval + indptr_S = getindptr(S) + indval_S = getindval(S) nzval_S = S.nzval - w = copy(colptr_T) - @inbounds for j = 1:mT, p = colptr_S[j]:(colptr_S[j+1]-1) - ind = rowval_S[p] + w = copy(indptr_T) + @inbounds for j = 1:nspdim(T), p = indptr_S[j]:(indptr_S[j+1]-1) + ind = indval_S[p] q = w[ind] w[ind] += 1 - rowval_T[q] = j + indval_T[q] = j nzval_T[q] = nzval_S[p] end - return SparseMatrixCSC(mT, nT, colptr_T, rowval_T, nzval_T) + return sparse(size(T,1), size(T,2), CompressedSparseStore(indptr_T, indval_T, nzval_T), format=sptype(S)) end -function transpose{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}) - (nT, mT) = size(S) +function transpose{Tv,Ti}(S::CompressedSparseMatrix{Tv,Ti}) nnzS = nnz(S) - rowval_S = S.rowval + indval_S = getindval(S) - rowval_T = Array(Ti, nnzS) + indval_T = Array(Ti, nnzS) nzval_T = Array(Tv, nnzS) - colptr_T = zeros(Ti, nT+1) - colptr_T[1] = 1 + indptr_T = zeros(Ti, nspdim(S)+1) + indptr_T[1] = 1 @inbounds for i=1:nnz(S) - colptr_T[rowval_S[i]+1] += 1 + indptr_T[indval_S[i]+1] += 1 end - colptr_T = cumsum(colptr_T) + indptr_T = cumsum(indptr_T) - T = SparseMatrixCSC(mT, nT, colptr_T, rowval_T, nzval_T) + T = sparse(size(S,2), size(S,1), CompressedSparseStore(indptr_T, indval_T, nzval_T), format=sptype(S)) return transpose!(S, T) end -# NOTE: When calling ctranspose!(S,T), the colptr in the result matrix T must be set up +# NOTE: When calling ctranspose!(S,T), the indptr in the result matrix T must be set up # by counting the nonzeros in every row of S. -function ctranspose!{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, T::SparseMatrixCSC{Tv,Ti}) - (mT, nT) = size(T) - colptr_T = T.colptr - rowval_T = T.rowval +function ctranspose!{Tc<:CompressedSparseMatrix}(S::Tc, T::Tc) + indptr_T = getindptr(T) + indval_T = getindval(T) nzval_T = T.nzval nnzS = nnz(S) - colptr_S = S.colptr - rowval_S = S.rowval + indptr_S = getindptr(S) + indval_S = getindval(S) nzval_S = S.nzval - w = copy(colptr_T) - @inbounds for j = 1:mT, p = colptr_S[j]:(colptr_S[j+1]-1) - ind = rowval_S[p] + w = copy(indptr_T) + @inbounds for j = 1:nspdim(T), p = indptr_S[j]:(indptr_S[j+1]-1) + ind = indval_S[p] q = w[ind] w[ind] += 1 - rowval_T[q] = j + indval_T[q] = j nzval_T[q] = conj(nzval_S[p]) end - - return SparseMatrixCSC(mT, nT, colptr_T, rowval_T, nzval_T) + return sparse(size(T,1), size(T,2), CompressedSparseStore(indptr_T, indval_T, nzval_T), format=sptype(S)) end -function ctranspose{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}) - (nT, mT) = size(S) +function ctranspose{Tv,Ti}(S::CompressedSparseMatrix{Tv,Ti}) nnzS = nnz(S) - rowval_S = S.rowval + indval_S = getindval(S) - rowval_T = Array(Ti, nnzS) + indval_T = Array(Ti, nnzS) nzval_T = Array(Tv, nnzS) - colptr_T = zeros(Ti, nT+1) - colptr_T[1] = 1 + indptr_T = zeros(Ti, nspdim(S)+1) + indptr_T[1] = 1 @inbounds for i=1:nnz(S) - colptr_T[rowval_S[i]+1] += 1 + indptr_T[indval_S[i]+1] += 1 end - colptr_T = cumsum(colptr_T) + indptr_T = cumsum(indptr_T) - T = SparseMatrixCSC(mT, nT, colptr_T, rowval_T, nzval_T) + T = sparse(size(S,2), size(S,1), CompressedSparseStore(indptr_T, indval_T, nzval_T), format=sptype(S)) return ctranspose!(S, T) end @@ -214,8 +226,8 @@ end # A trivial example is speye(n, n) in which every node is a root. function etree{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, postorder::Bool) m,n = size(A) - Ap = A.colptr - Ai = A.rowval + Ap = getindptr(A) + Ai = getindval(A) parent = zeros(Ti, n) ancestor = zeros(Ti, n) for k in 1:n, p in Ap[k]:(Ap[k+1] - 1) @@ -263,7 +275,7 @@ etree(A::SparseMatrixCSC) = etree(A, false) # find nonzero pattern of Cholesky L[k,1:k-1] using etree and triu(A[:,k]) # based on cs_ereach p. 43, "Direct Methods for Sparse Linear Systems" function ereach{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, k::Integer, parent::Vector{Ti}) - m,n = size(A); Ap = A.colptr; Ai = A.rowval + m,n = size(A); Ap = getindptr(A); Ai = getindval(A) s = Ti[]; sizehint(s, n) # to be used as a stack visited = falses(n) visited[k] = true @@ -281,12 +293,12 @@ end # based on cs_permute p. 21, "Direct Methods for Sparse Linear Systems" function csc_permute{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti}, q::Vector{Ti}) - m,n = size(A); Ap = A.colptr; Ai = A.rowval; Ax = A.nzval + m,n = size(A); Ap = getindptr(A); Ai = getindval(A); Ax = A.nzval if length(pinv) != m || length(q) != n error("dimension mismatch, size(A) = $(size(A)), length(pinv) = $(length(pinv)) and length(q) = $(length(q))") end if !isperm(pinv) || !isperm(q) error("both pinv and q must be permutations") end - C = copy(A); Cp = C.colptr; Ci = C.rowval; Cx = C.nzval + C = copy(A); Cp = getindptr(C); Ci = getindval(C); Cx = C.nzval nz = zero(Ti) for k in 1:n Cp[k] = nz @@ -304,10 +316,10 @@ end # based on cs_symperm p. 21, "Direct Methods for Sparse Linear Systems" # form A[p,p] for a symmetric A stored in the upper triangle function symperm{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti}) - m,n = size(A); Ap = A.colptr; Ai = A.rowval; Ax = A.nzval + m,n = size(A); Ap = getindptr(A); Ai = getindval(A); Ax = A.nzval isperm(pinv) || error("perm must be a permutation") m == n == length(pinv) || error("dimension mismatch") - C = copy(A); Cp = C.colptr; Ci = C.rowval; Cx = C.nzval + C = copy(A); Cp = getindptr(C); Ci = getindval(C); Cx = C.nzval w = zeros(Ti,n) for j in 1:n # count entries in each column of C j2 = pinv[j] @@ -337,23 +349,26 @@ end function fkeep!{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, f, other) nzorig = nnz(A) nz = 1 + indptrA = getindptr(A) + indvalA = getindval(A) + nzvalA = A.nzval for j = 1:A.n - p = A.colptr[j] # record current position - A.colptr[j] = nz # set new position - while p < A.colptr[j+1] - if f(A.rowval[p], j, A.nzval[p], other) - A.nzval[nz] = A.nzval[p] - A.rowval[nz] = A.rowval[p] + p = indptrA[j] # record current position + indptrA[j] = nz # set new position + while p < indptrA[j+1] + if f(indvalA[p], j, nzvalA[p], other) + nzvalA[nz] = nzvalA[p] + indvalA[nz] = indvalA[p] nz += 1 end p += 1 end end - A.colptr[A.n + 1] = nz + indptrA[A.n + 1] = nz nz -= 1 if nz < nzorig - resize!(A.nzval, nz) - resize!(A.rowval, nz) + resize!(nzvalA, nz) + resize!(indvalA, nz) end A end diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index cba0df7bb0a92..9df155a125218 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -10,34 +10,89 @@ type SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti} rowval::Vector{Ti} # Row values of nonzeros nzval::Vector{Tv} # Nonzero values end +type SparseMatrixCSR{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti} + m::Int # Number of rows + n::Int # Number of columns + rowptr::Vector{Ti} # Row i is in rowptr[i]:(rowptr[i+1]-1) + colval::Vector{Ti} # Column values of nonzeros + nzval::Vector{Tv} # Nonzero values +end +typealias CompressedSparseMatrix{Tv,Ti} Union(SparseMatrixCSC{Tv,Ti}, SparseMatrixCSR{Tv,Ti}) + +immutable CompressedSparseStore{Tv,Ti} + indptr::Vector{Ti} + indval::Vector{Ti} + nzval::Vector{Tv} +end + +immutable TripletSparseStore{Tv,Ti} + rows::AbstractVector{Ti} + cols::AbstractVector{Ti} + vals::AbstractVector{Tv} +end + +SparseMatrixCSC{Tv,Ti<:Integer}(m::Integer, n::Integer, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}) = + SparseMatrixCSC{Tv,Ti}(int(m), int(n), colptr, rowval, nzval) + +SparseMatrixCSR{Tv,Ti<:Integer}(m::Integer, n::Integer, rowptr::Vector{Ti}, colval::Vector{Ti}, nzval::Vector{Tv}) = + SparseMatrixCSR{Tv,Ti}(int(m), int(n), rowptr, colval, nzval) + +sptype{T<:CompressedSparseMatrix}(::Type{T}) = (T <: SparseMatrixCSC) ? CSC : CSR +sptype{T<:CompressedSparseMatrix}(::T) = sptype(T) +types{Tv,Ti}(::CompressedSparseMatrix{Tv,Ti}) = (Tv,Ti) +types{T<:CompressedSparseMatrix}(::Type{T}) = types(T.types[1]) +types{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}) = (Tv,Ti) +types{Tv,Ti}(::Type{SparseMatrixCSR{Tv,Ti}}) = (Tv,Ti) -SparseMatrixCSC{Tv,Ti}(m::Integer, n::Integer, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}) = - SparseMatrixCSC(int(m), int(n), colptr, rowval, nzval) +spaxis{T<:CompressedSparseMatrix}(::T) = (T <: SparseMatrixCSC) ? 2 : 1 +nspaxis{T<:CompressedSparseMatrix}(S::T) = spaxis(S) $ 3 +spdim(S::CompressedSparseMatrix) = size(S, spaxis(S)) +nspdim(S::CompressedSparseMatrix) = size(S, nspaxis(S)) -size(S::SparseMatrixCSC) = (S.m, S.n) -nnz(S::SparseMatrixCSC) = int(S.colptr[end]-1) -countnz(S::SparseMatrixCSC) = countnz(S.nzval) +getindptr{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}) = S.colptr +getindval{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}) = S.rowval +setindptr!{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, colptr::Vector{Ti}) = (S.colptr = colptr) +setindval!{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, rowval::Vector{Ti}) = (S.rowval = rowval) -function Base.showarray(io::IO, S::SparseMatrixCSC; +getindptr{Tv,Ti}(S::SparseMatrixCSR{Tv,Ti}) = S.rowptr +getindval{Tv,Ti}(S::SparseMatrixCSR{Tv,Ti}) = S.colval +setindptr!{Tv,Ti}(S::SparseMatrixCSR{Tv,Ti}, rowptr::Vector{Ti}) = (S.rowptr = rowptr) +setindval!{Tv,Ti}(S::SparseMatrixCSR{Tv,Ti}, colval::Vector{Ti}) = (S.colval = colval) + +valsptype(format) = (format === CSC) || (format === CSR) || error("unknown sparse matrix format $format") + +function sparse{Tv,Ti}(m::Integer, n::Integer, data::CompressedSparseStore{Tv,Ti}; format=CSC) + valsptype(format) + ((format === CSC) ? SparseMatrixCSC : SparseMatrixCSR)(m, n, data.indptr, data.indval, data.nzval) +end + +size(S::CompressedSparseMatrix) = (S.m, S.n) +size(S::CompressedSparseMatrix, d::Integer) = (d == 1) ? S.m : (d == 2) ? S.n : 1 +nnz(S::CompressedSparseMatrix) = int(getindptr(S)[end]-1) +countnz(S::CompressedSparseMatrix) = countnz(S.nzval) + +summary(a::CompressedSparseMatrix) = string(dims2string(size(a)), " $(sptype(a)) Sparse Matrix with $(nnz(a)) $(eltype(a)) entries:") +writemime(io::IO, ::MIME"text/plain", v::CompressedSparseMatrix) = with_output_limit(()->showarray(io, v, header=true, repr=false)) + +function showarray(io::IO, S::CompressedSparseMatrix; header::Bool=true, limit::Bool=Base._limit_output, rows = Base.tty_size()[1], repr=false) # TODO: repr? + header && print(io, summary(S)) - if header - print(io, S.m, "x", S.n, " sparse matrix with ", nnz(S), " ", eltype(S), " entries:") - end - - if limit - half_screen_rows = div(rows - 8, 2) - else - half_screen_rows = typemax(Int) - end + half_screen_rows = limit ? div(rows - 8, 2) : typemax(Int) pad = ndigits(max(S.m,S.n)) k = 0 sep = "\n\t" - for col = 1:S.n, k = S.colptr[col] : (S.colptr[col+1]-1) + + caxis = spaxis(S) + ncaxis = nspaxis(S) + rowcol = Array(Int,2) + for cidx = 1:spdim(S), k = getindptr(S)[cidx] : (getindptr(S)[cidx+1]-1) if k < half_screen_rows || k > nnz(S)-half_screen_rows - print(io, sep, '[', rpad(S.rowval[k], pad), ", ", lpad(col, pad), "] = ") + rowcol[caxis] = cidx + rowcol[ncaxis] = getindval(S)[k] + print(io, sep, '[', rpad(rowcol[1], pad), ", ", lpad(rowcol[2], pad), "] = ") showcompact(io, S.nzval[k]) elseif k == half_screen_rows print(io, sep, '\u22ee') @@ -48,136 +103,147 @@ end ## Reinterpret and Reshape -function reinterpret{T,Tv,Ti}(::Type{T}, a::SparseMatrixCSC{Tv,Ti}) - if sizeof(T) != sizeof(Tv) - throw(ArgumentError("SparseMatrixCSC reinterpret is only supported for element types of the same size")) - end +function reinterpret{T,Tv,Ti}(::Type{T}, a::CompressedSparseMatrix{Tv,Ti}) + (sizeof(T) == sizeof(Tv)) || throw(ArgumentError("SparseMatrix reinterpret is only supported for element types of the same size")) + mA, nA = size(a) - colptr = copy(a.colptr) - rowval = copy(a.rowval) + iptr = copy(getindptr(a)) + ival = copy(getindval(a)) nzval = reinterpret(T, a.nzval) - return SparseMatrixCSC{T,Ti}(mA, nA, colptr, rowval, nzval) + sparse(mA, nA, CompressedSparseStore(iptr, ival, nzval), format=sptype(a)) end -function sparse_compute_reshaped_colptr_and_rowval(colptrS, rowvalS, mS, nS, colptrA, rowvalA, mA, nA) - colptrS[1] = 1 +function sparse_compute_reshaped_indptr_and_indval{Ti}(indptrS::Vector{Ti}, indvalS::Vector{Ti}, mS::Int, nS::Int, indptrA::Vector{Ti}, indvalA::Vector{Ti}, mA::Int, nA::Int) + indptrS[1] = 1 colA = 1 colS = 1 ptr = 1 while colA <= nA - while ptr <= colptrA[colA+1]-1 - rowA = rowvalA[ptr] + while ptr <= indptrA[colA+1]-1 + rowA = indvalA[ptr] i = (colA - 1) * mA + rowA - 1 colSn = div(i, mS) + 1 rowS = mod(i, mS) + 1 while colS < colSn - colptrS[colS+1] = ptr + indptrS[colS+1] = ptr colS += 1 end - rowvalS[ptr] = rowS + indvalS[ptr] = rowS ptr += 1 end colA += 1 end while colS <= nS - colptrS[colS+1] = ptr + indptrS[colS+1] = ptr colS += 1 end end -function reinterpret{T,Tv,Ti,N}(::Type{T}, a::SparseMatrixCSC{Tv,Ti}, dims::NTuple{N,Int}) - if sizeof(T) != sizeof(Tv) - throw(ArgumentError("SparseMatrixCSC reinterpret is only supported for element types of the same size")) - end - if prod(dims) != length(a) - throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(length(a))")) - end - mS,nS = dims - mA,nA = size(a) +function reinterpret{T,Tv,Ti,N}(::Type{T}, a::CompressedSparseMatrix{Tv,Ti}, dims::NTuple{N,Int}) + (sizeof(T) == sizeof(Tv)) || throw(ArgumentError("SparseMatrix reinterpret is only supported for element types of the same size")) + (prod(dims) == length(a)) || throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(length(a))")) + + caxis = spaxis(a) + ncaxis = nspaxis(a) + + mS,nS = (dims[ncaxis], dims[caxis]) + mA,nA = (size(a,ncaxis), size(a,caxis)) numnz = nnz(a) - colptr = Array(Ti, nS+1) - rowval = Array(Ti, numnz) + iptr = Array(Ti, nS+1) + ival = Array(Ti, numnz) nzval = reinterpret(T, a.nzval) - sparse_compute_reshaped_colptr_and_rowval(colptr, rowval, mS, nS, a.colptr, a.rowval, mA, nA) + sparse_compute_reshaped_indptr_and_indval(iptr, ival, mS, nS, getindptr(a), getindval(a), mA, nA) - return SparseMatrixCSC{T,Ti}(mS, nS, colptr, rowval, nzval) + return sparse(dims[1], dims[2], CompressedSparseStore{T,Ti}(iptr, ival, nzval), format=sptype(a)) end -function reshape{Tv,Ti}(a::SparseMatrixCSC{Tv,Ti}, dims::NTuple{2,Int}) - if prod(dims) != length(a) - throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(length(a))")) - end - mS,nS = dims - mA,nA = size(a) +function reshape{Tv,Ti}(a::CompressedSparseMatrix{Tv,Ti}, dims::NTuple{2,Int}) + (prod(dims) == length(a)) || throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(length(a))")) + + caxis = spaxis(a) + ncaxis = nspaxis(a) + + mS,nS = (dims[ncaxis], dims[caxis]) + mA,nA = (size(a,ncaxis), size(a,caxis)) numnz = nnz(a) - colptr = Array(Ti, nS+1) - rowval = Array(Ti, numnz) + iptr = Array(Ti, nS+1) + ival = Array(Ti, numnz) nzval = a.nzval - sparse_compute_reshaped_colptr_and_rowval(colptr, rowval, mS, nS, a.colptr, a.rowval, mA, nA) + sparse_compute_reshaped_indptr_and_indval(iptr, ival, mS, nS, getindptr(a), getindval(a), mA, nA) - return SparseMatrixCSC{Tv,Ti}(mS, nS, colptr, rowval, nzval) + return sparse(dims[1], dims[2], CompressedSparseStore{Tv,Ti}(iptr, ival, nzval), format=sptype(a)) end ## Constructors -copy(S::SparseMatrixCSC) = - SparseMatrixCSC(S.m, S.n, copy(S.colptr), copy(S.rowval), copy(S.nzval)) - -similar(S::SparseMatrixCSC, Tv::NonTupleType=eltype(S)) = SparseMatrixCSC(S.m, S.n, copy(S.colptr), copy(S.rowval), Array(Tv, length(S.nzval))) -similar{Tv,Ti,TvNew}(S::SparseMatrixCSC{Tv,Ti}, ::Type{TvNew}, ::Type{Ti}) = similar(S, TvNew) -similar{Tv,Ti,TvNew,TiNew}(S::SparseMatrixCSC{Tv,Ti}, ::Type{TvNew}, ::Type{TiNew}) = SparseMatrixCSC(S.m, S.n, convert(Array{TiNew},S.colptr), convert(Array{TiNew}, S.rowval), Array(TvNew, length(S.nzval))) -similar{Tv}(S::SparseMatrixCSC, ::Type{Tv}, d::NTuple{Integer}) = spzeros(Tv, d...) - -function convert{Tv,Ti,TvS,TiS}(::Type{SparseMatrixCSC{Tv,Ti}}, S::SparseMatrixCSC{TvS,TiS}) - if Tv == TvS && Ti == TiS - return S - else - return SparseMatrixCSC(S.m, S.n, - convert(Vector{Ti},S.colptr), - convert(Vector{Ti},S.rowval), - convert(Vector{Tv},S.nzval)) +copy(S::CompressedSparseMatrix) = sparse(S.m, S.n, CompressedSparseStore(copy(getindptr(S)), copy(getindval(S)), copy(S.nzval)), format=sptype(S)) + +similar(S::CompressedSparseMatrix, Tv::NonTupleType=eltype(S)) = + sparse(S.m, S.n, CompressedSparseStore(copy(getindptr(S)), copy(getindval(S)), Array(Tv, length(S.nzval))), format=sptype(S)) +similar{Tv,Ti,TvNew}(S::CompressedSparseMatrix{Tv}, ::Type{TvNew}, ::Type{Ti}) = similar(S, TvNew) +similar{Tv,Ti,TvNew,TiNew}(S::CompressedSparseMatrix{Tv,Ti}, ::Type{TvNew}, ::Type{TiNew}) = + sparse(S.m, S.n, CompressedSparseStore(convert(Array{TiNew},getindptr(S)), convert(Array{TiNew}, getindval(S)), Array(TvNew, length(S.nzval))), format=sptype(S)) +similar{Tv}(S::CompressedSparseMatrix, ::Type{Tv}, d::NTuple{Integer}) = spzeros(Tv, d..., format=sptype(S)) + +convert{T<:CompressedSparseMatrix}(::Type{T}, S::T) = S +convert{TvNew,Tv,Ti}(::Type{SparseMatrixCSC{TvNew}}, S::CompressedSparseMatrix{Tv,Ti}) = convert(SparseMatrixCSC{TvNew,Ti}, S) +convert{TvNew,Tv,Ti}(::Type{SparseMatrixCSR{TvNew}}, S::CompressedSparseMatrix{Tv,Ti}) = convert(SparseMatrixCSR{TvNew,Ti}, S) +function convert{T1<:CompressedSparseMatrix,T2<:CompressedSparseMatrix}(::Type{T1}, S::T2) + if sptype(T1) != sptype(T2) + S = S' + S = sparse(S.m, S.n, CompressedSparseStore(getindptr(S), getindval(S), S.nzval), format=sptype(T1)) end -end + (Tv,Ti) = types(T2) + (TvNew,TiNew) = types(T1) + (Tv == TvNew) && (Ti == TiNew) && (return S) -function convert{Tv,TvS,TiS}(::Type{SparseMatrixCSC{Tv}}, S::SparseMatrixCSC{TvS,TiS}) - if Tv == TvS - return S + if Ti === TiNew + iptr = getindptr(S) + ival = getindval(S) else - return SparseMatrixCSC(S.m, S.n, - S.colptr, - S.rowval, - convert(Vector{Tv},S.nzval)) + iptr = convert(Vector{TiNew}, getindptr(S)) + ival = convert(Vector{TiNew}, getindval(S)) end -end + nzval = (Tv === TvNew) ? S.nzval : convert(Vector{TvNew}, S.nzval) -function convert{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}, M::Matrix) + sparse(S.m, S.n, CompressedSparseStore(iptr, ival, nzval), format=sptype(S)) +end +function convert{T<:CompressedSparseMatrix}(::Type{T}, M::Matrix) m, n = size(M) (I, J, V) = findnz(M) + (Tv,Ti) = types(T) return sparse_IJ_sorted!(convert(Vector{Ti},I), convert(Vector{Ti},J), convert(Vector{Tv},V), - m, n) + m, n, sptype(T)) end -convert{T}(::Type{AbstractMatrix{T}}, A::SparseMatrixCSC) = convert(SparseMatrixCSC{T}, A) -convert(::Type{Matrix}, S::SparseMatrixCSC) = full(S) +convert(::Type{Matrix}, S::CompressedSparseMatrix) = full(S) -function full{Tv}(S::SparseMatrixCSC{Tv}) +function full{Tv}(S::CompressedSparseMatrix{Tv}) A = zeros(Tv, S.m, S.n) - for col = 1 : S.n, k = S.colptr[col] : (S.colptr[col+1]-1) - A[S.rowval[k], col] = S.nzval[k] + + caxis = spaxis(S) + cdim = spdim(S) + if sptype(S) === CSC + for cidx = 1:cdim, k = getindptr(S)[cidx] : (getindptr(S)[cidx+1]-1) + A[getindval(S)[k], cidx] = S.nzval[k] + end + else + for cidx = 1:cdim, k = getindptr(S)[cidx] : (getindptr(S)[cidx+1]-1) + A[cidx, getindval(S)[k]] = S.nzval[k] + end end return A end -float(S::SparseMatrixCSC) = SparseMatrixCSC(S.m, S.n, copy(S.colptr), copy(S.rowval), float(copy(S.nzval))) - -complex(S::SparseMatrixCSC) = SparseMatrixCSC(S.m, S.n, copy(S.colptr), copy(S.rowval), complex(copy(S.nzval))) - +float(S::CompressedSparseMatrix) = sparse(S.m, S.n, CompressedSparseStore(copy(getindptr(S)), copy(getindval(S)), float(copy(S.nzval))), format=sptype(S)) +complex(S::CompressedSparseMatrix) = sparse(S.m, S.n, CompressedSparseStore(copy(getindptr(S)), copy(getindval(S)), complex(copy(S.nzval))), format=sptype(S)) complex(A::SparseMatrixCSC, B::SparseMatrixCSC) = A + im*B +complex(A::SparseMatrixCSR, B::SparseMatrixCSR) = A + im*B # Construct a sparse vector @@ -196,7 +262,7 @@ function sparsevec(I::AbstractVector, V, m::Integer, combine::Function) I = I[p] m >= I[end] || throw(DimensionMismatch("indices cannot be larger than length of vector")) V = V[p] - sparse_IJ_sorted!(I, ones(eltype(I), nI), V, m, 1, combine) + sparse_IJ_sorted!(I, ones(eltype(I), nI), V, m, 1, combine,CSC) end function sparsevec(a::Vector) @@ -204,28 +270,41 @@ function sparsevec(a::Vector) I = find(a) J = ones(Int, n) V = a[I] - return sparse_IJ_sorted!(I,J,V,n,1,+) + return sparse_IJ_sorted!(I,J,V,n,1,+,CSC) end sparse(a::Vector) = sparsevec(a) ## Construct a sparse matrix -sparse{Tv}(A::Matrix{Tv}) = convert(SparseMatrixCSC{Tv,Int}, A) - -sparse(S::SparseMatrixCSC) = copy(S) - -sparse_IJ_sorted!(I,J,V,m,n) = sparse_IJ_sorted!(I,J,V,m,n,+) +function sparse{Tv}(A::Matrix{Tv}; format=CSC) + valsptype(format) + convert((format === CSC) ? SparseMatrixCSC{Tv,Int} : SparseMatrixCSR{Tv,Int}, A) +end -sparse_IJ_sorted!(I,J,V::AbstractVector{Bool},m,n) = sparse_IJ_sorted!(I,J,V,m,n,|) +sparse(S::CompressedSparseMatrix) = copy(S) + +sparse_IJ_sorted!(I,J,V,m,n,format) = sparse_IJ_sorted!(I,J,V,m,n,+,format) +sparse_IJ_sorted!(I,J,V::AbstractVector{Bool},m,n,format) = sparse_IJ_sorted!(I,J,V,m,n,|,format) +# TODO: remove sortperm +function sparse_IJ_sorted!{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector, + m::Integer, n::Integer, combine::Function, format) + (length(V) == 0) && (return spzeros(eltype(V),Ti,m,n, format=format)) + if format === CSC + #s = sortperm(J) + #(iptr, ival, V) = compress_IJ_sorted!(I[s], J[s], V, m, n, combine) + (iptr, ival, V) = compress_IJ_sorted!(I, J, V, m, n, combine) + else + s = sortperm(I) + (iptr, ival, V) = compress_IJ_sorted!(J[s], I[s], V, n, m, combine) + end + return sparse(m, n, CompressedSparseStore(iptr, ival, V), format=format) +end -function sparse_IJ_sorted!{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, - V::AbstractVector, +function compress_IJ_sorted!{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector, m::Integer, n::Integer, combine::Function) - m = m < 0 ? 0 : m n = n < 0 ? 0 : n - if length(V) == 0; return spzeros(eltype(V),Ti,m,n); end cols = zeros(Ti, n+1) cols[1] = 1 # For cumsum purposes @@ -252,7 +331,7 @@ function sparse_IJ_sorted!{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector end end - colptr = cumsum(cols) + iptr = cumsum(cols) # Allow up to 20% slack if ndups > 0.2*length(I) @@ -261,65 +340,66 @@ function sparse_IJ_sorted!{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector V = V[1:numnz] end - return SparseMatrixCSC(m, n, colptr, I, V) + return (iptr, I, V) end ## sparse() can take its inputs in unsorted order (the parent method is now in jlsparse.jl) dimlub(I) = length(I)==0 ? 0 : int(maximum(I)) #least upper bound on required sparse matrix dimension -sparse(I,J,v::Number) = sparse(I, J, fill(v,length(I)), dimlub(I), dimlub(J), +) - -sparse(I,J,V::AbstractVector) = sparse(I, J, V, dimlub(I), dimlub(J), +) - -sparse(I,J,v::Number,m,n) = sparse(I, J, fill(v,length(I)), int(m), int(n), +) - -sparse(I,J,V::AbstractVector,m,n) = sparse(I, J, V, int(m), int(n), +) - -sparse(I,J,V::AbstractVector{Bool},m,n) = sparse(I, J, V, int(m), int(n), |) - -sparse(I,J,v::Number,m,n,combine::Function) = sparse(I, J, fill(v,length(I)), int(m), int(n), combine) - -function find(S::SparseMatrixCSC) +sparse{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, v::Number; format=CSC) = + sparse(dimlub(I), dimlub(J), TripletSparseStore(I,J,fill(v,length(I))), +, format=format) +sparse{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector; format=CSC) = + sparse(dimlub(I), dimlub(J), TripletSparseStore(I, J, V), +, format=format) +sparse{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, v::Number, m, n; format=CSC) = + sparse(int(m), int(n), TripletSparseStore(I, J, fill(v,length(I))), +, format=format) +sparse{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector, m, n; format=CSC) = + sparse(int(m), int(n), TripletSparseStore(I, J, V), +, format=format) +sparse{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{Bool}, m, n; format=CSC) = + sparse(int(m), int(n), TripletSparseStore(I, J, V), |, format=format) +sparse{Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti}, v::Number, m, n, combine::Function; format=CSC) = + sparse(int(m), int(n), TripletSparseStore(I, J, fill(v,length(I))), combine, format=format) + +function find(S::CompressedSparseMatrix) sz = size(S) I, J = findn(S) return sub2ind(sz, I, J) end -function findn{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}) +function findn{Tv,Ti}(S::CompressedSparseMatrix{Tv,Ti}) numnz = nnz(S) I = Array(Ti, numnz) J = Array(Ti, numnz) count = 1 - for col = 1 : S.n, k = S.colptr[col] : (S.colptr[col+1]-1) + for cidx = 1:spdim(S), k = getindptr(S)[cidx] : (getindptr(S)[cidx+1]-1) if S.nzval[k] != 0 - I[count] = S.rowval[k] - J[count] = col + I[count] = getindval(S)[k] + J[count] = cidx count += 1 end end count -= 1 if numnz != count - I = I[1:count] - J = J[1:count] + deleteat!(I, (count+1):numnz) + deleteat!(J, (count+1):numnz) end - return (I, J) + return (sptype(S) === CSC) ? (I, J) : (J, I) end -function findnz{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}) +function findnz{Tv,Ti}(S::CompressedSparseMatrix{Tv,Ti}) numnz = nnz(S) I = Array(Ti, numnz) J = Array(Ti, numnz) V = Array(Tv, numnz) count = 1 - for col = 1 : S.n, k = S.colptr[col] : (S.colptr[col+1]-1) + for cidx = 1:spdim(S), k = getindptr(S)[cidx] : (getindptr(S)[cidx+1]-1) if S.nzval[k] != 0 - I[count] = S.rowval[k] - J[count] = col + I[count] = getindval(S)[k] + J[count] = cidx V[count] = S.nzval[k] count += 1 end @@ -327,22 +407,22 @@ function findnz{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}) count -= 1 if numnz != count - I = I[1:count] - J = J[1:count] - V = V[1:count] + deleteat!(I, (count+1):numnz) + deleteat!(J, (count+1):numnz) + deleteat!(V, (count+1):numnz) end - return (I, J, V) + return (sptype(S) === CSC) ? (I, J, V) : (J, I, V) end -nonzeros(S::SparseMatrixCSC) = S.nzval +nonzeros(S::CompressedSparseMatrix) = S.nzval function sprand{T}(m::Integer, n::Integer, density::FloatingPoint, - rng::Function,::Type{T}=eltype(rng(1))) + rng::Function,::Type{T}=eltype(rng(1)); format=CSC) 0 <= density <= 1 || throw(ArgumentError("$density not in [0,1]")) N = n*m - N == 0 && return spzeros(T,m,n) - N == 1 && return rand() <= density ? sparse(rng(1)) : spzeros(T,1,1) + N == 0 && return spzeros(T,m,n, format=format) + N == 1 && return rand() <= density ? sparse(reshape(rng(1),1,1), format=format) : spzeros(T,1,1, format=format) I, J = Array(Int, 0), Array(Int, 0) # indices of nonzero elements sizehint(I, int(N*density)) @@ -373,41 +453,38 @@ function sprand{T}(m::Integer, n::Integer, density::FloatingPoint, J[i] = j end end - return sparse_IJ_sorted!(I, J, rng(length(I)), m, n, +) # it will never need to combine + return sparse_IJ_sorted!(I, J, rng(length(I)), m, n, +, format) # it will never need to combine end -sprand(m::Integer, n::Integer, density::FloatingPoint) = sprand(m,n,density,rand,Float64) -sprandn(m::Integer, n::Integer, density::FloatingPoint) = sprand(m,n,density,randn,Float64) +sprand(m::Integer, n::Integer, density::FloatingPoint; format=CSC) = sprand(m,n,density,rand,Float64,format=format) +sprandn(m::Integer, n::Integer, density::FloatingPoint; format=CSC) = sprand(m,n,density,randn,Float64,format=format) truebools(n::Integer) = ones(Bool, n) -sprandbool(m::Integer, n::Integer, density::FloatingPoint) = sprand(m,n,density,truebools,Bool) +sprandbool(m::Integer, n::Integer, density::FloatingPoint; format=CSC) = sprand(m,n,density,truebools,Bool,format=format) -spones{T}(S::SparseMatrixCSC{T}) = - SparseMatrixCSC(S.m, S.n, copy(S.colptr), copy(S.rowval), ones(T, S.colptr[end]-1)) +spones{T}(S::CompressedSparseMatrix{T}) = sparse(S.m, S.n, CompressedSparseStore(copy(getindptr(S)), copy(getindval(S)), ones(T, getindptr(S)[end]-1)), format=sptype(S)) -spzeros(m::Integer, n::Integer) = spzeros(Float64, m, n) -spzeros(Tv::Type, m::Integer, n::Integer) = - SparseMatrixCSC(m, n, ones(Int, n+1), Array(Int, 0), Array(Tv, 0)) -spzeros(Tv::Type, Ti::Type, m::Integer, n::Integer) = - SparseMatrixCSC(m, n, ones(Ti, n+1), Array(Ti, 0), Array(Tv, 0)) +spzeros(m::Integer, n::Integer; format=CSC) = spzeros(Float64, m, n, format=format) +spzeros(Tv::Type, m::Integer, n::Integer; format=CSC) = sparse(m, n, CompressedSparseStore(ones(Int, ((format===CSC)?n:m)+1), Array(Int, 0), Array(Tv, 0)), format=format) +spzeros(Tv::Type, Ti::Type, m::Integer, n::Integer; format=CSC) = sparse(m, n, CompressedSparseStore(ones(Ti, ((format===CSC)?n:m)+1), Array(Ti, 0), Array(Tv, 0)), format=format) -speye(n::Integer) = speye(Float64, n) -speye(T::Type, n::Integer) = speye(T, n, n) -speye(m::Integer, n::Integer) = speye(Float64, m, n) -speye{T}(S::SparseMatrixCSC{T}) = speye(T, size(S, 1), size(S, 2)) -eye(S::SparseMatrixCSC) = speye(S) +speye(n::Integer; format=CSC) = speye(Float64, n, format=format) +speye(T::Type, n::Integer; format=CSC) = speye(T, n, n, format=format) +speye(m::Integer, n::Integer; format=CSC) = speye(Float64, m, n, format=format) +speye{T}(S::CompressedSparseMatrix{T}) = speye(T, size(S, 1), size(S, 2), format=sptype(S)) +eye(S::CompressedSparseMatrix) = speye(S) -function speye(T::Type, m::Integer, n::Integer) +function speye(T::Type, m::Integer, n::Integer; format=CSC) x = min(m,n) - rowval = [1:x] - colptr = [rowval, fill(int(x+1), n+1-x)] + ival = [1:x] + iptr = [ival, fill(int(x+1), ((format===CSC)?n:m)+1-x)] nzval = ones(T, x) - return SparseMatrixCSC(m, n, colptr, rowval, nzval) + return sparse(m, n, CompressedSparseStore(iptr, ival, nzval), format=format) end -function one{T}(S::SparseMatrixCSC{T}) +function one{T}(S::CompressedSparseMatrix{T}) m,n = size(S) if m != n; throw(DimensionMismatch("multiplicative identity only defined for square matrices")); end - speye(T, m) + speye(T, m, format=sptype(S)) end ## Unary arithmetic and boolean operators @@ -427,29 +504,30 @@ for (op, restype) in ((:iceil, Int), (:ceil, Nothing), (:asind, Nothing), (:atand, Nothing) ) @eval begin - function ($op){Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}) + function ($op){Tv,Ti}(A::CompressedSparseMatrix{Tv,Ti}) + cdim = spdim(A) nfilledA = nnz(A) - colptrB = Array(Ti, A.n+1) - rowvalB = Array(Ti, nfilledA) + indptrB = Array(Ti, cdim+1) + indvalB = Array(Ti, nfilledA) nzvalB = Array($(restype==Nothing ? (:Tv) : restype), nfilledA) k = 0 # number of additional zeros introduced by op(A) - for i = 1 : A.n - colptrB[i] = A.colptr[i] - k - for j = A.colptr[i] : A.colptr[i+1]-1 + for i = 1 : cdim + indptrB[i] = getindptr(A)[i] - k + for j = getindptr(A)[i] : getindptr(A)[i+1]-1 opAj = $(op)(A.nzval[j]) if opAj == 0 k += 1 else - rowvalB[j - k] = A.rowval[j] + indvalB[j - k] = getindval(A)[j] nzvalB[j - k] = opAj end end end - colptrB[end] = A.colptr[end] - k - deleteat!(rowvalB, colptrB[end]:nfilledA) - deleteat!(nzvalB, colptrB[end]:nfilledA) - return SparseMatrixCSC(A.m, A.n, colptrB, rowvalB, nzvalB) + indptrB[end] = getindptr(A)[end] - k + deleteat!(indvalB, indptrB[end]:nfilledA) + deleteat!(nzvalB, indptrB[end]:nfilledA) + return sparse(A.m, A.n, CompressedSparseStore(indptrB, indvalB, nzvalB), format=sptype(A)) end end # quote @@ -460,7 +538,7 @@ end # macro for op in (:-, :abs, :abs2, :log1p, :expm1) @eval begin - function ($op)(A::SparseMatrixCSC) + function ($op)(A::CompressedSparseMatrix) B = similar(A) nzvalB = B.nzval nzvalA = A.nzval @@ -481,13 +559,20 @@ for op in (:cos, :cosh, :acos, :sec, :csc, :cot, :acot, :sech, :exp, :exp2, :exp10) @eval begin - function ($op){Tv}(A::SparseMatrixCSC{Tv}) + function ($op){Tv}(A::CompressedSparseMatrix{Tv}) B = fill($(op)(zero(Tv)), size(A)) - for col = 1 : A.n - for j = A.colptr[col] : A.colptr[col+1]-1 - row = A.rowval[j] + + caxis = spaxis(A) + ncaxis = nspaxis(A) + cdim = spdim(A) + rowcol = Array(Int,2) + + for cidx = 1 : cdim + rowcol[caxis] = cidx + for j = getindptr(A)[cidx] : getindptr(A)[cidx+1]-1 + rowcol[ncaxis] = getindval(A)[j] nz = A.nzval[j] - B[row,col] = $(op)(nz) + B[rowcol[1],rowcol[2]] = $(op)(nz) end end return B @@ -501,49 +586,49 @@ end for (op, restype) in ( (:+, Nothing), (:-, Nothing), (:.*, Nothing), (:(.<), Bool) ) @eval begin - - function ($op){TvA,TiA,TvB,TiB}(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) + function ($op){T1<:CompressedSparseMatrix, T2<:CompressedSparseMatrix}(A::T1, B::T2) + (TvA,TiA) = types(T1) + (TvB,TiB) = types(T2) Tv = promote_type(TvA, TvB) Ti = promote_type(TiA, TiB) - A = convert(SparseMatrixCSC{Tv,Ti}, A) - B = convert(SparseMatrixCSC{Tv,Ti}, B) - return ($op)(A, B) + T = (sptype(A) === CSC) ? SparseMatrixCSC{Tv,Ti} : SparseMatrixCSR{Tv,Ti} + return ($op)(convert(T,A), convert(T,B)) end - function ($op){Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) - if size(A,1) != size(B,1) || size(A,2) != size(B,2) - throw(DimensionMismatch("")) - end + function ($op){T<:CompressedSparseMatrix}(A::T, B::T) + (size(A) == size(B)) || throw(DimensionMismatch("")) (m, n) = size(A) + cdim = spdim(A) + (Tv,Ti) = types(A) # TODO: Need better method to estimate result space nnzS = nnz(A) + nnz(B) - colptrS = Array(Ti, A.n+1) - rowvalS = Array(Ti, nnzS) + indptrS = Array(Ti, cdim+1) + indvalS = Array(Ti, nnzS) nzvalS = Array($(restype==Nothing ? (:Tv) : restype), nnzS) z = zero(Tv) - colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval - colptrB = B.colptr; rowvalB = B.rowval; nzvalB = B.nzval + indptrA = getindptr(A); indvalA = getindval(A); nzvalA = A.nzval + indptrB = getindptr(B); indvalB = getindval(B); nzvalB = B.nzval ptrS = 1 - colptrS[1] = 1 + indptrS[1] = 1 - for col = 1:n - ptrA::Int = colptrA[col] - stopA::Int = colptrA[col+1] - ptrB::Int = colptrB[col] - stopB::Int = colptrB[col+1] + for cidx = 1:cdim + ptrA::Int = indptrA[cidx] + stopA::Int = indptrA[cidx+1] + ptrB::Int = indptrB[cidx] + stopB::Int = indptrB[cidx+1] while ptrA < stopA && ptrB < stopB - rowA = rowvalA[ptrA] - rowB = rowvalB[ptrB] + rowA = indvalA[ptrA] + rowB = indvalB[ptrB] if rowA < rowB res = ($op)(nzvalA[ptrA], z) if res != z - rowvalS[ptrS] = rowA + indvalS[ptrS] = rowA nzvalS[ptrS] = res ptrS += 1 end @@ -551,7 +636,7 @@ for (op, restype) in ( (:+, Nothing), (:-, Nothing), (:.*, Nothing), elseif rowB < rowA res = ($op)(z, nzvalB[ptrB]) if res != z - rowvalS[ptrS] = rowB + indvalS[ptrS] = rowB nzvalS[ptrS] = res ptrS += 1 end @@ -559,7 +644,7 @@ for (op, restype) in ( (:+, Nothing), (:-, Nothing), (:.*, Nothing), else res = ($op)(nzvalA[ptrA], nzvalB[ptrB]) if res != z - rowvalS[ptrS] = rowA + indvalS[ptrS] = rowA nzvalS[ptrS] = res ptrS += 1 end @@ -571,8 +656,8 @@ for (op, restype) in ( (:+, Nothing), (:-, Nothing), (:.*, Nothing), while ptrA < stopA res = ($op)(nzvalA[ptrA], z) if res != z - rowA = rowvalA[ptrA] - rowvalS[ptrS] = rowA + rowA = indvalA[ptrA] + indvalS[ptrS] = rowA nzvalS[ptrS] = res ptrS += 1 end @@ -582,99 +667,101 @@ for (op, restype) in ( (:+, Nothing), (:-, Nothing), (:.*, Nothing), while ptrB < stopB res = ($op)(z, nzvalB[ptrB]) if res != z - rowB = rowvalB[ptrB] - rowvalS[ptrS] = rowB + rowB = indvalB[ptrB] + indvalS[ptrS] = rowB nzvalS[ptrS] = res ptrS += 1 end ptrB += 1 end - colptrS[col+1] = ptrS + indptrS[cidx+1] = ptrS end - deleteat!(rowvalS, colptrS[end]:length(rowvalS)) - deleteat!(nzvalS, colptrS[end]:length(nzvalS)) - return SparseMatrixCSC(m, n, colptrS, rowvalS, nzvalS) + deleteat!(indvalS, indptrS[end]:length(indvalS)) + deleteat!(nzvalS, indptrS[end]:length(nzvalS)) + return sparse(m, n, CompressedSparseStore(indptrS, indvalS, nzvalS), format=sptype(A)) end - end # quote end # macro -(.+)(A::SparseMatrixCSC, B::Number) = full(A) .+ B -( +)(A::SparseMatrixCSC, B::Array ) = full(A) + B -(.+)(A::Number, B::SparseMatrixCSC) = A .+ full(B) -( +)(A::Array , B::SparseMatrixCSC) = A + full(B) - -(.-)(A::SparseMatrixCSC, B::Number) = full(A) .- B -( -)(A::SparseMatrixCSC, B::Array ) = full(A) - B -(.-)(A::Number, B::SparseMatrixCSC) = A .- full(B) -( -)(A::Array , B::SparseMatrixCSC) = A - full(B) - -(.*)(A::SparseMatrixCSC, B::Number) = SparseMatrixCSC(A.m, A.n, copy(A.colptr), copy(A.rowval), A.nzval .* B) -(.*)(A::Number, B::SparseMatrixCSC) = SparseMatrixCSC(B.m, B.n, copy(B.colptr), copy(B.rowval), A .* B.nzval) -(.*)(A::SparseMatrixCSC, B::Array) = (.*)(A, sparse(B)) -(.*)(A::Array, B::SparseMatrixCSC) = (.*)(sparse(A), B) - -(./)(A::SparseMatrixCSC, B::Number) = SparseMatrixCSC(A.m, A.n, copy(A.colptr), copy(A.rowval), A.nzval ./ B) -(./)(A::Number, B::SparseMatrixCSC) = (./)(A, full(B)) -(./)(A::SparseMatrixCSC, B::Array) = (./)(full(A), B) -(./)(A::Array, B::SparseMatrixCSC) = (./)(A, full(B)) -(./)(A::SparseMatrixCSC, B::SparseMatrixCSC) = (./)(full(A), full(B)) - -(.\)(A::SparseMatrixCSC, B::Number) = (.\)(full(A), B) -(.\)(A::Number, B::SparseMatrixCSC) = SparseMatrixCSC(B.m, B.n, copy(B.colptr), copy(B.rowval), B.nzval .\ A) -(.\)(A::SparseMatrixCSC, B::Array) = (.\)(full(A), B) -(.\)(A::Array, B::SparseMatrixCSC) = (.\)(A, full(B)) -(.\)(A::SparseMatrixCSC, B::SparseMatrixCSC) = (.\)(full(A), full(B)) - -(.^)(A::SparseMatrixCSC, B::Number) = - B==0 ? sparse(ones(typeof(one(eltype(A)).^B), A.m, A.n)) : - SparseMatrixCSC(A.m, A.n, copy(A.colptr), copy(A.rowval), A.nzval .^ B) -(.^)(A::Number, B::SparseMatrixCSC) = (.^)(A, full(B)) -(.^)(A::SparseMatrixCSC, B::Array) = (.^)(full(A), B) -(.^)(A::Array, B::SparseMatrixCSC) = (.^)(A, full(B)) - -(.<)(A::SparseMatrixCSC, B::Number) = (.<)(full(A), B) -(.<)(A::Number, B::SparseMatrixCSC) = (.<)(A, full(B)) +(.+)(A::CompressedSparseMatrix, B::Number) = full(A) .+ B +( +)(A::CompressedSparseMatrix, B::Array ) = full(A) + B +(.+)(A::Number, B::CompressedSparseMatrix) = A .+ full(B) +( +)(A::Array , B::CompressedSparseMatrix) = A + full(B) + +(.-)(A::CompressedSparseMatrix, B::Number) = full(A) .- B +( -)(A::CompressedSparseMatrix, B::Array ) = full(A) - B +(.-)(A::Number, B::CompressedSparseMatrix) = A .- full(B) +( -)(A::Array, B::CompressedSparseMatrix) = A - full(B) + +(.*)(A::CompressedSparseMatrix, B::Number) = sparse(A.m, A.n, CompressedSparseStore(copy(getindptr(A)), copy(getindval(A)), A.nzval .* B), format=sptype(A)) +(.*)(A::Number, B::CompressedSparseMatrix) = sparse(B.m, B.n, CompressedSparseStore(copy(getindptr(B)), copy(getindval(B)), A .* B.nzval), format=sptype(B)) +(.*)(A::CompressedSparseMatrix, B::Array) = (.*)(A, sparse(B, format=sptype(A))) +(.*)(A::Array, B::CompressedSparseMatrix) = (.*)(sparse(A, format=sptype(B)), B) + +(./)(A::CompressedSparseMatrix, B::Number) = sparse(A.m, A.n, CompressedSparseStore(copy(getindptr(A)), copy(getindval(A)), A.nzval ./ B), format=sptype(A)) +(./)(A::Number, B::CompressedSparseMatrix) = (./)(A, full(B)) +(./)(A::CompressedSparseMatrix, B::Array) = (./)(full(A), B) +(./)(A::Array, B::CompressedSparseMatrix) = (./)(A, full(B)) +(./)(A::CompressedSparseMatrix, B::CompressedSparseMatrix) = (./)(full(A), full(B)) + +(.\)(A::CompressedSparseMatrix, B::Number) = (.\)(full(A), B) +(.\)(A::Number, B::CompressedSparseMatrix) = sparse(B.m, B.n, CompressedSparseStore(copy(getindptr(B)), copy(getindval(B)), B.nzval .\ A), format=sptype(B)) +(.\)(A::CompressedSparseMatrix, B::Array) = (.\)(full(A), B) +(.\)(A::Array, B::CompressedSparseMatrix) = (.\)(A, full(B)) +(.\)(A::CompressedSparseMatrix, B::CompressedSparseMatrix) = (.\)(full(A), full(B)) + +(.^)(A::CompressedSparseMatrix, B::Number) = + B==0 ? sparse(ones(typeof(one(eltype(A)).^B), A.m, A.n), format=sptype(A)) : + sparse(A.m, A.n, CompressedSparseStore(copy(getindptr(A)), copy(getindval(A)), A.nzval .^ B), format=sptype(A)) +(.^)(A::Number, B::CompressedSparseMatrix) = (.^)(A, full(B)) +(.^)(A::CompressedSparseMatrix, B::Array) = (.^)(full(A), B) +(.^)(A::Array, B::CompressedSparseMatrix) = (.^)(A, full(B)) + +(.<)(A::CompressedSparseMatrix, B::Number) = (.<)(full(A), B) +(.<)(A::Number, B::CompressedSparseMatrix) = (.<)(A, full(B)) # Reductions # TODO: Should the results of sparse reductions be sparse? -function reducedim{Tv,Ti}(f::Function, A::SparseMatrixCSC{Tv,Ti}, region, v0) - if region == 1 || region == (1,) +function reducedim{Tv,Ti}(f::Function, A::CompressedSparseMatrix{Tv,Ti}, region, v0) + cdim = spdim(A) + ncdim = nspdim(A) + + if (length(region) == 1) && (region[1] == spaxis(A)) + + S = fill(v0, ncdim, 1) + counts = zeros(Ti, ncdim) + for i = 1 : cdim, j = getindptr(A)[i] : getindptr(A)[i+1]-1 + ncidx = getindval(A)[j] + S[ncidx] = f(S[ncidx], A.nzval[j]) + counts[ncidx] += 1 + end + for i = 1:ncdim + if counts[i] != cdim; S[i] = f(S[i], zero(Tv)); end + end + return S + + elseif (length(region) == 1) && (region[1] == nspaxis(A)) - S = Array(Tv, 1, A.n) - for i = 1 : A.n + S = Array(Tv, 1, cdim) + for i = 1 : cdim Si = v0 - ccount = 0 - for j = A.colptr[i] : A.colptr[i+1]-1 + count = 0 + for j = getindptr(A)[i] : getindptr(A)[i+1]-1 Si = f(Si, A.nzval[j]) - ccount += 1 + count += 1 end - if ccount != A.m; Si = f(Si, zero(Tv)); end + if count != ncdim; Si = f(Si, zero(Tv)); end S[i] = Si end return S - elseif region == 2 || region == (2,) - - S = fill(v0, A.m, 1) - rcounts = zeros(Ti, A.m) - for i = 1 : A.n, j = A.colptr[i] : A.colptr[i+1]-1 - row = A.rowval[j] - S[row] = f(S[row], A.nzval[j]) - rcounts[row] += 1 - end - for i = 1:A.m - if rcounts[i] != A.n; S[i] = f(S[i], zero(Tv)); end - end - return S - elseif region == (1,2) S = v0 - for i = 1 : A.n, j = A.colptr[i] : A.colptr[i+1]-1 + for i = 1 : cdim, j = getindptr(A)[i] : getindptr(A)[i+1]-1 S = f(S, A.nzval[j]) end if nnz(A) != A.m*A.n; S = f(S, zero(Tv)); end @@ -686,29 +773,29 @@ function reducedim{Tv,Ti}(f::Function, A::SparseMatrixCSC{Tv,Ti}, region, v0) end end -function maximum{T}(A::SparseMatrixCSC{T}) +function maximum{T}(A::CompressedSparseMatrix{T}) isempty(A) && throw(ArgumentError("argument must not be empty")) m = maximum(A.nzval) nnz(A)!=length(A) ? max(m,zero(T)) : m end -maximum{T}(A::SparseMatrixCSC{T}, region) = +maximum{T}(A::CompressedSparseMatrix{T}, region) = isempty(A) ? similar(A, reduced_dims0(A,region)) : reducedim(Base.scalarmax,A,region,typemin(T)) -function minimum{T}(A::SparseMatrixCSC{T}) +function minimum{T}(A::CompressedSparseMatrix{T}) isempty(A) && throw(ArgumentError("argument must not be empty")) m = minimum(A.nzval) nnz(A)!=length(A) ? min(m,zero(T)) : m end -minimum{T}(A::SparseMatrixCSC{T}, region) = +minimum{T}(A::CompressedSparseMatrix{T}, region) = isempty(A) ? similar(A, reduced_dims0(A,region)) : reducedim(Base.scalarmin,A,region,typemax(T)) -sum{T}(A::SparseMatrixCSC{T}) = sum(A.nzval) -sum{T}(A::SparseMatrixCSC{T}, region) = reducedim(+,A,region,zero(T)) +sum{T}(A::CompressedSparseMatrix{T}) = sum(A.nzval) +sum{T}(A::CompressedSparseMatrix{T}, region) = reducedim(+,A,region,zero(T)) -prod{T}(A::SparseMatrixCSC{T}) = nnz(A)!=length(A) ? zero(T) : prod(A.nzval) -prod{T}(A::SparseMatrixCSC{T}, region) = reducedim(*,A,region,one(T)) +prod{T}(A::CompressedSparseMatrix{T}) = nnz(A)!=length(A) ? zero(T) : prod(A.nzval) +prod{T}(A::CompressedSparseMatrix{T}, region) = reducedim(*,A,region,one(T)) #all(A::SparseMatrixCSC{Bool}, region) = reducedim(all,A,region,true) #any(A::SparseMatrixCSC{Bool}, region) = reducedim(any,A,region,false) @@ -721,53 +808,59 @@ function rangesearch(haystack::Range, needle) (rem==0 && 1<=i+1<=length(haystack)) ? i+1 : 0 end -getindex(A::SparseMatrixCSC, i::Integer) = getindex(A, ind2sub(size(A),i)) -getindex(A::SparseMatrixCSC, I::(Integer,Integer)) = getindex(A, I[1], I[2]) +getindex{T<:CompressedSparseMatrix}(A::T, i::Integer) = getindex(A, ind2sub(size(A),i)) +getindex{T<:CompressedSparseMatrix}(A::T, I::(Integer,Integer)) = getindex(A, I[1], I[2]) -function getindex{T}(A::SparseMatrixCSC{T}, i0::Integer, i1::Integer) - if !(1 <= i0 <= A.m && 1 <= i1 <= A.n); throw(BoundsError()); end - r1 = int(A.colptr[i1]) - r2 = int(A.colptr[i1+1]-1) - (r1 > r2) && return zero(T) - r1 = searchsortedfirst(A.rowval, i0, r1, r2, Forward) - ((r1 > r2) || (A.rowval[r1] != i0)) ? zero(T) : A.nzval[r1] +function getindex_single{Tv,Ti}(iptr::Vector{Ti}, ival::Vector{Ti}, nzval::Vector{Tv}, ic::Int, inc::Int) + r1 = int(iptr[ic]) + r2 = int(iptr[ic+1]-1) + (r1 > r2) && return zero(Tv) + r1 = searchsortedfirst(ival, inc, r1, r2, Forward) + ((r1 > r2) || (ival[r1] != inc)) ? zero(Tv) : nzval[r1] end -getindex{T<:Integer}(A::SparseMatrixCSC, I::AbstractVector{T}, j::Integer) = getindex(A,I,[j]) -getindex{T<:Integer}(A::SparseMatrixCSC, i::Integer, J::AbstractVector{T}) = getindex(A,[i],J) +getindex{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, i0::Integer, i1::Integer) = + (1 <= i0 <= A.m && 1 <= i1 <= A.n) ? getindex_single(A.colptr, A.rowval, A.nzval, i1, i0) : throw(BoundsError()) +getindex{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti}, i0::Integer, i1::Integer) = + (1 <= i0 <= A.m && 1 <= i1 <= A.n) ? getindex_single(A.rowptr, A.colval, A.nzval, i0, i1) : throw(BoundsError()) -function getindex_cols{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, J::AbstractVector) - # for indexing whole columns - (m, n) = size(A) - nJ = length(J) +getindex{T<:Integer}(A::CompressedSparseMatrix, I::AbstractVector{T}, j::Integer) = getindex(A, I, [j]) +getindex{T<:Integer}(A::CompressedSparseMatrix, i::Integer, J::AbstractVector{T}) = getindex(A, [i], J) - colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval +function getindex_along_caxis{Tv,Ti}(iptr::Vector{Ti}, ival::Vector{Ti}, nzval::Vector{Tv}, X::AbstractVector) + nX = length(X) - colptrS = Array(Ti, nJ+1) - colptrS[1] = 1 + indptrS = Array(Ti, nX+1) + indptrS[1] = 1 nnzS = 0 - for j = 1:nJ - col = J[j] - nnzS += colptrA[col+1] - colptrA[col] - colptrS[j+1] = nnzS + 1 + for x = 1:nX + @inbounds cidx = X[x] + nnzS += iptr[cidx+1] - iptr[cidx] + @inbounds indptrS[x+1] = nnzS + 1 end - rowvalS = Array(Ti, nnzS) + indvalS = Array(Ti, nnzS) nzvalS = Array(Tv, nnzS) ptrS = 0 - for j = 1:nJ - col = J[j] - for k = colptrA[col]:colptrA[col+1]-1 + for x = 1:nX + @inbounds cidx = X[x] + + for ncidx = iptr[cidx]:iptr[cidx+1]-1 ptrS += 1 - rowvalS[ptrS] = rowvalA[k] - nzvalS[ptrS] = nzvalA[k] + indvalS[ptrS] = ival[ncidx] + nzvalS[ptrS] = nzval[ncidx] end end - return SparseMatrixCSC(m, nJ, colptrS, rowvalS, nzvalS) + (indptrS, indvalS, nzvalS) end +getindex_cols{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, J::AbstractVector) = + sparse(A.m, length(J), CompressedSparseStore(getindex_along_caxis(getindptr(A), getindval(A), A.nzval, J)...), format=CSC) +getindex_rows{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractVector) = + sparse(length(I), A.n, CompressedSparseStore(getindex_along_caxis(getindptr(A), getindval(A), A.nzval, I)...), format=CSR) + function getindex{Tv,Ti<:Integer}(A::SparseMatrixCSC{Tv,Ti}, I::Range, J::AbstractVector) # Ranges for indexing rows (m, n) = size(A) @@ -785,11 +878,11 @@ function getindex{Tv,Ti<:Integer}(A::SparseMatrixCSC{Tv,Ti}, I::Range, J::Abstra # Form the structure of the result and compute space for j = 1:nJ - col = J[j] - for k in colptrA[col]:colptrA[col+1]-1 + @inbounds col = J[j] + @simd for k in colptrA[col]:colptrA[col+1]-1 if rowvalA[k] in I; nnzS += 1 end # `in` is fast for ranges end - colptrS[j+1] = nnzS+1 + @inbounds colptrS[j+1] = nnzS+1 end # Populate the values in the result @@ -798,7 +891,7 @@ function getindex{Tv,Ti<:Integer}(A::SparseMatrixCSC{Tv,Ti}, I::Range, J::Abstra ptrS = 1 for j = 1:nJ - col = J[j] + @inbounds col = J[j] for k = colptrA[col]:colptrA[col+1]-1 rowA = rowvalA[k] i = rangesearch(I, rowA) @@ -813,26 +906,140 @@ function getindex{Tv,Ti<:Integer}(A::SparseMatrixCSC{Tv,Ti}, I::Range, J::Abstra return SparseMatrixCSC(nI, nJ, colptrS, rowvalS, nzvalS) end -function getindex_I_sorted{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::AbstractVector) - # Sorted vectors for indexing rows. - # Similar to getindex_general but without the transpose trick. +function getindex{Tv,Ti<:Integer}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractVector, J::Range) + # Ranges for indexing rows (m, n) = size(A) + # whole columns: + if J == 1:n + return getindex_rows(A, I) + end + nI = length(I) nJ = length(J) + rowptrA = A.rowptr; colvalA = A.colval; nzvalA = A.nzval + rowptrS = Array(Ti, nI+1) + rowptrS[1] = 1 + nnzS = 0 - colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval - colptrS = Array(Ti, nJ+1) - colptrS[1] = 1 + # Form the structure of the result and compute space + for i = 1:nI + @inbounds row = I[i] + @simd for k in rowptrA[row]:rowptrA[row+1]-1 + if colvalA[k] in J; nnzS += 1 end # `in` is fast for ranges + end + @inbounds rowptrS[i+1] = nnzS+1 + end + + # Populate the values in the result + colvalS = Array(Ti, nnzS) + nzvalS = Array(Tv, nnzS) + ptrS = 1 + + for i = 1:nI + @inbounds row = I[i] + for k = rowptrA[row]:rowptrA[row+1]-1 + colA = colvalA[k] + j = rangesearch(J, colA) + if j > 0 + colvalS[ptrS] = j + nzvalS[ptrS] = nzvalA[k] + ptrS += 1 + end + end + end + + return SparseMatrixCSR(nI, nJ, rowptrS, colvalS, nzvalS) +end + + +# TODO: See if growing arrays is faster than pre-computing structure +# and then populating nonzeros +# TODO: Use binary search in cases where nI >> nnz(A[:,j]) or nI << nnz(A[:,j]) +#= +function getindex_sorted{Tv,Ti}(m::Int, n::Int, iptr::Vector{Ti}, ival::Vector{Ti}, nzval::Vector{Tv}, I::AbstractVector, J::AbstractVector) + nI = length(I) + nJ = length(J) + + I_ref = falses(m) + I_ref[I] = true + + I_repeat = zeros(Int, m) + for i=1:nI; I_repeat[I[i]] += 1; end + + indptrS = Array(Ti, nJ+1) + indptrS[1] = 1 nnzS = 0 # Form the structure of the result and compute space for j = 1:nJ col = J[j] + + for k = iptr[col]:iptr[col+1]-1 + rowA = ival[k] + + if I_ref[rowA] + for r = 1:I_repeat[rowA] + nnzS += 1 + end + end + + end + indptrS[j+1] = nnzS+1 + end + + # Populate the values in the result + indvalS = Array(Ti, nnzS) + nzvalS = Array(Tv, nnzS) + ptrS = 1 + + fI = zeros(Ti, m) + for k=1:nI + Ik = I[k] + if fI[Ik] == 0; fI[Ik] = k; end + end + + for j = 1:nJ + col = J[j] + + for k = iptr[col]:iptr[col+1]-1 + rowA = ival[k] + + if I_ref[rowA] + for r = 1:I_repeat[rowA] + indvalS[ptrS] = fI[rowA] + r - 1 + nzvalS[ptrS] = nzval[k] + ptrS += 1 + end + end + + end + end + + return (indptrS, indvalS, nzvalS) +end +=# + +function getindex_sorted{Tv,Ti}(m::Int, n::Int, iptr::Vector{Ti}, ival::Vector{Ti}, nzval::Vector{Tv}, I::AbstractVector, J::AbstractVector) + # Sorted vectors for indexing rows. + # Similar to getindex_general but without the transpose trick. + nI = length(I) + nJ = length(J) + nV = length(ival) + (nV == length(nzval)) || throw(BoundsError()) + + indptrS = Array(Ti, nJ+1) + indptrS[1] = 1 + nnzS = 0 + + # Form the structure of the result and compute space + for j = 1:nJ + @inbounds col = J[j] ptrI::Int = 1 # runs through I - ptrA::Int = colptrA[col] - stopA::Int = colptrA[col+1] - while ptrI <= nI && ptrA < stopA - rowA = rowvalA[ptrA] + ptrA::Int = iptr[col] + stopA::Int = iptr[col+1] + (stopA <= (nV+1)) || throw(BoundsError()) + @inbounds while ptrI <= nI && ptrA < stopA + rowA = ival[ptrA] rowI = I[ptrI] if rowI > rowA @@ -844,22 +1051,22 @@ function getindex_I_sorted{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, ptrI += 1 end end - colptrS[j+1] = nnzS+1 + @inbounds indptrS[j+1] = nnzS+1 end # Populate the values in the result - rowvalS = Array(Ti, nnzS) + indvalS = Array(Ti, nnzS) nzvalS = Array(Tv, nnzS) ptrS = 1 - for j = 1:nJ + @inbounds for j = 1:nJ col = J[j] ptrI::Int = 1 # runs through I - ptrA::Int = colptrA[col] - stopA::Int = colptrA[col+1] + ptrA::Int = iptr[col] + stopA::Int = iptr[col+1] while ptrI <= nI && ptrA < stopA - rowA = rowvalA[ptrA] + rowA = ival[ptrA] rowI = I[ptrI] if rowI > rowA @@ -867,100 +1074,110 @@ function getindex_I_sorted{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, elseif rowI < rowA ptrI += 1 else - rowvalS[ptrS] = ptrI - nzvalS[ptrS] = nzvalA[ptrA] + indvalS[ptrS] = ptrI + nzvalS[ptrS] = nzval[ptrA] ptrS += 1 ptrI += 1 end end end - return SparseMatrixCSC(nI, nJ, colptrS, rowvalS, nzvalS) + return (indptrS, indvalS, nzvalS) end -function getindex_general{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::AbstractVector) +function getindex_I_sorted{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::AbstractVector) + (indptrS, indvalS, nzvalS) = getindex_sorted(A.m, A.n, getindptr(A), getindval(A), A.nzval, I, J) + SparseMatrixCSC(length(I), length(J), indptrS, indvalS, nzvalS) +end +function getindex_J_sorted{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractVector, J::AbstractVector) + (indptrS, indvalS, nzvalS) = getindex_sorted(A.n, A.m, getindptr(A), getindval(A), A.nzval, J, I) + SparseMatrixCSR(length(I), length(J), indptrS, indvalS, nzvalS) +end + +function getindex_general{Tv,Ti}(iptr::Vector{Ti}, ival::Vector{Ti}, nzval::Vector{Tv}, NC::Vector, C::AbstractVector) # Anything for indexing rows. # This sorts I first then does some trick with constructing the transpose. - (m, n) = size(A) - nI = length(I) - nJ = length(J) + nNC = length(NC) + nC = length(C) + nV = length(ival) + (nV == length(nzval)) || throw(BoundsError()) - colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval nnzS = 0 - pI = sortperm(I); I = I[pI] - fI = find(I) - W = zeros(Ti, nI + 1) # Keep row counts - W[1] = 1 # For cumsum later + pNC = sortperm(NC); NC = NC[pNC] + fNC = find(NC) + W = zeros(Ti, nNC + 1) # Keep counts for the non compressed axis + W[1] = 1 # For cumsum later # Form the structure of the result and compute space - for j = 1:nJ - col = J[j] - - ptrI::Int = 1 - - ptrA::Int = colptrA[col] - stopA::Int = colptrA[col+1] - - while ptrI <= nI && ptrA < stopA - rowA = rowvalA[ptrA] - rowI = I[ptrI] - - if rowI > rowA + for idx = 1:nC + @inbounds cidx = C[idx] + ptrNC::Int = 1 + ptrA::Int = iptr[cidx] + stopA::Int = iptr[cidx+1] + (stopA <= (nV+1)) || throw(BoundsError()) + @inbounds while ptrNC <= nNC && ptrA < stopA + idxA = ival[ptrA] + idxNC = NC[ptrNC] + if idxNC > idxA ptrA += 1 - elseif rowI < rowA - ptrI += 1 + elseif idxNC < idxA + ptrNC += 1 else - W[fI[pI[ptrI]]+1] += 1 + W[fNC[pNC[ptrNC]]+1] += 1 nnzS += 1 - ptrI += 1 + ptrNC += 1 end end - end - colptrS_T = cumsum(W) + indptrS_T = cumsum(W) # Populate the values in the result, but transposed - rowvalS_T = Array(Ti, nnzS) + indvalS_T = Array(Ti, nnzS) nzvalS_T = Array(Tv, nnzS) - for i=1:nI; W[i] = 0; end # Zero out W to store row positions - - for j = 1:nJ - col = J[j] - - ptrI::Int = 1 - - ptrA::Int = colptrA[col] - stopA::Int = colptrA[col+1] - - while ptrI <= nI && ptrA < stopA - rowA = rowvalA[ptrA] - rowI = I[ptrI] - - if rowI > rowA + @simd for i=1:nNC; @inbounds W[i] = 0; end # Zero out W to store positions + + @inbounds for idx = 1:nC + cidx = C[idx] + ptrNC::Int = 1 + ptrA::Int = iptr[cidx] + stopA::Int = iptr[cidx+1] + while ptrNC <= nNC && ptrA < stopA + idxA = ival[ptrA] + idxNC = NC[ptrNC] + if idxNC > idxA ptrA += 1 - elseif rowI < rowA - ptrI += 1 + elseif idxNC < idxA + ptrNC += 1 else - rowS = fI[pI[ptrI]] - k = colptrS_T[rowS] + W[rowS] - rowvalS_T[k] = j - nzvalS_T[k] = nzvalA[ptrA] - W[rowS] += 1 - ptrI += 1 + idxS = fNC[pNC[ptrNC]] + k = indptrS_T[idxS] + W[idxS] + indvalS_T[k] = idx + nzvalS_T[k] = nzval[ptrA] + W[idxS] += 1 + ptrNC += 1 end end - end + (indptrS_T, indvalS_T, nzvalS_T) +end + +function getindex_general{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, I::Vector, J::AbstractVector) + (indptrS_T, indvalS_T, nzvalS_T) = getindex_general(getindptr(A), getindval(A), A.nzval, I, J) # Transpose so that rows are in sorted order and return - S_T = SparseMatrixCSC(nJ, nI, colptrS_T, rowvalS_T, nzvalS_T) + S_T = SparseMatrixCSC(length(J), length(I), indptrS_T, indvalS_T, nzvalS_T) return S_T.' end +function getindex_general{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractVector, J::Vector) + (indptrS_T, indvalS_T, nzvalS_T) = getindex_general(getindptr(A), getindval(A), A.nzval, J, I) + # Transpose so that cols are in sorted order and return + S_T = SparseMatrixCSR(length(J), length(I), indptrS_T, indvalS_T, nzvalS_T) + return S_T.' +end -# the general case: +# S = A[I, J] function getindex{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::AbstractVector) - loop_overI_threshold = 0 # ~40 but doesn't matter much. if issorted(I) return getindex_I_sorted(A, I, J) else @@ -968,214 +1185,206 @@ function getindex{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::Abstra end end +function getindex{Tv,Ti}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractVector, J::AbstractVector) + if issorted(J) + return getindex_J_sorted(A, I, J) + else + return getindex_general(A, I, J) + end +end + # logical getindex getindex{Tv,Ti<:Integer}(A::SparseMatrixCSC{Tv,Ti}, I::Range{Bool}, J::AbstractVector{Bool}) = error("Cannot index with Range{Bool}") getindex{Tv,Ti<:Integer,T<:Integer}(A::SparseMatrixCSC{Tv,Ti}, I::Range{Bool}, J::AbstractVector{T}) = error("Cannot index with Range{Bool}") - +getindex{Tv,Ti<:Integer}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractVector{Bool}, J::Range{Bool}) = error("Cannot index with Range{Bool}") +getindex{Tv,Ti<:Integer,T<:Integer}(A::SparseMatrixCSR{Tv,Ti}, I::AbstractVector{T}, J::Range{Bool}) = error("Cannot index with Range{Bool}") + getindex{T<:Integer}(A::SparseMatrixCSC, I::Range{T}, J::AbstractVector{Bool}) = A[I,find(J)] -getindex(A::SparseMatrixCSC, I::Integer, J::AbstractVector{Bool}) = A[I,find(J)] -getindex(A::SparseMatrixCSC, I::AbstractVector{Bool}, J::Integer) = A[find(I),J] +getindex{T<:Integer}(A::SparseMatrixCSR, I::AbstractVector{Bool}, J::Range{T}) = A[find(I),J] + +getindex(A::CompressedSparseMatrix, I::Integer, J::AbstractVector{Bool}) = A[I,find(J)] +getindex(A::CompressedSparseMatrix, I::AbstractVector{Bool}, J::Integer) = A[find(I),J] getindex(A::SparseMatrixCSC, I::AbstractVector{Bool}, J::AbstractVector{Bool}) = A[find(I),find(J)] +getindex(A::SparseMatrixCSR, I::AbstractVector{Bool}, J::AbstractVector{Bool}) = A[find(I),find(J)] getindex{T<:Integer}(A::SparseMatrixCSC, I::AbstractVector{T}, J::AbstractVector{Bool}) = A[I,find(J)] +getindex{T<:Integer}(A::SparseMatrixCSR, I::AbstractVector{T}, J::AbstractVector{Bool}) = A[I,find(J)] getindex{T<:Integer}(A::SparseMatrixCSC, I::AbstractVector{Bool}, J::AbstractVector{T}) = A[find(I),J] +getindex{T<:Integer}(A::SparseMatrixCSR, I::AbstractVector{Bool}, J::AbstractVector{T}) = A[find(I),J] -function getindex{Tv}(A::SparseMatrixCSC{Tv}, I::AbstractArray{Bool}) - checkbounds(A, I) - n = sum(I) - - colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval - colptrB = Int[1,n+1] - rowvalB = Array(Int, n) - nzvalB = Array(Tv, n) - c = 1 - rowB = 1 - - for col in 1:A.n - r1 = colptrA[col] - r2 = colptrA[col+1]-1 - - for row in 1:A.m - if I[row, col] - while (r1 <= r2) && (rowvalA[r1] < row) - r1 += 1 - end - if (r1 <= r2) && (rowvalA[r1] == row) - nzvalB[c] = nzvalA[r1] - rowvalB[c] = rowB - c += 1 - end - rowB += 1 - (rowB > n) && break - end - end - (rowB > n) && break - end - colptrB[end] = c - n = length(nzvalB) - if n > (c-1) - deleteat!(nzvalB, c:n) - deleteat!(rowvalB, c:n) - end - SparseMatrixCSC(n, 1, colptrB, rowvalB, nzvalB) -end - -function getindex{Tv}(A::SparseMatrixCSC{Tv}, I::AbstractArray) - szA = size(A); colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval +function getindex{Tv,Ti}(A::CompressedSparseMatrix{Tv,Ti}, I::AbstractArray) + szA = size(A); indptrA = getindptr(A); indvalA = getindval(A); nzvalA = A.nzval n = length(I) outm = size(I,1) outn = size(I,2) szB = (outm, outn) - colptrB = zeros(Int, outn+1) - rowvalB = Array(Int, n) + Ts = sptype(A) + indptrB = zeros(Int, ((Ts===CSC)?outn:outm)+1) + indvalB = Array(Int, n) nzvalB = Array(Tv, n) - colB = 1 - rowB = 1 - colptrB[colB] = 1 + spidxB = nspidxB = 1 + indptrB[spidxB] = 1 idxB = 1 + i_nsp = nspaxis(A); i_sp = spaxis(A) for i in 1:n - row,col = ind2sub(szA, I[i]) - for r in colptrA[col]:(colptrA[col+1]-1) - if rowvalA[r] == row - rowB,colB = ind2sub(szB, i) - colptrB[colB+1] += 1 - rowvalB[idxB] = rowB + subA = ind2sub(szA, I[i]) + @inbounds nspidxA = subA[i_nsp] + @inbounds spidxA = subA[i_sp] + subB = ind2sub(szB, i) + @inbounds nspidxB = subB[i_nsp] + @inbounds spidxB = subB[i_sp] + + for r in indptrA[spidxA]:(indptrA[spidxA+1]-1) + if indvalA[r] == nspidxA + indptrB[spidxB+1] += 1 + indvalB[idxB] = nspidxB nzvalB[idxB] = nzvalA[r] idxB += 1 break end end end - colptrB = cumsum(colptrB) + indptrB = cumsum(indptrB) if n > (idxB-1) deleteat!(nzvalB, idxB:n) - deleteat!(rowvalB, idxB:n) + deleteat!(indvalB, idxB:n) end - SparseMatrixCSC(outm, outn, colptrB, rowvalB, nzvalB) + sparse(outm, outn, CompressedSparseStore(indptrB, indvalB, nzvalB), format=Ts) end +getindex{Tv,Ti}(A::CompressedSparseMatrix{Tv,Ti}, I::AbstractArray{Bool}) = getindex(A, find(I)) ## setindex! -setindex!(A::SparseMatrixCSC, v, i::Integer) = setindex!(A, v, ind2sub(size(A),i)...) +setindex!(A::CompressedSparseMatrix, v, i::Integer) = setindex!(A, v, ind2sub(size(A),i)...) -function setindex!{T,Ti}(A::SparseMatrixCSC{T,Ti}, v, i0::Integer, i1::Integer) +function setindex_single!{Tv,Ti}(iptr::Vector{Ti}, ival::Vector{Ti}, nzval::Vector{Tv}, v, i0::Integer, i1::Integer) i0 = convert(Ti, i0) i1 = convert(Ti, i1) - if !(1 <= i0 <= A.m && 1 <= i1 <= A.n); throw(BoundsError()); end - v = convert(T, v) - r1 = int(A.colptr[i1]) - r2 = int(A.colptr[i1+1]-1) + v = convert(Tv, v) + + r1 = int(iptr[i1]) + r2 = int(iptr[i1+1]-1) + N = length(iptr) if v == 0 #either do nothing or delete entry if it exists if r1 <= r2 - r1 = searchsortedfirst(A.rowval, i0, r1, r2, Forward) - if (r1 <= r2) && (A.rowval[r1] == i0) - deleteat!(A.rowval, r1) - deleteat!(A.nzval, r1) - for j = (i1+1):(A.n+1) - A.colptr[j] -= 1 + r1 = searchsortedfirst(ival, i0, r1, r2, Forward) + if (r1 <= r2) && (ival[r1] == i0) + deleteat!(ival, r1) + deleteat!(ival, r1) + for j = (i1+1):N + iptr[j] -= 1 end end end - return A - end - i = (r1 > r2) ? r1 : searchsortedfirst(A.rowval, i0, r1, r2, Forward) - - if (i <= r2) && (A.rowval[i] == i0) - A.nzval[i] = v else - insert!(A.rowval, i, i0) - insert!(A.nzval, i, v) - for j = (i1+1):(A.n+1) - A.colptr[j] += 1 + i = (r1 > r2) ? r1 : searchsortedfirst(ival, i0, r1, r2, Forward) + + if (i <= r2) && (ival[i] == i0) + nzval[i] = v + else + insert!(ival, i, i0) + insert!(nzval, i, v) + for j = (i1+1):N + iptr[j] += 1 + end end end - return A + nothing end -setindex!{T<:Integer}(A::SparseMatrixCSC, v::AbstractMatrix, i::Integer, J::AbstractVector{T}) = setindex!(A, v, [i], J) -setindex!{T<:Integer}(A::SparseMatrixCSC, v::AbstractMatrix, I::AbstractVector{T}, j::Integer) = setindex!(A, v, I, [j]) +function setindex!{Tv,Ti}(A::CompressedSparseMatrix{Tv,Ti}, v, i0::Integer, i1::Integer) + (1 <= i0 <= A.m && 1 <= i1 <= A.n) || throw(BoundsError()) + Ts = sptype(A) + setindex_single!(getindptr(A), getindval(A), A.nzval, v, (Ts===CSC)?i0:i1, (Ts===CSC)?i1:i0) + A +end + +setindex!{T<:Integer}(A::CompressedSparseMatrix, v::AbstractMatrix, i::Integer, J::AbstractVector{T}) = setindex!(A, v, [i], J) +setindex!{T<:Integer}(A::CompressedSparseMatrix, v::AbstractMatrix, I::AbstractVector{T}, j::Integer) = setindex!(A, v, I, [j]) -setindex!{T<:Integer}(A::SparseMatrixCSC, x::Number, i::Integer, J::AbstractVector{T}) = setindex!(A, x, [i], J) -setindex!{T<:Integer}(A::SparseMatrixCSC, x::Number, I::AbstractVector{T}, j::Integer) = setindex!(A, x, I, [j]) +setindex!{T<:Integer}(A::CompressedSparseMatrix, x::Number, i::Integer, J::AbstractVector{T}) = setindex!(A, x, [i], J) +setindex!{T<:Integer}(A::CompressedSparseMatrix, x::Number, I::AbstractVector{T}, j::Integer) = setindex!(A, x, I, [j]) -setindex!{Tv,T<:Integer}(A::SparseMatrixCSC{Tv}, x::Number, I::AbstractVector{T}, J::AbstractVector{T}) = +setindex!{Tv,T<:Integer}(A::CompressedSparseMatrix{Tv}, x::Number, I::AbstractVector{T}, J::AbstractVector{T}) = (0 == x) ? spdelete!(A, I, J) : spset!(A, convert(Tv,x), I, J) -function spset!{Tv,Ti<:Integer}(A::SparseMatrixCSC{Tv}, x::Tv, I::AbstractVector{Ti}, J::AbstractVector{Ti}) +function spset!{Tv,Ti<:Integer}(A::CompressedSparseMatrix{Tv}, x::Tv, I::AbstractVector{Ti}, J::AbstractVector{Ti}) !issorted(I) && (I = I[sortperm(I)]) !issorted(J) && (J = J[sortperm(J)]) - m, n = size(A) - ((I[end] > m) || (J[end] > n)) && throw(DimensionMismatch("")) + ((I[end] > size(A,1)) || (J[end] > size(A,2))) && throw(DimensionMismatch("")) nnzA = nnz(A) + length(I) * length(J) - colptrA = colptr = A.colptr - rowvalA = rowval = A.rowval + indptrA = iptr = getindptr(A) + indvalA = ival = getindval(A) nzvalA = nzval = A.nzval - rowidx = 1 + SPX = (sptype(A) === CSC) ? J : I + NSPX = (sptype(A) === CSC) ? I : J + + ncidx = 1 nadd = 0 - for col in 1:n - rrange = colptr[col]:(colptr[col+1]-1) - (nadd > 0) && (colptrA[col] = colptr[col] + nadd) + for cidx in 1:spdim(A) + rrange = iptr[cidx]:(iptr[cidx+1]-1) + (nadd > 0) && (indptrA[cidx] = iptr[cidx] + nadd) - if col in J + if cidx in SPX if isempty(rrange) # set new vals only - nincl = length(I) + nincl = length(NSPX) if nadd == 0 - colptrA = copy(colptr) - rowvalA = Array(Ti, nnzA); copy!(rowvalA, 1, rowval, 1, length(rowval)) + indptrA = copy(iptr) + indvalA = Array(Ti, nnzA); copy!(indvalA, 1, ival, 1, length(ival)) nzvalA = Array(Tv, nnzA); copy!(nzvalA, 1, nzval, 1, length(nzval)) end - r = rowidx:(rowidx+nincl-1) - rowvalA[r] = I + r = ncidx:(ncidx+nincl-1) + indvalA[r] = NSPX nzvalA[r] = x - rowidx += nincl + ncidx += nincl nadd += nincl else # set old + new vals old_ptr = rrange[1] old_stop = rrange[end] new_ptr = 1 - new_stop = length(I) + new_stop = length(NSPX) while true - old_row = rowval[old_ptr] - new_row = I[new_ptr] + old_row = ival[old_ptr] + new_row = NSPX[new_ptr] if old_row < new_row - rowvalA[rowidx] = old_row - nzvalA[rowidx] = nzval[old_ptr] - rowidx += 1 + indvalA[ncidx] = old_row + nzvalA[ncidx] = nzval[old_ptr] + ncidx += 1 old_ptr += 1 else if old_row == new_row old_ptr += 1 else if nadd == 0 - colptrA = copy(colptr) - rowvalA = Array(Ti, nnzA); copy!(rowvalA, 1, rowval, 1, length(rowval)) + indptrA = copy(iptr) + indvalA = Array(Ti, nnzA); copy!(indvalA, 1, ival, 1, length(ival)) nzvalA = Array(Tv, nnzA); copy!(nzvalA, 1, nzval, 1, length(nzval)) end nadd += 1 end - rowvalA[rowidx] = new_row - nzvalA[rowidx] = x - rowidx += 1 + indvalA[ncidx] = new_row + nzvalA[ncidx] = x + ncidx += 1 new_ptr += 1 end if old_ptr > old_stop if new_ptr <= new_stop if nadd == 0 - colptrA = copy(colptr) - rowvalA = Array(Ti, nnzA); copy!(rowvalA, 1, rowval, 1, length(rowval)) + indptrA = copy(iptr) + indvalA = Array(Ti, nnzA); copy!(indvalA, 1, ival, 1, length(ival)) nzvalA = Array(Tv, nnzA); copy!(nzvalA, 1, nzval, 1, length(nzval)) end - r = rowidx:(rowidx+(new_stop-new_ptr)) - rowvalA[r] = I[new_ptr:new_stop] + r = ncidx:(ncidx+(new_stop-new_ptr)) + indvalA[r] = NSPX[new_ptr:new_stop] nzvalA[r] = x - rowidx += length(r) + ncidx += length(r) nadd += length(r) end break @@ -1183,107 +1392,109 @@ function spset!{Tv,Ti<:Integer}(A::SparseMatrixCSC{Tv}, x::Tv, I::AbstractVector if new_ptr > new_stop nincl = old_stop-old_ptr+1 - copy!(rowvalA, rowidx, rowval, old_ptr, nincl) - copy!(nzvalA, rowidx, nzval, old_ptr, nincl) - rowidx += nincl + copy!(indvalA, ncidx, ival, old_ptr, nincl) + copy!(nzvalA, ncidx, nzval, old_ptr, nincl) + ncidx += nincl break end end end elseif !isempty(rrange) # set old vals only nincl = length(rrange) - copy!(rowvalA, rowidx, rowval, rrange[1], nincl) - copy!(nzvalA, rowidx, nzval, rrange[1], nincl) - rowidx += nincl + copy!(indvalA, ncidx, ival, rrange[1], nincl) + copy!(nzvalA, ncidx, nzval, rrange[1], nincl) + ncidx += nincl end end if nadd > 0 - colptrA[n+1] = rowidx - deleteat!(rowvalA, rowidx:nnzA) - deleteat!(nzvalA, rowidx:nnzA) + indptrA[end] = ncidx + deleteat!(indvalA, ncidx:nnzA) + deleteat!(nzvalA, ncidx:nnzA) - A.colptr = colptrA - A.rowval = rowvalA + setindptr!(A, indptrA) + setindval!(A, indvalA) A.nzval = nzvalA end return A end -function spdelete!{Tv,Ti<:Integer}(A::SparseMatrixCSC{Tv}, I::AbstractVector{Ti}, J::AbstractVector{Ti}) - m, n = size(A) +function spdelete!{Tv,Ti<:Integer}(A::CompressedSparseMatrix{Tv}, I::AbstractVector{Ti}, J::AbstractVector{Ti}) nnzA = nnz(A) (nnzA == 0) && (return A) !issorted(I) && (I = I[sortperm(I)]) !issorted(J) && (J = J[sortperm(J)]) - ((I[end] > m) || (J[end] > n)) && throw(DimensionMismatch("")) + ((I[end] > size(A,1)) || (J[end] > size(A,2))) && throw(DimensionMismatch("")) - colptr = colptrA = A.colptr - rowval = rowvalA = A.rowval + iptr = indptrA = getindptr(A) + ival = indvalA = getindval(A) nzval = nzvalA = A.nzval - rowidx = 1 + + SPX = (sptype(A) === CSC) ? J : I + NSPX = (sptype(A) === CSC) ? I : J + + ncidx = 1 ndel = 0 - for col in 1:n - rrange = colptr[col]:(colptr[col+1]-1) - (ndel > 0) && (colptrA[col] = colptr[col] - ndel) - if isempty(rrange) || !(col in J) + for cidx in 1:spdim(A) + rrange = iptr[cidx]:(iptr[cidx+1]-1) + (ndel > 0) && (indptrA[cidx] = iptr[cidx] - ndel) + if isempty(rrange) || !(cidx in SPX) nincl = length(rrange) if(ndel > 0) && !isempty(rrange) - copy!(rowvalA, rowidx, rowval, rrange[1], nincl) - copy!(nzvalA, rowidx, nzval, rrange[1], nincl) + copy!(indvalA, ncidx, ival, rrange[1], nincl) + copy!(nzvalA, ncidx, nzval, rrange[1], nincl) end - rowidx += nincl + ncidx += nincl else for ridx in rrange - if rowval[ridx] in I + if ival[ridx] in NSPX if ndel == 0 - colptrA = copy(colptr) - rowvalA = copy(rowval) + indptrA = copy(iptr) + indvalA = copy(ival) nzvalA = copy(nzval) end ndel += 1 else if ndel > 0 - rowvalA[rowidx] = rowval[ridx] - nzvalA[rowidx] = nzval[ridx] + indvalA[ncidx] = ival[ridx] + nzvalA[ncidx] = nzval[ridx] end - rowidx += 1 + ncidx += 1 end end end end if ndel > 0 - colptrA[n+1] = rowidx - deleteat!(rowvalA, rowidx:nnzA) - deleteat!(nzvalA, rowidx:nnzA) + indptrA[end] = ncidx + deleteat!(indvalA, ncidx:nnzA) + deleteat!(nzvalA, ncidx:nnzA) - A.colptr = colptrA - A.rowval = rowvalA + setindptr!(A, indptrA) + setindval!(A, indvalA) A.nzval = nzvalA end return A end -setindex!{Tv,Ti,T<:Integer}(A::SparseMatrixCSC{Tv,Ti}, S::Matrix, I::AbstractVector{T}, J::AbstractVector{T}) = - setindex!(A, convert(SparseMatrixCSC{Tv,Ti}, S), I, J) +setindex!{Tv,Ti,T<:Integer}(A::CompressedSparseMatrix{Tv,Ti}, S::Matrix, I::AbstractVector{T}, J::AbstractVector{T}) = + setindex!(A, convert(typeof(A), S), I, J) -setindex!{Tv,Ti,T<:Integer}(A::SparseMatrixCSC{Tv,Ti}, v::AbstractVector, I::AbstractVector{T}, j::Integer) = setindex!(A, v, I, [j]) -setindex!{Tv,Ti,T<:Integer}(A::SparseMatrixCSC{Tv,Ti}, v::AbstractVector, i::Integer, J::AbstractVector{T}) = setindex!(A, v, [i], J) -setindex!{Tv,Ti,T<:Integer}(A::SparseMatrixCSC{Tv,Ti}, v::AbstractVector, I::AbstractVector{T}, J::AbstractVector{T}) = +setindex!{Tv,Ti,T<:Integer}(A::CompressedSparseMatrix{Tv,Ti}, v::AbstractVector, I::AbstractVector{T}, j::Integer) = setindex!(A, v, I, [j]) +setindex!{Tv,Ti,T<:Integer}(A::CompressedSparseMatrix{Tv,Ti}, v::AbstractVector, i::Integer, J::AbstractVector{T}) = setindex!(A, v, [i], J) +setindex!{Tv,Ti,T<:Integer}(A::CompressedSparseMatrix{Tv,Ti}, v::AbstractVector, I::AbstractVector{T}, J::AbstractVector{T}) = setindex!(A, reshape(v, length(I), length(J)), I, J) # A[I,J] = B -function setindex!{Tv,Ti,T<:Integer}(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}, I::AbstractVector{T}, J::AbstractVector{T}) - if size(B,1) != length(I) || size(B,2) != length(J) - throw(DimensionMismatch("")) - end +function setindex!{Tv,Ti,T<:Integer}(A::CompressedSparseMatrix{Tv,Ti}, B::CompressedSparseMatrix{Tv,Ti}, I::AbstractVector{T}, J::AbstractVector{T}) + (sptype(A) == sptype(B)) || (B = convert(typeof(A), B)) + Ts = sptype(A) + (size(B,1) != length(I) || size(B,2) != length(J)) && throw(DimensionMismatch("")) - issortedI = issorted(I) - issortedJ = issorted(J) + issortedI = issorted(I); issortedJ = issorted(J) if !issortedI && !issortedJ pI = sortperm(I); I = I[pI] @@ -1297,65 +1508,70 @@ function setindex!{Tv,Ti,T<:Integer}(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixC B = B[:, pJ] end - m, n = size(A) - mB, nB = size(B) + (indptrS, indvalS, nzvalS) = setindex_general!(getindptr(A), getindval(A), A.nzval, nspdim(A), + getindptr(B), getindval(B), B.nzval, + (Ts===CSC)?I:J, (Ts===CSC)?J:I) + setindptr!(A, indptrS) + setindval!(A, indvalS) + A.nzval = nzvalS + return A +end - nI = length(I) +function setindex_general!{Tv,Ti,T<:Integer}(indptrA::Vector{Ti}, indvalA::Vector{Ti}, nzvalA::Vector{Tv}, nspsz::Int, + indptrB::Vector{Ti}, indvalB::Vector{Ti}, nzvalB::Vector{Tv}, + I::AbstractVector{T}, J::AbstractVector{T}) nJ = length(J) - colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval - colptrB = B.colptr; rowvalB = B.rowval; nzvalB = B.nzval - - nnzS = nnz(A) + nnz(B) - colptrS = Array(Ti, n+1) - rowvalS = Array(Ti, nnzS) + nnzS = indptrA[end] + indptrB[end] - 2 + indptrS = Array(Ti, length(indptrA)) + indvalS = Array(Ti, nnzS) nzvalS = Array(Tv, nnzS) - colptrS[1] = 1 + indptrS[1] = 1 colB = 1 asgn_col = J[colB] - I_asgn = falses(m) + I_asgn = falses(nspsz) I_asgn[I] = true ptrS = 1 - for col = 1:n + for col = 1:(length(indptrA)-1) # Copy column of A if it is not being assigned into if colB > nJ || col != J[colB] - colptrS[col+1] = colptrS[col] + (colptrA[col+1]-colptrA[col]) + indptrS[col+1] = indptrS[col] + (indptrA[col+1]-indptrA[col]) - for k = colptrA[col]:colptrA[col+1]-1 - rowvalS[ptrS] = rowvalA[k] + for k = indptrA[col]:indptrA[col+1]-1 + indvalS[ptrS] = indvalA[k] nzvalS[ptrS] = nzvalA[k] ptrS += 1 end continue end - ptrA::Int = colptrA[col] - stopA::Int = colptrA[col+1] - ptrB::Int = colptrB[colB] - stopB::Int = colptrB[colB+1] + ptrA::Int = indptrA[col] + stopA::Int = indptrA[col+1] + ptrB::Int = indptrB[colB] + stopB::Int = indptrB[colB+1] while ptrA < stopA && ptrB < stopB - rowA = rowvalA[ptrA] - rowB = I[rowvalB[ptrB]] + rowA = indvalA[ptrA] + rowB = I[indvalB[ptrB]] if rowA < rowB if ~I_asgn[rowA] - rowvalS[ptrS] = rowA + indvalS[ptrS] = rowA nzvalS[ptrS] = nzvalA[ptrA] ptrS += 1 end ptrA += 1 elseif rowB < rowA - rowvalS[ptrS] = rowB + indvalS[ptrS] = rowB nzvalS[ptrS] = nzvalB[ptrB] ptrS += 1 ptrB += 1 else - rowvalS[ptrS] = rowB + indvalS[ptrS] = rowB nzvalS[ptrS] = nzvalB[ptrB] ptrS += 1 ptrB += 1 @@ -1364,9 +1580,9 @@ function setindex!{Tv,Ti,T<:Integer}(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixC end while ptrA < stopA - rowA = rowvalA[ptrA] + rowA = indvalA[ptrA] if ~I_asgn[rowA] - rowvalS[ptrS] = rowA + indvalS[ptrS] = rowA nzvalS[ptrS] = nzvalA[ptrA] ptrS += 1 end @@ -1374,139 +1590,166 @@ function setindex!{Tv,Ti,T<:Integer}(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixC end while ptrB < stopB - rowB = I[rowvalB[ptrB]] - rowvalS[ptrS] = rowB + rowB = I[indvalB[ptrB]] + indvalS[ptrS] = rowB nzvalS[ptrS] = nzvalB[ptrB] ptrS += 1 ptrB += 1 end - colptrS[col+1] = ptrS + indptrS[col+1] = ptrS colB += 1 end - deleteat!(rowvalS, colptrS[end]:length(rowvalS)) - deleteat!(nzvalS, colptrS[end]:length(nzvalS)) - - A.colptr = colptrS - A.rowval = rowvalS - A.nzval = nzvalS - return A + deleteat!(indvalS, indptrS[end]:length(indvalS)) + deleteat!(nzvalS, indptrS[end]:length(nzvalS)) + (indptrS, indvalS, nzvalS) end + # Logical setindex! -setindex!(A::SparseMatrixCSC, x::Matrix, I::Integer, J::AbstractVector{Bool}) = setindex!(A, sparse(x), I, find(J)) -setindex!(A::SparseMatrixCSC, x::Matrix, I::AbstractVector{Bool}, J::Integer) = setindex!(A, sparse(x), find(I), J) -setindex!(A::SparseMatrixCSC, x::Matrix, I::AbstractVector{Bool}, J::AbstractVector{Bool}) = setindex!(A, sparse(x), find(I), find(J)) -setindex!{T<:Integer}(A::SparseMatrixCSC, x::Matrix, I::AbstractVector{T}, J::AbstractVector{Bool}) = setindex!(A, sparse(x), I, find(J)) -setindex!{T<:Integer}(A::SparseMatrixCSC, x::Matrix, I::AbstractVector{Bool}, J::AbstractVector{T}) = setindex!(A, sparse(x), find(I),J) +setindex!(A::CompressedSparseMatrix, x::Matrix, I::Integer, J::AbstractVector{Bool}) = setindex!(A, sparse(x, format=sptype(A)), I, find(J)) +setindex!(A::CompressedSparseMatrix, x::Matrix, I::AbstractVector{Bool}, J::Integer) = setindex!(A, sparse(x, format=sptype(A)), find(I), J) +setindex!(A::CompressedSparseMatrix, x::Matrix, I::AbstractVector{Bool}, J::AbstractVector{Bool}) = setindex!(A, sparse(x, format=sptype(A)), find(I), find(J)) +setindex!{T<:Integer}(A::CompressedSparseMatrix, x::Matrix, I::AbstractVector{T}, J::AbstractVector{Bool}) = setindex!(A, sparse(x, format=sptype(A)), I, find(J)) +setindex!{T<:Integer}(A::CompressedSparseMatrix, x::Matrix, I::AbstractVector{Bool}, J::AbstractVector{T}) = setindex!(A, sparse(x, format=sptype(A)), find(I),J) -setindex!(A::Matrix, x::SparseMatrixCSC, I::Integer, J::AbstractVector{Bool}) = setindex!(A, full(x), I, find(J)) -setindex!(A::Matrix, x::SparseMatrixCSC, I::AbstractVector{Bool}, J::Integer) = setindex!(A, full(x), find(I), J) -setindex!(A::Matrix, x::SparseMatrixCSC, I::AbstractVector{Bool}, J::AbstractVector{Bool}) = setindex!(A, full(x), find(I), find(J)) -setindex!{T<:Integer}(A::Matrix, x::SparseMatrixCSC, I::AbstractVector{T}, J::AbstractVector{Bool}) = setindex!(A, full(x), I, find(J)) -setindex!{T<:Integer}(A::Matrix, x::SparseMatrixCSC, I::AbstractVector{Bool}, J::AbstractVector{T}) = setindex!(A, full(x), find(I), J) +setindex!(A::Matrix, x::CompressedSparseMatrix, I::Integer, J::AbstractVector{Bool}) = setindex!(A, full(x), I, find(J)) +setindex!(A::Matrix, x::CompressedSparseMatrix, I::AbstractVector{Bool}, J::Integer) = setindex!(A, full(x), find(I), J) +setindex!(A::Matrix, x::CompressedSparseMatrix, I::AbstractVector{Bool}, J::AbstractVector{Bool}) = setindex!(A, full(x), find(I), find(J)) +setindex!{T<:Integer}(A::Matrix, x::CompressedSparseMatrix, I::AbstractVector{T}, J::AbstractVector{Bool}) = setindex!(A, full(x), I, find(J)) +setindex!{T<:Integer}(A::Matrix, x::CompressedSparseMatrix, I::AbstractVector{Bool}, J::AbstractVector{T}) = setindex!(A, full(x), find(I), J) -function setindex!{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, x, I::AbstractArray{Bool,2}) +function setindex!(A::CompressedSparseMatrix, x, I::AbstractArray{Bool,2}) checkbounds(A, I) - n = sum(I) + setindex!(A, x, find(I)) +end + +sortidx_colfirst{T<:Real}(I::AbstractVector{T}, sz::NTuple{2,Int}) = sortperm(I) +function sortidx_rowfirst{T<:Real}(I::AbstractVector{T}, sz::NTuple{2,Int}) + nr = sz[1] + NI = I .- 1 + l = length(NI) + @inbounds for i in 1:l + NI[i] = I[i] - nr * div(NI[i], nr) + end + sortperm(NI) +end + +function setindex!{Tv,Ti,T<:Real}(A::CompressedSparseMatrix{Tv,Ti}, x, I::AbstractVector{T}) + n = length(I) (n == 0) && (return A) - colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval - colptrB = colptrA; rowvalB = rowvalA; nzvalB = nzvalA + indptrB = indptrA = getindptr(A); indvalB = indvalA = getindval(A); nzvalB = nzvalA = A.nzval; + szA = size(A); caxis = spaxis(A); ncaxis = nspaxis(A) nadd = ndel = 0 - bidx = xidx = 1 - r1 = r2 = 0 - - for col in 1:A.n - r1 = int(colptrA[col]) - r2 = int(colptrA[col+1]-1) - - for row in 1:A.m - if I[row, col] - v = isa(x, AbstractArray) ? x[xidx] : x - xidx += 1 - - if r1 <= r2 - copylen = searchsortedfirst(rowvalA, row, r1, r2, Forward) - r1 - if (copylen > 0) - if (nadd > 0) || (ndel > 0) - copy!(rowvalB, bidx, rowvalA, r1, copylen) - copy!(nzvalB, bidx, nzvalA, r1, copylen) - end - bidx += copylen - r1 += copylen - end - end + bidx = aidx = 0 + + S = isa(A, SparseMatrixCSC) ? sortidx_colfirst(I, szA) : sortidx_rowfirst(I, szA) + sxidx = r1 = r2 = 0 + + n = length(I) + lastcidx = 0 + for xidx in 1:n + sxidx = S[xidx] + (sxidx < n) && (I[sxidx] == I[sxidx+1]) && continue + + subA = ind2sub(szA, I[sxidx]) + ncidxA = subA[ncaxis] + cidxA = subA[caxis] + v = isa(x, AbstractArray) ? x[sxidx] : x - # 0: no change, 1: update, 2: delete, 3: add new - mode = ((r1 <= r2) && (rowvalA[r1] == row)) ? ((v == 0) ? 2 : 1) : ((v == 0) ? 0 : 3) + if cidxA > lastcidx + r1 = int(indptrA[cidxA]) + r2 = int(indptrA[cidxA+1] - 1) - if (mode > 1) && (nadd == 0) && (ndel == 0) - # copy storage to take changes - colptrB = copy(colptrA) - memreq = (x == 0) ? 0 : n - rowvalB = Array(Ti, length(rowvalA)+memreq); copy!(rowvalB, 1, rowvalA, 1, r1-1) - nzvalB = Array(Tv, length(nzvalA)+memreq); copy!(nzvalB, 1, nzvalA, 1, r1-1) + # copy from last position till current column + if (nadd > 0) || (ndel > 0) + ndiff = nadd - ndel + @inbounds for i in (lastcidx+1):cidxA + indptrB[i] = indptrA[i] + ndiff end - if mode == 1 - rowvalB[bidx] = row - nzvalB[bidx] = v - bidx += 1 - r1 += 1 - elseif mode == 2 - r1 += 1 - ndel += 1 - elseif mode == 3 - rowvalB[bidx] = row - nzvalB[bidx] = v - bidx += 1 - nadd += 1 + copylen = r1 - aidx + if copylen > 0 + copy!(indvalB, bidx, indvalA, aidx, copylen) + copy!(nzvalB, bidx, nzvalA, aidx, copylen) + aidx += copylen + bidx += copylen end - (xidx > n) && break - end # if I[row, col] - end # for row in 1:A.m - - if ((nadd != 0) || (ndel != 0)) - l = r2-r1+1 - if l > 0 - copy!(rowvalB, bidx, rowvalA, r1, l) - copy!(nzvalB, bidx, nzvalA, r1, l) - bidx += l + else + aidx = bidx = r1 end - colptrB[col+1] = bidx - - if (xidx > n) && (length(colptrB) > (col+1)) - diff = nadd - ndel - colptrB[(col+2):end] = colptrA[(col+2):end] .+ diff - r1 = colptrA[col+1] - r2 = colptrA[end]-1 - l = r2-r1+1 - if l > 0 - copy!(rowvalB, bidx, rowvalA, r1, l) - copy!(nzvalB, bidx, nzvalA, r1, l) - bidx += l + lastcidx = cidxA + end + + if r1 <= r2 + copylen = searchsortedfirst(indvalA, ncidxA, r1, r2, Forward) - r1 + if copylen > 0 + if (nadd > 0) || (ndel > 0) + copy!(indvalB, bidx, indvalA, r1, copylen) + copy!(nzvalB, bidx, nzvalA, r1, copylen) end + bidx += copylen + r1 += copylen + aidx += copylen end - else - bidx = colptrA[col+1] end - (xidx > n) && break - end # for col in 1:A.n - if (nadd != 0) || (ndel != 0) + # 0: no change, 1: update, 2: delete, 3: add new + mode = ((r1 <= r2) && (indvalA[r1] == ncidxA)) ? ((v == 0) ? 2 : 1) : ((v == 0) ? 0 : 3) + + if (mode > 1) && (nadd == 0) && (ndel == 0) + # copy storage to take changes + indptrB = copy(indptrA) + memreq = (x == 0) ? 0 : n + indvalB = Array(Ti, length(indvalA)+memreq); copy!(indvalB, 1, indvalA, 1, r1-1) + nzvalB = Array(Tv, length(nzvalA)+memreq); copy!(nzvalB, 1, nzvalA, 1, r1-1) + end + if mode == 1 + indvalB[bidx] = ncidxA + nzvalB[bidx] = v + bidx += 1 + aidx += 1 + r1 += 1 + elseif mode == 2 + r1 += 1 + aidx += 1 + ndel += 1 + elseif mode == 3 + indvalB[bidx] = ncidxA + nzvalB[bidx] = v + bidx += 1 + nadd += 1 + end + end + + # copy the rest + if (nadd > 0) || (ndel > 0) + ndiff = nadd - ndel + @inbounds for i in (lastcidx+1):length(indptrA) + indptrB[i] = indptrA[i] + ndiff + end + r1 = indptrA[end]-1 + copylen = r1 - aidx + 1 + if copylen > 0 + copy!(indvalB, bidx, indvalA, aidx, copylen) + copy!(nzvalB, bidx, nzvalA, aidx, copylen) + aidx += copylen + bidx += copylen + end + n = length(nzvalB) if n > (bidx-1) deleteat!(nzvalB, bidx:n) - deleteat!(rowvalB, bidx:n) + deleteat!(indvalB, bidx:n) end - A.nzval = nzvalB; A.rowval = rowvalB; A.colptr = colptrB + A.nzval = nzvalB; setindval!(A, indvalB); setindptr!(A, indptrB) end A end - +#= function setindex!{Tv,Ti,T<:Real}(A::SparseMatrixCSC{Tv,Ti}, x, I::AbstractVector{T}) n = length(I) (n == 0) && (return A) @@ -1609,15 +1852,14 @@ function setindex!{Tv,Ti,T<:Real}(A::SparseMatrixCSC{Tv,Ti}, x, I::AbstractVecto end A end - - +=# # Sparse concatenation -function vcat(X::SparseMatrixCSC...) +function cat_along_ncaxis(X::CompressedSparseMatrix...) num = length(X) - mX = [ size(x, 1) for x in X ] - nX = [ size(x, 2) for x in X ] + mX = [ nspdim(x) for x in X ] + nX = [ spdim(x) for x in X ] n = nX[1] for i = 2 : num if nX[i] != n; throw(DimensionMismatch("")); end @@ -1625,37 +1867,46 @@ function vcat(X::SparseMatrixCSC...) m = sum(mX) Tv = promote_type(map(x->eltype(x.nzval), X)...) - Ti = promote_type(map(x->eltype(x.rowval), X)...) + Ti = promote_type(map(x->eltype(getindval(x)), X)...) - colptr = Array(Ti, n + 1) + iptr = Array(Ti, n + 1) nnzX = [ nnz(x) for x in X ] nnz_res = sum(nnzX) - rowval = Array(Ti, nnz_res) + ival = Array(Ti, nnz_res) nzval = Array(Tv, nnz_res) - colptr[1] = 1 + iptr[1] = 1 for c = 1 : n mX_sofar = 0 - rr1 = colptr[c] + rr1 = iptr[c] for i = 1 : num - rX1 = X[i].colptr[c] - rX2 = X[i].colptr[c + 1] - 1 + rX1 = getindptr(X[i])[c] + rX2 = getindptr(X[i])[c + 1] - 1 rr2 = rr1 + (rX2 - rX1) - rowval[rr1 : rr2] = X[i].rowval[rX1 : rX2] .+ mX_sofar + ival[rr1 : rr2] = getindval(X[i])[rX1 : rX2] .+ mX_sofar nzval[rr1 : rr2] = X[i].nzval[rX1 : rX2] mX_sofar += mX[i] rr1 = rr2 + 1 end - colptr[c + 1] = rr1 + iptr[c + 1] = rr1 end - SparseMatrixCSC(m, n, colptr, rowval, nzval) + (m, n, iptr, ival, nzval) end -function hcat(X::SparseMatrixCSC...) +function vcat(X::SparseMatrixCSC...) + (ncdim, cdim, iptr, ival, nzval) = cat_along_ncaxis(X...) + SparseMatrixCSC(ncdim, cdim, iptr, ival, nzval) +end +function hcat(X::SparseMatrixCSR...) + (ncdim, cdim, iptr, ival, nzval) = cat_along_ncaxis(X...) + SparseMatrixCSR(cdim, ncdim, iptr, ival, nzval) +end + +function cat_along_caxis(X::CompressedSparseMatrix...) num = length(X) - mX = [ size(x, 1) for x in X ] - nX = [ size(x, 2) for x in X ] + mX = [ nspdim(x) for x in X ] + nX = [ spdim(x) for x in X ] m = mX[1] for i = 2 : num if mX[i] != m; throw(DimensionMismatch("")); end @@ -1663,31 +1914,41 @@ function hcat(X::SparseMatrixCSC...) n = sum(nX) Tv = promote_type(map(x->eltype(x.nzval), X)...) - Ti = promote_type(map(x->eltype(x.rowval), X)...) + Ti = promote_type(map(x->eltype(getindval(x)), X)...) - colptr = Array(Ti, n + 1) + iptr = Array(Ti, n + 1) nnzX = [ nnz(x) for x in X ] nnz_res = sum(nnzX) - rowval = Array(Ti, nnz_res) + ival = Array(Ti, nnz_res) nzval = Array(Tv, nnz_res) nnz_sofar = 0 nX_sofar = 0 for i = 1 : num - colptr[(1 : nX[i] + 1) + nX_sofar] = X[i].colptr .+ nnz_sofar - rowval[(1 : nnzX[i]) + nnz_sofar] = X[i].rowval + iptr[(1 : nX[i] + 1) + nX_sofar] = getindptr(X[i]) .+ nnz_sofar + ival[(1 : nnzX[i]) + nnz_sofar] = getindval(X[i]) nzval[(1 : nnzX[i]) + nnz_sofar] = X[i].nzval nnz_sofar += nnzX[i] nX_sofar += nX[i] end - SparseMatrixCSC(m, n, colptr, rowval, nzval) + (m, n, iptr, ival, nzval) +end + +function hcat(X::SparseMatrixCSC...) + (ncdim, cdim, iptr, ival, nzval) = cat_along_caxis(X...) + SparseMatrixCSC(ncdim, cdim, iptr, ival, nzval) +end + +function vcat(X::SparseMatrixCSR...) + (ncdim, cdim, iptr, ival, nzval) = cat_along_caxis(X...) + SparseMatrixCSR(cdim, ncdim, iptr, ival, nzval) end -function hvcat(rows::(Int...), X::SparseMatrixCSC...) +function hvcat{T<:CompressedSparseMatrix}(rows::(Int...), X::T...) nbr = length(rows) # number of block rows - tmp_rows = Array(SparseMatrixCSC, nbr) + tmp_rows = Array(T, nbr) k = 0 for i = 1 : nbr tmp_rows[i] = hcat(X[(1 : rows[i]) + k]...) @@ -1696,59 +1957,66 @@ function hvcat(rows::(Int...), X::SparseMatrixCSC...) vcat(tmp_rows...) end -function blkdiag(X::SparseMatrixCSC...) +function blkdiag{T<:CompressedSparseMatrix}(X::T...) num = length(X) - mX = [ size(x, 1) for x in X ] - nX = [ size(x, 2) for x in X ] + mX = [ nspdim(x) for x in X ] + nX = [ spdim(x) for x in X ] m = sum(mX) n = sum(nX) Tv = promote_type(map(x->eltype(x.nzval), X)...) - Ti = promote_type(map(x->eltype(x.rowval), X)...) + Ti = promote_type(map(x->eltype(getindval(x)), X)...) - colptr = Array(Ti, n + 1) + iptr = Array(Ti, n + 1) nnzX = [ nnz(x) for x in X ] nnz_res = sum(nnzX) - rowval = Array(Ti, nnz_res) + ival = Array(Ti, nnz_res) nzval = Array(Tv, nnz_res) nnz_sofar = 0 nX_sofar = 0 mX_sofar = 0 for i = 1 : num - colptr[(1 : nX[i] + 1) + nX_sofar] = X[i].colptr .+ nnz_sofar - rowval[(1 : nnzX[i]) + nnz_sofar] = X[i].rowval .+ mX_sofar + iptr[(1 : nX[i] + 1) + nX_sofar] = getindptr(X[i]) .+ nnz_sofar + ival[(1 : nnzX[i]) + nnz_sofar] = getindval(X[i]) .+ mX_sofar nzval[(1 : nnzX[i]) + nnz_sofar] = X[i].nzval nnz_sofar += nnzX[i] nX_sofar += nX[i] mX_sofar += mX[i] end - SparseMatrixCSC(m, n, colptr, rowval, nzval) + if T <: SparseMatrixCSC + SparseMatrixCSC(m, n, iptr, ival, nzval) + else + SparseMatrixCSR(n, m, iptr, ival, nzval) + end end ## Structure query functions -function issym(A::SparseMatrixCSC) +function issym(A::CompressedSparseMatrix) m, n = size(A) if m != n; return false; end return countnz(A - A.') == 0 end -function ishermitian(A::SparseMatrixCSC) +function ishermitian(A::CompressedSparseMatrix) m, n = size(A) if m != n; return false; end return countnz(A - A') == 0 end -function istriu{Tv}(A::SparseMatrixCSC{Tv}) - for col = 1:min(A.n,A.m-1) - l1 = A.colptr[col+1]-1 - for i = 0 : (l1 - A.colptr[col]) - if A.rowval[l1-i] <= col +function istric{Tv}(A::CompressedSparseMatrix{Tv}) + iptr = getindptr(A) + ival = getindval(A) + nzval = A.nzval + for spidx = 1:min(spdim(A),nspdim(A)-1) + l1 = iptr[spidx+1]-1 + for i = 0 : (l1 - iptr[spidx]) + if ival[l1-i] <= spidx break end - if A.nzval[l1-i] != 0 + if nzval[l1-i] != 0 return false end end @@ -1756,13 +2024,16 @@ function istriu{Tv}(A::SparseMatrixCSC{Tv}) return true end -function istril{Tv}(A::SparseMatrixCSC{Tv}) - for col = 2:A.n - for i = A.colptr[col] : (A.colptr[col+1]-1) - if A.rowval[i] >= col +function istrinc{Tv}(A::CompressedSparseMatrix{Tv}) + iptr = getindptr(A) + ival = getindval(A) + nzval = A.nzval + for spidx = 2:spdim(A) + for i = iptr[spidx] : (iptr[spidx+1]-1) + if ival[i] >= spidx break end - if A.nzval[i] != 0 + if nzval[i] != 0 return false end end @@ -1770,6 +2041,11 @@ function istril{Tv}(A::SparseMatrixCSC{Tv}) return true end +istriu{Tv}(A::SparseMatrixCSC{Tv}) = istric(A) +istril{Tv}(A::SparseMatrixCSC{Tv}) = istrinc(A) +istriu{Tv}(A::SparseMatrixCSR{Tv}) = istrinc(A) +istril{Tv}(A::SparseMatrixCSR{Tv}) = istric(A) + # Create a sparse diagonal matrix by specifying multiple diagonals # packed into a tuple, alongside their diagonal offsets and matrix shape @@ -1823,7 +2099,7 @@ spdiagm(B::AbstractVector, d::Number, m::Integer, n::Integer) = spdiagm((B,), (d spdiagm(B::AbstractVector, d::Number=0) = spdiagm((B,), (d,)) -## expand a colptr or rowptr into a dense index vector +## expand a indptr or rowptr into a dense index vector function expandptr{T<:Integer}(V::Vector{T}) if V[1] != 1 throw(ArgumentError("first index must be one")) end res = similar(V, (int64(V[end]-1),)) @@ -1834,25 +2110,25 @@ end ## diag and related using an iterator type SpDiagIterator{Tv,Ti} - A::SparseMatrixCSC{Tv,Ti} + A::CompressedSparseMatrix{Tv,Ti} n::Int end -SpDiagIterator(A::SparseMatrixCSC) = SpDiagIterator(A,minimum(size(A))) +SpDiagIterator(A::CompressedSparseMatrix) = SpDiagIterator(A,minimum(size(A))) length(d::SpDiagIterator) = d.n start(d::SpDiagIterator) = 1 done(d::SpDiagIterator, j) = j > d.n function next{Tv}(d::SpDiagIterator{Tv}, j) - A = d.A - r1 = int(A.colptr[j]) - r2 = int(A.colptr[j+1]-1) + p = getindptr(d.A); i = getindval(d.A); + r1 = int(p[j]) + r2 = int(p[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) + r1 = searchsortedfirst(i, j, r1, r2, Forward) + (((r1 > r2) || (i[r1] != j)) ? zero(Tv) : d.A.nzval[r1], j+1) end -function trace{Tv}(A::SparseMatrixCSC{Tv}) +function trace{Tv}(A::CompressedSparseMatrix{Tv}) if size(A,1) != size(A,2) throw(DimensionMismatch("expected square matrix")) end @@ -1863,48 +2139,48 @@ function trace{Tv}(A::SparseMatrixCSC{Tv}) s end -diag(A::SparseMatrixCSC) = [d for d in SpDiagIterator(A)] +diag(A::CompressedSparseMatrix) = [d for d in SpDiagIterator(A)] -function diagm{Tv,Ti}(v::SparseMatrixCSC{Tv,Ti}) +function diagm{Tv,Ti}(v::CompressedSparseMatrix{Tv,Ti}) if (size(v,1) != 1 && size(v,2) != 1) throw(DimensionMismatch("input should be nx1 or 1xn")) end n = length(v) numnz = nnz(v) - colptr = Array(Ti, n+1) - rowval = Array(Ti, numnz) + iptr = Array(Ti, n+1) + ival = Array(Ti, numnz) nzval = Array(Tv, numnz) if size(v,1) == 1 - copy!(colptr, 1, v.colptr, 1, n+1) + copy!(iptr, 1, getindptr(v), 1, n+1) ptr = 1 for col = 1:n - if colptr[col] != colptr[col+1] - rowval[ptr] = col + if iptr[col] != iptr[col+1] + ival[ptr] = col nzval[ptr] = v.nzval[ptr] ptr += 1 end end else - copy!(rowval, 1, v.rowval, 1, numnz) + copy!(ival, 1, getindval(v), 1, numnz) copy!(nzval, 1, v.nzval, 1, numnz) - colptr[1] = 1 + iptr[1] = 1 ptr = 1 col = 1 while col <= n && ptr <= numnz - while rowval[ptr] > col - colptr[col+1] = colptr[col] + while ival[ptr] > col + iptr[col+1] = iptr[col] col += 1 end - colptr[col+1] = colptr[col] + 1 + iptr[col+1] = iptr[col] + 1 ptr += 1 col += 1 end if col <= n - colptr[(col+1):(n+1)] = colptr[col] + iptr[(col+1):(n+1)] = iptr[col] end end - return SparseMatrixCSC{Tv,Ti}(n, n, colptr, rowval, nzval) + sparse(n, n, CompressedSparseStore(iptr, ival, nzval), format=sptype(v)) end diff --git a/test/sparse.jl b/test/sparse.jl index 5205120164e18..b94806dc6b003 100644 --- a/test/sparse.jl +++ b/test/sparse.jl @@ -1,5 +1,6 @@ # check sparse matrix construction -@test isequal(full(sparse(complex(ones(5,5),ones(5,5)))), complex(ones(5,5),ones(5,5))) +@test isequal(full(sparse(complex(ones(5,5),ones(5,5)), format=CSC)), complex(ones(5,5),ones(5,5))) +@test isequal(full(sparse(complex(ones(5,5),ones(5,5)), format=CSR)), complex(ones(5,5),ones(5,5))) # check matrix operations se33 = speye(3) @@ -332,31 +333,50 @@ let A = spzeros(Int, 10, 20) @test A[4:8,8:16] == 15 * ones(Int, 5, 9) end -let ASZ = 1000, TSZ = 800 - A = sprand(ASZ, 2*ASZ, 0.0001) - B = copy(A) - nA = countnz(A) - x = A[1:TSZ, 1:(2*TSZ)] - nx = countnz(x) - A[1:TSZ, 1:(2*TSZ)] = 0 - nB = countnz(A) - @test nB == (nA - nx) - A[1:TSZ, 1:(2*TSZ)] = x - @test countnz(A) == nA - @test A == B - A[1:TSZ, 1:(2*TSZ)] = 10 - @test countnz(A) == nB + 2*TSZ*TSZ - A[1:TSZ, 1:(2*TSZ)] = x - @test countnz(A) == nA - @test A == B +for Ts in [CSC, CSR] + let ASZ = 1000, TSZ = 800 + A = sprand(ASZ, 2*ASZ, 0.0001, format=Ts) + B = copy(A) + nA = countnz(A) + x = A[1:TSZ, 1:(2*TSZ)] + nx = countnz(x) + A[1:TSZ, 1:(2*TSZ)] = 0 + nB = countnz(A) + @test nB == (nA - nx) + A[1:TSZ, 1:(2*TSZ)] = x + @test countnz(A) == nA + @test A == B + A[1:TSZ, 1:(2*TSZ)] = 10 + @test countnz(A) == nB + 2*TSZ*TSZ + A[1:TSZ, 1:(2*TSZ)] = x + @test countnz(A) == nA + @test A == B + end end -let A = speye(Int, 5), I=[1:10], X=reshape([trues(10), falses(15)],5,5) +for Ts in [CSC, CSR] + A = speye(Int, 5, format=Ts) + I = [1:10] + X = reshape([trues(10), falses(15)],5,5) @test A[I] == A[X] == reshape([1,0,0,0,0,0,1,0,0,0], 10, 1) A[I] = [1:10] @test A[I] == A[X] == reshape([1:10], 10, 1) end +# hcat, vcat, hvcat test +for Ts in [CSC, CSR] + (m,n) = (3,2) + spm = {} + for i in 1:m + push!(spm, sprand(4, 2, 0.5, format=Ts)) + end + for i in 1:n + push!(spm, sprand(2, 3, 0.5, format=Ts)) + end + + @test vcat(hcat(spm[1:m]...), hcat(spm[(m+1):(m+n)]...)) == hvcat((m,n), spm...) +end + let S = sprand(50, 30, 0.5, x->int(rand(x)*100)), I = sprandbool(50, 30, 0.2) FS = full(S) FI = full(I)