Merge pull request #307 from SciML/defaultalg

Rework default algorithm to be fully type stable
SciML · May 31, 2023 · 9718ea2 · 9718ea2
2 parents cb31d58 + 1eccc20
commit 9718ea2
Show file tree

Hide file tree

Showing 11 changed files with 919 additions and 290 deletions.
diff --git a/Project.toml b/Project.toml
@@ -27,7 +27,7 @@ SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
-ArrayInterface = "7.3"
+ArrayInterface = "7.4.7"
 DocStringExtensions = "0.8, 0.9"
 EnumX = "1"
 FastLapackInterface = "1"

diff --git a/docs/src/solvers/solvers.md b/docs/src/solvers/solvers.md
@@ -89,6 +89,9 @@ customized per-package, details given below describe a subset of important array
         When `A` is sparse, a similar polyalgorithm is used. For indefinite matrices, the `LDLt`
         factorization does not use pivoting during the numerical factorization and therefore the
         procedure can fail even for invertible matrices.
+  - CholeskyFactorization
+  - BunchKaufmanFactorization
+  - CHOLMODFactorization
 
 ### LinearSolve.jl
 

diff --git a/ext/LinearSolveHYPREExt.jl b/ext/LinearSolveHYPREExt.jl
@@ -88,8 +88,8 @@ function SciMLBase.init(prob::LinearProblem, alg::HYPREAlgorithm,
 
     cache = LinearCache{
                         typeof(A), typeof(b), typeof(u0), typeof(p), typeof(alg), Tc,
-                        typeof(Pl), typeof(Pr), typeof(reltol), __issquare(assumptions),
-                        __conditioning(assumptions)
+                        typeof(Pl), typeof(Pr), typeof(reltol),
+                        typeof(__issquare(assumptions))
                         }(A, b, u0, p, alg, cacheval, isfresh, Pl, Pr, abstol, reltol,
                           maxiters,
                           verbose, assumptions)

diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl
@@ -77,6 +77,7 @@ isopenblas() = IS_OPENBLAS[]
 
 import PrecompileTools
 
+#=
 PrecompileTools.@compile_workload begin
     A = rand(4, 4)
     b = rand(4)
@@ -104,12 +105,14 @@ PrecompileTools.@compile_workload begin
     sol = solve(prob) # in case sparspak is used as default
     sol = solve(prob, SparspakFactorization())
 end
+=#
 
 export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization,
        GenericLUFactorization, SimpleLUFactorization, RFLUFactorization,
        NormalCholeskyFactorization, NormalBunchKaufmanFactorization,
        UMFPACKFactorization, KLUFactorization, FastLUFactorization, FastQRFactorization,
-       SparspakFactorization, DiagonalFactorization
+       SparspakFactorization, DiagonalFactorization, CholeskyFactorization,
+       BunchKaufmanFactorization, CHOLMODFactorization, LDLtFactorization
 
 export LinearSolveFunction, DirectLdiv!
 

diff --git a/src/common.jl b/src/common.jl
@@ -53,17 +53,19 @@ end
 
 Sets the operator `A` assumptions used as part of the default algorithm
 """
-struct OperatorAssumptions{issq, condition} end
+struct OperatorAssumptions{T}
+    issq::T
+    condition::OperatorCondition.T
+end
+
 function OperatorAssumptions(issquare = nothing;
                              condition::OperatorCondition.T = OperatorCondition.IllConditioned)
-    issq = something(_unwrap_val(issquare), Nothing)
-    condition = _unwrap_val(condition)
-    OperatorAssumptions{issq, condition}()
+    OperatorAssumptions{typeof(issquare)}(issquare, condition)
 end
-__issquare(::OperatorAssumptions{issq, condition}) where {issq, condition} = issq
-__conditioning(::OperatorAssumptions{issq, condition}) where {issq, condition} = condition
+__issquare(assump::OperatorAssumptions) = assump.issq
+__conditioning(assump::OperatorAssumptions) = assump.condition
 
-mutable struct LinearCache{TA, Tb, Tu, Tp, Talg, Tc, Tl, Tr, Ttol, issq, condition}
+mutable struct LinearCache{TA, Tb, Tu, Tp, Talg, Tc, Tl, Tr, Ttol, issq}
     A::TA
     b::Tb
     u::Tu
@@ -77,12 +79,15 @@ mutable struct LinearCache{TA, Tb, Tu, Tp, Talg, Tc, Tl, Tr, Ttol, issq, conditi
     reltol::Ttol
     maxiters::Int
     verbose::Bool
-    assumptions::OperatorAssumptions{issq, condition}
+    assumptions::OperatorAssumptions{issq}
 end
 
 function Base.setproperty!(cache::LinearCache, name::Symbol, x)
     if name === :A
         setfield!(cache, :isfresh, true)
+    elseif name === :cacheval && cache.alg isa DefaultLinearSolver
+        @assert cache.cacheval isa DefaultLinearSolverInit
+        return setfield!(cache.cacheval, Symbol(cache.alg.alg), x)
     end
     setfield!(cache, name, x)
 end
@@ -116,11 +121,20 @@ function SciMLBase.init(prob::LinearProblem, alg::SciMLLinearSolveAlgorithm,
                         kwargs...)
     @unpack A, b, u0, p = prob
 
-    A = alias_A ? A : deepcopy(A)
+    A = if alias_A
+        A
+    elseif A isa Array || A isa SparseMatrixCSC
+        copy(A)
+    else
+        deepcopy(A)
+    end
+
     b = if b isa SparseArrays.AbstractSparseArray && !(A isa Diagonal)
         Array(b) # the solution to a linear solve will always be dense!
     elseif alias_b
         b
+    elseif b isa Array || b isa SparseMatrixCSC
+        copy(b)
     else
         deepcopy(b)
     end
@@ -147,8 +161,7 @@ function SciMLBase.init(prob::LinearProblem, alg::SciMLLinearSolveAlgorithm,
                         typeof(Pl),
                         typeof(Pr),
                         typeof(reltol),
-                        __issquare(assumptions),
-                        __conditioning(assumptions)
+                        typeof(assumptions.issq)
                         }(A,
                           b,
                           u0,