diff --git a/src/cg.jl b/src/cg.jl index 1345a6232..50517e427 100644 --- a/src/cg.jl +++ b/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)) diff --git a/src/cgls.jl b/src/cgls.jl index e36d5acbd..c9355e52d 100644 --- a/src/cgls.jl +++ b/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 diff --git a/src/cr.jl b/src/cr.jl index d1d912464..d816ce480 100644 --- a/src/cr.jl +++ b/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 diff --git a/src/crls.jl b/src/crls.jl index bf43fa79b..40cca6f2e 100644 --- a/src/crls.jl +++ b/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 diff --git a/src/krylov_utils.jl b/src/krylov_utils.jl index ca73bc121..239e56b75 100644 --- a/src/krylov_utils.jl +++ b/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 diff --git a/src/lsmr.jl b/src/lsmr.jl index 085d941db..1e28c3b22 100644 --- a/src/lsmr.jl +++ b/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) diff --git a/src/lsqr.jl b/src/lsqr.jl index fe7acc37c..94622971e 100644 --- a/src/lsqr.jl +++ b/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) diff --git a/test/test_aux.jl b/test/test_aux.jl index 6c43142c0..9f5481e1d 100644 --- a/test/test_aux.jl +++ b/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