finishing touches

Author: ChrisRackauckas

Project.toml

@@ -25,7 +25,6 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
@@ -44,6 +43,7 @@ Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
 Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
+Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
 
 [extensions]
 LinearSolveBandedMatricesExt = "BandedMatrices"
@@ -60,6 +60,8 @@ LinearSolveMetalExt = "Metal"
 LinearSolvePardisoExt = "Pardiso"
 LinearSolveRecursiveArrayToolsExt = "RecursiveArrayTools"
 LinearSolveRecursiveFactorizationExt = "RecursiveFactorization"
+LinearSolveSparseArraysExt = "SparseArrays"
+LinearSolveSparspakExt = "Sparspak"
 
 [compat]
 AllocCheck = "0.2"
@@ -141,10 +143,11 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "Test", "IterativeSolvers", "InteractiveUtils", "JET", "KrylovKit", "KrylovPreconditioners", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff", "HYPRE", "MPI", "BlockDiagonals", "Enzyme", "FiniteDiff", "BandedMatrices", "FastAlmostBandedMatrices", "StaticArrays", "AllocCheck", "StableRNGs", "Zygote", "RecursiveFactorization"]
+test = ["Aqua", "Test", "IterativeSolvers", "InteractiveUtils", "JET", "KrylovKit", "KrylovPreconditioners", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff", "HYPRE", "MPI", "BlockDiagonals", "Enzyme", "FiniteDiff", "BandedMatrices", "FastAlmostBandedMatrices", "StaticArrays", "AllocCheck", "StableRNGs", "Zygote", "RecursiveFactorization", "Sparspak"]

ext/LinearSolveSparseArraysExt.jl

@@ -1,4 +1,4 @@
-module LinearSolveSparseArrays
+module LinearSolveSparseArraysExt
 
 using LinearSolve, LinearAlgebra
 using SparseArrays
@@ -11,12 +11,44 @@ using SparseArrays: AbstractSparseMatrixCSC, nonzeros, rowvals, getcolptr
 include("../src/KLU/klu.jl")
 
 LinearSolve.issparsematrixcsc(A::AbstractSparseMatrixCSC) = true
+LinearSolve.issparsematrix(A::AbstractSparseArray) = true
+LinearSolve.make_SparseMatrixCSC(A::AbstractSparseArray) = SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A))
+LinearSolve.makeempty_SparaseMatrixCSC(A::AbstractSparseArray) = SparseMatrixCSC(0, 0, [1], Int[], eltype(A)[])
+
+function LinearSolve.init_cacheval(alg::RFLUFactorization,
+        A::Union{AbstractSparseArray, LinearSolve.SciMLOperators.AbstractSciMLOperator}, b, u, Pl, Pr,
+        maxiters::Int,
+        abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions)
+    nothing, nothing
+end
+
+
+function LinearSolve.init_cacheval(
+        alg::QRFactorization, A::Symmetric{<:Number, <:SparseMatrixCSC}, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol, verbose::Bool,
+        assumptions::OperatorAssumptions)
+    return nothing
+end
 
 function LinearSolve.handle_sparsematrixcsc_lu(A::AbstractSparseMatrixCSC)
     lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)),
     check = false)
 end
 
+function LinearSolve.defaultalg(
+        A::Symmetric{<:Number, <:SparseMatrixCSC}, b, ::OperatorAssumptions{Bool})
+        LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.CHOLMODFactorization)
+end
+
+function LinearSolve.defaultalg(A::AbstractSparseMatrixCSC{Tv, Ti}, b,
+        assump::OperatorAssumptions{Bool}) where {Tv, Ti}
+    if assump.issq
+        DefaultLinearSolver(DefaultAlgorithmChoice.SparspakFactorization)
+    else
+        error("Generic number sparse factorization for non-square is not currently handled")
+    end
+end
+
 function LinearSolve.init_cacheval(alg::GenericFactorization,
         A::Union{Hermitian{T, <:SparseMatrixCSC},
             Symmetric{T, <:SparseMatrixCSC}}, b, u, Pl, Pr,
@@ -46,7 +78,7 @@ function LinearSolve.init_cacheval(alg::UMFPACKFactorization, A::AbstractSparseA
         rowvals(A), nonzeros(A)))
 end
 
-function SciMLBase.solve!(cache::LinearCache, alg::UMFPACKFactorization; kwargs...)
+function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::UMFPACKFactorization; kwargs...)
     A = cache.A
     A = convert(AbstractMatrix, A)
     if cache.isfresh
@@ -101,7 +133,7 @@ function LinearSolve.init_cacheval(alg::KLUFactorization, A::AbstractSparseArray
 end
 
 # TODO: guard this against errors
-function SciMLBase.solve!(cache::LinearCache, alg::KLUFactorization; kwargs...)
+function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::KLUFactorization; kwargs...)
     A = cache.A
     A = convert(AbstractMatrix, A)
     if cache.isfresh
@@ -146,11 +178,11 @@ function LinearSolve.init_cacheval(alg::CHOLMODFactorization,
 end
 
 function LinearSolve.init_cacheval(alg::NormalCholeskyFactorization,
-        A::Union{AbstractSparseArray, GPUArraysCore.AnyGPUArray,
+        A::Union{AbstractSparseArray, LinearSolve.GPUArraysCore.AnyGPUArray,
             Symmetric{<:Number, <:AbstractSparseArray}}, b, u, Pl, Pr,
         maxiters::Int, abstol, reltol, verbose::Bool,
         assumptions::OperatorAssumptions)
-    ArrayInterface.cholesky_instance(convert(AbstractMatrix, A))
+        LinearSolve.ArrayInterface.cholesky_instance(convert(AbstractMatrix, A))
 end
 
 # Specialize QR for the non-square case
@@ -170,16 +202,16 @@ function LinearSolve._ldiv!(x::AbstractVector,
 end
 
 # Ambiguity removal
-function LinearSolve._ldiv!(::SVector,
+function LinearSolve._ldiv!(::LinearSolve.SVector,
         A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR,
             LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse},
         b::AbstractVector)
     (A \ b)
 end
-function LinearSolve._ldiv!(::SVector,
+function LinearSolve._ldiv!(::LinearSolve.SVector,
         A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR,
             LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse},
-        b::SVector)
+        b::LinearSolve.SVector)
     (A \ b)
 end
 

ext/LinearSolveSparsepakExt.jl renamed to ext/LinearSolveSparspakExt.jl

@@ -1,9 +1,9 @@
-module LinearSolveSparsepakExt
+module LinearSolveSparspakExt
 
 using LinearSolve, LinearAlgebra
-using SparseArrays
-using SparseArrays: AbstractSparseMatrixCSC, nonzeros, rowvals, getcolptr
 using Sparspak
+using Sparspak.SparseCSCInterface.SparseArrays
+using SparseArrays: AbstractSparseMatrixCSC, nonzeros, rowvals, getcolptr
 
 const PREALLOCATED_SPARSEPAK = sparspaklu(SparseMatrixCSC(0, 0, [1], Int[], Float64[]),
     factorize = false)
@@ -15,7 +15,7 @@ function LinearSolve.init_cacheval(::SparspakFactorization, A::SparseMatrixCSC{F
     PREALLOCATED_SPARSEPAK
 end
 
-function init_cacheval(::SparspakFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol,
+function LinearSolve.init_cacheval(::SparspakFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol,
         reltol,
         verbose::Bool, assumptions::OperatorAssumptions)
     A = convert(AbstractMatrix, A)
@@ -30,7 +30,7 @@ function init_cacheval(::SparspakFactorization, A, b, u, Pl, Pr, maxiters::Int,
     end
 end
 
-function SciMLBase.solve!(cache::LinearCache, alg::SparspakFactorization; kwargs...)
+function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::SparspakFactorization; kwargs...)
     A = cache.A
     if cache.isfresh
         if cache.cacheval !== nothing && alg.reuse_symbolic
@@ -48,4 +48,12 @@ function SciMLBase.solve!(cache::LinearCache, alg::SparspakFactorization; kwargs
     SciMLBase.build_linear_solution(alg, y, nothing, cache)
 end
 
+LinearSolve.PrecompileTools.@compile_workload begin
+    A = sprand(4, 4, 0.3) + I
+    b = rand(4)
+    prob = LinearProblem(A * A', b)
+    sol = solve(prob) # in case sparspak is used as default
+    sol = solve(prob, SparspakFactorization())
+end
+
 end

src/LinearSolve.jl

@@ -91,6 +91,9 @@ function _fast_sym_givens! end
 
 issparsematrixcsc(A) = false
 handle_sparsematrixcsc_lu(A) = lu(A)
+issparsematrix(A) = false 
+make_SparseMatrixCSC(A) = nothing
+makeempty_SparaseMatrixCSC(A) = nothing
 
 EnumX.@enumx DefaultAlgorithmChoice begin
     LUFactorization
@@ -207,14 +210,6 @@ PrecompileTools.@compile_workload begin
     sol = solve(prob, KrylovJL_GMRES())
 end
 
-PrecompileTools.@compile_workload begin
-    A = sprand(4, 4, 0.3) + I
-    b = rand(4)
-    prob = LinearProblem(A * A', b)
-    sol = solve(prob) # in case sparspak is used as default
-    sol = solve(prob, SparspakFactorization())
-end
-
 ALREADY_WARNED_CUDSS = Ref{Bool}(false)
 error_no_cudss_lu(A) = nothing
 cudss_loaded(A) = false

src/common.jl

@@ -194,19 +194,19 @@ function SciMLBase.init(prob::LinearProblem, alg::SciMLLinearSolveAlgorithm,
         A
     elseif A isa Array
         copy(A)
-    elseif A isa AbstractSparseMatrixCSC
-        SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A))
+    elseif issparsematrixcsc(A)
+        make_SparseMatrixCSC(A)
     else
         deepcopy(A)
     end
 
-    b = if b isa SparseArrays.AbstractSparseArray && !(A isa Diagonal)
+    b = if issparsematrix(b) && !(A isa Diagonal)
         Array(b) # the solution to a linear solve will always be dense!
     elseif alias_b || b isa SVector
         b
     elseif b isa Array
         copy(b)
-    elseif b isa AbstractSparseMatrixCSC
+    elseif issparsematrixcsc(b)
         SparseMatrixCSC(size(b)..., getcolptr(b), rowvals(b), nonzeros(b))
     else
         deepcopy(b)

src/default.jl

@@ -92,20 +92,6 @@ function defaultalg(A::Symmetric{<:Number, <:Array}, b, ::OperatorAssumptions{Bo
     DefaultLinearSolver(DefaultAlgorithmChoice.BunchKaufmanFactorization)
 end
 
-function defaultalg(
-        A::Symmetric{<:Number, <:SparseMatrixCSC}, b, ::OperatorAssumptions{Bool})
-    DefaultLinearSolver(DefaultAlgorithmChoice.CHOLMODFactorization)
-end
-
-function defaultalg(A::AbstractSparseMatrixCSC{Tv, Ti}, b,
-        assump::OperatorAssumptions{Bool}) where {Tv, Ti}
-    if assump.issq
-        DefaultLinearSolver(DefaultAlgorithmChoice.SparspakFactorization)
-    else
-        error("Generic number sparse factorization for non-square is not currently handled")
-    end
-end
-
 function defaultalg(A::GPUArraysCore.AnyGPUArray, b, assump::OperatorAssumptions{Bool})
     if assump.condition === OperatorCondition.IllConditioned || !assump.issq
         DefaultLinearSolver(DefaultAlgorithmChoice.QRFactorization)
@@ -308,7 +294,7 @@ cache.cacheval = NamedTuple(LUFactorization = cache of LUFactorization, ...)
     caches = map(first.(EnumX.symbol_map(DefaultAlgorithmChoice.T))) do alg
         if alg === :KrylovJL_GMRES || alg === :KrylovJL_CRAIGMR || alg === :KrylovJL_LSMR
             quote
-                if A isa Matrix || A isa SparseMatrixCSC
+                if A isa Matrix || issparsematrixcsc(A)
                     nothing
                 else
                     init_cacheval($(algchoice_to_alg(alg)), A, b, u, Pl, Pr, maxiters,

src/extension_algs.jl

@@ -99,11 +99,12 @@ struct RFLUFactorization{P, T} <: AbstractDenseFactorization
             throwerror &&
                 error("RFLUFactorization requires that RecursiveFactorization.jl is loaded, i.e. `using RecursiveFactorization`")
         end
+        new{P, T}()
     end
 end
 
-function RFLUFactorization(; pivot = Val(true), thread = Val(true))
-    RFLUFactorization(pivot, thread)
+function RFLUFactorization(; pivot = Val(true), thread = Val(true), throwerror = true)
+    RFLUFactorization(pivot, thread; throwerror)
 end
 
 """

src/factorization.jl

@@ -45,13 +45,6 @@ function LinearSolve.init_cacheval(
     PREALLOCATED_LU, ipiv
 end
 
-function LinearSolve.init_cacheval(alg::RFLUFactorization,
-        A::Union{AbstractSparseArray, AbstractSciMLOperator}, b, u, Pl, Pr,
-        maxiters::Int,
-        abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions)
-    nothing, nothing
-end
-
 function LinearSolve.init_cacheval(alg::RFLUFactorization,
         A::Union{Diagonal, SymTridiagonal, Tridiagonal}, b, u, Pl, Pr,
         maxiters::Int,
@@ -111,7 +104,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::LUFactorization; kwargs...)
     A = convert(AbstractMatrix, A)
     if cache.isfresh
         cacheval = @get_cacheval(cache, :LUFactorization)
-        if A isa AbstractSparseMatrix && alg.reuse_symbolic
+        if issparsematrix(A) && alg.reuse_symbolic
             # Caches the symbolic factorization: https://github.com/JuliaLang/julia/pull/33738
             # If SparseMatrixCSC, check if the pattern has changed
             if alg.check_pattern && pattern_changed(cacheval, A)
@@ -248,13 +241,6 @@ function init_cacheval(alg::QRFactorization, A::Symmetric{<:Number, <:Array}, b,
     return qr(convert(AbstractMatrix, A), alg.pivot)
 end
 
-function init_cacheval(
-        alg::QRFactorization, A::Symmetric{<:Number, <:SparseMatrixCSC}, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Bool,
-        assumptions::OperatorAssumptions)
-    return nothing
-end
-
 const PREALLOCATED_QR_ColumnNorm = ArrayInterface.qr_instance(rand(1, 1), ColumnNorm())
 
 function init_cacheval(alg::QRFactorization{ColumnNorm}, A::Matrix{Float64}, b, u, Pl, Pr,

src/iterative_wrappers.jl

@@ -195,16 +195,16 @@ function init_cacheval(alg::KrylovJL, A, b, u, Pl, Pr, maxiters::Int, abstol, re
                      alg.KrylovAlg === Krylov.fgmres! ||
                      alg.KrylovAlg === Krylov.gpmr! ||
                      alg.KrylovAlg === Krylov.fom!)
-            if A isa SparseMatrixCSC
-                KS(SparseMatrixCSC(0, 0, [1], Int[], eltype(A)[]), eltype(b)[], 1)
+            if issparsematrixcsc(A)
+                KS(makeempty_SparaseMatrixCSC(A), eltype(b)[], 1)
             elseif A isa Matrix
                 KS(Matrix{eltype(A)}(undef, 0, 0), eltype(b)[], 1)
             else
                 KS(A, b, 1)
             end
         else
-            if A isa SparseMatrixCSC
-                KS(SparseMatrixCSC(0, 0, [1], Int[], eltype(A)[]), eltype(b)[])
+            if issparsematrixcsc(A)
+                KS(makeempty_SparaseMatrixCSC(A), eltype(b)[])
             elseif A isa Matrix
                 KS(Matrix{eltype(A)}(undef, 0, 0), eltype(b)[])
             else

test/basictests.jl

@@ -1,5 +1,5 @@
 using LinearSolve, LinearAlgebra, SparseArrays, MultiFloats, ForwardDiff
-using SciMLOperators, RecursiveFactorization
+using SciMLOperators, RecursiveFactorization, Sparspak
 using IterativeSolvers, KrylovKit, MKL_jll, KrylovPreconditioners
 using Test
 import Random