Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

MIT-licensed SparseMatrixCSC transposition methods #14631

Merged
merged 1 commit into from
Jan 13, 2016
+114 −63
Merged

MIT-licensed SparseMatrixCSC transposition methods #14631

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 0 additions & 63 deletions base/sparse/csparse.jl
View file
Open in desktop
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
114 changes: 114 additions & 0 deletions base/sparse/sparsematrix.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down