test cases

Author: vpuri3

src/LinearSolve.jl

@@ -48,8 +48,11 @@ isopenblas() = IS_OPENBLAS[]
 
 export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization,
        GenericLUFactorization, SimpleLUFactorization, RFLUFactorization,
-       UMFPACKFactorization, KLUFactorization,
-       FunctionCall
+       UMFPACKFactorization, KLUFactorization
+
+export FunctionCall, LdivBang2Args, LDivBang3Args,
+       ApplyLDivBang2Args, ApplyLDivBang3Args
+
 export KrylovJL, KrylovJL_CG, KrylovJL_GMRES, KrylovJL_BICGSTAB, KrylovJL_MINRES,
        IterativeSolversJL, IterativeSolversJL_CG, IterativeSolversJL_GMRES,
        IterativeSolversJL_BICGSTAB, IterativeSolversJL_MINRES,

src/function_call.jl

@@ -1,26 +1,46 @@
+
+""" user passes in inverse function, and arg symbols """
 struct FunctionCall{F,A} <: SciMLLinearSolveAlgorithm
     func!::F
-    argsyms::A
+    argnames::A
 
-    function FunctionCall(func!::Function, argsyms::Tuple)
-        new{typeof(func!), typeof(argsyms)}(func!, argsyms)
+    function FunctionCall(func!::Function, argnames::Tuple)
+#       @assert hasfield(::LinearCache, argnames)
+#       @assert isdefined
+        new{typeof(func!), typeof(argnames)}(func!, argnames)
     end
 end
 
-LdivBang2Args() = FunctionCall(ldiv!, (:A, :u))
-LdivBang3Args() = FunctionCall(ldiv!, (:u, :A, :b))
-
 function (f::FunctionCall)(cache::LinearCache)
-    @unpack func!, argsyms = f
-    args = [getproperty(cache,argsym) for argsym in argsyms]
+    @unpack func!, argnames = f
+    args = [getproperty(cache,argname) for argname in argnames]
     func!(args...)
 end
 
 function SciMLBase.solve(cache::LinearCache, alg::FunctionCall,
                          args...; kwargs...)
     @unpack u, b = cache
-    copy!(u, b)
     alg(cache)
+    return SciMLBase.build_linear_solution(alg,cache.u,nothing,cache)
+end
+
+##
 
+""" apply ldiv!(A, u) """
+struct ApplyLDivBang2Args <: SciMLLinearSolveAlgorithm end
+function SciMLBase.solve(cache::LinearCache, alg::ApplyLDivBang2Args,
+                         args...; kwargs...)
+    @unpack A, b, u = cache
+    copy!(u, b)
+    LinearAlgebra.ldiv!(A, u)
+    return SciMLBase.build_linear_solution(alg,cache.u,nothing,cache)
+end
+
+""" apply ldiv!(u, A, b) """
+struct ApplyLDivBang3Args <: SciMLLinearSolveAlgorithm end
+function SciMLBase.solve(cache::LinearCache, alg::ApplyLDivBang3Args,
+                         args...; kwargs...)
+    @unpack A, b, u = cache
+    LinearAlgebra.ldiv!(u, A, b)
     return SciMLBase.build_linear_solution(alg,cache.u,nothing,cache)
 end

test/basictests.jl

@@ -37,6 +37,64 @@ 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)),
+                    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)),
+                    ApplyLDivBang2Args(),
+                    ApplyLDivBang3Args(),
+                   )
+            test_interface(alg, prob1, prob2)
+        end
+    end
+
+end
+#=
+
 @testset "Default Linear Solver" begin
     test_interface(nothing, prob1, prob2)
 
@@ -286,5 +344,6 @@ end
     @test sol13.u ≈ sol23.u
     @test sol13.u ≈ sol33.u
 end
+=#
 
 end # testset