add methods defaultalg to support scimlops

Author: vpuri3

src/common.jl

@@ -25,6 +25,7 @@ end
 """
 $(SIGNATURES)
 """
+# TODO - this should modify OperatorAssumption??
 function set_A(cache::LinearCache, A)
     @set! cache.A = A
     @set! cache.isfresh = true
@@ -94,7 +95,7 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
                         verbose::Bool = false,
                         Pl = IdentityOperator{size(prob.A, 1)}(),
                         Pr = IdentityOperator{size(prob.A, 2)}(),
-                        assumptions = OperatorAssumptions(),
+                        assumptions = OperatorAssumptions(issquare(prob.A)),
                         kwargs...)
     @unpack A, b, u0, p = prob
 

src/default.jl

@@ -36,6 +36,7 @@ end
 function defaultalg(A::Diagonal, b, ::OperatorAssumptions{true})
     DiagonalFactorization()
 end
+# TODO - Diagonal matrices are always square
 function defaultalg(A::Diagonal, b, ::OperatorAssumptions{false})
     DiagonalFactorization()
 end
@@ -76,13 +77,21 @@ function defaultalg(A, b::GPUArraysCore.AbstractGPUArray, ::OperatorAssumptions{
 end
 
 function defaultalg(A::SciMLBase.AbstractSciMLOperator, b,
-                    assumptions::OperatorAssumptions)
+                    assumptions::OperatorAssumptions{true})
+    if has_ldiv!(A)
+        return DirectLdiv!()
+    end
+
     KrylovJL_GMRES()
 end
 
 # Ambiguity handling
 function defaultalg(A::SciMLBase.AbstractSciMLOperator, b,
                     assumptions::OperatorAssumptions{Nothing})
+    if has_ldiv!(A)
+        return DirectLdiv!()
+    end
+
     KrylovJL_GMRES()
 end
 

test/basictests.jl

@@ -27,20 +27,20 @@ function test_interface(alg, prob1, prob2)
     b2 = prob2.b
     x2 = prob2.u0
 
-    y = solve(prob1, alg; cache_kwargs...)
-    @test A1 * y ≈ b1
+    sol = solve(prob1, alg; cache_kwargs...)
+    @test A1 * sol.u ≈ b1
 
     cache = SciMLBase.init(prob1, alg; cache_kwargs...) # initialize cache
-    y = solve(cache)
-    @test A1 * y ≈ b1
+    sol = solve(cache)
+    @test A1 * sol.u ≈ b1
 
-    cache = LinearSolve.set_A(cache, copy(A2))
-    y = solve(cache; cache_kwargs...)
-    @test A2 * y ≈ b1
+    cache = LinearSolve.set_A(cache, deepcopy(A2))
+    sol = solve(cache; cache_kwargs...)
+    @test A2 * sol.u ≈ b1
 
     cache = LinearSolve.set_b(cache, b2)
-    y = solve(cache; cache_kwargs...)
-    @test A2 * y ≈ b2
+    sol = solve(cache; cache_kwargs...)
+    @test A2 * sol.u ≈ b2
 
     return
 end
@@ -359,26 +359,58 @@ end
         b2 = rand(n)
         x2 = zero(b1)
 
-        function sol_func(A, b, u, p, newA, Pl, Pr, solverdata; verbose = true, kwargs...)
-            if verbose == true
-                println("out-of-place solve")
+        @testset "LinearSolveFunction" begin
+            function sol_func(A, b, u, p, newA, Pl, Pr, solverdata; verbose = true, kwargs...)
+                if verbose == true
+                    println("out-of-place solve")
+                end
+                u = A \ b
             end
-            u = A \ b
-        end
 
-        function sol_func!(A, b, u, p, newA, Pl, Pr, solverdata; verbose = true, kwargs...)
-            if verbose == true
-                println("in-place solve")
+            function sol_func!(A, b, u, p, newA, Pl, Pr, solverdata; verbose = true, kwargs...)
+                if verbose == true
+                    println("in-place solve")
+                end
+                ldiv!(u, A, b)
+            end
+
+            prob1 = LinearProblem(A1, b1; u0 = x1)
+            prob2 = LinearProblem(A1, b1; u0 = x1)
+
+            for alg in (LinearSolveFunction(sol_func),
+                        LinearSolveFunction(sol_func!))
+                test_interface(alg, prob1, prob2)
             end
-            ldiv!(u, A, b)
         end
 
-        prob1 = LinearProblem(A1, b1; u0 = x1)
-        prob2 = LinearProblem(A1, b1; u0 = x1)
+        @testset "DirectLdiv!" begin
+            function get_operator(A, u)
+                F = lu(A)
 
-        for alg in (LinearSolveFunction(sol_func),
-                    LinearSolveFunction(sol_func!))
-            test_interface(alg, prob1, prob2)
+                function f(du, u, p, t)
+                    println("using FunctionOperator mul!")
+                    mul!(du, A, u)
+                end
+
+                function fi(du, u, p, t)
+                    println("using FunctionOperator ldiv!")
+                    ldiv!(du, F, u)
+                end
+
+                FunctionOperator(f, u, u; isinplace=true, op_inverse=fi)
+            end
+
+            op1 = get_operator(A1, x1*0)
+            op2 = get_operator(A2, x2*0)
+
+            prob1 = LinearProblem(op1, b1; u0 = x1)
+            prob2 = LinearProblem(op2, b2; u0 = x2)
+
+            @test LinearSolve.defaultalg(op1, x1) isa DirectLdiv!
+            @test LinearSolve.defaultalg(op2, x2) isa DirectLdiv!
+
+            test_interface(DirectLdiv!(), prob1, prob2)
+            test_interface(nothing, prob1, prob2)
         end
     end
 end # testset