From c45b61fe5ca04597df7d61522a4f00c13f25ab0f Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 4 Sep 2023 14:44:02 +0200 Subject: [PATCH 1/6] Code quality cleanup (#438) --- src/linalg.jl | 2 +- test/ambiguous.jl | 23 +++++++++++++++++++---- test/linalg.jl | 12 ++++++++++++ 3 files changed, 32 insertions(+), 5 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 8dd9cda3..f315300f 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -560,7 +560,7 @@ function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, tf end return C end -function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, _, xA::AdjOrTrans{<:Any,<:SparseMatrixCSCUnion}, B::AbstractVecOrMat) +function LinearAlgebra.generic_trimatmul!(C::StridedVecOrMat, uploc, isunitc, ::Function, xA::AdjOrTrans{<:Any,<:SparseMatrixCSCUnion}, B::AbstractVecOrMat) A = parent(xA) nrowC = size(C, 1) ncol = checksquare(A) diff --git a/test/ambiguous.jl b/test/ambiguous.jl index e636944d..00c892f5 100644 --- a/test/ambiguous.jl +++ b/test/ambiguous.jl @@ -25,13 +25,13 @@ using Test, LinearAlgebra, SparseArrays, Aqua @testset "code quality" begin @testset "Method ambiguity" begin - # Aqua.test_ambiguities([SparseArrays, Base, Core]) + Aqua.test_ambiguities([SparseArrays, Base, Core]) end @testset "Unbound type parameters" begin @test_broken Aqua.detect_unbound_args_recursively(SparseArrays) == [] end @testset "Undefined exports" begin - Aqua.test_undefined_exports(SparseArrays) == [] + Aqua.test_undefined_exports(SparseArrays) end @testset "Compare Project.toml and test/Project.toml" begin Aqua.test_project_extras(SparseArrays) @@ -50,8 +50,23 @@ using Test, LinearAlgebra, SparseArrays, Aqua end end -@testset "detect_ambiguities" begin - @test_nowarn detect_ambiguities(SparseArrays; recursive=true, ambiguous_bottom=false) +let ambig = detect_ambiguities(SparseArrays; recursive=true) + @test isempty(ambig) + ambig = Set{Any}(((m1.sig, m2.sig) for (m1, m2) in ambig)) + expect = [] + good = true + while !isempty(ambig) + sigs = pop!(ambig) + i = findfirst(==(sigs), expect) + if i === nothing + println(stderr, "push!(expect, (", sigs[1], ", ", sigs[2], "))") + good = false + continue + end + deleteat!(expect, i) + end + @test isempty(expect) + @test good end ## This was the older version that was disabled diff --git a/test/linalg.jl b/test/linalg.jl index c1f5880f..7b8aba54 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -162,6 +162,18 @@ begin vMAW = tr(wr(view(MA, :, 1:n))) @test vAW * B ≈ vMAW * B end + a = sprand(rng, ComplexF64, n, n, 0.01) + ma = Matrix(a) + @testset "triangular multiply with conjugate matrices" for tr in (x -> adjoint(transpose(x)), x -> transpose(adjoint(x))), + wr in (UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular) + AW = tr(wr(a)) + MAW = tr(wr(ma)) + @test AW * B ≈ MAW * B + # and for SparseMatrixCSCView - a view of all rows and unit range of cols + vAW = tr(wr(view(a, :, 1:n))) + vMAW = tr(wr(view(ma, :, 1:n))) + @test vAW * B ≈ vMAW * B + end A = A - Diagonal(diag(A)) + 2I # avoid rounding errors by division MA = Matrix(A) @testset "triangular solver for $tr($wr)" for tr in (identity, adjoint, transpose), From ddfed4c7b763c092f20c6a431398c4bbff12ef9c Mon Sep 17 00:00:00 2001 From: Jameson Nash Date: Tue, 1 Aug 2023 14:34:18 -0400 Subject: [PATCH 2/6] cat: ensure vararg is more inferrable Ensures that Union{} is not part of the output possibilities after type-piracy of Base.cat methods. Refs https://github.com/JuliaLang/julia/issues/50550 (cherry picked from commit c402d09cf05492179fad2def5632e354a81f5b30) --- src/sparsevector.jl | 48 +++++++++++++++++++++++++++------------------ 1 file changed, 29 insertions(+), 19 deletions(-) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 22b1badb..cd75e5b1 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1194,13 +1194,14 @@ anysparse() = false anysparse(X) = X isa AbstractArray && issparse(X) anysparse(X, Xs...) = anysparse(X) || anysparse(Xs...) -function hcat(X::Union{Vector, AbstractSparseVector}...) +const _SparseVecConcatGroup = Union{Vector, AbstractSparseVector} +function hcat(X::_SparseVecConcatGroup...) if anysparse(X...) X = map(sparse, X) end return cat(X...; dims=Val(2)) end -function vcat(X::Union{Vector, AbstractSparseVector}...) +function vcat(X::_SparseVecConcatGroup...) if anysparse(X...) X = map(sparse, X) end @@ -1213,30 +1214,30 @@ end const _SparseConcatGroup = Union{AbstractVecOrMat{<:Number},Number} # `@constprop :aggressive` allows `dims` to be propagated as constant improving return type inference -Base.@constprop :aggressive function Base._cat(dims, X::_SparseConcatGroup...) - T = promote_eltype(X...) - if anysparse(X...) - X = (_sparse(first(X)), map(_makesparse, Base.tail(X))...) +Base.@constprop :aggressive function Base._cat(dims, X1::_SparseConcatGroup, X::_SparseConcatGroup...) + T = promote_eltype(X1, X...) + if anysparse(X1) || anysparse(X...) + X1, X = _sparse(X1), map(_makesparse, X) end - return Base._cat_t(dims, T, X...) + return Base._cat_t(dims, T, X1, X...) end -function hcat(X::_SparseConcatGroup...) - if anysparse(X...) - X = (_sparse(first(X)), map(_makesparse, Base.tail(X))...) +function hcat(X1::_SparseConcatGroup, X::_SparseConcatGroup...) + if anysparse(X1) || anysparse(X...) + X1, X = _sparse(X1), map(_makesparse, X) end - return cat(X..., dims=Val(2)) + return cat(X1, X..., dims=Val(2)) end -function vcat(X::_SparseConcatGroup...) - if anysparse(X...) - X = (_sparse(first(X)), map(_makesparse, Base.tail(X))...) +function vcat(X1::_SparseConcatGroup, X::_SparseConcatGroup...) + if anysparse(X1) || anysparse(X...) + X1, X = _sparse(X1), map(_makesparse, X) end - return cat(X..., dims=Val(1)) + return cat(X1, X..., dims=Val(1)) end -function hvcat(rows::Tuple{Vararg{Int}}, X::_SparseConcatGroup...) - if anysparse(X...) - vcat(_hvcat_rows(rows, X...)...) +function hvcat(rows::Tuple{Vararg{Int}}, X1::_SparseConcatGroup, X::_SparseConcatGroup...) + if anysparse(X1) || anysparse(X...) + vcat(_hvcat_rows(rows, X1, X...)...) else - Base.typed_hvcat(Base.promote_eltypeof(X...), rows, X...) + Base.typed_hvcat(Base.promote_eltypeof(X1, X...), rows, X1, X...) end end function _hvcat_rows((row1, rows...)::Tuple{Vararg{Int}}, X::_SparseConcatGroup...) @@ -1254,6 +1255,15 @@ function _hvcat_rows((row1, rows...)::Tuple{Vararg{Int}}, X::_SparseConcatGroup. end _hvcat_rows(::Tuple{}, X::_SparseConcatGroup...) = () +# disambiguation for type-piracy problems created above +hcat(n1::Number, ns::Vararg{Number}) = invoke(hcat, Tuple{Vararg{Number}}, n1, ns...) +vcat(n1::Number, ns::Vararg{Number}) = invoke(vcat, Tuple{Vararg{Number}}, n1, ns...) +hcat(n1::Type{N}, ns::Vararg{N}) where {N<:Number} = invoke(hcat, Tuple{Vararg{Number}}, n1, ns...) +vcat(n1::Type{N}, ns::Vararg{N}) where {N<:Number} = invoke(vcat, Tuple{Vararg{Number}}, n1, ns...) +hvcat(rows::Tuple{Vararg{Int}}, n1::Number, ns::Vararg{Number}) = invoke(hvcat, Tuple{typeof(rows), Vararg{Number}}, rows, n1, ns...) +hvcat(rows::Tuple{Vararg{Int}}, n1::N, ns::Vararg{N}) where {N<:Number} = invoke(hvcat, Tuple{typeof(rows), Vararg{N}}, rows, n1, ns...) + + # make sure UniformScaling objects are converted to sparse matrices for concatenation promote_to_array_type(A::Tuple{Vararg{Union{_SparseConcatGroup,UniformScaling}}}) = anysparse(A...) ? SparseMatrixCSC : Matrix promote_to_arrays_(n::Int, ::Type{SparseMatrixCSC}, J::UniformScaling) = sparse(J, n, n) From 4a6d98c3f4578ba08c20e7fa8567d32909234035 Mon Sep 17 00:00:00 2001 From: Jameson Nash Date: Tue, 1 Aug 2023 14:36:07 -0400 Subject: [PATCH 3/6] fix inference of SparseVector cat This little gadget creates a closure over Type{T} instead of DataType. Fix https://github.com/JuliaLang/julia/issues/50623 (cherry picked from commit 68afc6e335942052ffa7d4f477d599ddaa02c1f0) --- src/sparsevector.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index cd75e5b1..7e66bc06 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1108,7 +1108,9 @@ function hcat(Xin::AbstractSparseVector...) X = map(_unsafe_unfix, Xin) Tv = promote_type(map(eltype, X)...) Ti = promote_type(map(indtype, X)...) - r = _absspvec_hcat(map(x -> convert(SparseVector{Tv,Ti}, x), X)...) + r = (function (::Type{SV}) where SV + _absspvec_hcat(map(x -> convert(SV, x), X)...) + end)(SparseVector{Tv,Ti}) return @if_move_fixed Xin... r end function _absspvec_hcat(X::AbstractSparseVector{Tv,Ti}...) where {Tv,Ti} @@ -1144,7 +1146,9 @@ function vcat(Xin::AbstractSparseVector...) X = map(_unsafe_unfix, Xin) Tv = promote_type(map(eltype, X)...) Ti = promote_type(map(indtype, X)...) - r = _absspvec_vcat(map(x -> convert(SparseVector{Tv,Ti}, x), X)...) + r = (function (::Type{SV}) where SV + _absspvec_vcat(map(x -> convert(SV, x), X)...) + end)(SparseVector{Tv,Ti}) return @if_move_fixed Xin... r end function _absspvec_vcat(X::AbstractSparseVector{Tv,Ti}...) where {Tv,Ti} From 2a84a1fbfc0bfc2804e7819c1a2bb686511905a7 Mon Sep 17 00:00:00 2001 From: Jameson Nash Date: Fri, 1 Sep 2023 03:36:29 -0400 Subject: [PATCH 4/6] faster cat performance (#432) The generic `cat` does more shape analysis, that typed_cat does not always do (before either may then dispatch to _cat_t), so we can make this faster by calling it instead. Secondly, we can make `issparse` non-recursive once it hits a base case where all trailing elements are the same. This makes it much better at handling large splat, since we do not need to check each recursively smaller type down to the base case using generic code, and just generate one const method specialized on the full length instead. Fix JuliaLang/julia#51011 (cherry picked from commit 4e6776a825f2a26c3c580f9b77ba230a6598d7dd) --- src/sparsevector.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 7e66bc06..9fc8bec8 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -1197,6 +1197,7 @@ _makesparse(x::AbstractMatrix) = convert(SparseMatrixCSC, issparse(x) ? x : spar anysparse() = false anysparse(X) = X isa AbstractArray && issparse(X) anysparse(X, Xs...) = anysparse(X) || anysparse(Xs...) +anysparse(X::T, Xs::T...) where {T} = anysparse(X) const _SparseVecConcatGroup = Union{Vector, AbstractSparseVector} function hcat(X::_SparseVecConcatGroup...) @@ -1229,13 +1230,13 @@ function hcat(X1::_SparseConcatGroup, X::_SparseConcatGroup...) if anysparse(X1) || anysparse(X...) X1, X = _sparse(X1), map(_makesparse, X) end - return cat(X1, X..., dims=Val(2)) + return Base.typed_hcat(Base.promote_eltype(X1, X...), X1, X...) end function vcat(X1::_SparseConcatGroup, X::_SparseConcatGroup...) if anysparse(X1) || anysparse(X...) X1, X = _sparse(X1), map(_makesparse, X) end - return cat(X1, X..., dims=Val(1)) + return Base.typed_vcat(Base.promote_eltype(X1, X...), X1, X...) end function hvcat(rows::Tuple{Vararg{Int}}, X1::_SparseConcatGroup, X::_SparseConcatGroup...) if anysparse(X1) || anysparse(X...) From baa418f9e182b14d24551048d31b0d42e50fcc90 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 24 Aug 2023 17:41:12 +0530 Subject: [PATCH 5/6] Respect `IOContext` while displaying a `SparseMatrixCSC` (#423) * Respect IOContext in printing column indices * Add docstring to ColumnIndices --- src/sparsematrix.jl | 36 +++++++++++++--------- test/sparsematrix_constructors_indexing.jl | 9 ++++++ 2 files changed, 30 insertions(+), 15 deletions(-) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 10bd525a..1703c870 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -344,6 +344,26 @@ function Base.print_array(io::IO, S::AbstractSparseMatrixCSCInclAdjointAndTransp end end +""" + ColumnIndices(S::AbstractSparseMatrixCSC) + +Return the column indices of the stored values in `S`. +This is an internal type that is used in displaying sparse matrices, +and is not a part of the public interface. +""" +struct ColumnIndices{Ti,S<:AbstractSparseMatrixCSC{<:Any,Ti}} <: AbstractVector{Ti} + arr :: S +end + +size(C::ColumnIndices) = (nnz(C.arr),) +# returns the column index of the n-th non-zero value from the column pointer +@inline function getindex(C::ColumnIndices, i::Int) + @boundscheck checkbounds(C, i) + colptr = getcolptr(C.arr) + ind = searchsortedlast(colptr, i) + eltype(C)(ind) +end + # always show matrices as `sparse(I, J, K)` function Base.show(io::IO, _S::AbstractSparseMatrixCSCInclAdjointAndTranspose) _checkbuffers(_S) @@ -358,21 +378,7 @@ function Base.show(io::IO, _S::AbstractSparseMatrixCSCInclAdjointAndTranspose) print(io, "transpose(") end print(io, "sparse(", I, ", ") - if length(I) == 0 - print(io, eltype(getcolptr(S)), "[]") - else - print(io, "[") - il = nnz(S) - 1 - for col in 1:size(S, 2), - k in getcolptr(S)[col] : (getcolptr(S)[col+1]-1) - print(io, col) - if il > 0 - print(io, ", ") - il -= 1 - end - end - print(io, "]") - end + show(io, ColumnIndices(S)) print(io, ", ", K, ", ", m, ", ", n, ")") if _S isa Adjoint || _S isa Transpose print(io, ")") diff --git a/test/sparsematrix_constructors_indexing.jl b/test/sparsematrix_constructors_indexing.jl index 1ef19a1d..5aad5628 100644 --- a/test/sparsematrix_constructors_indexing.jl +++ b/test/sparsematrix_constructors_indexing.jl @@ -1574,6 +1574,15 @@ end _show_with_braille_patterns(ioc, _filled_sparse(8, 16)) @test String(take!(io)) == "⎡⣿⣿⎤\n" * "⎣⣿⣿⎦" + + # respect IOContext while displaying J + I, J, V = shuffle(1:50), shuffle(1:50), [1:50;] + S = sparse(I, J, V) + I, J, V = I[sortperm(J)], sort(J), V[sortperm(J)] + @test repr(S) == "sparse($I, $J, $V, $(size(S,1)), $(size(S,2)))" + limctxt(x) = repr(x, context=:limit=>true) + expstr = "sparse($(limctxt(I)), $(limctxt(J)), $(limctxt(V)), $(size(S,1)), $(size(S,2)))" + @test limctxt(S) == expstr end @testset "issparse for specialized matrix types" begin From a27210928fd7c2e2e5f7c64b69800012c1e7dcfa Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 10 Oct 2023 15:32:11 +0200 Subject: [PATCH 6/6] Adapt to Documenter v1 (#444) * Adapt to Documenter v1 --- docs/Project.toml | 3 +++ docs/make.jl | 3 +-- docs/src/solvers.md | 13 ++++++++----- src/solvers/cholmod.jl | 7 ++++--- 4 files changed, 16 insertions(+), 10 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index dfa65cd1..1814eb33 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,2 +1,5 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +Documenter = "1" diff --git a/docs/make.jl b/docs/make.jl index de3381c5..53843712 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,8 +10,7 @@ makedocs( "SparseArrays" => "index.md", "Sparse Linear Algebra" => "solvers.md", ]; - # strict = true, - strict = Symbol[:doctest], + warnonly = [:missing_docs, :cross_references], ) deploydocs(repo = "github.com/JuliaSparse/SparseArrays.jl.git") diff --git a/docs/src/solvers.md b/docs/src/solvers.md index 91214c5c..e633be9d 100644 --- a/docs/src/solvers.md +++ b/docs/src/solvers.md @@ -14,11 +14,14 @@ Sparse matrix solvers call functions from [SuiteSparse](http://suitesparse.com). Other solvers such as [Pardiso.jl](https://github.com/JuliaSparse/Pardiso.jl/) are available as external packages. [Arpack.jl](https://julialinearalgebra.github.io/Arpack.jl/stable/) provides `eigs` and `svds` for iterative solution of eigensystems and singular value decompositions. -These factorizations are described in more detail in [`Linear Algebra`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/) section of the manual: -1. [`cholesky`](@ref) -2. [`ldlt`](@ref) -3. [`lu`](@ref) -4. [`qr`](@ref) +These factorizations are described in more detail in the +[`Linear Algebra`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/) +section of the manual: + +1. [`cholesky`](@ref SparseArrays.CHOLMOD.cholesky) +2. [`ldlt`](@ref SparseArrays.CHOLMOD.ldlt) +3. [`lu`](@ref SparseArrays.UMFPACK.lu) +4. [`qr`](@ref SparseArrays.SPQR.qr) ```@docs SparseArrays.CHOLMOD.cholesky diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 9c832ef7..731eaa3f 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -18,7 +18,8 @@ using LinearAlgebra using LinearAlgebra: RealHermSymComplexHerm, AdjOrTrans import LinearAlgebra: (\), AdjointFactorization, cholesky, cholesky!, det, diag, ishermitian, isposdef, - issuccess, issymmetric, ldlt, ldlt!, logdet, lowrankdowndate! + issuccess, issymmetric, ldlt, ldlt!, logdet, + lowrankdowndate, lowrankdowndate!, lowrankupdate, lowrankupdate! using SparseArrays using SparseArrays: getcolptr, AbstractSparseVecOrMat @@ -1549,7 +1550,7 @@ factor will be `L*L' == P*A*P' + C'*C` `update`: `Cint(1)` for `A + CC'`, `Cint(0)` for `A - CC'` """ -lowrankdowndate! +lowrankupdowndate! #Helper functions for rank updates lowrank_reorder(V::AbstractArray,p) = Sparse(sparse(V[p,:])) @@ -1598,7 +1599,7 @@ lowrankupdate(F::Factor{Tv}, V::AbstractArray{Tv}) where {Tv<:VTypes} = lowrankupdate!(copy(F), V) """ - lowrankupdate(F::CHOLMOD.Factor, C::AbstractArray) -> FF::CHOLMOD.Factor + lowrankdowndate(F::CHOLMOD.Factor, C::AbstractArray) -> FF::CHOLMOD.Factor Get an `LDLt` Factorization of `A + C*C'` given an `LDLt` or `LLt` factorization `F` of `A`.