Merge branch 'main' into do/csc-interface-baseline

Author: termi-official

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],
     )

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)

docs/src/solvers.md

@@ -18,7 +18,7 @@ using LinearAlgebra
 using LinearAlgebra: RealHermSymComplexHerm, AdjOrTrans
 import LinearAlgebra: (\), AdjointFactorization,
                  cholesky, cholesky!, det, diag, ishermitian, isposdef,
-                 issuccess, issymmetric, ldlt, ldlt!, logdet,
+                 issuccess, issymmetric, ldiv!, ldlt, ldlt!, logdet,
                  lowrankdowndate, lowrankdowndate!, lowrankupdate, lowrankupdate!
 
 using SparseArrays
@@ -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)
 
 
@@ -913,6 +913,32 @@ function Base.convert(::Type{Dense{Tnew}}, A::Dense{T}) where {Tnew, T}
 end
 Base.convert(::Type{Dense{T}}, A::Dense{T}) where T = A
 
+# Just calling Dense(x) or Dense(b) will allocate new
+# `cholmod_dense_struct`s in CHOLMOD. Instead, we want to reuse
+# the existing memory. We can do this by creating new
+# `cholmod_dense_struct`s and filling them manually.
+function wrap_dense_and_ptr(x::StridedVecOrMat{T}) where {T <: VTypes}
+    dense_x = cholmod_dense_struct()
+    dense_x.nrow = size(x, 1)
+    dense_x.ncol = size(x, 2)
+    dense_x.nzmax = length(x)
+    dense_x.d = stride(x, 2)
+    dense_x.x = pointer(x)
+    dense_x.z = C_NULL
+    dense_x.xtype = xtyp(eltype(x))
+    dense_x.dtype = dtyp(eltype(x))
+    return dense_x, pointer_from_objref(dense_x)
+end
+# We need to use a special handling for the case of `Dense`
+# input arrays since the `pointer` refers to the pointer to the
+# `cholmod_dense`, not to the array values themselves as for
+# standard arrays.
+function wrap_dense_and_ptr(x::Dense{T}) where {T <: VTypes}
+    dense_x_ptr = x.ptr
+    dense_x = unsafe_load(dense_x_ptr)
+    return dense_x, pointer_from_objref(dense_x)
+end
+
 # This constructor assumes zero based colptr and rowval
 function Sparse(m::Integer, n::Integer,
         colptr0::Vector{Ti}, rowval0::Vector{Ti},
@@ -1055,8 +1081,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 +1102,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 +1402,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 +1439,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 +1449,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 +1501,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 +1621,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 +1667,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 +1729,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 +1801,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 +1816,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)
 )
 
@@ -1905,6 +1939,66 @@ const AbstractSparseVecOrMatInclAdjAndTrans = Union{AbstractSparseVecOrMat, AdjO
     throw(ArgumentError("self-adjoint sparse system solve not implemented for sparse rhs B," *
         " consider to convert B to a dense array"))
 
+# in-place ldiv!
+for TI in IndexTypes
+    @eval function ldiv!(x::StridedVecOrMat{T},
+                         L::Factor{T, $TI},
+                         b::StridedVecOrMat{T}) where {T<:VTypes}
+        if x === b
+            throw(ArgumentError("output array must not be aliased with input array"))
+        end
+        if size(L, 1) != size(b, 1)
+            throw(DimensionMismatch("Factorization and RHS should have the same number of rows. " *
+                "Factorization has $(size(L, 2)) rows, but RHS has $(size(b, 1)) rows."))
+        end
+        if size(L, 2) != size(x, 1)
+            throw(DimensionMismatch("Factorization and solution should match sizes. " *
+                "Factorization has $(size(L, 1)) columns, but solution has $(size(x, 1)) rows."))
+        end
+        if size(x, 2) != size(b, 2)
+            throw(DimensionMismatch("Solution and RHS should have the same number of columns. " *
+                "Solution has $(size(x, 2)) columns, but RHS has $(size(b, 2)) columns."))
+        end
+        if !issuccess(L)
+            s = unsafe_load(pointer(L))
+            if s.is_ll == 1
+                throw(LinearAlgebra.PosDefException(s.minor))
+            else
+                throw(LinearAlgebra.ZeroPivotException(s.minor))
+            end
+        end
+
+        # Just calling Dense(x) or Dense(b) will allocate new
+        # `cholmod_dense_struct`s in CHOLMOD. Instead, we want to reuse
+        # the existing memory. We can do this by creating new
+        # `cholmod_dense_struct`s and filling them manually.
+        dense_x, dense_x_ptr = wrap_dense_and_ptr(x)
+        dense_b, dense_b_ptr = wrap_dense_and_ptr(b)
+
+        X_Handle = Ptr{cholmod_dense_struct}(dense_x_ptr)
+        Y_Handle = Ptr{cholmod_dense_struct}(C_NULL)
+        E_Handle = Ptr{cholmod_dense_struct}(C_NULL)
+        status = GC.@preserve x dense_x b dense_b begin
+            $(cholname(:solve2, TI))(
+                CHOLMOD_A, L,
+                Ref(dense_b), C_NULL,
+                Ref(X_Handle), C_NULL,
+                Ref(Y_Handle),
+                Ref(E_Handle),
+                getcommon($TI))
+        end
+        if Y_Handle != C_NULL
+            free!(Y_Handle)
+        end
+        if E_Handle != C_NULL
+            free!(E_Handle)
+        end
+        @assert !iszero(status)
+
+        return x
+    end
+end
+
 ## Other convenience methods
 function diag(F::Factor{Tv, Ti}) where {Tv, Ti}
     f = unsafe_load(typedpointer(F))