Replace the LGPL-licensed SparseMatrixCSC fkeep! method and children …

…tril!, triu!, droptol!, and dropzeros[!] with MIT-licensed versions. See JuliaLang#13001 and JuliaLang#14631.
Sacha0 · Jan 17, 2016 · 0d970aa · 0d970aa
1 parent a1fc070
commit 0d970aa
Show file tree

Hide file tree

Showing 2 changed files with 67 additions and 61 deletions.
diff --git a/base/sparse/csparse.jl b/base/sparse/csparse.jl
@@ -299,64 +299,4 @@ function symperm{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti})
         end
     end
     (C.').' # double transpose to order the columns
-end
-
-# Based on Direct Methods for Sparse Linear Systems, T. A. Davis, SIAM, Philadelphia, Sept. 2006.
-# Section 2.7: Removing entries from a matrix
-# http://www.cise.ufl.edu/research/sparse/CSparse/
-function fkeep!{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, f, other)
-    nzorig = nnz(A)
-    nz = 1
-    colptr = A.colptr
-    rowval = A.rowval
-    nzval = A.nzval
-    @inbounds for j = 1:A.n
-        p = colptr[j]                 # record current position
-        colptr[j] = nz                # set new position
-        while p < colptr[j+1]
-            if f(rowval[p], j, nzval[p], other)
-                nzval[nz] = nzval[p]
-                rowval[nz] = rowval[p]
-                nz += 1
-            end
-            p += 1
-        end
-    end
-    colptr[A.n + 1] = nz
-    nz -= 1
-    if nz < nzorig
-        resize!(nzval, nz)
-        resize!(rowval, nz)
-    end
-    A
-end
-
-
-immutable DropTolFun <: Func{4} end
-call(::DropTolFun, i,j,x,other) = abs(x)>other
-immutable DropZerosFun <: Func{4} end
-call(::DropZerosFun, i,j,x,other) = x!=0
-immutable TriuFun <: Func{4} end
-call(::TriuFun, i,j,x,other) = j>=i + other
-immutable TrilFun <: Func{4} end
-call(::TrilFun, i,j,x,other) = i>=j - other
-
-droptol!(A::SparseMatrixCSC, tol) = fkeep!(A, DropTolFun(), tol)
-dropzeros!(A::SparseMatrixCSC) = fkeep!(A, DropZerosFun(), nothing)
-dropzeros(A::SparseMatrixCSC) = dropzeros!(copy(A))
-
-function triu!(A::SparseMatrixCSC, k::Integer=0)
-    m,n = size(A)
-    if (k > 0 && k > n) || (k < 0 && -k > m)
-        throw(BoundsError())
-    end
-    fkeep!(A, TriuFun(), k)
-end
-
-function tril!(A::SparseMatrixCSC, k::Integer=0)
-    m,n = size(A)
-    if (k > 0 && k > n) || (k < 0 && -k > m)
-        throw(BoundsError())
-    end
-    fkeep!(A, TrilFun(), k)
-end
+end
diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl
@@ -462,6 +462,72 @@ 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)
 
+
+## fkeep! and children tril!, triu!, droptol!, dropzeros[!]
+
+"""
+    fkeep!{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, f, other)
+
+Keep elements of `A` for which test `f` returns `true`. `f`'s signature should be
+
+    f{Tv,Ti}(i::Ti, j::Ti, x::Tv, other::Any) -> Bool
+
+  where `i` and `j` are an element's row and column indices, `x` is the element's value,
+  and `other` is passed in from the call to `fkeep!`. This method makes a single sweep
+  through `A`, requiring `O(A.n, nnz(A))`-time and 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 fkeep!{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, f, other)
+    # Attach source matrix
+    Am, An = A.m, A.n
+    Acolptr = A.colptr
+    Arowval = A.rowval
+    Anzval = A.nzval
+    Annz = Acolptr[end]-1
+
+    # Sweep through columns, rewriting kept elements in their new positions and updating
+    # the column pointers accordingly as we go.
+    Awritepos = 1
+    oldAcolptrAj = 1
+    @inbounds for Aj in 1:An
+        for Ak in oldAcolptrAj:(Acolptr[Aj+1]-1)
+            Ai = Arowval[Ak]
+            Ax = Anzval[Ak]
+            # If this element should be kept, rewrite in new position
+            if f(Ai, Aj, Ax, other)
+                if Awritepos != Ak
+                    Arowval[Awritepos] = Ai
+                    Anzval[Awritepos] = Ax
+                end
+                Awritepos += 1
+            end
+        end
+        oldAcolptrAj = Acolptr[Aj+1]
+        Acolptr[Aj+1] = Awritepos
+    end
+
+    A
+end
+
+immutable TrilFunc <: Base.Func{4} end
+immutable TriuFunc <: Base.Func{4} end
+call{Tv,Ti}(::TrilFunc, i::Ti, j::Ti, x::Tv, k::Integer) = i + k >= j
+call{Tv,Ti}(::TriuFunc, i::Ti, j::Ti, x::Tv, k::Integer) = j >= i + k
+tril!(A::SparseMatrixCSC, k::Integer = 0) = fkeep!(A, TrilFunc(), k)
+triu!(A::SparseMatrixCSC, k::Integer = 0) = fkeep!(A, TriuFunc(), k)
+
+immutable DroptolFunc <: Base.Func{4} end
+call{Tv,Ti}(::DroptolFunc, i::Ti, j::Ti, x::Tv, tol::Real) = abs(x) >= tol
+droptol!(A::SparseMatrixCSC, tol) = fkeep!(A, DroptolFunc(), tol)
+
+immutable DropzerosFunc <: Base.Func{4} end
+call{Tv,Ti}(::DropzerosFunc, i::Ti, j::Ti, x::Tv, other) = x != 0
+dropzeros!(A::SparseMatrixCSC) = fkeep!(A, DropzerosFunc(), Void)
+dropzeros(A::SparseMatrixCSC) = dropzeros!(copy(A))
+
+
 ## Find methods
 
 function find(S::SparseMatrixCSC)