make sure to init the default algorithm and move first factorization

Author: ChrisRackauckas

src/common.jl

@@ -65,7 +65,7 @@ function set_cacheval(cache::LinearCache, alg_cache)
     return cache
 end
 
-init_cacheval(alg::Union{SciMLLinearSolveAlgorithm,Nothing}, A, b, u) = nothing
+init_cacheval(alg::SciMLLinearSolveAlgorithm, A, b, u) = nothing
 
 SciMLBase.init(prob::LinearProblem, args...; kwargs...) = SciMLBase.init(prob,nothing,args...;kwargs...)
 

src/default.jl

@@ -56,3 +56,58 @@ function SciMLBase.solve(cache::LinearCache, alg::Nothing,
         SciMLBase.solve(cache, alg, args...; kwargs...)
     end
 end
+
+function init_cacheval(alg::Nothing, A, b, u)
+    if A isa DiffEqArrayOperator
+        A = A.A
+    end
+
+    # Special case on Arrays: avoid BLAS for RecursiveFactorization.jl when
+    # it makes sense according to the benchmarks, which is dependent on
+    # whether MKL or OpenBLAS is being used
+    if A isa Matrix
+        if eltype(A) <: Union{Float32,Float64,ComplexF32,ComplexF64} &&
+                    ArrayInterface.can_setindex(b) && (size(A,1) <= 100 ||
+                                              (isopenblas() && size(A,1) <= 500)
+                                             )
+            alg = RFLUFactorization()
+            init_cacheval(alg, A, b, u)
+        else
+            alg = LUFactorization()
+            init_cacheval(alg, A, b, u)
+        end
+
+    # These few cases ensure the choice is optimal without the
+    # dynamic dispatching of factorize
+    elseif A isa Tridiagonal
+        alg = GenericFactorization(;fact_alg=lu!)
+        init_cacheval(alg, A, b, u)
+    elseif A isa SymTridiagonal
+        alg = GenericFactorization(;fact_alg=ldlt!)
+        init_cacheval(alg, A, b, u)
+    elseif A isa SparseMatrixCSC
+        alg = UMFPACKFactorization()
+        init_cacheval(alg, A, b, u)
+
+    # This catches the cases where a factorization overload could exist
+    # For example, BlockBandedMatrix
+    elseif ArrayInterface.isstructured(A)
+        alg = GenericFactorization()
+        init_cacheval(alg, A, b, u)
+
+    # This catches the case where A is a CuMatrix
+    # Which does not have LU fully defined
+    elseif !(A isa AbstractDiffEqOperator)
+        alg = QRFactorization()
+        init_cacheval(alg, A, b, u)
+
+    # Not factorizable operator, default to only using A*x
+    # IterativeSolvers is faster on CPU but not GPU-compatible
+    elseif cache.u isa Array
+        alg = IterativeSolversJL_GMRES()
+        init_cacheval(alg, A, b, u)
+    else
+        alg = KrylovJL_GMRES()
+        init_cacheval(alg, A, b, u)
+    end
+end