default

Author: vpuri3

Project.toml

@@ -9,9 +9,11 @@ IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]

src/LinearSolve.jl

@@ -1,8 +1,10 @@
 module LinearSolve
 
-using ArrayInterface: lu_instance
+using ArrayInterface
+using RecursiveFactorization
 using Base: cache_dependencies, Bool
 using LinearAlgebra
+using SparseArrays
 using SciMLBase: AbstractDiffEqOperator, AbstractLinearAlgorithm
 using Setfield
 using UnPack
@@ -20,15 +22,15 @@ abstract type AbstractFactorization <: SciMLLinearSolveAlgorithm end
 abstract type AbstractKrylovSubspaceMethod <: SciMLLinearSolveAlgorithm end
 
 include("common.jl")
-include("default.jl")
 include("factorization.jl")
 include("wrappers.jl")
+include("default.jl")
 
-export DefaultLinSolve
-export LUFactorization, SVDFactorization, QRFactorization
+export LUFactorization, SVDFactorization, QRFactorization, DefaultFactorization
 export KrylovJL, KrylovJL_CG, KrylovJL_GMRES, KrylovJL_BICGSTAB,
        KrylovJL_MINRES, 
        IterativeSolversJL, IterativeSolversJL_CG, IterativeSolversJL_GMRES,
        IterativeSolversJL_BICGSTAB, IterativeSolversJL_MINRES
+export DefaultLinSolve
 
 end

src/default.jl

@@ -1,15 +1,51 @@
-#
-mutable struct DefaultLinSolve <: SciMLLinearSolveAlgorithm
-    iterable
+## Default algorithm
+
+struct DefaultLinSolve{Ta} <: SciMLLinearSolveAlgorithm
+    linalg::Ta
+    ifopenblas::Union{Bool,Nothing}
+    isset::Bool # true => do nothing, false => find alg
 end
 
-DefaultLinSolve() = DefaultLinSolve(nothing, nothing, nothing)
+DefaultLinSolve() = DefaultLinSolve(nothing, nothing, true)
+
+function isopenblas()
+    @static if VERSION < v"1.7beta"
+        blas = BLAS.vendor()
+        blas == :openblas64 || blas == :openblas
+    else
+        occursin("openblas", BLAS.get_config().loaded_libs[1].libname)
+    end
+end
+
+function SciMLBase.solve(cache::LinearCache, alg::DefaultLinSolve,
+                         args...; kwargs...)
+    @unpack A = cache
+
+    if alg.isset
+      linalg = if A isa Matrix
+          if ArrayInterface.can_setindex(x) && (size(A,1) <= 100 ||
+                                                (p.openblas && size(A,1) <= 500)
+                                               )
+              DefaultFactorization(;fact_alg=:(RecursiveFactorization.lu!))
+          else
+              LUFactorization()
+          end
+      elseif A isa Union{Tridiagonal,} # ForwardSensitivityJacobian
+          DefaultFactorization(;fact_alg=lu!)
+      elseif A isa Union{SymTridiagonal}
+          DefaultFactorization(;fact_alg=ldlt!)
+      elseif A isa SparseMatrixCSC
+          LUFactorization()
+      elseif ArrayInterface.isstructured(A)
+          DefaultFactorization() # change fact_alg=LinearAlgebra.factorize
+      elseif !(A isa AbstractDiffEqOperator)
+          QRFactorization()
+      else
+          IterativeSolversJL_GMRES()
+      end
 
-function SciMLBase.solve(cache::LinearCache,
-                         alg::DefaultLinSolve,
-                         args...;kwargs...)
-    @unpack iterable = alg
-    @unpack A, b, u, cacheval = cache
+      @set! alg.linalg = linalg
+    end
 
-    return A \ u
+    SciMLBase.solve(cache, alg.linalg, args...; kwargs...)
 end

src/factorization.jl

@@ -74,13 +74,20 @@ end
 
 ## DefaultFactorization
 
-struct DefaultFactorization <: AbstractFactorization
+struct DefaultFactorization{F} <: AbstractFactorization
+    fact_alg::F
 end
 
+DefaultFactorization(
+                     # https://github.com/SciML/SciMLBase.jl/pull/119
+#                    ;fact_alg = LinearAlgebra.factorize
+                     ;fact_alg = lu!
+                    ) = DefaultFactorization(fact_alg)
+
 function init_cacheval(alg::DefaultFactorization, A, b, u)
     A isa Union{AbstractMatrix,AbstractDiffEqOperator} ||
         error("DefaultFactorization is not defined for $(typeof(A))")
 
-    fact = factorize(A)
+    fact = alg.fact_alg(A)
     return fact
 end