Skip to content

Commit

Permalink
Replace the LGPL-licensed SparseMatrixCSC transposition methods in ba…
Browse files Browse the repository at this point in the history 
…se/sparse/csparse.jl ([c|f]transpose[!]) with MIT-licensed versions. See #13001.
  • Loading branch information
@Sacha0
Sacha0 committed Jan 11, 2016
1 parent bcfc967 commit e2efd22
Show file tree
Hide file tree
Showing 2 changed files with 95 additions and 63 deletions.
63 changes: 0 additions & 63 deletions base/sparse/csparse.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
95 changes: 95 additions & 0 deletions base/sparse/sparsematrix.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -350,6 +350,101 @@ 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
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)"))
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)"))
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)"))
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)"))
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)"))

# 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 additional 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)
Expand Down

0 comments on commit e2efd22

Please sign in to comment.