From 91892f3f16db50616594a19ce17c13ef08314a9f Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 07:45:08 -0400 Subject: [PATCH 01/10] Refactor finite DMRG into a complete-sweep iterator --- src/algorithms/groundstate/dmrg.jl | 135 +++++++++++++++++------------ 1 file changed, 79 insertions(+), 56 deletions(-) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index b2400fdfa..2b3e8c9c7 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, 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 = it.state) + 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(string(nameof(typeof(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...) From 8e65db7e36347eea718c722bb9afb7faba704ec7 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 18:04:14 -0400 Subject: [PATCH 02/10] Refactor finite DMRG approximation into a complete-sweep iterator Co-Authored-By: Claude Opus 5.5 --- src/algorithms/approximate/fvomps.jl | 112 +++++++++++++-------------- 1 file changed, 56 insertions(+), 56 deletions(-) diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl index 465b40f1b..5a17d7a22 100644 --- a/src/algorithms/approximate/fvomps.jl +++ b/src/algorithms/approximate/fvomps.jl @@ -1,74 +1,74 @@ -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(string(nameof(typeof(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 + ϵ = approximate_sweep!(state.mps, state.operator, it.alg, state.envs, state.allocator) + ψ, 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.allocator) + return (ψ, envs, ϵ), it.state +end - ψ.AC[pos] = AC′ - end +function approximate_sweep!(ψ, Oϕ, alg::DMRG2, envs, allocator) + ϵ = 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)) - # 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 ϵ +end - return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, localchange = ϵ, numiter = iter) +function approximate_sweep!(ψ, Oϕ, alg::DMRG, envs, allocator) + ϵ = 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′)) + + ψ.AC[pos] = AC′ + end + return ϵ end From 9195042d6f0a0e5908163679b067ea6d60f9b113 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 19:14:38 -0400 Subject: [PATCH 03/10] Refactor dynamical DMRG into a complete-sweep iterator Co-Authored-By: Claude Opus 5.5 --- src/algorithms/propagator/corvector.jl | 150 +++++++++++++++---------- 1 file changed, 90 insertions(+), 60 deletions(-) diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl index bd72f892c..c73518f96 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -69,6 +69,49 @@ 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 + +function Base.iterate(it::IterativeSolver{<:DynamicalDMRG}, state::DDMRGState) + ϵ = ddmrg_sweep!(state, it.alg) + it.state = DDMRGState( + state.mps, state.target, state.z, state.operator, state.envs, state.iter + 1, ϵ, + state.allocator, + ) + return (state.mps, ϵ), it.state +end + +function _propagator_sweeps!(alg::DynamicalDMRG, state::DDMRGState) + log = IterLog("DDMRG") + 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 +120,33 @@ function propagator( h_envs = environments(init, H, init) # environments for h mixedenvs = environments(init, A) # environments for - ϵ = 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 ddmrg_sweep!(state::DDMRGState, alg::DynamicalDMRG{NaiveInvert}) + (; mps, target, z, operator, allocator) = state + init, A, H = mps, target, operator + h_envs, mixedenvs = state.envs + ϵ = 0.0 - 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 [1:(length(A) - 1); length(A):-1:2] + 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 ϵ end """ @@ -154,39 +189,8 @@ function propagator( H2, envs2 = squaredenvs(init, H, envs1) # environments for h^2 mixedenvs = environments(init, A) # environments for - ϵ = 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 +204,32 @@ function propagator( return v, init end +function ddmrg_sweep!(state::DDMRGState, alg::DynamicalDMRG{Jeckelmann}) + (; mps, target, z, allocator) = state + init, A = mps, target + H, H2 = state.operator + envs1, envs2, mixedenvs = state.envs + ω = real(z) + η = imag(z) + ϵ = 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 + + return ϵ +end + function squaredenvs( state::AbstractFiniteMPS, H, envs = environments(state, H, state) ) From f9b1b582f39f6f0737e9ddd8df0b46ccd82459de Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 19:16:39 -0400 Subject: [PATCH 04/10] Refactor multiline IDMRG approximation into a sweep iterator Co-Authored-By: Claude Opus 5.5 --- src/algorithms/approximate/idmrg.jl | 337 ++++++++++++++-------------- 1 file changed, 170 insertions(+), 167 deletions(-) diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl index 8b8079379..45fa513ab 100644 --- a/src/algorithms/approximate/idmrg.jl +++ b/src/algorithms/approximate/idmrg.jl @@ -1,199 +1,202 @@ +# Internal state of the multiline IDMRG/IDMRG2 approximation, where `ϵ` is the change of the +# center bond tensor over the last sweep. IDMRG2 reuses one `truncation_errors` matrix across sweeps. +struct IDMRGApproximateState{S, O, E, V, A} + mps::S + operator::O + envs::E + iter::Int + ϵ::Float64 + truncation_errors::V + allocator::A +end + 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(string(nameof(typeof(alg)))) + ϵ_truncs = alg isa IDMRG2 ? + PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ))) : nothing + state = IDMRGApproximateState(ψ, toapprox, envs, 0, 2 * alg.tol, ϵ_truncs, 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::IDMRGApproximateState) + ϵ = approximate_sweep!( + state.mps, state.operator, it.alg, state.envs, state.allocator, state.truncation_errors ) - 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(ψ))) + it.state = IDMRGApproximateState( + state.mps, state.operator, state.envs, state.iter + 1, ϵ, + state.truncation_errors, state.allocator, + ) + 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 approximate_sweep!( + ψ::MultilineMPS, toapprox, alg::IDMRG, envs, allocator, ::Nothing + ) + 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 - # 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] + # 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) - 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 + return norm(C_current - ψ.C[:, 0]) +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) +function approximate_sweep!( + ψ::MultilineMPS, toapprox, alg::IDMRG2, envs, allocator, ϵ_truncs + ) + 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 - normalize!(envs, ψ, toapprox) + transfer_leftenv!(envs, ψ, toapprox, site + 1) + transfer_rightenv!(envs, ψ, toapprox, site) + end - # 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 + # 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] - 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 + # 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 - # 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) + + # update error + return 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 end From bf1c9989f9e08b0717df5f7a334caadef5601ab5 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 19:16:40 -0400 Subject: [PATCH 05/10] Refactor IDMRG leading_boundary into a sweep iterator Co-Authored-By: Claude Opus 5.5 --- src/algorithms/statmech/idmrg.jl | 327 ++++++++++++++++--------------- 1 file changed, 169 insertions(+), 158 deletions(-) diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 8e5a9b0f0..e8161aaf2 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -1,201 +1,212 @@ +# Internal state of the IDMRG/IDMRG2 leading boundary search, where `ϵ` is the change of the +# center bond tensor over the last sweep. IDMRG2 reuses one `truncation_errors` matrix across sweeps. +struct IDMRGBoundaryState{S, O, E, V, A} + mps::S + operator::O + envs::E + iter::Int + ϵ::Float64 + truncation_errors::V + allocator::A +end + 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(string(nameof(typeof(alg)))) + ϵ_truncs = alg isa IDMRG2 ? + PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ))) : nothing + state = IDMRGBoundaryState(ψ, operator, envs, 0, 2 * alg.tol, ϵ_truncs, 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::IDMRGBoundaryState) + iter = state.iter + 1 + alg_eigsolve = adapt_solver(it.alg_eigsolve; iter, g_global = state.ϵ) + ϵ = leading_boundary_sweep!( + state.mps, state.operator, it.alg, state.envs, alg_eigsolve, state.allocator, + state.truncation_errors, ) - 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(ψ))) + it.state = IDMRGBoundaryState( + state.mps, state.operator, state.envs, iter, ϵ, state.truncation_errors, state.allocator, + ) + 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 leading_boundary_sweep!( + ψ::InfiniteMultilineMPS, operator, alg::IDMRG, envs, alg_eigsolve, allocator, ::Nothing + ) + C_current = ψ.C[:, 0] - normalize!(envs, ψ, operator, ψ) + # 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) - # 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) + ac = ψ.AC[row, col] + (col == width(ψ)) && (ac = copy(ac)) # needed in next sweep + ψ.AL[row, col], ψ.C[row, col] = left_orth!(ac) + end - 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) + transfer_leftenv!(envs, ψ, operator, ψ, col + 1) + end - ψ.AL[row + 1, site] = al - ψ.C[row + 1, site] = complex(c) - ψ.AR[row + 1, site + 1] = _transpose_front(ar) + # 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) - ψ.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 + 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 - # TODO: decide if we should compare at the half-sweep level? - # C_current = ψ.C[:, site] + transfer_rightenv!(envs, ψ, operator, ψ, col - 1) + 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) + return norm(C_current - ψ.C[:, 0]) +end - 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) +function leading_boundary_sweep!( + ψ::InfiniteMultilineMPS, operator, alg::IDMRG2, envs, alg_eigsolve, allocator, ϵ_truncs + ) + 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 - ψ.AL[row + 1, site] = al - ψ.C[row + 1, site] = complex(c) - ψ.AR[row + 1, site + 1] = _transpose_front(ar) - end + transfer_leftenv!(envs, ψ, operator, ψ, site + 1) + transfer_rightenv!(envs, ψ, operator, ψ, site) + end - transfer_leftenv!(envs, ψ, operator, ψ, site + 1) - transfer_rightenv!(envs, ψ, operator, ψ, site) - end + normalize!(envs, ψ, operator, ψ) - normalize!(envs, ψ, operator, ψ) + # 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) - # 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) + 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, end] = 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) - ψ.AL[row + 1, end] = al - ψ.C[row + 1, end] = complex(c) - ψ.AR[row + 1, 1] = _transpose_front(ar) + ψ.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 - ψ.AR[row + 1, end] = _transpose_front( - ψ.C[row + 1, end - 1] \ _transpose_tail(al * c) - ) - ψ.AC[row + 1, 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, ψ, 1) - transfer_rightenv!(envs, ψ, operator, ψ, 0) + transfer_leftenv!(envs, ψ, operator, ψ, 1) + transfer_rightenv!(envs, ψ, operator, ψ, 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 + # 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) - 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 + 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 + + transfer_leftenv!(envs, ψ, operator, ψ, site + 1) + transfer_rightenv!(envs, ψ, operator, ψ, site) end - alg_gauge = adapt_solver(alg.alg_gauge; iter, g_global = ϵ) - ψ = MultilineMPS(map(identity, ψ.AR); alg_gauge.tol, alg_gauge.maxiter) + normalize!(envs, ψ, operator, ψ) - recalculate!(envs, ψ, operator, ψ) - return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, bondresidual = ϵ, truncation_errors = parent(ϵ_truncs), numiter = iter) + # 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) + + 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) + + ψ.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, ψ, operator, ψ, 1) + transfer_rightenv!(envs, ψ, operator, ψ, 0) + + # update error + return 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 end From 60bd1e51a26dd96666ece3d36b26293ba8a12412 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 19:22:19 -0400 Subject: [PATCH 06/10] Restructure finite TDVP sweeps into per-site local updates Co-Authored-By: Claude Opus 5.5 --- src/algorithms/timestep/tdvp.jl | 190 +++++++++++++++++--------------- 1 file changed, 101 insertions(+), 89 deletions(-) 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) From 56f8d9ec1bdedd5c06dea6f0d70d4d35cdb30bc6 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 19:58:00 -0400 Subject: [PATCH 07/10] Drive time_evolve with a time-step iterator Co-Authored-By: Claude Opus 5.5 --- src/algorithms/groundstate/dmrg.jl | 2 +- src/algorithms/groundstate/vumps.jl | 2 +- src/algorithms/statmech/vomps.jl | 2 +- src/algorithms/timestep/time_evolve.jl | 88 ++++++++++++++++++-------- src/states/ortho.jl | 4 +- 5 files changed, 65 insertions(+), 33 deletions(-) diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 2b3e8c9c7..453fbe02c 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -340,7 +340,7 @@ function sweep!(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state, direction, ite ) end -function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state = it.state) +function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::DMRGState = it.state) iter = state.iter + 1 timeroutput = state.timeroutput state = @timeit timeroutput "sweep" begin diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index 9e7b144ef..487e4b134 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -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 = it.state) 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/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl index 55c006cef..9d17f15ef 100644 --- a/src/algorithms/statmech/vomps.jl +++ b/src/algorithms/statmech/vomps.jl @@ -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/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 9d9e87f40..489a3a3e0 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(string(nameof(typeof(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..6cbc8d930 100644 --- a/src/states/ortho.jl +++ b/src/states/ortho.jl @@ -250,7 +250,7 @@ function uniform_leftorth!( end end -function Base.iterate(it::IterativeSolver{LeftCanonical}, state = it.state) +function Base.iterate(it::IterativeSolver{LeftCanonical}, state::NamedTuple = it.state) timeroutput = state.timeroutput C₀ = @timeit timeroutput "gauge_eigsolve" gauge_eigsolve_step!(it, state) C₁ = @timeit timeroutput "gauge_orth" gauge_orth_step!(it, state) @@ -318,7 +318,7 @@ function uniform_rightorth!( end end -function Base.iterate(it::IterativeSolver{RightCanonical}, state = it.state) +function Base.iterate(it::IterativeSolver{RightCanonical}, state::NamedTuple = it.state) timeroutput = state.timeroutput C₀ = @timeit timeroutput "gauge_eigsolve" gauge_eigsolve_step!(it, state) C₁ = @timeit timeroutput "gauge_orth" gauge_orth_step!(it, state) From b5ca3adc861b99bf329b7ca5fc97223f6ca9fffb Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 9 Oct 2026 13:52:16 -0400 Subject: [PATCH 08/10] Address review: uniform half-sweep iterators and generic IterLog All sweep-based iterators now share DMRG's shape: `Base.iterate` runs `sweep!(it, state, Val(:right), iter)` and `sweep!(it, state, Val(:left), iter)`, with finite sweeps taking their sites from `_sweep_ranges`. The IDMRG convergence measure is shared as `_center_change`. `IterLog(::Algorithm)` names logs after the algorithm, and the redundant `= it.state` defaults on `Base.iterate` are dropped. Co-Authored-By: Claude Opus 5.5 --- src/algorithms/approximate/fvomps.jl | 38 +++++++++------ src/algorithms/approximate/idmrg.jl | 65 ++++++++++++++------------ src/algorithms/approximate/vomps.jl | 2 +- src/algorithms/groundstate/dmrg.jl | 6 +-- src/algorithms/groundstate/idmrg.jl | 4 +- src/algorithms/groundstate/vumps.jl | 4 +- src/algorithms/propagator/corvector.jl | 47 ++++++++++++++----- src/algorithms/statmech/idmrg.jl | 59 ++++++++++------------- src/algorithms/statmech/vomps.jl | 2 +- src/algorithms/timestep/time_evolve.jl | 2 +- src/states/ortho.jl | 11 +++-- src/utility/iterativesolvers.jl | 3 ++ 12 files changed, 137 insertions(+), 106 deletions(-) diff --git a/src/algorithms/approximate/fvomps.jl b/src/algorithms/approximate/fvomps.jl index 5a17d7a22..ed17ce1a0 100644 --- a/src/algorithms/approximate/fvomps.jl +++ b/src/algorithms/approximate/fvomps.jl @@ -14,7 +14,7 @@ function approximate!( envs = environments(ψ, _environment_args(Oϕ)...) ) allocator = default_allocator(ψ, SerialScheduler()) - log = IterLog(string(nameof(typeof(alg)))) + log = IterLog(alg) it = IterativeSolver(alg, ApproximateState(ψ, Oϕ, envs, 0, 2 * alg.tol, allocator)) with_verbosity(; alg.verbosity) do @@ -38,19 +38,25 @@ end function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::ApproximateState) iter = state.iter + 1 - ϵ = approximate_sweep!(state.mps, state.operator, it.alg, state.envs, state.allocator) + 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.allocator) - return (ψ, envs, ϵ), it.state + it.state = ApproximateState(ψ, state.operator, envs, iter, state.ϵ, state.allocator) + return (ψ, envs, state.ϵ), it.state end -function approximate_sweep!(ψ, Oϕ, alg::DMRG2, envs, allocator) - ϵ = 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)) +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)) AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1]) ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2)) @@ -58,17 +64,19 @@ function approximate_sweep!(ψ, Oϕ, alg::DMRG2, envs, allocator) ψ.AC[pos] = (al, complex(c)) ψ.AC[pos + 1] = (complex(c), _transpose_front(ar)) end - return ϵ + return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator) end -function approximate_sweep!(ψ, Oϕ, alg::DMRG, envs, allocator) - ϵ = 0.0 - for pos in [1:(length(ψ) - 1); length(ψ):-1:2] - AC′ = AC_projection(pos, ψ, Oϕ, envs; alg.backend, allocator) +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 ϵ + return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator) end diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl index 45fa513ab..fec1aa5cb 100644 --- a/src/algorithms/approximate/idmrg.jl +++ b/src/algorithms/approximate/idmrg.jl @@ -16,7 +16,7 @@ function approximate!( ) allocator = default_allocator(ψ, SerialScheduler()) alg isa IDMRG2 && width(ψ) < 2 && throw(ArgumentError("unit cell should be >= 2")) - log = IterLog(string(nameof(typeof(alg)))) + log = IterLog(alg) ϵ_truncs = alg isa IDMRG2 ? PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ))) : nothing state = IDMRGApproximateState(ψ, toapprox, envs, 0, 2 * alg.tol, ϵ_truncs, allocator) @@ -53,22 +53,31 @@ function approximate!( end function Base.iterate(it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state::IDMRGApproximateState) - ϵ = approximate_sweep!( - state.mps, state.operator, it.alg, state.envs, state.allocator, state.truncation_errors - ) + iter = state.iter + 1 + C_current = state.mps.C[:, 0] + state = sweep!(it, state, Val(:right), iter) + state = sweep!(it, state, Val(:left), iter) + ϵ = _center_change(it.alg, C_current, state.mps.C[:, 0]) it.state = IDMRGApproximateState( - state.mps, state.operator, state.envs, state.iter + 1, ϵ, - state.truncation_errors, state.allocator, + state.mps, state.operator, state.envs, iter, ϵ, state.truncation_errors, state.allocator, ) return (it.state.mps, it.state.envs, it.state.ϵ), it.state end -function approximate_sweep!( - ψ::MultilineMPS, toapprox, alg::IDMRG, envs, allocator, ::Nothing - ) - C_current = ψ.C[:, 0] +# change of the center bond tensors over a sweep; for IDMRG2 the bond dimension may have +# changed, so both are compared in their common subspace +_center_change(::IDMRG, C_old, C_new) = norm(C_old - C_new) +function _center_change(::IDMRG2, C_old, C_new) + return sum(zip(C_old, C_new)) 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 +end - # left to right sweep +function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGApproximateState, ::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( @@ -80,8 +89,11 @@ function approximate_sweep!( end transfer_leftenv!(envs, ψ, toapprox, col + 1) end + return state +end - # right to left sweep +function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGApproximateState, ::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( @@ -95,16 +107,12 @@ function approximate_sweep!( transfer_rightenv!(envs, ψ, toapprox, col - 1) end normalize!(envs, ψ, toapprox) - - return norm(C_current - ψ.C[:, 0]) + return state end -function approximate_sweep!( - ψ::MultilineMPS, toapprox, alg::IDMRG2, envs, allocator, ϵ_truncs - ) - C_current = ψ.C[:, 0] - - # sweep from left to right +function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGApproximateState, ::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( @@ -142,15 +150,17 @@ function approximate_sweep!( ψ.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 - # update environments transfer_leftenv!(envs, ψ, toapprox, 1) transfer_rightenv!(envs, ψ, toapprox, 0) normalize!(envs, ψ, toapprox) + return state +end - # sweep from right to left +function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGApproximateState, ::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( @@ -191,12 +201,5 @@ function approximate_sweep!( transfer_rightenv!(envs, ψ, toapprox, 0) normalize!(envs, ψ, toapprox) - - # update error - return 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 + 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 453fbe02c..ad0ade0a1 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -320,7 +320,7 @@ function DMRGState(ψ, H, alg::Union{DMRG, DMRG2}, envs, allocator, timeroutput) ) end -function sweep!(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state, direction, iter) +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.ϵ @@ -340,7 +340,7 @@ function sweep!(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state, direction, ite ) end -function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::DMRGState = it.state) +function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::DMRGState) iter = state.iter + 1 timeroutput = state.timeroutput state = @timeit timeroutput "sweep" begin @@ -364,7 +364,7 @@ sweep_converged(alg, state::DMRGState) = function find_groundstate_sweep!( ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs, allocator, timeroutput ) - log = IterLog(string(nameof(typeof(alg)))) + log = IterLog(alg) it = IterativeSolver(alg, DMRGState(ψ, H, alg, envs, allocator, timeroutput)) with_verbosity(; alg.verbosity) do diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index a0c314148..0fad1c656 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -99,7 +99,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,7 +147,7 @@ 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{<:Any, <: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 diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index 487e4b134..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::VUMPSState = 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 c73518f96..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, ψ) @@ -82,17 +84,26 @@ struct DDMRGState{S, T, Z, O, E, A} 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) - ϵ = ddmrg_sweep!(state, it.alg) + 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, state.iter + 1, ϵ, + state.mps, state.target, state.z, state.operator, state.envs, iter, state.ϵ, state.allocator, ) - return (state.mps, ϵ), it.state + return (state.mps, state.ϵ), it.state end function _propagator_sweeps!(alg::DynamicalDMRG, state::DDMRGState) - log = IterLog("DDMRG") + log = IterLog(alg) it = IterativeSolver(alg, state) with_verbosity(; alg.verbosity) do @@ -126,13 +137,17 @@ function propagator( return dot(A, init), init end -function ddmrg_sweep!(state::DDMRGState, alg::DynamicalDMRG{NaiveInvert}) +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 - ϵ = 0.0 + fwd, bwd = _sweep_ranges(alg, A) + ϵ = state.ϵ - for i in [1:(length(A) - 1); length(A):-1:2] + 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) @@ -146,7 +161,9 @@ function ddmrg_sweep!(state::DDMRGState, alg::DynamicalDMRG{NaiveInvert}) @warn "propagator ($i) failed to converge: normres = $(convhist.normres)" end - return ϵ + return DDMRGState( + init, A, z, operator, state.envs, state.iter, ϵ, allocator, + ) end """ @@ -204,16 +221,20 @@ function propagator( return v, init end -function ddmrg_sweep!(state::DDMRGState, alg::DynamicalDMRG{Jeckelmann}) +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) - ϵ = 0.0 + fwd, bwd = _sweep_ranges(alg, A) + ϵ = state.ϵ - for i in [1:(length(A) - 1); length(A):-1:2] + 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) @@ -227,7 +248,9 @@ function ddmrg_sweep!(state::DDMRGState, alg::DynamicalDMRG{Jeckelmann}) @warn "propagator ($i) failed to converge: normres $(convhist.normres)" end - return ϵ + return DDMRGState( + init, A, z, state.operator, state.envs, state.iter, ϵ, allocator, + ) end function squaredenvs( diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index e8161aaf2..076d9c935 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -16,7 +16,7 @@ function leading_boundary( ) allocator = default_allocator(ψ, SerialScheduler()) alg isa IDMRG2 && width(ψ) < 2 && throw(ArgumentError("unit cell should be >= 2")) - log = IterLog(string(nameof(typeof(alg)))) + log = IterLog(alg) ϵ_truncs = alg isa IDMRG2 ? PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ))) : nothing state = IDMRGBoundaryState(ψ, operator, envs, 0, 2 * alg.tol, ϵ_truncs, allocator) @@ -54,23 +54,19 @@ _boundary_objective(::IDMRG2, ψ, operator, envs) = nothing function Base.iterate(it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state::IDMRGBoundaryState) iter = state.iter + 1 - alg_eigsolve = adapt_solver(it.alg_eigsolve; iter, g_global = state.ϵ) - ϵ = leading_boundary_sweep!( - state.mps, state.operator, it.alg, state.envs, alg_eigsolve, state.allocator, - state.truncation_errors, - ) + C_current = state.mps.C[:, 0] + state = sweep!(it, state, Val(:right), iter) + state = sweep!(it, state, Val(:left), iter) + ϵ = _center_change(it.alg, C_current, state.mps.C[:, 0]) it.state = IDMRGBoundaryState( state.mps, state.operator, state.envs, iter, ϵ, state.truncation_errors, state.allocator, ) return (it.state.mps, it.state.envs, it.state.ϵ), it.state end -function leading_boundary_sweep!( - ψ::InfiniteMultilineMPS, operator, alg::IDMRG, envs, alg_eigsolve, allocator, ::Nothing - ) - C_current = ψ.C[:, 0] - - # left to right sweep +function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGBoundaryState, ::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) @@ -83,8 +79,12 @@ function leading_boundary_sweep!( transfer_leftenv!(envs, ψ, operator, ψ, col + 1) end + return state +end - # right to left sweep +function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGBoundaryState, ::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) @@ -96,18 +96,14 @@ function leading_boundary_sweep!( transfer_rightenv!(envs, ψ, operator, ψ, col - 1) end - normalize!(envs, ψ, operator, ψ) - - return norm(C_current - ψ.C[:, 0]) + return state end -function leading_boundary_sweep!( - ψ::InfiniteMultilineMPS, operator, alg::IDMRG2, envs, alg_eigsolve, allocator, ϵ_truncs - ) - C_current = ψ.C[:, 0] - - # sweep from left to right +function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGBoundaryState, ::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) @@ -150,13 +146,15 @@ function leading_boundary_sweep!( ψ.AL[row + 1, 1] = ψ.AC[row + 1, 1] / ψ.C[row + 1, 1] end - # TODO: decide if we should compare at the half-sweep level? - # C_current = ψ.C[:, site] - transfer_leftenv!(envs, ψ, operator, ψ, 1) transfer_rightenv!(envs, ψ, operator, ψ, 0) + return state +end - # sweep from right to left +function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGBoundaryState, ::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) @@ -201,12 +199,5 @@ function leading_boundary_sweep!( transfer_leftenv!(envs, ψ, operator, ψ, 1) transfer_rightenv!(envs, ψ, operator, ψ, 0) - - # update error - return 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 + return state end diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl index 9d17f15ef..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 = ϵ) diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 489a3a3e0..61c64bd69 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -80,7 +80,7 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time end function _time_evolve(alg, state::TimeEvolveState; verbosity::Int) - log = IterLog(string(nameof(typeof(alg)))) + 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 diff --git a/src/states/ortho.jl b/src/states/ortho.jl index 6cbc8d930..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::NamedTuple = 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::NamedTuple = 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)))) From 146d8597c70df1c1720b203d6d532ab67057e3ca Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 9 Oct 2026 14:32:00 -0400 Subject: [PATCH 09/10] Add changelog entry for iterator-based sweep solvers Co-Authored-By: Claude Opus 5.5 --- docs/src/changelog.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 7efbec564..cdd0fdce1 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -27,6 +27,12 @@ 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)) + ### Deprecated ### Removed From e3bbaea131cbbb4634fb0857a987f9df84fa3183 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 9 Oct 2026 15:29:59 -0400 Subject: [PATCH 10/10] Share IDMRG state and bond-change measure across solvers `bond_change` replaces the three copies of the center-bond convergence measure (groundstate, multiline approximation and leading boundary IDMRG/IDMRG2). It compares bond tensors of different dimension on their common subspace through block views instead of isometry products, and sums the changes over the rows of a multiline MPS, as IDMRG2 already did. The multiline approximation and leading boundary searches now use the groundstate `IDMRGState`, dispatching on its MPS and operator types. Co-Authored-By: Claude Opus 5.5 --- docs/src/changelog.md | 4 +++ src/algorithms/approximate/idmrg.jl | 43 +++++++---------------------- src/algorithms/groundstate/idmrg.jl | 21 +++++--------- src/algorithms/statmech/idmrg.jl | 31 +++++++-------------- src/utility/utility.jl | 24 ++++++++++++++++ 5 files changed, 55 insertions(+), 68 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index cdd0fdce1..09b1d49d1 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -32,6 +32,10 @@ When releasing a new version, move the "Unreleased" changes to a new version sec 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 diff --git a/src/algorithms/approximate/idmrg.jl b/src/algorithms/approximate/idmrg.jl index fec1aa5cb..9fe55fbf7 100644 --- a/src/algorithms/approximate/idmrg.jl +++ b/src/algorithms/approximate/idmrg.jl @@ -1,15 +1,3 @@ -# Internal state of the multiline IDMRG/IDMRG2 approximation, where `ϵ` is the change of the -# center bond tensor over the last sweep. IDMRG2 reuses one `truncation_errors` matrix across sweeps. -struct IDMRGApproximateState{S, O, E, V, A} - mps::S - operator::O - envs::E - iter::Int - ϵ::Float64 - truncation_errors::V - allocator::A -end - function approximate!( ψ::MultilineMPS, toapprox::Tuple{<:MultilineMPO, <:MultilineMPS}, alg::Union{IDMRG, IDMRG2}, envs = environments(ψ, toapprox...) @@ -19,7 +7,7 @@ function approximate!( log = IterLog(alg) ϵ_truncs = alg isa IDMRG2 ? PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ))) : nothing - state = IDMRGApproximateState(ψ, toapprox, envs, 0, 2 * alg.tol, ϵ_truncs, allocator) + state = IDMRGState(ψ, toapprox, envs, 0, 2 * alg.tol, ϵ_truncs, nothing, NoTimerOutput(), allocator) it = IterativeSolver(alg, state) with_verbosity(; alg.verbosity) do @@ -52,31 +40,20 @@ function approximate!( return ψ, envs, info end -function Base.iterate(it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state::IDMRGApproximateState) +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) - ϵ = _center_change(it.alg, C_current, state.mps.C[:, 0]) - it.state = IDMRGApproximateState( - state.mps, state.operator, state.envs, iter, ϵ, state.truncation_errors, state.allocator, + ϵ = 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, ) return (it.state.mps, it.state.envs, it.state.ϵ), it.state end -# change of the center bond tensors over a sweep; for IDMRG2 the bond dimension may have -# changed, so both are compared in their common subspace -_center_change(::IDMRG, C_old, C_new) = norm(C_old - C_new) -function _center_change(::IDMRG2, C_old, C_new) - return sum(zip(C_old, C_new)) 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 -end - -function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGApproximateState, ::Val{:right}, iter) +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) @@ -92,7 +69,7 @@ function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGApproximateState, ::Va return state end -function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGApproximateState, ::Val{:left}, iter) +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) @@ -110,7 +87,7 @@ function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGApproximateState, ::Va return state end -function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGApproximateState, ::Val{:right}, iter) +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) @@ -158,7 +135,7 @@ function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGApproximateState, ::V return state end -function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGApproximateState, ::Val{:left}, iter) +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)) diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 0fad1c656..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 @@ -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::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/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 076d9c935..fe5e0073f 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -1,15 +1,3 @@ -# Internal state of the IDMRG/IDMRG2 leading boundary search, where `ϵ` is the change of the -# center bond tensor over the last sweep. IDMRG2 reuses one `truncation_errors` matrix across sweeps. -struct IDMRGBoundaryState{S, O, E, V, A} - mps::S - operator::O - envs::E - iter::Int - ϵ::Float64 - truncation_errors::V - allocator::A -end - function leading_boundary( ψ::InfiniteMultilineMPS, operator, alg::Union{IDMRG, IDMRG2}, envs = environments(ψ, operator, ψ) @@ -19,7 +7,7 @@ function leading_boundary( log = IterLog(alg) ϵ_truncs = alg isa IDMRG2 ? PeriodicMatrix(zeros(real(scalartype(ψ)), length(ψ), width(ψ))) : nothing - state = IDMRGBoundaryState(ψ, operator, envs, 0, 2 * alg.tol, ϵ_truncs, allocator) + state = IDMRGState(ψ, operator, envs, 0, 2 * alg.tol, ϵ_truncs, nothing, NoTimerOutput(), allocator) it = IterativeSolver(alg, state) with_verbosity(; alg.verbosity) do @@ -52,19 +40,20 @@ end _boundary_objective(::IDMRG, ψ, operator, envs) = leading_eigenvalue(ψ, operator, envs) _boundary_objective(::IDMRG2, ψ, operator, envs) = nothing -function Base.iterate(it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state::IDMRGBoundaryState) +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) - ϵ = _center_change(it.alg, C_current, state.mps.C[:, 0]) - it.state = IDMRGBoundaryState( - state.mps, state.operator, state.envs, iter, ϵ, state.truncation_errors, state.allocator, + ϵ = 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, ) return (it.state.mps, it.state.envs, it.state.ϵ), it.state end -function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGBoundaryState, ::Val{:right}, iter) +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(ψ) @@ -82,7 +71,7 @@ function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGBoundaryState, ::Val{: return state end -function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGBoundaryState, ::Val{:left}, iter) +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 @@ -100,7 +89,7 @@ function sweep!(it::IterativeSolver{<:IDMRG}, state::IDMRGBoundaryState, ::Val{: return state end -function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGBoundaryState, ::Val{:right}, iter) +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 @@ -151,7 +140,7 @@ function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGBoundaryState, ::Val{ return state end -function sweep!(it::IterativeSolver{<:IDMRG2}, state::IDMRGBoundaryState, ::Val{:left}, iter) +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 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})