impose function syntax

Author: vpuri3

src/LinearSolve.jl

@@ -26,13 +26,13 @@ using Reexport
 abstract type SciMLLinearSolveAlgorithm <: SciMLBase.AbstractLinearAlgorithm end
 abstract type AbstractFactorization <: SciMLLinearSolveAlgorithm end
 abstract type AbstractKrylovSubspaceMethod <: SciMLLinearSolveAlgorithm end
-abstract type AbstractFunctionCall <: SciMLLinearSolveAlgorithm end
+abstract type AbstractSolveFunction <: SciMLLinearSolveAlgorithm end
 
 # Traits
 
 needs_concrete_A(alg::AbstractFactorization) = true
 needs_concrete_A(alg::AbstractKrylovSubspaceMethod) = false
-needs_concrete_A(alg::AbstractFunctionCall) = false
+needs_concrete_A(alg::AbstractSolveFunction) = false
 
 # Code
 
@@ -41,7 +41,7 @@ include("factorization.jl")
 include("simplelu.jl")
 include("iterative_wrappers.jl")
 include("preconditioners.jl")
-include("function_call.jl")
+include("solve_function.jl")
 include("default.jl")
 include("init.jl")
 
@@ -52,7 +52,7 @@ export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization,
        GenericLUFactorization, SimpleLUFactorization, RFLUFactorization,
        UMFPACKFactorization, KLUFactorization
 
-export FunctionCall, ApplyLDivBang, ApplyLDivBang2Args, ApplyLDivBang3Args
+export LinearSolveFunction
 
 export KrylovJL, KrylovJL_CG, KrylovJL_GMRES, KrylovJL_BICGSTAB, KrylovJL_MINRES,
        IterativeSolversJL, IterativeSolversJL_CG, IterativeSolversJL_GMRES,

src/function_call.jl

@@ -0,0 +1,19 @@
+#
+function DEFAULT_LINEAR_SOLVE(A,b,u,p,newA,Pl,Pr,solverdata;kwargs...)
+    solve(LinearProblem(A, b; u0=u); p=p, kwargs...).u
+end
+
+Base.@kwdef struct LinearSolveFunction{F} <: AbstractSolveFunction
+    solve_func::F = DEFAULT_LINEAR_SOLVE
+end
+
+function SciMLBase.solve(cache::LinearCache, alg::LinearSolveFunction,
+                         args...; kwargs...)
+    @unpack A,b,u,p,isfresh,Pl,Pr,cacheval = cache
+    @unpack solve_func = alg
+
+    u = solve_func(A,b,u,p,isfresh,Pl,Pr,cacheval;kwargs...)
+    cache = set_u(cache, u)
+
+    return SciMLBase.build_linear_solution(alg,cache.u,nothing,cache)
+end

src/solve_function.jl

@@ -37,65 +37,6 @@ end
 
 @testset "LinearSolve" begin
 
-@testset "Apply Function" begin
-
-    @testset "Diagonal Type" begin
-        A1 = rand(n) |> Diagonal; b1 = rand(n); x1 = zero(b1)
-        A2 = rand(n) |> Diagonal; b2 = rand(n); x2 = zero(b1)
-                                                  
-        prob1 = LinearProblem(A1, b1; u0=x1)
-        prob2 = LinearProblem(A1, b1; u0=x1)
-
-        for alg in (
-                    FunctionCall(LinearAlgebra.ldiv!, (:u, :A, :b)),
-                    ApplyLDivBang(),
-                    ApplyLDivBang2Args(),
-                    ApplyLDivBang3Args(),
-                   )
-            test_interface(alg, prob1, prob2)
-        end
-    end
-
-    @testset "Custom Type" begin
-
-        struct MyDiag
-            d
-        end
-
-        # overloads
-        (D::MyDiag)(du, u, p, t) = mul!(du, D, u)
-        Base.:*(D::MyDiag, u) = Diagonal(D.d) * u
-
-        Base.copy(D::MyDiag) = copy(D.d) |> MyDiag
-
-        LinearAlgebra.mul!(y, D::MyDiag, x) = mul!(y, Diagonal(D.d), x)
-        LinearAlgebra.ldiv!(y, D::MyDiag, x) = ldiv!(y, Diagonal(D.d), x)
-        LinearAlgebra.ldiv!(D::MyDiag, x) = ldiv!(Diagonal(D.d), x)
-
-        # custom inverse function
-        function my_inv!(D::MyDiag, u, b)
-            @. u = b / D.d 
-        end
-    
-        A1 = rand(n) |> MyDiag; b1 = rand(n); x1 = zero(b1)
-        A2 = rand(n) |> MyDiag; b2 = rand(n); x2 = zero(b1)
-                                                  
-        prob1 = LinearProblem(A1, b1; u0=x1)
-        prob2 = LinearProblem(A1, b1; u0=x1)
-
-        for alg in (
-                    FunctionCall(LinearAlgebra.ldiv!, (:u, :A, :b)),
-                    FunctionCall(my_inv!, (:A, :u, :b)),
-                    ApplyLDivBang(),
-                    ApplyLDivBang2Args(),
-                    ApplyLDivBang3Args(),
-                   )
-            test_interface(alg, prob1, prob2)
-        end
-    end
-end
-#=
-
 @testset "Default Linear Solver" begin
     test_interface(nothing, prob1, prob2)
 
@@ -345,6 +286,36 @@ end
     @test sol13.u ≈ sol23.u
     @test sol13.u ≈ sol33.u
 end
-=#
+
+@testset "Solve Function" begin
+
+    A1 = rand(n) |> Diagonal; b1 = rand(n); x1 = zero(b1)
+    A2 = rand(n) |> Diagonal; 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")
+        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")
+        end
+        ldiv!(u,A,b)
+    end
+
+    prob1 = LinearProblem(A1, b1; u0=x1)
+    prob2 = LinearProblem(A1, b1; u0=x1)
+
+    for alg in (
+                LinearSolveFunction(),
+                LinearSolveFunction(sol_func),
+                LinearSolveFunction(sol_func!),
+               )
+        test_interface(alg, prob1, prob2)
+    end
+end
 
 end # testset