update

Author: vpuri3

src/LinearSolve.jl

@@ -5,6 +5,7 @@ using RecursiveFactorization
 using Base: cache_dependencies, Bool
 import Base: eltype, adjoint, inv
 using LinearAlgebra
+using IterativeSolvers:Identity
 using SparseArrays
 using SciMLBase: AbstractDiffEqOperator, AbstractLinearAlgorithm
 using Setfield

src/common.jl

@@ -60,20 +60,20 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
                         reltol=√eps(eltype(prob.A)),
                         maxiters=length(prob.b),
                         verbose=false,
-                        scale=one(eltype(prob.A)),
+                        Pl = nothing,
+                        Pr = nothing,
                         kwargs...,
                        )
     @unpack A, b, u0, p = prob
 
-    u0 = (u0 === nothing) ? zero(b) : u0
+    u0 = (u0 !== nothing) ? u0 : zero(b)   
+    Pl = (Pl !== nothing) ? Pl : Identity()
+    Pr = (Pr !== nothing) ? Pr : Identity()
 
     cacheval = init_cacheval(alg, A, b, u0)
     isfresh = cacheval === nothing
     Tc = isfresh ? Any : typeof(cacheval)
 
-    Pl = scaling_preconditioner(scale, true)
-    Pr = scaling_preconditioner(scale, false)
-
     A = alias_A ? A : deepcopy(A)
     b = alias_b ? b : deepcopy(b)
 

src/wrappers.jl

@@ -1,7 +1,7 @@
 
 ## Preconditioners
 
-scaling_preconditioner(s, isleft) = isleft ? I * s : I * (1/s)
+scaling_preconditioner(s) = I * s , I * (1/s)
 
 struct ComposePreconditioner{Ti,To}
     inner::Ti
@@ -41,6 +41,23 @@ function LinearAlgebra.mul!(y, A::InvComposePreconditioner, x)
     ldiv!(y, P, x)
 end
 
+function get_preconditioner(Pi, Po)
+
+    ifPi = Pi !== Identity()
+    ifPo = Po !== Identity()
+
+    P = 
+    if ifPi & ifPo
+        ComposePreconditioner(Pi, Po)
+    elseif ifPi | ifPo
+        ifPi ? Pi : Po
+    else
+        Identity()
+    end
+
+    return P
+end
+
 ## Krylov.jl
 
 struct KrylovJL{F,Tl,Tr,I,A,K} <: AbstractKrylovSubspaceMethod
@@ -53,10 +70,14 @@ struct KrylovJL{F,Tl,Tr,I,A,K} <: AbstractKrylovSubspaceMethod
     kwargs::K
 end
 
-function KrylovJL(args...; KrylovAlg = Krylov.gmres!, Pl=I, Pr=I,
+function KrylovJL(args...; KrylovAlg = Krylov.gmres!,
+                  Pl=nothing, Pr=nothing,
                   gmres_restart=0, window=0,
                   kwargs...)
 
+    Pl = (Pl === nothing) ? Identity() : Pl
+    Pr = (Pr === nothing) ? Identity() : Pr
+
     return KrylovJL(KrylovAlg, Pl, Pr, gmres_restart, window,
                     args, kwargs)
 end
@@ -144,16 +165,11 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
         cache = set_cacheval(cache, solver)
     end
 
-    M = I # left precond
-    N = I # right precond
+    M = get_preconditioner(alg.Pl, cache.Pl)
+    N = get_preconditioner(alg.Pr, cache.Pr)
 
-    if (cache.Pl != I) | (alg.Pl != I)
-        M = InvComposePreconditioner(alg.Pl, cache.Pl)
-    end
-
-    if (cache.Pr != I) | (alg.Pr != I)
-        N = InvComposePreconditioner(alg.Pr, cache.Pr)
-    end
+    M = (M === Identity()) ? I : inv(M)
+    N = (N === Identity()) ? I : inv(N)
 
     atol    = cache.abstol
     rtol    = cache.reltol
@@ -165,7 +181,7 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
               alg.kwargs...)
 
     if cache.cacheval isa Krylov.CgSolver
-        N != I  &&
+        N !== I  &&
             @warn "$(alg.KrylovAlg) doesn't support right preconditioning."
         Krylov.solve!(args...; M=M,
                       kwargs...)
@@ -176,7 +192,7 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
         Krylov.solve!(args...; M=M, N=N,
                       kwargs...)
     elseif cache.cacheval isa Krylov.MinresSolver
-        N != I  &&
+        N !== I  &&
             @warn "$(alg.KrylovAlg) doesn't support right preconditioning."
         Krylov.solve!(args...; M=M,
                       kwargs...)
@@ -200,9 +216,12 @@ end
 
 function IterativeSolversJL(args...;
                             generate_iterator = IterativeSolvers.gmres_iterable!,
-                            Pl=IterativeSolvers.Identity(),
-                            Pr=IterativeSolvers.Identity(),
+                            Pl=nothing, Pr=nothing,
                             gmres_restart=0, kwargs...)
+
+    Pl = (Pl === nothing) ? Identity() : Pl
+    Pr = (Pr === nothing) ? Identity() : Pr
+
     return IterativeSolversJL(generate_iterator, Pl, Pr, gmres_restart,
                               args, kwargs)
 end
@@ -227,16 +246,8 @@ IterativeSolversJL_MINRES(args...;kwargs...) =
 function init_cacheval(alg::IterativeSolversJL, cache::LinearCache)
     @unpack A, b, u = cache
 
-    Pl = IterativeSolvers.Identity()
-    Pr = IterativeSolvers.Identity()
-
-    if (cache.Pl != I) | (alg.Pl != IterativeSolvers.Identity())
-        Pl = ComposePreconditioner(alg.Pl, cache.Pl)
-    end
-
-    if (cache.Pr != I) | (alg.Pr != IterativeSolvers.Identity())
-        Pr = ComposePreconditioner(alg.Pr, cache.Pr)
-    end
+    Pl = get_preconditioner(alg.Pl, cache.Pl)
+    Pr = get_preconditioner(alg.Pr, cache.Pr)
 
     abstol  = cache.abstol
     reltol  = cache.reltol
@@ -249,15 +260,15 @@ function init_cacheval(alg::IterativeSolversJL, cache::LinearCache)
               alg.kwargs...)
 
     iterable = if alg.generate_iterator === IterativeSolvers.cg_iterator!
-        Pr != IterativeSolvers.Identity() &&
+        Pr !== Identity() &&
           @warn "$(alg.generate_iterator) doesn't support right preconditioning"
         alg.generate_iterator(u, A, b, Pl;
                               kwargs...)
     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.bicgstabl_iterator!
-        Pr != IterativeSolvers.Identity() &&
+        Pr !== Identity() &&
           @warn "$(alg.generate_iterator) doesn't support right preconditioning"
         alg.generate_iterator(u, A, b, alg.args...; Pl=Pl,
                               abstol=abstol, reltol=reltol,

test/runtests.jl

@@ -133,8 +133,7 @@ end
         x = rand(n,n)
         y = rand(n,n)
 
-        Pl = LinearSolve.scaling_preconditioner(s, true)
-        Pr = LinearSolve.scaling_preconditioner(s, false)
+        Pl, Pr = LinearSolve.scaling_preconditioner(s)
 
         mul!(y, Pl, x); @test y ≈ s * x
         mul!(y, Pr, x); @test y ≈ s \ x
@@ -154,13 +153,13 @@ end
         x = rand(n,n)
         y = rand(n,n)
 
-        P1 = LinearSolve.scaling_preconditioner(s1, true)
-        P2 = LinearSolve.scaling_preconditioner(s2, true)
+        P1, _ = LinearSolve.scaling_preconditioner(s1)
+        P2, _ = LinearSolve.scaling_preconditioner(s2)
 
         P  = LinearSolve.ComposePreconditioner(P1,P2)
         Pi = LinearSolve.InvComposePreconditioner(P)
 
-        @test Pi  == LinearSolve.InvComposePreconditioner(P1,P2)
+        @test Pi  == LinearSolve.InvComposePreconditioner(P1, P2)
         @test Pi  == inv(P)
         @test P   == inv(Pi)
         @test Pi' == inv(P')