diff --git a/src/LMTR_alg.jl b/src/LMTR_alg.jl
index 621bedb4..fe7c5e4d 100644
--- a/src/LMTR_alg.jl
+++ b/src/LMTR_alg.jl
@@ -1,6 +1,94 @@
-export LMTR
+export LMTR, LMTRSolver, solve!
+
+import SolverCore.solve!
+
+mutable struct LMTRSolver{
+  T <: Real,
+  G <: ShiftedProximableFunction,
+  V <: AbstractVector{T},
+  N,
+  ST <: AbstractOptimizationSolver,
+  PB <: AbstractRegularizedNLPModel,
+} <: AbstractOptimizationSolver
+  xk::V
+  ∇fk::V
+  mν∇fk::V
+  Fk::V
+  Fkn::V
+  ψ::G
+  χ::N
+  xkn::V
+  s::V
+  has_bnds::Bool
+  l_bound::V
+  u_bound::V
+  l_bound_m_x::V
+  u_bound_m_x::V
+  subsolver::ST
+  subpb::PB
+  substats::GenericExecutionStats{T, V, V, T}
+end
+
+function LMTRSolver(
+  reg_nls::AbstractRegularizedNLPModel{T, V};
+  subsolver = R2Solver,
+  χ = NormLinf(one(T))
+) where{T, V}
+  x0 = reg_nls.model.meta.x0
+  l_bound = reg_nls.model.meta.lvar
+  u_bound = reg_nls.model.meta.uvar
+
+  xk = similar(x0)
+  ∇fk = similar(x0)
+  mν∇fk = similar(x0)
+  Fk = similar(x0, reg_nls.model.nls_meta.nequ)
+  Fkn = similar(Fk)
+  xkn = similar(x0)
+  s = similar(x0)
+  has_bnds = any(l_bound .!= T(-Inf)) || any(u_bound .!= T(Inf)) || subsolver == TRDHSolver
+  if has_bnds
+    l_bound_m_x = similar(xk)
+    u_bound_m_x = similar(xk)
+    @. l_bound_m_x = l_bound - x0
+    @. u_bound_m_x = u_bound - x0
+  else
+    l_bound_m_x = similar(xk, 0)
+    u_bound_m_x = similar(xk, 0)
+  end
+
+  ψ =
+    has_bnds ? shifted(reg_nls.h, xk, max.(-one(T), l_bound_m_x), min.(one(T), u_bound_m_x), reg_nls.selected) :
+    shifted(reg_nls.h, xk, one(T), χ)
+  
+  Jk = jac_op_residual(reg_nls.model, xk)
+  sub_nlp = LMModel(Jk, Fk, T(0), xk)
+  subpb = RegularizedNLPModel(sub_nlp, ψ)
+  substats = RegularizedExecutionStats(subpb)
+  subsolver = subsolver(subpb)
+
+  return LMTRSolver{T, typeof(ψ), V, typeof(χ), typeof(subsolver), typeof(subpb)}(
+    xk,
+    ∇fk,
+    mν∇fk,
+    Fk,
+    Fkn,
+    ψ,
+    χ,
+    xkn,
+    s,
+    has_bnds,
+    l_bound,
+    u_bound,
+    l_bound_m_x,
+    u_bound_m_x,
+    subsolver,
+    subpb,
+    substats
+  )
+end
 
 """
+    LMTR(reg_nls; kwargs...)
     LMTR(nls, h, χ, options; kwargs...)
 
 A trust-region Levenberg-Marquardt method for the problem
@@ -17,233 +105,309 @@ At each iteration, a step s is computed as an approximate solution of
 where F(x) and J(x) are the residual and its Jacobian at x, respectively, ψ(s; x) = h(x + s),
 ‖⋅‖ is a user-defined norm and Δ > 0 is a trust-region radius.
 
-### Arguments
-
-* `nls::AbstractNLSModel`: a smooth nonlinear least-squares problem
-* `h`: a regularizer such as those defined in ProximalOperators
-* `χ`: a norm used to define the trust region in the form of a regularizer
-* `options::ROSolverOptions`: a structure containing algorithmic parameters
-
-### Keyword arguments
-
-* `x0::AbstractVector`: an initial guess (default: `nls.meta.x0`)
-* `subsolver_logger::AbstractLogger`: a logger to pass to the subproblem solver
-* `subsolver`: the procedure used to compute a step (`PG`, `R2` or `TRDH`)
-* `subsolver_options::ROSolverOptions`: default options to pass to the subsolver.
-* `selected::AbstractVector{<:Integer}`: (default `1:nls.meta.nvar`).
-
-### Return values
-
-* `xk`: the final iterate
-* `Fobj_hist`: an array with the history of values of the smooth objective
-* `Hobj_hist`: an array with the history of values of the nonsmooth objective
-* `Complex_hist`: an array with the history of number of inner iterations.
+For advanced usage, first define a solver "LMSolver" to preallocate the memory used in the algorithm, and then call `solve!`:
+
+    solver = LMTRSolver(reg_nls; χ = NormLinf(one(T)), subsolver = R2Solver)
+    solve!(solver, reg_nls)
+
+    stats = RegularizedExecutionStats(reg_nls)
+    solve!(solver, reg_nls, stats)
+  
+# Arguments
+* `reg_nls::AbstractRegularizedNLPModel{T, V}`: the problem to solve, see `RegularizedProblems.jl`, `NLPModels.jl`.
+
+# Keyword arguments
+- `x::V = nlp.meta.x0`: the initial guess;
+- `atol::T = √eps(T)`: absolute tolerance;
+- `sub_atol::T = atol`: subsolver absolute tolerance;
+- `rtol::T = √eps(T)`: relative tolerance;
+- `neg_tol::T = zero(T): negative tolerance;
+- `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited);
+- `max_time::Float64 = 30.0`: maximum time limit in seconds;
+- `max_iter::Int = 10000`: maximum number of iterations;
+- `verbose::Int = 0`: if > 0, display iteration details every `verbose` iteration;
+- `Δk::T = eps(T)`: initial value of the trust-region radius;
+- `η1::T = √√eps(T)`: successful iteration threshold;
+- `η2::T = T(0.9)`: very successful iteration threshold;
+- `γ::T = T(3)`: trust-region radius parameter multiplier, Δ := Δ*γ when the iteration is very successful and Δ := Δ/γ when the iteration is unsuccessful;
+- `α::T = 1/eps(T)`: TODO
+- `β::T = 1/eps(T)`: TODO
+- `χ =  NormLinf(1)`: norm used to define the trust-region;`
+- `subsolver::S = R2Solver`: subsolver used to solve the subproblem that appears at each iteration.
+
+The algorithm stops either when `√(ξₖ/νₖ) < atol + rtol*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure.
+
+# Output
+The value returned is a `GenericExecutionStats`, see `SolverCore.jl`.
+
+# Callback
+$(callback_docstring)
 """
 function LMTR(
   nls::AbstractNLSModel,
   h::H,
   χ::X,
   options::ROSolverOptions;
-  x0::AbstractVector = nls.meta.x0,
-  subsolver_logger::Logging.AbstractLogger = Logging.NullLogger(),
-  subsolver = R2,
-  subsolver_options = ROSolverOptions(ϵa = options.ϵa),
-  selected::AbstractVector{<:Integer} = 1:(nls.meta.nvar),
+  kwargs...
 ) where {H, X}
-  start_time = time()
-  elapsed_time = 0.0
-  # initialize passed options
-  ϵ = options.ϵa
-  ϵ_subsolver = subsolver_options.ϵa
-  ϵr = options.ϵr
-  Δk = options.Δk
-  verbose = options.verbose
-  maxIter = options.maxIter
-  maxTime = options.maxTime
-  η1 = options.η1
-  η2 = options.η2
-  γ = options.γ
-  α = options.α
-  θ = options.θ
-  β = options.β
-
-  # store initial values of the subsolver_options fields that will be modified
-  ν_subsolver = subsolver_options.ν
-  ϵa_subsolver = subsolver_options.ϵa
-  Δk_subsolver = subsolver_options.Δk
-
-  local l_bound, u_bound
-  treats_bounds = has_bounds(nls) || subsolver == TRDH
-  if treats_bounds
-    l_bound = nls.meta.lvar
-    u_bound = nls.meta.uvar
-  end
+  kwargs_dict = Dict(kwargs...)
+  selected = pop!(kwargs_dict, :selected, 1:(nls.meta.nvar))
+  reg_nls = RegularizedNLPModel(nls, h, selected)
+  x0 = pop!(kwargs_dict, :x0, nls.meta.x0)
+  return LMTR(
+    reg_nls;
+    x = x0,
+    χ = χ,
+    atol = options.ϵa,
+    rtol = options.ϵr,
+    neg_tol = options.neg_tol,
+    verbose = options.verbose,
+    max_iter = options.maxIter,
+    max_time = options.maxTime,
+    Δk = options.Δk,
+    η1 = options.η1,
+    η2 = options.η2,
+    γ = options.γ,
+    α = options.α,
+    β = options.β,
+    kwargs_dict...,
+  )
+end
 
-  if verbose == 0
-    ptf = Inf
-  elseif verbose == 1
-    ptf = round(maxIter / 10)
-  elseif verbose == 2
-    ptf = round(maxIter / 100)
-  else
-    ptf = 1
+function LMTR(reg_nls::AbstractRegularizedNLPModel{T}; kwargs...) where{T}
+  kwargs_dict = Dict(kwargs...)
+  subsolver = pop!(kwargs_dict, :subsolver, R2Solver)
+  χ = pop!(kwargs_dict, :χ, NormLinf(one(T)))
+  solver = LMTRSolver(reg_nls, χ = χ, subsolver = subsolver)
+  stats = RegularizedExecutionStats(reg_nls)
+  solve!(solver, reg_nls, stats; kwargs_dict...)
+end
+
+function SolverCore.solve!(
+  solver::LMTRSolver{T, G, V},
+  reg_nls::AbstractRegularizedNLPModel{T, V},
+  stats::GenericExecutionStats{T, V};
+  callback = (args...) -> nothing,
+  x::V = reg_nls.model.meta.x0,
+  atol::T = √eps(T),
+  sub_atol::T = atol,
+  rtol::T = √eps(T),
+  neg_tol::T = zero(T),
+  verbose::Int = 0,
+  max_iter::Int = 10000,
+  max_time::Float64 = 30.0,
+  max_eval::Int = -1,
+  Δk::T = T(1),
+  η1::T = √√eps(T),
+  η2::T = T(0.9),
+  γ::T = T(3),
+  α::T = 1 / eps(T),
+  β::T = 1 / eps(T)
+) where {T, G, V}
+  reset!(stats)
+
+  # Retrieve workspace
+  selected = reg_nls.selected
+  h = reg_nls.h
+  nls = reg_nls.model
+
+  xk = solver.xk .= x
+
+  # Make sure ψ has the correct shift 
+  shift!(solver.ψ, xk)
+
+  Fk = solver.Fk
+  Fkn = solver.Fkn
+  ∇fk = solver.∇fk
+  mν∇fk = solver.mν∇fk
+  ψ = solver.ψ
+  χ = solver.χ
+  xkn = solver.xkn
+  s = solver.s
+  has_bnds = solver.has_bnds
+  if has_bnds
+    l_bound_m_x, u_bound_m_x = solver.l_bound_m_x, solver.u_bound_m_x
+    l_bound, u_bound = solver.l_bound, solver.u_bound
+    update_bounds!(l_bound_m_x, u_bound_m_x, false, l_bound, u_bound, xk, Δk)
+    set_bounds!(ψ, l_bound_m_x, u_bound_m_x)
+    set_bounds!(solver.subsolver.ψ, l_bound_m_x, u_bound_m_x)
   end
 
   # initialize parameters
-  xk = copy(x0)
-  hk = h(xk[selected])
+  improper = false
+  hk = @views h(xk[selected])
   if hk == Inf
-    verbose > 0 && @info "LMTR: finding initial guess where nonsmooth term is finite"
-    prox!(xk, h, x0, one(eltype(x0)))
-    hk = h(xk[selected])
+    verbose > 0 && @info "LM: finding initial guess where nonsmooth term is finite"
+    prox!(xk, h, xk, one(T))
+    hk = @views h(xk[selected])
     hk < Inf || error("prox computation must be erroneous")
-    verbose > 0 && @debug "LMTR: found point where h has value" hk
+    verbose > 0 && @debug "LM: found point where h has value" hk
   end
-  hk == -Inf && error("nonsmooth term is not proper")
-
-  xkn = similar(xk)
-  s = zero(xk)
-  ψ =
-    treats_bounds ? shifted(h, xk, max.(-Δk, l_bound - xk), min.(Δk, u_bound - xk), selected) :
-    shifted(h, xk, Δk, χ)
-
-  Fobj_hist = zeros(maxIter)
-  Hobj_hist = zeros(maxIter)
-  Complex_hist = zeros(Int, maxIter)
-  Grad_hist = zeros(Int, maxIter)
-  Resid_hist = zeros(Int, maxIter)
+  improper = (hk == -Inf)
+  improper == true && @warn "LM: Improper term detected"
+  improper == true && return stats
 
   if verbose > 0
-    #! format: off
-    @info @sprintf "%6s %8s %8s %8s %7s %7s %8s %7s %7s %7s %7s %1s" "outer" "inner" "f(x)" "h(x)" "√ξ1" "√ξ" "ρ" "Δ" "‖x‖" "‖s‖" "ν" "TR"
-    #! format: on
+    @info log_header(
+      [:outer, :inner, :fx, :hx, :xi, :ρ, :Δ, :normx, :norms, :ν, :arrow],
+      [Int, Int, T, T, T, T, T, T, T, T, Char],
+      hdr_override = Dict{Symbol, String}(
+        :fx => "f(x)",
+        :hx => "h(x)",
+        :xi => "√(ξ1/ν)",
+        :normx => "‖x‖",
+        :norms => "‖s‖",
+        :arrow => "LMTR",
+      ),
+      colsep = 1,
+    )
   end
 
-  local ξ1
-  k = 0
+  local ξ1::T
+  local ρk::T = zero(T)
 
-  Fk = residual(nls, xk)
-  Fkn = similar(Fk)
+  residual!(nls, xk, Fk)
+  jtprod_residual!(nls, xk, Fk, ∇fk)
   fk = dot(Fk, Fk) / 2
-  Jk = jac_op_residual(nls, xk)
-  ∇fk = Jk' * Fk
-  JdFk = similar(Fk)   # temporary storage
-  Jt_Fk = similar(∇fk)   # temporary storage
 
-  σmax, found_σ = opnorm(Jk)
+  σmax, found_σ = opnorm(solver.subpb.model.J)
   found_σ || error("operator norm computation failed")
-  α⁻¹Δ⁻¹ = 1 / (α * Δk)
-  ν = 1 / (α⁻¹Δ⁻¹ + σmax^2 * (α⁻¹Δ⁻¹ + 1)) # ‖J'J‖ = ‖J‖²
-
-  mν∇fk = -∇fk * ν
-
-  optimal = false
-  tired = k ≥ maxIter || elapsed_time > maxTime
-
-  while !(optimal || tired)
-    k = k + 1
-    elapsed_time = time() - start_time
-    Fobj_hist[k] = fk
-    Hobj_hist[k] = hk
-    Grad_hist[k] = nls.counters.neval_jtprod_residual + nls.counters.neval_jprod_residual
-    Resid_hist[k] = nls.counters.neval_residual
-
-    # model for first prox-gradient iteration
-    φ1(d) = begin
-      jtprod_residual!(nls, xk, Fk, Jt_Fk)
-      dot(Fk, Fk) / 2 + dot(Jt_Fk, d)
-    end
+  ν = α * Δk / (1 + σmax^2 * (α * Δk + 1))
+  @. mν∇fk = -∇fk * ν
+  sqrt_ξ1_νInv = one(T)
 
-    mk1(d) = φ1(d) + ψ(d)
+  set_iter!(stats, 0)
+  start_time = time()
+  set_time!(stats, 0.0)
+  set_objective!(stats, fk + hk)
+  set_solver_specific!(stats, :smooth_obj, fk)
+  set_solver_specific!(stats, :nonsmooth_obj, hk)
 
-    # TODO: reuse residual computation
-    # model for subsequent prox-gradient iterations
-    φ(d) = begin
-      jprod_residual!(nls, xk, d, JdFk)
-      JdFk .+= Fk
-      dot(JdFk, JdFk) / 2
-    end
+  φ1 = let Fk = Fk, ∇fk = ∇fk
+    d -> dot(Fk, Fk) / 2 + dot(∇fk, d) # ∇fk = Jk^T Fk
+  end
 
-    ∇φ!(g, d) = begin
-      jprod_residual!(nls, xk, d, JdFk)
-      JdFk .+= Fk
-      jtprod_residual!(nls, xk, JdFk, g)
-      g
-    end
+  mk1 = let φ1 = φ1, ψ = ψ
+    d -> φ1(d) + ψ(d)
+  end
 
-    mk(d) = φ(d) + ψ(d)
+  mk = let ψ = ψ, solver = solver
+    d -> obj(solver.subpb.model, d) + ψ(d)
+  end
 
-    # Take first proximal gradient step s1 and see if current xk is nearly stationary.
-    # s1 minimizes φ1(d) + ‖d‖² / 2 / ν + ψ(d) ⟺ s1 ∈ prox{νψ}(-ν∇φ1(0))
-    prox!(s, ψ, mν∇fk, ν)
-    ξ1 = fk + hk - mk1(s) + max(1, abs(fk + hk)) * 10 * eps()
-    ξ1 > 0 || error("LMTR: first prox-gradient step should produce a decrease but ξ1 = $(ξ1)")
+  # Take first proximal gradient step s1 and see if current xk is nearly stationary.
+  # s1 minimizes φ1(d) + ‖d‖² / 2 / ν + ψ(d) ⟺ s1 ∈ prox{νψ}(-ν∇φ1(0))
+  prox!(s, ψ, mν∇fk, ν)
+  ξ1 = fk + hk - mk1(s) + max(1, abs(fk + hk)) * 10 * eps()
+  sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν)
+  solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol)
+  (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) &&
+    error("LMTR: prox-gradient step should produce a decrease but ξ1 = $(ξ1)")
+  atol += rtol * sqrt_ξ1_νInv # make stopping test absolute and relative
+  sub_atol += rtol * sqrt_ξ1_νInv
+
+  set_status!(
+    stats,
+    get_status(
+      reg_nls,
+      elapsed_time = stats.elapsed_time,
+      iter = stats.iter,
+      optimal = solved,
+      improper = improper,
+      max_eval = max_eval,
+      max_time = max_time,
+      max_iter = max_iter,
+    ),
+  )
+
+  callback(nls, solver, stats)
+
+  done = stats.status != :unknown
+
+  while !done
 
-    if ξ1 ≥ 0 && k == 1
-      ϵ_increment = ϵr * sqrt(ξ1)
-      ϵ += ϵ_increment  # make stopping test absolute and relative
-      ϵ_subsolver += ϵ_increment
-    end
+    ∆_effective = min(β * χ(s), Δk)
 
-    if sqrt(ξ1) < ϵ
-      # the current xk is approximately first-order stationary
-      optimal = true
-      continue
+    if has_bnds
+      @. l_bound_m_x = l_bound - xk
+      @. u_bound_m_x = u_bound - xk
+      @. l_bound_m_x .= max.(l_bound_m_x, -∆_effective)
+      @. u_bound_m_x .= min.(u_bound_m_x, ∆_effective)
+      set_bounds!(ψ, l_bound_m_x, u_bound_m_x)
+      set_bounds!(solver.subsolver.ψ, l_bound_m_x, u_bound_m_x)
+    else
+      set_radius!(solver.subsolver.ψ, ∆_effective)
+      set_radius!(ψ, ∆_effective)
     end
 
-    subsolver_options.ϵa = k == 1 ? 1.0e-5 : max(ϵ_subsolver, min(1.0e-1, ξ1 / 10))
-    ∆_effective = min(β * χ(s), Δk)
-    treats_bounds ?
-    set_bounds!(ψ, max.(-∆_effective, l_bound - xk), min.(∆_effective, u_bound - xk)) :
-    set_radius!(ψ, ∆_effective)
-    subsolver_options.Δk = ∆_effective / 10
-    subsolver_options.ν = ν
-    subsolver_args = subsolver == TRDH ? (SpectralGradient(1 / ν, nls.meta.nvar),) : ()
-    s, iter, _ = with_logger(subsolver_logger) do
-      subsolver(φ, ∇φ!, ψ, subsolver_args..., subsolver_options, s)
+    if isa(solver.subsolver, TRDHSolver) 
+      solver.subsolver.D.d[1] = 1/ν
+      solve!(
+        solver.subsolver,
+        solver.subpb,
+        solver.substats,
+        x = s, 
+        atol = stats.iter == 0 ? 1.0e-5 : max(sub_atol, min(1.0e-1, ξ1 / 10)),
+        Δk = ∆_effective / 10
+      )
+    else
+      solve!(
+        solver.subsolver,
+        solver.subpb,
+        solver.substats,
+        x = s, 
+        atol = stats.iter == 0 ? 1.0e-5 : max(sub_atol, min(1.0e-1, ξ1 / 10)),
+        ν = ν,
+      )
     end
-    # restore initial values of subsolver_options here so that it is not modified
-    # if there is an error
-    subsolver_options.ν = ν_subsolver
-    subsolver_options.ϵa = ϵa_subsolver
-    subsolver_options.Δk = Δk_subsolver
 
-    Complex_hist[k] = iter
+    s .= solver.substats.solution
 
     sNorm = χ(s)
     xkn .= xk .+ s
     residual!(nls, xkn, Fkn)
     fkn = dot(Fkn, Fkn) / 2
-    hkn = h(xkn[selected])
-    hkn == -Inf && error("nonsmooth term is not proper")
-
-    Δobj = fk + hk - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
-    ξ = fk + hk - mk(s) + max(1, abs(fk + hk)) * 10 * eps() # TODO: isn't mk(s) returned by subsolver?
-
-    @debug "computed step" s norm(s, Inf) Δk
+    hkn = @views h(xkn[selected])
+    mks = mk(s)
 
+    ξ = fk + hk - mks + max(1, abs(hk)) * 10 * eps()
     if (ξ ≤ 0 || isnan(ξ))
       error("LMTR: failed to compute a step: ξ = $ξ")
     end
 
+    Δobj = fk + hk - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
     ρk = Δobj / ξ
 
-    TR_stat = (η2 ≤ ρk < Inf) ? "↗" : (ρk < η1 ? "↘" : "=")
-
-    if (verbose > 0) && (k % ptf == 0)
-      #! format: off
-      @info @sprintf "%6d %8d %8.1e %8.1e %7.1e %7.1e %8.1e %7.1e %7.1e %7.1e %7.1e %1s" k iter fk hk sqrt(ξ1) sqrt(ξ) ρk ∆_effective χ(xk) sNorm ν TR_stat
-      #! format: on
-    end
-
-    if η2 ≤ ρk < Inf
-      Δk = max(Δk, γ * sNorm)
-      !treats_bounds && set_radius!(ψ, Δk)
-    end
-
+    verbose > 0 &&
+      stats.iter % verbose == 0 &&
+      @info log_row(
+        Any[
+          stats.iter,
+          solver.substats.iter,
+          fk,
+          hk,
+          sqrt_ξ1_νInv,
+          ρk,
+          ∆_effective,
+          χ(xk),
+          sNorm,
+          ν,
+          (η2 ≤ ρk < Inf) ? '↗' : (ρk < η1 ? '↘' : '='),
+        ],
+        colsep = 1,
+      )
+    
     if η1 ≤ ρk < Inf
+
       xk .= xkn
-      treats_bounds && set_bounds!(ψ, max.(-Δk, l_bound - xk), min.(Δk, u_bound - xk))
+      if has_bnds
+        @. l_bound_m_x = l_bound - xk
+        @. u_bound_m_x = u_bound - xk
+        @. l_bound_m_x .= max.(l_bound_m_x, -Δk)
+        @. u_bound_m_x .= min.(u_bound_m_x, Δk)
+        set_bounds!(ψ, l_bound_m_x, u_bound_m_x)
+        set_bounds!(solver.subsolver.ψ, l_bound_m_x, u_bound_m_x)
+      end
 
       #update functions
       Fk .= Fkn
@@ -251,58 +415,93 @@ function LMTR(
       hk = hkn
 
       shift!(ψ, xk)
-      Jk = jac_op_residual(nls, xk)
       jtprod_residual!(nls, xk, Fk, ∇fk)
-      σmax, found_σ = opnorm(Jk)
+      
+      σmax, found_σ = opnorm(solver.subpb.model.J)
       found_σ || error("operator norm computation failed")
-      α⁻¹Δ⁻¹ = 1 / (α * Δk)
-      ν = 1 / (α⁻¹Δ⁻¹ + σmax^2 * (α⁻¹Δ⁻¹ + 1)) # ‖J'J‖ = ‖J‖²
-      @. mν∇fk = -∇fk * ν
+    end
+
+    if η2 ≤ ρk < Inf
+      Δk = max(Δk, γ * sNorm)
+      if !has_bnds 
+        set_radius!(ψ, Δk)
+        set_radius!(solver.subsolver.ψ, Δk)
+      end
     end
 
     if ρk < η1 || ρk == Inf
       Δk = Δk / 2
-      treats_bounds ? set_bounds!(ψ, max.(-Δk, l_bound - xk), min.(Δk, u_bound - xk)) :
-      set_radius!(ψ, Δk)
+      if has_bnds
+        @. l_bound_m_x = l_bound - xk
+        @. u_bound_m_x = u_bound - xk
+        @. l_bound_m_x .= max.(l_bound_m_x, -Δk)
+        @. u_bound_m_x .= min.(u_bound_m_x, Δk)
+        set_bounds!(ψ, l_bound_m_x, u_bound_m_x)
+        set_bounds!(solver.subsolver.ψ, l_bound_m_x, u_bound_m_x)
+      else
+        set_radius!(solver.subsolver.ψ, Δk)
+        set_radius!(ψ, Δk)
+      end
     end
 
-    tired = k ≥ maxIter || elapsed_time > maxTime
-  end
+    set_objective!(stats, fk + hk)
+    set_solver_specific!(stats, :smooth_obj, fk)
+    set_solver_specific!(stats, :nonsmooth_obj, hk)
+    set_iter!(stats, stats.iter + 1)
+    set_time!(stats, time() - start_time)
 
-  if verbose > 0
-    if k == 1
-      @info @sprintf "%6d %8s %8.1e %8.1e" k "" fk hk
-    elseif optimal
-      #! format: off
-      @info @sprintf "%6d %8d %8.1e %8.1e %7.1e %7.1e %8s %7.1e %7.1e %7.1e %7.1e" k 1 fk hk sqrt(ξ1) sqrt(ξ1) "" Δk χ(xk) χ(s) ν
-      #! format: on
-      @info "LMTR: terminating with √ξ1 = $(sqrt(ξ1))"
-    end
+    ν = α * Δk / (1 + σmax^2 * (α * Δk + 1))
+    @. mν∇fk = -∇fk * ν
+
+    prox!(s, ψ, mν∇fk, ν)
+    mks = mk1(s)
+
+    ξ1 = fk + hk - mks + max(1, abs(hk)) * 10 * eps()
+    sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν)
+    solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol)
+    (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) &&
+      error("LM: prox-gradient step should produce a decrease but ξ1 = $(ξ1)")
+
+    set_status!(
+      stats,
+      get_status(
+        reg_nls,
+        elapsed_time = stats.elapsed_time,
+        iter = stats.iter,
+        optimal = solved,
+        improper = improper,
+        max_eval = max_eval,
+        max_time = max_time,
+        max_iter = max_iter,
+      ),
+    )
+
+    callback(nls, solver, stats)
+
+    done = stats.status != :unknown
   end
 
-  status = if optimal
-    :first_order
-  elseif elapsed_time > maxTime
-    :max_time
-  elseif tired
-    :max_iter
-  else
-    :exception
+  if verbose > 0 && stats.status == :first_order
+    @info log_row(
+      Any[
+        stats.iter,
+        0,
+        fk,
+        hk,
+        sqrt_ξ1_νInv,
+        ρk,
+        Δk,
+        χ(xk),
+        χ(s),
+        ν,
+        "",
+      ],
+      colsep = 1,
+    )
+    @info "LMTR: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv)"
   end
 
-  stats = GenericExecutionStats(nls)
-  set_status!(stats, status)
   set_solution!(stats, xk)
-  set_objective!(stats, fk + hk)
-  set_residuals!(stats, zero(eltype(xk)), ξ1 ≥ 0 ? sqrt(ξ1) : ξ1)
-  set_iter!(stats, k)
-  set_time!(stats, elapsed_time)
-  set_solver_specific!(stats, :radius, Δk)
-  set_solver_specific!(stats, :Fhist, Fobj_hist[1:k])
-  set_solver_specific!(stats, :Hhist, Hobj_hist[1:k])
-  set_solver_specific!(stats, :NonSmooth, h)
-  set_solver_specific!(stats, :SubsolverCounter, Complex_hist[1:k])
-  set_solver_specific!(stats, :NLSGradHist, Grad_hist[1:k])
-  set_solver_specific!(stats, :ResidHist, Resid_hist[1:k])
+  set_residuals!(stats, zero(T), sqrt_ξ1_νInv)
   return stats
-end
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 62dfbaed..4daaa1f5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -100,20 +100,9 @@ for (h, h_name) ∈ ((NormL1(λ), "l1"),)
     p = randperm(bpdn_nls.meta.nvar)[1:nz]
     x0[p[1:nz]] = sign.(randn(nz))  # initial guess with nz nonzeros (necessary for h = B0)
     LMTR_out = LMTR(bpdn_nls, h, NormL2(1.0), options, x0 = x0)
-    @test typeof(LMTR_out.solution) == typeof(bpdn_nls.meta.x0)
-    @test length(LMTR_out.solution) == bpdn_nls.meta.nvar
-    @test typeof(LMTR_out.solver_specific[:Fhist]) == typeof(LMTR_out.solution)
-    @test typeof(LMTR_out.solver_specific[:Hhist]) == typeof(LMTR_out.solution)
-    @test typeof(LMTR_out.solver_specific[:SubsolverCounter]) == Array{Int, 1}
+    @test typeof(LMTR_out.solution) == typeof(bpdn.meta.x0)
+    @test length(LMTR_out.solution) == bpdn.meta.nvar
     @test typeof(LMTR_out.dual_feas) == eltype(LMTR_out.solution)
-    @test length(LMTR_out.solver_specific[:Fhist]) == length(LMTR_out.solver_specific[:Hhist])
-    @test length(LMTR_out.solver_specific[:Fhist]) ==
-          length(LMTR_out.solver_specific[:SubsolverCounter])
-    @test length(LMTR_out.solver_specific[:Fhist]) == length(LMTR_out.solver_specific[:NLSGradHist])
-    @test LMTR_out.solver_specific[:NLSGradHist][end] ==
-          bpdn_nls.counters.neval_jprod_residual + bpdn_nls.counters.neval_jtprod_residual - 1
-    @test obj(bpdn_nls, LMTR_out.solution) == LMTR_out.solver_specific[:Fhist][end]
-    @test h(LMTR_out.solution) == LMTR_out.solver_specific[:Hhist][end]
     @test LMTR_out.status == :first_order
   end
 end
diff --git a/test/test_allocs.jl b/test/test_allocs.jl
index 10e477ea..795d4528 100644
--- a/test/test_allocs.jl
+++ b/test/test_allocs.jl
@@ -91,3 +91,18 @@ end
     end
   end
 end
+
+@testset "NLS allocs" begin
+  for (h, h_name) ∈ ((NormL0(λ), "l0"),)
+    for (solver, solver_name) ∈ ((:LMTRSolver, "LMTR"), )
+      @testset "$(solver_name)" begin
+        solver_name == "LMTR" && continue #FIXME
+        reg_nlp = RegularizedNLPModel(bpdn_nls, h)
+        solver = eval(solver)(reg_nlp)
+        stats = RegularizedExecutionStats(reg_nlp)
+        @test @wrappedallocs(solve!(solver, reg_nlp, stats, Δk = 1.0, atol = 1e-6, rtol = 1e-6)) == 0
+        @test stats.status == :first_order
+      end
+    end
+  end
+end
\ No newline at end of file
diff --git a/test/test_bounds.jl b/test/test_bounds.jl
index 110546a2..de2472d5 100644
--- a/test/test_bounds.jl
+++ b/test/test_bounds.jl
@@ -66,23 +66,11 @@ for (h, h_name) ∈ ((NormL0(λ), "l0"), (NormL1(λ), "l1"))
       NormLinf(1.0),
       options,
       x0 = x0,
-      subsolver = TRDH,
-      subsolver_options = subsolver_options,
+      subsolver = TRDHSolver,
     )
-    @test typeof(LMTR_out.solution) == typeof(bpdn_nls2.meta.x0)
-    @test length(LMTR_out.solution) == bpdn_nls2.meta.nvar
-    @test typeof(LMTR_out.solver_specific[:Fhist]) == typeof(LMTR_out.solution)
-    @test typeof(LMTR_out.solver_specific[:Hhist]) == typeof(LMTR_out.solution)
-    @test typeof(LMTR_out.solver_specific[:SubsolverCounter]) == Array{Int, 1}
+    @test typeof(LMTR_out.solution) == typeof(bpdn.meta.x0)
+    @test length(LMTR_out.solution) == bpdn.meta.nvar
     @test typeof(LMTR_out.dual_feas) == eltype(LMTR_out.solution)
-    @test length(LMTR_out.solver_specific[:Fhist]) == length(LMTR_out.solver_specific[:Hhist])
-    @test length(LMTR_out.solver_specific[:Fhist]) ==
-          length(LMTR_out.solver_specific[:SubsolverCounter])
-    @test length(LMTR_out.solver_specific[:Fhist]) == length(LMTR_out.solver_specific[:NLSGradHist])
-    @test LMTR_out.solver_specific[:NLSGradHist][end] ==
-          bpdn_nls2.counters.neval_jprod_residual + bpdn_nls2.counters.neval_jtprod_residual - 1
-    @test obj(bpdn_nls2, LMTR_out.solution) == LMTR_out.solver_specific[:Fhist][end]
-    @test h(LMTR_out.solution) == LMTR_out.solver_specific[:Hhist][end]
     @test LMTR_out.status == :first_order
   end
 end