Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Replace the LGPL-licensed SparseMatrixCSC transposition methods in ba…
…se/sparse/csparse.jl ([c|f]transpose[!]) with MIT-licensed versions. See JuliaLang#13001.
- Loading branch information
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -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)")) | ||
This comment has been minimized.
Sorry, something went wrong.
This comment has been minimized.
Sorry, something went wrong. 
tkelman
|
||
|
|
||
| # 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) | ||
|
|
||
Should I wrap these longer lines somehow, to observe the 92-character line-length guideline? If so, what is the preferred format?