From ab0bdf2313d0056498b95c8ffd234f68ad90f76a Mon Sep 17 00:00:00 2001 From: Sacha Verweij Date: Sun, 10 Jan 2016 16:34:07 -0800 Subject: [PATCH] Replace the LGPL-licensed SparseMatrixCSC transposition methods in base/sparse/csparse.jl ([c|f]transpose[!]) with MIT-licensed versions. See #13001. --- base/sparse/csparse.jl | 63 -------------------- base/sparse/sparsematrix.jl | 114 ++++++++++++++++++++++++++++++++++++ 2 files changed, 114 insertions(+), 63 deletions(-) diff --git a/base/sparse/csparse.jl b/base/sparse/csparse.jl index 0a62afda659f5..a9b251eea5b5b 100644 --- a/base/sparse/csparse.jl +++ b/base/sparse/csparse.jl @@ -138,69 +138,6 @@ function sparse{Tv,Ti<:Integer}(I::AbstractVector{Ti}, return SparseMatrixCSC(nrow, ncol, RpT, RiT, RxT) end -## Transpose and apply f - -# Based on Direct Methods for Sparse Linear Systems, T. A. Davis, SIAM, Philadelphia, Sept. 2006. -# Section 2.5: Transpose -# http://www.cise.ufl.edu/research/sparse/CSparse/ -function ftranspose!{Tv,Ti}(T::SparseMatrixCSC{Tv,Ti}, S::SparseMatrixCSC{Tv,Ti}, f) - (mS, nS) = size(S) - nnzS = nnz(S) - colptr_S = S.colptr - rowval_S = S.rowval - nzval_S = S.nzval - - (mT, nT) = size(T) - colptr_T = T.colptr - rowval_T = T.rowval - nzval_T = T.nzval - - fill!(colptr_T, 0) - colptr_T[1] = 1 - for i=1:nnzS - @inbounds colptr_T[rowval_S[i]+1] += 1 - end - cumsum!(colptr_T, colptr_T) - - w = copy(colptr_T) - @inbounds for j = 1:nS, p = colptr_S[j]:(colptr_S[j+1]-1) - ind = rowval_S[p] - q = w[ind] - w[ind] += 1 - rowval_T[q] = j - nzval_T[q] = f(nzval_S[p]) - end - - return T -end - -function ftranspose{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, f) - (nT, mT) = size(S) - nnzS = nnz(S) - colptr_T = Array(Ti, nT+1) - rowval_T = Array(Ti, nnzS) - nzval_T = Array(Tv, nnzS) - - T = SparseMatrixCSC(mT, nT, colptr_T, rowval_T, nzval_T) - return ftranspose!(T, S, f) -end - -function transpose!{Tv,Ti}(T::SparseMatrixCSC{Tv,Ti}, S::SparseMatrixCSC{Tv,Ti}) - ftranspose!(T, S, IdFun()) -end - -function transpose{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}) - ftranspose(S, IdFun()) -end - -function ctranspose!{Tv,Ti}(T::SparseMatrixCSC{Tv,Ti}, S::SparseMatrixCSC{Tv,Ti}) - ftranspose!(T, S, ConjFun()) -end - -function ctranspose{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}) - ftranspose(S, ConjFun()) -end - # Compute the elimination tree of A using triu(A) returning the parent vector. # A root node is indicated by 0. This tree may actually be a forest in that # there may be more than one root, indicating complete separability. diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index fb1a5cd516a74..b4b976d4b89b6 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -350,6 +350,120 @@ function sparse(B::Bidiagonal) return sparse([1:m;1:m-1],[1:m;2:m],[B.dv;B.ev], Int(m), Int(m)) # upper bidiagonal end +## Transposition methods + +# qftranspose! is the parent method on which the others are built. +""" + qftranspose!{Tv,Ti}(C::SparseMatrixCSC{Tv,Ti}, A::SparseMatrixCSC{Tv,Ti}, q::AbstractVector, f) + +Column-permute and transpose `A` (`(AQ)^T`), applying `f` to each element of `A` in the + process, and store the result in preallocated `C`. Permutation vector `q` defines the + column-permutation `Q`. The number of columns of `C` (`C.n`) must match the number of + rows of `A` (`A.m`). The number of rows of `C` (`C.m`) must match the number of columns + of `A` (`A.n`). The length of `C`'s internal row-index (`length(C.rowval)`) and + entry-value (`length(C.nzval)`) arrays must be at least the number of allocated entries + in `A` (`nnz(A)`). The length of the permutation vector `q` (`length(q)`) must match the + number of columns of `A` (`A.n`). + +This method implements the HALFPERM algorithm described in F. Gustavson, "Two fast + algorithms for sparse matrices: multiplication and permuted transposition," ACM TOMS + 4(3), 250-269 (1978). The algorithm runs in `O(A.m, A.n, nnz(A))` time and requires no + space beyond that passed in. + +Performance note: As of January 2016, `f` should be a functor for this method to perform + well. This caveat may disappear when the work in `jb/functions` lands. +""" +function qftranspose!{Tv,Ti}(C::SparseMatrixCSC{Tv,Ti}, A::SparseMatrixCSC{Tv,Ti}, q::AbstractVector, f) + # Attach source matrix + Am, An = A.m, A.n + Acolptr = A.colptr + Arowval = A.rowval + Anzval = A.nzval + Annz = Acolptr[end]-1 + + # Attach destination matrix + Cm, Cn = C.m, C.n + Ccolptr = C.colptr + Crowval = C.rowval + Cnzval = C.nzval + + # Check compatibility of source and destination + if !(Cm == An) + throw(DimensionMismatch("the number of rows of the first argument, C.m = $(Cm)," + * "must match the number of columns of the second argument, A.n = $(An)") ) + elseif !(Cn == Am) + throw(DimensionMismatch("the number of columns of the first argument, C.n = $(Cn)," + * "must match the number of rows of the second argument, A.m = $(Am)") ) + elseif !(length(q) == A.n) + throw(DimensionMismatch("the length of the permtuation vector, length(q) = " + * "$(length(q)), must match the number of columns of the second argument," + * "A.n = $(An)") ) + elseif !(length(Crowval) >= Annz) + throw(ArgumentError("the first argument's row-index array's length," + * "length(C.rowval) = $(length(Crowval)), must be at least the number of" + * "allocated entries in the second argument, nnz(A) = $(Annz)") ) + elseif !(length(Cnzval) >= Annz) + throw(ArgumentError("the first argument's entry-value array's length," + * "length(C.nzval) = $(length(Cnzval)), must be at least the number of" + * "allocated entries in the second argument, nnz(A) = $(Annz)") ) + end + + # Compute the column counts of C and store them shifted forward by one in Ccolptr + Ccolptr[1:end] = 0 + @inbounds for k in 1:Annz + Ccolptr[Arowval[k]+1] += 1 + end + + # From these column counts, compute C's column pointers + # and store them shifted forward by one in Ccolptr + countsum = 1 + @inbounds for k in 2:(Cn+1) + overwritten = Ccolptr[k] + Ccolptr[k] = countsum + countsum += overwritten + end + + # Distribution-sort the row indices and nonzero values into Crowval and Cnzval, + # tracking write positions in Ccolptr + @inbounds for Aj in 1:An + qAj = q[Aj] + for Ak in Acolptr[qAj]:(Acolptr[qAj+1]-1) + Ai = Arowval[Ak] + Ck = Ccolptr[Ai+1] + Crowval[Ck] = qAj + Cnzval[Ck] = f(Anzval[Ak]) + Ccolptr[Ai+1] += 1 + end + end + + # Tracking write positions in Ccolptr as in the last block fixes the colptr shift, + # but the first colptr remains incorrect + Ccolptr[1] = 1 + + C +end +transpose!{Tv,Ti}(C::SparseMatrixCSC{Tv,Ti}, A::SparseMatrixCSC{Tv,Ti}) = qftranspose!(C, A, 1:A.n, Base.IdFun()) +ctranspose!{Tv,Ti}(C::SparseMatrixCSC{Tv,Ti}, A::SparseMatrixCSC{Tv,Ti}) = qftranspose!(C, A, 1:A.n, Base.ConjFun()) +"See `qftranspose!`" ftranspose!{Tv,Ti}(C::SparseMatrixCSC{Tv,Ti}, A::SparseMatrixCSC{Tv,Ti}, f) = qftranspose!(C, A, 1:A.n, f) + +""" + qftranspose{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, q::AbstractVector, f) + +Return-allocating version of `qftranspose!`. See `qftranspose!` for documentation. +""" +function qftranspose{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, q::AbstractVector, f) + Cm, Cn, Cnnz = A.n, A.m, nnz(A) + Ccolptr = zeros(Ti, Cn+1) + Crowval = Array{Ti}(Cnnz) + Cnzval = Array{Tv}(Cnnz) + qftranspose!(SparseMatrixCSC(Cm, Cn, Ccolptr, Crowval, Cnzval), A, q, f) +end +transpose{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}) = qftranspose(A, 1:A.n, Base.IdFun()) +ctranspose{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}) = qftranspose(A, 1:A.n, Base.ConjFun()) +"See `qftranspose`" ftranspose{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, f) = qftranspose(A, 1:A.n, f) + +## Find methods + function find(S::SparseMatrixCSC) sz = size(S) I, J = findn(S)