Iterator-based sweep solvers

docs/src/changelog.md

@@ -27,6 +27,16 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 
 ### Changed
 
+- Internal: the finite `DMRG`/`DMRG2` ground-state search and approximation, `DynamicalDMRG`,
+  multiline `IDMRG`/`IDMRG2` approximation and `leading_boundary`, and `time_evolve` now run as
+  iterators whose iterations are a forward and a backward half-sweep (or a single time step),
+  sharing one structure across solvers. Finite `TDVP`/`TDVP2` steps are composed of per-site local
+  updates. Results and logging are unchanged. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
+- The `IDMRG`/`IDMRG2` convergence measure (`bondresidual`) is computed by one shared,
+  lower-allocation routine. For multiline `IDMRG` with several rows, the changes of the rows are
+  now summed, as `IDMRG2` already did, instead of combined in a 2-norm, which makes the stopping
+  criterion marginally stricter. ([#537](https://github.com/QuantumKitHub/MPSKit.jl/pull/537))
+
 ### Deprecated
 
 ### Removed

src/algorithms/approximate/fvomps.jl

@@ -1,74 +1,82 @@
-function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG2, envs = environments(ψ, _environment_args(Oϕ)...))
+# Internal state of the finite DMRG/DMRG2 approximation, where `ϵ` is the largest relative
+# local change of the last sweep
+struct ApproximateState{S, O, E, A}
+    mps::S
+    operator::O
+    envs::E
+    iter::Int
+    ϵ::Float64
+    allocator::A
+end
+
+function approximate!(
+        ψ::AbstractFiniteMPS, Oϕ, alg::Union{DMRG, DMRG2},
+        envs = environments(ψ, _environment_args(Oϕ)...)
+    )
     allocator = default_allocator(ψ, SerialScheduler())
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
-    log = IterLog("DMRG2")
+    log = IterLog(alg)
+    it = IterativeSolver(alg, ApproximateState(ψ, Oϕ, envs, 0, 2 * alg.tol, allocator))
 
     with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            ϵ = 0.0
-            for pos in [1:(length(ψ) - 1); (length(ψ) - 2):-1:1]
-                AC2′ = AC2_projection(pos, ψ, Oϕ, envs; alg.backend, allocator)
-                al, c, ar, = svd_trunc!(AC2′, inner_alg_gauge(alg))
-
-                AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
-                ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2))
-
-                ψ.AC[pos] = (al, complex(c))
-                ψ.AC[pos + 1] = (complex(c), _transpose_front(ar))
-            end
-
-            # finalize
-            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ), typeof(envs)}
-
+        @log_initialization loginit!(log, it.ϵ)
+        for (_, _, ϵ) in Iterators.take(it, alg.maxiter)
             if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
+                @log_convergence logfinish!(log, it.iter, ϵ)
                 break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
+            elseif it.iter == alg.maxiter
+                @log_nonconvergence logcancel!(log, it.iter, ϵ)
             else
-                @log_iteration logiter!(log, iter, ϵ)
+                @log_iteration logiter!(log, it.iter, ϵ)
             end
         end
     end
 
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, localchange = ϵ, numiter = iter)
+    state = it.state
+    info = AlgorithmInfo(; converged = state.ϵ <= alg.tol, localchange = state.ϵ, numiter = state.iter)
+    return state.mps, state.envs, info
 end
 
-function approximate!(ψ::AbstractFiniteMPS, Oϕ, alg::DMRG, envs = environments(ψ, _environment_args(Oϕ)...))
-    allocator = default_allocator(ψ, SerialScheduler())
-    ϵ::Float64 = 2 * alg.tol
-    iter = 0
-    log = IterLog("DMRG")
-
-    with_verbosity(; alg.verbosity) do
-        @log_initialization loginit!(log, ϵ)
-        for outer iter in 1:(alg.maxiter)
-            ϵ = 0.0
-            for pos in [1:(length(ψ) - 1); length(ψ):-1:2]
-                AC′ = AC_projection(pos, ψ, Oϕ, envs; alg.backend, allocator)
-                AC = ψ.AC[pos]
-                ϵ = max(ϵ, norm(AC′ - AC) / norm(AC′))
+function Base.iterate(it::IterativeSolver{<:Union{DMRG, DMRG2}}, state::ApproximateState)
+    iter = state.iter + 1
+    state = ApproximateState(
+        state.mps, state.operator, state.envs, state.iter, 0.0, state.allocator
+    )
+    state = sweep!(it, state, Val(:right), iter)
+    state = sweep!(it, state, Val(:left), iter)
+    ψ, envs = it.finalize(
+        iter, state.mps, state.operator, state.envs
+    )::Tuple{typeof(state.mps), typeof(state.envs)}
+    it.state = ApproximateState(ψ, state.operator, envs, iter, state.ϵ, state.allocator)
+    return (ψ, envs, state.ϵ), it.state
+end
 
-                ψ.AC[pos] = AC′
-            end
+function sweep!(it::IterativeSolver{<:DMRG2}, state::ApproximateState, direction, iter)
+    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    sites = direction === Val(:right) ? fwd : bwd
+    ψ, Oϕ, envs, ϵ = state.mps, state.operator, state.envs, state.ϵ
+    for pos in sites
+        AC2′ = AC2_projection(pos, ψ, Oϕ, envs; it.alg.backend, state.allocator)
+        al, c, ar, = svd_trunc!(AC2′, inner_alg_gauge(it.alg))
 
-            # finalize
-            ψ, envs = alg.finalize(iter, ψ, Oϕ, envs)::Tuple{typeof(ψ), typeof(envs)}
+        AC2 = ψ.AC[pos] * _transpose_tail(ψ.AR[pos + 1])
+        ϵ = max(ϵ, norm(al * c * ar - AC2) / norm(AC2))
 
-            if ϵ <= alg.tol
-                @log_convergence logfinish!(log, iter, ϵ)
-                break
-            end
-            if iter == alg.maxiter
-                @log_nonconvergence logcancel!(log, iter, ϵ)
-            else
-                @log_iteration logiter!(log, iter, ϵ)
-            end
-        end
+        ψ.AC[pos] = (al, complex(c))
+        ψ.AC[pos + 1] = (complex(c), _transpose_front(ar))
     end
+    return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator)
+end
 
-    return ψ, envs, AlgorithmInfo(; converged = ϵ <= alg.tol, localchange = ϵ, numiter = iter)
+function sweep!(it::IterativeSolver{<:DMRG}, state::ApproximateState, direction, iter)
+    fwd, bwd = _sweep_ranges(it.alg, state.mps)
+    sites = direction === Val(:right) ? fwd : bwd
+    ψ, Oϕ, envs, ϵ = state.mps, state.operator, state.envs, state.ϵ
+    for pos in sites
+        AC′ = AC_projection(pos, ψ, Oϕ, envs; it.alg.backend, state.allocator)
+        AC = ψ.AC[pos]
+        ϵ = max(ϵ, norm(AC′ - AC) / norm(AC′))
+
+        ψ.AC[pos] = AC′
+    end
+    return ApproximateState(ψ, Oϕ, envs, state.iter, ϵ, state.allocator)
 end