diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl index 48754751..2a924a14 100644 --- a/src/algorithms/groundstate/dmrg.jl +++ b/src/algorithms/groundstate/dmrg.jl @@ -20,30 +20,33 @@ Single site DMRG algorithm for finding groundstates. end function find_groundstate!(Ψ::AbstractFiniteMPS, H, alg::DMRG, envs=environments(Ψ, H)) - tol = alg.tol - maxiter = alg.maxiter - iter = 0 - delta::Float64 = 2 * tol - - while iter < maxiter && delta > tol - delta = 0.0 - - for pos in [1:(length(Ψ) - 1); length(Ψ):-1:2] - h = ∂∂AC(pos, Ψ, H, envs) - _, vecs = eigsolve(h, Ψ.AC[pos], 1, :SR, alg.eigalg) - delta = max(delta, calc_galerkin(Ψ, pos, envs)) - Ψ.AC[pos] = vecs[1] + t₀ = Base.time_ns() + ϵ::Float64 = 2 * alg.tol + + for iter in 1:(alg.maxiter) + global ϵ = 0.0 + Δt = @elapsed begin + for pos in [1:(length(Ψ) - 1); length(Ψ):-1:2] + h = ∂∂AC(pos, Ψ, H, envs) + _, vecs = eigsolve(h, Ψ.AC[pos], 1, :SR, alg.eigalg) + ϵ = max(ϵ, calc_galerkin(Ψ, pos, envs)) + Ψ.AC[pos] = vecs[1] + end + + Ψ, envs = alg.finalize(iter, Ψ, H, envs)::Tuple{typeof(Ψ),typeof(envs)} end - alg.verbose && @info "Iteraton $(iter) error $(delta)" - flush(stdout) + alg.verbose && + @info "DMRG iteration:" iter ϵ λ = sum(expectation_value(Ψ, H, envs)) Δt - iter += 1 - - Ψ, envs = alg.finalize(iter, Ψ, H, envs)::Tuple{typeof(Ψ),typeof(envs)} + ϵ <= alg.tol && break + iter == alg.maxiter && + @warn "DMRG maximum iterations" iter ϵ λ = sum(expectation_value(Ψ, H, envs)) end - return Ψ, envs, delta + Δt = (Base.time_ns() - t₀) / 1.0e9 + alg.verbose && @info "DMRG summary:" ϵ λ = sum(expectation_value(Ψ, H, envs)) Δt + return Ψ, envs, ϵ end """ @@ -70,53 +73,62 @@ end end function find_groundstate!(Ψ::AbstractFiniteMPS, H, alg::DMRG2, envs=environments(Ψ, H)) - tol = alg.tol - maxiter = alg.maxiter - iter = 0 - delta::Float64 = 2 * tol - - while iter < maxiter && delta > tol - delta = 0.0 - - #left to right sweep - for pos in 1:(length(Ψ) - 1) - @plansor ac2[-1 -2; -3 -4] := Ψ.AC[pos][-1 -2; 1] * Ψ.AR[pos + 1][1 -4; -3] - - _, vecs = eigsolve(∂∂AC2(pos, Ψ, H, envs), ac2, 1, :SR, alg.eigalg) - newA2center = first(vecs) - - al, c, ar, ϵ = tsvd(newA2center; trunc=alg.trscheme, alg=TensorKit.SVD()) - normalize!(c) - v = @plansor ac2[1 2; 3 4] * conj(al[1 2; 5]) * conj(c[5; 6]) * conj(ar[6; 3 4]) - delta = max(delta, abs(1 - abs(v))) - - Ψ.AC[pos] = (al, complex(c)) - Ψ.AC[pos + 1] = (complex(c), _transpose_front(ar)) + t₀ = Base.time_ns() + ϵ::Float64 = 2 * alg.tol + + for iter in 1:(alg.maxiter) + ϵ = 0.0 + Δt = @elapsed begin + #left to right sweep + for pos in 1:(length(Ψ) - 1) + @plansor ac2[-1 -2; -3 -4] := Ψ.AC[pos][-1 -2; 1] * Ψ.AR[pos + 1][1 -4; -3] + + _, vecs = eigsolve(∂∂AC2(pos, Ψ, H, envs), ac2, 1, :SR, alg.eigalg) + newA2center = first(vecs) + + al, c, ar, = tsvd!(newA2center; trunc=alg.trscheme, alg=TensorKit.SVD()) + normalize!(c) + v = @plansor ac2[1 2; 3 4] * + conj(al[1 2; 5]) * + conj(c[5; 6]) * + conj(ar[6; 3 4]) + ϵ = max(ϵ, abs(1 - abs(v))) + + Ψ.AC[pos] = (al, complex(c)) + Ψ.AC[pos + 1] = (complex(c), _transpose_front(ar)) + end + + for pos in (length(Ψ) - 2):-1:1 + @plansor ac2[-1 -2; -3 -4] := Ψ.AL[pos][-1 -2; 1] * Ψ.AC[pos + 1][1 -4; -3] + + _, vecs = eigsolve(∂∂AC2(pos, Ψ, H, envs), ac2, 1, :SR, alg.eigalg) + newA2center = first(vecs) + + al, c, ar, = tsvd!(newA2center; trunc=alg.trscheme, alg=TensorKit.SVD()) + normalize!(c) + v = @plansor ac2[1 2; 3 4] * + conj(al[1 2; 5]) * + conj(c[5; 6]) * + conj(ar[6; 3 4]) + ϵ = max(ϵ, abs(1 - abs(v))) + + Ψ.AC[pos + 1] = (complex(c), _transpose_front(ar)) + Ψ.AC[pos] = (al, complex(c)) + end + + Ψ, envs = alg.finalize(iter, Ψ, H, envs)::Tuple{typeof(Ψ),typeof(envs)} end - for pos in (length(Ψ) - 2):-1:1 - @plansor ac2[-1 -2; -3 -4] := Ψ.AL[pos][-1 -2; 1] * Ψ.AC[pos + 1][1 -4; -3] - - _, vecs = eigsolve(∂∂AC2(pos, Ψ, H, envs), ac2, 1, :SR, alg.eigalg) - newA2center = first(vecs) - - al, c, ar, ϵ = tsvd(newA2center; trunc=alg.trscheme, alg=TensorKit.SVD()) - normalize!(c) - v = @plansor ac2[1 2; 3 4] * conj(al[1 2; 5]) * conj(c[5; 6]) * conj(ar[6; 3 4]) - delta = max(delta, abs(1 - abs(v))) - - Ψ.AC[pos + 1] = (complex(c), _transpose_front(ar)) - Ψ.AC[pos] = (al, complex(c)) - end + alg.verbose && @info "DMRG iteration:" iter ϵ λ = expectation_value(Ψ, H, envs) Δt - alg.verbose && @info "Iteraton $(iter) error $(delta)" - flush(stdout) - #finalize - Ψ, envs = alg.finalize(iter, Ψ, H, envs)::Tuple{typeof(Ψ),typeof(envs)} - iter += 1 + ϵ <= alg.tol && break + iter == alg.maxiter && + @warn "DMRG maximum iterations" iter ϵ λ = expectation_value(Ψ, H, envs) end - return Ψ, envs, delta + Δt = (Base.time_ns() - t₀) / 1.0e9 + alg.verbose && @info "DMRG summary:" ϵ λ = sum(expectation_value(Ψ, H, envs)) Δt + return Ψ, envs, ϵ end function find_groundstate(Ψ, H, alg::Union{DMRG,DMRG2}, envs...) diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl index 2af61d88..9914102b 100644 --- a/src/algorithms/groundstate/vumps.jl +++ b/src/algorithms/groundstate/vumps.jl @@ -35,11 +35,11 @@ https://arxiv.org/abs/1701.07035. gauge_tolfactor::Float64 = Defaults.gauge_tolfactor end -function updatetols(alg::VUMPS, iter, ε) +function updatetols(alg::VUMPS, iter, ϵ) if alg.dynamical_tols - tol_eigs = between(alg.tol_min, ε * alg.eigs_tolfactor / sqrt(iter), alg.tol_max) - tol_envs = between(alg.tol_min, ε * alg.envs_tolfactor / sqrt(iter), alg.tol_max) - tol_gauge = between(alg.tol_min, ε * alg.gauge_tolfactor / sqrt(iter), alg.tol_max) + tol_eigs = between(alg.tol_min, ϵ * alg.eigs_tolfactor / sqrt(iter), alg.tol_max) + tol_envs = between(alg.tol_min, ϵ * alg.envs_tolfactor / sqrt(iter), alg.tol_max) + tol_gauge = between(alg.tol_min, ϵ * alg.gauge_tolfactor / sqrt(iter), alg.tol_max) else # preserve legacy behavior tol_eigs = alg.tol_galerkin / 10 tol_envs = Defaults.tol @@ -56,11 +56,11 @@ end function find_groundstate(Ψ::InfiniteMPS, H, alg::VUMPS, envs=environments(Ψ, H)) t₀ = Base.time_ns() - ε::Float64 = calc_galerkin(Ψ, envs) + ϵ::Float64 = calc_galerkin(Ψ, envs) temp_ACs = similar.(Ψ.AC) for iter in 1:(alg.maxiter) - tol_eigs, tol_envs, tol_gauge = updatetols(alg, iter, ε) + tol_eigs, tol_envs, tol_gauge = updatetols(alg, iter, ϵ) Δt = @elapsed begin eigalg = Arnoldi(; tol=tol_eigs) @@ -83,21 +83,20 @@ function find_groundstate(Ψ::InfiniteMPS, H, alg::VUMPS, envs=environments(Ψ, Ψ, envs = alg.finalize(iter, Ψ, H, envs)::Tuple{typeof(Ψ),typeof(envs)} - ε = calc_galerkin(Ψ, envs) + ϵ = calc_galerkin(Ψ, envs) end alg.verbose && - @info "VUMPS iteration:" iter ε λ = sum(expectation_value(Ψ, H, envs)) Δt + @info "VUMPS iteration:" iter ϵ λ = sum(expectation_value(Ψ, H, envs)) Δt - ε <= alg.tol_galerkin && break + ϵ <= alg.tol_galerkin && break iter == alg.maxiter && - @warn "VUMPS maximum iterations" iter ε λ = sum(expectation_value(Ψ, H, envs)) Δt + @warn "VUMPS maximum iterations" iter ϵ λ = sum(expectation_value(Ψ, H, envs)) end - alg.verbose && @info "VUMPS summary:" ε λ = sum(expectation_value(Ψ, H, envs)) Δt = ( - (Base.time_ns() - t₀) / 1.0e9 - ) - return Ψ, envs, ε + Δt = (Base.time_ns() - t₀) / 1.0e9 + alg.verbose && @info "VUMPS summary:" ϵ λ = sum(expectation_value(Ψ, H, envs)) Δt + return Ψ, envs, ϵ end function _vumps_localupdate!(AC′, loc, Ψ, H, envs, eigalg, factalg=QRpos())