iterativesolvers.jl wrapper

Author: vpuri3

src/wrappers.jl

@@ -10,6 +10,10 @@ function KrylovJL(args...; KrylovAlg = Krylov.gmres!, kwargs...)
     return KrylovJL(KrylovAlg, args, kwargs)
 end
 
+KrylovJL_CG(args...;kwargs...) = KrylovJL(Krylov.cg!, args...; kwargs...)
+KrylovJL_GMRES(args...;kwargs...) = KrylovJL(Krylov.gmres!, args...; kwargs...)
+KrylovJL_BICGSTAB(args...;kwargs...) = KrylovJL(Krylov.bicgstab!, args...; kwargs...)
+
 function init_cacheval(alg::KrylovJL, A, b, u)
     cacheval = if alg.KrylovAlg === Krylov.cg!
         Krylov.CgSolver(A,b)
@@ -36,34 +40,48 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
     return cache.u
 end
 
-KrylovJL_CG(args...;kwargs...) = KrylovJL(Krylov.cg!, args...; kwargs...)
-KrylovJL_GMRES(args...;kwargs...) = KrylovJL(Krylov.gmres!, args...; kwargs...)
-KrylovJL_BICGSTAB(args...;kwargs...) = KrylovJL(Krylov.bicgstab!, args...; kwargs...)
-
 ## IterativeSolvers.jl
 
 struct IterativeSolversJL{F,A,K} <: SciMLLinearSolveAlgorithm
-    solver::F
+    generate_iterator::F
     args::A
     kwargs::K
 end
 
-## KrylovKit.jl
-
-struct KrylovKitJL{F,A,K} <: SciMLLinearSolveAlgorithm
-    solver::F
-    args::A
-    kwargs::K
+function IterativeSolversJL(args...;
+                            generate_iterator = IterativeSolvers.gmres_iterable!,
+                            kwargs...)
+    return IterativeSolversJL(generate_iterator, args, kwargs)
 end
 
-function KrylovKitJL(args...; solver = KrylovKit.CG(), kwargs...)
-    return KrylovKitJL(solver, args, kwargs)
+#IterativeSolversJL_CG(args...; kwargs...)
+#    = IterativeSolversJL(IterativeSolvers.cg_iterator!, args...; kwargs...)
+#IterativeSolversJL_GMRES(args...;kwargs...)
+#    = IterativeSolversJL(IterativeSolvers.gmres_iterable!, args...; kwargs...)
+#IterativeSolversJL_BICGSTAB(args...;kwargs...)
+#    = IterativeSolversJL(IterativeSolvers.bicgstabl_iterator!, args...;kwargs...)
+
+function init_cacheval(alg::IterativeSolversJL, A, b, u)
+    cacheval = if alg.generate_iterator === IterativeSolvers.cg_iterator!
+        alg.generate_iterator(u, A, b)
+    elseif alg.generate_iterator === IterativeSolvers.gmres_iterable!
+        alg.generate_iterator(u, A, b)
+    elseif alg.generate_iterator === IterativeSolvers.bicgstabl_iterator!
+        alg.generate_iterator(u, A, b)
+    else
+        alg.generate_iterator(u, A, b)
+    end
+    return cacheval
 end
 
-function SciMLBase.solve(cache::LinearCache, alg::KrylovKitJL,args...;kwargs...)
-    @unpack A, b, u = cache
-    @unpack solver = alg
-    u = KrylovKit.linsolve(A, b, u, solver, args...; kwargs...)[1] #no precond?!
-    return u
+function SciMLBase.solve(cache::LinearCache, alg::IterativeSolversJL; kwargs...)
+    if cache.isfresh
+        solver = init_cacheval(alg, cache.A, cache.b, cache.u)
+        cache = set_cacheval(cache, solver)
+    end
+
+    for resi in cache.cacheval end
+
+    return cache.u
 end
 

test/runtests.jl

@@ -3,13 +3,13 @@ using Test
 
 @testset "LinearSolve.jl" begin
     using LinearAlgebra
-    n = 2
+    n = 32
 
     A = Matrix(I,n,n)
     b = ones(n)
-    A1 = A/1; b1 = ones(n); x1 = zero(b)
-    A2 = A/2; b2 = ones(n); x2 = zero(b)
-    A3 = A/3; b3 = ones(n); x3 = zero(b)
+    A1 = A/1; b1 = rand(n); x1 = zero(b)
+    A2 = A/2; b2 = rand(n); x2 = zero(b)
+    A3 = A/3; b3 = rand(n); x3 = zero(b)
 
     prob1 = LinearProblem(A1, b1; u0=x1)
     prob2 = LinearProblem(A2, b2; u0=x2)
@@ -23,9 +23,7 @@ using Test
 #               :DefaultLinSolve,
 
                 :KrylovJL, :KrylovJL_CG, :KrylovJL_GMRES, :KrylovJL_BICGSTAB,
-#               :IterativeSolversJL
-#               :KrylovKitJL,
-
+                :IterativeSolversJL,#:IterativeSolversJL_GMRES, :IterativeSolversJL_BICGSTAB,
                )
         @eval begin
             y = solve($prob1, $alg())