Merge pull request #404 from mschauer/idrs

Add IDR(s) from IterativeSolvers

Author: ChrisRackauckas

Project.toml

@@ -67,7 +67,7 @@ EnzymeCore = "0.5, 0.6"
 FastLapackInterface = "1, 2"
 GPUArraysCore = "0.1"
 HYPRE = "1.4.0"
-IterativeSolvers = "0.9.2"
+IterativeSolvers = "0.9.3"
 KLU = "0.3.0, 0.4"
 KernelAbstractions = "0.9"
 Krylov = "0.9"

docs/src/solvers/solvers.md

@@ -221,6 +221,7 @@ IterativeSolversJL_CG
 IterativeSolversJL_GMRES
 IterativeSolversJL_BICGSTAB
 IterativeSolversJL_MINRES
+IterativeSolversJL_IDRS
 IterativeSolversJL
 ```
 

ext/LinearSolveIterativeSolversExt.jl

@@ -27,6 +27,12 @@ function LinearSolve.IterativeSolversJL_GMRES(args...; kwargs...)
         generate_iterator = IterativeSolvers.gmres_iterable!,
         kwargs...)
 end
+function LinearSolve.IterativeSolversJL_IDRS(args...; kwargs...)
+    IterativeSolversJL(args...;
+        generate_iterator = IterativeSolvers.idrs_iterable!,
+        kwargs...)
+end
+
 function LinearSolve.IterativeSolversJL_BICGSTAB(args...; kwargs...)
     IterativeSolversJL(args...;
         generate_iterator = IterativeSolvers.bicgstabl_iterator!,
@@ -47,6 +53,7 @@ function LinearSolve.init_cacheval(alg::IterativeSolversJL, A, b, u, Pl, Pr, max
     reltol,
     verbose::Bool, assumptions::OperatorAssumptions)
     restart = (alg.gmres_restart == 0) ? min(20, size(A, 1)) : alg.gmres_restart
+    s = :idrs_s in keys(alg.kwargs) ? alg.kwargs.idrs_s : 4 # shadow space
 
     kwargs = (abstol = abstol, reltol = reltol, maxiter = maxiters,
         alg.kwargs...)
@@ -59,6 +66,14 @@ function LinearSolve.init_cacheval(alg::IterativeSolversJL, A, b, u, Pl, Pr, max
     elseif alg.generate_iterator === IterativeSolvers.gmres_iterable!
         alg.generate_iterator(u, A, b; Pl = Pl, Pr = Pr, restart = restart,
             kwargs...)
+    elseif alg.generate_iterator === IterativeSolvers.idrs_iterable!
+        !!LinearSolve._isidentity_struct(Pr) &&
+            @warn "$(alg.generate_iterator) doesn't support right preconditioning"
+        history = IterativeSolvers.ConvergenceHistory(partial=true)
+        history[:abstol] = abstol
+        history[:reltol] = reltol
+        IterativeSolvers.idrs_iterable!(history, u, A, b, s, Pl, abstol, reltol, maxiters; 
+            alg.kwargs...)
     elseif alg.generate_iterator === IterativeSolvers.bicgstabl_iterator!
         !!LinearSolve._isidentity_struct(Pr) &&
             @warn "$(alg.generate_iterator) doesn't support right preconditioning"
@@ -95,7 +110,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::IterativeSolversJL; kwargs...
     end
     cache.verbose && println()
 
-    resid = cache.cacheval.residual
+    resid = cache.cacheval isa IterativeSolvers.IDRSIterable ? cache.cacheval.R : cache.cacheval.residual
     if resid isa IterativeSolvers.Residual
         resid = resid.current
     end

src/LinearSolve.jl

@@ -217,7 +217,7 @@ export LinearSolveFunction, DirectLdiv!
 export KrylovJL, KrylovJL_CG, KrylovJL_MINRES, KrylovJL_GMRES,
     KrylovJL_BICGSTAB, KrylovJL_LSMR, KrylovJL_CRAIGMR,
     IterativeSolversJL, IterativeSolversJL_CG, IterativeSolversJL_GMRES,
-    IterativeSolversJL_BICGSTAB, IterativeSolversJL_MINRES,
+    IterativeSolversJL_BICGSTAB, IterativeSolversJL_MINRES, IterativeSolversJL_IDRS,
     KrylovKitJL, KrylovKitJL_CG, KrylovKitJL_GMRES
 
 export SimpleGMRES

src/extension_algs.jl

@@ -309,6 +309,21 @@ A wrapper over the IterativeSolvers.jl GMRES.
 """
 function IterativeSolversJL_GMRES end
 
+"""
+```julia
+IterativeSolversJL_IDRS(args...; Pl = nothing, kwargs...)
+```
+
+A wrapper over the IterativeSolvers.jl IDR(S).
+
+
+!!! note
+    
+    Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers`
+
+"""
+function IterativeSolversJL_IDRS end
+
 """
 ```julia
 IterativeSolversJL_BICGSTAB(args...; Pl = nothing, Pr = nothing, kwargs...)

test/basictests.jl

@@ -267,7 +267,8 @@ end
             kwargs = (; gmres_restart = 5)
             for alg in (("Default", IterativeSolversJL(kwargs...)),
                 ("CG", IterativeSolversJL_CG(kwargs...)),
-                ("GMRES", IterativeSolversJL_GMRES(kwargs...))
+                ("GMRES", IterativeSolversJL_GMRES(kwargs...)),
+                ("IDRS", IterativeSolversJL_IDRS(kwargs...))
                 #           ("BICGSTAB",IterativeSolversJL_BICGSTAB(kwargs...)),
                 #            ("MINRES",IterativeSolversJL_MINRES(kwargs...)),
             )