From 3fa0e3ab2e1fc15cb07a3609fed4ef6e1af5616c Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Wed, 27 Apr 2016 11:27:34 +0200 Subject: [PATCH] remove csparse functions The functions and tests for these function have been moved to the package SuiteSparse.jl in the JuliaSparse organization. --- LICENSE.md | 1 - base/sparse.jl | 34 +++++- base/sparse/csparse.jl | 178 -------------------------------- contrib/add_license_to_files.jl | 1 - doc/stdlib/arrays.rst | 6 +- test/sparsedir/sparse.jl | 28 +---- 6 files changed, 40 insertions(+), 208 deletions(-) delete mode 100644 base/sparse/csparse.jl diff --git a/LICENSE.md b/LICENSE.md index 1f00b208885a30..f15371de45da16 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 4a16213bda9711..92d32d28102b86 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,37 @@ if Base.USE_GPL_LIBS include("sparse/spqr.jl") end +# point users to SuiteSparse +const SUITESPARSE_END_ERR_SRING = 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. +""" +etree{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, postorder::Bool) = + error("etree(A[, post])", SUITESPARSE_END_ERR_SRING) +etree(A::SparseMatrixCSC) = etree(A, false) + +ereach{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, k::Integer, parent::Vector{Ti}) = + error("ereach(A, k, parent)", SUITESPARSE_END_ERR_SRING) + +csc_permute{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti}, q::Vector{Ti})= + error("csc_permute(A, pinv, q)", SUITESPARSE_END_ERR_SRING) + +""" + 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. +""" +symperm{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti}) = + error("symperm(A, pinv)", SuiteSparse_END_ERR_SRING) + end diff --git a/base/sparse/csparse.jl b/base/sparse/csparse.jl deleted file mode 100644 index 53a153c8aa77b7..00000000000000 --- 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 53bb003f5267f2..cd751b79839038 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 0a6e60cd8e4396..ab5cd36e1032fb 100644 --- a/doc/stdlib/arrays.rst +++ b/doc/stdlib/arrays.rst @@ -227,7 +227,7 @@ Constructors .. Docstring generated from Julia source - Create an array with the same data as the given array, but with different dimensions. An implementation for a particular type of array may choose whether the data is copied or shared. + Create an array with the same data as the given array, but with different dimensions. .. function:: similar(array, [element_type=eltype(array)], [dims=size(array)]) @@ -971,12 +971,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 CSparse.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 CSparse.jl package. + .. function:: nonzeros(A) .. Docstring generated from Julia source diff --git a/test/sparsedir/sparse.jl b/test/sparsedir/sparse.jl index 9252c212e497de..c2ad9a7343e33b 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))