Change default to use KrylovJL_GMRES on non Array

Because this is GPU-compatible. IterativeSolvers.jl does some low rank Q updating which might make it more efficient on CPU but precludes it from being used on the GPU, so if that's the case then we use IterativeSolvers for the speed but more generally fall back to Krylov.jl

Author: ChrisRackauckas

src/default.jl

@@ -7,6 +7,9 @@ function SciMLBase.solve(cache::LinearCache, alg::Nothing,
         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 ArrayInterface.can_setindex(cache.b) && (size(A,1) <= 100 ||
                                               (isopenblas() && size(A,1) <= 500)
@@ -17,6 +20,9 @@ function SciMLBase.solve(cache::LinearCache, alg::Nothing,
             alg = LUFactorization()
             SciMLBase.solve(cache, alg, args...; kwargs...)
         end
+    
+    # These few cases ensure the choice is optimal without the
+    # dynamic dispatching of factorize
     elseif A isa Tridiagonal
         alg = GenericFactorization(;fact_alg=lu!)
         SciMLBase.solve(cache, alg, args...; kwargs...)
@@ -26,14 +32,26 @@ function SciMLBase.solve(cache::LinearCache, alg::Nothing,
     elseif A isa SparseMatrixCSC
         alg = LUFactorization()
         SciMLBase.solve(cache, alg, args...; kwargs...)
+    
+    # This catches the cases where a factorization overload could exist
+    # For example, BlockBandedMatrix
     elseif ArrayInterface.isstructured(A)
         alg = GenericFactorization()
         SciMLBase.solve(cache, alg, args...; kwargs...)
+    
+    # This catches the case where A is a CuMatrix
+    # Which does not have LU fully defined
     elseif !(A isa AbstractDiffEqOperator)
         alg = QRFactorization()
         SciMLBase.solve(cache, alg, args...; kwargs...)
-    else
+    
+    # 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()
         SciMLBase.solve(cache, alg, args...; kwargs...)
+    else
+        alg = KrylovJL_GMRES()
+        SciMLBase.solve(cache, alg, args...; kwargs...)
     end
 end