diff --git a/LICENSE.md b/LICENSE.md index 1f00b208885a3..f15371de45da1 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -65,7 +65,6 @@ their own licenses: The following components of Julia's standard library have separate licenses: - base/fftw.jl (see [FFTW](http://fftw.org/doc/License-and-Copyright.html)) -- base/sparse/csparse.jl (LGPL-2.1+) - base/linalg/umfpack.jl (see [SUITESPARSE](http://faculty.cse.tamu.edu/davis/suitesparse.html)) - base/linalg/cholmod.jl (see [SUITESPARSE](http://faculty.cse.tamu.edu/davis/suitesparse.html)) diff --git a/base/sparse.jl b/base/sparse.jl index 9c7a4cc446558..5c5fedd8929e9 100644 --- a/base/sparse.jl +++ b/base/sparse.jl @@ -36,7 +36,6 @@ export AbstractSparseArray, AbstractSparseMatrix, AbstractSparseVector, include("sparse/abstractsparse.jl") include("sparse/sparsematrix.jl") include("sparse/sparsevector.jl") -include("sparse/csparse.jl") include("sparse/linalg.jl") if Base.USE_GPL_LIBS @@ -45,4 +44,56 @@ if Base.USE_GPL_LIBS include("sparse/spqr.jl") end +# point users to SuiteSparse +const SUITESPARSE_END_STRING = string(" has been moved to the package SuiteSparse.jl.\n", + "Run Pkg.add(\"SuiteSparse\") to install SuiteSparse on Julia v0.5-") + +""" + etree(A[, post]) +Compute the elimination tree of a symmetric sparse matrix `A` from `triu(A)` +and, optionally, its post-ordering permutation. +Note: This function has been moved to the SuiteSparse.jl package. +""" +function etree{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, postorder::Bool) + if isdefined(Main, :SuiteSparse) + Main.SuiteSparse.etree(A, postorder) + else + error("etree(A[, post])", SUITESPARSE_END_STRING) + end +end + +etree(A::SparseMatrixCSC) = etree(A, false) + +function ereach{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, k::Integer, parent::Vector{Ti}) + if isdefined(Main, :SuiteSparse) + Main.SuiteSparse.ereach(A, k, parent) + else + error("ereach(A, k, parent)", SUITESPARSE_END_STRING) + end +end + +function csc_permute{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti}, q::Vector{Ti}) + if isdefined(Main, :SuiteSparse) + Main.SuiteSparse.csc_permute(A, pinv, q) + else + error("csc_permute(A, pinv, q)", SUITESPARSE_END_STRING) + end +end + +""" + symperm(A, p) +Return the symmetric permutation of `A`, which is `A[p,p]`. `A` should be +symmetric, sparse, and only contain nonzeros in the upper triangular part of the +matrix is stored. This algorithm ignores the lower triangular part of the +matrix. Only the upper triangular part of the result is returned. +Note: This function has been moved to the SuiteSparse.jl package. +""" +function symperm{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti}) + if isdefined(Main, :SuiteSparse) + Main.SuiteSparse.symperm(A, pinv) + else + error("symperm(A, pinv)", SUITESPARSE_END_STRING) + end +end + end diff --git a/base/sparse/csparse.jl b/base/sparse/csparse.jl deleted file mode 100644 index 53a153c8aa77b..0000000000000 --- a/base/sparse/csparse.jl +++ /dev/null @@ -1,178 +0,0 @@ -# These functions are based on C functions in the CSparse library by Tim Davis. -# These are pure Julia implementations, and do not link to the CSparse library. -# CSparse can be downloaded from http://www.cise.ufl.edu/research/sparse/CSparse/CSparse.tar.gz -# CSparse is Copyright (c) 2006-2007, Timothy A. Davis and released under -# Lesser GNU Public License, version 2.1 or later. A copy of the license can be -# downloaded from http://www.gnu.org/licenses/lgpl-2.1.html - -# Because these functions are based on code covered by LGPL-2.1+ the same license -# must apply to the code in this file which is -# Copyright (c) 2013-2014 Viral Shah, Douglas Bates and other contributors - -# Based on Direct Methods for Sparse Linear Systems, T. A. Davis, SIAM, Philadelphia, Sept. 2006. -# Section 2.4: Triplet form -# http://www.cise.ufl.edu/research/sparse/CSparse/ - - -# 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. -# A trivial example is speye(n, n) in which every node is a root. - -""" - etree(A[, post]) - -Compute the elimination tree of a symmetric sparse matrix `A` from `triu(A)` -and, optionally, its post-ordering permutation. -""" -function etree{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, postorder::Bool) - m,n = size(A) - Ap = A.colptr - Ai = A.rowval - parent = zeros(Ti, n) - ancestor = zeros(Ti, n) - for k in 1:n, p in Ap[k]:(Ap[k+1] - 1) - i = Ai[p] - while i != 0 && i < k - inext = ancestor[i] # inext = ancestor of i - ancestor[i] = k # path compression - if (inext == 0) parent[i] = k end # no anc., parent is k - i = inext - end - end - if !postorder return parent end - head = zeros(Ti,n) # empty linked lists - next = zeros(Ti,n) - for j in n:-1:1 # traverse in reverse order - if (parent[j] == 0); continue; end # j is a root - next[j] = head[parent[j]] # add j to list of its parent - head[parent[j]] = j - end - stack = Ti[] - sizehint!(stack, n) - post = zeros(Ti,n) - k = 1 - for j in 1:n - if (parent[j] != 0) continue end # skip j if it is not a root - push!(stack, j) # place j on the stack - while (!isempty(stack)) # while (stack is not empty) - p = stack[end] # p = top of stack - i = head[p] # i = youngest child of p - if (i == 0) - pop!(stack) - post[k] = p # node p is the kth postordered node - k += 1 - else - head[p] = next[i] # remove i from children of p - push!(stack, i) - end - end - end - parent, post -end - -etree(A::SparseMatrixCSC) = etree(A, false) - -# find nonzero pattern of Cholesky L[k,1:k-1] using etree and triu(A[:,k]) -# based on cs_ereach p. 43, "Direct Methods for Sparse Linear Systems" -function ereach{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, k::Integer, parent::Vector{Ti}) - m,n = size(A); Ap = A.colptr; Ai = A.rowval - s = Ti[]; sizehint!(s, n) # to be used as a stack - visited = falses(n) - visited[k] = true - for p in Ap[k]:(Ap[k+1] - 1) - i = Ai[p] # A[i,k] is nonzero - if i > k continue end # only use upper triangular part of A - while !visited[i] # traverse up etree - push!(s,i) # L[k,i] is nonzero - visited[i] = true - i = parent[i] - end - end - s -end - -# based on cs_permute p. 21, "Direct Methods for Sparse Linear Systems" -function csc_permute{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti}, q::Vector{Ti}) - m, n = size(A) - Ap = A.colptr - Ai = A.rowval - Ax = A.nzval - lpinv = length(pinv) - if m != lpinv - throw(DimensionMismatch( - "the number of rows of sparse matrix A must equal the length of pinv, $m != $lpinv")) - end - lq = length(q) - if n != lq - throw(DimensionMismatch( - "the number of columns of sparse matrix A must equal the length of q, $n != $lq")) - end - if !isperm(pinv) || !isperm(q) - throw(ArgumentError("both pinv and q must be permutations")) - end - C = copy(A); Cp = C.colptr; Ci = C.rowval; Cx = C.nzval - nz = one(Ti) - for k in 1:n - Cp[k] = nz - j = q[k] - for t = Ap[j]:(Ap[j+1]-1) - Cx[nz] = Ax[t] - Ci[nz] = pinv[Ai[t]] - nz += one(Ti) - end - end - Cp[n + 1] = nz - (C.').' # double transpose to order the columns -end - - -# based on cs_symperm p. 21, "Direct Methods for Sparse Linear Systems" -# form A[p,p] for a symmetric A stored in the upper triangle -""" - symperm(A, p) - -Return the symmetric permutation of `A`, which is `A[p,p]`. `A` should be -symmetric, sparse, and only contain nonzeros in the upper triangular part of the -matrix is stored. This algorithm ignores the lower triangular part of the -matrix. Only the upper triangular part of the result is returned. -""" -function symperm{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti}) - m, n = size(A) - if m != n - throw(DimensionMismatch("sparse matrix A must be square")) - end - Ap = A.colptr - Ai = A.rowval - Ax = A.nzval - if !isperm(pinv) - throw(ArgumentError("pinv must be a permutation")) - end - lpinv = length(pinv) - if n != lpinv - throw(DimensionMismatch( - "dimensions of sparse matrix A must equal the length of pinv, $((m,n)) != $lpinv")) - end - C = copy(A); Cp = C.colptr; Ci = C.rowval; Cx = C.nzval - w = zeros(Ti,n) - for j in 1:n # count entries in each column of C - j2 = pinv[j] - for p in Ap[j]:(Ap[j+1]-1) - (i = Ai[p]) > j || (w[max(pinv[i],j2)] += one(Ti)) - end - end - Cp[:] = cumsum(vcat(one(Ti),w)) - copy!(w,Cp[1:n]) # needed to be consistent with cs_cumsum - for j in 1:n - j2 = pinv[j] - for p = Ap[j]:(Ap[j+1]-1) - (i = Ai[p]) > j && continue - i2 = pinv[i] - ind = max(i2,j2) - Ci[q = w[ind]] = min(i2,j2) - w[ind] += 1 - Cx[q] = Ax[p] - end - end - (C.').' # double transpose to order the columns -end \ No newline at end of file diff --git a/contrib/add_license_to_files.jl b/contrib/add_license_to_files.jl index 383fff8c059f5..2e4fcb2aee13c 100644 --- a/contrib/add_license_to_files.jl +++ b/contrib/add_license_to_files.jl @@ -34,7 +34,6 @@ const skipfiles = [ # files to check - already copyright # see: https://github.com/JuliaLang/julia/pull/11073#issuecomment-98099389 "../base/special/trig.jl", - "../base/sparse/csparse.jl", "../base/linalg/givens.jl", # "../src/abi_llvm.cpp", diff --git a/doc/stdlib/arrays.rst b/doc/stdlib/arrays.rst index 45f1140a3b6f5..4571aeddc6707 100644 --- a/doc/stdlib/arrays.rst +++ b/doc/stdlib/arrays.rst @@ -958,12 +958,16 @@ dense counterparts. The following functions are specific to sparse arrays. Compute the elimination tree of a symmetric sparse matrix ``A`` from ``triu(A)`` and, optionally, its post-ordering permutation. + Note: This function has been moved to the SuiteSparse.jl package. + .. function:: symperm(A, p) .. Docstring generated from Julia source Return the symmetric permutation of ``A``\ , which is ``A[p,p]``\ . ``A`` should be symmetric, sparse, and only contain nonzeros in the upper triangular part of the matrix is stored. This algorithm ignores the lower triangular part of the matrix. Only the upper triangular part of the result is returned. + Note: This function has been moved to the SuiteSparse.jl package. + .. function:: nonzeros(A) .. Docstring generated from Julia source diff --git a/test/sparsedir/sparse.jl b/test/sparsedir/sparse.jl index 0d31455fbd7ec..600f32708faae 100644 --- a/test/sparsedir/sparse.jl +++ b/test/sparsedir/sparse.jl @@ -343,20 +343,6 @@ end @test full(spdiagm((ones(2), ones(2)), (0, -1), 3, 3)) == [1.0 0.0 0.0; 1.0 1.0 0.0; 0.0 1.0 0.0] -# elimination tree -## upper triangle of the pattern test matrix from Figure 4.2 of -## "Direct Methods for Sparse Linear Systems" by Tim Davis, SIAM, 2006 -rowval = Int32[1,2,2,3,4,5,1,4,6,1,7,2,5,8,6,9,3,4,6,8,10,3,5,7,8,10,11] -colval = Int32[1,2,3,3,4,5,6,6,6,7,7,8,8,8,9,9,10,10,10,10,10,11,11,11,11,11,11] -A = sparse(rowval, colval, ones(length(rowval))) -p = etree(A) -P,post = etree(A, true) -@test P == p -@test P == Int32[6,3,8,6,8,7,9,10,10,11,0] -@test post == Int32[2,3,5,8,1,4,6,7,9,10,11] -@test isperm(post) - - # issue #4986, reinterpret sfe22 = speye(Float64, 2) mfe22 = eye(Float64, 2) @@ -1019,13 +1005,9 @@ A = sparse(tril(rand(5,5))) @test istril(A) @test !istril(sparse(ones(5,5))) -# symperm -srand(1234321) -A = triu(sprand(10,10,0.2)) # symperm operates on upper triangle -perm = randperm(10) -@test symperm(A,perm).colptr == [1,1,2,4,5,5,6,8,8,9,10] - # droptol +srand(1234321) +A = triu(sprand(10,10,0.2)) @test Base.droptol!(A,0.01).colptr == [1,1,1,2,2,3,4,6,6,7,9] @test isequal(Base.droptol!(sparse([1], [1], [1]), 1), SparseMatrixCSC(1,1,Int[1,1],Int[],Int[])) @@ -1224,12 +1206,6 @@ if Base.USE_GPL_LIBS end @test_throws DimensionMismatch Base.SparseArrays.normestinv(sprand(3,5,.9)) -# csc_permute -A = sprand(10,10,0.2) -p = randperm(10) -q = randperm(10) -@test Base.SparseArrays.csc_permute(A, invperm(p), q) == full(A)[p, q] - # issue #13008 @test_throws ArgumentError sparse(collect(1:100), collect(1:100), fill(5,100), 5, 5) @test_throws ArgumentError sparse(Int[], collect(1:5), collect(1:5))