From 21d8f307fa193573408ae4e2c370425bdf92d1bc Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Sun, 26 Jul 2026 13:05:33 -0400
Subject: [PATCH 1/5] fix: operator-free two-site MPO derivative contraction

The `MPO_AC2_Hamiltonian{<:MPSBondTensor, Nothing, Nothing, <:MPSBondTensor}`
action (the operator-free / overlap two-site derivative) was missing a
semicolon in the `rightenv` index, declaring a `(2, 0)` index pattern for a
`(1, 1)` bond tensor. This threw a `TensorOperations.IndexError` for the planar
bond-tensor case (e.g. `approximate(x, b, DMRG2())` on planar states); the
ComplexSpace path happened to avoid the broken method.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/algorithms/derivatives/mpo_derivatives.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl
index 46146377b..1e7bde691 100644
--- a/src/algorithms/derivatives/mpo_derivatives.jl
+++ b/src/algorithms/derivatives/mpo_derivatives.jl
@@ -138,7 +138,7 @@ function (h::MPO_AC2_Hamiltonian{<:MPSBondTensor, Nothing, Nothing, <:MPSBondTen
     )
     backend, allocator = h.backend, h.allocator
     @plansor backend = backend allocator = allocator begin
-        y[-1 -2; -3 -4] ≔ h.leftenv[-1; 1] * x[1 -2; 2 -4] * h.rightenv[2 -3]
+        y[-1 -2; -3 -4] ≔ h.leftenv[-1; 1] * x[1 -2; 2 -4] * h.rightenv[2; -3]
     end
     return y isa AbstractBlockTensorMap ? only(y) : y
 end

From 62523f41cf4a733f6a2bac2a240728fa87a126c8 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Sun, 26 Jul 2026 13:05:56 -0400
Subject: [PATCH 2/5] feat: optional truncation-error floor for DynamicTol

Add a `truncation_factor` field to `DynamicTol` and let `adapt_solver` floor the
retuned tolerance at `truncation_factor * eps_trunc`, i.e.
`tol = clamp(max(truncation_factor*eps_trunc, tol_factor*g_global/sqrt(iter)), tol_min, tol_max)`.
The default `truncation_factor = 0` keeps existing behavior (DMRG's gauge /
eigensolver uses are unchanged); a positive value is used by truncating two-site
sweeps to avoid solving a bond far below what the SVD truncation will discard.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/utility/dynamictols.jl | 32 +++++++++++++++++++++++---------
 1 file changed, 23 insertions(+), 9 deletions(-)

diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl
index 5ef4cc9f1..8c39e44cd 100644
--- a/src/utility/dynamictols.jl
+++ b/src/utility/dynamictols.jl
@@ -65,26 +65,40 @@ struct DynamicTol{A} <: Algorithm
     "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}
+    "factor on the local truncation error, setting a tolerance floor (0 disables it)"
+    truncation_factor::Float64
+
+    function DynamicTol(
+            alg::A, tol_min::Real, tol_max::Real, tol_factor::Real, truncation_factor::Real = 0.0
+        ) 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)
+        truncation_factor >= 0 ||
+            throw(ArgumentError("truncation_factor must be non-negative"))
+        return new{A}(alg, tol_min, tol_max, tol_factor, truncation_factor)
     end
 end
-function DynamicTol(alg; tol_min = 1.0e-6, tol_max = 1.0e-2, tol_factor = 0.1)
-    return DynamicTol(alg, tol_min, tol_max, tol_factor)
+function DynamicTol(alg; tol_min = 1.0e-6, tol_max = 1.0e-2, tol_factor = 0.1, truncation_factor = 0.0)
+    return DynamicTol(alg, tol_min, tol_max, tol_factor, truncation_factor)
 end
 
 """
-    adapt_solver(alg::DynamicTol; iter, g_global, ...)
+    adapt_solver(alg::DynamicTol; iter, g_global, eps_trunc, ...)
 
 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:
+global gradient / convergence-error scalar `g_global`, damped by the iteration count, but never
+below the truncation-error floor `truncation_factor · eps_trunc` (relevant for truncating two-site
+sweeps; with the default `truncation_factor = 0` the floor is inactive):
 
-    tol = clamp(tol_factor · g_global / √iter, tol_min, tol_max)
+    tol = clamp(max(truncation_factor · eps_trunc, tol_factor · g_global / √iter), tol_min, tol_max)
 """
-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)
+function adapt_solver(
+        alg::DynamicTol; iter::Integer = 1, g_global::Real = 0.0, eps_trunc::Real = 0.0, kwargs...
+    )
+    tol = clamp(
+        max(alg.truncation_factor * eps_trunc, alg.tol_factor * g_global / sqrt(max(iter, 1))),
+        alg.tol_min, alg.tol_max
+    )
     return _updatetol(alg.alg, tol)
 end
 

From 6ee774a9ea739e1170d08ef55abcfd01d3793250 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Wed, 14 Oct 2026 19:25:31 -0400
Subject: [PATCH 3/5] refactor: shared sweep geometry traits and
 setproperties-based state updates

Finite sweeping algorithms declare `sweep_sites(alg)` (1 for site updates, 2 for bond
updates), from which the shared `num_updates` and default `sweep_ranges` are derived.
`DMRG`/`DMRG2` and their `approximate` sweeps use these; `TDVP`, `TDVP2` and `BUG` keep their
own `sweep_ranges` methods. Iterator states are updated with `Accessors.setproperties`
instead of repeating every field.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/algorithms/approximate/fvomps.jl | 14 +++++-------
 src/algorithms/approximate/idmrg.jl  |  5 +----
 src/algorithms/groundstate/dmrg.jl   | 33 ++++++++--------------------
 src/algorithms/statmech/idmrg.jl     |  5 +----
 src/algorithms/timestep/bug.jl       |  2 +-
 src/algorithms/timestep/tdvp.jl      |  6 ++---
 src/utility/iterativesolvers.jl      | 12 ++++++++++
 7 files changed, 33 insertions(+), 44 deletions(-)

diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl
index ed17ce1a0..b3ced70b3 100644
--- a/src/algorithms/approximate/fvomps.jl
+++ b/src/algorithms/approximate/fvomps.jl
@@ -38,20 +38,18 @@ end
 
 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 = setproperties(state; ϵ = 0.0)
     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)
+    it.state = setproperties(state; mps = ψ, envs, iter)
     return (ψ, envs, state.ϵ), it.state
 end
 
 function sweep!(it::IterativeSolver{<:DMRG2}, state::ApproximateState, direction, iter)
-    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    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
@@ -64,11 +62,11 @@ function sweep!(it::IterativeSolver{<:DMRG2}, state::ApproximateState, direction
         ψ.AC[pos] = (al, complex(c))
         ψ.AC[pos + 1] = (complex(c), _transpose_front(ar))
     end
-    return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator)
+    return setproperties(state; mps = ψ, ϵ)
 end
 
 function sweep!(it::IterativeSolver{<:DMRG}, state::ApproximateState, direction, iter)
-    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    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
@@ -78,5 +76,5 @@ function sweep!(it::IterativeSolver{<:DMRG}, state::ApproximateState, direction,
 
         ψ.AC[pos] = AC′
     end
-    return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator)
+    return setproperties(state; mps = ψ, ϵ)
 end
diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl
index 9fe55fbf7..0b74cf1f8 100644
--- a/src/algorithms/approximate/idmrg.jl
+++ b/src/algorithms/approximate/idmrg.jl
@@ -46,10 +46,7 @@ function Base.iterate(it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state::IDMRGS
     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,
-    )
+    it.state = setproperties(state; iter, ϵ)
     return (it.state.mps, it.state.envs, it.state.ϵ), it.state
 end
 
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index 1b9575e3a..64f67aac0 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -243,20 +243,11 @@ function local_update!(
     return ψ, ϵ_local, ϵ_trunc, decay_rate
 end
 
-# 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_sites(::DMRG) = 1
+sweep_sites(::DMRG2) = 2
 
-_sweep_ranges(::DMRG, ψ) = (1:(length(ψ) - 1), length(ψ):-1:2)
-_sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 2):-1:1)
-
-# per-bond truncation errors of the last cut from the per-update slots: single-site slot `pos`
-# holds the backward cut of bond `pos - 1`, which is the last one made there
-_bond_truncation_errors(::DMRG, ϵ_truncs) = ϵ_truncs[2:end]
-_bond_truncation_errors(::DMRG2, ϵ_truncs) = copy(ϵ_truncs)
+# single-site slot `pos` holds the backward cut of bond `pos - 1`, which is the last one made there
+bond_truncation_errors(alg, ϵ_truncs) = sweep_sites(alg) == 1 ? ϵ_truncs[2:end] : copy(ϵ_truncs)
 
 inner_alg_gauge(alg::Union{DMRG, DMRG2}) = alg_gauge(alg.alg_gauge)
 
@@ -310,7 +301,7 @@ end
 
 function DMRGState(ψ, H, alg::Union{DMRG, DMRG2}, envs, allocator, timeroutput)
     Tr = real(scalartype(ψ))
-    n = _num_updates(alg, ψ)
+    n = num_updates(alg, ψ)
     local_errors = ones(Tr, n)
     return DMRGState(
         ψ, H, envs, 0, maximum(local_errors), local_errors, zeros(Tr, n),
@@ -319,7 +310,7 @@ function DMRGState(ψ, H, alg::Union{DMRG, DMRG2}, envs, allocator, timeroutput)
 end
 
 function sweep!(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::DMRGState, direction, iter)
-    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    fwd, bwd = sweep_ranges(it.alg, state.mps)
     sites = direction === Val(:right) ? fwd : bwd
     ψ, ϵ = state.mps, state.ϵ
     for pos in sites
@@ -332,10 +323,7 @@ function sweep!(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::DMRGState, dir
         # 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,
-    )
+    return setproperties(state; mps = ψ, ϵ)
 end
 
 function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::DMRGState)
@@ -348,10 +336,7 @@ function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::DMRGStat
     ψ, 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,
-    )
+    it.state = setproperties(state; mps = ψ, envs, iter)
     return (ψ, envs, state.ϵ), it.state
 end
 
@@ -384,7 +369,7 @@ function find_groundstate_sweep!(
     state = it.state
     info = AlgorithmInfo(;
         converged = sweep_converged(alg, state), galerkin = state.ϵ,
-        truncation_errors = _bond_truncation_errors(alg, state.truncation_errors),
+        truncation_errors = bond_truncation_errors(alg, state.truncation_errors),
         numiter = state.iter,
     )
     return state.mps, state.envs, info
diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl
index fbfcea38c..8e4958c1c 100644
--- a/src/algorithms/statmech/idmrg.jl
+++ b/src/algorithms/statmech/idmrg.jl
@@ -47,10 +47,7 @@ function Base.iterate(it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state::IDMRGS
     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,
-    )
+    it.state = setproperties(state; iter, ϵ)
     return (it.state.mps, it.state.envs, it.state.ϵ), it.state
 end
 
diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl
index f9e79e199..b9c5f9dba 100644
--- a/src/algorithms/timestep/bug.jl
+++ b/src/algorithms/timestep/bug.jl
@@ -168,7 +168,7 @@ function local_update!(site, direction::Val, ψ, alg::BUG, state::TimeStepState)
 end
 
 # like single-site TDVP, BUG visits every site and only evolves the center at the far end
-_sweep_ranges(::BUG, ψ) = (1:length(ψ), length(ψ):-1:1)
+sweep_ranges(::BUG, ψ) = (1:length(ψ), length(ψ):-1:1)
 
 _promotes_to_complex(ψ, H, dt, ::BUG; imaginary_evolution) =
     scalartype(ψ) <: Real && !(imaginary_evolution && isreal(dt) && scalartype(H) <: Real)
diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl
index 1c0c126a3..af9331c4c 100644
--- a/src/algorithms/timestep/tdvp.jl
+++ b/src/algorithms/timestep/tdvp.jl
@@ -205,14 +205,14 @@ end
 
 # Single-site integrators visit every site, finishing a half-sweep by evolving only the center
 # tensor at its far end.
-_sweep_ranges(::TDVP, ψ) = (1:length(ψ), length(ψ):-1:1)
+sweep_ranges(::TDVP, ψ) = (1:length(ψ), length(ψ):-1:1)
 
 # bond cut by the update at `site`, or `nothing` at the far end of a single-site half-sweep
 _swept_bond(alg, site, ψ, ::Val{:right}) = site < length(ψ) ? site : nothing
 _swept_bond(alg, site, ψ, ::Val{:left}) = site > 1 ? site - 1 : nothing
 
 function sweep!(alg, state::TimeStepState, direction)
-    fwd, bwd = _sweep_ranges(alg, state.mps)
+    fwd, bwd = sweep_ranges(alg, state.mps)
     sites = direction === Val(:right) ? fwd : bwd
     ψ = state.mps
     @timeit state.timeroutput "half-sweep" for site in sites
@@ -310,7 +310,7 @@ Used as the `algorithm` argument of [`timestep`](@ref), [`timestep!`](@ref) and
 end
 
 # two-site TDVP visits every bond, and each update cuts the bond it acts on
-_sweep_ranges(::TDVP2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 1):-1:1)
+sweep_ranges(::TDVP2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 1):-1:1)
 _swept_bond(::TDVP2, pos, ψ, ::Val{:right}) = pos
 _swept_bond(::TDVP2, pos, ψ, ::Val{:left}) = pos
 
diff --git a/src/utility/iterativesolvers.jl b/src/utility/iterativesolvers.jl
index ce1aea330..c9a53ff56 100644
--- a/src/utility/iterativesolvers.jl
+++ b/src/utility/iterativesolvers.jl
@@ -22,3 +22,15 @@ 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))))
+
+# Finite sweeps update `sweep_sites(alg)` neighbouring sites at once: 1 for a site, 2 for a bond.
+function sweep_sites end
+
+# number of distinct local updates in a sweep: all sites, or all bonds
+num_updates(alg, ψ) = length(ψ) - sweep_sites(alg) + 1
+
+# forward and backward positions of a sweep, without updating the turning point twice
+function sweep_ranges(alg, ψ)
+    L = length(ψ)
+    return sweep_sites(alg) == 1 ? (1:(L - 1), L:-1:2) : (1:(L - 1), (L - 2):-1:1)
+end

From 4df3c4e3846736baf03e36b67068335ba4691b8a Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Mon, 12 Oct 2026 16:43:05 -0400
Subject: [PATCH 4/5] feat: MPS linear solver (linear_solve) category
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Add a first-class `linear_solve`/`linear_solve!` for finite MPS, solving
`(a0 + a1*A)*x = b` with `A` an MPO/MPOHamiltonian and `x`, `b` MPS. The
signature is state-first (`linear_solve(x0, A, b, [alg], [envs]; a0=0, a1=1, ...)`)
and returns `(x, envs, info)`, matching find_groundstate/approximate; the shift
follows KrylovKit's convention. It is deliberately not named `linsolve`: that
would clash with the `linsolve` KrylovKit exports for anyone using both
packages, and its argument order differs from KrylovKit's.

Algorithms are DMRG-style sweeps `DMRGSolve` (single-site) and `DMRGSolve2`
(two-site, rank-adaptive), each parameterized by a local `LinsolveFormulation`:
`Galerkin` (solve the effective system directly) and `LeastSquares` (normal
equations via squared environments). Convergence uses the local Galerkin
residual `||(a0 + a1*A_eff)*AC - b_eff|| / ||b||`, formed as a local tensor
(machine-precision floor, no catastrophic cancellation). The two-site stop is
truncation-aware. Both algorithms carry a `backend` and serve their local
updates from a per-sweep allocator, like DMRG.

The sweeps follow the iterator structure of `DMRG`/`DMRG2`: an
`IterativeSolver` over a `LinearSolveState`, whose iterations are a forward and
a backward `sweep!` of per-site `local_update!`s over `sweep_ranges`, with
`sweep_converged` as the stopping test. `info.galerkin` reports the relative local
residual, the linear-system counterpart of DMRG's Galerkin error, and
`DMRGSolve2` also its `truncation_errors`.

`LeastSquares` forms the `A†b` term of the normal equations from the mixed
sandwich <x|A|b> instead of materializing `A*b` as an MPS, avoiding a state of
bond dimension chi_A*chi_b and making the formulation usable on `WindowMPS`.
The JordanMPO effective operator assumes `below === above`, so the mixed case
goes through the generic sparse MPO derivative.

Local solves default to adaptive tolerances via `Defaults.alg_linsolve()` (a
structure-selected CG/GMRES wrapped in `DynamicTol`, with a truncation floor for
two-site); pass a plain solver to disable adaptation. KrylovKit has no
`linsolve` for `MINRES`, so the hermitian-indefinite case also uses GMRES. A
local solve only warns about non-convergence when it misses `tol * norm(b)`.

Tests cover the resolvent `(z - H)x = V|gs>` against the analytic
`1/(z - E0)` and `propagator` (Galerkin and LeastSquares, single- and two-site),
a positive-definite solve via the keyword interface, adaptive tolerances, a
general RHS checked against a dense solve, and a `WindowMPS` resolvent.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/MPSKit.jl                                |   5 +
 src/algorithms/linear_solve/flinear_solve.jl | 166 ++++++++++++++++
 src/algorithms/linear_solve/linear_solve.jl  | 187 +++++++++++++++++++
 src/utility/algorithminfo.jl                 |  13 +-
 src/utility/defaults.jl                      |  13 +-
 test/algorithms/linear_solve.jl              | 124 ++++++++++++
 6 files changed, 502 insertions(+), 6 deletions(-)
 create mode 100644 src/algorithms/linear_solve/flinear_solve.jl
 create mode 100644 src/algorithms/linear_solve/linear_solve.jl
 create mode 100644 test/algorithms/linear_solve.jl

diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index ea209f449..468e15a46 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -43,6 +43,8 @@ export NoiseSchedule, FunctionalSchedule, ExponentialDecay, Warmup, DMRG3S
 export Zipup
 export propagator
 export DynamicalDMRG, NaiveInvert, Jeckelmann
+export linear_solve, linear_solve!
+export DMRGSolve, DMRGSolve2, Galerkin, LeastSquares
 export exact_diagonalization, fidelity_susceptibility
 
 # toolbox:
@@ -209,6 +211,9 @@ include("algorithms/ED.jl")
 
 include("algorithms/unionalg.jl")
 
+include("algorithms/linear_solve/linear_solve.jl")
+include("algorithms/linear_solve/flinear_solve.jl")
+
 include("utility/show.jl")
 
 function __init__()
diff --git a/src/algorithms/linear_solve/flinear_solve.jl b/src/algorithms/linear_solve/flinear_solve.jl
new file mode 100644
index 000000000..c69fa4bd5
--- /dev/null
+++ b/src/algorithms/linear_solve/flinear_solve.jl
@@ -0,0 +1,166 @@
+# Finite-MPS sweeps for `linear_solve`. Convergence is measured, as in `find_groundstate!`, by the
+# local residual of the original system `(a₀ + a₁·A)·x = b`, evaluated before each local solve.
+
+sweep_sites(::DMRGSolve) = 1
+sweep_sites(::DMRGSolve2) = 2
+
+# environments beyond `⟨x|A|x⟩`: `LeastSquares` also needs `A²` and the mixed sandwich `⟨x|A|b⟩`
+linear_solve_environments(::Galerkin, x, A, b, openvs) = (; openvs, rhsenvs = environments(x, b))
+function linear_solve_environments(::LeastSquares, x, A, b, openvs)
+    Asq, sqenvs = squaredenvs(x, A, openvs)
+    return (; openvs, rhsenvs = environments(x, b), Asq, sqenvs, abenvs = environments(x, A, b))
+end
+
+# `⟨∂x|A|b⟩` through the generic MPO derivative, since the JordanMPO one assumes `below === above`
+function mixed_projection(alg, pos, x, A, b, envs; backend, allocator)
+    W = A isa WindowMPOHamiltonian ? A.finite_ham : A
+    GL = leftenv(envs, pos, x; backend, allocator)
+    if sweep_sites(alg) == 1
+        GR = rightenv(envs, pos, x; backend, allocator)
+        return MPO_AC_Hamiltonian(GL, W[pos], GR, backend, allocator) * b.AC[pos]
+    else
+        GR = rightenv(envs, pos + 1, x; backend, allocator)
+        return MPO_AC2_Hamiltonian(GL, W[pos], W[pos + 1], GR, backend, allocator) * AC2(b, pos)
+    end
+end
+
+function local_hamiltonian(alg, pos, x, A, envs; kwargs...)
+    return sweep_sites(alg) == 1 ? AC_hamiltonian(pos, x, A, x, envs; kwargs...) :
+        AC2_hamiltonian(pos, x, A, x, envs; kwargs...)
+end
+
+# the local system `(α₀ + α₁·M)·v = rhs` actually handed to the Krylov solver
+local_linear_system(::Galerkin, alg, pos, x, A, b, ctx, A_eff, b_eff, a₀, a₁; kwargs...) =
+    (A_eff, b_eff, a₀, a₁)
+function local_linear_system(::LeastSquares, alg, pos, x, A, b, ctx, A_eff, b_eff, a₀, a₁; kwargs...)
+    Asq_eff = local_hamiltonian(alg, pos, x, ctx.Asq, ctx.sqenvs; kwargs...)
+    N_eff = LinearCombination((A_eff, Asq_eff), (2 * real(conj(a₀) * a₁), abs2(a₁)))
+    Ab_eff = mixed_projection(alg, pos, x, A, b, ctx.abenvs; kwargs...)
+    return N_eff, conj(a₀) * b_eff + conj(a₁) * Ab_eff, abs2(a₀), one(a₁)
+end
+
+# returns the updated local tensor and the original-system residual before the update
+function local_linear_solve(alg, pos, x, state; kind = :ACAR, adapt_kwargs...)
+    A, b, ctx, (a₀, a₁) = state.operator, state.rhs, state.context, state.shift
+    backend, allocator = alg.backend, state.allocator
+    onesite = sweep_sites(alg) == 1
+    A_eff = local_hamiltonian(alg, pos, x, A, ctx.openvs; backend, allocator)
+    b_eff = onesite ? AC_projection(pos, x, b, ctx.rhsenvs; backend, allocator) :
+        AC2_projection(pos, x, b, ctx.rhsenvs; backend, allocator)
+    v = onesite ? x.AC[pos] : AC2(x, pos; kind)
+    # `‖(a₀ + a₁·A)·x − b‖`, reusing the buffer of `A·x`
+    res = norm(add!!(add!!(A_eff * v, v, a₀, a₁), b_eff, -1))
+    M, rhs, α₀, α₁ = local_linear_system(
+        alg.formulation, alg, pos, x, A, b, ctx, A_eff, b_eff, a₀, a₁; backend, allocator
+    )
+    v′, info = KrylovKit.linsolve(M, rhs, v, adapt_solver(alg.solver; adapt_kwargs...), α₀, α₁)
+    # only warn when the local solve also misses the accuracy the sweep needs
+    info.converged == 0 && info.normres > alg.tol * state.normb &&
+        @warn "linear_solve: local solve at $pos did not converge (normres = $(info.normres))"
+    return v′, res
+end
+
+function local_update!(pos, direction, x, alg::DMRGSolve, state, g_global, iter)
+    AC′, res = local_linear_solve(alg, pos, x, state; iter, g_global)
+    x.AC[pos] = AC′
+    return x, res, 0.0
+end
+function local_update!(pos, direction, x, alg::DMRGSolve2, state, g_global, iter)
+    kind = direction === Val(:right) ? :ACAR : :ALAC
+    eps_trunc = state.truncation_errors[pos]
+    AC2′, res = local_linear_solve(alg, pos, x, state; kind, iter, g_global, eps_trunc)
+    x, ϵ_trunc = gauge2!(x, pos, direction, AC2′, alg.alg_gauge; normalize = false)
+    return x, res, ϵ_trunc
+end
+
+# Internal state of the finite linear-solve sweeps, where `ϵ` is the largest local residual
+# relative to `‖b‖`
+struct LinearSolveState{S, O, B, E, C, T, A}
+    mps::S
+    operator::O
+    rhs::B
+    envs::E
+    context::C
+    shift::T
+    normb::Float64
+    iter::Int
+    ϵ::Float64
+    residuals::Vector{Float64}
+    truncation_errors::Vector{Float64}
+    allocator::A
+end
+
+function LinearSolveState(x, A, b, alg::Union{DMRGSolve, DMRGSolve2}, envs, a₀, a₁)
+    n = num_updates(alg, x)
+    context = linear_solve_environments(alg.formulation, x, A, b, envs)
+    allocator = default_allocator(x, SerialScheduler())
+    # guard the relative residual against a zero right-hand side
+    normb = norm(b)
+    iszero(normb) && (normb = one(normb))
+    return LinearSolveState(
+        x, A, b, envs, context, (a₀, a₁), normb, 0, 1.0, ones(n), zeros(n), allocator
+    )
+end
+
+function sweep!(
+        it::IterativeSolver{<:Union{DMRGSolve, DMRGSolve2}}, state::LinearSolveState,
+        direction, iter
+    )
+    fwd, bwd = sweep_ranges(it.alg, state.mps)
+    x = state.mps
+    for pos in (direction === Val(:right) ? fwd : bwd)
+        # the next local solver uses the residuals accumulated so far, as in `DMRG`
+        g_global = maximum(state.residuals) * state.normb
+        x, res, state.truncation_errors[pos] = local_update!(
+            pos, direction, x, it.alg, state, g_global, iter
+        )
+        state.residuals[pos] = res / state.normb
+    end
+    return setproperties(state; mps = x, ϵ = maximum(state.residuals))
+end
+
+function Base.iterate(
+        it::IterativeSolver{<:Union{DMRGSolve, DMRGSolve2}}, state::LinearSolveState
+    )
+    iter = state.iter + 1
+    state = sweep!(it, state, Val(:right), iter)
+    state = sweep!(it, state, Val(:left), iter)
+    x, envs = it.finalize(
+        iter, state.mps, state.operator, state.envs
+    )::Tuple{typeof(state.mps), typeof(state.envs)}
+    it.state = setproperties(state; mps = x, envs, iter)
+    return (x, envs, state.ϵ), it.state
+end
+
+# the residual cannot drop below the level set by the discarded weight
+sweep_converged(alg, state::LinearSolveState) =
+    state.ϵ <= max(alg.tol, maximum(state.truncation_errors) / state.normb)
+
+function linear_solve!(
+        x::AbstractFiniteMPS, A, b, alg::Union{DMRGSolve, DMRGSolve2},
+        envs = environments(x, A, x); a₀ = 0, a₁ = 1
+    )
+    log = IterLog(alg)
+    it = IterativeSolver(alg, LinearSolveState(x, A, b, alg, envs, a₀, a₁))
+
+    with_verbosity(; alg.verbosity) do
+        @log_initialization loginit!(log, it.ϵ)
+        for (_, _, ϵ) in Iterators.take(it, alg.maxiter)
+            if sweep_converged(alg, it.state)
+                @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
+
+    state = it.state
+    info = AlgorithmInfo(;
+        converged = sweep_converged(alg, state), galerkin = state.ϵ, numiter = state.iter,
+        truncation_errors = alg isa DMRGSolve2 ? copy(state.truncation_errors) : nothing,
+    )
+    return state.mps, state.envs, info
+end
diff --git a/src/algorithms/linear_solve/linear_solve.jl b/src/algorithms/linear_solve/linear_solve.jl
new file mode 100644
index 000000000..2b62f9ef7
--- /dev/null
+++ b/src/algorithms/linear_solve/linear_solve.jl
@@ -0,0 +1,187 @@
+# Flavours
+# --------
+"""
+$(TYPEDEF)
+
+Abstract supertype selecting how the local linear subproblem is posed in an MPS
+[`linear_solve`](@ref) sweep. See [`Galerkin`](@ref) and [`LeastSquares`](@ref).
+"""
+abstract type LinsolveFormulation end
+
+"""
+$(TYPEDEF)
+
+Local formulation imposing the Galerkin condition: the residual `(a₀ + a₁·A)·x − b` is projected
+orthogonal to the local tangent space, yielding the effective system `(a₀ + a₁·A_eff)·x = b_eff`
+which is solved with the configured Krylov `solver`. Use a `CG` solver for hermitian
+positive-definite `A`, and `GMRES`/`BiCGStab` for general (non-hermitian or indefinite) `A`.
+"""
+struct Galerkin <: LinsolveFormulation end
+
+"""
+$(TYPEDEF)
+
+Local formulation that minimizes `‖(a₀ + a₁·A)·x − b‖²`, i.e. solves the normal equations
+`M†M·x = M†b` with `M = a₀ + a₁·A`, built from squared operator environments. The normal operator
+is unconditionally positive-definite (solve it with `CG`) at the cost of squaring the condition
+number.
+
+Currently assumes `A` is hermitian (the resolvent / dynamical-DMRG case), for which
+`M†M = |a₀|² + 2·Re(ā₀·a₁)·A + |a₁|²·A²` and `M†b = ā₀·b + ā₁·A·b`. This is a generalization of the
+[`Jeckelmann`](@ref) dynamical-DMRG functional.
+"""
+struct LeastSquares <: LinsolveFormulation end
+
+# Algorithms
+# ----------
+"""
+$(TYPEDEF)
+
+Single-site DMRG-style sweeping algorithm for the MPS linear solver [`linear_solve`](@ref). It keeps
+the bond dimension of the initial guess fixed and, at each site, solves the local linear problem
+selected by `formulation` with the local Krylov `solver`.
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
+@kwdef struct DMRGSolve{L <: LinsolveFormulation, S, F, B} <: Algorithm
+    "formulation of the local subproblem, either [`Galerkin`](@ref) or [`LeastSquares`](@ref)"
+    formulation::L = Galerkin()
+    "local linear solver; a plain KrylovKit solver, or one wrapped in `DynamicTol` for per-bond adaptive tolerances (the default)"
+    solver::S = Defaults.alg_linsolve()
+    "tolerance for convergence criterium"
+    tol::Float64 = Defaults.tol
+    "maximal amount of iterations"
+    maxiter::Int = Defaults.maxiter
+    "setting for how much information is displayed"
+    verbosity::Int = Defaults.verbosity
+    "callback function applied after each iteration, of signature `finalize(iter, x, A, envs) -> x, envs`"
+    finalize::F = Defaults._finalize
+    "backend for tensor contractions and index manipulations"
+    backend::B = Defaults.backend()
+end
+
+"""
+$(TYPEDEF)
+
+Two-site DMRG-style sweeping algorithm for the MPS linear solver [`linear_solve`](@ref). Each bond
+update solves the local two-site linear problem selected by `formulation` and truncates the enlarged
+bond back down with `alg_gauge` (a truncated SVD built from `trunc`), making the bond dimension
+adaptive.
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
+struct DMRGSolve2{L <: LinsolveFormulation, S, G, F, B} <: Algorithm
+    "formulation of the local subproblem, either [`Galerkin`](@ref) or [`LeastSquares`](@ref)"
+    formulation::L
+    "local linear solver; a plain KrylovKit solver, or one wrapped in `DynamicTol` for per-bond adaptive tolerances"
+    solver::S
+    "tolerance for convergence criterium"
+    tol::Float64
+    "maximal amount of iterations"
+    maxiter::Int
+    "setting for how much information is displayed"
+    verbosity::Int
+    "truncated SVD used for the post-update gauge"
+    alg_gauge::G
+    "callback function applied after each iteration, of signature `finalize(iter, x, A, envs) -> x, envs`"
+    finalize::F
+    "backend for tensor contractions and index manipulations"
+    backend::B
+end
+function DMRGSolve2(;
+        formulation = Galerkin(), solver = Defaults.alg_linsolve(), tol = Defaults.tol,
+        maxiter = Defaults.maxiter, verbosity = Defaults.verbosity,
+        alg_svd = Defaults.alg_svd(), trunc, finalize = Defaults._finalize,
+        backend = Defaults.backend()
+    )
+    alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc)
+    return DMRGSolve2(formulation, solver, tol, maxiter, verbosity, alg_gauge, finalize, backend)
+end
+
+# Interface
+# ---------
+"""
+    linear_solve(x₀, A, b, [algorithm], [environments]; kwargs...) -> (x, environments, info)
+    linear_solve!(x₀, A, b, [algorithm], [environments]; kwargs...) -> (x, environments, info)
+
+Solve the linear system `(a₀ + a₁·A)·x = b` for the MPS `x`, where `A` is an MPO / MPOHamiltonian
+and `b` is an MPS. The `a₀ + a₁·A` shift follows `KrylovKit.linsolve` and is applied implicitly;
+the defaults `a₀ = 0`, `a₁ = 1` give the plain system `A·x = b`. `linear_solve!` solves in-place,
+overwriting `x₀`.
+
+Currently only finite MPS are supported.
+
+# Arguments
+- `x₀::AbstractMPS`: initial guess (and, for `linear_solve!`, the state overwritten with the
+  solution)
+- `A`: operator of the linear system (MPO / MPOHamiltonian)
+- `b::AbstractMPS`: right-hand side
+- `algorithm`: linear-solve algorithm, see [`DMRGSolve`](@ref) and [`DMRGSolve2`](@ref)
+- `[environments]`: MPS environment manager for the operator sandwich `⟨x|A|x⟩`
+
+# Keyword Arguments
+The keyword-based call (no explicit `algorithm`) selects an algorithm from the structure flags:
+- `tol::Float64`: tolerance for convergence criterium
+- `maxiter::Int`: maximum amount of iterations
+- `verbosity::Int`: display progress information
+- `a₀`, `a₁`: shift/scale scalars of the system `(a₀ + a₁·A)·x = b`
+- `isposdef::Bool`: declare `(a₀ + a₁·A)` positive-definite, so the local solver is `CG` instead of `GMRES`
+- `trunc`: if supplied, a truncated two-site sweep ([`DMRGSolve2`](@ref)) is prepended to
+  adapt the bond dimension before the single-site algorithm polishes the result
+
+# Returns
+- `x::AbstractMPS`: the (bond-dimension-limited) solution
+- `environments`: operator environments corresponding to `x`
+- `info::AlgorithmInfo`: how the algorithm terminated. `info.galerkin` is the Galerkin error,
+  the largest local residual `‖(a₀ + a₁·A)·x − b‖` of the last sweep relative to `‖b‖`, and `info.converged` whether
+  it met the stopping test. [`DMRGSolve2`](@ref) also reports `info.truncation_errors`. See
+  [`AlgorithmInfo`](@ref).
+"""
+function linear_solve end
+@doc (@doc linear_solve) function linear_solve! end
+
+# scalar-promote an initial state so it can hold a complex solution when the shift is complex
+function promote_to_shift(x::AbstractMPS, a₀, a₁)
+    T = promote_type(scalartype(x), typeof(a₀), typeof(a₁))
+    return T <: Complex ? complex(x) : x
+end
+
+# a loose two-site pass to grow the bond dimension (if `trunc` is given), then a single-site polish
+function linear_solve_algorithm(; solver, tol, maxiter, verbosity, trunc, backend = Defaults.backend())
+    alg = DMRGSolve(; solver, tol, maxiter, verbosity, backend)
+    isnothing(trunc) && return alg
+    return DMRGSolve2(;
+        solver, tol = min(1.0e-2, 100tol), maxiter, verbosity, trunc, backend
+    ) & alg
+end
+
+function linear_solve(
+        x₀::AbstractMPS, A, b::AbstractMPS;
+        tol = Defaults.tol, maxiter = Defaults.maxiter, verbosity = Defaults.verbosity,
+        a₀ = 0, a₁ = 1, isposdef = false, trunc = nothing
+    )
+    x₀ isa AbstractFiniteMPS ||
+        throw(ArgumentError("`linear_solve` currently only supports finite MPS"))
+    solver = Defaults.alg_linsolve(; isposdef, tol, maxiter)
+    alg = linear_solve_algorithm(; solver, tol, maxiter, verbosity, trunc)
+    return linear_solve(x₀, A, b, alg; a₀, a₁)
+end
+
+# `envs...` lets `linear_solve!` build the environments from the promoted copy when none are supplied
+function linear_solve(x₀, A, b, alg::Union{DMRGSolve, DMRGSolve2}, envs...; a₀ = 0, a₁ = 1)
+    return linear_solve!(promote_to_shift(copy(x₀), a₀, a₁), A, b, alg, envs...; a₀, a₁)
+end
+
+# sequential chaining of algorithms (e.g. a two-site pass then a single-site polish)
+function linear_solve(x₀, A, b, alg::UnionAlg, envs...; a₀ = 0, a₁ = 1)
+    x, newenvs, = linear_solve(x₀, A, b, alg.alg1, envs...; a₀, a₁)
+    return linear_solve!(x, A, b, alg.alg2, newenvs; a₀, a₁)
+end
+function linear_solve!(x₀, A, b, alg::UnionAlg, envs...; a₀ = 0, a₁ = 1)
+    x, newenvs, = linear_solve!(x₀, A, b, alg.alg1, envs...; a₀, a₁)
+    return linear_solve!(x, A, b, alg.alg2, newenvs; a₀, a₁)
+end
diff --git a/src/utility/algorithminfo.jl b/src/utility/algorithminfo.jl
index 1da88d34e..13d7b93e7 100644
--- a/src/utility/algorithminfo.jl
+++ b/src/utility/algorithminfo.jl
@@ -2,8 +2,8 @@
 $(TYPEDEF)
 
 Information about how an algorithm arrived at its result, returned as the last value by
-[`find_groundstate`](@ref), [`leading_boundary`](@ref), [`approximate`](@ref), [`timestep`](@ref)
-and [`time_evolve`](@ref).
+[`find_groundstate`](@ref), [`leading_boundary`](@ref), [`approximate`](@ref), [`linear_solve`](@ref),
+[`timestep`](@ref) and [`time_evolve`](@ref).
 
 Algorithms in MPSKit produce genuinely different measures, and not every algorithm even has access
 to the same information. To avoid reporting two different quantities under one name, the information is
@@ -28,9 +28,12 @@ it reports. Nothing prevents an algorithm from adding its own.
 The quantity that was compared against the algorithm's `tol` is stored under a name that says
 which measure it is:
 
-  - `galerkin`: the Galerkin error, i.e. the maximum over sites of the local update projected onto
-    the orthogonal complement of the current tensor. Reported by [`DMRG`](@ref), [`DMRG2`](@ref),
-    [`VUMPS`](@ref) and [`VOMPS`](@ref) when solving for a state.
+  - `galerkin`: the Galerkin error, i.e. the maximum over sites of the residual projected onto the
+    local variational space. For an eigenproblem this is the local update projected onto the
+    orthogonal complement of the current tensor, as reported by [`DMRG`](@ref), [`DMRG2`](@ref),
+    [`VUMPS`](@ref) and [`VOMPS`](@ref) when solving for a state. For a linear system it is the
+    local residual `‖(a₀ + a₁·A)·x − b‖` relative to `‖b‖`, as reported by [`DMRGSolve`](@ref) and
+    [`DMRGSolve2`](@ref) in [`linear_solve`](@ref).
   - `gradientnorm`: the norm of the Riemannian (Grassmann) gradient, as supplied by the optimiser.
     Reported by [`GradientGrassmann`](@ref).
   - `bondresidual`: the change in the center bond tensor over a sweep. This is a fixed-point
diff --git a/src/utility/defaults.jl b/src/utility/defaults.jl
index c6c100c29..4b94e0fca 100644
--- a/src/utility/defaults.jl
+++ b/src/utility/defaults.jl
@@ -5,7 +5,7 @@ Some default values and settings for MPSKit.
 """
 module Defaults
 
-import KrylovKit: GMRES, Arnoldi, Lanczos
+import KrylovKit: GMRES, CG, Arnoldi, Lanczos
 using OhMyThreads
 using ..MPSKit: DynamicTol, AdaptiveKrylov
 using Preferences: @load_preference, @set_preferences!
@@ -75,6 +75,17 @@ function alg_eigsolve(;
     return dynamic_tols ? DynamicTol(alg, tol_min, tol_max, tol_factor) : alg
 end
 
+function alg_linsolve(;
+        isposdef = false, tol = tol, maxiter = maxiter,
+        krylovdim = krylovdim, verbosity = 0, dynamic_tols = dynamic_tols,
+        tol_min = tol_min, tol_max = tol_max, tol_factor = eigs_tolfactor, truncation_factor = 1.0e-1
+    )
+    # KrylovKit has no `linsolve` for `MINRES`, so the hermitian-indefinite case also uses GMRES
+    base = isposdef ? CG(; tol, maxiter, verbosity) :
+        GMRES(; tol, maxiter, krylovdim, verbosity)
+    return dynamic_tols ? DynamicTol(base, tol_min, tol_max, tol_factor, truncation_factor) : base
+end
+
 function alg_environments(;
         tol = tol, maxiter = maxiter, verbosity = 0, krylovdim = krylovdim, eager = true,
         dynamic_tols = dynamic_tols, tol_min = tol_min, tol_max = tol_max,
diff --git a/test/algorithms/linear_solve.jl b/test/algorithms/linear_solve.jl
new file mode 100644
index 000000000..663122057
--- /dev/null
+++ b/test/algorithms/linear_solve.jl
@@ -0,0 +1,124 @@
+println("
+-----------------------------
+|   Linear solver tests      |
+-----------------------------
+")
+
+using .TestSetup
+using Test, TestExtras
+using MPSKit
+using TensorKit
+using TensorKit: ℙ
+# KrylovKit is not a direct test dependency; reach the solver types through MPSKit's re-export
+const GMRES = MPSKit.KrylovKit.GMRES
+const CG = MPSKit.KrylovKit.CG
+const BiCGStab = MPSKit.KrylovKit.BiCGStab
+const DynamicTol = MPSKit.DynamicTol
+
+verbosity_conv = 1
+
+# helper: relative residual of (a₀ + a₁·A)·x = b, computed densely (robust for small systems;
+# MPS-level `+`/`-` does not reliably cancel across states of differing bond structure)
+function rel_residual(A, x, b, a₀, a₁)
+    Am = convert(TensorMap, A)
+    xv = convert(TensorMap, x)
+    bv = convert(TensorMap, b)
+    return norm(a₀ * xv + a₁ * (Am * xv) - bv) / norm(bv)
+end
+
+@testset "linear_solve FiniteMPS" verbose = true begin
+    L = 8
+    H = force_planar(-transverse_field_ising(; L, g = -4))
+    gs, = find_groundstate(FiniteMPS(L, ℙ^2, ℙ^16), H; verbosity = verbosity_conv)
+    E₀ = real(expectation_value(gs, H))
+
+    # resolvent of an eigenstate: (z − H)⁻¹|gs⟩ = 1/(z − E₀)|gs⟩, so ⟨gs|x⟩ = 1/(z − E₀).
+    z = E₀ + 0.5 + 0.3im
+    predicted = 1 / (z - E₀)
+
+    @testset "resolvent, single-site, formulation $flav" for flav in (Galerkin(), LeastSquares())
+        solver = flav isa LeastSquares ? CG(; tol = 1.0e-12) : GMRES(; tol = 1.0e-12)
+        alg = DMRGSolve(; formulation = flav, solver, tol = 1.0e-10, verbosity = 0)
+        # solve (z − H) x = gs  ⟺  a₀ = z, a₁ = −1
+        x, envs, ϵ = linear_solve(complex(copy(gs)), H, gs, alg; a₀ = z, a₁ = -1)
+        @test dot(gs, x) ≈ predicted atol = 1.0e-6
+        @test rel_residual(H, x, gs, z, -1) < 1.0e-6
+    end
+
+    @testset "resolvent, two-site (adaptive χ), formulation $flav" for flav in
+        (Galerkin(), LeastSquares())
+        solver = flav isa LeastSquares ? CG(; tol = 1.0e-12) : GMRES(; tol = 1.0e-12)
+        alg = DMRGSolve2(;
+            formulation = flav, solver, trunc = truncrank(16), tol = 1.0e-10, verbosity = 0
+        )
+        x, = linear_solve(complex(copy(gs)), H, gs, alg; a₀ = z, a₁ = -1)
+        @test dot(gs, x) ≈ predicted atol = 1.0e-5
+        @test rel_residual(H, x, gs, z, -1) < 1.0e-5
+    end
+
+    @testset "matches propagator (correction vector)" begin
+        alg_ls = DMRGSolve(; solver = GMRES(; tol = 1.0e-12), tol = 1.0e-10, verbosity = 0)
+        x, = linear_solve(complex(copy(gs)), H, gs, alg_ls; a₀ = z, a₁ = -1)
+        for f in (NaiveInvert(), Jeckelmann())
+            g_prop, = propagator(gs, z, H, DynamicalDMRG(; flavour = f, tol = 1.0e-10, verbosity = 0))
+            @test dot(gs, x) ≈ g_prop atol = 1.0e-5
+        end
+    end
+
+    @testset "adaptive tolerances (incl. BiCGStab)" begin
+        # DynamicTol works for any solver, including short-recurrence BiCGStab (retunes only `tol`)
+        alg_bicg = DMRGSolve(;
+            solver = DynamicTol(BiCGStab(; tol = 1.0e-12); tol_min = 1.0e-14, tol_max = 1.0e-4),
+            tol = 1.0e-9, verbosity = 0
+        )
+        x, = linear_solve(complex(copy(gs)), H, gs, alg_bicg; a₀ = z, a₁ = -1)
+        @test dot(gs, x) ≈ predicted atol = 1.0e-5
+        @test rel_residual(H, x, gs, z, -1) < 1.0e-5
+
+        # adaptive-by-default keyword path (DynamicTol-wrapped GMRES)
+        x2, envs, ϵ = linear_solve(complex(copy(gs)), H, gs; a₀ = z, a₁ = -1, tol = 1.0e-9, verbosity = 0)
+        @test dot(gs, x2) ≈ predicted atol = 1.0e-5
+        @test rel_residual(H, x2, gs, z, -1) < 1.0e-5
+    end
+
+    @testset "positive-definite solve via keyword interface" begin
+        # s ≫ spectrum ⇒ (s − H) is positive-definite; the CG path is selected by `isposdef`
+        s = E₀ + 50.0
+        # eigenstate RHS: exact solution 1/(s − E₀) · gs, real and positive
+        x, envs, ϵ = linear_solve(copy(gs), H, gs; a₀ = s, a₁ = -1, isposdef = true, tol = 1.0e-10, verbosity = 0)
+        @test dot(gs, x) ≈ 1 / (s - E₀) atol = 1.0e-8
+        @test rel_residual(H, x, gs, s, -1) < 1.0e-8
+    end
+
+    @testset "general RHS, plain A·x = b against dense" begin
+        # shift H to a well-conditioned, invertible operator A = s·I − H (s ≫ spectrum)
+        s = E₀ + 50.0
+        b = normalize!(FiniteMPS(randn, ComplexF64, L, ℙ^2, ℙ^8))
+        alg = DMRGSolve2(;
+            solver = CG(; tol = 1.0e-12), trunc = truncrank(64), tol = 1.0e-9, verbosity = 0
+        ) & DMRGSolve(; solver = CG(; tol = 1.0e-12), tol = 1.0e-9, verbosity = 0)
+        x, = linear_solve(complex(copy(b)), H, b, alg; a₀ = s, a₁ = -1)
+        @test rel_residual(H, x, b, s, -1) < 1.0e-4
+    end
+end
+
+@testset "linear_solve WindowMPS" verbose = true begin
+    # `LeastSquares` builds the `A†b` term from the mixed sandwich ⟨x|A|b⟩, which is the only way to
+    # reach it for a window operator (there is no `*(::WindowMPOHamiltonian, ::WindowMPS)`)
+    N = 10
+    H = transverse_field_ising(; g = -4)
+    Ω, = find_groundstate(InfiniteMPS(ℂ^2, ℂ^10), H, VUMPS(; verbosity = 0))
+    XΩ = WindowMPS(Ω, N)
+    H_w = WindowMPOHamiltonian(H, 1:N)
+
+    E₀ = expectation_value(XΩ, H_w)
+    z = E₀ + 0.5 + 0.3im
+    predicted = 1 / (z - E₀)
+
+    @testset "resolvent, formulation $flav" for flav in (Galerkin(), LeastSquares())
+        solver = flav isa LeastSquares ? CG(; tol = 1.0e-12) : GMRES(; tol = 1.0e-12)
+        alg = DMRGSolve(; formulation = flav, solver, tol = 1.0e-10, verbosity = 0)
+        x, = linear_solve(XΩ, H_w, XΩ, alg; a₀ = z, a₁ = -1)
+        @test dot(XΩ, x) ≈ predicted atol = 1.0e-6
+    end
+end

From 0cea74c26c71e3c0e249689be2454d90411758d5 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Mon, 12 Oct 2026 16:43:15 -0400
Subject: [PATCH 5/5] refactor: propagator/DynamicalDMRG on top of linear_solve

`corvector.jl` carried its own DMRG sweep iterator, which had drifted from the one in the
`linear_solve` category on three counts: it converged on the per-site update norm
||AC' - AC|| (which cannot tell converged from stalled) rather than on the
residual of the linear system, it used a fixed-tolerance GMRES rather than the
structure-selected adaptive solver, and it was single-site only. `propagator` is
now a single method for both flavours whose body is a `linear_solve!` call.

Both flavours keep their variational problem unchanged. `NaiveInvert` is a
`Galerkin` solve of `(z - H)x = |psi0>`. `Jeckelmann` still minimizes functional
(14) of Jeckelmann2002 and reconstructs G(z) through equation (11): its operator
is assembled as `LinearCombination((H^2, H), (1, -2w))` over a hand-built
`LazyLincoCache` -- hand-built so that `squaredenvs`' boundary conditions, which
matter for `WindowMPS`, survive -- and its stationarity condition
`((w - H)^2 + eta^2) psi = -eta |psi0>` is scaled by `-1/eta` so that the
right-hand side is `|psi0>` itself. No scaled copy of the right-hand side is
needed and `x` is the identical vector; the shift comes out real, so a real
initial guess stays real.

Consequent changes to the defaults of `DynamicalDMRG`:

- `tol` defaults to `Defaults.tol` rather than `10 * Defaults.tol`, and now
  measures the relative residual of the linear system.
- `solver` defaults to `Defaults.alg_linsolve()` (adaptive tolerances).
- a new `trunc` field prepends a truncated two-site sweep, so the bond
  dimension no longer has to be fixed by the initial guess.
- the `backend` field is forwarded to the underlying `DMRGSolve`/`DMRGSolve2`.
- iteration logs are named after the underlying `DMRGSolve`/`DMRGSolve2`
  instead of `DDMRG`.
- the `Jeckelmann` flavour throws for real `z` instead of returning `NaN`; the
  equation-(11) reconstruction divides by `imag(z)`.

The docstrings claimed that the second return value is `1/(z - H)|psi0>` for
both flavours. That is only true for `NaiveInvert`; for `Jeckelmann` it is the
minimizer of the functional, i.e. the imaginary part of the propagator. Say so.

Adds a `linear_solve` section to the algorithms manual: the category was exported
but appeared nowhere in the docs.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 docs/src/changelog.md                         |   5 +
 docs/src/man/algorithms.md                    |  28 +-
 src/algorithms/propagator/corvector.jl        | 283 ++++++------------
 .../{linear_solve.jl => propagator.jl}        |  69 ++++-
 test/groundstate/dynamical_dmrg.jl            |  62 ----
 5 files changed, 188 insertions(+), 259 deletions(-)
 rename test/algorithms/{linear_solve.jl => propagator.jl} (67%)
 delete mode 100644 test/groundstate/dynamical_dmrg.jl

diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 1e8ea3fa9..226a27b5d 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -29,6 +29,9 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - `changebonds(::InfiniteMPOHamiltonian, ::SvdCut)` compresses infinite Hamiltonians with the
   almost-Schmidt decomposition of [Parker et al.](https://arxiv.org/abs/1909.06341), and
   `left_canonicalize!`/`right_canonicalize!` bring an `InfiniteMPOHamiltonian` into canonical form. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539))
+- `linear_solve`/`linear_solve!` for solving `(a₀ + a₁·A)·x = b` with finite MPS `x` and `b`, using the
+  DMRG-style sweeps `DMRGSolve` and `DMRGSolve2`; see [`linear_solve`](@ref) for details.
+- `DynamicalDMRG` gained a `trunc` field for a bond-adaptive `propagator`; see [`DynamicalDMRG`](@ref).
 
 ### Changed
 
@@ -46,6 +49,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   its allocator across steps, rather than copying the state at every step. A `finalize` callback
   that keeps the state of an earlier step therefore sees it updated by later steps, as with
   `time_evolve!`. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
+- `propagator`/`DynamicalDMRG` are now implemented on top of `linear_solve`, which changes the default
+  convergence criterion, `tol` and local solver; `Jeckelmann` now throws for real `z`.
 
 ### Deprecated
 
diff --git a/docs/src/man/algorithms.md b/docs/src/man/algorithms.md
index df90b7fd8..dfad8dfca 100644
--- a/docs/src/man/algorithms.md
+++ b/docs/src/man/algorithms.md
@@ -444,6 +444,28 @@ state, by a new state.
 approximate
 ```
 
+## `linear_solve`
+
+Solving a linear system ``(a₀ + a₁ A) x = b`` for an MPS ``x``, with ``A`` an MPO and ``b`` an MPS.
+The shift ``a₀ + a₁ A`` follows the convention of `KrylovKit.linsolve` and is applied implicitly, so
+no shifted operator is ever formed.
+The sweeps are DMRG-like: [`DMRGSolve`](@ref) keeps the bond dimension of the initial guess fixed,
+while [`DMRGSolve2`](@ref) updates two sites at a time and is therefore rank-adaptive.
+Each is parameterized by how the local subproblem is posed — [`Galerkin`](@ref) solves the effective
+system directly, [`LeastSquares`](@ref) minimizes ``\|(a₀ + a₁ A) x - b\|²`` through the normal
+equations.
+
+The resolvent is the canonical application: ``\frac{1}{z - H}|ψ₀⟩`` is `linear_solve(ψ₀, H, ψ₀; a₀ = z, a₁ = -1)`.
+
+```@docs; canonical=false
+linear_solve
+linear_solve!
+DMRGSolve
+DMRGSolve2
+Galerkin
+LeastSquares
+```
+
 ## Varia
 
 What follows is a medley of lesser known (or used) algorithms and don't entirely fit under
@@ -451,9 +473,9 @@ one of the above categories.
 
 ### Dynamical DMRG
 
-Dynamical DMRG has been described in other papers and is a way to find the propagator. The
-basic idea is that to calculate ``G(z) = ⟨ V | (H-z)^{-1} | V ⟩ `` , one can variationally
-find ``(H-z) |W ⟩ = | V ⟩ `` and then the propagator simply equals ``G(z) = ⟨ V | W ⟩``.
+Dynamical DMRG has been described in other papers and is a way to find the propagator.
+The basic idea is that to calculate ``G(z) = ⟨ V | (H-z)^{-1} | V ⟩ `` , one can variationally find ``(H-z) |W ⟩ = | V ⟩ `` and then the propagator simply equals ``G(z) = ⟨ V | W ⟩``.
+This is a thin wrapper around [`linear_solve`](@ref), see there for the sweep, the convergence criterion and the local solver options.
 
 ```@docs; canonical=false
 propagator
diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl
index 0704743e6..6f2f63a91 100644
--- a/src/algorithms/propagator/corvector.jl
+++ b/src/algorithms/propagator/corvector.jl
@@ -8,54 +8,6 @@ abstract type DDMRG_Flavour end
 """
 $(TYPEDEF)
 
-A dynamical DMRG method for calculating dynamical properties and excited states, based on a
-variational principle for dynamical correlation functions.
-
-# Fields
-
-$(TYPEDFIELDS)
-
-# See also
-
-Used as the `algorithm` argument of [`propagator`](@ref).
-
-# References
-
-* [Jeckelmann. Phys. Rev. B 66 (2002)](@cite jeckelmann2002)
-"""
-@kwdef struct DynamicalDMRG{F <: DDMRG_Flavour, S, B} <: Algorithm
-    "flavour of the algorithm to use, either of type [`NaiveInvert`](@ref) or [`Jeckelmann`](@ref)"
-    flavour::F = NaiveInvert()
-    "algorithm used for the linear solvers"
-    solver::S = Defaults.linearsolver
-    "convergence tolerance on the largest change of a center tensor over a sweep"
-    tol::Float64 = Defaults.tol * 10
-    "maximal amount of iterations"
-    maxiter::Int = Defaults.maxiter
-    "setting for how much information is displayed"
-    verbosity::Int = Defaults.verbosity
-    "backend for tensor contractions and index manipulations"
-    backend::B = Defaults.backend()
-end
-
-IterativeLoggers.IterLog(::DynamicalDMRG) = IterLog("DDMRG")
-
-"""
-    propagator(ψ₀::AbstractFiniteMPS, z::Number, H::MPOHamiltonian, alg::DynamicalDMRG; init = copy(ψ₀)) -> (g, ψ)
-
-Calculate the action of the propagator ``\\frac{1}{z - H}|ψ₀⟩`` using the dynamical DMRG
-algorithm.
-
-# Returns
-
-- `g`: approximation of the propagator matrix element ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩``
-- `ψ`: MPS approximation of ``\\frac{1}{z - H}|ψ₀⟩``
-"""
-function propagator end
-
-"""
-$(TYPEDEF)
-
 An alternative approach to the dynamical DMRG algorithm, without quadratic terms but with a
 less controlled approximation.
 This algorithm minimizes the following cost function
@@ -71,100 +23,6 @@ Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac
 """
 struct NaiveInvert <: DDMRG_Flavour end
 
-# Internal state of the dynamical DMRG sweeps
-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)
-    )
-    allocator = default_allocator(A, SerialScheduler())
-    h_envs = environments(init, H, init) # environments for h
-    mixedenvs = environments(init, A) # environments for <init | A>
-
-    state = DDMRGState(init, A, z, H, (h_envs, mixedenvs), 0, 2 * alg.tol, allocator)
-    _propagator_sweeps!(alg, state)
-
-    return dot(A, init), init
-end
-
-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.ϵ
-
-    for i in (direction === Val(:right) ? fwd : bwd)
-        tos = AC_projection(i, init, A, mixedenvs; alg.backend, allocator)
-
-        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))
-
-        ϵ = max(ϵ, norm(AC′ - AC))
-        init.AC[i] = AC′
-
-        convhist.converged == 0 &&
-            @warn "propagator ($i) failed to converge: normres = $(convhist.normres)"
-    end
-
-    return DDMRGState(
-        init, A, z, operator, state.envs, state.iter, ϵ, allocator,
-    )
-end
-
 """
 $(TYPEDEF)
 
@@ -178,10 +36,12 @@ which attains its minimum at
 ((ω - H)^2 + η^2)|ψ⟩ = -η|ψ₀⟩
 ```
 
-Together with equation (11) from that same paper we can determine the full propagator
-``\\frac{1}{z - H}|ψ₀⟩``.
+The solution of that equation is the imaginary part of the propagator; together with equation (11)
+from that same paper it determines the full ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩``. Because of that
+reconstruction step this flavour requires ``η = \\mathrm{Im}(z) ≠ 0``.
 
-Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac{1}{z - H}|ψ₀⟩``.
+Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and the minimizer ``|ψ⟩`` of the functional
+above (*not* ``\\frac{1}{z - H}|ψ₀⟩`` itself).
 
 # See also
 
@@ -193,63 +53,106 @@ Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac
 """
 struct Jeckelmann <: DDMRG_Flavour end
 
-function propagator(
-        A::AbstractFiniteMPS, z::Number, H,
-        alg::DynamicalDMRG{Jeckelmann}; init = copy(A)
-    )
-    allocator = default_allocator(A, SerialScheduler())
-    ω = real(z)
-    η = imag(z)
+"""
+$(TYPEDEF)
 
-    envs1 = environments(init, H, init) # environments for h
-    H2, envs2 = squaredenvs(init, H, envs1) # environments for h^2
-    mixedenvs = environments(init, A) # environments for <init | A>
+A dynamical DMRG method for calculating dynamical properties and excited states, based on a
+variational principle for dynamical correlation functions.
 
-    state = DDMRGState(init, A, z, (H, H2), (envs1, envs2, mixedenvs), 0, 2 * alg.tol, allocator)
-    _propagator_sweeps!(alg, state)
+This is a thin wrapper around [`linear_solve`](@ref): the sweep, the convergence criterion and the
+adaptive local tolerances are those of [`DMRGSolve`](@ref) / [`DMRGSolve2`](@ref).
 
-    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)
-    b = zero(a)
-    for i in 1:length(cb)
-        b += @plansor cb[i][1 2; 3] * init.C[end][3; 4] *
-            rightenv(envs1, length(A), A)[i][4 2; 5] * conj(A.C[end][1; 5])
-    end
+# Fields
+
+$(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`propagator`](@ref).
+
+# References
 
-    v = b / η - ω / η * a + 1im * a
-    return v, init
+* [Jeckelmann. Phys. Rev. B 66 (2002)](@cite jeckelmann2002)
+"""
+@kwdef struct DynamicalDMRG{F <: DDMRG_Flavour, S, T, B} <: Algorithm
+    "flavour of the algorithm to use, either of type [`NaiveInvert`](@ref) or [`Jeckelmann`](@ref)"
+    flavour::F = NaiveInvert()
+    "local linear solver; a plain KrylovKit solver, or one wrapped in `DynamicTol` for per-bond adaptive tolerances (the default)"
+    solver::S = Defaults.alg_linsolve()
+    "tolerance for convergence criterium, measured as the relative residual of the linear system"
+    tol::Float64 = Defaults.tol
+    "maximal amount of iterations"
+    maxiter::Int = Defaults.maxiter
+    "setting for how much information is displayed"
+    verbosity::Int = Defaults.verbosity
+    "if supplied, a truncated two-site sweep ([`DMRGSolve2`](@ref)) is prepended to adapt the bond dimension"
+    trunc::T = nothing
+    "backend for tensor contractions and index manipulations"
+    backend::B = Defaults.backend()
 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)
+# `Jeckelmann` solves `((ω - H)² + η²)·x = -η|ψ₀⟩`, rescaled by `-1/η` so that `b = |ψ₀⟩`
+propagator_shift(::NaiveInvert, z) = (z, -one(z))
+function propagator_shift(::Jeckelmann, 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)
+    iszero(η) && throw(
+        ArgumentError(
+            "`Jeckelmann` requires `imag(z) != 0`; use `NaiveInvert` flavour for real `z`"
+        )
+    )
+    return (-abs2(z) / η, -inv(η))
+end
 
-        ϵ = max(ϵ, norm(AC′ - init.AC[i]))
-        init.AC[i] = AC′
+propagator_operator(::NaiveInvert, H, x, envs, z) = (H, envs)
+function propagator_operator(::Jeckelmann, H, x, envs, z)
+    H², envs² = squaredenvs(x, H, envs)
+    A = LinearCombination((H², H), (one(real(z)), -2 * real(z)))
+    return A, LazyLincoCache(A, (envs², envs))
+end
 
-        convhist.converged == 0 &&
-            @warn "propagator ($i) failed to converge: normres $(convhist.normres)"
+# `G(z)` from the solution vector
+propagator_value(::NaiveInvert, x, ψ₀, z, H, envs) = dot(ψ₀, x)
+function propagator_value(::Jeckelmann, x, ψ₀, z, H, envs)
+    ω, η = real(z), imag(z)
+    # equation (11) of Jeckelmann2002: the real part of `G(z)` follows from ⟨ψ₀|H|x⟩
+    a = dot(ψ₀, x)
+    cb = leftenv(envs, 1, ψ₀) * TransferMatrix(x.AL, H[1:length(ψ₀.AL)], ψ₀.AL)
+    b = zero(a)
+    for i in 1:length(cb)
+        b += @plansor cb[i][1 2; 3] * x.C[end][3; 4] *
+            rightenv(envs, length(ψ₀), ψ₀)[i][4 2; 5] * conj(ψ₀.C[end][1; 5])
     end
+    return b / η - ω / η * a + 1im * a
+end
+
+"""
+    propagator(ψ₀::AbstractFiniteMPS, z::Number, H, alg::DynamicalDMRG; init = ψ₀) -> (g, ψ)
+
+Calculate the action of the propagator ``\\frac{1}{z - H}|ψ₀⟩`` using the dynamical DMRG
+algorithm.
+
+# Returns
 
-    return DDMRGState(
-        init, A, z, state.operator, state.envs, state.iter, ϵ, allocator,
+- `g`: approximation of the propagator matrix element ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩``
+- `ψ`: for [`NaiveInvert`](@ref), the MPS approximation of ``\\frac{1}{z - H}|ψ₀⟩``; for
+  [`Jeckelmann`](@ref), the vector its functional optimizes,
+  ``-η[(ω - H)^2 + η^2]^{-1}|ψ₀⟩``, i.e. the imaginary part of the propagator, from which `g` is
+  reconstructed.
+
+`init` is used as the initial guess and is left untouched. The underlying variational problem is
+solved with [`linear_solve`](@ref); for full control over the sweep, call that directly, e.g.
+`linear_solve(ψ₀, H, ψ₀; a₀ = z, a₁ = -1)` for the [`NaiveInvert`](@ref) flavour.
+"""
+function propagator(
+        ψ₀::AbstractFiniteMPS, z::Number, H, alg::DynamicalDMRG; init = ψ₀
     )
+    a₀, a₁ = propagator_shift(alg.flavour, z)
+    x = promote_to_shift(copy(init), a₀, a₁)
+    envs = environments(x, H, x)
+    A, Aenvs = propagator_operator(alg.flavour, H, x, envs, z)
+    sweep = linear_solve_algorithm(; alg.solver, alg.tol, alg.maxiter, alg.verbosity, alg.trunc, alg.backend)
+    x, = linear_solve!(x, A, ψ₀, sweep, Aenvs; a₀, a₁)
+    return propagator_value(alg.flavour, x, ψ₀, z, H, envs), x
 end
 
 function squaredenvs(
diff --git a/test/algorithms/linear_solve.jl b/test/algorithms/propagator.jl
similarity index 67%
rename from test/algorithms/linear_solve.jl
rename to test/algorithms/propagator.jl
index 663122057..eb0173229 100644
--- a/test/algorithms/linear_solve.jl
+++ b/test/algorithms/propagator.jl
@@ -1,6 +1,6 @@
 println("
 -----------------------------
-|   Linear solver tests      |
+|   Propagator tests         |
 -----------------------------
 ")
 
@@ -40,7 +40,8 @@ end
         solver = flav isa LeastSquares ? CG(; tol = 1.0e-12) : GMRES(; tol = 1.0e-12)
         alg = DMRGSolve(; formulation = flav, solver, tol = 1.0e-10, verbosity = 0)
         # solve (z − H) x = gs  ⟺  a₀ = z, a₁ = −1
-        x, envs, ϵ = linear_solve(complex(copy(gs)), H, gs, alg; a₀ = z, a₁ = -1)
+        x, envs, info = linear_solve(complex(copy(gs)), H, gs, alg; a₀ = z, a₁ = -1)
+        @test info.converged && convergence_measure(info) == info.galerkin
         @test dot(gs, x) ≈ predicted atol = 1.0e-6
         @test rel_residual(H, x, gs, z, -1) < 1.0e-6
     end
@@ -76,7 +77,7 @@ end
         @test rel_residual(H, x, gs, z, -1) < 1.0e-5
 
         # adaptive-by-default keyword path (DynamicTol-wrapped GMRES)
-        x2, envs, ϵ = linear_solve(complex(copy(gs)), H, gs; a₀ = z, a₁ = -1, tol = 1.0e-9, verbosity = 0)
+        x2, = linear_solve(complex(copy(gs)), H, gs; a₀ = z, a₁ = -1, tol = 1.0e-9, verbosity = 0)
         @test dot(gs, x2) ≈ predicted atol = 1.0e-5
         @test rel_residual(H, x2, gs, z, -1) < 1.0e-5
     end
@@ -85,7 +86,7 @@ end
         # s ≫ spectrum ⇒ (s − H) is positive-definite; the CG path is selected by `isposdef`
         s = E₀ + 50.0
         # eigenstate RHS: exact solution 1/(s − E₀) · gs, real and positive
-        x, envs, ϵ = linear_solve(copy(gs), H, gs; a₀ = s, a₁ = -1, isposdef = true, tol = 1.0e-10, verbosity = 0)
+        x, = linear_solve(copy(gs), H, gs; a₀ = s, a₁ = -1, isposdef = true, tol = 1.0e-10, verbosity = 0)
         @test dot(gs, x) ≈ 1 / (s - E₀) atol = 1.0e-8
         @test rel_residual(H, x, gs, s, -1) < 1.0e-8
     end
@@ -122,3 +123,63 @@ end
         @test dot(XΩ, x) ≈ predicted atol = 1.0e-6
     end
 end
+
+@testset "Dynamical DMRG FiniteMPS" verbose = true begin
+    L = 10
+    H = force_planar(-transverse_field_ising(; L, g = -4))
+    gs, = find_groundstate(FiniteMPS(L, ℙ^2, ℙ^10), H; verbosity = verbosity_conv)
+    E₀ = expectation_value(gs, H)
+
+    vals = (-0.5:0.2:0.5) .+ E₀
+    eta = 0.3im
+
+    predicted = [1 / (v + eta - E₀) for v in vals]
+
+    @testset "Flavour $f" for f in (Jeckelmann(), NaiveInvert())
+        alg = DynamicalDMRG(; flavour = f, verbosity = 0, tol = 1.0e-8)
+        data = map(vals) do v
+            result, = propagator(gs, v + eta, H, alg)
+            return result
+        end
+        @test data ≈ predicted atol = 1.0e-8
+    end
+
+    # `trunc` prepends a two-site sweep, so the bond dimension of the initial guess need not
+    # already be large enough: start from χ = 2 and let it grow back to the χ = 10 of the solution
+    @testset "trunc (adaptive χ), flavour $f" for f in (Jeckelmann(), NaiveInvert())
+        alg = DynamicalDMRG(;
+            flavour = f, verbosity = 0, tol = 1.0e-8, trunc = truncrank(16)
+        )
+        init = FiniteMPS(L, ℙ^2, ℙ^2)
+        data = map(vals) do v
+            result, = propagator(gs, v + eta, H, alg; init)
+            return result
+        end
+        @test data ≈ predicted atol = 1.0e-6
+    end
+end
+
+@testset "Dynamical DMRG WindowMPS" verbose = true begin
+    N = 20
+
+    H = transverse_field_ising(g = -4)
+    Ω = InfiniteMPS(ComplexSpace(2), ComplexSpace(20))
+
+    (Ω, _) = find_groundstate(Ω, H, VUMPS(verbosity = verbosity_conv))
+    XΩ = WindowMPS(Ω, N)
+    H_w = WindowMPOHamiltonian(H, 1:N)
+
+    gs_en = expectation_value(XΩ, H_w)
+    vals = range(gs_en - 1.0, gs_en + 1.0, length = 5)
+    eta = 0.3im
+    predicted = [1 / (v + eta - gs_en) for v in vals]
+
+    @testset "Flavour $f" for f in (Jeckelmann(), NaiveInvert())
+        alg = DynamicalDMRG(; flavour = f, verbosity = 0, tol = 1.0e-8)
+        data = map(vals) do v
+            result, = propagator(XΩ, v + eta, H_w, alg)
+            return result
+        end
+        @test data ≈ predicted atol = 1.0e-8
+    end
+end
diff --git a/test/groundstate/dynamical_dmrg.jl b/test/groundstate/dynamical_dmrg.jl
deleted file mode 100644
index 8e07757f4..000000000
--- a/test/groundstate/dynamical_dmrg.jl
+++ /dev/null
@@ -1,62 +0,0 @@
-println("
------------------------------
-|   Dynamical DMRG tests     |
------------------------------
-")
-
-using .TestSetup
-using Test, TestExtras
-using MPSKit
-using TensorKit
-using TensorKit: ℙ
-
-verbosity_conv = 1
-
-@testset "Dynamical DMRG FiniteMPS" verbose = true begin
-    L = 10
-    H = force_planar(-transverse_field_ising(; L, g = -4))
-    gs, = find_groundstate(FiniteMPS(L, ℙ^2, ℙ^10), H; verbosity = verbosity_conv)
-    E₀ = expectation_value(gs, H)
-
-    vals = (-0.5:0.2:0.5) .+ E₀
-    eta = 0.3im
-
-    predicted = [1 / (v + eta - E₀) for v in vals]
-
-    @testset "Flavour $f" for f in (Jeckelmann(), NaiveInvert())
-        alg = DynamicalDMRG(; flavour = f, verbosity = 0, tol = 1.0e-8)
-        data = map(vals) do v
-            result, = propagator(gs, v + eta, H, alg)
-            return result
-        end
-        @test data ≈ predicted atol = 1.0e-8
-    end
-end
-
-
-@testset "Dynamical DMRG WindowMPS" verbose = true begin
-    N = 20
-
-    H = transverse_field_ising(g = -4)
-    Ω = InfiniteMPS(ComplexSpace(2), ComplexSpace(20))
-
-    (Ω, _) = find_groundstate(Ω, H, VUMPS(verbosity = verbosity_conv))
-    XΩ = WindowMPS(Ω, N)
-    H_w = WindowMPOHamiltonian(H, 1:N)
-
-
-    gs_en = expectation_value(XΩ, H_w)
-    vals = range(gs_en - 1.0, gs_en + 1.0, length = 5)
-    eta = 0.3im
-    predicted = [1 / (v + eta - gs_en) for v in vals]
-
-
-    @testset "Flavour $f" for f in (Jeckelmann(), NaiveInvert())
-        alg = DynamicalDMRG(; flavour = f, verbosity = 0, tol = 1.0e-8)
-        data = map(vals) do v
-            result, = propagator(XΩ, v + eta, H_w, alg)
-            return result
-        end
-        @test data ≈ predicted atol = 1.0e-8
-    end
-end