Merge branch 'main' into vp-precond

Author: vpuri3

.github/workflows/CI.yml

@@ -11,6 +11,7 @@ jobs:
       matrix:
         version:
           - '1'
+          - '1.6'
         os:
           - ubuntu-latest
         arch:

Project.toml

@@ -1,7 +1,7 @@
 name = "LinearSolve"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 authors = ["Jonathan <[email protected]> and contributors"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
@@ -26,7 +26,7 @@ Reexport = "1"
 SciMLBase = "1.18.6"
 Setfield = "0.7, 0.8"
 UnPack = "1"
-julia = "1"
+julia = "1.6"
 
 [extras]
 DiffEqProblemLibrary = "a077e3f3-b75c-5d7f-a0c6-6bc4c8ec64a9"

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