diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 7efbec564..09b1d49d1 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -27,6 +27,16 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 
 ### Changed
 
+- Internal: the finite `DMRG`/`DMRG2` ground-state search and approximation, `DynamicalDMRG`,
+  multiline `IDMRG`/`IDMRG2` approximation and `leading_boundary`, and `time_evolve` now run as
+  iterators whose iterations are a forward and a backward half-sweep (or a single time step),
+  sharing one structure across solvers. Finite `TDVP`/`TDVP2` steps are composed of per-site local
+  updates. Results and logging are unchanged. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
+- The `IDMRG`/`IDMRG2` convergence measure (`bondresidual`) is computed by one shared,
+  lower-allocation routine. For multiline `IDMRG` with several rows, the changes of the rows are
+  now summed, as `IDMRG2` already did, instead of combined in a 2-norm, which makes the stopping
+  criterion marginally stricter. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
+
 ### Deprecated
 
 ### Removed
diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl
index 465b40f1b..ed17ce1a0 100644
--- a/src/algorithms/approximate/fvomps.jl
+++ b/src/algorithms/approximate/fvomps.jl
@@ -1,74 +1,82 @@
-function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environments(ψ, _environment_args(Oϕ)...))
+# Internal state of the finite DMRG/DMRG2 approximation, where `ϵ` is the largest relative
+# local change of the last sweep
+struct ApproximateState{S, O, E, A}
+    mps::S
+    operator::O
+    envs::E
+    iter::Int
+    ϵ::Float64
+    allocator::A
+end
+
+function approximate!(
+        ψ::AbstractFiniteMPS, Oϕ, alg::Union{DMRG, DMRG2},
+        envs = environments(ψ, _environment_args(Oϕ)...)
+    )
     allocator = default_allocator(ψ, SerialScheduler())
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
-    log = IterLog("DMRG2")
+    log = IterLog(alg)
+    it = IterativeSolver(alg, ApproximateState(ψ, Oϕ, envs, 0, 2 * alg.tol, allocator))
 
     with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            ϵ = 0.0
-            for pos in [1:(length(ψ) - 1); (length(ψ) - 2):-1:1]
-                AC2′ = AC2_projection(pos, ψ, Oϕ, envs; alg.backend, allocator)
-                al, c, ar, = svd_trunc!(AC2′, inner_alg_gauge(alg))
-
-                AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
-                ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2))
-
-                ψ.AC[pos] = (al, complex(c))
-                ψ.AC[pos + 1] = (complex(c), _transpose_front(ar))
-            end
-
-            # finalize
-            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ), typeof(envs)}
-
+        @log_initialization loginit!(log, it.ϵ)
+        for (_, _, ϵ) in Iterators.take(it, alg.maxiter)
             if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
+                @log_convergence logfinish!(log, it.iter, ϵ)
                 break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
+            elseif it.iter == alg.maxiter
+                @log_nonconvergence logcancel!(log, it.iter, ϵ)
             else
-                @log_iteration logiter!(log, iter, ϵ)
+                @log_iteration logiter!(log, it.iter, ϵ)
             end
         end
     end
 
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, localchange = ϵ, numiter = iter)
+    state = it.state
+    info = AlgorithmInfo(; converged = state.ϵ <= alg.tol, localchange = state.ϵ, numiter = state.iter)
+    return state.mps, state.envs, info
 end
 
-function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs = environments(ψ, _environment_args(Oϕ)...))
-    allocator = default_allocator(ψ, SerialScheduler())
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
-    log = IterLog("DMRG")
-
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            ϵ = 0.0
-            for pos in [1:(length(ψ) - 1); length(ψ):-1:2]
-                AC′ = AC_projection(pos, ψ, Oϕ, envs; alg.backend, allocator)
-                AC = ψ.AC[pos]
-                ϵ = max(ϵ, norm(AC′ - AC) / norm(AC′))
+function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::ApproximateState)
+    iter = state.iter + 1
+    state = ApproximateState(
+        state.mps, state.operator, state.envs, state.iter, 0.0, state.allocator
+    )
+    state = sweep!(it, state, Val(:right), iter)
+    state = sweep!(it, state, Val(:left), iter)
+    ψ, envs = it.finalize(
+        iter, state.mps, state.operator, state.envs
+    )::Tuple{typeof(state.mps), typeof(state.envs)}
+    it.state = ApproximateState(ψ, state.operator, envs, iter, state.ϵ, state.allocator)
+    return (ψ, envs, state.ϵ), it.state
+end
 
-                ψ.AC[pos] = AC′
-            end
+function sweep!(it::IterativeSolver{<:DMRG2}, state::ApproximateState, direction, iter)
+    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    sites = direction === Val(:right) ? fwd : bwd
+    ψ, Oϕ, envs, ϵ = state.mps, state.operator, state.envs, state.ϵ
+    for pos in sites
+        AC2′ = AC2_projection(pos, ψ, Oϕ, envs; it.alg.backend, state.allocator)
+        al, c, ar, = svd_trunc!(AC2′, inner_alg_gauge(it.alg))
 
-            # finalize
-            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ), typeof(envs)}
+        AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
+        ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2))
 
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
-        end
+        ψ.AC[pos] = (al, complex(c))
+        ψ.AC[pos + 1] = (complex(c), _transpose_front(ar))
     end
+    return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator)
+end
 
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, localchange = ϵ, numiter = iter)
+function sweep!(it::IterativeSolver{<:DMRG}, state::ApproximateState, direction, iter)
+    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    sites = direction === Val(:right) ? fwd : bwd
+    ψ, Oϕ, envs, ϵ = state.mps, state.operator, state.envs, state.ϵ
+    for pos in sites
+        AC′ = AC_projection(pos, ψ, Oϕ, envs; it.alg.backend, state.allocator)
+        AC = ψ.AC[pos]
+        ϵ = max(ϵ, norm(AC′ - AC) / norm(AC′))
+
+        ψ.AC[pos] = AC′
+    end
+    return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator)
 end
diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl
index 8b8079379..9fe55fbf7 100644
--- a/src/algorithms/approximate/idmrg.jl
+++ b/src/algorithms/approximate/idmrg.jl
@@ -1,199 +1,182 @@
 function approximate!(
-        ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO, <:MultilineMPS}, alg::IDMRG,
-        envs = environments(ψ, toapprox...)
+        ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO, <:MultilineMPS},
+        alg::Union{IDMRG, IDMRG2}, envs = environments(ψ, toapprox...)
     )
     allocator = default_allocator(ψ, SerialScheduler())
-    log = IterLog("IDMRG")
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
+    alg isa IDMRG2 && width(ψ) < 2 && throw(ArgumentError("unit cell should be >= 2"))
+    log = IterLog(alg)
+    ϵ_truncs = alg isa IDMRG2 ?
+        PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ))) : nothing
+    state = IDMRGState(ψ, toapprox, envs, 0, 2 * alg.tol, ϵ_truncs, nothing, NoTimerOutput(), allocator)
+    it = IterativeSolver(alg, state)
 
     with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            C_current = ψ.C[:, 0]
-
-            # left to right sweep
-            for col in 1:width(ψ)
-                for row in 1:size(ψ, 1)
-                    ψ.AC[row + 1, col] = AC_projection(
-                        CartesianIndex(row, col), ψ, toapprox, envs;
-                        alg.backend, allocator
-                    )
-                    normalize!(ψ.AC[row + 1, col])
-                    ψ.AL[row + 1, col], ψ.C[row + 1, col] = left_orth!(ψ.AC[row + 1, col])
-                end
-                transfer_leftenv!(envs, ψ, toapprox, col + 1)
-            end
-
-            # right to left sweep
-            for col in reverse(1:width(ψ))
-                for row in 1:size(ψ, 1)
-                    ψ.AC[row + 1, col] = AC_projection(
-                        CartesianIndex(row, col), ψ, toapprox, envs;
-                        alg.backend, allocator
-                    )
-                    normalize!(ψ.AC[row + 1, col])
-                    ψ.C[row + 1, col - 1], temp = right_orth!(_transpose_tail(ψ.AC[row + 1, col]))
-                    ψ.AR[row + 1, col] = _transpose_front(temp)
-                end
-                transfer_rightenv!(envs, ψ, toapprox, col - 1)
-            end
-            normalize!(envs, ψ, toapprox)
-
-            ϵ = norm(C_current - ψ.C[:, 0])
-
+        @log_initialization loginit!(log, it.ϵ)
+        for _ in Iterators.take(it, alg.maxiter)
+            ϵ = it.ϵ
             if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
+                @log_convergence logfinish!(log, it.iter, ϵ)
                 break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
+            elseif it.iter == alg.maxiter
+                @log_nonconvergence logcancel!(log, it.iter, ϵ)
             else
-                @log_iteration logiter!(log, iter, ϵ)
+                @log_iteration logiter!(log, it.iter, ϵ)
             end
         end
     end
 
+    (; iter, ϵ) = it.state
+
     # TODO: immediately compute in-place
     alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
-    ψ′ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
+    ψ′ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
     copy!(ψ, ψ′) # ensure output destination is unchanged
 
     recalculate!(envs, ψ, toapprox)
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, numiter = iter)
+    info = AlgorithmInfo(;
+        converged = ϵ <= alg.tol, bondresidual = ϵ,
+        truncation_errors = isnothing(ϵ_truncs) ? nothing : parent(ϵ_truncs), numiter = iter,
+    )
+    return ψ, envs, info
 end
 
-function approximate!(
-        ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO, <:MultilineMPS},
-        alg::IDMRG2, envs = environments(ψ, toapprox...)
+function Base.iterate(it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state::IDMRGState{<:MultilineMPS, <:Tuple})
+    iter = state.iter + 1
+    C_current = state.mps.C[:, 0]
+    state = sweep!(it, state, Val(:right), iter)
+    state = sweep!(it, state, Val(:left), iter)
+    ϵ = bond_change(C_current, state.mps.C[:, 0])
+    it.state = IDMRGState(
+        state.mps, state.operator, state.envs, iter, ϵ, state.truncation_errors, state.energy,
+        state.timeroutput, state.allocator,
     )
-    allocator = default_allocator(ψ, SerialScheduler())
-    width(ψ) < 2 && throw(ArgumentError("unit cell should be >= 2"))
-    ϵ::Float64 = 2 * alg.tol
-    log = IterLog("IDMRG2")
-    O, ϕ = toapprox
-    iter = 0
-    ϵ_truncs = PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ)))
+    return (it.state.mps, it.state.envs, it.state.ϵ), it.state
+end
 
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            C_current = ψ.C[:, 0]
-
-            # sweep from left to right
-            for site in 1:(width(ψ) - 1)
-                for row in 1:size(ψ, 1)
-                    AC2′ = AC2_projection(
-                        CartesianIndex(row, site), ψ, toapprox, envs;
-                        kind = :ACAR, alg.backend, allocator
-                    )
-                    al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
-                    normalize!(c)
-
-                    ψ.AL[row + 1, site] = al
-                    ψ.C[row + 1, site] = complex(c)
-                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
-                    ψ.AC[row + 1, site + 1] = _transpose_front(c * ar)
-                end
-
-                transfer_leftenv!(envs, ψ, toapprox, site + 1)
-                transfer_rightenv!(envs, ψ, toapprox, site)
-            end
+function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGState{<:MultilineMPS, <:Tuple}, ::Val{:right}, iter)
+    alg, ψ, toapprox, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    for col in 1:width(ψ)
+        for row in 1:size(ψ, 1)
+            ψ.AC[row + 1, col] = AC_projection(
+                CartesianIndex(row, col), ψ, toapprox, envs;
+                alg.backend, allocator
+            )
+            normalize!(ψ.AC[row + 1, col])
+            ψ.AL[row + 1, col], ψ.C[row + 1, col] = left_orth!(ψ.AC[row + 1, col])
+        end
+        transfer_leftenv!(envs, ψ, toapprox, col + 1)
+    end
+    return state
+end
 
-            # update the edge
-            ψ.AL[:, end] .= ψ.AC[:, end] ./ ψ.C[:, end]
-            ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
-            for row in 1:size(ψ, 1)
-                AC2′ = AC2_projection(
-                    CartesianIndex(row, width(ψ)), ψ, toapprox, envs;
-                    kind = :ALAC, alg.backend, allocator
-                )
-                al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
-                normalize!(c)
-
-                ψ.AL[row + 1, end] = al
-                ψ.C[row + 1, end] = complex(c)
-                ψ.AR[row + 1, 1] = _transpose_front(ar)
-
-                ψ.AC[row + 1, end] = _mul_tail(al, c)
-                ψ.AC[row + 1, 1] = _transpose_front(c * ar)
-                ψ.AL[row + 1, 1] = ψ.AC[row + 1, 1] / ψ.C[row + 1, 1]
+function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGState{<:MultilineMPS, <:Tuple}, ::Val{:left}, iter)
+    alg, ψ, toapprox, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    for col in reverse(1:width(ψ))
+        for row in 1:size(ψ, 1)
+            ψ.AC[row + 1, col] = AC_projection(
+                CartesianIndex(row, col), ψ, toapprox, envs;
+                alg.backend, allocator
+            )
+            normalize!(ψ.AC[row + 1, col])
+            ψ.C[row + 1, col - 1], temp = right_orth!(_transpose_tail(ψ.AC[row + 1, col]))
+            ψ.AR[row + 1, col] = _transpose_front(temp)
+        end
+        transfer_rightenv!(envs, ψ, toapprox, col - 1)
+    end
+    normalize!(envs, ψ, toapprox)
+    return state
+end
 
-            end
-            # update environments
-            transfer_leftenv!(envs, ψ, toapprox, 1)
-            transfer_rightenv!(envs, ψ, toapprox, 0)
-
-            normalize!(envs, ψ, toapprox)
-
-            # sweep from right to left
-            for site in reverse(1:(width(ψ) - 1))
-                for row in 1:size(ψ, 1)
-                    AC2′ = AC2_projection(
-                        CartesianIndex(row, site), ψ, toapprox, envs;
-                        kind = :ALAC, alg.backend, allocator
-                    )
-                    al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
-                    normalize!(c)
-
-                    ψ.AL[row + 1, site] = al
-                    ψ.C[row + 1, site] = complex(c)
-                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
-                end
-
-                transfer_leftenv!(envs, ψ, toapprox, site + 1)
-                transfer_rightenv!(envs, ψ, toapprox, site)
-            end
+function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGState{<:MultilineMPS, <:Tuple}, ::Val{:right}, iter)
+    alg, ψ, toapprox, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    ϵ_truncs = state.truncation_errors
+    for site in 1:(width(ψ) - 1)
+        for row in 1:size(ψ, 1)
+            AC2′ = AC2_projection(
+                CartesianIndex(row, site), ψ, toapprox, envs;
+                kind = :ACAR, alg.backend, allocator
+            )
+            al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
+            normalize!(c)
+
+            ψ.AL[row + 1, site] = al
+            ψ.C[row + 1, site] = complex(c)
+            ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+            ψ.AC[row + 1, site + 1] = _transpose_front(c * ar)
+        end
 
-            # update the edge
-            ψ.AC[:, end] .= _mul_front.(ψ.C[:, end - 1], ψ.AR[:, end])
-            ψ.AR[:, 1] .= _transpose_front.(ψ.C[:, end] .\ _transpose_tail.(ψ.AC[:, 1]))
-            for row in 1:size(ψ, 1)
-                AC2′ = AC2_projection(
-                    CartesianIndex(row, 0), ψ, toapprox, envs;
-                    kind = :ACAR, alg.backend, allocator
-                )
-                al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
-                normalize!(c)
-
-                ψ.AL[row + 1, end] = al
-                ψ.C[row + 1, end] = complex(c)
-                ψ.AR[row + 1, 1] = _transpose_front(ar)
-
-                ψ.AR[row + 1, end] = _transpose_front(ψ.C[row + 1, end - 1] \ _transpose_tail(al * c))
-                ψ.AC[row + 1, 1] = _transpose_front(c * ar)
-            end
-            transfer_leftenv!(envs, ψ, toapprox, 1)
-            transfer_rightenv!(envs, ψ, toapprox, 0)
+        transfer_leftenv!(envs, ψ, toapprox, site + 1)
+        transfer_rightenv!(envs, ψ, toapprox, site)
+    end
 
-            normalize!(envs, ψ, toapprox)
+    # update the edge
+    ψ.AL[:, end] .= ψ.AC[:, end] ./ ψ.C[:, end]
+    ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
+    for row in 1:size(ψ, 1)
+        AC2′ = AC2_projection(
+            CartesianIndex(row, width(ψ)), ψ, toapprox, envs;
+            kind = :ALAC, alg.backend, allocator
+        )
+        al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
+        normalize!(c)
+
+        ψ.AL[row + 1, end] = al
+        ψ.C[row + 1, end] = complex(c)
+        ψ.AR[row + 1, 1] = _transpose_front(ar)
+
+        ψ.AC[row + 1, end] = _mul_tail(al, c)
+        ψ.AC[row + 1, 1] = _transpose_front(c * ar)
+        ψ.AL[row + 1, 1] = ψ.AC[row + 1, 1] / ψ.C[row + 1, 1]
+    end
+    transfer_leftenv!(envs, ψ, toapprox, 1)
+    transfer_rightenv!(envs, ψ, toapprox, 0)
 
-            # update error
-            ϵ = sum(zip(C_current, ψ.C[:, 0])) do (c1, c2)
-                smallest = infimum(_firstspace(c1), _firstspace(c2))
-                e1 = isometry(_firstspace(c1), smallest)
-                e2 = isometry(_firstspace(c2), smallest)
-                return norm(e2' * c2 * e2 - e1' * c1 * e1)
-            end
+    normalize!(envs, ψ, toapprox)
+    return state
+end
 
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
+function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGState{<:MultilineMPS, <:Tuple}, ::Val{:left}, iter)
+    alg, ψ, toapprox, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    ϵ_truncs = state.truncation_errors
+    for site in reverse(1:(width(ψ) - 1))
+        for row in 1:size(ψ, 1)
+            AC2′ = AC2_projection(
+                CartesianIndex(row, site), ψ, toapprox, envs;
+                kind = :ALAC, alg.backend, allocator
+            )
+            al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
+            normalize!(c)
+
+            ψ.AL[row + 1, site] = al
+            ψ.C[row + 1, site] = complex(c)
+            ψ.AR[row + 1, site + 1] = _transpose_front(ar)
         end
+
+        transfer_leftenv!(envs, ψ, toapprox, site + 1)
+        transfer_rightenv!(envs, ψ, toapprox, site)
     end
 
-    # TODO: immediately compute in-place
-    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
+    # update the edge
+    ψ.AC[:, end] .= _mul_front.(ψ.C[:, end - 1], ψ.AR[:, end])
+    ψ.AR[:, 1] .= _transpose_front.(ψ.C[:, end] .\ _transpose_tail.(ψ.AC[:, 1]))
+    for row in 1:size(ψ, 1)
+        AC2′ = AC2_projection(
+            CartesianIndex(row, 0), ψ, toapprox, envs;
+            kind = :ACAR, alg.backend, allocator
+        )
+        al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(AC2′; trunc = alg.trunc, alg = alg.alg_svd)
+        normalize!(c)
+
+        ψ.AL[row + 1, end] = al
+        ψ.C[row + 1, end] = complex(c)
+        ψ.AR[row + 1, 1] = _transpose_front(ar)
+
+        ψ.AR[row + 1, end] = _transpose_front(ψ.C[row + 1, end - 1] \ _transpose_tail(al * c))
+        ψ.AC[row + 1, 1] = _transpose_front(c * ar)
+    end
+    transfer_leftenv!(envs, ψ, toapprox, 1)
+    transfer_rightenv!(envs, ψ, toapprox, 0)
 
-    recalculate!(envs, ψ, toapprox)
-    info = AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, truncation_errors = parent(ϵ_truncs), numiter = iter)
-    return ψ, envs, info
+    normalize!(envs, ψ, toapprox)
+    return state
 end
diff --git a/src/algorithms/approximate/vomps.jl b/src/algorithms/approximate/vomps.jl
index 475913e08..361c846cc 100644
--- a/src/algorithms/approximate/vomps.jl
+++ b/src/algorithms/approximate/vomps.jl
@@ -21,7 +21,7 @@ function approximate(
 end
 
 function _approximate_vomps(mps, toapprox, alg::VOMPS, envs)
-    log = IterLog("VOMPS")
+    log = IterLog(alg)
     iter = 0
     ϵ = calc_galerkin(mps, toapprox..., envs; alg.backend)
     alg_environments = adapt_solver(alg.alg_environments; iter, g_global = ϵ)
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index b2400fdfa..ad0ade0a1 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -295,78 +295,101 @@ function find_groundstate!(
     return find_groundstate_sweep!(ψ, H, alg, envs, allocator, timeroutput)
 end
 
-function find_groundstate_sweep!(
-        ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs, allocator, timeroutput
-    )
-    name = string(nameof(typeof(alg)))
-    log = IterLog(name)
+# A solver iteration is one complete forward/backward sweep. The vectors retain the
+# local history used by the adaptive solvers between sweeps.
+struct DMRGState{S, O, E, R, V, D, T, A}
+    mps::S
+    operator::O
+    envs::E
+    iter::Int
+    ϵ::R
+    local_errors::V
+    truncation_errors::V
+    decay_rates::D
+    timeroutput::T
+    allocator::A
+end
 
+function DMRGState(ψ, H, alg::Union{DMRG, DMRG2}, envs, allocator, 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, ψ)
-    iter = 0
+    local_errors = ones(Tr, n)
+    return DMRGState(
+        ψ, H, envs, 0, maximum(local_errors), local_errors, zeros(Tr, n),
+        zeros(n), timeroutput, allocator,
+    )
+end
 
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ_global, expectation_value(ψ, H, envs))
-        for outer iter in 1:(alg.maxiter)
-            @timeit timeroutput "sweep" begin
-                # 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],
-                        iter, timeroutput, allocator
-                    )
-                    ϵ_global = maximum(ϵ_locals)
-                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],
-                        iter, timeroutput, allocator
-                    )
-                    ϵ_global = maximum(ϵ_locals)
-                end
-            end
-            ϵ_global = maximum(ϵ_locals)
+function sweep!(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::DMRGState, direction, iter)
+    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    sites = direction === Val(:right) ? fwd : bwd
+    ψ, ϵ = state.mps, state.ϵ
+    for pos in sites
+        ψ, state.local_errors[pos], state.truncation_errors[pos], state.decay_rates[pos] =
+            local_update!(
+            pos, direction, ψ, state.operator, it.alg, state.envs,
+            ϵ, state.truncation_errors[pos], state.decay_rates[pos],
+            iter, state.timeroutput, state.allocator,
+        )
+        # The next local solver uses the errors accumulated so far in this sweep.
+        ϵ = maximum(state.local_errors)
+    end
+    return DMRGState(
+        ψ, state.operator, state.envs, state.iter, ϵ, state.local_errors,
+        state.truncation_errors, state.decay_rates, state.timeroutput, state.allocator,
+    )
+end
+
+function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::DMRGState)
+    iter = state.iter + 1
+    timeroutput = state.timeroutput
+    state = @timeit timeroutput "sweep" begin
+        state = sweep!(it, state, Val(:right), iter)
+        sweep!(it, state, Val(:left), iter)
+    end
+    ψ, envs = @timeit timeroutput "finalize" it.finalize(
+        iter, state.mps, state.operator, state.envs
+    )::Tuple{typeof(state.mps), typeof(state.envs)}
+    it.state = DMRGState(
+        ψ, state.operator, envs, iter, state.ϵ, state.local_errors,
+        state.truncation_errors, state.decay_rates, timeroutput, state.allocator,
+    )
+    return (ψ, envs, state.ϵ), it.state
+end
 
-            ψ, envs = @timeit timeroutput "finalize" alg.finalize(
-                iter, ψ, H, envs
-            )::Tuple{typeof(ψ), typeof(envs)}
+# Truncation sets the attainable floor of the Galerkin error.
+sweep_converged(alg, state::DMRGState) =
+    state.ϵ <= max(alg.tol, maximum(state.truncation_errors))
 
-            # 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))
+function find_groundstate_sweep!(
+        ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs, allocator, timeroutput
+    )
+    log = IterLog(alg)
+    it = IterativeSolver(alg, DMRGState(ψ, H, alg, envs, allocator, timeroutput))
+
+    with_verbosity(; alg.verbosity) do
+        @log_initialization loginit!(log, it.ϵ, expectation_value(ψ, H, envs))
+        for (ψ, envs, ϵ) in Iterators.take(it, alg.maxiter)
+            if sweep_converged(alg, it.state)
                 @info TimerReport(timeroutput) _group = :mpskit_timing
-                @log_convergence logfinish!(log, iter, ϵ_global, expectation_value(ψ, H, envs))
+                @log_convergence logfinish!(log, it.iter, ϵ, expectation_value(ψ, H, envs))
                 break
-            end
-            if iter == alg.maxiter
+            elseif it.iter == alg.maxiter
                 @info TimerReport(timeroutput) _group = :mpskit_timing
-                @log_nonconvergence logcancel!(log, iter, ϵ_global, expectation_value(ψ, H, envs))
+                @log_nonconvergence logcancel!(log, it.iter, ϵ, expectation_value(ψ, H, envs))
             else
-                @log_iteration logiter!(log, iter, ϵ_global, expectation_value(ψ, H, envs))
+                @log_iteration logiter!(log, it.iter, ϵ, expectation_value(ψ, H, envs))
             end
         end
     end
 
+    state = it.state
     info = AlgorithmInfo(;
-        converged = ϵ_global <= max(alg.tol, maximum(ϵ_truncs)), galerkin = ϵ_global,
-        truncation_errors = _bond_truncation_errors(alg, ϵ_truncs), numiter = iter
+        converged = sweep_converged(alg, state), galerkin = state.ϵ,
+        truncation_errors = _bond_truncation_errors(alg, state.truncation_errors),
+        numiter = state.iter,
     )
-    return ψ, envs, info
+    return state.mps, state.envs, info
 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 a0c314148..214183bc4 100644
--- a/src/algorithms/groundstate/idmrg.jl
+++ b/src/algorithms/groundstate/idmrg.jl
@@ -70,14 +70,17 @@ Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_bounda
     backend::B = Defaults.backend()
 end
 
-# Internal state of the IDMRG algorithm
+# Internal state of the IDMRG/IDMRG2 algorithms, shared by the ground state search
+# (`mps::InfiniteMPS`), the leading boundary search (`mps::MultilineMPS`) and approximation
+# (`mps::MultilineMPS`, `operator::Tuple`). `ϵ` is the change of the center bond tensor over the
+# last sweep, and `energy` is only tracked by the ground state search.
 struct IDMRGState{S, O, E, V, T, TO, A}
     mps::S
     operator::O
     envs::E
     iter::Int
     ϵ::Float64 # TODO: Could be any <:Real
-    truncation_errors::V # per bond, of the most recent sweep only
+    truncation_errors::V
     energy::T
     timeroutput::TO
     allocator::A
@@ -99,7 +102,7 @@ end
 function _find_groundstate_idmrg(mps, operator, alg::alg_type, envs) where {alg_type <: Union{<:IDMRG, <:IDMRG2}}
     (length(mps) ≤ 1 && alg isa IDMRG2) && throw(ArgumentError("unit cell should be >= 2"))
     name = alg isa IDMRG ? "IDMRG" : "IDMRG2"
-    log = IterLog(name)
+    log = IterLog(alg)
     timeroutput = alg.verbosity > 3 ? TimerOutput(name) : NoTimerOutput()
     mps = copy(mps)
     iter = 0
@@ -147,24 +150,14 @@ function _find_groundstate_idmrg(mps, operator, alg::alg_type, envs) where {alg_
 end
 
 function Base.iterate(
-        it::IterativeSolver{alg_type}, state::IDMRGState{<:Any, <:Any, <:Any, <:Any, T} = it.state
+        it::IterativeSolver{alg_type}, state::IDMRGState{<:InfiniteMPS, <:Any, <:Any, <:Any, T}
     ) where {alg_type <: Union{<:IDMRG, <:IDMRG2}, T}
     timeroutput = state.timeroutput
     ϵ_truncs = zero(state.truncation_errors) # fresh each sweep, filled by the sweep itself
     mps, envs, C_old, E_new = @timeit timeroutput "localupdate" localupdate_step!(it, state, ϵ_truncs)
 
     # error criterion
-    C = mps.C[0]
-    space_C_old = _firstspace(C_old)
-    space_C = _firstspace(C)
-    if space_C != space_C_old
-        smallest = infimum(space_C_old, space_C)
-        e1 = isometry(space_C_old, smallest)
-        e2 = isometry(space_C, smallest)
-        ϵ = norm(e2' * C * e2 - e1' * C_old * e1)
-    else
-        ϵ = norm(C - C_old)
-    end
+    ϵ = bond_change(C_old, mps.C[0])
 
     # New energy
     ΔE = (E_new - state.energy) / 2
diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl
index 9e7b144ef..2d1d8fd5e 100644
--- a/src/algorithms/groundstate/vumps.jl
+++ b/src/algorithms/groundstate/vumps.jl
@@ -67,7 +67,7 @@ function dominant_eigsolve(
         operator, mps, alg::VUMPS, envs = environments(mps, operator, mps, alg.alg_environments);
         which
     )
-    log = IterLog("VUMPS")
+    log = IterLog(alg)
     timeroutput = alg.verbosity > 3 ? TimerOutput("VUMPS") : NoTimerOutput()
     iter = 0
 
@@ -106,7 +106,7 @@ function dominant_eigsolve(
     return result
 end
 
-function Base.iterate(it::IterativeSolver{<:VUMPS}, state = it.state)
+function Base.iterate(it::IterativeSolver{<:VUMPS}, state::VUMPSState)
     timeroutput = state.timeroutput
     ACs = @timeit timeroutput "localupdate (parallel)" localupdate_step!(it, state)
     mps = @timeit timeroutput "gauge" gauge_step!(it, state, ACs)
diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl
index bd72f892c..5450e5bbe 100644
--- a/src/algorithms/propagator/corvector.jl
+++ b/src/algorithms/propagator/corvector.jl
@@ -38,6 +38,8 @@ Used as the `algorithm` argument of [`propagator`](@ref).
     backend::B = Defaults.backend()
 end
 
+IterativeLoggers.IterLog(::DynamicalDMRG) = IterLog("DDMRG")
+
 """
     propagator(ψ₀::AbstractFiniteMPS, z::Number, H::MPOHamiltonian, alg::DynamicalDMRG; init = copy(ψ₀)) -> (g, ψ)
 
@@ -69,6 +71,58 @@ Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac
 """
 struct NaiveInvert <: DDMRG_Flavour end
 
+# Internal state of the dynamical DMRG sweeps, which update `mps` in place towards the
+# propagator applied to `target`; `ϵ` is the largest change of a center tensor over the last sweep
+struct DDMRGState{S, T, Z, O, E, A}
+    mps::S
+    target::T
+    z::Z
+    operator::O
+    envs::E
+    iter::Int
+    ϵ::Float64
+    allocator::A
+end
+
+# dynamical DMRG sweeps over the sites like single-site DMRG
+_sweep_ranges(::DynamicalDMRG, ψ) = (1:(length(ψ) - 1), length(ψ):-1:2)
+
+function Base.iterate(it::IterativeSolver{<:DynamicalDMRG}, state::DDMRGState)
+    iter = state.iter + 1
+    state = DDMRGState(
+        state.mps, state.target, state.z, state.operator, state.envs, state.iter, 0.0,
+        state.allocator,
+    )
+    state = sweep!(it, state, Val(:right), iter)
+    state = sweep!(it, state, Val(:left), iter)
+    it.state = DDMRGState(
+        state.mps, state.target, state.z, state.operator, state.envs, iter, state.ϵ,
+        state.allocator,
+    )
+    return (state.mps, state.ϵ), it.state
+end
+
+function _propagator_sweeps!(alg::DynamicalDMRG, state::DDMRGState)
+    log = IterLog(alg)
+    it = IterativeSolver(alg, state)
+
+    with_verbosity(; alg.verbosity) do
+        @log_initialization loginit!(log, it.ϵ)
+        for (_, ϵ) in Iterators.take(it, alg.maxiter)
+            if ϵ <= alg.tol
+                @log_convergence logfinish!(log, it.iter, ϵ)
+                break
+            elseif it.iter == alg.maxiter
+                @log_nonconvergence logcancel!(log, it.iter, ϵ)
+            else
+                @log_iteration logiter!(log, it.iter, ϵ)
+            end
+        end
+    end
+
+    return it.state
+end
+
 function propagator(
         A::AbstractFiniteMPS, z::Number, H,
         alg::DynamicalDMRG{NaiveInvert}; init = copy(A)
@@ -77,41 +131,39 @@ function propagator(
     h_envs = environments(init, H, init) # environments for h
     mixedenvs = environments(init, A) # environments for <init | A>
 
-    ϵ = 2 * alg.tol
-    log = IterLog("DDMRG")
+    state = DDMRGState(init, A, z, H, (h_envs, mixedenvs), 0, 2 * alg.tol, allocator)
+    _propagator_sweeps!(alg, state)
 
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for iter in 1:(alg.maxiter)
-            ϵ = 0.0
+    return dot(A, init), init
+end
 
-            for i in [1:(length(A) - 1); length(A):-1:2]
-                tos = AC_projection(i, init, A, mixedenvs; alg.backend, allocator)
+function sweep!(
+        it::IterativeSolver{<:DynamicalDMRG{NaiveInvert}}, state::DDMRGState, direction, iter
+    )
+    alg = it.alg
+    (; mps, target, z, operator, allocator) = state
+    init, A, H = mps, target, operator
+    h_envs, mixedenvs = state.envs
+    fwd, bwd = _sweep_ranges(alg, A)
+    ϵ = state.ϵ
 
-                H_AC = AC_hamiltonian(i, init, H, init, h_envs; alg.backend, allocator)
-                AC = init.AC[i]
-                AC′, convhist = linsolve(H_AC, -tos, AC, alg.solver, -z, one(z))
+    for i in (direction === Val(:right) ? fwd : bwd)
+        tos = AC_projection(i, init, A, mixedenvs; alg.backend, allocator)
 
-                ϵ = max(ϵ, norm(AC′ - AC))
-                init.AC[i] = AC′
+        H_AC = AC_hamiltonian(i, init, H, init, h_envs; alg.backend, allocator)
+        AC = init.AC[i]
+        AC′, convhist = linsolve(H_AC, -tos, AC, alg.solver, -z, one(z))
 
-                convhist.converged == 0 &&
-                    @warn "propagator ($i) failed to converge: normres = $(convhist.normres)"
-            end
+        ϵ = max(ϵ, norm(AC′ - AC))
+        init.AC[i] = AC′
 
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
-        end
+        convhist.converged == 0 &&
+            @warn "propagator ($i) failed to converge: normres = $(convhist.normres)"
     end
 
-    return dot(A, init), init
+    return DDMRGState(
+        init, A, z, operator, state.envs, state.iter, ϵ, allocator,
+    )
 end
 
 """
@@ -154,39 +206,8 @@ function propagator(
     H2, envs2 = squaredenvs(init, H, envs1) # environments for h^2
     mixedenvs = environments(init, A) # environments for <init | A>
 
-    ϵ = 2 * alg.tol
-    log = IterLog("DDMRG")
-
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for iter in 1:(alg.maxiter)
-            ϵ = 0.0
-
-            for i in [1:(length(A) - 1); length(A):-1:2]
-                tos = AC_projection(i, init, A, mixedenvs; alg.backend, allocator)
-                H1_AC = AC_hamiltonian(i, init, H, init, envs1; alg.backend, allocator)
-                H2_AC = AC_hamiltonian(i, init, H2, init, envs2; alg.backend, allocator)
-                H_AC = LinearCombination((H1_AC, H2_AC), (-2 * ω, 1))
-                AC′, convhist = linsolve(H_AC, -η * tos, init.AC[i], alg.solver, abs2(z), 1)
-
-                ϵ = max(ϵ, norm(AC′ - init.AC[i]))
-                init.AC[i] = AC′
-
-                convhist.converged == 0 &&
-                    @warn "propagator ($i) failed to converge: normres $(convhist.normres)"
-            end
-
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
-        end
-    end
+    state = DDMRGState(init, A, z, (H, H2), (envs1, envs2, mixedenvs), 0, 2 * alg.tol, allocator)
+    _propagator_sweeps!(alg, state)
 
     a = dot(AC_projection(1, init, A, mixedenvs; alg.backend, allocator), init.AC[1])
     cb = leftenv(envs1, 1, A) * TransferMatrix(init.AL, H[1:length(A.AL)], A.AL)
@@ -200,6 +221,38 @@ function propagator(
     return v, init
 end
 
+function sweep!(
+        it::IterativeSolver{<:DynamicalDMRG{Jeckelmann}}, state::DDMRGState, direction, iter
+    )
+    alg = it.alg
+    (; mps, target, z, allocator) = state
+    init, A = mps, target
+    H, H2 = state.operator
+    envs1, envs2, mixedenvs = state.envs
+    ω = real(z)
+    η = imag(z)
+    fwd, bwd = _sweep_ranges(alg, A)
+    ϵ = state.ϵ
+
+    for i in (direction === Val(:right) ? fwd : bwd)
+        tos = AC_projection(i, init, A, mixedenvs; alg.backend, allocator)
+        H1_AC = AC_hamiltonian(i, init, H, init, envs1; alg.backend, allocator)
+        H2_AC = AC_hamiltonian(i, init, H2, init, envs2; alg.backend, allocator)
+        H_AC = LinearCombination((H1_AC, H2_AC), (-2 * ω, 1))
+        AC′, convhist = linsolve(H_AC, -η * tos, init.AC[i], alg.solver, abs2(z), 1)
+
+        ϵ = max(ϵ, norm(AC′ - init.AC[i]))
+        init.AC[i] = AC′
+
+        convhist.converged == 0 &&
+            @warn "propagator ($i) failed to converge: normres $(convhist.normres)"
+    end
+
+    return DDMRGState(
+        init, A, z, state.operator, state.envs, state.iter, ϵ, allocator,
+    )
+end
+
 function squaredenvs(
         state::AbstractFiniteMPS, H, envs = environments(state, H, state)
     )
diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl
index 8e5a9b0f0..fe5e0073f 100644
--- a/src/algorithms/statmech/idmrg.jl
+++ b/src/algorithms/statmech/idmrg.jl
@@ -1,201 +1,192 @@
 function leading_boundary(
-        ψ::InfiniteMultilineMPS, operator, alg::IDMRG, envs = environments(ψ, operator, ψ)
+        ψ::InfiniteMultilineMPS, operator, alg::Union{IDMRG, IDMRG2},
+        envs = environments(ψ, operator, ψ)
     )
     allocator = default_allocator(ψ, SerialScheduler())
-    log = IterLog("IDMRG")
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
+    alg isa IDMRG2 && width(ψ) < 2 && throw(ArgumentError("unit cell should be >= 2"))
+    log = IterLog(alg)
+    ϵ_truncs = alg isa IDMRG2 ?
+        PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ))) : nothing
+    state = IDMRGState(ψ, operator, envs, 0, 2 * alg.tol, ϵ_truncs, nothing, NoTimerOutput(), allocator)
+    it = IterativeSolver(alg, state)
 
     with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ, leading_eigenvalue(ψ, operator, envs))
-        for outer iter in 1:(alg.maxiter)
-            alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = ϵ)
-            C_current = ψ.C[:, 0]
-
-            # left to right sweep
-            for col in 1:width(ψ)
-                Hac = AC_hamiltonian(col, ψ, operator, ψ, envs; alg.backend, allocator)
-                _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve)
-
-                for row in 1:size(ψ, 1)
-                    ac = ψ.AC[row, col]
-                    (col == width(ψ)) && (ac = copy(ac)) # needed in next sweep
-                    ψ.AL[row, col], ψ.C[row, col] = left_orth!(ac)
-                end
-
-                transfer_leftenv!(envs, ψ, operator, ψ, col + 1)
-            end
-
-            # right to left sweep
-            for col in width(ψ):-1:1
-                Hac = AC_hamiltonian(col, ψ, operator, ψ, envs; alg.backend, allocator)
-                _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve)
-
-                for row in 1:size(ψ, 1)
-                    ψ.C[row, col - 1], temp = right_orth!(_transpose_tail(ψ.AC[row, col]; copy = true))
-                    ψ.AR[row, col] = _transpose_front(temp)
-                end
-
-                transfer_rightenv!(envs, ψ, operator, ψ, col - 1)
-            end
-
-            normalize!(envs, ψ, operator, ψ)
-
-            ϵ = norm(C_current - ψ.C[:, 0])
-
+        @log_initialization loginit!(log, it.ϵ, _boundary_objective(alg, ψ, operator, envs))
+        for (ψ, envs, ϵ) in Iterators.take(it, alg.maxiter)
             if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ, leading_eigenvalue(ψ, operator, envs))
+                @log_convergence logfinish!(log, it.iter, ϵ, _boundary_objective(alg, ψ, operator, envs))
                 break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ, leading_eigenvalue(ψ, operator, envs))
+            elseif it.iter == alg.maxiter
+                @log_nonconvergence logcancel!(log, it.iter, ϵ, _boundary_objective(alg, ψ, operator, envs))
             else
-                @log_iteration logiter!(log, iter, ϵ, leading_eigenvalue(ψ, operator, envs))
+                @log_iteration logiter!(log, it.iter, ϵ, _boundary_objective(alg, ψ, operator, envs))
             end
         end
     end
 
+    (; iter, ϵ) = it.state
     alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
-    ψ = MultilineMPS(map(x -> x, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
+    ψ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
 
     recalculate!(envs, ψ, operator, ψ)
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, numiter = iter)
+    info = AlgorithmInfo(;
+        converged = ϵ <= alg.tol, bondresidual = ϵ,
+        truncation_errors = isnothing(ϵ_truncs) ? nothing : parent(ϵ_truncs), numiter = iter,
+    )
+    return ψ, envs, info
 end
 
-function leading_boundary(
-        ψ::InfiniteMultilineMPS, operator, alg::IDMRG2, envs = environments(ψ, operator, ψ)
+# only single-site IDMRG reports the leading eigenvalue in its log
+_boundary_objective(::IDMRG, ψ, operator, envs) = leading_eigenvalue(ψ, operator, envs)
+_boundary_objective(::IDMRG2, ψ, operator, envs) = nothing
+
+function Base.iterate(it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state::IDMRGState{<:MultilineMPS})
+    iter = state.iter + 1
+    C_current = state.mps.C[:, 0]
+    state = sweep!(it, state, Val(:right), iter)
+    state = sweep!(it, state, Val(:left), iter)
+    ϵ = bond_change(C_current, state.mps.C[:, 0])
+    it.state = IDMRGState(
+        state.mps, state.operator, state.envs, iter, ϵ, state.truncation_errors, state.energy,
+        state.timeroutput, state.allocator,
     )
-    allocator = default_allocator(ψ, SerialScheduler())
-    width(ψ) < 2 && throw(ArgumentError("unit cell should be >= 2"))
-    ϵ::Float64 = 2 * alg.tol
-    log = IterLog("IDMRG2")
-    iter = 0
-    ϵ_truncs = PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ)))
+    return (it.state.mps, it.state.envs, it.state.ϵ), it.state
+end
 
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = ϵ)
-            C_current = ψ.C[:, 0]
-
-            # sweep from left to right
-            for site in 1:(width(ψ) - 1)
-                ac2 = AC2(ψ, site; kind = :ACAR)
-                h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
-                _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
-
-                for row in 1:size(ψ, 1)
-                    al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
-                    normalize!(c)
-
-                    ψ.AL[row + 1, site] = al
-                    ψ.C[row + 1, site] = complex(c)
-                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
-                    ψ.AC[row + 1, site + 1] = _transpose_front(c * ar)
-                end
-
-                transfer_leftenv!(envs, ψ, operator, ψ, site + 1)
-                transfer_rightenv!(envs, ψ, operator, ψ, site)
-            end
+function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGState{<:MultilineMPS}, ::Val{:right}, iter)
+    alg, ψ, operator, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = state.ϵ)
+    for col in 1:width(ψ)
+        Hac = AC_hamiltonian(col, ψ, operator, ψ, envs; alg.backend, allocator)
+        _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve)
+
+        for row in 1:size(ψ, 1)
+            ac = ψ.AC[row, col]
+            (col == width(ψ)) && (ac = copy(ac)) # needed in next sweep
+            ψ.AL[row, col], ψ.C[row, col] = left_orth!(ac)
+        end
 
-            normalize!(envs, ψ, operator, ψ)
+        transfer_leftenv!(envs, ψ, operator, ψ, col + 1)
+    end
+    return state
+end
 
-            # update the edge
-            site = width(ψ)
-            ψ.AL[:, end] .= ψ.AC[:, end] ./ ψ.C[:, end]
-            ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
-            ac2 = AC2(ψ, site; kind = :ALAC)
-            h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
-            _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGState{<:MultilineMPS}, ::Val{:left}, iter)
+    alg, ψ, operator, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = state.ϵ)
+    for col in width(ψ):-1:1
+        Hac = AC_hamiltonian(col, ψ, operator, ψ, envs; alg.backend, allocator)
+        _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve)
 
-            for row in 1:size(ψ, 1)
-                al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
-                normalize!(c)
+        for row in 1:size(ψ, 1)
+            ψ.C[row, col - 1], temp = right_orth!(_transpose_tail(ψ.AC[row, col]; copy = true))
+            ψ.AR[row, col] = _transpose_front(temp)
+        end
 
-                ψ.AL[row + 1, site] = al
-                ψ.C[row + 1, site] = complex(c)
-                ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+        transfer_rightenv!(envs, ψ, operator, ψ, col - 1)
+    end
+    normalize!(envs, ψ, operator, ψ)
+    return state
+end
 
-                ψ.AC[row + 1, site] = _mul_tail(al, c)
-                ψ.AC[row + 1, 1] = _transpose_front(c * ar)
-                ψ.AL[row + 1, 1] = ψ.AC[row + 1, 1] / ψ.C[row + 1, 1]
-            end
+function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGState{<:MultilineMPS}, ::Val{:right}, iter)
+    alg, ψ, operator, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = state.ϵ)
+    ϵ_truncs = state.truncation_errors
+    for site in 1:(width(ψ) - 1)
+        ac2 = AC2(ψ, site; kind = :ACAR)
+        h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
+        _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+
+        for row in 1:size(ψ, 1)
+            al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
+            normalize!(c)
+
+            ψ.AL[row + 1, site] = al
+            ψ.C[row + 1, site] = complex(c)
+            ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+            ψ.AC[row + 1, site + 1] = _transpose_front(c * ar)
+        end
 
-            # TODO: decide if we should compare at the half-sweep level?
-            # C_current = ψ.C[:, site]
+        transfer_leftenv!(envs, ψ, operator, ψ, site + 1)
+        transfer_rightenv!(envs, ψ, operator, ψ, site)
+    end
 
-            transfer_leftenv!(envs, ψ, operator, ψ, 1)
-            transfer_rightenv!(envs, ψ, operator, ψ, 0)
+    normalize!(envs, ψ, operator, ψ)
 
-            # sweep from right to left
-            for site in reverse(1:(width(ψ) - 1))
-                ac2 = AC2(ψ, site; kind = :ALAC)
-                h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
-                _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+    # update the edge
+    site = width(ψ)
+    ψ.AL[:, end] .= ψ.AC[:, end] ./ ψ.C[:, end]
+    ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
+    ac2 = AC2(ψ, site; kind = :ALAC)
+    h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
+    _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
 
-                for row in 1:size(ψ, 1)
-                    al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
-                    normalize!(c)
+    for row in 1:size(ψ, 1)
+        al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
+        normalize!(c)
 
-                    ψ.AL[row + 1, site] = al
-                    ψ.C[row + 1, site] = complex(c)
-                    ψ.AR[row + 1, site + 1] = _transpose_front(ar)
-                end
+        ψ.AL[row + 1, site] = al
+        ψ.C[row + 1, site] = complex(c)
+        ψ.AR[row + 1, site + 1] = _transpose_front(ar)
 
-                transfer_leftenv!(envs, ψ, operator, ψ, site + 1)
-                transfer_rightenv!(envs, ψ, operator, ψ, site)
-            end
+        ψ.AC[row + 1, site] = _mul_tail(al, c)
+        ψ.AC[row + 1, 1] = _transpose_front(c * ar)
+        ψ.AL[row + 1, 1] = ψ.AC[row + 1, 1] / ψ.C[row + 1, 1]
+    end
 
-            normalize!(envs, ψ, operator, ψ)
+    transfer_leftenv!(envs, ψ, operator, ψ, 1)
+    transfer_rightenv!(envs, ψ, operator, ψ, 0)
+    return state
+end
 
-            # update the edge
-            ψ.AC[:, end] .= _mul_front.(ψ.C[:, end - 1], ψ.AR[:, end])
-            ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
-            ψ.AR[:, 1] .= _transpose_front.(ψ.C[:, end] .\ _transpose_tail.(ψ.AC[:, 1]))
-            ac2 = AC2(ψ, 0; kind = :ACAR)
-            h = AC2_hamiltonian(0, ψ, operator, ψ, envs; alg.backend, allocator)
-            _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGState{<:MultilineMPS}, ::Val{:left}, iter)
+    alg, ψ, operator, envs, allocator = it.alg, state.mps, state.operator, state.envs, state.allocator
+    alg_eigsolve = adapt_solver(alg.alg_eigsolve; iter, g_global = state.ϵ)
+    ϵ_truncs = state.truncation_errors
+    for site in reverse(1:(width(ψ) - 1))
+        ac2 = AC2(ψ, site; kind = :ALAC)
+        h = AC2_hamiltonian(site, ψ, operator, ψ, envs; alg.backend, allocator)
+        _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
+
+        for row in 1:size(ψ, 1)
+            al, c, ar, ϵ_truncs[row + 1, site] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
+            normalize!(c)
+
+            ψ.AL[row + 1, site] = al
+            ψ.C[row + 1, site] = complex(c)
+            ψ.AR[row + 1, site + 1] = _transpose_front(ar)
+        end
 
-            for row in 1:size(ψ, 1)
-                al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
-                normalize!(c)
+        transfer_leftenv!(envs, ψ, operator, ψ, site + 1)
+        transfer_rightenv!(envs, ψ, operator, ψ, site)
+    end
 
-                ψ.AL[row + 1, end] = al
-                ψ.C[row + 1, end] = complex(c)
-                ψ.AR[row + 1, 1] = _transpose_front(ar)
+    normalize!(envs, ψ, operator, ψ)
 
-                ψ.AR[row + 1, end] = _transpose_front(
-                    ψ.C[row + 1, end - 1] \ _transpose_tail(al * c)
-                )
-                ψ.AC[row + 1, 1] = _transpose_front(c * ar)
-            end
+    # update the edge
+    ψ.AC[:, end] .= _mul_front.(ψ.C[:, end - 1], ψ.AR[:, end])
+    ψ.AC[:, 1] .= _mul_tail.(ψ.AL[:, 1], ψ.C[:, 1])
+    ψ.AR[:, 1] .= _transpose_front.(ψ.C[:, end] .\ _transpose_tail.(ψ.AC[:, 1]))
+    ac2 = AC2(ψ, 0; kind = :ACAR)
+    h = AC2_hamiltonian(0, ψ, operator, ψ, envs; alg.backend, allocator)
+    _, ac2′ = fixedpoint(h, ac2, :LM, alg_eigsolve)
 
-            transfer_leftenv!(envs, ψ, operator, ψ, 1)
-            transfer_rightenv!(envs, ψ, operator, ψ, 0)
+    for row in 1:size(ψ, 1)
+        al, c, ar, ϵ_truncs[row + 1, end] = svd_trunc!(ac2′[row]; trunc = alg.trunc, alg = alg.alg_svd)
+        normalize!(c)
 
-            # update error
-            ϵ = sum(zip(C_current, ψ.C[:, 0])) do (c1, c2)
-                smallest = infimum(_firstspace(c1), _firstspace(c2))
-                e1 = isometry(_firstspace(c1), smallest)
-                e2 = isometry(_firstspace(c2), smallest)
-                return norm(e2' * c2 * e2 - e1' * c1 * e1)
-            end
+        ψ.AL[row + 1, end] = al
+        ψ.C[row + 1, end] = complex(c)
+        ψ.AR[row + 1, 1] = _transpose_front(ar)
 
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
-        end
+        ψ.AR[row + 1, end] = _transpose_front(
+            ψ.C[row + 1, end - 1] \ _transpose_tail(al * c)
+        )
+        ψ.AC[row + 1, 1] = _transpose_front(c * ar)
     end
 
-    alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ)
-    ψ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter)
-
-    recalculate!(envs, ψ, operator, ψ)
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, truncation_errors = parent(ϵ_truncs), numiter = iter)
+    transfer_leftenv!(envs, ψ, operator, ψ, 1)
+    transfer_rightenv!(envs, ψ, operator, ψ, 0)
+    return state
 end
diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl
index 55c006cef..7ea87bb6b 100644
--- a/src/algorithms/statmech/vomps.jl
+++ b/src/algorithms/statmech/vomps.jl
@@ -61,7 +61,7 @@ function dominant_eigsolve(
         which
     )
     @assert which === :LM "VOMPS only supports the LM eigenvalue problem"
-    log = IterLog("VOMPS")
+    log = IterLog(alg)
     iter = 0
     ϵ = calc_galerkin(mps, operator, mps, envs; alg.backend)
     alg_environments = adapt_solver(alg.alg_environments; iter, g_global = ϵ)
@@ -90,7 +90,7 @@ function dominant_eigsolve(
     end
 end
 
-function Base.iterate(it::IterativeSolver{<:VOMPS}, state)
+function Base.iterate(it::IterativeSolver{<:VOMPS}, state::VOMPSState)
     ACs = localupdate_step!(it, state)
     mps = gauge_step!(it, state, ACs)
     envs = envs_step!(it, state, mps)
diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl
index af4707f2d..c87d74684 100644
--- a/src/algorithms/timestep/tdvp.jl
+++ b/src/algorithms/timestep/tdvp.jl
@@ -150,70 +150,74 @@ function timestep!(
     )
 end
 
-function _timestep_finite!(
-        ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs, allocator;
-        imaginary_evolution::Bool, normalize::Bool
+# Start times of the forward center update and of the backward update that follows it, for the
+# half-sweep in `direction` of a step `t → t + dt`
+_half_sweep_times(::Val{:right}, t, dt) = (t, t + dt / 2)
+_half_sweep_times(::Val{:left}, t, dt) = (t + dt / 2, t + dt)
+
+function local_update!(
+        site, direction::Val, ψ, H, alg::TDVP, envs, t, dt, allocator;
+        imaginary_evolution, normalize
     )
-    ϵ_truncs = zeros(real(scalartype(ψ)), length(ψ) - 1)
+    t_AC, t_C = _half_sweep_times(direction, t, dt)
 
-    # sweep left to right
-    for i in 1:(length(ψ) - 1)
-        # 1. optionally expand the bond ahead of the local update (CBE)
-        isnothing(alg.alg_expand) ||
-            changebond!(i, Val(:right), ψ, H, alg.alg_expand, envs; normalize, allocator)
-
-        # 2. evolve the (possibly expanded) center tensor forward
-        Hac = AC_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator)
-        AC = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator; imaginary_evolution)
-
-        # 3. gauge: split AC -> AL[i], C[i] (QR center-move, or truncated SVD cutting the
-        #    enlarged bond back down) and move the center to i+1. By default the norm is
-        #    preserved; `normalize` renormalizes.
-        _, ϵ_truncs[i] = left_gauge!(ψ, i, AC, alg.alg_gauge; normalize)
-
-        # 4. evolve the bond tensor backward
-        Hc = C_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator)
-        ψ.C[i] = integrate(
-            Hc, ψ.C[i], t + dt / 2, -dt / 2, alg.integrator;
-            imaginary_evolution
-        )
+    # at the far end of the sweep there is no bond ahead: only evolve the center tensor
+    if site == _sweep_end(ψ, direction)
+        Hac = AC_hamiltonian(site, ψ, H, ψ, envs; alg.backend, allocator)
+        ψ.AC[site] = integrate(Hac, ψ.AC[site], t_AC, dt / 2, alg.integrator; imaginary_evolution)
+        return ψ, zero(real(scalartype(ψ)))
     end
 
-    # edge case
-    Hac = AC_hamiltonian(length(ψ), ψ, H, ψ, envs; alg.backend, allocator)
-    ψ.AC[end] = integrate(Hac, ψ.AC[end], t, dt / 2, alg.integrator; imaginary_evolution)
-
-    # sweep right to left
-    for i in length(ψ):-1:2
-        # 1. optionally expand the bond ahead of the local update (CBE)
-        isnothing(alg.alg_expand) ||
-            changebond!(i, Val(:left), ψ, H, alg.alg_expand, envs; normalize, allocator)
-
-        # 2. evolve the (possibly expanded) center tensor forward
-        Hac = AC_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator)
-        AC = integrate(
-            Hac, ψ.AC[i], t + dt / 2, dt / 2, alg.integrator;
-            imaginary_evolution
-        )
+    # 1. optionally expand the bond ahead of the local update (CBE)
+    isnothing(alg.alg_expand) ||
+        changebond!(site, direction, ψ, H, alg.alg_expand, envs; normalize, allocator)
 
-        # 3. gauge: split AC -> C[i-1], AR[i] and move the center to i-1 (norm preserved by
-        #    default; `normalize` renormalizes)
-        _, ϵ_truncs[i - 1] = right_gauge!(ψ, i, AC, alg.alg_gauge; normalize)
+    # 2. evolve the (possibly expanded) center tensor forward
+    Hac = AC_hamiltonian(site, ψ, H, ψ, envs; alg.backend, allocator)
+    AC = integrate(Hac, ψ.AC[site], t_AC, dt / 2, alg.integrator; imaginary_evolution)
 
-        # 4. evolve the bond tensor backward
-        Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator)
-        ψ.C[i - 1] = integrate(
-            Hc, ψ.C[i - 1], t + dt, -dt / 2, alg.integrator;
-            imaginary_evolution
-        )
+    # 3. gauge: split AC onto the bond ahead (QR center-move, or truncated SVD cutting the
+    #    enlarged bond back down) and move the center across it. By default the norm is
+    #    preserved; `normalize` renormalizes.
+    if direction === Val(:right)
+        _, ϵ = left_gauge!(ψ, site, AC, alg.alg_gauge; normalize)
+        bond = site
+    else
+        _, ϵ = right_gauge!(ψ, site, AC, alg.alg_gauge; normalize)
+        bond = site - 1
     end
 
-    # edge case
-    Hac = AC_hamiltonian(1, ψ, H, ψ, envs; alg.backend, allocator)
-    ψ.AC[1] = integrate(
-        Hac, ψ.AC[1], t + dt / 2, dt / 2, alg.integrator;
-        imaginary_evolution
+    # 4. evolve the bond tensor backward
+    Hc = C_hamiltonian(bond, ψ, H, ψ, envs; alg.backend, allocator)
+    ψ.C[bond] = integrate(Hc, ψ.C[bond], t_C, -dt / 2, alg.integrator; imaginary_evolution)
+
+    return ψ, ϵ
+end
+
+function _timestep_finite!(
+        ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs, allocator;
+        imaginary_evolution::Bool, normalize::Bool
     )
+    L = length(ψ)
+    ϵ_truncs = zeros(real(scalartype(ψ)), L - 1)
+
+    # left→right half-sweep: `t → t + dt / 2`
+    for site in 1:L
+        ψ, ϵ = local_update!(
+            site, Val(:right), ψ, H, alg, envs, t, dt, allocator;
+            imaginary_evolution, normalize
+        )
+        site < L && (ϵ_truncs[site] = ϵ)
+    end
+
+    # right→left half-sweep: `t + dt / 2 → t + dt`
+    for site in L:-1:1
+        ψ, ϵ = local_update!(
+            site, Val(:left), ψ, H, alg, envs, t, dt, allocator;
+            imaginary_evolution, normalize
+        )
+        site > 1 && (ϵ_truncs[site - 1] = ϵ)
+    end
 
     return ψ, envs, AlgorithmInfo(; truncation_errors = ϵ_truncs)
 end
@@ -271,48 +275,56 @@ function timestep!(
     )
 end
 
+function local_update!(
+        pos, direction::Val, ψ, H, alg::TDVP2, envs, t, dt, allocator;
+        imaginary_evolution, normalize
+    )
+    t_AC2, t_AC = _half_sweep_times(direction, t, dt)
+
+    # 1. evolve the two-site center tensor at `(pos, pos + 1)` forward
+    ac2 = if direction === Val(:right)
+        _transpose_front(ψ.AC[pos]) * _transpose_tail(ψ.AR[pos + 1])
+    else
+        _transpose_front(ψ.AL[pos]) * _transpose_tail(ψ.AC[pos + 1])
+    end
+    Hac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs; alg.backend, allocator)
+    ac2′ = integrate(Hac2, ac2, t_AC2, dt / 2, alg.integrator; imaginary_evolution)
+
+    # 2. gauge: the two-site center always has to be split back up, so this is always a
+    #    truncated SVD, and the norm of the discarded singular values is the truncation error
+    alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trunc)
+    _, ϵ = gauge2!(ψ, pos, direction, ac2′, alg_gauge; normalize)
+
+    # 3. evolve the new single-site center backward, except at the far end of the sweep
+    if direction === Val(:right) ? pos != length(ψ) - 1 : pos != 1
+        site = direction === Val(:right) ? pos + 1 : pos
+        Hac = AC_hamiltonian(site, ψ, H, ψ, envs; alg.backend, allocator)
+        ψ.AC[site] = integrate(Hac, ψ.AC[site], t_AC, -dt / 2, alg.integrator; imaginary_evolution)
+    end
+
+    return ψ, ϵ
+end
+
 function _timestep2_finite!(
         ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, envs, allocator;
         imaginary_evolution::Bool, normalize::Bool
     )
-    # the two-site center always has to be split back up, so the gauge is always a truncated SVD
-    alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trunc)
-
     ϵ_truncs = zeros(real(scalartype(ψ)), length(ψ) - 1)
 
-    # sweep left to right
-    for i in 1:(length(ψ) - 1)
-        ac2 = _transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1])
-        Hac2 = AC2_hamiltonian(i, ψ, H, ψ, envs; alg.backend, allocator)
-        ac2′ = integrate(Hac2, ac2, t, dt / 2, alg.integrator; imaginary_evolution)
-
-        # the norm of the discarded singular values is the truncation error
-        _, ϵ_truncs[i] = gauge2!(ψ, i, Val(:right), ac2′, alg_gauge; normalize)
-
-        if i != (length(ψ) - 1)
-            Hac = AC_hamiltonian(i + 1, ψ, H, ψ, envs; alg.backend, allocator)
-            ψ.AC[i + 1] = integrate(
-                Hac, ψ.AC[i + 1], t + dt / 2, -dt / 2, alg.integrator;
-                imaginary_evolution
-            )
-        end
+    # left→right half-sweep: `t → t + dt / 2`
+    for pos in 1:(length(ψ) - 1)
+        ψ, ϵ_truncs[pos] = local_update!(
+            pos, Val(:right), ψ, H, alg, envs, t, dt, allocator;
+            imaginary_evolution, normalize
+        )
     end
 
-    # sweep right to left
-    for i in length(ψ):-1:2
-        ac2 = _transpose_front(ψ.AL[i - 1]) * _transpose_tail(ψ.AC[i])
-        Hac2 = AC2_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator)
-        ac2′ = integrate(Hac2, ac2, t + dt / 2, dt / 2, alg.integrator; imaginary_evolution)
-
-        _, ϵ_truncs[i - 1] = gauge2!(ψ, i - 1, Val(:left), ac2′, alg_gauge; normalize)
-
-        if i != 2
-            Hac = AC_hamiltonian(i - 1, ψ, H, ψ, envs; alg.backend, allocator)
-            ψ.AC[i - 1] = integrate(
-                Hac, ψ.AC[i - 1], t + dt, -dt / 2, alg.integrator;
-                imaginary_evolution
-            )
-        end
+    # right→left half-sweep: `t + dt / 2 → t + dt`
+    for pos in (length(ψ) - 1):-1:1
+        ψ, ϵ_truncs[pos] = local_update!(
+            pos, Val(:left), ψ, H, alg, envs, t, dt, allocator;
+            imaginary_evolution, normalize
+        )
     end
 
     return ψ, envs, AlgorithmInfo(; truncation_errors = ϵ_truncs)
diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl
index 9d9e87f40..61c64bd69 100644
--- a/src/algorithms/timestep/time_evolve.jl
+++ b/src/algorithms/timestep/time_evolve.jl
@@ -36,43 +36,75 @@ evolution at `verbosity ≥ 2`.
 """
 function time_evolve end, function time_evolve! end
 
+# Internal state of `time_evolve`: the first `iter` steps of `t_span` have been taken, each by
+# `stepper` (`timestep` or `timestep!`) called with `kwargs`
+struct TimeEvolveState{S, O, E, T, F, K}
+    mps::S
+    operator::O
+    envs::E
+    iter::Int
+    t_span::T
+    stepper::F
+    kwargs::K
+end
+
+# one iteration is a single time step followed by `finalize`, which receives the step's start time
+function Base.iterate(it::IterativeSolver{<:Any, <:TimeEvolveState}, state::TimeEvolveState)
+    state.iter >= length(state.t_span) - 1 && return nothing
+    iter = state.iter + 1
+    t = state.t_span[iter]
+    dt = state.t_span[iter + 1] - t
+
+    ψ, envs, info = state.stepper(
+        state.mps, state.operator, t, dt, it.alg, state.envs; state.kwargs...
+    )
+    ψ, envs = it.finalize(t, ψ, state.operator, envs)::Tuple{typeof(ψ), typeof(envs)}
+
+    it.state = TimeEvolveState(
+        ψ, state.operator, envs, iter, state.t_span, state.stepper, state.kwargs
+    )
+    return (ψ, envs, info), it.state
+end
+
 for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time_evolve!))
     @eval function $time_evolve(
             ψ, H, t_span::AbstractVector{<:Number}, alg,
             envs = environments(ψ, H, ψ);
             verbosity::Int = 0, imaginary_evolution::Bool = false, normalize::Bool = false
         )
-        log = IterLog(string(nameof(typeof(alg))))
-        truncation_errors = []
-        ϵ_max = 0.0
-        with_verbosity(; verbosity) do
-            @log_initialization loginit!(log, 0.0, first(t_span))
-            for iter in 1:(length(t_span) - 1)
-                t = t_span[iter]
-                dt = t_span[iter + 1] - t
-
-                ψ, envs, info_step = $timestep(
-                    ψ, H, t, dt, alg, envs; imaginary_evolution, normalize
-                )
-                ψ, envs = alg.finalize(t, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)}
-
-                # the log shows the largest per-bond error, or zero for a non-truncating algorithm
-                ϵ_step = 0.0
-                if haskey(info_step, :truncation_errors)
-                    push!(truncation_errors, info_step.truncation_errors)
-                    ϵ_step = Float64(maximum(info_step.truncation_errors; init = 0.0))
-                end
-                ϵ_max = max(ϵ_max, ϵ_step)
-                @log_iteration logiter!(log, iter, ϵ_step, t)
+        state = TimeEvolveState(
+            ψ, H, envs, 0, t_span, $timestep, (; imaginary_evolution, normalize)
+        )
+        return _time_evolve(alg, state; verbosity)
+    end
+end
+
+function _time_evolve(alg, state::TimeEvolveState; verbosity::Int)
+    log = IterLog(alg)
+    # the state type is left abstract, since the first step may promote a real state to complex
+    it = IterativeSolver{typeof(alg), TimeEvolveState}(alg, state)
+    t_span = state.t_span
+    truncation_errors = []
+    ϵ_max = 0.0
+    with_verbosity(; verbosity) do
+        @log_initialization loginit!(log, 0.0, first(t_span))
+        for (_, _, info_step) in it
+            # the log shows the largest per-bond error, or zero for a non-truncating algorithm
+            ϵ_step = 0.0
+            if haskey(info_step, :truncation_errors)
+                push!(truncation_errors, info_step.truncation_errors)
+                ϵ_step = Float64(maximum(info_step.truncation_errors; init = 0.0))
             end
-            @log_convergence logfinish!(log, length(t_span), ϵ_max, t_span[end])
+            ϵ_max = max(ϵ_max, ϵ_step)
+            @log_iteration logiter!(log, it.iter, ϵ_step, t_span[it.iter])
         end
-        info = AlgorithmInfo(;
-            numiter = length(t_span) - 1,
-            truncation_errors = isempty(truncation_errors) ? nothing : copy(truncation_errors)
-        )
-        return ψ, envs, info
+        @log_convergence logfinish!(log, length(t_span), ϵ_max, t_span[end])
     end
+    info = AlgorithmInfo(;
+        numiter = length(t_span) - 1,
+        truncation_errors = isempty(truncation_errors) ? nothing : copy(truncation_errors)
+    )
+    return it.state.mps, it.state.envs, info
 end
 
 """
diff --git a/src/states/ortho.jl b/src/states/ortho.jl
index d7a53a82d..966428cf4 100644
--- a/src/states/ortho.jl
+++ b/src/states/ortho.jl
@@ -61,6 +61,9 @@ Used as the `alg` argument of [`gaugefix!`](@ref).
     eig_miniter::Int = 10
 end
 
+IterativeLoggers.IterLog(::LeftCanonical) = IterLog("LC")
+IterativeLoggers.IterLog(::RightCanonical) = IterLog("RC")
+
 """
 $(TYPEDEF)
 
@@ -228,7 +231,7 @@ function uniform_leftorth!(
     C[end] = normalize!(C₀)
     return with_verbosity(; alg.verbosity) do
         # initialize algorithm and temporary variables
-        log = IterLog("LC")
+        log = IterLog(alg)
         A_tail = _transpose_tail.(A) # pre-transpose A
         CA_tail = similar.(A_tail)  # pre-allocate workspace
         state = (; AL, C, A, A_tail, CA_tail, iter = 0, ϵ = Inf, timeroutput, backend, allocator)
@@ -250,7 +253,7 @@ function uniform_leftorth!(
     end
 end
 
-function Base.iterate(it::IterativeSolver{LeftCanonical}, state = it.state)
+function Base.iterate(it::IterativeSolver{LeftCanonical}, state::NamedTuple)
     timeroutput = state.timeroutput
     C₀ = @timeit timeroutput "gauge_eigsolve" gauge_eigsolve_step!(it, state)
     C₁ = @timeit timeroutput "gauge_orth" gauge_orth_step!(it, state)
@@ -297,7 +300,7 @@ function uniform_rightorth!(
     C[end] = normalize!(C₀)
     return with_verbosity(; alg.verbosity) do
         # initialize algorithm and temporary variables
-        log = IterLog("RC")
+        log = IterLog(alg)
         AC_tail = _similar_tail.(A) # pre-allocate workspace
         state = (; AR, C, A, AC_tail, iter = 0, ϵ = Inf, timeroutput, backend, allocator)
         it = IterativeSolver(alg, state)
@@ -318,7 +321,7 @@ function uniform_rightorth!(
     end
 end
 
-function Base.iterate(it::IterativeSolver{RightCanonical}, state = it.state)
+function Base.iterate(it::IterativeSolver{RightCanonical}, state::NamedTuple)
     timeroutput = state.timeroutput
     C₀ = @timeit timeroutput "gauge_eigsolve" gauge_eigsolve_step!(it, state)
     C₁ = @timeit timeroutput "gauge_orth" gauge_orth_step!(it, state)
diff --git a/src/utility/iterativesolvers.jl b/src/utility/iterativesolvers.jl
index 8eb863d12..ce1aea330 100644
--- a/src/utility/iterativesolvers.jl
+++ b/src/utility/iterativesolvers.jl
@@ -19,3 +19,6 @@ function Base.getproperty(it::IterativeSolver{A, B}, name::Symbol) where {A, B}
 end
 
 Base.iterate(it::IterativeSolver) = iterate(it, it.state)
+
+# iteration logs are named after their algorithm unless the algorithm specifies otherwise
+IterativeLoggers.IterLog(alg::Algorithm) = IterLog(string(nameof(typeof(alg))))
diff --git a/src/utility/utility.jl b/src/utility/utility.jl
index d527e06c3..88dce8c20 100644
--- a/src/utility/utility.jl
+++ b/src/utility/utility.jl
@@ -40,6 +40,30 @@ end
 _firstspace(t::AbstractTensorMap) = space(t, 1)
 _lastspace(t::AbstractTensorMap) = space(t, numind(t))
 
+"""
+    bond_change(C_old, C_new)
+
+Norm of the change of a bond tensor. If the bond space changed, both tensors are compared on
+their common subspace, i.e. after restricting them with `isometry(V, infimum(V_old, V_new))`.
+For vectors of bond tensors, such as a column of a multiline MPS, the changes are summed.
+"""
+function bond_change(C_old::AbstractTensorMap, C_new::AbstractTensorMap)
+    V_old, V_new = _firstspace(C_old), _firstspace(C_new)
+    V_old == V_new && return norm(C_new - C_old)
+    # the isometry onto a subspace selects the leading rows and columns of every block
+    V = infimum(V_old, V_new)
+    ϵ² = zero(real(scalartype(C_new)))
+    for c in sectors(V)
+        d = dim(V, c)
+        b_old = view(block(C_old, c), 1:d, 1:d)
+        b_new = view(block(C_new, c), 1:d, 1:d)
+        ϵ² += dim(c) * norm(b_new - b_old)^2
+    end
+    return sqrt(ϵ²)
+end
+bond_change(C_old::AbstractVector, C_new::AbstractVector) =
+    sum(splat(bond_change), zip(C_old, C_new))
+
 """
     similar_scalartype(T::Type{<:AbstractTensorMap}, S::Type{<:Number})