diff --git a/examples/heisenberg_su/main.jl b/examples/heisenberg_su/main.jl
index a77405d3f..d028c459c 100644
--- a/examples/heisenberg_su/main.jl
+++ b/examples/heisenberg_su/main.jl
@@ -70,12 +70,11 @@ fix a truncation error (if that can be reached by remaining below `Dbond`):
 
 dts = [1.0e-2, 1.0e-3, 4.0e-4]
 tols = [1.0e-6, 1.0e-8, 1.0e-8]
-maxiter = 10000
+nstep = 10000
 trunc_peps = truncerror(; atol = 1.0e-10) & truncrank(Dbond)
-
+alg = SimpleUpdate(; trunc = trunc_peps, bipartite = true)
 for (dt, tol) in zip(dts, tols)
-    alg = SimpleUpdate(dt, tol, maxiter, trunc_peps)
-    global peps, wts, = simpleupdate(peps, H, alg, wts; bipartite = true)
+    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts; tol, check_interval = 500)
 end
 
 md"""
diff --git a/examples/hubbard_su/main.jl b/examples/hubbard_su/main.jl
index 68c94cc2f..786dab680 100644
--- a/examples/hubbard_su/main.jl
+++ b/examples/hubbard_su/main.jl
@@ -66,26 +66,24 @@ maxiter = 20000
 
 for (dt, tol, Dbond) in zip(dts, tols, Ds)
     trunc = truncerror(; atol = 1.0e-10) & truncrank(Dbond)
-    alg = SimpleUpdate(dt, tol, maxiter, trunc)
-    global peps, wts, = simpleupdate(
-        peps, H, alg, wts; bipartite = false, check_interval = 2000
-    )
+    alg = SimpleUpdate(; trunc, bipartite = false)
+    global peps, wts, = time_evolve(peps, H, dt, maxiter, alg, wts; tol, check_interval = 2000)
 end
 
 md"""
 ## Computing the ground-state energy
 
 In order to compute the energy expectation value with evolved PEPS, we need to converge a
-CTMRG environment on it. We first converge an environment with a small enviroment dimension
-and then use that to initialize another run with bigger environment dimension. We'll use
-`trunc=truncrank(χ)` for that such that the dimension is increased during the second CTMRG
-run:
+CTMRG environment on it. We first converge an environment with a small enviroment dimension,
+which is initialized using the simple update bond weights. Next we use it to initialize 
+another run with bigger environment dimension. The dynamic adjustment of environment dimension
+is achieved by using `trunc=truncrank(χ)` with different `χ`s in the CTMRG runs:
 """
 
 χenv₀, χenv = 6, 16
 env_space = Vect[fℤ₂](0 => χenv₀ / 2, 1 => χenv₀ / 2)
 normalize!.(peps.A, Inf)
-env = CTMRGEnv(rand, Float64, peps, env_space)
+env = CTMRGEnv(wts, peps)
 for χ in [χenv₀, χenv]
     global env, = leading_boundary(
         env, peps; alg = :sequential, tol = 1.0e-8, maxiter = 50, trunc = truncrank(χ)
diff --git a/examples/j1j2_su/main.jl b/examples/j1j2_su/main.jl
index 5dac25897..e027f5bcc 100644
--- a/examples/j1j2_su/main.jl
+++ b/examples/j1j2_su/main.jl
@@ -48,15 +48,16 @@ Therefore, we shall gradually increase $J_2 / J_1$ from 0.1 to 0.5, each time in
 on the previously evolved PEPS:
 """
 
-dt, tol, maxiter = 1.0e-2, 1.0e-8, 30000
+dt, tol, nstep = 1.0e-2, 1.0e-8, 30000
 check_interval = 4000
 trunc_peps = truncerror(; atol = 1.0e-10) & truncrank(Dbond)
-alg = SimpleUpdate(dt, tol, maxiter, trunc_peps)
+alg = SimpleUpdate(; trunc = trunc_peps)
 for J2 in 0.1:0.1:0.5
-    H = real( ## convert Hamiltonian `LocalOperator` to real floats
+    ## convert Hamiltonian `LocalOperator` to real floats
+    H = real(
         j1_j2_model(ComplexF64, symm, InfiniteSquare(Nr, Nc); J1, J2, sublattice = false),
     )
-    global peps, wts, = simpleupdate(peps, H, alg, wts; check_interval)
+    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts; tol, check_interval)
 end
 
 md"""
@@ -69,8 +70,7 @@ tols = [1.0e-9, 1.0e-9]
 J2 = 0.5
 H = real(j1_j2_model(ComplexF64, symm, InfiniteSquare(Nr, Nc); J1, J2, sublattice = false))
 for (dt, tol) in zip(dts, tols)
-    alg′ = SimpleUpdate(dt, tol, maxiter, trunc_peps)
-    global peps, wts, = simpleupdate(peps, H, alg′, wts; check_interval)
+    global peps, wts, = time_evolve(peps, H, dt, nstep, alg, wts; tol)
 end
 
 md"""
diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index b3730aa46..2fb5131be 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -71,6 +71,7 @@ include("algorithms/truncation/fullenv_truncation.jl")
 include("algorithms/truncation/bond_truncation.jl")
 
 include("algorithms/time_evolution/evoltools.jl")
+include("algorithms/time_evolution/time_evolve.jl")
 include("algorithms/time_evolution/simpleupdate.jl")
 include("algorithms/time_evolution/simpleupdate3site.jl")
 
@@ -100,7 +101,8 @@ export fixedpoint
 
 export absorb_weight
 export ALSTruncation, FullEnvTruncation
-export su_iter, su3site_iter, simpleupdate, SimpleUpdate
+export SimpleUpdate
+export TimeEvolver, timestep, time_evolve
 
 export InfiniteSquareNetwork
 export InfinitePartitionFunction
diff --git a/src/algorithms/time_evolution/evoltools.jl b/src/algorithms/time_evolution/evoltools.jl
index 5c5466d52..627e38bf1 100644
--- a/src/algorithms/time_evolution/evoltools.jl
+++ b/src/algorithms/time_evolution/evoltools.jl
@@ -1,4 +1,22 @@
-const InfiniteState = Union{InfinitePEPS, InfinitePEPO}
+"""
+Process the Trotter time step `dt` according to the intended usage.
+"""
+function _get_dt(
+        state::InfiniteState, dt::Number, imaginary_time::Bool
+    )
+    # PEPS update: exp(-H dt)|ψ⟩
+    # PEPO update (purified): exp(-H dt/2)|ρ⟩
+    # PEPO update (not purified): exp(-H dt/2) ρ exp(-H dt/2)
+    dt′ = (state isa InfinitePEPS) ? dt : (dt / 2)
+    if (state isa InfinitePEPO)
+        @assert size(state)[3] == 1
+    end
+    if !imaginary_time
+        @assert (state isa InfinitePEPS) "Real time evolution of InfinitePEPO (Heisenberg picture) is not implemented."
+        dt′ = 1.0im * dt′
+    end
+    return dt′
+end
 
 function MPSKit.infinite_temperature_density_matrix(H::LocalOperator)
     T = scalartype(H)
diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl
index d466483c7..17e3fc46d 100644
--- a/src/algorithms/time_evolution/simpleupdate.jl
+++ b/src/algorithms/time_evolution/simpleupdate.jl
@@ -1,38 +1,96 @@
 """
 $(TYPEDEF)
 
-Algorithm struct for simple update (SU) of infinite PEPS with bond weights.
-Each SU run is converged when the singular value difference becomes smaller than `tol`.
+Algorithm struct for simple update (SU) of InfinitePEPS or InfinitePEPO.
 
 ## Fields
 
 $(TYPEDFIELDS)
 """
-struct SimpleUpdate
-    dt::Number
-    tol::Float64
-    maxiter::Int
+@kwdef struct SimpleUpdate <: TimeEvolution
+    "Truncation strategy for bonds updated by Trotter gates"
     trunc::TruncationStrategy
+    "When true (or false), the Trotter gate is `exp(-H dt)` (or `exp(-iH dt)`)"
+    imaginary_time::Bool = true
+    "When true, force the usage of 3-site simple update"
+    force_3site::Bool = false
+    "(Only applicable to InfinitePEPS) When true, assume bipartite unit cell structure"
+    bipartite::Bool = false
+    "(Only applicable to InfinitePEPO) 
+    When true, the PEPO is regarded as a purified PEPS, and updated as
+    `|ρ(t + dt)⟩ = exp(-H dt/2) |ρ(t)⟩`.
+    When false, the PEPO is updated as 
+    `ρ(t + dt) = exp(-H dt/2) ρ(t) exp(-H dt/2)`."
+    purified::Bool = true
 end
-# TODO: add kwarg constructor and SU Defaults
 
 """
-$(SIGNATURES)
+Internal state of simple update algorithm
+"""
+struct SUState{S <: InfiniteState, E <: SUWeight, N <: Number}
+    "number of performed iterations"
+    iter::Int
+    "evolved time"
+    t::N
+    "PEPS/PEPO"
+    psi::S
+    "SUWeight environment"
+    env::E
+end
 
-Simple update of the x-bond between `[r,c]` and `[r,c+1]`.
+"""
+    TimeEvolver(
+        psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, 
+        alg::SimpleUpdate, env0::SUWeight; t0::Number = 0.0
+    )
+
+Initialize a TimeEvolver with Hamiltonian `H` and simple update `alg`, 
+starting from the initial state `psi0` and SUWeight environment `env0`.
 
+- The initial time is specified by `t0`.
+"""
+function TimeEvolver(
+        psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int,
+        alg::SimpleUpdate, env0::SUWeight; t0::Number = 0.0
+    )
+    # sanity checks
+    _timeevol_sanity_check(psi0, physicalspace(H), alg)
+    # create Trotter gates
+    nnonly = is_nearest_neighbour(H)
+    use_3site = alg.force_3site || !nnonly
+    if alg.bipartite
+        @assert !use_3site "3-site simple update is incompatible with bipartite lattice."
+    end
+    dt′ = _get_dt(psi0, dt, alg.imaginary_time)
+    gate = if use_3site
+        [
+            _get_gatempos_se(H, dt′),
+            _get_gatempos_se(rotl90(H), dt′),
+            _get_gatempos_se(rot180(H), dt′),
+            _get_gatempos_se(rotr90(H), dt′),
+        ]
+    else
+        get_expham(H, dt′)
+    end
+    state = SUState(0, t0, psi0, env0)
+    return TimeEvolver(alg, dt, nstep, gate, state)
+end
+
+"""
+Simple update of the x-bond between `[r,c]` and `[r,c+1]`.
 ```
         |           |
     -- T[r,c] -- T[r,c+1] --
         |           |
 ```
+When `gate_ax = 1` (or `2`), the gate will be applied to 
+the codomain (or domain) physicsl legs of `state`.
 """
 function _su_xbond!(
         state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trunc::TruncationStrategy; gate_ax::Int = 1
     ) where {T <: Number, S <: ElementarySpace}
-    Nr, Nc = size(state)[1:2]
-    @assert 1 <= row <= Nr && 1 <= col <= Nc
+    Nr, Nc, = size(state)
     cp1 = _next(col, Nc)
     # absorb environment weights
     A, B = state.A[row, col], state.A[row, cp1]
@@ -65,13 +123,14 @@ Simple update of the y-bond between `[r,c]` and `[r-1,c]`.
     -- T[r,c] ---
         |
 ```
+When `gate_ax = 1` (or `2`), the gate will be applied to 
+the codomain (or domain) physicsl legs of `state`.
 """
 function _su_ybond!(
         state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::SUWeight,
         row::Int, col::Int, trunc::TruncationStrategy; gate_ax::Int = 1
     ) where {T <: Number, S <: ElementarySpace}
-    Nr, Nc = size(state)[1:2]
-    @assert 1 <= row <= Nr && 1 <= col <= Nc
+    Nr, Nc, = size(state)
     rm1 = _prev(row, Nr)
     # absorb environment weights
     A, B = state.A[row, col], state.A[rm1, col]
@@ -96,36 +155,22 @@ function _su_ybond!(
 end
 
 """
-    su_iter(state::InfinitePEPS, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight; bipartite::Bool = false)
-    su_iter(densitymatrix::InfinitePEPO, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight; gate_bothsides::Bool = true)
-
-One round of simple update, which applies the nearest neighbor `gate` on an InfinitePEPS `state` or InfinitePEPO `densitymatrix`.
-When the input `state` has a unit cell size of (2, 2), one can set `bipartite = true` to enforce the bipartite structure. 
+One iteration of simple update
 """
 function su_iter(
-        state::InfiniteState, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        bipartite::Bool = false, gate_bothsides::Bool = true
+        state::InfiniteState, gate::LocalOperator, alg::SimpleUpdate, env::SUWeight
     )
-    @assert size(gate.lattice) == size(state)[1:2]
-    if state isa InfinitePEPS
-        gate_bothsides = false
-    else
-        @assert size(state, 3) == 1
-    end
-    Nr, Nc = size(state)[1:2]
-    bipartite && (@assert Nr == Nc == 2)
-    (Nr >= 2 && Nc >= 2) || throw(
-        ArgumentError("`state` unit cell size for simple update should be no smaller than (2, 2)."),
-    )
-    state2, env2 = deepcopy(state), deepcopy(env)
-    gate_axs = gate_bothsides ? (1:2) : (1:1)
+    Nr, Nc, = size(state)
+    state2, env2, ϵ = deepcopy(state), deepcopy(env), 0.0
+    gate_axs = alg.purified ? (1:1) : (1:2)
     for r in 1:Nr, c in 1:Nc
         term = get_gateterm(gate, (CartesianIndex(r, c), CartesianIndex(r, c + 1)))
         trunc = truncation_strategy(alg.trunc, 1, r, c)
         for gate_ax in gate_axs
-            _su_xbond!(state2, term, env2, r, c, trunc; gate_ax)
+            ϵ′ = _su_xbond!(state2, term, env2, r, c, trunc; gate_ax)
+            ϵ = max(ϵ, ϵ′)
         end
-        if bipartite
+        if alg.bipartite
             rp1, cp1 = _next(r, Nr), _next(c, Nc)
             state2.A[rp1, cp1] = deepcopy(state2.A[r, c])
             state2.A[rp1, c] = deepcopy(state2.A[r, cp1])
@@ -134,94 +179,134 @@ function su_iter(
         term = get_gateterm(gate, (CartesianIndex(r, c), CartesianIndex(r - 1, c)))
         trunc = truncation_strategy(alg.trunc, 2, r, c)
         for gate_ax in gate_axs
-            _su_ybond!(state2, term, env2, r, c, trunc; gate_ax)
+            ϵ′ = _su_ybond!(state2, term, env2, r, c, trunc; gate_ax)
+            ϵ = max(ϵ, ϵ′)
         end
-        if bipartite
+        if alg.bipartite
             rm1, cm1 = _prev(r, Nr), _prev(c, Nc)
             state2.A[rm1, cm1] = deepcopy(state2.A[r, c])
             state2.A[r, cm1] = deepcopy(state2.A[rm1, c])
             env2.data[2, rm1, cm1] = deepcopy(env2.data[2, r, c])
         end
     end
-    return state2, env2
+    return state2, env2, ϵ
+end
+
+function Base.iterate(it::TimeEvolver{<:SimpleUpdate}, state = it.state)
+    iter, t = state.iter, state.t
+    (iter == it.nstep) && return nothing
+    psi, env, ϵ = su_iter(state.psi, it.gate, it.alg, state.env)
+    # update internal state
+    iter += 1
+    t += it.dt
+    it.state = SUState(iter, t, psi, env)
+    info = (; t, ϵ)
+    return (psi, env, info), it.state
+end
+
+"""
+    timestep(
+        it::TimeEvolver{<:SimpleUpdate}, psi::InfiniteState, env::SUWeight
+    ) -> (psi, env, info)
+
+Given the TimeEvolver iterator `it`, perform one step of time evolution
+on the input state `psi` and its environment `env`.
+"""
+function MPSKit.timestep(
+        it::TimeEvolver{<:SimpleUpdate}, psi::InfiniteState, env::SUWeight
+    )
+    _timeevol_sanity_check(psi, physicalspace(it.state.psi), it.alg)
+    state = SUState(it.state.iter, it.state.t, psi, env)
+    result = iterate(it, state)
+    if result === nothing
+        @warn "TimeEvolver `it` has already reached the end."
+        return nothing
+    else
+        return first(result)
+    end
 end
 
 """
-Perform simple update with Hamiltonian `ham` containing up to nearest neighbor interaction terms. 
+    time_evolve(
+        it::TimeEvolver{<:SimpleUpdate}; 
+        tol::Float64 = 0.0, check_interval::Int = 500
+    ) -> (psi, env, info)
+
+Perform time evolution to the end of TimeEvolver iterator `it`,
+or until convergence of SUWeight set by a positive `tol`.
+
+- Setting `tol > 0` enables convergence check (for imaginary time evolution of InfinitePEPS only).
+    For other usages it should not be changed.
+- `check_interval` sets the number of iterations between outputs of information.
 """
-function _simpleupdate2site(
-        state::InfiniteState, ham::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        bipartite::Bool = false, check_interval::Int = 500, gate_bothsides::Bool = true
+function MPSKit.time_evolve(
+        it::TimeEvolver{<:SimpleUpdate};
+        tol::Float64 = 0.0, check_interval::Int = 500
     )
     time_start = time()
-    # exponentiating the 2-site Hamiltonian gate
-    if state isa InfinitePEPS
-        gate_bothsides = false
+    check_convergence = (tol > 0)
+    if check_convergence
+        @assert (it.state.psi isa InfinitePEPS) && it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking."
     end
-    dt = gate_bothsides ? (alg.dt / 2) : alg.dt
-    gate = get_expham(ham, dt)
-    wtdiff = 1.0
-    env0 = deepcopy(env)
-    for count in 1:(alg.maxiter)
-        time0 = time()
-        state, env = su_iter(state, gate, alg, env; bipartite, gate_bothsides)
-        wtdiff = compare_weights(env, env0)
-        converge = (wtdiff < alg.tol)
-        cancel = (count == alg.maxiter)
-        env0 = deepcopy(env)
+    env0, time0 = it.state.env, time()
+    for (psi, env, info) in it
+        iter = it.state.iter
+        diff = compare_weights(env0, env)
+        stop = (iter == it.nstep) || (diff < tol)
+        showinfo = (check_interval > 0) &&
+            ((iter % check_interval == 0) || (iter == 1) || stop)
         time1 = time()
-        if ((count == 1) || (count % check_interval == 0) || converge || cancel)
+        if showinfo
             @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
-            label = (converge ? "conv" : (cancel ? "cancel" : "iter"))
-            message = @sprintf(
-                "SU %s %-7d:  dt = %.0e,  weight diff = %.3e,  time = %.3f sec\n",
-                label, count, alg.dt, wtdiff, time1 - ((converge || cancel) ? time_start : time0)
+            @info @sprintf(
+                "SU iter %-7d: dt = %s, |Δλ| = %.3e. Time = %.3f s/it",
+                # using `string` since `dt` can be complex
+                iter, string(it.dt), diff, time1 - time0
             )
-            cancel ? (@warn message) : (@info message)
         end
-        converge && break
+        if check_convergence
+            if (iter == it.nstep) && (diff >= tol)
+                @warn "SU: bond weights have not converged."
+            end
+            if diff < tol
+                @info "SU: bond weights have converged."
+            end
+        end
+        if stop
+            time_end = time()
+            @info @sprintf("Simple update finished. Total time elasped: %.2f s", time_end - time_start)
+            return psi, env, info
+        else
+            env0 = env
+        end
+        time0 = time()
     end
-    return state, env, wtdiff
+    return
 end
 
 """
-    simpleupdate(
-        state::InfinitePEPS, ham::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        bipartite::Bool = false, force_3site::Bool = false, check_interval::Int = 500
-    )
-    simpleupdate(
-        densitymatrix::InfinitePEPO, ham::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        gate_bothsides::Bool = true, force_3site::Bool = false, check_interval::Int = 500
-    )
+    time_evolve(
+        psi0::Union{InfinitePEPS, InfinitePEPO}, H::LocalOperator, 
+        dt::Number, nstep::Int, alg::SimpleUpdate, env0::SUWeight;
+        tol::Float64 = 0.0, t0::Number = 0.0, check_interval::Int = 500
+    ) -> (psi, env, info)
 
-Perform a simple update on an InfinitePEPS `state` or an InfinitePEPO `densitymatrix`
-using the Hamiltonian `ham`, which can contain up to next-nearest-neighbor interaction terms.
+Perform time evolution on the initial state `psi0` and initial environment `env0`
+with Hamiltonian `H`, using SimpleUpdate algorithm `alg`, time step `dt` for 
+`nstep` number of steps. 
 
-## Keyword Arguments
-
-- `bipartite::Bool=false`: If `true`, enforces the bipartite structure of the PEPS. 
-- `gate_bothsides::Bool=true`: (Effective only for PEPO evolution) If true, apply the Trotter gate to both side of the PEPO. Otherwise, the gate is only applied to the physical codomain legs. 
-- `force_3site::Bool=false`: Forces the use of the 3-site simple update algorithm, even if the Hamiltonian contains only nearest-neighbor terms.
-- `check_interval::Int=500`: Specifies the number of evolution steps between printing progress information.
-
-## Notes
-
-- Setting `bipartite = true` is allowed only for PEPS evolution with up to next-nearest neighbor terms, and requires the input `peps` to have a unit cell size of (2, 2). 
+- Setting `tol > 0` enables convergence check (for imaginary time evolution of InfinitePEPS only).
+    For other usages it should not be changed.
+- Use `t0` to specify the initial time of the evolution.
+- `check_interval` sets the interval to output information. Output during the evolution can be turned off by setting `check_interval <= 0`.
+- `info` is a NamedTuple containing information of the evolution, 
+    including the time `info.t` evolved since `psi0`.
 """
-function simpleupdate(
-        state::InfiniteState, ham::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        bipartite::Bool = false, gate_bothsides::Bool = true,
-        force_3site::Bool = false, check_interval::Int = 500
+function MPSKit.time_evolve(
+        psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int,
+        alg::SimpleUpdate, env0::SUWeight;
+        tol::Float64 = 0.0, t0::Number = 0.0, check_interval::Int = 500
     )
-    # determine if Hamiltonian contains nearest neighbor terms only
-    nnonly = is_nearest_neighbour(ham)
-    use_3site = force_3site || !nnonly
-    @assert !(state isa InfinitePEPO && bipartite) "Evolution of PEPO with bipartite structure is not implemented."
-    @assert !(bipartite && use_3site) "3-site simple update is incompatible with bipartite lattice."
-    # TODO: check SiteDependentTruncation is compatible with bipartite structure
-    if use_3site
-        return _simpleupdate3site(state, ham, alg, env; gate_bothsides, check_interval)
-    else
-        return _simpleupdate2site(state, ham, alg, env; bipartite, gate_bothsides, check_interval)
-    end
+    it = TimeEvolver(psi0, H, dt, nstep, alg, env0; t0)
+    return time_evolve(it; tol, check_interval)
 end
diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl
index e1dae4e16..7b80acde3 100644
--- a/src/algorithms/time_evolution/simpleupdate3site.jl
+++ b/src/algorithms/time_evolution/simpleupdate3site.jl
@@ -462,7 +462,7 @@ end
 function _su3site_se!(
         state::InfiniteState, gs::Vector{T}, env::SUWeight,
         row::Int, col::Int, truncs::Vector{E};
-        gate_bothsides::Bool = true
+        purified::Bool = true
     ) where {T <: AbstractTensorMap, E <: TruncationStrategy}
     Nr, Nc = size(state)
     @assert 1 <= row <= Nr && 1 <= col <= Nc
@@ -477,7 +477,7 @@ function _su3site_se!(
     # weights in the cluster
     wt_idxs = ((1, row, col), (2, row, cp1))
     # apply gate MPOs
-    gate_axs = gate_bothsides ? (1:2) : (1:1)
+    gate_axs = purified ? (1:1) : (1:2)
     ϵs = nothing
     for gate_ax in gate_axs
         _apply_gatempo!(Ms, gs; gate_ax)
@@ -500,23 +500,13 @@ function _su3site_se!(
         M = absorb_weight(M, env, coord[1], coord[2], openaxs; inv = true)
         state.A[CartesianIndex(coord)] = normalize(M, Inf)
     end
-    return ϵs
+    return maximum(ϵs)
 end
 
-"""
-    su3site_iter(state::InfinitePEPS, gatempos, alg::SimpleUpdate, env::SUWeight)
-    su3site_iter(densitymatrix::InfinitePEPO, gatempos, alg::SimpleUpdate, env::SUWeight; gate_bothsides::Bool = true)
-
-One round of 3-site simple update, which applies the Trotter gate MPOs `gatempos`
-on an InfinitePEPS `state` or InfinitePEPO `densitymatrix`.
-"""
-function su3site_iter(
-        state::InfiniteState, gatempos::Vector{G}, alg::SimpleUpdate, env::SUWeight;
-        gate_bothsides::Bool = true
+function su_iter(
+        state::InfiniteState, gatempos::Vector{G}, alg::SimpleUpdate, env::SUWeight
     ) where {G <: AbstractMatrix}
-    if state isa InfinitePEPS
-        gate_bothsides = false
-    else
+    if state isa InfinitePEPO
         @assert size(state, 3) == 1
     end
     Nr, Nc = size(state)[1:2]
@@ -525,7 +515,7 @@ function su3site_iter(
             "iPEPS unit cell size for simple update should be no smaller than (2, 2)."
         ),
     )
-    state2, env2 = deepcopy(state), deepcopy(env)
+    state2, env2, ϵ = deepcopy(state), deepcopy(env), 0.0
     trunc = alg.trunc
     for i in 1:4
         Nr, Nc = size(state2)[1:2]
@@ -535,53 +525,10 @@ function su3site_iter(
                 truncation_strategy(trunc, 1, r, c)
                 truncation_strategy(trunc, 2, r, _next(c, Nc))
             ]
-            _su3site_se!(state2, gs, env2, r, c, truncs; gate_bothsides)
+            ϵ = _su3site_se!(state2, gs, env2, r, c, truncs; alg.purified)
         end
         state2, env2 = rotl90(state2), rotl90(env2)
         trunc = rotl90(trunc)
     end
-    return state2, env2
-end
-
-"""
-Perform 3-site simple update for Hamiltonian `ham`.
-"""
-function _simpleupdate3site(
-        state::InfiniteState, ham::LocalOperator, alg::SimpleUpdate, env::SUWeight;
-        check_interval::Int = 500, gate_bothsides::Bool = true
-    )
-    time_start = time()
-    # convert Hamiltonian to 3-site exponentiated gate MPOs
-    if state isa InfinitePEPS
-        gate_bothsides = false
-    end
-    dt = gate_bothsides ? (alg.dt / 2) : alg.dt
-    gatempos = [
-        _get_gatempos_se(ham, dt),
-        _get_gatempos_se(rotl90(ham), dt),
-        _get_gatempos_se(rot180(ham), dt),
-        _get_gatempos_se(rotr90(ham), dt),
-    ]
-    wtdiff = 1.0
-    env0 = deepcopy(env)
-    for count in 1:(alg.maxiter)
-        time0 = time()
-        state, env = su3site_iter(state, gatempos, alg, env; gate_bothsides)
-        wtdiff = compare_weights(env, env0)
-        converge = (wtdiff < alg.tol)
-        cancel = (count == alg.maxiter)
-        env0 = deepcopy(env)
-        time1 = time()
-        if ((count == 1) || (count % check_interval == 0) || converge || cancel)
-            @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))"
-            label = (converge ? "conv" : (cancel ? "cancel" : "iter"))
-            message = @sprintf(
-                "SU %s %-7d:  dt = %.0e,  weight diff = %.3e,  time = %.3f sec\n",
-                label, count, alg.dt, wtdiff, time1 - ((converge || cancel) ? time_start : time0)
-            )
-            cancel ? (@warn message) : (@info message)
-        end
-        converge && break
-    end
-    return state, env, wtdiff
+    return state2, env2, ϵ
 end
diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl
new file mode 100644
index 000000000..807caaa02
--- /dev/null
+++ b/src/algorithms/time_evolution/time_evolve.jl
@@ -0,0 +1,47 @@
+"""
+$(TYPEDEF)
+
+Abstract super type for time evolution algorithms of InfinitePEPS or InfinitePEPO.
+"""
+abstract type TimeEvolution end
+
+"""
+    mutable struct TimeEvolver{TE <: TimeEvolution, G, S, N <: Number}
+
+Iterator for Trotter-based time evolution of InfinitePEPS or InfinitePEPO.
+
+## Fields
+
+$(TYPEDFIELDS)
+"""
+mutable struct TimeEvolver{TE <: TimeEvolution, G, S, N <: Number}
+    "Time evolution algorithm (currently supported: `SimpleUpdate`)"
+    alg::TE
+    "Trotter time step"
+    dt::N
+    "The number of iteration steps"
+    nstep::Int
+    "Trotter gates"
+    gate::G
+    "Internal state of the iterator, including the number of
+    already performed iterations, evolved time, PEPS/PEPO and its environment"
+    state::S
+end
+
+Base.iterate(it::TimeEvolver) = iterate(it, it.state)
+
+function _timeevol_sanity_check(
+        ψ₀::InfiniteState, Pspaces::M, alg::A
+    ) where {A <: TimeEvolution, M <: AbstractMatrix{<:ElementarySpace}}
+    Nr, Nc, = size(ψ₀)
+    @assert (Nr >= 2 && Nc >= 2) "Unit cell size for simple update should be no smaller than (2, 2)."
+    @assert Pspaces == physicalspace(ψ₀) "Physical spaces of `ψ₀` do not match `Pspaces`."
+    if hasfield(typeof(alg), :purified) && !alg.purified
+        @assert ψ₀ isa InfinitePEPO "alg.purified = false is only applicable to PEPO."
+    end
+    if hasfield(typeof(alg), :bipartite) && alg.bipartite
+        @assert Nr == Nc == 2 "`bipartite = true` requires 2 x 2 unit cell size."
+        @assert ψ₀ isa InfinitePEPS "Evolution of PEPO with bipartite structure is not implemented."
+    end
+    return nothing
+end
diff --git a/src/algorithms/truncation/bond_truncation.jl b/src/algorithms/truncation/bond_truncation.jl
index 25fcc190f..b2f54ac1c 100644
--- a/src/algorithms/truncation/bond_truncation.jl
+++ b/src/algorithms/truncation/bond_truncation.jl
@@ -15,7 +15,7 @@ The truncation algorithm can be constructed from the following keyword arguments
 
 * `trunc::TruncationStrategy`: SVD truncation strategy when initilizing the truncated tensors connected by the bond.
 * `maxiter::Int=50` : Maximal number of ALS iterations.
-* `tol::Float64=1e-15` : ALS converges when fidelity change between two FET iterations is smaller than `tol`.
+* `tol::Float64=1e-15` : ALS converges when fidelity change between two iterations is smaller than `tol`.
 * `check_interval::Int=0` : Set number of iterations to print information. Output is suppressed when `check_interval <= 0`. 
 """
 @kwdef struct ALSTruncation
diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl
index e16b4cfca..c1aaa1c3a 100644
--- a/src/operators/infinitepepo.jl
+++ b/src/operators/infinitepepo.jl
@@ -27,6 +27,8 @@ struct InfinitePEPO{T <: PEPOTensor}
     end
 end
 
+const InfiniteState = Union{InfinitePEPS, InfinitePEPO}
+
 ## Constructors
 
 """
diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl
index 967b02c17..1a531032a 100644
--- a/test/bondenv/benv_fu.jl
+++ b/test/bondenv/benv_fu.jl
@@ -13,8 +13,9 @@ function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0)
     Vphy = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1 // 2) => 1, (1, -1 // 2) => 1)
     peps = InfinitePEPS(rand, ComplexF64, Vphy, Vphy; unitcell = (Nr, Nc))
     wts = SUWeight(peps)
-    alg = SimpleUpdate(1.0e-2, 1.0e-8, 10000, truncerror(; atol = 1.0e-10) & truncrank(4))
-    peps, = simpleupdate(peps, H, alg, wts; bipartite = false, check_interval = 2000)
+    alg = SimpleUpdate(; trunc = truncerror(; atol = 1.0e-10) & truncrank(4))
+    evolver = TimeEvolver(peps, H, 1.0e-2, 10000, alg, wts)
+    peps, = time_evolve(evolver; tol = 1.0e-8, check_interval = 2000)
     normalize!.(peps.A, Inf)
     return peps
 end
diff --git a/test/examples/heisenberg.jl b/test/examples/heisenberg.jl
index 6fab3f4f2..df526e9f4 100644
--- a/test/examples/heisenberg.jl
+++ b/test/examples/heisenberg.jl
@@ -67,12 +67,12 @@ end
     # simple update
     dts = [1.0e-2, 1.0e-3, 1.0e-3, 1.0e-4]
     tols = [1.0e-7, 1.0e-8, 1.0e-8, 1.0e-8]
-    maxiter = 5000
+    nstep = 5000
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         Dbond2 = (n == 2) ? Dbond + 2 : Dbond
         trunc = truncerror(; atol = 1.0e-10) & truncrank(Dbond2)
-        alg = SimpleUpdate(dt, tol, maxiter, trunc)
-        peps, wts, = simpleupdate(peps, ham, alg, wts; bipartite = false)
+        alg = SimpleUpdate(; trunc, bipartite = false)
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; tol)
     end
 
     # measure physical quantities with CTMRG
diff --git a/test/examples/j1j2_finiteT.jl b/test/examples/j1j2_finiteT.jl
index 4df82ef25..991521c35 100644
--- a/test/examples/j1j2_finiteT.jl
+++ b/test/examples/j1j2_finiteT.jl
@@ -1,12 +1,9 @@
 using Test
-using Random
 using LinearAlgebra
 using TensorKit
 import MPSKitModels: σˣ, σᶻ
 using PEPSKit
 
-Random.seed!(10235876)
-
 # Benchmark energy from high-temperature expansion
 # at β = 0.3, 0.6
 # Physical Review B 86, 045139 (2012) Fig. 15-16
@@ -14,7 +11,7 @@ bm = [-0.1235, -0.213]
 
 function converge_env(state, χ::Int)
     trunc1 = truncrank(χ) & truncerror(; atol = 1.0e-12)
-    env0 = CTMRGEnv(randn, Float64, state, Vect[SU2Irrep](0 => 1))
+    env0 = CTMRGEnv(ones, Float64, state, Vect[SU2Irrep](0 => 1))
     env, = leading_boundary(env0, state; alg = :sequential, trunc = trunc1, tol = 1.0e-10)
     return env
 end
@@ -31,24 +28,27 @@ trunc_pepo = truncrank(7) & truncerror(; atol = 1.0e-12)
 check_interval = 100
 
 # PEPO approach
-dt, maxiter = 1.0e-3, 600
-alg = SimpleUpdate(dt, 0.0, maxiter, trunc_pepo)
-pepo, wts, = simpleupdate(pepo0, ham, alg, wts0; check_interval, gate_bothsides = true)
+dt, nstep = 1.0e-3, 600
+alg = SimpleUpdate(; trunc = trunc_pepo, purified = false)
+evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0)
+pepo, wts, info = time_evolve(evolver; check_interval)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 energy = expectation_value(pepo, ham, env) / (Nr * Nc)
-@info "β = $(dt * maxiter): tr(ρH) = $(energy)"
+@info "β = $(dt * nstep): tr(ρH) = $(energy)"
+@test dt * nstep ≈ info.t
 @test energy ≈ bm[2] atol = 5.0e-3
 
 # PEPS (purified PEPO) approach
-dt, maxiter = 1.0e-3, 300
-alg = SimpleUpdate(dt, 0.0, maxiter, trunc_pepo)
-pepo, wts, = simpleupdate(pepo0, ham, alg, wts0; check_interval, gate_bothsides = false)
+alg = SimpleUpdate(; trunc = trunc_pepo, purified = true)
+evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0)
+pepo, wts, info = time_evolve(evolver; check_interval)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 energy = expectation_value(pepo, ham, env) / (Nr * Nc)
-@info "β = $(dt * maxiter): tr(ρH) = $(energy)"
+@info "β = $(dt * nstep) / 2: tr(ρH) = $(energy)"
 @test energy ≈ bm[1] atol = 5.0e-3
 
 env = converge_env(InfinitePEPS(pepo), 16)
 energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc)
-@info "β = 2 × $(dt * maxiter): ⟨ρ|H|ρ⟩ = $(energy)"
+@info "β = $(dt * nstep): ⟨ρ|H|ρ⟩ = $(energy)"
+@test dt * nstep ≈ info.t
 @test energy ≈ bm[2] atol = 5.0e-3
diff --git a/test/examples/tf_ising_finiteT.jl b/test/examples/tf_ising_finiteT.jl
index d2cce0ea5..584b21bf0 100644
--- a/test/examples/tf_ising_finiteT.jl
+++ b/test/examples/tf_ising_finiteT.jl
@@ -1,12 +1,9 @@
 using Test
-using Random
 using LinearAlgebra
 using TensorKit
 import MPSKitModels: σˣ, σᶻ
 using PEPSKit
 
-Random.seed!(10235876)
-
 # Benchmark data of [σx, σz] from HOTRG
 # Physical Review B 86, 045139 (2012) Fig. 15-16
 bm_β = [0.5632, 0.0]
@@ -25,7 +22,7 @@ end
 
 function converge_env(state, χ::Int)
     trunc1 = truncrank(4) & truncerror(; atol = 1.0e-12)
-    env0 = CTMRGEnv(randn, Float64, state, ℂ^4)
+    env0 = CTMRGEnv(ones, Float64, state, ℂ^1)
     env, = leading_boundary(env0, state; alg = :sequential, trunc = trunc1, tol = 1.0e-10)
     trunc2 = truncrank(χ) & truncerror(; atol = 1.0e-12)
     env, = leading_boundary(env, state; alg = :sequential, trunc = trunc2, tol = 1.0e-10)
@@ -54,30 +51,33 @@ wts0 = SUWeight(pepo0)
 
 trunc_pepo = truncrank(8) & truncerror(; atol = 1.0e-12)
 
-dt, maxiter = 1.0e-3, 400
-β = dt * maxiter
-alg = SimpleUpdate(dt, 0.0, maxiter, trunc_pepo)
+dt, nstep = 1.0e-3, 400
+β = dt * nstep
 
 # when g = 2, β = 0.4 and 2β = 0.8 belong to two phases (without and with nonzero σᶻ)
 
-# PEPO approach
-## results at β, or T = 2.5
-pepo, wts, = simpleupdate(pepo0, ham, alg, wts0; gate_bothsides = true)
+# PEPO approach: results at β, or T = 2.5
+alg = SimpleUpdate(; trunc = trunc_pepo, purified = false)
+pepo, wts, info = time_evolve(pepo0, ham, dt, nstep, alg, wts0)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 result_β = measure_mag(pepo, env)
-@info "Magnetization at T = $(1 / β)" result_β
+@info "tr(σ(x,z)ρ) at T = $(1 / β)" result_β
+@test β ≈ info.t
 @test isapprox(abs.(result_β), bm_β, rtol = 1.0e-2)
 
-## results at 2β, or T = 1.25
-pepo, wts, = simpleupdate(pepo, ham, alg, wts; gate_bothsides = true)
+# continue to get results at 2β, or T = 1.25
+pepo, wts, info = time_evolve(pepo, ham, dt, nstep, alg, wts; t0 = β)
 env = converge_env(InfinitePartitionFunction(pepo), 16)
 result_2β = measure_mag(pepo, env)
-@info "Magnetization at T = $(1 / (2β))" result_2β
+@info "tr(σ(x,z)ρ) at T = $(1 / (2β))" result_2β
+@test 2 * β ≈ info.t
 @test isapprox(abs.(result_2β), bm_2β, rtol = 1.0e-4)
 
-# purification approach (should match 2β result)
-pepo, = simpleupdate(pepo0, ham, alg, wts0; gate_bothsides = false)
+# Purification approach: results at 2β, or T = 1.25
+alg = SimpleUpdate(; trunc = trunc_pepo, purified = true)
+pepo, wts, info = time_evolve(pepo0, ham, dt, 2 * nstep, alg, wts0)
 env = converge_env(InfinitePEPS(pepo), 8)
 result_2β′ = measure_mag(pepo, env; purified = true)
-@info "Magnetization at T = $(1 / (2β)) (purification approach)" result_2β′
+@info "⟨ρ|σ(x,z)|ρ⟩ at T = $(1 / (2β)) (purification approach)" result_2β′
+@test 2 * β ≈ info.t
 @test isapprox(abs.(result_2β′), bm_2β, rtol = 1.0e-2)
diff --git a/test/runtests.jl b/test/runtests.jl
index f18007942..38d4fdd4f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -61,6 +61,9 @@ end
         end
     end
     if GROUP == "ALL" || GROUP == "TIMEEVOL"
+        @time @safetestset "`timestep` function" begin
+            include("timeevol/timestep.jl")
+        end
         @time @safetestset "Cluster truncation with projectors" begin
             include("timeevol/cluster_projectors.jl")
         end
diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl
index 36359e1a4..279907c16 100644
--- a/test/timeevol/cluster_projectors.jl
+++ b/test/timeevol/cluster_projectors.jl
@@ -117,11 +117,11 @@ end
     # usual 2-site simple update, and measure energy
     dts = [1.0e-2, 1.0e-2, 5.0e-3]
     tols = [1.0e-8, 1.0e-8, 1.0e-8]
-    maxiter = 5000
+    nstep = 5000
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
         trunc = truncerror(; atol = 1.0e-10) & truncrank(n == 1 ? 4 : 2)
-        alg = SimpleUpdate(dt, tol, maxiter, trunc)
-        peps, wts, = simpleupdate(peps, ham, alg, wts; bipartite = true, check_interval = 1000)
+        alg = SimpleUpdate(; trunc, bipartite = true)
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; tol, check_interval = 1000)
     end
     normalize!.(peps.A, Inf)
     env = CTMRGEnv(wts, peps)
@@ -132,11 +132,9 @@ end
     dts = [1.0e-2, 5.0e-3]
     tols = [1.0e-8, 1.0e-8]
     trunc = truncerror(; atol = 1.0e-10) & truncrank(2)
+    alg = SimpleUpdate(; trunc, force_3site = true)
     for (n, (dt, tol)) in enumerate(zip(dts, tols))
-        alg = SimpleUpdate(dt, tol, maxiter, trunc)
-        peps, wts, = simpleupdate(
-            peps, ham, alg, wts; check_interval = 1000, force_3site = true
-        )
+        peps, wts, = time_evolve(peps, ham, dt, nstep, alg, wts; tol, check_interval = 1000)
     end
     normalize!.(peps.A, Inf)
     env, = leading_boundary(env, peps; tol = ctmrg_tol, trunc = trunc_env)
diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl
index 6f067592d..fbc86fb86 100644
--- a/test/timeevol/sitedep_truncation.jl
+++ b/test/timeevol/sitedep_truncation.jl
@@ -27,8 +27,8 @@ end
     # set trunc to be compatible with bipartite structure
     bonddims = stack([[6 4; 4 6], [5 7; 7 5]]; dims = 1)
     trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
-    alg = SimpleUpdate(1.0e-2, 1.0e-14, 4, trunc)
-    peps, env, = simpleupdate(peps0, ham, alg, env0; bipartite = true)
+    alg = SimpleUpdate(; trunc, bipartite = true)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # check bipartite structure is preserved
@@ -53,13 +53,14 @@ end
     bonddims = rand(2:8, 2, Nr, Nc)
     @show bonddims
     trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims))
-    alg = SimpleUpdate(1.0e-2, 1.0e-14, 2, trunc)
     # 2-site SU
-    peps, env, = simpleupdate(peps0, ham, alg, env0; bipartite = false)
+    alg = SimpleUpdate(; trunc)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
     # 3-site SU
-    peps, env, = simpleupdate(peps0, ham, alg, env0; bipartite = false, force_3site = true)
+    alg = SimpleUpdate(; trunc, force_3site = true)
+    peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0)
     @test get_bonddims(peps) == bonddims
     @test get_bonddims(env) == bonddims
 end
diff --git a/test/timeevol/timestep.jl b/test/timeevol/timestep.jl
new file mode 100644
index 000000000..39153a898
--- /dev/null
+++ b/test/timeevol/timestep.jl
@@ -0,0 +1,33 @@
+using Test
+using Random
+using TensorKit
+using PEPSKit
+
+@testset "SimpleUpdate timestep" begin
+    Nr, Nc = 2, 2
+    H = real(heisenberg_XYZ(ComplexF64, Trivial, InfiniteSquare(Nr, Nc); Jx = 1, Jy = 1, Jz = 1))
+    Pspace, Vspace = ℂ^2, ℂ^4
+    ψ0 = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc))
+    env0 = SUWeight(ψ0)
+    alg = SimpleUpdate(; trunc = truncerror(; atol = 1.0e-10) & truncrank(4))
+    dt, nstep = 1.0e-2, 50
+    # manual timestep
+    evolver = TimeEvolver(ψ0, H, dt, nstep, alg, env0)
+    ψ1, env1, info1 = deepcopy(ψ0), deepcopy(env0), nothing
+    for iter in 0:(nstep - 1)
+        ψ1, env1, info1 = timestep(evolver, ψ1, env1)
+    end
+    # time_evolve
+    ψ2, env2, info2 = time_evolve(ψ0, H, dt, nstep, alg, env0)
+    # for-loop syntax
+    ## manually reset internal state of evolver
+    evolver.state = PEPSKit.SUState(0, 0.0, ψ0, env0)
+    ψ3, env3, info3 = nothing, nothing, nothing
+    for state in evolver
+        ψ3, env3, info3 = state
+    end
+    # results should be *exactly* the same
+    @test ψ1 == ψ2 == ψ3
+    @test env1 == env2 == env3
+    @test info1 == info2 == info3
+end