From 3d82d716d5eda988a0903d7ae2984779c76f72d3 Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Thu, 25 Nov 2021 11:02:06 -0500
Subject: [PATCH 01/17] initial commit

---
 Project.toml       |  1 +
 src/LinearSolve.jl |  1 +
 src/common.jl      |  6 ++++++
 src/wrappers.jl    | 46 +++++++++++++++++++++++++++++++++++-----------
 4 files changed, 43 insertions(+), 11 deletions(-)

diff --git a/Project.toml b/Project.toml
index ca80d41a4..f6d4107a6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.1.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl
index 3e05fa6a7..16c510a78 100644
--- a/src/LinearSolve.jl
+++ b/src/LinearSolve.jl
@@ -8,6 +8,7 @@ using SparseArrays
 using SciMLBase: AbstractDiffEqOperator, AbstractLinearAlgorithm
 using Setfield
 using UnPack
+using FastBroadcast
 
 # wrap
 import Krylov
diff --git a/src/common.jl b/src/common.jl
index 3ffb3a327..ff63bbc96 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -44,6 +44,12 @@ function set_cacheval(cache::LinearCache, alg_cache)
     return cache
 end
 
+function set_prec(cache, Pl, Pr)
+    @set! cache.Pl = Pl
+    @set! cache.Pr = Pr
+    return cache
+end
+
 init_cacheval(alg::Union{SciMLLinearSolveAlgorithm,Nothing}, A, b, u) = nothing
 
 SciMLBase.init(prob::LinearProblem, args...; kwargs...) = SciMLBase.init(prob,nothing,args...;kwargs...)
diff --git a/src/wrappers.jl b/src/wrappers.jl
index ad24d81d2..4654ae9e2 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -1,7 +1,4 @@
 
-#TODO: composed preconditioners, preconditioner setter for cache,
-#   detailed tests for wrappers
-
 ## Preconditioners
 
 struct ScaleVector{T}
@@ -9,10 +6,24 @@ struct ScaleVector{T}
     isleft::Bool
 end
 
-function LinearAlgebra.ldiv!(v::ScaleVector, x)
+function LinearAlgebra.ldiv!(P::ScaleVector, x)
+    P.s == one(eltype(P.s)) && return x
+
+    if P.isleft
+        @. x = x * P.s # @.. doesnt speed up scalar computation
+    else
+        @. x = x / P.s
+    end
 end
 
-function LinearAlgebra.ldiv!(y, v::ScaleVector, x)
+function LinearAlgebra.ldiv!(y, P::ScaleVector, x)
+    P.s == one(eltype(P.s)) && return y = x
+
+    if P.isleft
+        @. y = x / P.s
+    else
+        @. y = x * P.s
+    end
 end
 
 struct ComposePreconditioner{Ti,To}
@@ -21,12 +32,19 @@ struct ComposePreconditioner{Ti,To}
     isleft::Bool
 end
 
-function LinearAlgebra.ldiv!(v::ComposePreconditioner, x)
-    @unpack inner, outer, isleft = v
+function LinearAlgebra.ldiv!(P::ComposePreconditioner, x)
+    @unpack inner, outer, isleft = P
+    if isleft
+        ldiv!(outer, x)
+        ldiv!(inner, x)
+    else
+        ldiv!(inner, x)
+        ldiv!(outer, x)
+    end
 end
 
-function LinearAlgebra.ldiv!(y, v::ComposePreconditioner, x)
-    @unpack inner, outer, isleft = v
+function LinearAlgebra.ldiv!(y, P::ComposePreconditioner, x)
+    @unpack inner, outer, isleft = P
 end
 
 ## Krylov.jl
@@ -132,6 +150,9 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
         cache = set_cacheval(cache, solver)
     end
 
+#   Pl = ComposePreconditioner(alg.Pl, cache.Pl, true)
+#   Pr = ComposePreconditioner(alg.Pr, cache.Pr, false)
+
     atol    = cache.abstol
     rtol    = cache.reltol
     itmax   = cache.maxiters
@@ -204,8 +225,11 @@ IterativeSolversJL_MINRES(args...;kwargs...) =
 function init_cacheval(alg::IterativeSolversJL, cache::LinearCache)
     @unpack A, b, u = cache
 
-    Pl = (alg.Pl == LinearAlgebra.I) ? IterativeSolvers.Identity() : alg.Pl
-    Pr = (alg.Pr == LinearAlgebra.I) ? IterativeSolvers.Identity() : alg.Pr
+    Pl = alg.Pl
+    Pr = alg.Pr
+
+#   Pl = ComposePreconditioner(alg.Pl, cache.Pl, true)
+#   Pr = ComposePreconditioner(alg.Pr, cache.Pr, false)
 
     abstol  = cache.abstol
     reltol  = cache.reltol

From 6b6d2dd07a899b4f22d028ecebc116699d4cd959 Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Thu, 25 Nov 2021 13:16:15 -0500
Subject: [PATCH 02/17] commiting updates - still WIP

---
 src/LinearSolve.jl |  2 +-
 src/common.jl      |  4 +--
 src/wrappers.jl    | 86 +++++++++++++++++++++++++++++-----------------
 3 files changed, 57 insertions(+), 35 deletions(-)

diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl
index 16c510a78..2e5b0a21b 100644
--- a/src/LinearSolve.jl
+++ b/src/LinearSolve.jl
@@ -2,7 +2,7 @@ module LinearSolve
 
 using ArrayInterface
 using RecursiveFactorization
-using Base: cache_dependencies, Bool
+using Base: cache_dependencies, Bool, eltype
 using LinearAlgebra
 using SparseArrays
 using SciMLBase: AbstractDiffEqOperator, AbstractLinearAlgorithm
diff --git a/src/common.jl b/src/common.jl
index ff63bbc96..402a2c678 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -70,8 +70,8 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
     isfresh = cacheval === nothing
     Tc = isfresh ? Any : typeof(cacheval)
 
-    Pl = LinearAlgebra.I
-    Pr = LinearAlgebra.I
+    Pl = ScaleVector(one(eltype(A)), true)
+    Pr = ScaleVector(one(eltype(A)), false)
 
     A = alias_A ? A : deepcopy(A)
     b = alias_b ? b : deepcopy(b)
diff --git a/src/wrappers.jl b/src/wrappers.jl
index 4654ae9e2..48ec8274f 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -1,50 +1,72 @@
 
 ## Preconditioners
 
+"""
+ P  *  x = x .* (1/P.s)
+
+ Pi *  x = x .* (P.s)
+"""
 struct ScaleVector{T}
     s::T
     isleft::Bool
 end
 
-function LinearAlgebra.ldiv!(P::ScaleVector, x)
-    P.s == one(eltype(P.s)) && return x
+function LinearAlgebra.mul!(y, A::ScaleVector, x)
+end
+
+function LinearAlgebra.mul!(C, A::ScaleVector, B, α, β)
+end
+
+function LinearAlgebra.ldiv!(A::ScaleVector, x)
+    A.s == one(eltype(A.s)) && return x
 
-    if P.isleft
-        @. x = x * P.s # @.. doesnt speed up scalar computation
+    if A.isleft
+        @. x = x * A.s
     else
-        @. x = x / P.s
+        @. x = x / A.s
     end
 end
 
-function LinearAlgebra.ldiv!(y, P::ScaleVector, x)
-    P.s == one(eltype(P.s)) && return y = x
+function LinearAlgebra.ldiv!(y, A::ScaleVector, x)
+    P.s == one(eltype(A.s)) && return y = x
 
-    if P.isleft
-        @. y = x / P.s
+    if A.isleft
+        @. y = x * A.s
     else
-        @. y = x * P.s
+        @. y = x / A.s
     end
 end
 
+Base.eltype(A::ScaleVector) = eltype(A.s)
+
+"""
+ C  * x = P  * Q  * x
+
+ Ci * x = Qi * Pi * x
+"""
 struct ComposePreconditioner{Ti,To}
     inner::Ti
     outer::To
-    isleft::Bool
 end
 
-function LinearAlgebra.ldiv!(P::ComposePreconditioner, x)
-    @unpack inner, outer, isleft = P
-    if isleft
-        ldiv!(outer, x)
-        ldiv!(inner, x)
-    else
-        ldiv!(inner, x)
-        ldiv!(outer, x)
-    end
+function LinearAlgebra.ldiv!(A::ComposePreconditioner, x)
+    @unpack inner, outer, isleft = A
+
+    ldiv!(inner, x)
+    ldiv!(outer, x)
+end
+
+function LinearAlgebra.ldiv!(y, A::ComposePreconditioner, x)
+    @unpack inner, outer = A
+
+    ldiv!(y, inner, x)
+    ldiv!(outer, y)
+end
+
+function LinearAlgebra.mul!(y, A::ComposePreconditioner, x)
 end
 
-function LinearAlgebra.ldiv!(y, P::ComposePreconditioner, x)
-    @unpack inner, outer, isleft = P
+function LinearAlgebra.mul!(C, A::ComposePreconditioner, B, α, β)
 end
 
 ## Krylov.jl
@@ -150,8 +172,8 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
         cache = set_cacheval(cache, solver)
     end
 
-#   Pl = ComposePreconditioner(alg.Pl, cache.Pl, true)
-#   Pr = ComposePreconditioner(alg.Pr, cache.Pr, false)
+    M = alg.Pl #ComposePreconditioner(alg.Pl, cache.Pl) # left precond
+    N = alg.Pr #ComposePreconditioner(alg.Pr, cache.Pr) # right 
 
     atol    = cache.abstol
     rtol    = cache.reltol
@@ -163,20 +185,20 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
               alg.kwargs...)
 
     if cache.cacheval isa Krylov.CgSolver
-        alg.Pr != LinearAlgebra.I  &&
+        N != LinearAlgebra.I  &&
             @warn "$(alg.KrylovAlg) doesn't support right preconditioning."
-        Krylov.solve!(args...; M=alg.Pl,
+        Krylov.solve!(args...; M=M,
                       kwargs...)
     elseif cache.cacheval isa Krylov.GmresSolver
-        Krylov.solve!(args...; M=alg.Pl, N=alg.Pr,
+        Krylov.solve!(args...; M=M, N=N,
                       kwargs...)
     elseif cache.cacheval isa Krylov.BicgstabSolver
-        Krylov.solve!(args...; M=alg.Pl, N=alg.Pr,
+        Krylov.solve!(args...; M=M, N=N,
                       kwargs...)
     elseif cache.cacheval isa Krylov.MinresSolver
-        alg.Pr != LinearAlgebra.I  &&
+        N != LinearAlgebra.I  &&
             @warn "$(alg.KrylovAlg) doesn't support right preconditioning."
-        Krylov.solve!(args...; M=alg.Pl,
+        Krylov.solve!(args...; M=M,
                       kwargs...)
     else
         Krylov.solve!(args...; kwargs...)
@@ -228,8 +250,8 @@ function init_cacheval(alg::IterativeSolversJL, cache::LinearCache)
     Pl = alg.Pl
     Pr = alg.Pr
 
-#   Pl = ComposePreconditioner(alg.Pl, cache.Pl, true)
-#   Pr = ComposePreconditioner(alg.Pr, cache.Pr, false)
+#   Pl = ComposePreconditioner(alg.Pl, cache.Pl)
+#   Pr = ComposePreconditioner(alg.Pr, cache.Pr)
 
     abstol  = cache.abstol
     reltol  = cache.reltol

From 73d8808747747f0c5339f698a279edd198438948 Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Thu, 25 Nov 2021 15:55:00 -0500
Subject: [PATCH 03/17] need mul! for krylov.jl

---
 src/wrappers.jl | 32 +++++++++++++++++++++++++-------
 1 file changed, 25 insertions(+), 7 deletions(-)

diff --git a/src/wrappers.jl b/src/wrappers.jl
index 48ec8274f..d25978825 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -11,10 +11,22 @@ struct ScaleVector{T}
     isleft::Bool
 end
 
+# y = A x
 function LinearAlgebra.mul!(y, A::ScaleVector, x)
+    A.s == one(eltype(A.s)) && return y = x
+
+    if A.isleft
+        @. x = x / A.s
+    else
+        @. x = x * A.s
+    end
 end
 
+# A B α + C β
 function LinearAlgebra.mul!(C, A::ScaleVector, B, α, β)
+    
+    tmp = zero(B)
+    C = β * C + α * mul!(tmp, A, B)
 end
 
 function LinearAlgebra.ldiv!(A::ScaleVector, x)
@@ -49,8 +61,20 @@ struct ComposePreconditioner{Ti,To}
     outer::To
 end
 
+# y = A x
+function LinearAlgebra.mul!(y, A::ComposePreconditioner, x)
+    @unpack inner, outer = A
+    mul!(y, inner, x)
+    y = outer * y
+end
+
+# A B α + C β
+function LinearAlgebra.mul!(C, A::ComposePreconditioner, B, α, β)
+    @unpack inner, outer = A
+end
+
 function LinearAlgebra.ldiv!(A::ComposePreconditioner, x)
-    @unpack inner, outer, isleft = A
+    @unpack inner, outer = A
 
     ldiv!(inner, x)
     ldiv!(outer, x)
@@ -63,12 +87,6 @@ function LinearAlgebra.ldiv!(y, A::ComposePreconditioner, x)
     ldiv!(outer, y)
 end
 
-function LinearAlgebra.mul!(y, A::ComposePreconditioner, x)
-end
-
-function LinearAlgebra.mul!(C, A::ComposePreconditioner, B, α, β)
-end
-
 ## Krylov.jl
 
 struct KrylovJL{F,Tl,Tr,I,A,K} <: AbstractKrylovSubspaceMethod

From 7cc8f3b876a36a750a7592e03dc2563fdffab625 Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Fri, 26 Nov 2021 16:04:15 -0500
Subject: [PATCH 04/17] commiting updates

---
 src/LinearSolve.jl |  3 +-
 src/wrappers.jl    | 71 +++++++++++++++++++++++++---------------------
 test/runtests.jl   | 52 +++++++++++++++++++++++++++++++++
 3 files changed, 93 insertions(+), 33 deletions(-)

diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl
index 2e5b0a21b..26d547538 100644
--- a/src/LinearSolve.jl
+++ b/src/LinearSolve.jl
@@ -2,7 +2,8 @@ module LinearSolve
 
 using ArrayInterface
 using RecursiveFactorization
-using Base: cache_dependencies, Bool, eltype
+using Base: cache_dependencies, Bool
+import Base: eltype, *
 using LinearAlgebra
 using SparseArrays
 using SciMLBase: AbstractDiffEqOperator, AbstractLinearAlgorithm
diff --git a/src/wrappers.jl b/src/wrappers.jl
index d25978825..a9d8d3701 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -2,58 +2,59 @@
 ## Preconditioners
 
 """
+ LEFT 
  P  *  x = x .* (1/P.s)
-
  Pi *  x = x .* (P.s)
+
+ Right
+ P  *  x = x .* (P.s)
+ Pi *  x = x .* (1/P.s)
 """
 struct ScaleVector{T}
     s::T
     isleft::Bool
 end
 
+Base.eltype(A::ScaleVector) = eltype(A.s)
+
+#function Base.*(A::ScaleVector, x)
+#    y = similar(x)
+#    mul!(y, A, x)
+#end
+
 # y = A x
 function LinearAlgebra.mul!(y, A::ScaleVector, x)
     A.s == one(eltype(A.s)) && return y = x
 
-    if A.isleft
-        @. x = x / A.s
-    else
-        @. x = x * A.s
-    end
+    s = A.isleft ? 1/A.s : A.s
+    mul!(y, s, x)
+
 end
 
 # A B α + C β
 function LinearAlgebra.mul!(C, A::ScaleVector, B, α, β)
-    
-    tmp = zero(B)
-    C = β * C + α * mul!(tmp, A, B)
+    A.s == one(eltype(A.s)) && return @. C = α * B + β * C
+
+    s = A.isleft ? 1/A.s : A.s
+    mul!(C, s, B, α, β)
 end
 
 function LinearAlgebra.ldiv!(A::ScaleVector, x)
     A.s == one(eltype(A.s)) && return x
 
-    if A.isleft
-        @. x = x * A.s
-    else
-        @. x = x / A.s
-    end
+    s = A.isleft ? A.s : 1/A.s
+    @. x = x * s
 end
 
 function LinearAlgebra.ldiv!(y, A::ScaleVector, x)
-    P.s == one(eltype(A.s)) && return y = x
+    A.s == one(eltype(A.s)) && return y = x
 
-    if A.isleft
-        @. y = x * A.s
-    else
-        @. y = x / A.s
-    end
+    s = A.isleft ? A.s : 1/A.s
+    mul!(y, s, x)
 end
 
-Base.eltype(A::ScaleVector) = eltype(A.s)
-
 """
  C  * x = P  * Q  * x
-
  Ci * x = Qi * Pi * x
 """
 struct ComposePreconditioner{Ti,To}
@@ -61,16 +62,22 @@ struct ComposePreconditioner{Ti,To}
     outer::To
 end
 
+Base.eltype(A::ComposePreconditioner) = Float64 #eltype(A.inner)
+
 # y = A x
 function LinearAlgebra.mul!(y, A::ComposePreconditioner, x)
     @unpack inner, outer = A
-    mul!(y, inner, x)
-    y = outer * y
+    tmp = similar(y)
+    mul!(tmp, outer, x)
+    mul!(y, inner, tmp)
 end
 
 # A B α + C β
 function LinearAlgebra.mul!(C, A::ComposePreconditioner, B, α, β)
     @unpack inner, outer = A
+    tmp = similar(B)
+    mul!(tmp, inner, B)
+    mul!(C, outer, tmp, α, β)
 end
 
 function LinearAlgebra.ldiv!(A::ComposePreconditioner, x)
@@ -190,8 +197,11 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
         cache = set_cacheval(cache, solver)
     end
 
-    M = alg.Pl #ComposePreconditioner(alg.Pl, cache.Pl) # left precond
-    N = alg.Pr #ComposePreconditioner(alg.Pr, cache.Pr) # right 
+    M = alg.Pl
+    N = alg.Pr
+
+#   M = ComposePreconditioner(alg.Pl, cache.Pl) # left precond
+#   N = ComposePreconditioner(alg.Pr, cache.Pr) # right 
 
     atol    = cache.abstol
     rtol    = cache.reltol
@@ -265,11 +275,8 @@ IterativeSolversJL_MINRES(args...;kwargs...) =
 function init_cacheval(alg::IterativeSolversJL, cache::LinearCache)
     @unpack A, b, u = cache
 
-    Pl = alg.Pl
-    Pr = alg.Pr
-
-#   Pl = ComposePreconditioner(alg.Pl, cache.Pl)
-#   Pr = ComposePreconditioner(alg.Pr, cache.Pr)
+    Pl = ComposePreconditioner(alg.Pl, cache.Pl)
+    Pr = ComposePreconditioner(alg.Pr, cache.Pr)
 
     abstol  = cache.abstol
     reltol  = cache.reltol
diff --git a/test/runtests.jl b/test/runtests.jl
index dd2b99f4b..efee8f8c6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -34,6 +34,8 @@ function test_interface(alg, prob1, prob2)
     return
 end
 
+@testset "LinearSolve" begin
+
 @testset "Default Linear Solver" begin
     test_interface(nothing, prob1, prob2)
 
@@ -123,3 +125,53 @@ end
         end
     end
 end
+
+@testset "Preconditioners" begin
+    @testset "ScaleVector" begin
+        s = rand()
+        α = rand()
+        β = rand()
+
+        x = rand(n,n)
+        y = rand(n,n)
+
+        Pl = LinearSolve.ScaleVector(s, true)
+        Pr = LinearSolve.ScaleVector(s, false)
+
+        mul!(y, Pl, x)
+        mul!(y, Pr, x)
+
+        mul!(y, Pl, x, α, β)
+        mul!(y, Pr, x, α, β)
+
+        ldiv!(Pl, x)
+        ldiv!(Pr, x)
+
+        ldiv!(y, Pl, x)
+        ldiv!(y, Pr, x)
+
+    end
+
+    @testset "ComposePreconditioenr" begin
+        s = rand()
+        α = rand()
+        β = rand()
+
+        x = rand(n,n)
+        y = rand(n,n)
+
+        Pi = LinearSolve.ScaleVector(s, true)
+        Po = LinearSolve.ScaleVector(s, false)
+
+        P = LinearSolve.ComposePreconditioner(Pi,Po)
+
+        mul!(y, P, x)
+        mul!(y, P, x, α, β)
+
+        ldiv!(P, x)
+        ldiv!(y, P, x)
+
+    end
+end
+
+end # testset

From e020ce06c00a7483f7498ffd9902005dfd773940 Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Fri, 26 Nov 2021 16:58:22 -0500
Subject: [PATCH 05/17] rm fastbroadcast

---
 Project.toml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index f6d4107a6..ca80d41a4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,7 +5,6 @@ version = "0.1.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"

From 3575cb22bd277569045f795035081e9b84f5fada Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Fri, 26 Nov 2021 16:59:48 -0500
Subject: [PATCH 06/17] implemented adjoint, inv

---
 src/LinearSolve.jl |  3 +--
 src/wrappers.jl    | 16 +++++++++-------
 2 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl
index 26d547538..9ce4d735a 100644
--- a/src/LinearSolve.jl
+++ b/src/LinearSolve.jl
@@ -3,13 +3,12 @@ module LinearSolve
 using ArrayInterface
 using RecursiveFactorization
 using Base: cache_dependencies, Bool
-import Base: eltype, *
+import Base: eltype, adjoint, inv
 using LinearAlgebra
 using SparseArrays
 using SciMLBase: AbstractDiffEqOperator, AbstractLinearAlgorithm
 using Setfield
 using UnPack
-using FastBroadcast
 
 # wrap
 import Krylov
diff --git a/src/wrappers.jl b/src/wrappers.jl
index a9d8d3701..9f7f9b82c 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -6,7 +6,7 @@
  P  *  x = x .* (1/P.s)
  Pi *  x = x .* (P.s)
 
- Right
+ RIGHT
  P  *  x = x .* (P.s)
  Pi *  x = x .* (1/P.s)
 """
@@ -16,11 +16,8 @@ struct ScaleVector{T}
 end
 
 Base.eltype(A::ScaleVector) = eltype(A.s)
-
-#function Base.*(A::ScaleVector, x)
-#    y = similar(x)
-#    mul!(y, A, x)
-#end
+Base.adjoint(A::ScaleVector) = copy(A)
+Base.inv(A::ScaleVector) = ScaleVector(1/A.s, A.isleft)
 
 # y = A x
 function LinearAlgebra.mul!(y, A::ScaleVector, x)
@@ -28,7 +25,6 @@ function LinearAlgebra.mul!(y, A::ScaleVector, x)
 
     s = A.isleft ? 1/A.s : A.s
     mul!(y, s, x)
-
 end
 
 # A B α + C β
@@ -63,6 +59,8 @@ struct ComposePreconditioner{Ti,To}
 end
 
 Base.eltype(A::ComposePreconditioner) = Float64 #eltype(A.inner)
+Base.adjoint(A::ComposePreconditioner) = ComposePreconditioner(A.outer', A.inner')
+Base.inv(A::ComposePreconditioner) = ComposePreconditioner(inv(A.outer), inv(A.inner))
 
 # y = A x
 function LinearAlgebra.mul!(y, A::ComposePreconditioner, x)
@@ -200,6 +198,10 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
     M = alg.Pl
     N = alg.Pr
 
+    """
+    TODO - pass in inv(Pl), inv(Pr) to Krylov.jl
+    """
+
 #   M = ComposePreconditioner(alg.Pl, cache.Pl) # left precond
 #   N = ComposePreconditioner(alg.Pr, cache.Pr) # right 
 

From 1793e51e57af527db392035880a913887235d01c Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Fri, 26 Nov 2021 17:42:34 -0500
Subject: [PATCH 07/17] scalevector --> default_preconditioner, start with lazy
 invcomposepreconditioner

---
 src/common.jl    |  4 +--
 src/wrappers.jl  | 84 ++++++++++++++++++------------------------------
 test/runtests.jl |  8 ++---
 3 files changed, 38 insertions(+), 58 deletions(-)

diff --git a/src/common.jl b/src/common.jl
index 402a2c678..ba789c34f 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -70,8 +70,8 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
     isfresh = cacheval === nothing
     Tc = isfresh ? Any : typeof(cacheval)
 
-    Pl = ScaleVector(one(eltype(A)), true)
-    Pr = ScaleVector(one(eltype(A)), false)
+    Pl = default_preconditioner(one(eltype(A)), true)
+    Pr = default_preconditioner(one(eltype(A)), false)
 
     A = alias_A ? A : deepcopy(A)
     b = alias_b ? b : deepcopy(b)
diff --git a/src/wrappers.jl b/src/wrappers.jl
index 9f7f9b82c..94837afe3 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -1,53 +1,7 @@
 
 ## Preconditioners
 
-"""
- LEFT 
- P  *  x = x .* (1/P.s)
- Pi *  x = x .* (P.s)
-
- RIGHT
- P  *  x = x .* (P.s)
- Pi *  x = x .* (1/P.s)
-"""
-struct ScaleVector{T}
-    s::T
-    isleft::Bool
-end
-
-Base.eltype(A::ScaleVector) = eltype(A.s)
-Base.adjoint(A::ScaleVector) = copy(A)
-Base.inv(A::ScaleVector) = ScaleVector(1/A.s, A.isleft)
-
-# y = A x
-function LinearAlgebra.mul!(y, A::ScaleVector, x)
-    A.s == one(eltype(A.s)) && return y = x
-
-    s = A.isleft ? 1/A.s : A.s
-    mul!(y, s, x)
-end
-
-# A B α + C β
-function LinearAlgebra.mul!(C, A::ScaleVector, B, α, β)
-    A.s == one(eltype(A.s)) && return @. C = α * B + β * C
-
-    s = A.isleft ? 1/A.s : A.s
-    mul!(C, s, B, α, β)
-end
-
-function LinearAlgebra.ldiv!(A::ScaleVector, x)
-    A.s == one(eltype(A.s)) && return x
-
-    s = A.isleft ? A.s : 1/A.s
-    @. x = x * s
-end
-
-function LinearAlgebra.ldiv!(y, A::ScaleVector, x)
-    A.s == one(eltype(A.s)) && return y = x
-
-    s = A.isleft ? A.s : 1/A.s
-    mul!(y, s, x)
-end
+default_preconditioner(s, isleft) = isleft ? I * s : I * (1/s)
 
 """
  C  * x = P  * Q  * x
@@ -60,9 +14,8 @@ end
 
 Base.eltype(A::ComposePreconditioner) = Float64 #eltype(A.inner)
 Base.adjoint(A::ComposePreconditioner) = ComposePreconditioner(A.outer', A.inner')
-Base.inv(A::ComposePreconditioner) = ComposePreconditioner(inv(A.outer), inv(A.inner))
+Base.inv(A::ComposePreconditioner) = InvComposePreconditioner(A)
 
-# y = A x
 function LinearAlgebra.mul!(y, A::ComposePreconditioner, x)
     @unpack inner, outer = A
     tmp = similar(y)
@@ -70,7 +23,6 @@ function LinearAlgebra.mul!(y, A::ComposePreconditioner, x)
     mul!(y, inner, tmp)
 end
 
-# A B α + C β
 function LinearAlgebra.mul!(C, A::ComposePreconditioner, B, α, β)
     @unpack inner, outer = A
     tmp = similar(B)
@@ -92,6 +44,34 @@ function LinearAlgebra.ldiv!(y, A::ComposePreconditioner, x)
     ldiv!(outer, y)
 end
 
+struct InvComposePreconditioner{Tp <: ComposePreconditioner}
+    P::Tp
+end
+
+InvComposePreconditioner(inner, outer) = InvComposePreconditioner(ComposePreconditioner(inner, outer))
+
+Base.eltype(A::InvComposePreconditioner) = Float64 #eltype(A.inner)
+#Base.adjoint(A::InvComposePreconditioner) = ComposePreconditioner(A.outer', A.inner')
+Base.inv(A::InvComposePreconditioner) = ComposePreconditioner(A.P)
+
+function LinearAlgebra.mul!(y, A::InvComposePreconditioner, x)
+    @unpack P = A
+    ldiv!(y, P, x)
+end
+
+function LinearAlgebra.mul!(C, A::InvComposePreconditioner, B, α, β)
+    @unpack P = A
+end
+
+function LinearAlgebra.ldiv!(A::InvComposePreconditioner, x)
+    @unpack P = A
+end
+
+function LinearAlgebra.ldiv!(y, A::InvComposePreconditioner, x)
+    @unpack P = A
+    mul!(y, P, x)
+end
+
 ## Krylov.jl
 
 struct KrylovJL{F,Tl,Tr,I,A,K} <: AbstractKrylovSubspaceMethod
@@ -202,8 +182,8 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
     TODO - pass in inv(Pl), inv(Pr) to Krylov.jl
     """
 
-#   M = ComposePreconditioner(alg.Pl, cache.Pl) # left precond
-#   N = ComposePreconditioner(alg.Pr, cache.Pr) # right 
+#   M = InvComposePreconditioner(alg.Pl, cache.Pl) # left precond
+#   N = InvComposePreconditioner(alg.Pr, cache.Pr) # right 
 
     atol    = cache.abstol
     rtol    = cache.reltol
diff --git a/test/runtests.jl b/test/runtests.jl
index efee8f8c6..a6882e411 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -135,8 +135,8 @@ end
         x = rand(n,n)
         y = rand(n,n)
 
-        Pl = LinearSolve.ScaleVector(s, true)
-        Pr = LinearSolve.ScaleVector(s, false)
+        Pl = LinearSolve.default_preconditioner(s, true)
+        Pr = LinearSolve.default_preconditioner(s, false)
 
         mul!(y, Pl, x)
         mul!(y, Pr, x)
@@ -160,8 +160,8 @@ end
         x = rand(n,n)
         y = rand(n,n)
 
-        Pi = LinearSolve.ScaleVector(s, true)
-        Po = LinearSolve.ScaleVector(s, false)
+        Pi = LinearSolve.default_preconditioner(s, true)
+        Po = LinearSolve.default_preconditioner(s, false)
 
         P = LinearSolve.ComposePreconditioner(Pi,Po)
 

From 27a7d88ffd1bca29b4ca8d9f1765d4c40180f0c2 Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Fri, 26 Nov 2021 22:37:57 -0500
Subject: [PATCH 08/17] InvComposePreconditioner done

---
 src/wrappers.jl  | 22 ++++++++++------------
 test/runtests.jl | 32 ++++++++++++++++++++++++++++++++
 2 files changed, 42 insertions(+), 12 deletions(-)

diff --git a/src/wrappers.jl b/src/wrappers.jl
index 94837afe3..f4a7eda40 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -50,9 +50,9 @@ end
 
 InvComposePreconditioner(inner, outer) = InvComposePreconditioner(ComposePreconditioner(inner, outer))
 
-Base.eltype(A::InvComposePreconditioner) = Float64 #eltype(A.inner)
-#Base.adjoint(A::InvComposePreconditioner) = ComposePreconditioner(A.outer', A.inner')
-Base.inv(A::InvComposePreconditioner) = ComposePreconditioner(A.P)
+Base.eltype(A::InvComposePreconditioner) = Base.eltype(A.P)
+Base.adjoint(A::InvComposePreconditioner) = InvComposePreconditioner(A.P')
+Base.inv(A::InvComposePreconditioner) = deepcopy(A.P)
 
 function LinearAlgebra.mul!(y, A::InvComposePreconditioner, x)
     @unpack P = A
@@ -61,10 +61,15 @@ end
 
 function LinearAlgebra.mul!(C, A::InvComposePreconditioner, B, α, β)
     @unpack P = A
+    tmp = copy(B)
+    ldiv!(tmp, P, B)
+    mul!(C, I, tmp, α, β)
 end
 
 function LinearAlgebra.ldiv!(A::InvComposePreconditioner, x)
     @unpack P = A
+    y = copy(x)
+    mul!(x, P, y)
 end
 
 function LinearAlgebra.ldiv!(y, A::InvComposePreconditioner, x)
@@ -175,15 +180,8 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
         cache = set_cacheval(cache, solver)
     end
 
-    M = alg.Pl
-    N = alg.Pr
-
-    """
-    TODO - pass in inv(Pl), inv(Pr) to Krylov.jl
-    """
-
-#   M = InvComposePreconditioner(alg.Pl, cache.Pl) # left precond
-#   N = InvComposePreconditioner(alg.Pr, cache.Pr) # right 
+    M = InvComposePreconditioner(alg.Pl, cache.Pl) # left precond
+    N = InvComposePreconditioner(alg.Pr, cache.Pr) # right 
 
     atol    = cache.abstol
     rtol    = cache.reltol
diff --git a/test/runtests.jl b/test/runtests.jl
index a6882e411..e711e70b3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -172,6 +172,38 @@ end
         ldiv!(y, P, x)
 
     end
+
+    @testset "InvComposePreconditioenr" begin
+        s = rand()
+        α = rand()
+        β = rand()
+
+        x = rand(n,n)
+        y = rand(n,n)
+
+        P1 = LinearSolve.default_preconditioner(s, true)
+        P2 = LinearSolve.default_preconditioner(s, false)
+
+        P  = LinearSolve.ComposePreconditioner(P1,P2)
+        Pi = LinearSolve.InvComposePreconditioner(P)
+
+        @test Pi == LinearSolve.InvComposePreconditioner(P1,P2)
+        @test Pi == inv(P)
+        @test P  == inv(Pi)
+
+        mul!(y, P, x)
+        mul!(y, P, x, α, β)
+
+        ldiv!(P, x)
+        ldiv!(y, P, x)
+
+        mul!(y, Pi, x)
+        mul!(y, Pi, x, α, β)
+
+        ldiv!(Pi, x)
+        ldiv!(y, Pi, x)
+
+    end
 end
 
 end # testset

From 72f4445dbce84edebd8223b41e7c2465cd8b0ac2 Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Sat, 27 Nov 2021 07:52:28 -0500
Subject: [PATCH 09/17] comment out unused mul, ldiv

---
 src/wrappers.jl  |  4 ++++
 test/runtests.jl | 33 +++++++--------------------------
 2 files changed, 11 insertions(+), 26 deletions(-)

diff --git a/src/wrappers.jl b/src/wrappers.jl
index f4a7eda40..f2e64180b 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -16,6 +16,7 @@ Base.eltype(A::ComposePreconditioner) = Float64 #eltype(A.inner)
 Base.adjoint(A::ComposePreconditioner) = ComposePreconditioner(A.outer', A.inner')
 Base.inv(A::ComposePreconditioner) = InvComposePreconditioner(A)
 
+#= unused
 function LinearAlgebra.mul!(y, A::ComposePreconditioner, x)
     @unpack inner, outer = A
     tmp = similar(y)
@@ -29,6 +30,7 @@ function LinearAlgebra.mul!(C, A::ComposePreconditioner, B, α, β)
     mul!(tmp, inner, B)
     mul!(C, outer, tmp, α, β)
 end
+=#
 
 function LinearAlgebra.ldiv!(A::ComposePreconditioner, x)
     @unpack inner, outer = A
@@ -59,6 +61,7 @@ function LinearAlgebra.mul!(y, A::InvComposePreconditioner, x)
     ldiv!(y, P, x)
 end
 
+#=
 function LinearAlgebra.mul!(C, A::InvComposePreconditioner, B, α, β)
     @unpack P = A
     tmp = copy(B)
@@ -76,6 +79,7 @@ function LinearAlgebra.ldiv!(y, A::InvComposePreconditioner, x)
     @unpack P = A
     mul!(y, P, x)
 end
+=#
 
 ## Krylov.jl
 
diff --git a/test/runtests.jl b/test/runtests.jl
index e711e70b3..5772866ed 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -160,27 +160,6 @@ end
         x = rand(n,n)
         y = rand(n,n)
 
-        Pi = LinearSolve.default_preconditioner(s, true)
-        Po = LinearSolve.default_preconditioner(s, false)
-
-        P = LinearSolve.ComposePreconditioner(Pi,Po)
-
-        mul!(y, P, x)
-        mul!(y, P, x, α, β)
-
-        ldiv!(P, x)
-        ldiv!(y, P, x)
-
-    end
-
-    @testset "InvComposePreconditioenr" begin
-        s = rand()
-        α = rand()
-        β = rand()
-
-        x = rand(n,n)
-        y = rand(n,n)
-
         P1 = LinearSolve.default_preconditioner(s, true)
         P2 = LinearSolve.default_preconditioner(s, false)
 
@@ -191,17 +170,19 @@ end
         @test Pi == inv(P)
         @test P  == inv(Pi)
 
-        mul!(y, P, x)
-        mul!(y, P, x, α, β)
+        # ComposePreconditioner
+#       mul!(y, P, x)
+#       mul!(y, P, x, α, β)
 
         ldiv!(P, x)
         ldiv!(y, P, x)
 
+        # InvComposePreconditioner
         mul!(y, Pi, x)
-        mul!(y, Pi, x, α, β)
+#       mul!(y, Pi, x, α, β)
 
-        ldiv!(Pi, x)
-        ldiv!(y, Pi, x)
+#       ldiv!(Pi, x)
+#       ldiv!(y, Pi, x)
 
     end
 end

From c5bd4cb12a5569399ae4e8e41c08b6a3d4756113 Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Sat, 27 Nov 2021 07:58:42 -0500
Subject: [PATCH 10/17] rm unused methods

---
 src/wrappers.jl  | 36 ------------------------------------
 test/runtests.jl |  7 -------
 2 files changed, 43 deletions(-)

diff --git a/src/wrappers.jl b/src/wrappers.jl
index f2e64180b..0adf9b916 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -16,22 +16,6 @@ Base.eltype(A::ComposePreconditioner) = Float64 #eltype(A.inner)
 Base.adjoint(A::ComposePreconditioner) = ComposePreconditioner(A.outer', A.inner')
 Base.inv(A::ComposePreconditioner) = InvComposePreconditioner(A)
 
-#= unused
-function LinearAlgebra.mul!(y, A::ComposePreconditioner, x)
-    @unpack inner, outer = A
-    tmp = similar(y)
-    mul!(tmp, outer, x)
-    mul!(y, inner, tmp)
-end
-
-function LinearAlgebra.mul!(C, A::ComposePreconditioner, B, α, β)
-    @unpack inner, outer = A
-    tmp = similar(B)
-    mul!(tmp, inner, B)
-    mul!(C, outer, tmp, α, β)
-end
-=#
-
 function LinearAlgebra.ldiv!(A::ComposePreconditioner, x)
     @unpack inner, outer = A
 
@@ -61,26 +45,6 @@ function LinearAlgebra.mul!(y, A::InvComposePreconditioner, x)
     ldiv!(y, P, x)
 end
 
-#=
-function LinearAlgebra.mul!(C, A::InvComposePreconditioner, B, α, β)
-    @unpack P = A
-    tmp = copy(B)
-    ldiv!(tmp, P, B)
-    mul!(C, I, tmp, α, β)
-end
-
-function LinearAlgebra.ldiv!(A::InvComposePreconditioner, x)
-    @unpack P = A
-    y = copy(x)
-    mul!(x, P, y)
-end
-
-function LinearAlgebra.ldiv!(y, A::InvComposePreconditioner, x)
-    @unpack P = A
-    mul!(y, P, x)
-end
-=#
-
 ## Krylov.jl
 
 struct KrylovJL{F,Tl,Tr,I,A,K} <: AbstractKrylovSubspaceMethod
diff --git a/test/runtests.jl b/test/runtests.jl
index 5772866ed..3c5eb90e1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -171,18 +171,11 @@ end
         @test P  == inv(Pi)
 
         # ComposePreconditioner
-#       mul!(y, P, x)
-#       mul!(y, P, x, α, β)
-
         ldiv!(P, x)
         ldiv!(y, P, x)
 
         # InvComposePreconditioner
         mul!(y, Pi, x)
-#       mul!(y, Pi, x, α, β)
-
-#       ldiv!(Pi, x)
-#       ldiv!(y, Pi, x)
 
     end
 end

From 36e4dd784231344ecaf7624a6bb439f8ce267086 Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Mon, 29 Nov 2021 04:38:01 -0500
Subject: [PATCH 11/17] default_precond kwarg

---
 src/common.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/common.jl b/src/common.jl
index ba789c34f..c4731565d 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -60,6 +60,7 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
                         reltol=√eps(eltype(prob.A)),
                         maxiters=length(prob.b),
                         verbose=false,
+                        preconditioner_scale=one(eltype(prob.A)),
                         kwargs...,
                        )
     @unpack A, b, u0, p = prob
@@ -70,8 +71,8 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
     isfresh = cacheval === nothing
     Tc = isfresh ? Any : typeof(cacheval)
 
-    Pl = default_preconditioner(one(eltype(A)), true)
-    Pr = default_preconditioner(one(eltype(A)), false)
+    Pl = default_preconditioner(preconditioner_scale, true)
+    Pr = default_preconditioner(preconditioner_scale, false)
 
     A = alias_A ? A : deepcopy(A)
     b = alias_b ? b : deepcopy(b)

From 32aa9b458c3adfa661028bcd27f9e06475f584f0 Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Mon, 29 Nov 2021 04:38:42 -0500
Subject: [PATCH 12/17] precond handling

---
 src/wrappers.jl | 30 +++++++++++++++++++++++-------
 1 file changed, 23 insertions(+), 7 deletions(-)

diff --git a/src/wrappers.jl b/src/wrappers.jl
index 0adf9b916..f04bf90ef 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -12,7 +12,7 @@ struct ComposePreconditioner{Ti,To}
     outer::To
 end
 
-Base.eltype(A::ComposePreconditioner) = Float64 #eltype(A.inner)
+Base.eltype(A::ComposePreconditioner) = promote_type(eltype(A.inner), eltype(A.outer))
 Base.adjoint(A::ComposePreconditioner) = ComposePreconditioner(A.outer', A.inner')
 Base.inv(A::ComposePreconditioner) = InvComposePreconditioner(A)
 
@@ -148,8 +148,16 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
         cache = set_cacheval(cache, solver)
     end
 
-    M = InvComposePreconditioner(alg.Pl, cache.Pl) # left precond
-    N = InvComposePreconditioner(alg.Pr, cache.Pr) # right 
+    M = I # left precond
+    N = I # right precond
+
+    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
 
     atol    = cache.abstol
     rtol    = cache.reltol
@@ -161,7 +169,7 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
               alg.kwargs...)
 
     if cache.cacheval isa Krylov.CgSolver
-        N != LinearAlgebra.I  &&
+        N != I  &&
             @warn "$(alg.KrylovAlg) doesn't support right preconditioning."
         Krylov.solve!(args...; M=M,
                       kwargs...)
@@ -172,7 +180,7 @@ function SciMLBase.solve(cache::LinearCache, alg::KrylovJL; kwargs...)
         Krylov.solve!(args...; M=M, N=N,
                       kwargs...)
     elseif cache.cacheval isa Krylov.MinresSolver
-        N != LinearAlgebra.I  &&
+        N != I  &&
             @warn "$(alg.KrylovAlg) doesn't support right preconditioning."
         Krylov.solve!(args...; M=M,
                       kwargs...)
@@ -223,8 +231,16 @@ IterativeSolversJL_MINRES(args...;kwargs...) =
 function init_cacheval(alg::IterativeSolversJL, cache::LinearCache)
     @unpack A, b, u = cache
 
-    Pl = ComposePreconditioner(alg.Pl, cache.Pl)
-    Pr = ComposePreconditioner(alg.Pr, cache.Pr)
+    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
 
     abstol  = cache.abstol
     reltol  = cache.reltol

From 65022e7af06007050b350f5a5ed30dbfac06f76a Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Mon, 29 Nov 2021 04:39:08 -0500
Subject: [PATCH 13/17] add proper tests

---
 test/runtests.jl | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 3c5eb90e1..f25f91ca9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -166,9 +166,10 @@ end
         P  = LinearSolve.ComposePreconditioner(P1,P2)
         Pi = LinearSolve.InvComposePreconditioner(P)
 
-        @test Pi == LinearSolve.InvComposePreconditioner(P1,P2)
-        @test Pi == inv(P)
-        @test P  == inv(Pi)
+        @test Pi  == LinearSolve.InvComposePreconditioner(P1,P2)
+        @test Pi  == inv(P)
+        @test P   == inv(Pi)
+        @test Pi' == inv(P')
 
         # ComposePreconditioner
         ldiv!(P, x)

From ef35dac09566ffa92aacb77cb0901e61db4dd05a Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Sat, 4 Dec 2021 21:09:25 -0500
Subject: [PATCH 14/17] Krylov.jl eltype issue fixed

---
 Project.toml    | 2 +-
 src/wrappers.jl | 4 ----
 2 files changed, 1 insertion(+), 5 deletions(-)

diff --git a/Project.toml b/Project.toml
index ca80d41a4..2032c270f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,7 +19,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 [compat]
 ArrayInterface = "3"
 IterativeSolvers = "0.9.2"
-Krylov = "0.7"
+Krylov = "0.7.9"
 KrylovKit = "0.5"
 RecursiveFactorization = "0.2"
 Reexport = "1"
diff --git a/src/wrappers.jl b/src/wrappers.jl
index f04bf90ef..e5c64ae8d 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -3,10 +3,6 @@
 
 default_preconditioner(s, isleft) = isleft ? I * s : I * (1/s)
 
-"""
- C  * x = P  * Q  * x
- Ci * x = Qi * Pi * x
-"""
 struct ComposePreconditioner{Ti,To}
     inner::Ti
     outer::To

From 238c1fe225149f1acf1f034d0ddef72f596fba7c Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Sat, 4 Dec 2021 21:34:36 -0500
Subject: [PATCH 15/17] tests

---
 test/runtests.jl | 38 ++++++++++++++++----------------------
 1 file changed, 16 insertions(+), 22 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index f25f91ca9..95dce8686 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -127,41 +127,35 @@ end
 end
 
 @testset "Preconditioners" begin
-    @testset "ScaleVector" begin
+    @testset "scaling_preconditioner" begin
         s = rand()
-        α = rand()
-        β = rand()
 
         x = rand(n,n)
         y = rand(n,n)
 
-        Pl = LinearSolve.default_preconditioner(s, true)
-        Pr = LinearSolve.default_preconditioner(s, false)
+        Pl = LinearSolve.scaling_preconditioner(s, true)
+        Pr = LinearSolve.scaling_preconditioner(s, false)
 
-        mul!(y, Pl, x)
-        mul!(y, Pr, x)
+        mul!(y, Pl, x); @test y ≈ s * x
+        mul!(y, Pr, x); @test y ≈ s \ x
 
-        mul!(y, Pl, x, α, β)
-        mul!(y, Pr, x, α, β)
+        y .= x; ldiv!(Pl, x); @test x ≈ s \ y
+        y .= x; ldiv!(Pr, x); @test x ≈ s * y
 
-        ldiv!(Pl, x)
-        ldiv!(Pr, x)
-
-        ldiv!(y, Pl, x)
-        ldiv!(y, Pr, x)
+        ldiv!(y, Pl, x); @test y ≈ s \ x
+        ldiv!(y, Pr, x); @test y ≈ s * x
 
     end
 
     @testset "ComposePreconditioenr" begin
-        s = rand()
-        α = rand()
-        β = rand()
+        s1 = rand()
+        s2 = rand()
 
         x = rand(n,n)
         y = rand(n,n)
 
-        P1 = LinearSolve.default_preconditioner(s, true)
-        P2 = LinearSolve.default_preconditioner(s, false)
+        P1 = LinearSolve.scaling_preconditioner(s1, true)
+        P2 = LinearSolve.scaling_preconditioner(s2, true)
 
         P  = LinearSolve.ComposePreconditioner(P1,P2)
         Pi = LinearSolve.InvComposePreconditioner(P)
@@ -172,11 +166,11 @@ end
         @test Pi' == inv(P')
 
         # ComposePreconditioner
-        ldiv!(P, x)
-        ldiv!(y, P, x)
+        ldiv!(y, P, x);      @test y ≈ ldiv!(P2, ldiv!(P1, x))
+        y .= x; ldiv!(P, x); @test x ≈ ldiv!(P2, ldiv!(P1, y))
 
         # InvComposePreconditioner
-        mul!(y, Pi, x)
+        mul!(y, Pi, x); @test y ≈ ldiv!(P2, ldiv!(P1, x))
 
     end
 end

From c6c75aeebdd93222678d857aa26c575c9ec0f5d0 Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Sat, 4 Dec 2021 21:35:30 -0500
Subject: [PATCH 16/17] rename default_preconditioner to scaling_preconditioner

---
 src/common.jl   | 6 +++---
 src/wrappers.jl | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/common.jl b/src/common.jl
index c4731565d..6f0b374e1 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -60,7 +60,7 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
                         reltol=√eps(eltype(prob.A)),
                         maxiters=length(prob.b),
                         verbose=false,
-                        preconditioner_scale=one(eltype(prob.A)),
+                        scale=one(eltype(prob.A)),
                         kwargs...,
                        )
     @unpack A, b, u0, p = prob
@@ -71,8 +71,8 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
     isfresh = cacheval === nothing
     Tc = isfresh ? Any : typeof(cacheval)
 
-    Pl = default_preconditioner(preconditioner_scale, true)
-    Pr = default_preconditioner(preconditioner_scale, false)
+    Pl = scaling_preconditioner(scale, true)
+    Pr = scaling_preconditioner(scale, false)
 
     A = alias_A ? A : deepcopy(A)
     b = alias_b ? b : deepcopy(b)
diff --git a/src/wrappers.jl b/src/wrappers.jl
index e5c64ae8d..f780a0083 100644
--- a/src/wrappers.jl
+++ b/src/wrappers.jl
@@ -1,7 +1,7 @@
 
 ## Preconditioners
 
-default_preconditioner(s, isleft) = isleft ? I * s : I * (1/s)
+scaling_preconditioner(s, isleft) = isleft ? I * s : I * (1/s)
 
 struct ComposePreconditioner{Ti,To}
     inner::Ti

From 730b5017bc1d44f2a631d20c6ac137a4d4a10b7b Mon Sep 17 00:00:00 2001
From: Vedant Puri <vedantpuri@gmail.com>
Date: Wed, 8 Dec 2021 05:22:31 -0500
Subject: [PATCH 17/17] update

---
 src/LinearSolve.jl |  1 +
 src/common.jl      | 10 +++----
 src/wrappers.jl    | 65 +++++++++++++++++++++++++++-------------------
 test/runtests.jl   |  9 +++----
 4 files changed, 48 insertions(+), 37 deletions(-)

diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl
index 9ce4d735a..d28f1e58a 100644
--- a/src/LinearSolve.jl
+++ b/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
diff --git a/src/common.jl b/src/common.jl
index 6f0b374e1..ed8bd6a87 100644
--- a/src/common.jl
+++ b/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)
 
diff --git a/src/wrappers.jl b/src/wrappers.jl
index f780a0083..d9378d06c 100644
--- a/src/wrappers.jl
+++ 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,7 +260,7 @@ 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...)
@@ -257,7 +268,7 @@ function init_cacheval(alg::IterativeSolversJL, cache::LinearCache)
         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,
diff --git a/test/runtests.jl b/test/runtests.jl
index 95dce8686..9077b0399 100644
--- a/test/runtests.jl
+++ b/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')