diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index b7387f3da..2ef8cf900 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -52,7 +52,7 @@ export r_LL, l_LL, r_RR, l_RR, r_RL, r_LR, l_RL, l_LR # TODO: rename
 
 # unexported
 using Compat: @compat
-@compat public DynamicTols
+@compat public DynamicTols, DynamicTol, AdaptiveKrylov
 @compat public VERBOSE_NONE, VERBOSE_WARN, VERBOSE_CONV, VERBOSE_ITER, VERBOSE_ALL
 @compat public IterLog, loginit!, logiter!, logfinish!, logcancel!
 
diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl
index 5a9ec4d29..d5ffc6cc4 100644
--- a/src/algorithms/approximate/fvomps.jl
+++ b/src/algorithms/approximate/fvomps.jl
@@ -8,7 +8,7 @@ function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environment
             ϵ = 0.0
             for pos in [1:(length(ψ) - 1); (length(ψ) - 2):-1:1]
                 AC2′ = AC2_projection(pos, ψ, Oϕ, envs)
-                al, c, ar = svd_trunc!(AC2′; trunc = alg.trscheme, alg = alg.alg_svd)
+                al, c, ar, = svd_trunc!(AC2′, alg.alg_gauge)
 
                 AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
                 ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2))
diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl
index f437a55fa..de53ad77a 100644
--- a/src/algorithms/approximate/idmrg.jl
+++ b/src/algorithms/approximate/idmrg.jl
@@ -52,7 +52,7 @@ function approximate!(
     end
 
     # TODO: immediately compute in-place
-    alg_gauge = updatetol(alg.alg_gauge, iter, ϵ)
+    alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
     ψ′ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
     copy!(ψ, ψ′) # ensure output destination is unchanged
 
@@ -179,7 +179,7 @@ function approximate!(
     end
 
     # TODO: immediately compute in-place
-    alg_gauge = updatetol(alg.alg_gauge, iter, ϵ)
+    alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
     ψ′ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
     copy!(ψ, ψ′) # ensure output destination is unchanged
 
diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl
index 455265256..70f3335b0 100644
--- a/src/algorithms/approximate/vomps.jl
+++ b/src/algorithms/approximate/vomps.jl
@@ -17,7 +17,7 @@ function approximate(
     log = IterLog("VOMPS")
     iter = 0
     ϵ = calc_galerkin(mps, toapprox..., envs)
-    alg_environments = updatetol(alg.alg_environments, iter, ϵ)
+    alg_environments = adapt_solver(alg.alg_environments; iter, g_global = ϵ)
     recalculate!(envs, mps, toapprox..., alg_environments)
 
     state = VOMPSState(mps, toapprox, envs, iter, ϵ)
@@ -65,7 +65,7 @@ end
 function localupdate_step!(
         it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, ::SerialScheduler
     )
-    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
+    alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ)
     alg_orth = alg_gauge.alg_orth
 
     ACs = similar(state.mps.AC)
@@ -89,7 +89,7 @@ end
 function localupdate_step!(
         it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, scheduler
     )
-    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
+    alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ)
     alg_orth = alg_gauge.alg_orth
 
     ACs = similar(state.mps.AC)
@@ -119,6 +119,6 @@ function localupdate_step!(
 end
 
 function envs_step!(it::IterativeSolver{<:VOMPS}, state::VOMPSState{<:Any, <:Tuple}, mps)
-    alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ)
+    alg_environments = adapt_solver(it.alg_environments; iter = state.iter, g_global = state.ϵ)
     return recalculate!(state.envs, mps, state.operator..., alg_environments)
 end
diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl
index d88b3aa3b..87d54ee5b 100644
--- a/src/algorithms/changebonds/optimalexpand.jl
+++ b/src/algorithms/changebonds/optimalexpand.jl
@@ -114,7 +114,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Op
     # embed `left` into the enlarged domain (zero weight in the new directions), leaving the state
     # unchanged
     nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1))
-    nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
+    nal, nc, _ = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
     nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re))
 
     normalize && normalize!(nc)
diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl
index 2ab840ffe..85bbfa81c 100644
--- a/src/algorithms/changebonds/randexpand.jl
+++ b/src/algorithms/changebonds/randexpand.jl
@@ -98,7 +98,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Ra
     ar_re = Vᴴ * NR
     # embed `left` into the enlarged domain (zero weight in the new directions)
     nal_space = codomain(left) ← (only(domain(left)) ⊕ space(Vᴴ, 1))
-    nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
+    nal, nc, _ = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
     nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re))
 
     normalize && normalize!(nc)
diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl
index 9e89b3b81..deade493f 100644
--- a/src/algorithms/changebonds/sketchedexpand.jl
+++ b/src/algorithms/changebonds/sketchedexpand.jl
@@ -53,7 +53,7 @@ end
 function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::SketchedExpand, envs; normalize::Bool = true)
     left = ψ.AC[site]
     right = ψ.AR[site + 1]
-    AL, _ = left_gauge(left)        # local left-isometric form
+    AL, _, _ = left_gauge(left)        # local left-isometric form
     ARtt = _transpose_tail(right)   # AR is already right-isometric
 
     # nothing to add when either complement is empty (e.g. edge bonds)
@@ -81,7 +81,7 @@ function changebond!(site::Int, ::Val{:right}, ψ::AbstractFiniteMPS, H, alg::Sk
     # embed `left` into the enlarged domain (zero weight in the new directions), leaving the state
     # unchanged
     nal_space = codomain(left) ← (only(domain(left)) ⊕ space(ar_re, 1))
-    nal, nc = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
+    nal, nc, _ = left_gauge(absorb!(zerovector!(similar(left, nal_space)), left))
     nar = _transpose_front(catcodomain(_transpose_tail(right), ar_re))
 
     normalize && normalize!(nc)
diff --git a/src/algorithms/excitation/dmrgexcitation.jl b/src/algorithms/excitation/dmrgexcitation.jl
index 8bdb934c6..72c47a401 100644
--- a/src/algorithms/excitation/dmrgexcitation.jl
+++ b/src/algorithms/excitation/dmrgexcitation.jl
@@ -18,12 +18,21 @@ $(TYPEDFIELDS)
     weight::Float64 = 10.0
 end
 
+# Initialize excited state by perturbing the current eigenvector to avoid local minima
+function _perturbed_state(ψ::FiniteMPS; atol = 1.0e-2)
+    return FiniteMPS(
+        map(1:length(ψ)) do i
+            A = i == length(ψ) ? ψ.AC[i] : ψ.AL[i]
+            noise = randomize!(similar(A))
+            return add!(noise, A, 1, atol / norm(noise))
+        end
+    )
+end
+
 function excitations(
         H::FiniteMPOHamiltonian, alg::FiniteExcited,
         states::Tuple{T, Vararg{T}};
-        init = FiniteMPS(
-            [ copy(first(states).AC[i]) for i in 1:length(first(states)) ]
-        ), num = 1
+        init = _perturbed_state(first(states)), num = 1
     ) where {T <: FiniteMPS}
     num == 0 && return (scalartype(T)[], T[])
 
diff --git a/src/algorithms/fixedpoint.jl b/src/algorithms/fixedpoint.jl
index 2248d416d..26b692e11 100644
--- a/src/algorithms/fixedpoint.jl
+++ b/src/algorithms/fixedpoint.jl
@@ -1,32 +1,23 @@
 # wrapper around KrylovKit.jl's eigsolve function
 
+function fixedpoint(A, x₀, which::Symbol; kwargs...)
+    alg = KrylovKit.eigselector(A, scalartype(x₀); kwargs...)
+    return fixedpoint(A, x₀, which, alg)
+end
+
 """
-    fixedpoint(A, x₀, which::Symbol; kwargs...) -> val, vec
-    fixedpoint(A, x₀, which::Symbol, alg) -> val, vec
+    fixedpoint(A, x₀, which::Symbol; kwargs...) -> val, vec, info
+    fixedpoint(A, x₀, which::Symbol, alg) -> val, vec, info
 
-Compute the fixed point of a linear operator `A` using the specified eigensolver `alg`. The
-fixedpoint is assumed to be unique.
+Compute the fixed point of a given linear operator `A` with initial guess `x₀`.
+The dominant eigenvector is assumed to be unique.
 """
 function fixedpoint(A, x₀, which::Symbol, alg::Lanczos)
     vals, vecs, info = eigsolve(A, x₀, 1, which, alg)
-
-    info.converged == 0 &&
-        @warnv 1 "fixed point not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
-
-    return vals[1], vecs[1]
+    return vals[1], vecs[1], info
 end
-
 function fixedpoint(A, x₀, which::Symbol, alg::Arnoldi)
     TT, vecs, vals, info = schursolve(A, x₀, 1, which, alg)
-
-    info.converged == 0 &&
-        @warnv 1 "fixed point not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
     size(TT, 2) > 1 && TT[2, 1] != 0 && @warnv 1 "non-unique fixed point detected"
-
-    return vals[1], vecs[1]
-end
-
-function fixedpoint(A, x₀, which::Symbol; kwargs...)
-    alg = KrylovKit.eigselector(A, scalartype(x₀); kwargs...)
-    return fixedpoint(A, x₀, which, alg)
+    return vals[1], vecs[1], info
 end
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index 61f3097a0..de8ac4bd8 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -53,8 +53,10 @@ function DMRG(;
         alg_expand = nothing, trscheme = notrunc(),
         alg_svd = Defaults.alg_svd(), alg_orth = Defaults.alg_orth()
     )
-    alg_eigsolve′ = alg_eigsolve isa NamedTuple ? Defaults.alg_eigsolve(; alg_eigsolve...) :
-        alg_eigsolve
+    # single-site DMRG defaults to the per-bond adaptive controller (`AdaptiveKrylov`); pass
+    # `alg_eigsolve = (; adaptive = false, ...)` to opt out (the splat overrides the default).
+    alg_eigsolve′ = alg_eigsolve isa NamedTuple ?
+        Defaults.alg_eigsolve(; adaptive = true, alg_eigsolve...) : alg_eigsolve
     # a no-truncation `trscheme` selects a (bond-preserving) QR gauge, anything else a truncated SVD
     alg_gauge = trscheme isa MatrixAlgebraKit.NoTruncation ? alg_orth :
         MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
@@ -64,83 +66,38 @@ function DMRG(;
     return DMRG(tol, maxiter, verbosity, alg_eigsolve′, finalize, alg_expand, alg_gauge)
 end
 
-function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG, envs = environments(ψ, H, ψ))
-    ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ))
-    ϵ = maximum(ϵs)
-    log = IterLog("DMRG")
-    timeroutput = TimerOutput("DMRG")
-    alg.verbosity > 3 || disable_timer!(timeroutput)
 
-    LoggingExtras.withlevel(; alg.verbosity) do
-        @infov 2 loginit!(log, ϵ, expectation_value(ψ, H, envs))
-        for iter in 1:(alg.maxiter)
-            alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ)
+function local_update!(
+        site, direction,
+        ψ, O, alg::DMRG, envs,
+        ϵ_global, ϵ_trunc, decay_rate,
+        timeroutput
+    )
+    ϵ_local = calc_galerkin(site, ψ, O, ψ, envs)
+
+    # 1. expand
+    isnothing(alg.alg_expand) ||
+        @timeit timeroutput "expand" changebond!(site, direction, ψ, O, alg.alg_expand, envs)
+
+    # 2. local update
+    alg_eigsolve = adapt_solver(alg.alg_eigsolve; decay_rate, g_local = ϵ_local, g_global = ϵ_global, eps_trunc = ϵ_trunc)
+    ac_old = ψ.AC[site]
+    λ, AC′, info = @timeit timeroutput "AC_eigsolve" begin
+        H_effective = AC_hamiltonian(site, ψ, O, ψ, envs)
+        fixedpoint(H_effective, ac_old, :SR, alg_eigsolve)
+    end
 
-            zerovector!(ϵs)
-            @timeit timeroutput "sweep" begin
-                # left-to-right
-                for pos in 1:(length(ψ) - 1)
-                    local AC′
-                    # convergence: pre-expansion single-site Galerkin error
-                    ϵs[pos] = max(ϵs[pos], calc_galerkin(pos, ψ, H, ψ, envs))
-
-                    # 1. expand
-                    isnothing(alg.alg_expand) ||
-                        @timeit timeroutput "expand" changebond!(pos, Val(:right), ψ, H, alg.alg_expand, envs)
-
-                    # 2. local update
-                    @timeit timeroutput "AC_eigsolve" begin
-                        Hac = AC_hamiltonian(pos, ψ, H, ψ, envs)
-                        _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve)
-                    end
-
-                    # 3. gauge (QR center-move or truncated SVD, selected by `alg_gauge`)
-                    @timeit timeroutput "gauge" left_gauge!(ψ, pos, AC′, alg.alg_gauge; normalize = true)
-                end
+    # 3. gauge
+    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge!(ψ, site, direction, AC′, alg.alg_gauge; normalize = true)
 
-                # right-to-left
-                for pos in length(ψ):-1:2
-                    local AC′
-                    # convergence: pre-expansion single-site Galerkin error
-                    ϵs[pos] = max(ϵs[pos], calc_galerkin(pos, ψ, H, ψ, envs))
-
-                    # 1. expand
-                    isnothing(alg.alg_expand) ||
-                        @timeit timeroutput "expand" changebond!(pos, Val(:left), ψ, H, alg.alg_expand, envs)
-
-                    # 2. local update
-                    @timeit timeroutput "AC_eigsolve" begin
-                        Hac = AC_hamiltonian(pos, ψ, H, ψ, envs)
-                        _, AC′ = fixedpoint(Hac, ψ.AC[pos], :SR, alg_eigsolve)
-                    end
-
-                    # 3. gauge (QR center-move or truncated SVD, selected by `alg_gauge`)
-                    @timeit timeroutput "gauge" right_gauge!(ψ, pos, AC′, alg.alg_gauge; normalize = true)
-                end
-            end
-            ϵ = maximum(ϵs)
+    # 4. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1)
+    decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999)
 
-            ψ, envs = @timeit timeroutput "finalize" alg.finalize(
-                iter, ψ, H, envs
-            )::Tuple{typeof(ψ), typeof(envs)}
+    @debug "DMRG local update" site direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local
 
-            if ϵ <= alg.tol
-                @infov 4 timeroutput
-                @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, H, envs))
-                break
-            end
-            if iter == alg.maxiter
-                @infov 4 timeroutput
-                @warnv 1 logcancel!(log, iter, ϵ, expectation_value(ψ, H, envs))
-            else
-                @infov 3 logiter!(log, iter, ϵ, expectation_value(ψ, H, envs))
-            end
-        end
-    end
-    return ψ, envs, ϵ
+    return ψ, ϵ_local, ϵ_trunc, decay_rate
 end
 
-
 """
 $(TYPEDEF)
 
@@ -150,7 +107,7 @@ Two-site DMRG algorithm for finding the dominant eigenvector.
 
 $(TYPEDFIELDS)
 """
-struct DMRG2{A, S, F} <: Algorithm
+struct DMRG2{A, G, F} <: Algorithm
     "tolerance for convergence criterium"
     tol::Float64
 
@@ -163,11 +120,8 @@ struct DMRG2{A, S, F} <: Algorithm
     "algorithm used for the eigenvalue solvers"
     alg_eigsolve::A
 
-    "algorithm used for the singular value decomposition"
-    alg_svd::S
-
-    "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy) of the two-site update"
-    trscheme::TruncationStrategy
+    "factorization used for the post-update gauge: a truncated SVD (`alg_svd` with `trscheme`)"
+    alg_gauge::G
 
     "callback function applied after each iteration, of signature `finalize(iter, ψ, H, envs) -> ψ, envs`"
     finalize::F
@@ -178,84 +132,126 @@ function DMRG2(;
         alg_eigsolve = (;), alg_svd = Defaults.alg_svd(), trscheme,
         finalize = Defaults._finalize
     )
-    alg_eigsolve′ = alg_eigsolve isa NamedTuple ? Defaults.alg_eigsolve(; alg_eigsolve...) :
-        alg_eigsolve
-    return DMRG2(tol, maxiter, verbosity, alg_eigsolve′, alg_svd, trscheme, finalize)
+    # two-site DMRG defaults to the per-bond adaptive controller (`AdaptiveKrylov`); pass
+    # `alg_eigsolve = (; adaptive = false, ...)` to opt out (the splat overrides the default).
+    alg_eigsolve′ = alg_eigsolve isa NamedTuple ?
+        Defaults.alg_eigsolve(; adaptive = true, alg_eigsolve...) : alg_eigsolve
+    # two-site DMRG always truncates the enlarged bond back down, so the gauge is a truncated SVD
+    alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trscheme)
+    return DMRG2(tol, maxiter, verbosity, alg_eigsolve′, alg_gauge, finalize)
+end
+
+function local_update!(
+        pos, direction,
+        ψ, O, alg::DMRG2, envs,
+        ϵ_global, ϵ_trunc, decay_rate,
+        timeroutput
+    )
+    Heff = @timeit timeroutput "AC2_hamiltonian" AC2_hamiltonian(pos, ψ, O, ψ, envs)
+
+    kind = direction === Val(:right) ? :ACAR : :ALAC
+    ac2 = AC2(ψ, pos; kind)
+    AC2′ = normalize!(Heff * ac2)
+    project_complement!(AC2′, ψ.AL[pos])
+    ϵ_local = norm(AC2′)
+
+    # 1. local two-site update
+    alg_eigsolve = adapt_solver(alg.alg_eigsolve; decay_rate, g_local = ϵ_local, g_global = ϵ_global, eps_trunc = ϵ_trunc)
+    newA2center, info = @timeit timeroutput "AC2_eigsolve" begin
+        _, newA2center, info = fixedpoint(Heff, ac2, :SR, alg_eigsolve)
+        (newA2center, info)
+    end
+
+    # 2. gauge: truncated SVD split back into single-site tensors and install;
+    #           the discarded weight is the truncation error
+    ψ, ϵ_trunc = @timeit timeroutput "gauge" gauge2!(ψ, pos, direction, newA2center, alg.alg_gauge; normalize = true)
+
+    # 3. bookkeeping: measured contraction factor per matvec, kept a strict contraction in (0, 1)
+    decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999)
+
+    @debug "DMRG2 local update" pos direction numops = info.numops normres = first(info.normres) krylovdim = alg_eigsolve.krylovdim maxiter = alg_eigsolve.maxiter tol = alg_eigsolve.tol decay_rate ϵ_local ϵ_trunc
+
+    return ψ, ϵ_local, ϵ_trunc, decay_rate
 end
 
-function find_groundstate!(ψ::AbstractFiniteMPS, H, alg::DMRG2, envs = environments(ψ, H, ψ))
-    ϵs = map(pos -> calc_galerkin(pos, ψ, H, ψ, envs), 1:length(ψ))
-    ϵ = maximum(ϵs)
-    log = IterLog("DMRG2")
-    timeroutput = TimerOutput("DMRG2")
+# Per-algorithm sweep geometry: single-site DMRG updates all `length(ψ)` sites (endpoints once,
+# interior twice), whereas two-site DMRG2 updates the `length(ψ) - 1` bonds. `_num_updates`
+# gives the number of per-update bookkeeping slots and `_sweep_ranges` the forward/backward
+# index ranges; everything else in the sweep is shared.
+_num_updates(::DMRG, ψ) = length(ψ)
+_num_updates(::DMRG2, ψ) = length(ψ) - 1
+
+_sweep_ranges(::DMRG, ψ) = (1:(length(ψ) - 1), length(ψ):-1:2)
+_sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 2):-1:1)
+
+function find_groundstate!(
+        ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs = environments(ψ, H, ψ)
+    )
+    name = string(nameof(typeof(alg)))
+    log = IterLog(name)
+    timeroutput = TimerOutput(name)
     alg.verbosity > 3 || disable_timer!(timeroutput)
 
+    Tr = real(scalartype(ψ))
+    n = _num_updates(alg, ψ)
+    ϵ_locals = ones(Tr, n)      # local Galerkin errors (drive both the eigensolver tol and the stop test)
+    ϵ_global = maximum(ϵ_locals) # sweep-wide Galerkin error, the convergence measure
+    ϵ_truncs = zeros(Tr, n)     # local truncation error
+    decay_rates = zeros(n)      # local observed decay rate of eigensolver
+    fwd, bwd = _sweep_ranges(alg, ψ)
+
     LoggingExtras.withlevel(; alg.verbosity) do
+        @infov 2 loginit!(log, ϵ_global, expectation_value(ψ, H, envs))
         for iter in 1:(alg.maxiter)
-            alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ)
-            zerovector!(ϵs)
-
             @timeit timeroutput "sweep" begin
-                # left to right sweep
-                for pos in 1:(length(ψ) - 1)
-                    local ac2, newA2center, al, c, ar
-                    @timeit timeroutput "AC2_eigsolve" begin
-                        @plansor ac2[-1 -2; -3 -4] := ψ.AC[pos][-1 -2; 1] * ψ.AR[pos + 1][1 -4; -3]
-                        Hac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs)
-                        _, newA2center = fixedpoint(Hac2, ac2, :SR, alg_eigsolve)
-                    end
-                    @timeit timeroutput "svd_trunc" begin
-                        al, c, ar = svd_trunc!(newA2center; trunc = alg.trscheme, alg = alg.alg_svd)
-                        normalize!(c)
-                        v = @plansor ac2[1 2; 3 4] * conj(al[1 2; 5]) * conj(c[5; 6]) * conj(ar[6; 3 4])
-                        ϵs[pos] = max(ϵs[pos], abs(1 - abs(v)))
-                    end
-                    @timeit timeroutput "update_AC" begin
-                        ψ.AC[pos] = (al, complex(c))
-                        ψ.AC[pos + 1] = (complex(c), _transpose_front(ar))
-                    end
+                # left-to-right
+                for pos in fwd
+                    ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos] =
+                        local_update!(
+                        pos, Val(:right),
+                        ψ, H, alg, envs,
+                        ϵ_global, ϵ_truncs[pos], decay_rates[pos],
+                        timeroutput
+                    )
+                    ϵ_global = maximum(ϵ_locals)
                 end
 
-                # right to left sweep
-                for pos in (length(ψ) - 2):-1:1
-                    local ac2, newA2center, al, c, ar
-                    @timeit timeroutput "AC2_eigsolve" begin
-                        @plansor ac2[-1 -2; -3 -4] := ψ.AL[pos][-1 -2; 1] * ψ.AC[pos + 1][1 -4; -3]
-                        Hac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs)
-                        _, newA2center = fixedpoint(Hac2, ac2, :SR, alg_eigsolve)
-                    end
-                    @timeit timeroutput "svd_trunc" begin
-                        al, c, ar = svd_trunc!(newA2center; trunc = alg.trscheme, alg = alg.alg_svd)
-                        normalize!(c)
-                        v = @plansor ac2[1 2; 3 4] * conj(al[1 2; 5]) * conj(c[5; 6]) * conj(ar[6; 3 4])
-                        ϵs[pos] = max(ϵs[pos], abs(1 - abs(v)))
-                    end
-                    @timeit timeroutput "update_AC" begin
-                        ψ.AC[pos + 1] = (complex(c), _transpose_front(ar))
-                        ψ.AC[pos] = (al, complex(c))
-                    end
+                # right-to-left
+                for pos in bwd
+                    ψ, ϵ_locals[pos], ϵ_truncs[pos], decay_rates[pos] =
+                        local_update!(
+                        pos, Val(:left),
+                        ψ, H, alg, envs,
+                        ϵ_global, ϵ_truncs[pos], decay_rates[pos],
+                        timeroutput
+                    )
+                    ϵ_global = maximum(ϵ_locals)
                 end
             end
+            ϵ_global = maximum(ϵ_locals)
 
-            ϵ = maximum(ϵs)
             ψ, envs = @timeit timeroutput "finalize" alg.finalize(
                 iter, ψ, H, envs
             )::Tuple{typeof(ψ), typeof(envs)}
 
-            if ϵ <= alg.tol
+            # Truncation-aware convergence: the Galerkin gradient cannot drop below the level set by
+            # the discarded weight, so a truncating scheme converges once `ϵ_global` reaches the
+            # truncation error rather than the (unreachable) bare `tol`. With no truncation
+            # (`ϵ_truncs .= 0`, e.g. single-site/QR gauge) this reduces to the plain `ϵ_global ≤ tol`.
+            if ϵ_global <= max(alg.tol, maximum(ϵ_truncs))
                 @infov 4 timeroutput
-                @infov 2 logfinish!(log, iter, ϵ, expectation_value(ψ, H, envs))
+                @infov 2 logfinish!(log, iter, ϵ_global, expectation_value(ψ, H, envs))
                 break
             end
             if iter == alg.maxiter
                 @infov 4 timeroutput
-                @warnv 1 logcancel!(log, iter, ϵ, expectation_value(ψ, H, envs))
+                @warnv 1 logcancel!(log, iter, ϵ_global, expectation_value(ψ, H, envs))
             else
-                @infov 3 logiter!(log, iter, ϵ, expectation_value(ψ, H, envs))
+                @infov 3 logiter!(log, iter, ϵ_global, expectation_value(ψ, H, envs))
             end
         end
     end
-    return ψ, envs, ϵ
+    return ψ, envs, ϵ_global
 end
 
 function find_groundstate(ψ, H, alg::Union{DMRG, DMRG2}, envs...; kwargs...)
diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl
index ed83e664a..66143b3bf 100644
--- a/src/algorithms/groundstate/idmrg.jl
+++ b/src/algorithms/groundstate/idmrg.jl
@@ -110,7 +110,7 @@ function find_groundstate(mps, operator, alg::alg_type, envs = environments(mps,
             @infov 3 logiter!(log, it.iter, ϵ, ΔE)
         end
 
-        alg_gauge = updatetol(alg.alg_gauge, it.state.iter, it.state.ϵ)
+        alg_gauge = adapt_solver(alg.alg_gauge; iter = it.state.iter, g_global = it.state.ϵ)
         ψ′ = InfiniteMPS(it.state.mps.AR; alg_gauge.tol, alg_gauge.maxiter)
         envs = recalculate!(it.state.envs, ψ′, it.state.operator, ψ′)
         return ψ′, envs, it.state.ϵ
@@ -151,7 +151,7 @@ end
 function localupdate_step!(
         it::IterativeSolver{<:IDMRG}, state
     )
-    alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ)
+    alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = state.iter, g_global = state.ϵ)
     return _localupdate_sweep_idmrg!(
         state.mps, state.operator, state.envs, alg_eigsolve, state.timeroutput,
     )
@@ -160,7 +160,7 @@ end
 function localupdate_step!(
         it::IterativeSolver{<:IDMRG2}, state
     )
-    alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ)
+    alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = state.iter, g_global = state.ϵ)
     return _localupdate_sweep_idmrg2!(
         state.mps, state.operator, state.envs, alg_eigsolve,
         it.trscheme, it.alg_svd, state.timeroutput,
diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl
index aa1c413cb..1f1cf51f7 100644
--- a/src/algorithms/groundstate/vumps.jl
+++ b/src/algorithms/groundstate/vumps.jl
@@ -63,7 +63,7 @@ function dominant_eigsolve(
 
     mps = copy(mps)
     ϵ = calc_galerkin(mps, operator, mps, envs)
-    alg_environments = updatetol(alg.alg_environments, iter, ϵ)
+    alg_environments = adapt_solver(alg.alg_environments; iter, g_global = ϵ)
     recalculate!(envs, mps, operator, mps, alg_environments; timeroutput)
 
     state = VUMPSState(mps, operator, envs, iter, ϵ, which, timeroutput)
@@ -118,8 +118,8 @@ end
 function localupdate_step!(
         it::IterativeSolver{<:VUMPS}, state, scheduler = Defaults.scheduler[]
     )
-    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
-    alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ)
+    alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ)
+    alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = state.iter, g_global = state.ϵ)
     alg_orth = alg_gauge.alg_orth
 
     mps = state.mps
@@ -186,7 +186,7 @@ function _localupdate_vumps_step!(
 end
 
 function gauge_step!(it::IterativeSolver{<:VUMPS}, state, ACs::AbstractVector)
-    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
+    alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ)
     mps = gaugefix!(
         state.mps, ACs, state.mps.C[end];
         order = :R, timeroutput = state.timeroutput, alg_gauge...,
@@ -195,11 +195,11 @@ function gauge_step!(it::IterativeSolver{<:VUMPS}, state, ACs::AbstractVector)
     return mps
 end
 function gauge_step!(it::IterativeSolver{<:VUMPS}, state, ACs::AbstractMatrix)
-    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
+    alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ)
     return MultilineMPS(ACs, @view(state.mps.C[:, end]); alg_gauge.tol, alg_gauge.maxiter, alg_gauge.alg_orth)
 end
 
 function envs_step!(it::IterativeSolver{<:VUMPS}, state, mps)
-    alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ)
+    alg_environments = adapt_solver(it.alg_environments; iter = state.iter, g_global = state.ϵ)
     return recalculate!(state.envs, mps, state.operator, mps, alg_environments; state.timeroutput)
 end
diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl
index 614e9f6bf..208da3664 100644
--- a/src/algorithms/statmech/idmrg.jl
+++ b/src/algorithms/statmech/idmrg.jl
@@ -8,7 +8,7 @@ function leading_boundary(
     LoggingExtras.withlevel(; alg.verbosity) do
         @infov 2 loginit!(log, ϵ, expectation_value(ψ, operator, envs))
         for outer iter in 1:(alg.maxiter)
-            alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ)
+            alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = ϵ)
             C_current = ψ.C[:, 0]
 
             # left to right sweep
@@ -54,7 +54,7 @@ function leading_boundary(
         end
     end
 
-    alg_gauge = updatetol(alg.alg_gauge, iter, ϵ)
+    alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
     ψ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
 
     recalculate!(envs, ψ, operator, ψ)
@@ -72,7 +72,7 @@ function leading_boundary(
     LoggingExtras.withlevel(; alg.verbosity) do
         @infov 2 loginit!(log, ϵ)
         for outer iter in 1:(alg.maxiter)
-            alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ)
+            alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = ϵ)
             C_current = ψ.C[:, 0]
 
             # sweep from left to right
@@ -189,7 +189,7 @@ function leading_boundary(
         end
     end
 
-    alg_gauge = updatetol(alg.alg_gauge, iter, ϵ)
+    alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
     ψ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
 
     recalculate!(envs, ψ, operator, ψ)
diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl
index 9295f2b28..935d9111f 100644
--- a/src/algorithms/statmech/vomps.jl
+++ b/src/algorithms/statmech/vomps.jl
@@ -54,7 +54,7 @@ function dominant_eigsolve(
     log = IterLog("VOMPS")
     iter = 0
     ϵ = calc_galerkin(mps, operator, mps, envs)
-    alg_environments = updatetol(alg.alg_environments, iter, ϵ)
+    alg_environments = adapt_solver(alg.alg_environments; iter, g_global = ϵ)
     recalculate!(envs, mps, operator, mps, alg_environments)
 
     state = VOMPSState(mps, operator, envs, iter, ϵ)
@@ -100,7 +100,7 @@ end
 function localupdate_step!(
         it::IterativeSolver{<:VOMPS}, state, scheduler = Defaults.scheduler[]
     )
-    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
+    alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ)
     alg_orth = alg_gauge.alg_orth
     mps = state.mps
     ACs = similar(mps.AC)
@@ -135,15 +135,15 @@ function _localupdate_vomps_step!(
 end
 
 function gauge_step!(it::IterativeSolver{<:VOMPS}, state, ACs::AbstractVector)
-    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
+    alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ)
     return InfiniteMPS(ACs, state.mps.C[end]; alg_gauge...)
 end
 function gauge_step!(it::IterativeSolver{<:VOMPS}, state, ACs::AbstractMatrix)
-    alg_gauge = updatetol(it.alg_gauge, state.iter, state.ϵ)
+    alg_gauge = adapt_solver(it.alg_gauge; iter = state.iter, g_global = state.ϵ)
     return MultilineMPS(ACs, @view(state.mps.C[:, end]); alg_gauge...)
 end
 
 function envs_step!(it::IterativeSolver{<:VOMPS}, state, mps)
-    alg_environments = updatetol(it.alg_environments, state.iter, state.ϵ)
+    alg_environments = adapt_solver(it.alg_environments; iter = state.iter, g_global = state.ϵ)
     return recalculate!(state.envs, mps, state.operator, mps, alg_environments)
 end
diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index acea13fe3..d9a763402 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -170,7 +170,11 @@ end
 TensorKit.space(W::JordanMPOTensor) = space(W.tensors)
 Base.eltype(::Type{T}) where {T <: JordanMPOTensor} = eltype(fieldtype(T, :tensors))
 
-function Base.getproperty(W::JordanMPOTensor, sym::Symbol)
+Base.size(W::JordanMPOTensor) = size(W.tensors)
+Base.size(W::JordanMPOTensor, i::Int) = size(W.tensors, i)
+BlockTensorKit.eachspace(W::JordanMPOTensor) = eachspace(W.tensors)
+
+Base.@constprop :aggressive function Base.getproperty(W::JordanMPOTensor, sym::Symbol)
     sym === :A && return _jordan_A(W)
     sym === :B && return _jordan_B(W)
     sym === :C && return _jordan_C(W)
diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl
index 6f65dbbe5..f38989b9c 100644
--- a/src/states/finitemps.jl
+++ b/src/states/finitemps.jl
@@ -343,8 +343,15 @@ Base.@propagate_inbounds function Base.getindex(ψ::FiniteMPS, i::Int)
     end
 end
 
-# TODO: check where gauge center is to determine efficient kind
-AC2(psi::FiniteMPS, site::Int) = psi.AC[site] * _transpose_tail(psi.AR[site + 1])
+function AC2(psi::FiniteMPS, site::Int; kind = :ACAR)
+    if kind == :ACAR
+        return psi.AC[site] * _transpose_tail(psi.AR[site + 1])
+    elseif kind == :ALAC
+        return psi.AL[site] * _transpose_tail(psi.AC[site + 1])
+    else
+        throw(ArgumentError("Invalid kind: $kind"))
+    end
+end
 
 f_if_not_missing(f, x) = ismissing(x) ? x : f(x)
 _copy_if_not_missing(x) = f_if_not_missing(copy, x)
diff --git a/src/states/ortho.jl b/src/states/ortho.jl
index 0b4cbd6ca..d7901de08 100644
--- a/src/states/ortho.jl
+++ b/src/states/ortho.jl
@@ -246,7 +246,7 @@ end
 function gauge_eigsolve_step!(it::IterativeSolver{LeftCanonical}, state)
     (; AL, C, A, iter, ϵ) = state
     if iter ≥ it.eig_miniter
-        alg_eigsolve = updatetol(it.alg_eigsolve, 1, ϵ^2)
+        alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = 1, g_global = ϵ^2)
         _, vec = fixedpoint(flip(TransferMatrix(A, AL)), C[end], :LM, alg_eigsolve)
         _, C[end] = left_orth!(vec; alg = it.alg_orth)
     end
@@ -308,7 +308,7 @@ end
 function gauge_eigsolve_step!(it::IterativeSolver{RightCanonical}, state)
     (; AR, C, A, iter, ϵ) = state
     if iter ≥ it.eig_miniter
-        alg_eigsolve = updatetol(it.alg_eigsolve, 1, ϵ^2)
+        alg_eigsolve = adapt_solver(it.alg_eigsolve; iter = 1, g_global = ϵ^2)
         _, vec = fixedpoint(TransferMatrix(A, AR), C[end], :LM, alg_eigsolve)
         C[end], _ = right_orth!(vec; alg = it.alg_orth)
     end
diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl
index bb2338340..efa394b8a 100644
--- a/src/states/orthoview.jl
+++ b/src/states/orthoview.jl
@@ -60,7 +60,7 @@ function Base.getindex(v::CView{<:FiniteMPS, E}, i::Int)::E where {E}
         end
 
         for j in Iterators.reverse((i + 1):center)
-            v.parent.Cs[j], v.parent.ARs[j] = right_gauge(v.parent.ACs[j])
+            v.parent.Cs[j], v.parent.ARs[j], _ = right_gauge(v.parent.ACs[j])
             if j != i + 1 # last AC not needed
                 v.parent.ACs[j - 1] = _mul_tail(v.parent.ALs[j - 1], v.parent.Cs[j])
             end
@@ -75,7 +75,7 @@ function Base.getindex(v::CView{<:FiniteMPS, E}, i::Int)::E where {E}
         end
 
         for j in center:i
-            v.parent.ALs[j], v.parent.Cs[j + 1] = left_gauge(v.parent.ACs[j])
+            v.parent.ALs[j], v.parent.Cs[j + 1], _ = left_gauge(v.parent.ACs[j])
             if j != i # last AC not needed
                 v.parent.ACs[j + 1] = _mul_front(v.parent.Cs[j + 1], v.parent.ARs[j + 1])
             end
@@ -88,9 +88,9 @@ end
 function Base.setindex!(v::CView{<:FiniteMPS}, vec, i::Int)
     if ismissing(v.parent.Cs[i + 1])
         if !ismissing(v.parent.ALs[i])
-            v.parent.Cs[i + 1], v.parent.ARs[i + 1] = right_gauge(v.parent.AC[i + 1])
+            v.parent.Cs[i + 1], v.parent.ARs[i + 1], _ = right_gauge(v.parent.AC[i + 1])
         else
-            v.parent.ALs[i], v.parent.Cs[i + 1] = left_gauge(v.parent.AC[i])
+            v.parent.ALs[i], v.parent.Cs[i + 1], _ = left_gauge(v.parent.AC[i])
         end
     end
 
@@ -223,8 +223,8 @@ end
 # Gauging routines
 # ----------------
 @doc """
-    left_gauge(AC, [alg]) -> AL, C
-    right_gauge(AC, [alg]) -> C, AR
+    left_gauge(AC, [alg]) -> AL, C, ϵ
+    right_gauge(AC, [alg]) -> C, AR, ϵ
 
 Factor an updated center MPS tensor `AC` into left- or right-canonical form, `AC ≈ AL * C`
 (left, with `AL` left-isometric) or `AC ≈ C * AR` (right, with `AR` right-isometric), without
@@ -233,28 +233,38 @@ standard MPS-tensor form.
 
 `alg` selects the factorization and defaults to a (positive) QR/LQ center-move that preserves the virtual bond.
 Passing a [`TruncatedAlgorithm`](@extref MatrixAlgebraKit.TruncatedAlgorithm) instead performs a truncated SVD may shrink the bond.
+
+Also returns the truncation error `ϵ`: the 2-norm of the discarded singular values from the
+truncated SVD, or `0` for a norm-preserving QR/LQ gauge.
 """
 left_gauge
 @doc (@doc left_gauge) right_gauge
 
 left_gauge(AC) = left_gauge(AC, Defaults.alg_orth())
-left_gauge(AC, alg) = left_orth(AC; alg)
-left_gauge(AC, alg::MatrixAlgebraKit.TruncatedAlgorithm) =
-    left_orth(AC; alg = MatrixAlgebraKit.LeftOrthViaSVD(alg))
+function left_gauge(AC, alg)
+    AL, C = left_orth(AC; alg)
+    return AL, C, zero(real(scalartype(AC)))
+end
+function left_gauge(AC, alg::MatrixAlgebraKit.TruncatedAlgorithm)
+    U, S, Vᴴ, ϵ = svd_trunc(AC, alg)
+    C = LinearAlgebra.lmul!(S, Vᴴ) # C = S * Vᴴ, matching `LeftOrthViaSVD`
+    return U, C, ϵ
+end
 
 right_gauge(AC) = right_gauge(AC, Defaults.alg_orth())
 function right_gauge(AC, alg)
     C, AR = right_orth(_transpose_tail(AC); alg)
-    return C, _transpose_front(AR)
+    return C, _transpose_front(AR), zero(real(scalartype(AC)))
 end
 function right_gauge(AC, alg::MatrixAlgebraKit.TruncatedAlgorithm)
-    C, AR = right_orth(_transpose_tail(AC); alg = MatrixAlgebraKit.RightOrthViaSVD(alg))
-    return C, _transpose_front(AR)
+    U, S, Vᴴ, ϵ = svd_trunc(_transpose_tail(AC), alg)
+    C = LinearAlgebra.rmul!(U, S) # C = U * S, matching `RightOrthViaSVD`
+    return C, _transpose_front(Vᴴ), ϵ
 end
 
 @doc """
-    left_gauge!(ψ, pos, AC, [alg]; normalize = false) -> ψ
-    right_gauge!(ψ, pos, AC, [alg]; normalize = false) -> ψ
+    left_gauge!(ψ, pos, AC, [alg]; normalize = false) -> ψ, ϵ
+    right_gauge!(ψ, pos, AC, [alg]; normalize = false) -> ψ, ϵ
 
 Gauge an updated center tensor `AC` at site `pos` and install it into `ψ` in one step: factor
 `AC` with [`left_gauge`](@ref) / [`right_gauge`](@ref) and write the canonical tensors back,
@@ -264,19 +274,69 @@ be a [`TruncatedAlgorithm`](@extref MatrixAlgebraKit.TruncatedAlgorithm) to trun
 
 By default the factors are installed as-is. Pass `normalize = true` to renormalize the bond
 tensor, so `ψ` stays normalized after a local update that changed its norm.
+
+Also returns the truncation error `ϵ` of the factorization: the 2-norm of the discarded singular
+values from a truncated SVD gauge, or `0` for a norm-preserving QR gauge.
 """
 left_gauge!
 @doc (@doc left_gauge!) right_gauge!
 
 function left_gauge!(ψ::AbstractFiniteMPS, pos::Int, AC, alg = Defaults.alg_orth(); normalize::Bool = false)
-    AL, C = left_gauge(AC, alg)
+    AL, C, ϵ = left_gauge(AC, alg)
     normalize && normalize!(C)
     ψ.AC[pos] = (AL, C)
-    return ψ
+    return ψ, ϵ
 end
 function right_gauge!(ψ::AbstractFiniteMPS, pos::Int, AC, alg = Defaults.alg_orth(); normalize::Bool = false)
-    C, AR = right_gauge(AC, alg)
+    C, AR, ϵ = right_gauge(AC, alg)
     normalize && normalize!(C)
     ψ.AC[pos] = (C, AR)
-    return ψ
+    return ψ, ϵ
+end
+
+@doc """
+    gauge!(ψ, pos, direction, AC, [alg]; normalize = false) -> ψ, ϵ
+
+Direction-dispatching wrapper around [`left_gauge!`](@ref) / [`right_gauge!`](@ref): gauge an
+updated center tensor `AC` at site `pos` and install it into `ψ`, shifting the gauge center past
+`pos` to the right for `direction = Val(:right)` and to the left for `Val(:left)`. `alg` and
+`normalize` are forwarded unchanged, and the truncation error `ϵ` is returned as-is.
+"""
+function gauge!(ψ::AbstractFiniteMPS, pos::Int, ::Val{Dir}, AC, alg = Defaults.alg_orth(); kwargs...) where {Dir}
+    Dir === :right && return left_gauge!(ψ, pos, AC, alg; kwargs...)
+    Dir === :left && return right_gauge!(ψ, pos, AC, alg; kwargs...)
+    return throw(ArgumentError(lazy"invalid direction `$Dir`"))
+end
+
+@doc """
+    gauge2!(ψ, pos, direction, AC2, alg; normalize = false) -> ψ, ϵ
+
+Two-site analogue of [`gauge!`](@ref): factor an updated two-site center tensor `AC2` spanning
+sites `pos` and `pos+1` with the truncated SVD `alg` and install the resulting canonical tensors
+into `ψ` in one step, shifting the gauge center past the bond.
+(To the right for `direction = Val(:right)`, to the left for `Val(:left)`).
+
+Returns the truncation error `ϵ`, the 2-norm of the discarded singular values. Pass
+`normalize = true` to renormalize the bond tensor, so `ψ` stays normalized after a local update
+that changed its norm.
+"""
+function gauge2!(
+        ψ::AbstractFiniteMPS, pos::Int, ::Val{Dir}, AC2, alg;
+        normalize::Bool = false
+    ) where {Dir}
+    al, c, ar, ϵ = svd_trunc!(AC2, alg)
+    normalize && normalize!(c)
+    # the SVD center `c` (singular values) is real-valued; promote it to the state's scalartype so
+    # the installed tensors stay consistent, without needlessly forcing a real-valued state complex.
+    C = scalartype(ψ) <: Complex ? complex(c) : c
+    if Dir === :right
+        ψ.AC[pos] = (al, C)
+        ψ.AC[pos + 1] = (C, _transpose_front(ar))
+    elseif Dir === :left
+        ψ.AC[pos + 1] = (C, _transpose_front(ar))
+        ψ.AC[pos] = (al, C)
+    else
+        throw(ArgumentError(lazy"invalid direction `$Dir`"))
+    end
+    return ψ, ϵ
 end
diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl
index ddfd02afe..dbc221420 100644
--- a/src/utility/defaults.jl
+++ b/src/utility/defaults.jl
@@ -7,7 +7,7 @@ module Defaults
 
 import KrylovKit: GMRES, Arnoldi, Lanczos
 using OhMyThreads
-using ..MPSKit: DynamicTol
+using ..MPSKit: DynamicTol, AdaptiveKrylov
 using TensorKit: TensorKit
 using MatrixAlgebraKit: DefaultAlgorithm, Householder
 
@@ -54,8 +54,9 @@ end
 function alg_eigsolve(;
         ishermitian = true, tol = tol, maxiter = maxiter, verbosity = 0,
         eager = true, krylovdim = krylovdim, dynamic_tols = dynamic_tols, tol_min = tol_min,
-        tol_max = tol_max, tol_factor = eigs_tolfactor
+        tol_max = tol_max, tol_factor = eigs_tolfactor, adaptive = false, adaptive_kwargs...
     )
+    adaptive && return AdaptiveKrylov(; ishermitian, adaptive_kwargs...)
     alg = ishermitian ? Lanczos(; tol, maxiter, eager, krylovdim, verbosity) :
         Arnoldi(; tol, maxiter, eager, krylovdim, verbosity)
     return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg
diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl
index 7a67a9a01..fe52d3188 100644
--- a/src/utility/dynamictols.jl
+++ b/src/utility/dynamictols.jl
@@ -4,16 +4,35 @@ import ..MPSKit: Algorithm
 using Accessors
 using DocStringExtensions
 using MatrixAlgebraKit: DefaultAlgorithm
+using KrylovKit: KrylovKit, Lanczos, Arnoldi
 
-export updatetol, DynamicTol
+export adapt_solver, DynamicTol, AdaptiveKrylov
+# deprecated, kept for backwards compatibility (see the bottom of this file)
+export updatetol
+
+# ============================================================
+# Unified solver-adaptation interface
+# ============================================================
 
 @doc """
-    updatetol(alg, iter, ϵ)
+    adapt_solver(alg; iter, decay_rate, g_local, g_global, eps_trunc)
 
-Update the tolerance of the algorithm `alg` based on the current iteration `iter` and the current error `ϵ`.
-""" updatetol
+Resolve a (possibly adaptive) solver specification `alg` into a concrete solver, using
+whichever adaptation signals are supplied as keyword arguments. Numeric signals default to
+`0.0` (and `iter` to `1`) — i.e. "no contribution" — so each caller passes only what it has in
+scope:
 
-updatetol(alg, iter::Integer, ϵ::Real) = alg
+  - `iter`       — outer iteration count (drives the `1/√iter` damping)
+  - `decay_rate` — measured per-matvec contraction factor of the previous solve
+  - `g_local`    — local (per-bond) gradient / Galerkin norm
+  - `g_global`   — global (sweep-wide) gradient / convergence-error scalar
+  - `eps_trunc`  — local truncation error
+
+The generic fallback returns `alg` unchanged (for plain `Lanczos`/`Arnoldi`/… solvers, which
+are not adapted). `DynamicTol` and `AdaptiveKrylov` provide the adaptive implementations.
+""" adapt_solver
+
+adapt_solver(alg; kwargs...) = alg
 
 # Wrapper for dynamic tolerance adjustment
 # ----------------------------------------
@@ -21,13 +40,15 @@ updatetol(alg, iter::Integer, ϵ::Real) = alg
 """
 $(TYPEDEF)
 
-Algorithm wrapper with dynamically adjusted tolerances.
+Algorithm wrapper with dynamically adjusted tolerances. Only the wrapped solver's tolerance is
+retuned; its Krylov budget (if any) is left fixed — this is the simpler counterpart to
+[`AdaptiveKrylov`](@ref).
 
 ## Fields
 
 $(TYPEDFIELDS)
 
-See also [`updatetol`](@ref).
+See also [`adapt_solver`](@ref).
 """
 struct DynamicTol{A} <: Algorithm
     "parent algorithm"
@@ -39,12 +60,10 @@ struct DynamicTol{A} <: Algorithm
     "maximal value of the dynamic tolerance"
     tol_max::Float64
 
-    "tolerance factor for updating relative to current algorithm error"
+    "tolerance factor for updating relative to the current (global) gradient norm"
     tol_factor::Float64
 
-    function DynamicTol(
-            alg::A, tol_min::Real, tol_max::Real, tol_factor::Real
-        ) where {A}
+    function DynamicTol(alg::A, tol_min::Real, tol_max::Real, tol_factor::Real) where {A}
         0 <= tol_min <= tol_max ||
             throw(ArgumentError("tol_min must be between 0 and tol_max"))
         return new{A}(alg, tol_min, tol_max, tol_factor)
@@ -55,17 +74,16 @@ function DynamicTol(alg; tol_min = 1.0e-6, tol_max = 1.0e-2, tol_factor = 0.1)
 end
 
 """
-    updatetol(alg::DynamicTol, iter, ϵ)
+    adapt_solver(alg::DynamicTol; iter, g_global, ...)
+
+Tighten only the wrapped solver's tolerance (its Krylov budget, if any, is left fixed), from the
+global gradient / convergence-error scalar `g_global`, damped by the iteration count:
 
-Update the tolerance of the algorithm `alg` based on the current iteration `iter` and the current error `ϵ`,
-where the new tolerance is given by
-    
-    new_tol = clamp(ϵ * alg.tol_factor / sqrt(iter), alg.tol_min, alg.tol_max)
+    tol = clamp(tol_factor·g_global / √iter, tol_min, tol_max)
 """
-function updatetol(alg::DynamicTol, iter::Integer, ϵ::Real)
-    iter = max(iter, one(iter))
-    new_tol = clamp(ϵ * alg.tol_factor / sqrt(iter), alg.tol_min, alg.tol_max)
-    return _updatetol(alg.alg, new_tol)
+function adapt_solver(alg::DynamicTol; iter::Integer = 1, g_global::Real = 0.0, kwargs...)
+    tol = clamp(alg.tol_factor * g_global / sqrt(max(iter, 1)), alg.tol_min, alg.tol_max)
+    return _updatetol(alg.alg, tol)
 end
 
 # default implementation with Accessors.jl, but can be hooked into
@@ -77,4 +95,152 @@ function _updatetol(alg::DefaultAlgorithm, tol::Real)
     return DefaultAlgorithm(; kwargs...)
 end
 
+# Set several solver parameters at once (tol, krylovdim, maxiter). Used by the adaptive
+# controller, which retunes more than just the tolerance. `Lanczos`/`Arnoldi` expose these
+# as plain fields, so `Accessors.setproperties` rebuilds the immutable struct in one shot.
+function _set_eigsolve_params(alg::Union{Lanczos, Arnoldi}; kwargs...)
+    return Accessors.setproperties(alg, NamedTuple(kwargs))
+end
+
+# ============================================================
+# Adaptive Krylov controller (per-bond, stateful)
+# ============================================================
+
+"""
+$(TYPEDEF)
+
+Adaptive controller for local eigensolvers in DMRG-like algorithms.
+
+This controller attempts to balance minimizing the number of local function applications with
+minimizing the total number of iterations in order to obtain the globally fastest convergence.
+This is driven by the local and global gradient norm, the truncation error and the measured
+decay rate of previous iterations in an attempt to obtain fast convergence for gapped systems
+while avoiding stagnation for gapless ones.
+
+## Fields
+
+$(TYPEDFIELDS)
+
+See also [`adapt_solver`](@ref).
+"""
+struct AdaptiveKrylov{T, O <: KrylovKit.Orthogonalizer} <: Algorithm
+    "orthogonalizer passed to the instantiated `Lanczos`/`Arnoldi`"
+    orth::O
+
+    "minimal Krylov subspace dimension (conservative/cold-start default)"
+    krylovdim_min::Int
+    "maximal Krylov subspace dimension"
+    krylovdim_max::Int
+
+    "minimal number of restart iterations"
+    iter_min::Int
+    "maximal number of restart iterations"
+    iter_max::Int
+
+    "lower bound on the dynamically chosen tolerance"
+    tol_min::Float64
+    "upper bound on the dynamically chosen tolerance"
+    tol_max::Float64
+    "factor on the local truncation error, sets the inner-solve tolerance floor"
+    truncation_factor::Float64
+    "factor on the global gradient norm, sets the convergence-driven tolerance"
+    tol_factor::Float64
+
+    "whether to use the eager (early-stopping) solver mode"
+    eager::Bool
+    "verbosity of the instantiated eigensolver"
+    verbosity::Int
+
+    "hermitian flag (`Val(true)` → `Lanczos`, `Val(false)` → `Arnoldi`)"
+    ishermitian::Val{T}
+
+    function AdaptiveKrylov{T, O}(
+            orth::O, krylovdim_min, krylovdim_max, iter_min, iter_max,
+            tol_min, tol_max, truncation_factor, tol_factor, eager, verbosity, ishermitian::Val{T}
+        ) where {T, O <: KrylovKit.Orthogonalizer}
+        0 < krylovdim_min <= krylovdim_max ||
+            throw(ArgumentError("need 0 < krylovdim_min ≤ krylovdim_max, got ($krylovdim_min, $krylovdim_max)"))
+        1 <= iter_min <= iter_max ||
+            throw(ArgumentError("need 1 ≤ iter_min ≤ iter_max, got ($iter_min, $iter_max)"))
+        0 <= tol_min <= tol_max ||
+            throw(ArgumentError("need 0 ≤ tol_min ≤ tol_max, got ($tol_min, $tol_max)"))
+        truncation_factor >= 0 ||
+            throw(ArgumentError("truncation_factor must be non-negative, got $truncation_factor"))
+        tol_factor >= 0 ||
+            throw(ArgumentError("tol_factor must be non-negative, got $tol_factor"))
+        return new{T, O}(
+            orth, krylovdim_min, krylovdim_max, iter_min, iter_max,
+            tol_min, tol_max, truncation_factor, tol_factor, eager, verbosity, ishermitian
+        )
+    end
+end
+
+function AdaptiveKrylov(;
+        ishermitian::Bool = true, orth::KrylovKit.Orthogonalizer = KrylovKit.KrylovDefaults.orth,
+        krylovdim_min::Int = 3, krylovdim_max::Int = 20,
+        iter_min::Int = 1, iter_max::Int = 1,
+        tol_min::Real = 1.0e-12, tol_max::Real = 1.0e-2,
+        truncation_factor::Real = 1.0e-1, tol_factor::Real = 1.0e-1,
+        eager::Bool = true, verbosity::Int = 0
+    )
+    return AdaptiveKrylov{ishermitian, typeof(orth)}(
+        orth, krylovdim_min, krylovdim_max, iter_min, iter_max,
+        tol_min, tol_max, truncation_factor, tol_factor, eager, verbosity, Val(ishermitian)
+    )
+end
+
+"""
+    adapt_solver(alg::AdaptiveKrylov; decay_rate, g_local, g_global, eps_trunc, ...)
+
+Build a concrete `Lanczos`/`Arnoldi` for the current site. The target tolerance is the same as
+[`DynamicTol`](@ref)'s per-bond tolerance, but the Krylov `krylovdim`/`maxiter` are additionally
+predicted from the measured `decay_rate` and the local gradient norm `g_local`. A `decay_rate`
+of `0` (the default) signals a cold start and falls back to the minimal budget.
+"""
+function adapt_solver(
+        alg::AdaptiveKrylov{T}; decay_rate::Real = 0.0, g_local::Real = 0.0,
+        g_global::Real = 0.0, eps_trunc::Real = 0.0, kwargs...
+    ) where {T}
+    # 1. target tolerance: inexact inner solve depending on outer convergence,
+    #    never below the truncation error we already incur.
+    tol = clamp(
+        max(alg.truncation_factor * eps_trunc, alg.tol_factor * g_global), alg.tol_min, alg.tol_max
+    )
+
+    # the measured decay rate must be a genuine contraction in (0, 1); a non-positive or ≥ 1
+    # rate signals an uninitialized or hard/stalling bond → fall back to the largest budget.
+    ρ = clamp(decay_rate, 0.0, 1.0)
+    if !(0 < ρ < 1)
+        krylovdim = iszero(ρ) ? alg.krylovdim_min : alg.krylovdim_max
+        maxiter = iszero(ρ) ? alg.iter_min : alg.iter_max
+        return (T ? Lanczos : Arnoldi)(; alg.orth, krylovdim, maxiter, tol, alg.eager, alg.verbosity)
+    end
+
+    # 2. predicted matvecs from the *measured* decay and local convergence
+    #    estimate through Lanczos bounds: going from g_local to tol in factors of ρ
+    R = max(g_local, alg.tol_min) / tol
+    nmatvecs = round(Int, log(R) / log(inv(ρ)))
+
+    # 3. krylovdim and maxiter from matvec budget:
+    #    prioritize krylovdim and compensate with maxiter
+    #    thick restarts keep 3/5 of krylovdim, i.e. ~2/5 new vectors per cycle
+    krylovdim = clamp(nmatvecs, alg.krylovdim_min, alg.krylovdim_max)
+    maxiter = clamp(cld(5nmatvecs, 2krylovdim) - 3, alg.iter_min, alg.iter_max)
+
+    return (T ? Lanczos : Arnoldi)(; alg.orth, krylovdim, maxiter, tol, alg.eager, alg.verbosity)
+end
+
+# ============================================================
+# Deprecated entry points (subsumed by `adapt_solver`)
+# ============================================================
+
+Base.@deprecate updatetol(alg, iter::Integer, ϵ::Real) adapt_solver(alg; iter = iter, g_global = ϵ)
+
+@doc """
+    updatetol(alg, iter, ϵ)
+
+Deprecated: superseded by [`adapt_solver`](@ref). Equivalent to
+`adapt_solver(alg; iter, g_global = ϵ)`.
+""" updatetol
+
 end