Hotfix and add tests for default linear solver issues found in ODE

Tested locally

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "LinearSolve"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 authors = ["SciML"]
-version = "2.1.2"
+version = "2.1.3"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"

src/default.jl

@@ -25,12 +25,12 @@ defaultalg(A, b) = defaultalg(A, b, OperatorAssumptions(true))
 
 function defaultalg(A::Union{DiffEqArrayOperator, MatrixOperator}, b,
                     assump::OperatorAssumptions)
-    DefaultLinearSolver(defaultalg(A.A, b, assump))
+    defaultalg(A.A, b, assump)
 end
 
 function defaultalg(A, b, assump::OperatorAssumptions{Nothing})
     issq = issquare(A)
-    DefaultLinearSolver(defaultalg(A, b, OperatorAssumptions(issq, assump.condition)))
+    defaultalg(A, b, OperatorAssumptions(issq, assump.condition))
 end
 
 function defaultalg(A::Tridiagonal, b, assump::OperatorAssumptions)

src/iterative_wrappers.jl

@@ -161,6 +161,8 @@ function get_KrylovJL_solver(KrylovAlg)
         Krylov.GpmrSolver
     elseif (KrylovAlg === Krylov.fom!)
         Krylov.FomSolver
+    else
+        error("Invalid Krylov method detected")
     end
 
     return KS
@@ -241,7 +243,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; kwargs...)
     itmax = cache.maxiters
     verbose = cache.verbose ? 1 : 0
 
-    args = (cache.cacheval, cache.A, cache.b)
+    args = (@get_cacheval(cache, :KrylovJL_GMRES), cache.A, cache.b)
     kwargs = (atol = atol, rtol = rtol, itmax = itmax, verbose = verbose,
               ldiv = true, history = true, alg.kwargs...)
 

test/basictests.jl

@@ -403,7 +403,7 @@ end
         end
 
         @testset "DirectLdiv!" begin
-            function get_operator(A, u)
+            function get_operator(A, u; add_inverse = true)
                 function f(du, u, p, t)
                     println("using FunctionOperator mul!")
                     mul!(du, A, u)
@@ -414,22 +414,35 @@ end
                     ldiv!(du, A, u)
                 end
 
-                FunctionOperator(f, u, u; isinplace = true, op_inverse = fi)
+                if add_inverse
+                    FunctionOperator(f, u, u; isinplace = true, op_inverse = fi)
+                else
+                    FunctionOperator(f, u, u; isinplace = true)
+                end
             end
 
             op1 = get_operator(A1, x1 * 0)
             op2 = get_operator(A2, x2 * 0)
+            op3 = get_operator(A1, x1 * 0; add_inverse = false)
+            op4 = get_operator(A2, x2 * 0; add_inverse = false)
 
             prob1 = LinearProblem(op1, b1; u0 = x1)
             prob2 = LinearProblem(op2, b2; u0 = x2)
+            prob3 = LinearProblem(op1, b1; u0 = x1)
+            prob4 = LinearProblem(op2, b2; u0 = x2)
 
             @test LinearSolve.defaultalg(op1, x1).alg ===
                   LinearSolve.DefaultAlgorithmChoice.DirectLdiv!
             @test LinearSolve.defaultalg(op2, x2).alg ===
                   LinearSolve.DefaultAlgorithmChoice.DirectLdiv!
-
+            @test LinearSolve.defaultalg(op3, x1).alg ===
+                  LinearSolve.DefaultAlgorithmChoice.KrylovJL_GMRES
+            @test LinearSolve.defaultalg(op4, x2).alg ===
+                  LinearSolve.DefaultAlgorithmChoice.KrylovJL_GMRES
             test_interface(DirectLdiv!(), prob1, prob2)
             test_interface(nothing, prob1, prob2)
+            test_interface(KrylovJL_GMRES(), prob3, prob4)
+            test_interface(nothing, prob3, prob4)
         end
     end
 end # testset