diff --git a/docs/make.jl b/docs/make.jl index 53843712..df9e506c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -8,7 +8,6 @@ makedocs( sitename = "SparseArrays", pages = Any[ "SparseArrays" => "index.md", - "Sparse Linear Algebra" => "solvers.md", ]; warnonly = [:missing_docs, :cross_references], ) diff --git a/docs/src/index.md b/docs/src/index.md index 3868d041..dcfc729f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -206,6 +206,27 @@ section of the standard library reference. | [`sprandn(m,n,d)`](@ref) | [`randn(m,n)`](@ref) | Creates a *m*-by-*n* random matrix (of density *d*) with iid non-zero elements distributed according to the standard normal (Gaussian) distribution. | | [`sprandn(rng,m,n,d)`](@ref) | [`randn(rng,m,n)`](@ref) | Creates a *m*-by-*n* random matrix (of density *d*) with iid non-zero elements generated with the `rng` random number generator | +## [Sparse Linear Algebra](@id stdlib-sparse-linalg) + +Sparse matrix solvers call functions from [SuiteSparse](http://suitesparse.com). The following factorizations are available: + +| Type | Description | +|:----------------------|:--------------------------------------------- | +| `CHOLMOD.Factor` | Cholesky and LDLt factorizations | +| `UMFPACK.UmfpackLU` | LU factorization | +| `SPQR.QRSparse` | QR factorization | + +These factorizations are described in more detail in the [Sparse Linear Algebra API section](@ref stdlib-sparse-linalg-api): + +1. [`cholesky`](@ref SparseArrays.CHOLMOD.cholesky) +2. [`ldlt`](@ref SparseArrays.CHOLMOD.ldlt) +3. [`lu`](@ref SparseArrays.UMFPACK.lu) +4. [`qr`](@ref SparseArrays.SPQR.qr) + +```@meta +DocTestSetup = nothing +``` + # [SparseArrays API](@id stdlib-sparse-arrays) ```@docs @@ -245,6 +266,26 @@ SparseArrays.ftranspose! ```@meta DocTestSetup = nothing ``` + +# [Sparse Linear Algebra API](@id stdlib-sparse-linalg-api) + +```@docs +SparseArrays.CHOLMOD.cholesky +SparseArrays.CHOLMOD.cholesky! +SparseArrays.CHOLMOD.lowrankupdate +SparseArrays.CHOLMOD.lowrankupdate! +SparseArrays.CHOLMOD.lowrankdowndate +SparseArrays.CHOLMOD.lowrankdowndate! +SparseArrays.CHOLMOD.lowrankupdowndate! +SparseArrays.CHOLMOD.ldlt +SparseArrays.UMFPACK.lu +SparseArrays.SPQR.qr +``` + +```@meta +DocTestSetup = nothing +``` + # Noteworthy External Sparse Packages Several other Julia packages provide sparse matrix implementations that should be mentioned: @@ -264,3 +305,15 @@ Several other Julia packages provide sparse matrix implementations that should b 7. [ExtendableSparse.jl](https://github.com/j-fu/ExtendableSparse.jl) enables fast insertion into sparse matrices using a lazy approach to new stored indices. 8. [Finch.jl](https://github.com/willow-ahrens/Finch.jl) supports extensive multidimensional sparse array formats and operations through a mini tensor language and compiler, all in native Julia. Support for COO, CSF, CSR, CSC and more, as well as operations like broadcast, reduce, etc. and custom operations. + +External packages providing sparse direct solvers: +1. [KLU.jl](https://github.com/JuliaSparse/KLU.jl) +2. [Pardiso.jl](https://github.com/JuliaSparse/Pardiso.jl/) + +External packages providing solvers for iterative solution of eigensystems and singular value decompositions: +1. [ArnoldiMethods.jl](https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl) +2. [KrylovKit](https://github.com/Jutho/KrylovKit.jl) +3. [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl) + +External packages for working with graphs: +1. [Graphs.jl](https://github.com/JuliaGraphs/Graphs.jl) diff --git a/docs/src/solvers.md b/docs/src/solvers.md deleted file mode 100644 index e633be9d..00000000 --- a/docs/src/solvers.md +++ /dev/null @@ -1,44 +0,0 @@ -# Sparse Linear Algebra - -```@meta -DocTestSetup = :(using LinearAlgebra, SparseArrays) -``` - -Sparse matrix solvers call functions from [SuiteSparse](http://suitesparse.com). The following factorizations are available: - -| Type | Description | -|:--------------------------------- |:--------------------------------------------- | -| `CHOLMOD.Factor` | Cholesky factorization | -| `UMFPACK.UmfpackLU` | LU factorization | -| `SPQR.QRSparse` | QR factorization | - -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 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 -SparseArrays.CHOLMOD.cholesky! -SparseArrays.CHOLMOD.ldlt -SparseArrays.SPQR.qr -SparseArrays.UMFPACK.lu -``` - -```@docs -SparseArrays.CHOLMOD.lowrankupdate -SparseArrays.CHOLMOD.lowrankupdate! -SparseArrays.CHOLMOD.lowrankdowndate -SparseArrays.CHOLMOD.lowrankdowndate! -SparseArrays.CHOLMOD.lowrankupdowndate! -``` - -```@meta -DocTestSetup = nothing -``` diff --git a/src/solvers/cholmod.jl b/src/solvers/cholmod.jl index 117d801b..e155e72f 100644 --- a/src/solvers/cholmod.jl +++ b/src/solvers/cholmod.jl @@ -795,7 +795,7 @@ function ssmult(A::Sparse{Tv1, Ti1}, B::Sparse{Tv2, Ti2}, stype::Integer, A, B = convert.(Sparse{promote_type(Tv1, Tv2), promote_type(Ti1, Ti2)}, (A, B)) return ssmult(A, B, stype, values, sorted) end -function horzcat(A::Sparse{Tv1, Ti1}, B::Sparse{Tv2, Ti2}, values::Bool) where +function horzcat(A::Sparse{Tv1, Ti1}, B::Sparse{Tv2, Ti2}, values::Bool) where {Tv1<:VRealTypes, Tv2<:VRealTypes, Ti1, Ti2} A, B = convert.(Sparse{promote_type(Tv1, Tv2), promote_type(Ti1, Ti2)}, (A, B)) return horzcat(A, B, values) @@ -809,7 +809,7 @@ function sdmult!(A::Sparse{Tv1, Ti}, transpose::Bool, A, X = convert(Sparse{Tv3, Ti}, A), convert(Dense{Tv3}, X) return sdmult!(A, transpose, α, β, X, Y) end -function vertcat(A::Sparse{Tv1, Ti1}, B::Sparse{Tv2, Ti2}, values::Bool) where +function vertcat(A::Sparse{Tv1, Ti1}, B::Sparse{Tv2, Ti2}, values::Bool) where {Tv1<:VRealTypes, Ti1, Tv2<:VRealTypes, Ti2} A, B = convert.(Sparse{promote_type(Tv1, Tv2), promote_type(Ti1, Ti2)}, (A, B)) return vertcat(A, B, values) @@ -895,7 +895,7 @@ function Dense(A::StridedVecOrMatInclAdjAndTrans) return Dense{T}(A) end # Don't always promote to Float64 now that we have Float32 support. -Dense(A::StridedVecOrMatInclAdjAndTrans{T}) where +Dense(A::StridedVecOrMatInclAdjAndTrans{T}) where {T<:Union{Float16, ComplexF16, Float32, ComplexF32}} = Dense{promote_type(T, Float32)}(A) @@ -1055,8 +1055,8 @@ Sparse(A::Hermitian{Tv, SparseMatrixCSC{Tv,Ti}}) where {Tv, Ti} = Sparse{promote_type(Tv, Float64), Ti <: ITypes ? Ti : promote_type(Ti, Int)}( A.data, A.uplo == 'L' ? -1 : 1 ) -Sparse(A::Hermitian{Tv, SparseMatrixCSC{Tv,Ti}}) where - {Tv<:Union{Float16, Float32, ComplexF32, ComplexF16}, Ti} = +Sparse(A::Hermitian{Tv, SparseMatrixCSC{Tv,Ti}}) where + {Tv<:Union{Float16, Float32, ComplexF32, ComplexF16}, Ti} = Sparse{promote_type(Float32, Tv), Ti <: ITypes ? Ti : promote_type(Ti, Int)}( A.data, A.uplo == 'L' ? -1 : 1 ) @@ -1076,7 +1076,7 @@ function Base.convert(::Type{Sparse{Tnew, Inew}}, A::Sparse{Tv, Ti}) where {Tnew a = unsafe_load(typedpointer(A)) S = allocate_sparse(a.nrow, a.ncol, a.nzmax, Bool(a.sorted), Bool(a.packed), a.stype, Tnew, Inew) s = unsafe_load(typedpointer(S)) - + ap = unsafe_wrap(Array, a.p, (a.ncol + 1,), own = false) sp = unsafe_wrap(Array, s.p, (s.ncol + 1,), own = false) copyto!(sp, ap) @@ -1376,7 +1376,7 @@ end ## Multiplication (*)(A::Sparse, B::Sparse) = ssmult(A, B, 0, true, true) -(*)(A::Sparse, B::Dense) = sdmult!(A, false, 1., 0., B, +(*)(A::Sparse, B::Dense) = sdmult!(A, false, 1., 0., B, zeros(size(A, 1), size(B, 2), promote_type(eltype(A), eltype(B))) ) (*)(A::Sparse, B::VecOrMat) = (*)(A, Dense(B)) @@ -1413,7 +1413,7 @@ function *(adjA::Adjoint{<:Any,<:Sparse}, B::Sparse) end *(adjA::Adjoint{<:Any,<:Sparse}, B::Dense) = ( - A = parent(adjA); sdmult!(A, true, 1., 0., B, + A = parent(adjA); sdmult!(A, true, 1., 0., B, zeros(size(A, 2), size(B, 2), promote_type(eltype(A), eltype(B)))) ) *(adjA::Adjoint{<:Any,<:Sparse}, B::VecOrMat) = adjA * Dense(B) @@ -1423,25 +1423,33 @@ end ## Compute that symbolic factorization only function symbolic(A::Sparse{<:VTypes, Ti}; - perm::Union{Nothing,AbstractVector{<:Integer}}=nothing, - postorder::Bool=isnothing(perm)||isempty(perm), userperm_only::Bool=true) where Ti + perm::Union{Nothing,AbstractVector{<:Integer}}=nothing, + postorder::Bool=isnothing(perm)||isempty(perm), + userperm_only::Bool=true, + nested_dissection::Bool=false) where Ti sA = unsafe_load(pointer(A)) sA.stype == 0 && throw(ArgumentError("sparse matrix is not symmetric/Hermitian")) - @cholmod_param postorder = postorder begin - if perm === nothing || isempty(perm) # TODO: deprecate empty perm - return analyze(A) - else # user permutation provided - if userperm_only # use perm even if it is worse than AMD - @cholmod_param nmethods = 1 begin + # The default is to just use AMD. Use nested dissection only if explicitly asked for. + # https://github.com/JuliaSparse/SparseArrays.jl/issues/548 + # https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/26ababc7f3af725c5fb9168a1b94850eab74b666/CHOLMOD/Include/cholmod.h#L555-L574 + @cholmod_param nmethods = (nested_dissection ? 0 : 2) begin + @cholmod_param postorder = postorder begin + if perm === nothing || isempty(perm) # TODO: deprecate empty perm + return analyze(A) + else # user permutation provided + if userperm_only # use perm even if it is worse than AMD + @cholmod_param nmethods = 1 begin + return analyze_p(A, Ti[p-1 for p in perm]) + end + else return analyze_p(A, Ti[p-1 for p in perm]) end - else - return analyze_p(A, Ti[p-1 for p in perm]) end end end + end function cholesky!(F::Factor{Tv}, A::Sparse{Tv}; @@ -1467,7 +1475,7 @@ See also [`cholesky`](@ref). !!! note This method uses the CHOLMOD library from SuiteSparse, which only supports - real or complex types in single or double precision. + real or complex types in single or double precision. Input matrices not of those element types will be converted to these types as appropriate. """ @@ -1587,8 +1595,8 @@ true !!! note This method uses the CHOLMOD[^ACM887][^DavisHager2009] library from [SuiteSparse](https://github.com/DrTimothyAldenDavis/SuiteSparse). - CHOLMOD only supports real or complex types in single or double precision. - Input matrices not of those element types will be + CHOLMOD only supports real or complex types in single or double precision. + Input matrices not of those element types will be converted to these types as appropriate. Many other functions from CHOLMOD are wrapped but not exported from the @@ -1633,8 +1641,8 @@ have the type tag, it must still be symmetric or Hermitian. See also [`ldlt`](@ref). !!! note - This method uses the CHOLMOD library from [SuiteSparse](https://github.com/DrTimothyAldenDavis/SuiteSparse), - which only supports real or complex types in single or double precision. + This method uses the CHOLMOD library from [SuiteSparse](https://github.com/DrTimothyAldenDavis/SuiteSparse), + which only supports real or complex types in single or double precision. Input matrices not of those element types will be converted to these types as appropriate. """ @@ -1695,7 +1703,7 @@ it should be a permutation of `1:size(A,1)` giving the ordering to use !!! note This method uses the CHOLMOD[^ACM887][^DavisHager2009] library from [SuiteSparse](https://github.com/DrTimothyAldenDavis/SuiteSparse). - CHOLMOD only supports real or complex types in single or double precision. + CHOLMOD only supports real or complex types in single or double precision. Input matrices not of those element types will be converted to these types as appropriate. @@ -1767,7 +1775,7 @@ See also [`lowrankupdate!`](@ref), [`lowrankdowndate`](@ref), [`lowrankdowndate! """ lowrankupdate(F::Factor{Tv}, V::AbstractArray{Tv2}) where {Tv, Tv2} = lowrankupdate!( - change_xdtype(F, promote_type(Tv, Tv2)), + change_xdtype(F, promote_type(Tv, Tv2)), convert(AbstractArray{promote_type(Tv, Tv2)}, V) ) @@ -1782,7 +1790,7 @@ See also [`lowrankdowndate!`](@ref), [`lowrankupdate`](@ref), [`lowrankupdate!`] """ lowrankdowndate(F::Factor{Tv}, V::AbstractArray{Tv2}) where {Tv, Tv2} = lowrankdowndate!( - change_xdtype(F, promote_type(Tv, Tv2)), + change_xdtype(F, promote_type(Tv, Tv2)), convert(AbstractArray{promote_type(Tv, Tv2)}, V) )