add cglsshifts and lsqrshifts

Author: tmigot

src/PreProcess/PreProcessLSKARC.jl

@@ -0,0 +1,47 @@
+function preprocess(PData::PDataLSKARC, Jx, Fx, gNorm2, calls, max_calls, α)
+    ζ, ξ, maxtol, mintol = PData.ζ, PData.ξ, PData.maxtol, PData.mintol
+    nshifts = PData.nshifts
+    shifts = PData.shifts
+
+    m, n = size(Jx)
+    # Tolerance used in Assumption 2.6b in the paper ( ξ > 0, 0 < ζ ≤ 1 )
+    atol = PData.cgatol(ζ, ξ, maxtol, mintol, gNorm2)
+    rtol = PData.cgrtol(ζ, ξ, maxtol, mintol, gNorm2)
+
+    nshifts = length(shifts)
+    cb = (slv) -> begin
+        ind = setdiff(1:length(shifts), findall(.!slv.converged))
+        if length(ind) > 1
+            for i in ind
+                if (norm(slv.x[i]) / shifts[i] - α > 0)
+                    return true
+                end
+            end
+        end
+        return false
+    end
+    solver = PData.solver
+    solve!(
+        solver,
+        Jx,
+        Fx,
+        shifts,
+        itmax = min(max_calls - sum(calls), 2 * (m + n)),
+        atol = atol,
+        rtol = rtol,
+        verbose = 0,
+        callback = cb,
+    )
+
+    PData.indmin = 0
+    PData.positives .= solver.converged
+    for i = 1:nshifts
+        @. PData.xShift[i] = -solver.x[i]
+        PData.norm_dirs[i] = norm(solver.x[i]) # Get it from cg_lanczos?
+    end
+    PData.shifts .= shifts
+    PData.nshifts = nshifts
+    PData.OK = sum(solver.converged) != 0 # at least one system was solved
+
+    return PData
+end

src/SolveModel/SolveModelLSKARC.jl

@@ -0,0 +1,23 @@
+function solve_modelLSKARC(X::PDataLSKARC, Jx, Fx, gNorm2, calls, max_calls, α::T) where {T}
+    # target value should be close to satisfy αλ=||d||
+    start = findfirst(X.positives)
+    if isnothing(start)
+        start = 1
+    end
+    if VERSION < v"1.7.0"
+        positives = collect(start:length(X.positives))
+        target = [(abs(α * X.shifts[i] - X.norm_dirs[i])) for i in positives]
+    else
+        positives = start:length(X.positives)
+        target = ((abs(α * X.shifts[i] - X.norm_dirs[i])) for i in positives)
+    end
+
+    # pick the closest shift to the target within positive definite H+λI
+    indmin = argmin(target)
+    X.indmin = start + indmin - 1
+    p_imin = X.indmin
+    X.d .= X.xShift[p_imin]
+    X.λ = X.shifts[p_imin]
+
+    return X.d, X.λ
+end

src/solvers.jl

@@ -41,4 +41,6 @@ const solvers_nls_const = Dict(
     :ST_TROpGN => (HessGaussNewtonOp, PDataST, solve_modelST_TR, ()),
     :ST_TROpGNLS => (HessGaussNewtonOp, PDataNLSST, solve_modelNLSST_TR, ()),
     :ST_TROpLS => (HessOp, PDataNLSST, solve_modelNLSST_TR, ()),
+    :LSARCqKOpCgls => (HessGaussNewtonOp, PDataLSKARC, solve_modelLSKARC, [:shifts => 10.0 .^ (collect(-10.0:0.5:20.0)), :solver_method => :cgls]),
+    :LSARCqKOpLsqr => (HessGaussNewtonOp, PDataLSKARC, solve_modelLSKARC, [:shifts => 10.0 .^ (collect(-10.0:0.5:20.0)), :solver_method => :lsqr]),
 )

src/utils/pdata_struct.jl

@@ -115,6 +115,82 @@ function PDataKARC(
     )
 end
 
+"""
+    PDataLSKARC(::Type{S}, ::Type{T}, n)
+Return a structure used for the preprocessing of ARCqK methods.
+"""
+mutable struct PDataLSKARC{T} <: PDataIter{T}
+    d::Array{T,1}             # (H+λI)\g ; on first call = g
+    λ::T                      # "active" value of λ; on first call = 0
+    ζ::T                      # Inexact Newton order parameter: stop when ||∇q|| < ξ * ||g||^(1+ζ)
+    ξ::T                      # Inexact Newton order parameter: stop when ||∇q|| < ξ * ||g||^(1+ζ)
+    maxtol::T                 # Largest tolerance for Inexact Newton
+    mintol::T                 # Smallest tolerance for Inexact Newton
+    cgatol::Any
+    cgrtol::Any
+
+    indmin::Int               # index of best shift value  within "positive". On first call = 0
+
+    positives::Array{Bool,1}   # indices of the shift values yielding (H+λI)⪰0
+    xShift::Array{Array{T,1},1}        # solutions for each shifted system
+    shifts::Array{T,1}        # values of the shifts
+    nshifts::Int              # number of shifts
+    norm_dirs::Array{T,1}     # norms of xShifts
+    OK::Bool                  # preprocess success
+
+    solver::Union{CglsLanczosShiftSolver,LsqrShiftSolver}
+end
+
+function PDataLSKARC(
+    ::Type{S},
+    ::Type{T},
+    n;
+    ζ = T(0.5),
+    ξ = T(0.01),
+    maxtol = T(0.01),
+    mintol = sqrt(eps(T)),
+    cgatol = (ζ, ξ, maxtol, mintol, gNorm2) -> max(mintol, min(maxtol, ξ * gNorm2^(1 + ζ))),
+    cgrtol = (ζ, ξ, maxtol, mintol, gNorm2) -> max(mintol, min(maxtol, ξ * gNorm2^ζ)),
+    shifts = 10.0 .^ collect(-20.0:1.0:20.0),
+    solver_method = :cgls,
+    kwargs...,
+) where {S,T}
+    d = S(undef, n)
+    λ = zero(T)
+    indmin = 1
+    nshifts = length(shifts)
+    positives = Array{Bool,1}(undef, nshifts)
+    xShift = Array{S,1}(undef, nshifts)
+    for i = 1:nshifts
+        xShift[i] = S(undef, n)
+    end
+    norm_dirs = S(undef, nshifts)
+    OK = true
+    solver = if solver_method == :cgls
+        CglsLanczosShiftSolver(n, n, nshifts, S)
+    else
+        LsqrShiftSolver(n, n, nshifts, S)
+    end
+    return PDataLSKARC(
+        d,
+        λ,
+        ζ,
+        ξ,
+        maxtol,
+        mintol,
+        cgatol,
+        cgrtol,
+        indmin,
+        positives,
+        xShift,
+        T.(shifts),
+        nshifts,
+        norm_dirs,
+        OK,
+        solver,
+    )
+end
+
 """
     PDataTRK(::Type{S}, ::Type{T}, n)
 Return a structure used for the preprocessing of TRK methods.