Fix preconditioning in cg when optional radius argument is positive

src/cg.jl

@@ -210,9 +210,15 @@ kwargs_cg = (:M, :ldiv, :radius, :linesearch, :atol, :rtol, :itmax, :timemax, :v
       (zero_curvature || solved) && continue
 
       α = γ / pAp
-
+ 
       # Compute step size to boundary if applicable.
-      σ = radius > 0 ? maximum(to_boundary(n, x, p, radius, dNorm2=pNorm²)) : α
+      if radius == 0
+         σ = α
+      elseif MisI 
+         σ = maximum(to_boundary(n, x, p, z, radius, dNorm2=pNorm²))
+      else
+         σ = maximum(to_boundary(n, x, p, z, radius, M=M, ldiv=!ldiv)) 
+      end
 
       kdisplay(iter, verbose) && @printf(iostream, "  %8.1e  %8.1e  %8.1e  %.2fs\n", pAp, α, σ, ktimer(start_time))
 

src/cgls.jl

@@ -198,7 +198,7 @@ kwargs_cgls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose
       α = γ / δ
 
       # if a trust-region constraint is give, compute step to the boundary
-      σ = radius > 0 ? maximum(to_boundary(n, x, p, radius)) : α
+      σ = radius > 0 ? maximum(to_boundary(n, x, p, Mr, radius)) : α
       if (radius > 0) & (α > σ)
         α = σ
         on_boundary = true

src/cr.jl

@@ -236,10 +236,10 @@ kwargs_cr = (:M, :ldiv, :radius, :linesearch, :γ, :atol, :rtol, :itmax, :timema
         (verbose > 0) && @printf(iostream, "radius = %8.1e > 0 and ‖x‖ = %8.1e\n", radius, xNorm)
         # find t1 > 0 and t2 < 0 such that ‖x + ti * p‖² = radius²  (i = 1, 2)
         xNorm² = xNorm * xNorm
-        t = to_boundary(n, x, p, radius; flip = false, xNorm2 = xNorm², dNorm2 = pNorm²)
+        t = to_boundary(n, x, p, Mq, radius; flip = false, xNorm2 = xNorm², dNorm2 = pNorm²)
         t1 = maximum(t) # > 0
         t2 = minimum(t) # < 0
-        tr = maximum(to_boundary(n, x, r, radius; flip = false, xNorm2 = xNorm², dNorm2 = rNorm²))
+        tr = maximum(to_boundary(n, x, r, Mq, radius; flip = false, xNorm2 = xNorm², dNorm2 = rNorm²))
         (verbose > 0) && @printf(iostream, "t1 = %8.1e, t2 = %8.1e and tr = %8.1e\n", t1, t2, tr)
 
         if abspAp ≤ γ * pNorm * @knrm2(n, q) # pᴴAp ≃ 0

src/crls.jl

@@ -204,10 +204,10 @@ kwargs_crls = (:M, :ldiv, :radius, :λ, :atol, :rtol, :itmax, :timemax, :verbose
           p = Ar # p = Aᴴr
           pNorm² = ArNorm * ArNorm
           mul!(q, Aᴴ, s)
-          α = min(ArNorm^2 / γ, maximum(to_boundary(n, x, p, radius, flip = false, dNorm2 = pNorm²))) # the quadratic is minimal in the direction Aᴴr for α = ‖Ar‖²/γ
+          α = min(ArNorm^2 / γ, maximum(to_boundary(n, x, p, Ms, radius, flip = false, dNorm2 = pNorm²))) # the quadratic is minimal in the direction Aᴴr for α = ‖Ar‖²/γ
         else
           pNorm² = pNorm * pNorm
-          σ = maximum(to_boundary(n, x, p, radius, flip = false, dNorm2 = pNorm²))
+          σ = maximum(to_boundary(n, x, p, Ms, radius, flip = false, dNorm2 = pNorm²))
           if α ≥ σ
             α = σ
             on_boundary = true

src/krylov_utils.jl

@@ -405,21 +405,32 @@ If `flip` is set to `true`, `σ1` and `σ2` are computed such that
 
     ‖x - σi d‖ = radius, i = 1, 2.
 """
-function to_boundary(n :: Int, x :: AbstractVector{FC}, d :: AbstractVector{FC}, radius :: T; flip :: Bool=false, xNorm2 :: T=zero(T), dNorm2 :: T=zero(T)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}}
+function to_boundary(n :: Int, x :: AbstractVector{FC}, d :: AbstractVector{FC}, z::AbstractVector{FC}, radius :: T; flip :: Bool=false, xNorm2 :: T=zero(T), dNorm2 :: T=zero(T), M=I, ldiv :: Bool = false) where {T <: AbstractFloat, FC <: FloatOrComplex{T}}
   radius > 0 || error("radius must be positive")
 
-  # ‖d‖² σ² + (xᴴd + dᴴx) σ + (‖x‖² - Δ²).
-  rxd = @kdotr(n, x, d)
-  flip && (rxd = -rxd)
-  dNorm2 == zero(T) && (dNorm2 = @kdotr(n, d, d))
+  if M === I
+    # ‖d‖² σ² + (xᴴd + dᴴx) σ + (‖x‖² - Δ²).
+    rxd = @kdotr(n, x, d)
+    dNorm2 == zero(T) && (dNorm2 = @kdotr(n, d, d))
+    xNorm2 == zero(T) && (xNorm2 = @kdotr(n, x, x))
+  else
+    # (dᴴMd) σ² + (xᴴMd + dᴴMx) σ + (xᴴMx - Δ²).
+    mulorldiv!(z, M, x, ldiv)
+    rxd = dot(z, d)
+    xNorm2 = dot(z, x)
+    mulorldiv!(z, M, d, ldiv) 
+    dNorm2 = dot(z, d)
+  end 
   dNorm2 == zero(T) && error("zero direction")
-  xNorm2 == zero(T) && (xNorm2 = @kdotr(n, x, x))
+  flip && (rxd = -rxd)
+
   radius2 = radius * radius
   (xNorm2 ≤ radius2) || error(@sprintf("outside of the trust region: ‖x‖²=%7.1e, Δ²=%7.1e", xNorm2, radius2))
 
   # q₂ = ‖d‖², q₁ = xᴴd + dᴴx, q₀ = ‖x‖² - Δ²
   # ‖x‖² ≤ Δ² ⟹ (q₁)² - 4 * q₂ * q₀ ≥ 0
   roots = roots_quadratic(dNorm2, 2 * rxd, xNorm2 - radius2)
+
   return roots  # `σ1` and `σ2`
 end
 

src/lsmr.jl

@@ -349,7 +349,7 @@ kwargs_lsmr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim,
       # the step ϕ/ρ is not necessarily positive
       σ = ζ / (ρ * ρbar)
       if radius > 0
-        t1, t2 = to_boundary(n, x, hbar, radius)
+        t1, t2 = to_boundary(n, x, hbar, v, radius)
         tmax, tmin = max(t1, t2), min(t1, t2)
         on_boundary = σ > tmax || σ < tmin
         σ = σ > 0 ? min(σ, tmax) : max(σ, tmin)

src/lsqr.jl

@@ -346,7 +346,7 @@ kwargs_lsqr = (:M, :N, :ldiv, :sqd, :λ, :radius, :etol, :axtol, :btol, :conlim,
       # the step ϕ/ρ is not necessarily positive
       σ = ϕ / ρ
       if radius > 0
-        t1, t2 = to_boundary(n, x, w, radius)
+        t1, t2 = to_boundary(n, x, w, v, radius)
         tmax, tmin = max(t1, t2), min(t1, t2)
         on_boundary = σ > tmax || σ < tmin
         σ = σ > 0 ? min(σ, tmax) : max(σ, tmin)

test/test_aux.jl

@@ -106,13 +106,14 @@
     n = 5
     x = ones(n)
     d = ones(n); d[1:2:n] .= -1
-    @test_throws ErrorException Krylov.to_boundary(n, x, d, -1.0)
-    @test_throws ErrorException Krylov.to_boundary(n, x, d, 0.5)
-    @test_throws ErrorException Krylov.to_boundary(n, x, zeros(n), 1.0)
-    @test maximum(Krylov.to_boundary(n, x, d, 5.0)) ≈ 2.209975124224178
-    @test minimum(Krylov.to_boundary(n, x, d, 5.0)) ≈ -1.8099751242241782
-    @test maximum(Krylov.to_boundary(n, x, d, 5.0, flip=true)) ≈ 1.8099751242241782
-    @test minimum(Krylov.to_boundary(n, x, d, 5.0, flip=true)) ≈ -2.209975124224178
+    z = similar(d) # <-- placeholder for preconditioning storage
+    @test_throws ErrorException Krylov.to_boundary(n, x, d, z, -1.0)
+    @test_throws ErrorException Krylov.to_boundary(n, x, d, z, 0.5)
+    @test_throws ErrorException Krylov.to_boundary(n, x, zeros(n), z, 1.0)
+    @test maximum(Krylov.to_boundary(n, x, d, z, 5.0)) ≈ 2.209975124224178
+    @test minimum(Krylov.to_boundary(n, x, d, z, 5.0)) ≈ -1.8099751242241782
+    @test maximum(Krylov.to_boundary(n, x, d, z, 5.0, flip=true)) ≈ 1.8099751242241782
+    @test minimum(Krylov.to_boundary(n, x, d, z, 5.0, flip=true)) ≈ -2.209975124224178
   end
 
   @testset "kzeros" begin