From b031218253d5ef9f170f95386f0c9e6d4a5e8e5f Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 07:45:08 -0400 Subject: [PATCH 01/15] 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 55390b25553757e63abf8c1128f86a7347f96b8c Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 18:04:14 -0400 Subject: [PATCH 02/15] 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 a7212eb8376b5c3a614db59081d950c80fc263f0 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 19:14:38 -0400 Subject: [PATCH 03/15] 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 c4f0ab1f84bf133b661d8261baf8a948ca0f691e Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 19:16:39 -0400 Subject: [PATCH 04/15] 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 614baf4251c3f8dfdd299bcc3181b32b3c220d02 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 19:16:40 -0400 Subject: [PATCH 05/15] 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 2f9e576abcb75c2fbc1c8b1348d3ecd2de2d8fbe Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 19:22:19 -0400 Subject: [PATCH 06/15] 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 4e187e2f6f6edef0843df89e1d36cb27b1af9bb5 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 19:58:00 -0400 Subject: [PATCH 07/15] 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 e3ee57ae75866004f7dfbe888c9e11b17fe370f7 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 9 Oct 2026 13:52:16 -0400 Subject: [PATCH 08/15] 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 6c6cda7b06fd27f5322a6593c8777d795f1143ef Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 9 Oct 2026 14:32:00 -0400 Subject: [PATCH 09/15] Add changelog entry for iterator-based sweep solvers Co-Authored-By: Claude Opus 5.5 --- docs/src/changelog.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 64bb570b0..93f2679f9 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -33,6 +33,11 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - `left_canonicalize!`/`right_canonicalize!(::FiniteMPOHamiltonian)` now sweep over the whole chain instead of a single site, and return the gauge transformations alongside the operator, matching the infinite versions. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539)) +- 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 From 903404a0a0d60723eaa92b1b62b652ae0568c4ee Mon Sep 17 00:00:00 2001 From: lkdvos Date: Fri, 9 Oct 2026 15:29:59 -0400 Subject: [PATCH 10/15] 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 93f2679f9..50c973b51 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -38,6 +38,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}) From d0fc5d95951ccc720323a2c19b013ef91a353c35 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 12 Oct 2026 11:35:01 -0400 Subject: [PATCH 11/15] Share one per-step sweep across finite TDVP, TDVP2 and BUG MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A `TimeStepState` carries a single time step `t → t + dt` through its two half-sweeps, and one `sweep!(alg, state, direction)` drives the per-site `local_update!(site, direction, ψ, alg, state)` of each integrator, taking the sites from `_sweep_ranges` and recording truncation errors per bond. Co-Authored-By: Claude Opus 5.5 --- docs/src/changelog.md | 4 +- src/algorithms/timestep/bug.jl | 44 +++------- src/algorithms/timestep/tdvp.jl | 139 ++++++++++++++++---------------- 3 files changed, 86 insertions(+), 101 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 50c973b51..7d62ecdc3 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -36,8 +36,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - 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)) + sharing one structure across solvers. Finite `TDVP`, `TDVP2` and `BUG` steps share a single sweep + over 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 diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl index 24df2c962..a1a9e20fe 100644 --- a/src/algorithms/timestep/bug.jl +++ b/src/algorithms/timestep/bug.jl @@ -133,10 +133,13 @@ function _evolve_center(site, ψ, H, alg::BUG, envs, t, h, allocator; imaginary_ return integrate(Heff, ψ.AC[site], t, h, alg.integrator; imaginary_evolution) end -function local_update!( - site, direction::Val, ψ, H, alg::BUG, envs, t, h, allocator; - imaginary_evolution, normalize, timeroutput - ) +# a half-sweep evolves over `h = dt / 2`, starting at `t` (left→right) or `t + h` (right→left) +function local_update!(site, direction::Val, ψ, alg::BUG, state::TimeStepState) + (; operator, envs, allocator, timeroutput, imaginary_evolution, normalize) = state + H = operator + h = state.dt / 2 + t = direction === Val(:right) ? state.t : state.t + h + # at the far end of the sweep there is no bond ahead to cut: evolve and finalize if site == _sweep_end(ψ, direction) AC = @timeit timeroutput "AC_integrate" _evolve_center( @@ -171,35 +174,14 @@ function timestep!( imaginary_evolution::Bool = false, normalize::Bool = false, timeroutput = NoTimerOutput() ) - L = length(ψ) - h = dt / 2 - - # the sweep is serial, so a single allocator serves all local updates - allocator = default_allocator(ψ, SerialScheduler()) - - ϵ_truncs = zeros(real(scalartype(ψ)), L - 1) - - # left→right half-sweep (root = last site): `t → t + dt / 2` - @timeit timeroutput "half-sweep" for site in 1:L - ψ, ϵ = local_update!( - site, Val(:right), ψ, H, alg, envs, t, h, allocator; - imaginary_evolution, normalize, timeroutput - ) - site < L && (ϵ_truncs[site] = ϵ) - end - - # right→left half-sweep (root = first site): `t + dt / 2 → t + dt` - @timeit timeroutput "half-sweep" for site in L:-1:1 - ψ, ϵ = local_update!( - site, Val(:left), ψ, H, alg, envs, t + h, h, allocator; - imaginary_evolution, normalize, timeroutput - ) - site > 1 && (ϵ_truncs[site - 1] = ϵ) - end - - return ψ, envs, AlgorithmInfo(; truncation_errors = ϵ_truncs) + return _timestep_finite!( + ψ, H, t, dt, alg, envs; imaginary_evolution, normalize, timeroutput + ) end +# like single-site TDVP, BUG visits every site and only evolves the center at the far end +_sweep_ranges(::BUG, ψ) = (1:length(ψ), length(ψ):-1:1) + # copying version function timestep( ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index c87d74684..bac5c9888 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -143,22 +143,78 @@ function timestep!( envs::AbstractMPSEnvironments = environments(ψ, H, ψ); imaginary_evolution::Bool = false, normalize::Bool = false ) + return _timestep_finite!(ψ, H, t, dt, alg, envs; imaginary_evolution, normalize) +end + +# Internal state of a single finite time step `t → t + dt`, threaded through its two half-sweeps. +# `truncation_errors[i]` holds the truncation error of the last cut made at bond `i`. +struct TimeStepState{S, O, E, T, D, V, A, TO} + mps::S + operator::O + envs::E + t::T + dt::D + truncation_errors::V + allocator::A + timeroutput::TO + imaginary_evolution::Bool + normalize::Bool +end + +function TimeStepState( + ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput = NoTimerOutput() + ) # the sweep is serial, so a single allocator serves all local updates allocator = default_allocator(ψ, SerialScheduler()) - return _timestep_finite!( - ψ, H, t, dt, alg, envs, allocator; imaginary_evolution, normalize + truncation_errors = zeros(real(scalartype(ψ)), length(ψ) - 1) + return TimeStepState( + ψ, H, envs, t, dt, truncation_errors, allocator, timeroutput, imaginary_evolution, normalize + ) +end + +# Single-site integrators visit every site, finishing a half-sweep by evolving only the center +# tensor at its far end. +_sweep_ranges(::TDVP, ψ) = (1:length(ψ), length(ψ):-1:1) + +# bond cut by the update at `site`, or `nothing` at the far end of a single-site half-sweep +_swept_bond(alg, site, ψ, ::Val{:right}) = site < length(ψ) ? site : nothing +_swept_bond(alg, site, ψ, ::Val{:left}) = site > 1 ? site - 1 : nothing + +function sweep!(alg, state::TimeStepState, direction) + fwd, bwd = _sweep_ranges(alg, state.mps) + sites = direction === Val(:right) ? fwd : bwd + ψ = state.mps + @timeit state.timeroutput "half-sweep" for site in sites + ψ, ϵ = local_update!(site, direction, ψ, alg, state) + bond = _swept_bond(alg, site, ψ, direction) + isnothing(bond) || (state.truncation_errors[bond] = ϵ) + end + return TimeStepState( + ψ, state.operator, state.envs, state.t, state.dt, state.truncation_errors, + state.allocator, state.timeroutput, state.imaginary_evolution, state.normalize, ) end +# a step is a left→right half-sweep `t → t + dt / 2` followed by a right→left half-sweep +# `t + dt / 2 → t + dt` +function _timestep_finite!( + ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg, envs; + imaginary_evolution::Bool, normalize::Bool, timeroutput = NoTimerOutput() + ) + state = TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput) + state = sweep!(alg, state, Val(:right)) + state = sweep!(alg, state, Val(:left)) + return state.mps, state.envs, AlgorithmInfo(; truncation_errors = state.truncation_errors) +end + # 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 - ) +function local_update!(site, direction::Val, ψ, alg::TDVP, state::TimeStepState) + (; operator, envs, t, dt, allocator, imaginary_evolution, normalize) = state + H = operator t_AC, t_C = _half_sweep_times(direction, t, dt) # at the far end of the sweep there is no bond ahead: only evolve the center tensor @@ -194,34 +250,6 @@ function local_update!( 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 - """ $(TYPEDEF) @@ -268,17 +296,17 @@ function timestep!( envs::AbstractMPSEnvironments = environments(ψ, H, ψ); imaginary_evolution::Bool = false, normalize::Bool = false ) - # the sweep is serial, so a single allocator serves all local updates - allocator = default_allocator(ψ, SerialScheduler()) - return _timestep2_finite!( - ψ, H, t, dt, alg, envs, allocator; imaginary_evolution, normalize - ) + return _timestep_finite!(ψ, H, t, dt, alg, envs; imaginary_evolution, normalize) end -function local_update!( - pos, direction::Val, ψ, H, alg::TDVP2, envs, t, dt, allocator; - imaginary_evolution, normalize - ) +# two-site TDVP visits every bond, and each update cuts the bond it acts on +_sweep_ranges(::TDVP2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 1):-1:1) +_swept_bond(::TDVP2, pos, ψ, ::Val{:right}) = pos +_swept_bond(::TDVP2, pos, ψ, ::Val{:left}) = pos + +function local_update!(pos, direction::Val, ψ, alg::TDVP2, state::TimeStepState) + (; operator, envs, t, dt, allocator, imaginary_evolution, normalize) = state + H = operator t_AC2, t_AC = _half_sweep_times(direction, t, dt) # 1. evolve the two-site center tensor at `(pos, pos + 1)` forward @@ -305,31 +333,6 @@ function local_update!( return ψ, ϵ end -function _timestep2_finite!( - ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, envs, allocator; - imaginary_evolution::Bool, normalize::Bool - ) - ϵ_truncs = zeros(real(scalartype(ψ)), length(ψ) - 1) - - # 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 - - # 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) -end - # copying version function timestep( ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, From 03f6ede3aad0ba1ae5cf8bf07ef6ebc0a8386761 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 12 Oct 2026 13:11:01 -0400 Subject: [PATCH 12/15] Reuse the finite time-step state across time_evolve steps `time_evolve` with finite TDVP, TDVP2 or BUG now keeps its `TimeStepState`, and with it the sweep allocator, from one step to the next. The copying variant copies (or promotes to complex) the initial state once, through the same `_timestep_copy` used by the copying `timestep` methods, and takes every later step in place. Other algorithms still go through `timestep`/`timestep!`. Co-Authored-By: Claude Opus 5.5 --- docs/src/changelog.md | 5 +++ src/algorithms/timestep/bug.jl | 15 +++----- src/algorithms/timestep/tdvp.jl | 41 ++++++++++++++-------- src/algorithms/timestep/time_evolve.jl | 47 ++++++++++++++++++++++---- 4 files changed, 76 insertions(+), 32 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 7d62ecdc3..de07a4b81 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -42,6 +42,11 @@ When releasing a new version, move the "Unreleased" changes to a new version sec 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)) +- `time_evolve` with finite `TDVP`, `TDVP2` or `BUG` copies (or promotes to complex) the initial + state once and then takes every step in place, reusing one sweep state and its allocator across + steps, rather than copying the state at every step. A `finalize` callback that keeps the state + of an earlier step therefore sees it updated by later steps, as with `time_evolve!`. + ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537)) ### Deprecated diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl index a1a9e20fe..0458566e1 100644 --- a/src/algorithms/timestep/bug.jl +++ b/src/algorithms/timestep/bug.jl @@ -182,22 +182,15 @@ end # like single-site TDVP, BUG visits every site and only evolves the center at the far end _sweep_ranges(::BUG, ψ) = (1:length(ψ), length(ψ):-1:1) +_promotes_to_complex(ψ, H, ::BUG; imaginary_evolution) = + scalartype(ψ) <: Real && !(imaginary_evolution && scalartype(H) <: Real) + # copying version function timestep( ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, alg::BUG, envs::AbstractMPSEnvironments...; imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... ) - isreal = (scalartype(ψ) <: Real && !(imaginary_evolution && scalartype(H) <: Real)) - ψ′ = isreal ? complex(ψ) : copy(ψ) - if length(envs) != 0 && isreal - @warn "Currently cannot reuse real environments for complex evolution" - envs′ = environments(ψ′, H, ψ′) - elseif length(envs) == 1 - envs′ = only(envs) - else - @assert length(envs) == 0 "Invalid signature" - envs′ = environments(ψ′, H, ψ′) - end + ψ′, envs′ = _timestep_copy(ψ, H, envs, _promotes_to_complex(ψ, H, alg; imaginary_evolution)) return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...) end diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index bac5c9888..e17cd6be1 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -161,11 +161,12 @@ struct TimeStepState{S, O, E, T, D, V, A, TO} normalize::Bool end +# the sweep is serial, so a single allocator serves all local updates, and it can be carried over +# to later steps function TimeStepState( - ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput = NoTimerOutput() + ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput = NoTimerOutput(), + allocator = default_allocator(ψ, SerialScheduler()) ) - # the sweep is serial, so a single allocator serves all local updates - allocator = default_allocator(ψ, SerialScheduler()) truncation_errors = zeros(real(scalartype(ψ)), length(ψ) - 1) return TimeStepState( ψ, H, envs, t, dt, truncation_errors, allocator, timeroutput, imaginary_evolution, normalize @@ -202,10 +203,13 @@ function _timestep_finite!( imaginary_evolution::Bool, normalize::Bool, timeroutput = NoTimerOutput() ) state = TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput) - state = sweep!(alg, state, Val(:right)) - state = sweep!(alg, state, Val(:left)) + state = _timestep_finite!(alg, state) return state.mps, state.envs, AlgorithmInfo(; truncation_errors = state.truncation_errors) end +function _timestep_finite!(alg, state::TimeStepState) + state = sweep!(alg, state, Val(:right)) + return sweep!(alg, state, Val(:left)) +end # 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` @@ -333,15 +337,11 @@ function local_update!(pos, direction::Val, ψ, alg::TDVP2, state::TimeStepState return ψ, ϵ end -# copying version -function timestep( - ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, - alg::Union{TDVP, TDVP2}, envs::AbstractMPSEnvironments...; - imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... - ) - isreal = (scalartype(ψ) <: Real && !imaginary_evolution) - ψ′ = isreal ? complex(ψ) : copy(ψ) - if length(envs) != 0 && isreal +# Copy of `ψ` to step out of place, promoted to complex when the evolution leaves the reals, +# together with environments for that copy. +function _timestep_copy(ψ, H, envs, promote::Bool) + ψ′ = promote ? complex(ψ) : copy(ψ) + if length(envs) != 0 && promote @warn "Currently cannot reuse real environments for complex evolution" envs′ = environments(ψ′, H, ψ′) elseif length(envs) == 1 @@ -350,5 +350,18 @@ function timestep( @assert length(envs) == 0 "Invalid signature" envs′ = environments(ψ′, H, ψ′) end + return ψ′, envs′ +end + +_promotes_to_complex(ψ, H, ::Union{TDVP, TDVP2}; imaginary_evolution) = + scalartype(ψ) <: Real && !imaginary_evolution + +# copying version +function timestep( + ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, + alg::Union{TDVP, TDVP2}, envs::AbstractMPSEnvironments...; + imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... + ) + ψ′, envs′ = _timestep_copy(ψ, H, envs, _promotes_to_complex(ψ, H, alg; imaginary_evolution)) return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...) end diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 61c64bd69..6d8e72684 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -37,8 +37,9 @@ 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} +# `stepper` (`timestep` or `timestep!`) called with `kwargs`. Integrators that sweep a finite MPS +# carry their `TimeStepState` over in `step` instead, or `nothing` before the first step. +struct TimeEvolveState{S, O, E, T, F, K, P} mps::S operator::O envs::E @@ -46,6 +47,7 @@ struct TimeEvolveState{S, O, E, T, F, K} t_span::T stepper::F kwargs::K + step::P end # one iteration is a single time step followed by `finalize`, which receives the step's start time @@ -55,17 +57,48 @@ function Base.iterate(it::IterativeSolver{<:Any, <:TimeEvolveState}, state::Time 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, info, step = _evolve_step(it.alg, state, t, dt) ψ, 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 + ψ, state.operator, envs, iter, state.t_span, state.stepper, state.kwargs, step ) return (ψ, envs, info), it.state end +_evolve_step(alg, state::TimeEvolveState, t, dt) = _evolve_timestep(alg, state, t, dt) +function _evolve_timestep(alg, state::TimeEvolveState, t, dt) + ψ, envs, info = state.stepper( + state.mps, state.operator, t, dt, alg, state.envs; state.kwargs... + ) + return ψ, envs, info, nothing +end + +# A finite sweep keeps its `TimeStepState`, and with it its allocator, across steps. Only the +# first step of `time_evolve` copies the state (promoting it to complex if needed), after which +# every step is taken in place. +function _evolve_step(alg::Union{TDVP, TDVP2, BUG}, state::TimeEvolveState, t, dt) + state.mps isa AbstractFiniteMPS || return _evolve_timestep(alg, state, t, dt) + (; imaginary_evolution, normalize) = state.kwargs + H = state.operator + step = if isnothing(state.step) + ψ, envs = state.stepper === timestep ? + _timestep_copy( + state.mps, H, (state.envs,), + _promotes_to_complex(state.mps, H, alg; imaginary_evolution) + ) : (state.mps, state.envs) + TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize) + else + TimeStepState( + state.mps, H, state.envs, t, dt; imaginary_evolution, normalize, + state.step.allocator, + ) + end + step = _timestep_finite!(alg, step) + info = AlgorithmInfo(; truncation_errors = step.truncation_errors) + return step.mps, step.envs, info, step +end + for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time_evolve!)) @eval function $time_evolve( ψ, H, t_span::AbstractVector{<:Number}, alg, @@ -73,7 +106,7 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time verbosity::Int = 0, imaginary_evolution::Bool = false, normalize::Bool = false ) state = TimeEvolveState( - ψ, H, envs, 0, t_span, $timestep, (; imaginary_evolution, normalize) + ψ, H, envs, 0, t_span, $timestep, (; imaginary_evolution, normalize), nothing ) return _time_evolve(alg, state; verbosity) end From 2b34066612401f2243cf153283144666605b32a6 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 13 Oct 2026 10:29:41 -0400 Subject: [PATCH 13/15] Run infinite TDVP steps on the shared time-step state Split the infinite TDVP step into `localupdate_step!`, `gauge_step!` and an environment update on `TimeStepState`, mirroring VUMPS, and carry the step's allocator over between `time_evolve` steps. Co-Authored-By: Claude Opus 5.5 --- docs/src/changelog.md | 3 +- src/algorithms/timestep/tdvp.jl | 107 ++++++++++++++++--------- src/algorithms/timestep/time_evolve.jl | 25 +++++- 3 files changed, 92 insertions(+), 43 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index de07a4b81..3d8fc7fbe 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -37,7 +37,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec 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` and `BUG` steps share a single sweep - over per-site local updates. Results and logging are unchanged. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537)) + over per-site local updates, and infinite `TDVP` steps are split into local update, gauge and + environment steps on the same step state. 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 diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index e17cd6be1..0c790a056 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -72,6 +72,39 @@ function TDVP(; ) end +# Internal state of a single time step `t → t + dt`, threaded through the two half-sweeps of a +# finite MPS, or through the local update, gauge and environment steps of an infinite one. +# `truncation_errors[i]` holds the truncation error of the last cut made at bond `i`, or is +# `nothing` for an infinite MPS. +struct TimeStepState{S, O, E, T, D, V, A, TO} + mps::S + operator::O + envs::E + t::T + dt::D + truncation_errors::V + allocator::A + timeroutput::TO + imaginary_evolution::Bool + normalize::Bool +end + +# a finite sweep is serial, so by default a single allocator serves all local updates, and it +# can be carried over to later steps +function TimeStepState( + ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput = NoTimerOutput(), + allocator = default_allocator(ψ, SerialScheduler()) + ) + truncation_errors = _truncation_errors(ψ) + return TimeStepState( + ψ, H, envs, t, dt, truncation_errors, allocator, timeroutput, imaginary_evolution, normalize + ) +end + +# per-bond truncation errors of a finite sweep; an infinite step never truncates +_truncation_errors(ψ::AbstractFiniteMPS) = zeros(real(scalartype(ψ)), length(ψ) - 1) +_truncation_errors(ψ::InfiniteMPS) = nothing + function timestep( ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs::AbstractMPSEnvironments = environments(ψ, H, ψ); @@ -95,11 +128,34 @@ function _timestep_infinite( ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs, scheduler::Scheduler; leftorthflag, imaginary_evolution ) + # both sweeps together are a single unit of concurrent work, and share one allocator + allocator = default_allocator(ψ, scheduler) + state = TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize = false, allocator) + state = _timestep_infinite!(alg, state, scheduler; leftorthflag) + # infinite one-site TDVP runs at fixed bond dimension and never truncates, so it doesn't + # report `truncation_errors` (rather than report zeros that would look like a measurement) + # the gauge-fixing residual is controlled by `tolgauge`, not reported here + return state.mps, state.envs, AlgorithmInfo() +end + +# an infinite step evolves every center tensor at once, then restores the uniform gauge and +# recomputes the environments +function _timestep_infinite!(alg::TDVP, state::TimeStepState, scheduler::Scheduler; leftorthflag) + ACs, Cs = localupdate_step!(alg, state, scheduler) + ψ = gauge_step!(alg, state, ACs, Cs; leftorthflag) + recalculate!(state.envs, ψ, state.operator) + return TimeStepState( + ψ, state.operator, state.envs, state.t, state.dt, state.truncation_errors, + state.allocator, state.timeroutput, state.imaginary_evolution, state.normalize, + ) +end + +function localupdate_step!(alg::TDVP, state::TimeStepState, scheduler::Scheduler) + (; mps, operator, envs, t, dt, allocator, imaginary_evolution) = state + ψ, H = mps, operator temp_ACs = similar(ψ.AC) temp_Cs = similar(ψ.C) - # both sweeps together are a single unit of concurrent work, and share one allocator - allocator = default_allocator(ψ, scheduler) ac_sweep!() = tforeach(1:length(ψ); scheduler) do loc Hac = AC_hamiltonian(loc, ψ, H, ψ, envs; alg.backend, allocator) temp_ACs[loc] = integrate(Hac, ψ.AC[loc], t, dt, alg.integrator; imaginary_evolution) @@ -121,21 +177,19 @@ function _timestep_infinite( Threads.@spawn c_sweep!() end end + return temp_ACs, temp_Cs +end +function gauge_step!(alg::TDVP, state::TimeStepState, ACs, Cs; leftorthflag) + ψ = state.mps if leftorthflag - regauge!.(temp_ACs, temp_Cs) - ψ′ = InfiniteMPS(temp_ACs, ψ.C[end]; tol = alg.tolgauge, maxiter = alg.gaugemaxiter) + regauge!.(ACs, Cs) + return InfiniteMPS(ACs, ψ.C[end]; tol = alg.tolgauge, maxiter = alg.gaugemaxiter) else - circshift!(temp_Cs, 1) - regauge!.(temp_Cs, temp_ACs) - ψ′ = InfiniteMPS(ψ.C[0], temp_ACs; tol = alg.tolgauge, maxiter = alg.gaugemaxiter) + circshift!(Cs, 1) + regauge!.(Cs, ACs) + return InfiniteMPS(ψ.C[0], ACs; tol = alg.tolgauge, maxiter = alg.gaugemaxiter) end - - recalculate!(envs, ψ′, H) - # infinite one-site TDVP runs at fixed bond dimension and never truncates, so it doesn't - # report `truncation_errors` (rather than report zeros that would look like a measurement) - # the gauge-fixing residual is controlled by `tolgauge`, not reported here - return ψ′, envs, AlgorithmInfo() end function timestep!( @@ -146,33 +200,6 @@ function timestep!( return _timestep_finite!(ψ, H, t, dt, alg, envs; imaginary_evolution, normalize) end -# Internal state of a single finite time step `t → t + dt`, threaded through its two half-sweeps. -# `truncation_errors[i]` holds the truncation error of the last cut made at bond `i`. -struct TimeStepState{S, O, E, T, D, V, A, TO} - mps::S - operator::O - envs::E - t::T - dt::D - truncation_errors::V - allocator::A - timeroutput::TO - imaginary_evolution::Bool - normalize::Bool -end - -# the sweep is serial, so a single allocator serves all local updates, and it can be carried over -# to later steps -function TimeStepState( - ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput = NoTimerOutput(), - allocator = default_allocator(ψ, SerialScheduler()) - ) - truncation_errors = zeros(real(scalartype(ψ)), length(ψ) - 1) - return TimeStepState( - ψ, H, envs, t, dt, truncation_errors, allocator, timeroutput, imaginary_evolution, normalize - ) -end - # Single-site integrators visit every site, finishing a half-sweep by evolving only the center # tensor at its far end. _sweep_ranges(::TDVP, ψ) = (1:length(ψ), length(ψ):-1:1) diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 6d8e72684..6a819ff90 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -37,8 +37,8 @@ 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`. Integrators that sweep a finite MPS -# carry their `TimeStepState` over in `step` instead, or `nothing` before the first step. +# `stepper` (`timestep` or `timestep!`) called with `kwargs`. `TDVP`, `TDVP2` and `BUG` carry +# their `TimeStepState` over in `step` instead, or `nothing` before the first step. struct TimeEvolveState{S, O, E, T, F, K, P} mps::S operator::O @@ -99,6 +99,27 @@ function _evolve_step(alg::Union{TDVP, TDVP2, BUG}, state::TimeEvolveState, t, d return step.mps, step.envs, info, step end +# An infinite step builds a new state anyway, so only its allocator is carried over, as long as +# the scheduler still selects the same kind of allocator. +function _evolve_step(alg::TDVP, state::TimeEvolveState{<:InfiniteMPS}, t, dt) + state.stepper === timestep || return _evolve_timestep(alg, state, t, dt) + (; imaginary_evolution) = state.kwargs + ψ = state.mps + if scalartype(ψ) <: Real && (!imaginary_evolution || !isreal(dt)) + ψ = complex(ψ) + end + scheduler = Defaults.scheduler[] + allocator = default_allocator(ψ, scheduler) + if !isnothing(state.step) && state.step.allocator isa typeof(allocator) + allocator = state.step.allocator + end + step = TimeStepState( + ψ, state.operator, state.envs, t, dt; imaginary_evolution, normalize = false, allocator + ) + step = _timestep_infinite!(alg, step, scheduler; leftorthflag = true) + return step.mps, step.envs, AlgorithmInfo(), step +end + for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time_evolve!)) @eval function $time_evolve( ψ, H, t_span::AbstractVector{<:Number}, alg, From e526bb954192400abe442f919543b4a01281e1ba Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 13 Oct 2026 11:33:40 -0400 Subject: [PATCH 14/15] Add in-place infinite TDVP and fix two infinite TDVP bugs Add `timestep!(::InfiniteMPS, ..., ::TDVP)`, which restores the gauge in place, and route infinite TDVP through the shared copying `timestep` and the in-place `time_evolve` step, so that finite and infinite steps take the same path. `TimeStepState` now carries the scheduler with its allocator. Fix `InfiniteMPS(C, AR)`, which called `eltype` instead of `scalartype` and therefore failed for infinite TDVP with `leftorthflag = false`, and rebuild the environments when a real infinite state is promoted to complex, as for finite MPS, instead of writing complex data into real environments. A real state is now also promoted for imaginary-time evolution with a complex time step. Co-Authored-By: Claude Opus 5.5 --- docs/src/changelog.md | 19 +++-- src/algorithms/timestep/bug.jl | 7 +- src/algorithms/timestep/tdvp.jl | 99 ++++++++++++++++---------- src/algorithms/timestep/time_evolve.jl | 40 +++-------- src/states/infinitemps.jl | 2 +- test/states/infinitemps.jl | 6 ++ test/timeevolution/timestep.jl | 30 ++++++++ 7 files changed, 124 insertions(+), 79 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 3d8fc7fbe..81acee240 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -24,6 +24,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - `ClusterExpansion` time-evolution algorithm for `make_time_mpo`, constructing evolution MPOs for finite and infinite nearest-neighbor Hamiltonians through a nonperturbative cluster expansion that matches exact exponentials on clusters of up to `N` sites. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252)) +- In-place `timestep!` for an `InfiniteMPS` with `TDVP`, which also enables `time_evolve!` on + infinite states. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537)) - `changebonds(::InfiniteMPOHamiltonian, ::SvdCut)` compresses infinite Hamiltonians with the almost-Schmidt decomposition of [Parker et al.](https://arxiv.org/abs/1909.06341), and `left_canonicalize!`/`right_canonicalize!` bring an `InfiniteMPOHamiltonian` into canonical form. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539)) @@ -43,11 +45,11 @@ When releasing a new version, move the "Unreleased" changes to a new version sec 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)) -- `time_evolve` with finite `TDVP`, `TDVP2` or `BUG` copies (or promotes to complex) the initial - state once and then takes every step in place, reusing one sweep state and its allocator across - steps, rather than copying the state at every step. A `finalize` callback that keeps the state - of an earlier step therefore sees it updated by later steps, as with `time_evolve!`. - ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537)) +- `time_evolve` with finite `TDVP`, `TDVP2` or `BUG`, or infinite `TDVP`, copies (or promotes to + complex) the initial state once and then takes every step in place, reusing one step state and + its allocator across steps, rather than copying the state at every step. A `finalize` callback + that keeps the state of an earlier step therefore sees it updated by later steps, as with + `time_evolve!`. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537)) ### Deprecated @@ -59,6 +61,13 @@ When releasing a new version, move the "Unreleased" changes to a new version sec - `show` methods for finite MPS/MPOs no longer check unit space condition through `oneunit`, which is not defined for multifusion spaces. ([#538](https://github.com/QuantumKitHub/MPSKit.jl/pull/538)) - The canonical forms of `MPOHamiltonian`s, and therefore `SvdCut` compression, are now correct when interactions contain identity components or redundant channels. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539)) +- `InfiniteMPS(C, AR)` no longer throws a `MethodError`, which also fixes infinite `TDVP` with + `leftorthflag = false`. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537)) +- Real-time infinite `TDVP` from a real `InfiniteMPS` with a real Hamiltonian no longer fails with + an `InexactError`: the environments are now rebuilt for the complex state, as for finite MPS. + ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537)) +- `TDVP`, `TDVP2` and `BUG` promote a real state to complex for imaginary-time evolution with a + complex time step, rather than failing. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537)) ### Performance diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl index 0458566e1..5926ddff3 100644 --- a/src/algorithms/timestep/bug.jl +++ b/src/algorithms/timestep/bug.jl @@ -182,8 +182,8 @@ end # like single-site TDVP, BUG visits every site and only evolves the center at the far end _sweep_ranges(::BUG, ψ) = (1:length(ψ), length(ψ):-1:1) -_promotes_to_complex(ψ, H, ::BUG; imaginary_evolution) = - scalartype(ψ) <: Real && !(imaginary_evolution && scalartype(H) <: Real) +_promotes_to_complex(ψ, H, dt, ::BUG; imaginary_evolution) = + scalartype(ψ) <: Real && !(imaginary_evolution && isreal(dt) && scalartype(H) <: Real) # copying version function timestep( @@ -191,6 +191,7 @@ function timestep( alg::BUG, envs::AbstractMPSEnvironments...; imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... ) - ψ′, envs′ = _timestep_copy(ψ, H, envs, _promotes_to_complex(ψ, H, alg; imaginary_evolution)) + promote = _promotes_to_complex(ψ, H, timestep, alg; imaginary_evolution) + ψ′, envs′ = _timestep_copy(ψ, H, envs, promote) return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...) end diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 0c790a056..f0697b3ce 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -76,36 +76,64 @@ end # finite MPS, or through the local update, gauge and environment steps of an infinite one. # `truncation_errors[i]` holds the truncation error of the last cut made at bond `i`, or is # `nothing` for an infinite MPS. -struct TimeStepState{S, O, E, T, D, V, A, TO} +struct TimeStepState{S, O, E, T, D, V, SC, A, TO} mps::S operator::O envs::E t::T dt::D truncation_errors::V + scheduler::SC allocator::A timeroutput::TO imaginary_evolution::Bool normalize::Bool end -# a finite sweep is serial, so by default a single allocator serves all local updates, and it -# can be carried over to later steps +# a single allocator serves all local updates of a step, and it can be carried over to later +# steps together with the scheduler that selected it function TimeStepState( ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput = NoTimerOutput(), - allocator = default_allocator(ψ, SerialScheduler()) + scheduler = _step_scheduler(ψ), allocator = default_allocator(ψ, scheduler) ) truncation_errors = _truncation_errors(ψ) return TimeStepState( - ψ, H, envs, t, dt, truncation_errors, allocator, timeroutput, imaginary_evolution, normalize + ψ, H, envs, t, dt, truncation_errors, scheduler, allocator, timeroutput, + imaginary_evolution, normalize, ) end +# a finite sweep is serial, whereas an infinite step updates all sites concurrently +_step_scheduler(::AbstractFiniteMPS) = SerialScheduler() +_step_scheduler(::InfiniteMPS) = Defaults.scheduler[] + # per-bond truncation errors of a finite sweep; an infinite step never truncates _truncation_errors(ψ::AbstractFiniteMPS) = zeros(real(scalartype(ψ)), length(ψ) - 1) _truncation_errors(ψ::InfiniteMPS) = nothing -function timestep( +# infinite one-site TDVP runs at fixed bond dimension and never truncates, so it doesn't report +# `truncation_errors` (rather than report zeros that would look like a measurement) +_step_info(state::TimeStepState) = isnothing(state.truncation_errors) ? AlgorithmInfo() : + AlgorithmInfo(; truncation_errors = state.truncation_errors) + +function _with_mps(state::TimeStepState, ψ) + return TimeStepState( + ψ, state.operator, state.envs, state.t, state.dt, state.truncation_errors, + state.scheduler, state.allocator, state.timeroutput, state.imaginary_evolution, + state.normalize, + ) +end + +# a single time step on `state`, in place +_timestep!(alg, state::TimeStepState{<:AbstractFiniteMPS}) = _timestep_finite!(alg, state) +_timestep!(alg::TDVP, state::TimeStepState{<:InfiniteMPS}) = _timestep_infinite!(alg, state) + +# whether `alg` steps `ψ` in place through a `TimeStepState` +_steps_in_place(alg, ψ) = false +_steps_in_place(alg, ψ::AbstractFiniteMPS) = true +_steps_in_place(alg::TDVP, ψ::InfiniteMPS) = true + +function timestep!( ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs::AbstractMPSEnvironments = environments(ψ, H, ψ); leftorthflag = true, imaginary_evolution::Bool = false, normalize::Bool = false @@ -113,49 +141,38 @@ function timestep( # `normalize` is accepted for signature uniformity with the finite integrators, but an # `InfiniteMPS` is always normalized to norm-1-per-site by the gauge/reconstruction below # (a structural gauge requirement, not information erasure), so the flag has no effect here. - # convert state to complex if necessary - if scalartype(ψ) <: Real && (!imaginary_evolution || !isreal(dt)) - return timestep(complex(ψ), H, t, dt, alg, envs; leftorthflag, imaginary_evolution, normalize) - end - # the scheduler is read here rather than below, so that the allocator it selects is inferable - return _timestep_infinite( + return _timestep_infinite!( ψ, H, t, dt, alg, envs, Defaults.scheduler[]; leftorthflag, imaginary_evolution ) end -function _timestep_infinite( +function _timestep_infinite!( ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs, scheduler::Scheduler; leftorthflag, imaginary_evolution ) - # both sweeps together are a single unit of concurrent work, and share one allocator - allocator = default_allocator(ψ, scheduler) - state = TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize = false, allocator) - state = _timestep_infinite!(alg, state, scheduler; leftorthflag) - # infinite one-site TDVP runs at fixed bond dimension and never truncates, so it doesn't - # report `truncation_errors` (rather than report zeros that would look like a measurement) + state = TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize = false, scheduler) + state = _timestep_infinite!(alg, state; leftorthflag) # the gauge-fixing residual is controlled by `tolgauge`, not reported here - return state.mps, state.envs, AlgorithmInfo() + return state.mps, state.envs, _step_info(state) end # an infinite step evolves every center tensor at once, then restores the uniform gauge and # recomputes the environments -function _timestep_infinite!(alg::TDVP, state::TimeStepState, scheduler::Scheduler; leftorthflag) - ACs, Cs = localupdate_step!(alg, state, scheduler) +function _timestep_infinite!(alg::TDVP, state::TimeStepState; leftorthflag = true) + ACs, Cs = localupdate_step!(alg, state) ψ = gauge_step!(alg, state, ACs, Cs; leftorthflag) recalculate!(state.envs, ψ, state.operator) - return TimeStepState( - ψ, state.operator, state.envs, state.t, state.dt, state.truncation_errors, - state.allocator, state.timeroutput, state.imaginary_evolution, state.normalize, - ) + return _with_mps(state, ψ) end -function localupdate_step!(alg::TDVP, state::TimeStepState, scheduler::Scheduler) - (; mps, operator, envs, t, dt, allocator, imaginary_evolution) = state +function localupdate_step!(alg::TDVP, state::TimeStepState) + (; mps, operator, envs, t, dt, scheduler, allocator, imaginary_evolution) = state ψ, H = mps, operator temp_ACs = similar(ψ.AC) temp_Cs = similar(ψ.C) + # both sweeps together are a single unit of concurrent work, and share one allocator ac_sweep!() = tforeach(1:length(ψ); scheduler) do loc Hac = AC_hamiltonian(loc, ψ, H, ψ, envs; alg.backend, allocator) temp_ACs[loc] = integrate(Hac, ψ.AC[loc], t, dt, alg.integrator; imaginary_evolution) @@ -180,16 +197,22 @@ function localupdate_step!(alg::TDVP, state::TimeStepState, scheduler::Scheduler return temp_ACs, temp_Cs end +# the evolved center tensors are split into the new left (or right) isometries, from which the +# uniform gauge is restored in place function gauge_step!(alg::TDVP, state::TimeStepState, ACs, Cs; leftorthflag) ψ = state.mps if leftorthflag regauge!.(ACs, Cs) - return InfiniteMPS(ACs, ψ.C[end]; tol = alg.tolgauge, maxiter = alg.gaugemaxiter) + copyto!(ψ.AL, ACs) + gaugefix!(ψ, ψ.AL, ψ.C[end]; order = :R, tol = alg.tolgauge, maxiter = alg.gaugemaxiter) else circshift!(Cs, 1) regauge!.(Cs, ACs) - return InfiniteMPS(ψ.C[0], ACs; tol = alg.tolgauge, maxiter = alg.gaugemaxiter) + copyto!(ψ.AR, ACs) + gaugefix!(ψ, ψ.AR, ψ.C[0]; order = :L, tol = alg.tolgauge, maxiter = alg.gaugemaxiter) end + mul!.(ψ.AC, ψ.AL, ψ.C) + return ψ end function timestep!( @@ -217,10 +240,7 @@ function sweep!(alg, state::TimeStepState, direction) bond = _swept_bond(alg, site, ψ, direction) isnothing(bond) || (state.truncation_errors[bond] = ϵ) end - return TimeStepState( - ψ, state.operator, state.envs, state.t, state.dt, state.truncation_errors, - state.allocator, state.timeroutput, state.imaginary_evolution, state.normalize, - ) + return _with_mps(state, ψ) end # a step is a left→right half-sweep `t → t + dt / 2` followed by a right→left half-sweep @@ -231,7 +251,7 @@ function _timestep_finite!( ) state = TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput) state = _timestep_finite!(alg, state) - return state.mps, state.envs, AlgorithmInfo(; truncation_errors = state.truncation_errors) + return state.mps, state.envs, _step_info(state) end function _timestep_finite!(alg, state::TimeStepState) state = sweep!(alg, state, Val(:right)) @@ -380,15 +400,16 @@ function _timestep_copy(ψ, H, envs, promote::Bool) return ψ′, envs′ end -_promotes_to_complex(ψ, H, ::Union{TDVP, TDVP2}; imaginary_evolution) = - scalartype(ψ) <: Real && !imaginary_evolution +_promotes_to_complex(ψ, H, dt, ::Union{TDVP, TDVP2}; imaginary_evolution) = + scalartype(ψ) <: Real && !(imaginary_evolution && isreal(dt)) # copying version function timestep( - ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, + ψ::Union{AbstractFiniteMPS, InfiniteMPS}, H, time::Number, timestep::Number, alg::Union{TDVP, TDVP2}, envs::AbstractMPSEnvironments...; imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... ) - ψ′, envs′ = _timestep_copy(ψ, H, envs, _promotes_to_complex(ψ, H, alg; imaginary_evolution)) + promote = _promotes_to_complex(ψ, H, timestep, alg; imaginary_evolution) + ψ′, envs′ = _timestep_copy(ψ, H, envs, promote) return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...) end diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 6a819ff90..a149b2664 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -37,7 +37,7 @@ 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`. `TDVP`, `TDVP2` and `BUG` carry +# `stepper` (`timestep` or `timestep!`) called with `kwargs`. Integrators that step in place carry # their `TimeStepState` over in `step` instead, or `nothing` before the first step. struct TimeEvolveState{S, O, E, T, F, K, P} mps::S @@ -74,50 +74,28 @@ function _evolve_timestep(alg, state::TimeEvolveState, t, dt) return ψ, envs, info, nothing end -# A finite sweep keeps its `TimeStepState`, and with it its allocator, across steps. Only the -# first step of `time_evolve` copies the state (promoting it to complex if needed), after which -# every step is taken in place. +# An integrator that steps in place keeps its `TimeStepState`, and with it its scheduler and +# allocator, across steps. Only the first step of `time_evolve` copies the state (promoting it to +# complex if needed), after which every step is taken in place. function _evolve_step(alg::Union{TDVP, TDVP2, BUG}, state::TimeEvolveState, t, dt) - state.mps isa AbstractFiniteMPS || return _evolve_timestep(alg, state, t, dt) + _steps_in_place(alg, state.mps) || return _evolve_timestep(alg, state, t, dt) (; imaginary_evolution, normalize) = state.kwargs H = state.operator step = if isnothing(state.step) ψ, envs = state.stepper === timestep ? _timestep_copy( state.mps, H, (state.envs,), - _promotes_to_complex(state.mps, H, alg; imaginary_evolution) + _promotes_to_complex(state.mps, H, dt, alg; imaginary_evolution) ) : (state.mps, state.envs) TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize) else TimeStepState( state.mps, H, state.envs, t, dt; imaginary_evolution, normalize, - state.step.allocator, + state.step.scheduler, state.step.allocator, ) end - step = _timestep_finite!(alg, step) - info = AlgorithmInfo(; truncation_errors = step.truncation_errors) - return step.mps, step.envs, info, step -end - -# An infinite step builds a new state anyway, so only its allocator is carried over, as long as -# the scheduler still selects the same kind of allocator. -function _evolve_step(alg::TDVP, state::TimeEvolveState{<:InfiniteMPS}, t, dt) - state.stepper === timestep || return _evolve_timestep(alg, state, t, dt) - (; imaginary_evolution) = state.kwargs - ψ = state.mps - if scalartype(ψ) <: Real && (!imaginary_evolution || !isreal(dt)) - ψ = complex(ψ) - end - scheduler = Defaults.scheduler[] - allocator = default_allocator(ψ, scheduler) - if !isnothing(state.step) && state.step.allocator isa typeof(allocator) - allocator = state.step.allocator - end - step = TimeStepState( - ψ, state.operator, state.envs, t, dt; imaginary_evolution, normalize = false, allocator - ) - step = _timestep_infinite!(alg, step, scheduler; leftorthflag = true) - return step.mps, step.envs, AlgorithmInfo(), step + step = _timestep!(alg, step) + return step.mps, step.envs, _step_info(step), step end for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time_evolve!)) diff --git a/src/states/infinitemps.jl b/src/states/infinitemps.jl index 353b544f4..f91d06835 100644 --- a/src/states/infinitemps.jl +++ b/src/states/infinitemps.jl @@ -239,7 +239,7 @@ function InfiniteMPS(C₀::MPSBondTensor, AR::AbstractVector{<:GenericMPSTensor} # initialize tensor storage AC = similar.(AR) AL = similar.(AR) - T = TensorOperations.promote_contract(eltype(AR), eltype(C₀)) + T = TensorOperations.promote_contract(scalartype(AR), scalartype(C₀)) TC = TensorOperations.tensoradd_type(T, C₀, ((1,), (2,)), false) C = similar(AR, TC) ψ = InfiniteMPS{eltype(AL), TC}(AL, AR, C, AC) diff --git a/test/states/infinitemps.jl b/test/states/infinitemps.jl index 1226c395d..90f03c468 100644 --- a/test/states/infinitemps.jl +++ b/test/states/infinitemps.jl @@ -44,6 +44,12 @@ using Adapt @test TransferMatrix(ψ.AR[i], ψ.AL[i]) * r_RL(ψ, i) ≈ r_RL(ψ, i + 1) @test TransferMatrix(ψ.AR[i], ψ.AR[i]) * r_RR(ψ, i) ≈ r_RR(ψ, i + 1) end + + # rebuilding from the left or right isometries and a bond tensor gives back the same state + for ψ′ in (InfiniteMPS(ψ.AL, ψ.C[end]; tol), InfiniteMPS(ψ.C[0], ψ.AR; tol)) + @test scalartype(ψ′) == elt + @test abs(dot(ψ, ψ′)) ≈ 1 + end end @testset "InfiniteMPS copying" begin diff --git a/test/timeevolution/timestep.jl b/test/timeevolution/timestep.jl index 3baebcf91..9c42207c8 100644 --- a/test/timeevolution/timestep.jl +++ b/test/timeevolution/timestep.jl @@ -89,6 +89,31 @@ maxbond(ψ) = maximum(i -> dim(left_virtualspace(ψ, i)), 1:length(ψ)) @test E₀ ≈ E atol = 1.0e-2 end + @testset "Infinite TDVP in place" begin + ψ, = timestep(ψ₀, H, 0.0, dt, TDVP()) + ψ′ = copy(ψ₀) + ψ″, = timestep!(ψ′, H, 0.0, dt, TDVP()) + @test ψ″ === ψ′ + @test ψ″.AC == ψ.AC + # evolving in the right-orthogonal gauge gives the same state + ψr, = timestep(ψ₀, H, 0.0, dt, TDVP(); leftorthflag = false) + @test abs(dot(ψ, ψr)) ≈ 1 atol = 1.0e-8 + end + + @testset "Infinite real TDVP" begin + Hr = force_planar(transverse_field_ising(Float64; g = 0.7)) + ψr₀ = InfiniteMPS(rand, Float64, [ℙ^2], [ℙ^8]) + # real-time evolution promotes the state, and with it the environments, to complex + ψ, envs = timestep(ψr₀, Hr, 0.0, dt, TDVP()) + @test scalartype(ψ) <: Complex + ψc, = timestep(complex(ψr₀), Hr, 0.0, dt, TDVP()) + @test abs(dot(ψ, ψc)) ≈ 1 atol = 1.0e-8 + @test expectation_value(ψ, Hr, envs) ≈ expectation_value(ψc, Hr) + # imaginary-time evolution with a real step keeps it real + ψi, = timestep(ψr₀, Hr, 0.0, dt, TDVP(); imaginary_evolution = true) + @test scalartype(ψi) == Float64 + end + Hlazy = LazySum([3 * deepcopy(H), 1.55 * deepcopy(H), -0.1 * deepcopy(H)]) @testset "Infinite LazySum TDVP" begin @@ -329,5 +354,10 @@ end ψ, envs = time_evolve(ψ₀, H, t_span, TDVP(); verbosity = 2) E = expectation_value(ψ, H, envs) @test E₀ ≈ E atol = 1.0e-2 + + ψ′ = copy(ψ₀) + ψ″, = MPSKit.time_evolve!(ψ′, H, t_span, TDVP()) + @test ψ″ === ψ′ + @test ψ″.AC == ψ.AC end end From 34aef029cf76c339403dabc8850211ec9938afc6 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 13 Oct 2026 15:50:45 -0400 Subject: [PATCH 15/15] Share one time-step driver across the time-evolution integrators Finite `timestep!` and the copying `timestep` are now single methods that build the step state and call `_timestep!`, and `time_evolve` steps every integrator in place. The step state tracks the time of the state, so the local updates no longer recompute the half-sweep times. Also trims a few internal state comments, folds the in-place infinite TDVP checks into the scheduler testset, and drops the `bondresidual` changelog entry. Co-Authored-By: Claude Opus 5.5 --- docs/src/changelog.md | 4 -- src/algorithms/groundstate/dmrg.jl | 2 - src/algorithms/groundstate/idmrg.jl | 5 +- src/algorithms/propagator/corvector.jl | 3 +- src/algorithms/statmech/idmrg.jl | 1 + src/algorithms/timestep/bug.jl | 27 +------- src/algorithms/timestep/tdvp.jl | 94 ++++---------------------- src/algorithms/timestep/time_evolve.jl | 56 ++++++++++----- test/timeevolution/timestep.jl | 28 ++++---- 9 files changed, 71 insertions(+), 149 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 81acee240..1e8ea3fa9 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -41,10 +41,6 @@ When releasing a new version, move the "Unreleased" changes to a new version sec sharing one structure across solvers. Finite `TDVP`, `TDVP2` and `BUG` steps share a single sweep over per-site local updates, and infinite `TDVP` steps are split into local update, gauge and environment steps on the same step state. 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)) - `time_evolve` with finite `TDVP`, `TDVP2` or `BUG`, or infinite `TDVP`, copies (or promotes to complex) the initial state once and then takes every step in place, reusing one step state and its allocator across steps, rather than copying the state at every step. A `finalize` callback diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index ad0ade0a1..1b9575e3a 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -295,8 +295,6 @@ function find_groundstate!( return find_groundstate_sweep!(ψ, H, alg, envs, allocator, timeroutput) end -# 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 diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 214183bc4..afe939fc4 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -70,10 +70,7 @@ Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_bounda backend::B = Defaults.backend() end -# 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. +# Internal state of the IDMRG/IDMRG2 algorithms struct IDMRGState{S, O, E, V, T, TO, A} mps::S operator::O diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl index 5450e5bbe..0704743e6 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -71,8 +71,7 @@ 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 +# Internal state of the dynamical DMRG sweeps struct DDMRGState{S, T, Z, O, E, A} mps::S target::T diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index fe5e0073f..fbfcea38c 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -37,6 +37,7 @@ function leading_boundary( end # only single-site IDMRG reports the leading eigenvalue in its log +# this bypasses issues with intermediate IDMRG2 states not being correct InfiniteMPS _boundary_objective(::IDMRG, ψ, operator, envs) = leading_eigenvalue(ψ, operator, envs) _boundary_objective(::IDMRG2, ψ, operator, envs) = nothing diff --git a/src/algorithms/timestep/bug.jl b/src/algorithms/timestep/bug.jl index 5926ddff3..f9e79e199 100644 --- a/src/algorithms/timestep/bug.jl +++ b/src/algorithms/timestep/bug.jl @@ -133,12 +133,11 @@ function _evolve_center(site, ψ, H, alg::BUG, envs, t, h, allocator; imaginary_ return integrate(Heff, ψ.AC[site], t, h, alg.integrator; imaginary_evolution) end -# a half-sweep evolves over `h = dt / 2`, starting at `t` (left→right) or `t + h` (right→left) +# a half-sweep evolves over `h = dt / 2` function local_update!(site, direction::Val, ψ, alg::BUG, state::TimeStepState) - (; operator, envs, allocator, timeroutput, imaginary_evolution, normalize) = state + (; operator, envs, t, allocator, timeroutput, imaginary_evolution, normalize) = state H = operator h = state.dt / 2 - t = direction === Val(:right) ? state.t : state.t + h # at the far end of the sweep there is no bond ahead to cut: evolve and finalize if site == _sweep_end(ψ, direction) @@ -168,30 +167,8 @@ function local_update!(site, direction::Val, ψ, alg::BUG, state::TimeStepState) return ψ, ϵ end -function timestep!( - ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::BUG, - envs::AbstractMPSEnvironments = environments(ψ, H, ψ); - imaginary_evolution::Bool = false, normalize::Bool = false, - timeroutput = NoTimerOutput() - ) - return _timestep_finite!( - ψ, H, t, dt, alg, envs; imaginary_evolution, normalize, timeroutput - ) -end - # like single-site TDVP, BUG visits every site and only evolves the center at the far end _sweep_ranges(::BUG, ψ) = (1:length(ψ), length(ψ):-1:1) _promotes_to_complex(ψ, H, dt, ::BUG; imaginary_evolution) = scalartype(ψ) <: Real && !(imaginary_evolution && isreal(dt) && scalartype(H) <: Real) - -# copying version -function timestep( - ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, - alg::BUG, envs::AbstractMPSEnvironments...; - imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... - ) - promote = _promotes_to_complex(ψ, H, timestep, alg; imaginary_evolution) - ψ′, envs′ = _timestep_copy(ψ, H, envs, promote) - return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...) -end diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index f0697b3ce..1c0c126a3 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -72,10 +72,7 @@ function TDVP(; ) end -# Internal state of a single time step `t → t + dt`, threaded through the two half-sweeps of a -# finite MPS, or through the local update, gauge and environment steps of an infinite one. -# `truncation_errors[i]` holds the truncation error of the last cut made at bond `i`, or is -# `nothing` for an infinite MPS. +# Internal state of a single time step `t → t + dt` struct TimeStepState{S, O, E, T, D, V, SC, A, TO} mps::S operator::O @@ -116,23 +113,14 @@ _truncation_errors(ψ::InfiniteMPS) = nothing _step_info(state::TimeStepState) = isnothing(state.truncation_errors) ? AlgorithmInfo() : AlgorithmInfo(; truncation_errors = state.truncation_errors) -function _with_mps(state::TimeStepState, ψ) +function _with_mps(state::TimeStepState, ψ; t = state.t) return TimeStepState( - ψ, state.operator, state.envs, state.t, state.dt, state.truncation_errors, + ψ, state.operator, state.envs, t, state.dt, state.truncation_errors, state.scheduler, state.allocator, state.timeroutput, state.imaginary_evolution, state.normalize, ) end -# a single time step on `state`, in place -_timestep!(alg, state::TimeStepState{<:AbstractFiniteMPS}) = _timestep_finite!(alg, state) -_timestep!(alg::TDVP, state::TimeStepState{<:InfiniteMPS}) = _timestep_infinite!(alg, state) - -# whether `alg` steps `ψ` in place through a `TimeStepState` -_steps_in_place(alg, ψ) = false -_steps_in_place(alg, ψ::AbstractFiniteMPS) = true -_steps_in_place(alg::TDVP, ψ::InfiniteMPS) = true - function timestep!( ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs::AbstractMPSEnvironments = environments(ψ, H, ψ); @@ -152,18 +140,18 @@ function _timestep_infinite!( leftorthflag, imaginary_evolution ) state = TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize = false, scheduler) - state = _timestep_infinite!(alg, state; leftorthflag) + state = _timestep!(alg, state; leftorthflag) # the gauge-fixing residual is controlled by `tolgauge`, not reported here return state.mps, state.envs, _step_info(state) end # an infinite step evolves every center tensor at once, then restores the uniform gauge and # recomputes the environments -function _timestep_infinite!(alg::TDVP, state::TimeStepState; leftorthflag = true) +function _timestep!(alg::TDVP, state::TimeStepState{<:InfiniteMPS}; leftorthflag = true) ACs, Cs = localupdate_step!(alg, state) ψ = gauge_step!(alg, state, ACs, Cs; leftorthflag) recalculate!(state.envs, ψ, state.operator) - return _with_mps(state, ψ) + return _with_mps(state, ψ; t = state.t + state.dt) end function localupdate_step!(alg::TDVP, state::TimeStepState) @@ -215,14 +203,6 @@ function gauge_step!(alg::TDVP, state::TimeStepState, ACs, Cs; leftorthflag) return ψ end -function timestep!( - ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, - envs::AbstractMPSEnvironments = environments(ψ, H, ψ); - imaginary_evolution::Bool = false, normalize::Bool = false - ) - return _timestep_finite!(ψ, H, t, dt, alg, envs; imaginary_evolution, normalize) -end - # Single-site integrators visit every site, finishing a half-sweep by evolving only the center # tensor at its far end. _sweep_ranges(::TDVP, ψ) = (1:length(ψ), length(ψ):-1:1) @@ -240,33 +220,21 @@ function sweep!(alg, state::TimeStepState, direction) bond = _swept_bond(alg, site, ψ, direction) isnothing(bond) || (state.truncation_errors[bond] = ϵ) end - return _with_mps(state, ψ) + # each half-sweep evolves the state over `dt / 2` + return _with_mps(state, ψ; t = state.t + state.dt / 2) end # a step is a left→right half-sweep `t → t + dt / 2` followed by a right→left half-sweep # `t + dt / 2 → t + dt` -function _timestep_finite!( - ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg, envs; - imaginary_evolution::Bool, normalize::Bool, timeroutput = NoTimerOutput() - ) - state = TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput) - state = _timestep_finite!(alg, state) - return state.mps, state.envs, _step_info(state) -end -function _timestep_finite!(alg, state::TimeStepState) +function _timestep!(alg, state::TimeStepState{<:AbstractFiniteMPS}) state = sweep!(alg, state, Val(:right)) return sweep!(alg, state, Val(:left)) end -# 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, ψ, alg::TDVP, state::TimeStepState) (; operator, envs, t, dt, allocator, imaginary_evolution, normalize) = state H = operator - t_AC, t_C = _half_sweep_times(direction, t, dt) + t_AC, t_C = t, t + dt / 2 # at the far end of the sweep there is no bond ahead: only evolve the center tensor if site == _sweep_end(ψ, direction) @@ -283,9 +251,8 @@ function local_update!(site, direction::Val, ψ, alg::TDVP, state::TimeStepState Hac = AC_hamiltonian(site, ψ, H, ψ, envs; alg.backend, allocator) AC = integrate(Hac, ψ.AC[site], t_AC, 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. + # 3. gauge: split AC onto the bond ahead and move the center. + # By default the norm is preserved; `normalize` renormalizes. if direction === Val(:right) _, ϵ = left_gauge!(ψ, site, AC, alg.alg_gauge; normalize) bond = site @@ -342,14 +309,6 @@ Used as the `algorithm` argument of [`timestep`](@ref), [`timestep!`](@ref) and backend::B = Defaults.backend() end -function timestep!( - ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, - envs::AbstractMPSEnvironments = environments(ψ, H, ψ); - imaginary_evolution::Bool = false, normalize::Bool = false - ) - return _timestep_finite!(ψ, H, t, dt, alg, envs; imaginary_evolution, normalize) -end - # two-site TDVP visits every bond, and each update cuts the bond it acts on _sweep_ranges(::TDVP2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 1):-1:1) _swept_bond(::TDVP2, pos, ψ, ::Val{:right}) = pos @@ -358,7 +317,7 @@ _swept_bond(::TDVP2, pos, ψ, ::Val{:left}) = pos function local_update!(pos, direction::Val, ψ, alg::TDVP2, state::TimeStepState) (; operator, envs, t, dt, allocator, imaginary_evolution, normalize) = state H = operator - t_AC2, t_AC = _half_sweep_times(direction, t, dt) + t_AC2, t_AC = t, t + dt / 2 # 1. evolve the two-site center tensor at `(pos, pos + 1)` forward ac2 = if direction === Val(:right) @@ -384,32 +343,5 @@ function local_update!(pos, direction::Val, ψ, alg::TDVP2, state::TimeStepState return ψ, ϵ end -# Copy of `ψ` to step out of place, promoted to complex when the evolution leaves the reals, -# together with environments for that copy. -function _timestep_copy(ψ, H, envs, promote::Bool) - ψ′ = promote ? complex(ψ) : copy(ψ) - if length(envs) != 0 && promote - @warn "Currently cannot reuse real environments for complex evolution" - envs′ = environments(ψ′, H, ψ′) - elseif length(envs) == 1 - envs′ = only(envs) - else - @assert length(envs) == 0 "Invalid signature" - envs′ = environments(ψ′, H, ψ′) - end - return ψ′, envs′ -end - _promotes_to_complex(ψ, H, dt, ::Union{TDVP, TDVP2}; imaginary_evolution) = scalartype(ψ) <: Real && !(imaginary_evolution && isreal(dt)) - -# copying version -function timestep( - ψ::Union{AbstractFiniteMPS, InfiniteMPS}, H, time::Number, timestep::Number, - alg::Union{TDVP, TDVP2}, envs::AbstractMPSEnvironments...; - imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... - ) - promote = _promotes_to_complex(ψ, H, timestep, alg; imaginary_evolution) - ψ′, envs′ = _timestep_copy(ψ, H, envs, promote) - return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...) -end diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index a149b2664..866bb8230 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -36,9 +36,7 @@ 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`. Integrators that step in place carry -# their `TimeStepState` over in `step` instead, or `nothing` before the first step. +# Internal state of `time_evolve` struct TimeEvolveState{S, O, E, T, F, K, P} mps::S operator::O @@ -66,19 +64,8 @@ function Base.iterate(it::IterativeSolver{<:Any, <:TimeEvolveState}, state::Time return (ψ, envs, info), it.state end -_evolve_step(alg, state::TimeEvolveState, t, dt) = _evolve_timestep(alg, state, t, dt) -function _evolve_timestep(alg, state::TimeEvolveState, t, dt) - ψ, envs, info = state.stepper( - state.mps, state.operator, t, dt, alg, state.envs; state.kwargs... - ) - return ψ, envs, info, nothing -end - -# An integrator that steps in place keeps its `TimeStepState`, and with it its scheduler and -# allocator, across steps. Only the first step of `time_evolve` copies the state (promoting it to -# complex if needed), after which every step is taken in place. -function _evolve_step(alg::Union{TDVP, TDVP2, BUG}, state::TimeEvolveState, t, dt) - _steps_in_place(alg, state.mps) || return _evolve_timestep(alg, state, t, dt) +# only the first step copies the state, if `stepper` is `timestep` +function _evolve_step(alg, state::TimeEvolveState, t, dt) (; imaginary_evolution, normalize) = state.kwargs H = state.operator step = if isnothing(state.step) @@ -206,6 +193,43 @@ julia> round(real(expectation_value(ψ, 2 => X)); digits = 6) """ function timestep end, function timestep! end +function timestep!( + ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg, + envs::AbstractMPSEnvironments = environments(ψ, H, ψ); + imaginary_evolution::Bool = false, normalize::Bool = false, + timeroutput = NoTimerOutput() + ) + state = TimeStepState(ψ, H, envs, t, dt; imaginary_evolution, normalize, timeroutput) + state = _timestep!(alg, state) + return state.mps, state.envs, _step_info(state) +end + +function timestep( + ψ::Union{AbstractFiniteMPS, InfiniteMPS}, H, time::Number, timestep::Number, + alg, envs::AbstractMPSEnvironments...; + imaginary_evolution::Bool = false, normalize::Bool = false, kwargs... + ) + promote = _promotes_to_complex(ψ, H, timestep, alg; imaginary_evolution) + ψ′, envs′ = _timestep_copy(ψ, H, envs, promote) + return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, normalize, kwargs...) +end + +# Copy of `ψ` to step out of place, promoted to complex when the evolution leaves the reals, +# together with environments for that copy. +function _timestep_copy(ψ, H, envs, promote::Bool) + ψ′ = promote ? complex(ψ) : copy(ψ) + if length(envs) != 0 && promote + @warn "Currently cannot reuse real environments for complex evolution" + envs′ = environments(ψ′, H, ψ′) + elseif length(envs) == 1 + envs′ = only(envs) + else + @assert length(envs) == 0 "Invalid signature" + envs′ = environments(ψ′, H, ψ′) + end + return ψ′, envs′ +end + @doc """ make_time_mpo(H::MPOHamiltonian, dt::Number, alg; kwargs...) -> O::MPO diff --git a/test/timeevolution/timestep.jl b/test/timeevolution/timestep.jl index 9c42207c8..75d282957 100644 --- a/test/timeevolution/timestep.jl +++ b/test/timeevolution/timestep.jl @@ -82,22 +82,20 @@ maxbond(ψ) = maximum(i -> dim(left_virtualspace(ψ, i)), 1:length(ψ)) # the AC and C sweeps of the infinite integrator spawn over the unit cell, and additionally run # concurrently with each other, so they share one allocator: check both schedulers agree @testset "Infinite TDVP ($schedname)" for (schedname, scheduler) in SCHEDULERS - ψ, envs = with_scheduler(scheduler) do - return timestep(ψ₀, H, 0.0, dt, TDVP()) + with_scheduler(scheduler) do + ψ, envs = timestep(ψ₀, H, 0.0, dt, TDVP()) + E = expectation_value(ψ, H, envs) + @test E₀ ≈ E atol = 1.0e-2 + + ψ′ = copy(ψ₀) + ψ″, = timestep!(ψ′, H, 0.0, dt, TDVP()) + @test ψ″ === ψ′ + @test ψ″.AC == ψ.AC + + # evolving in the right-orthogonal gauge gives the same state + ψr, = timestep(ψ₀, H, 0.0, dt, TDVP(); leftorthflag = false) + @test abs(dot(ψ, ψr)) ≈ 1 atol = 1.0e-8 end - E = expectation_value(ψ, H, envs) - @test E₀ ≈ E atol = 1.0e-2 - end - - @testset "Infinite TDVP in place" begin - ψ, = timestep(ψ₀, H, 0.0, dt, TDVP()) - ψ′ = copy(ψ₀) - ψ″, = timestep!(ψ′, H, 0.0, dt, TDVP()) - @test ψ″ === ψ′ - @test ψ″.AC == ψ.AC - # evolving in the right-orthogonal gauge gives the same state - ψr, = timestep(ψ₀, H, 0.0, dt, TDVP(); leftorthflag = false) - @test abs(dot(ψ, ψr)) ≈ 1 atol = 1.0e-8 end @testset "Infinite real TDVP" begin