Combine pr

src/algorithms/changebonds/vumpssvd.jl

@@ -77,10 +77,68 @@
     return state, envs
 end
 
+function changebonds_parallel_n(state::InfiniteMPS, H, alg::VUMPSSvdCut,
+                                envs=environments(state, H))
+    #prepare the loops we'll have to run over, this is different depending on wether the length of the state is odd or even ! 
+    subsets = nothing
+    start_locs = nothing
+    if iseven(length(state))
+        subsets = ["even", "odd"]
+        start_locs = Dict("even" => 1:2:length(state), "odd" => 2:2:length(state))
+    else
+        subsets = ["even", "odd", "last"]
+        start_locs = Dict("even" => 1:2:(length(state) - 1),
+                          "odd" => 2:2:(length(state) - 1), "last" => length(state))
+    end
+
+    for subset in subsets
+        new_ALs = copy(state.AL)
+        @sync for loc in start_locs[subset]
+            Threads.@spawn begin
+                @plansor AC2[-1 -2; -3 -4] := state.AC[loc][-1 -2; 1] *
+                                              state.AR[loc + 1][1 -4; -3]
+
+                h_ac2 = ∂∂AC2(loc, state, H, envs)
+                (vals, vecs, _) = eigsolve(h_ac2, AC2, 1, :SR; tol=alg.tol_eigenval,
+                                           ishermitian=false)
+                nAC2 = vecs[1]
+
+                h_c = ∂∂C(loc + 1, state, H, envs)
+                (vals, vecs, _) = eigsolve(h_c, state.CR[loc + 1], 1, :SR;
+                                           tol=alg.tol_eigenval,
+                                           ishermitian=false)
+                nC2 = vecs[1]
+
+                #svd ac2, get new AL1 and S,V ---> AC
+                (AL1, S, V, eps) = tsvd(nAC2; trunc=alg.trscheme, alg=TensorKit.SVD())
+                @plansor AC[-1 -2; -3] := S[-1; 1] * V[1; -3 -2]
+
+                #find AL2 from AC and C as in vumps paper
+                QAC, _ = leftorth(AC; alg=QRpos())
+                QC, _ = leftorth(nC2; alg=QRpos())
+                dom_map = isometry(domain(QC), domain(QAC))
+
+                @plansor AL2[-1 -2; -3] := QAC[-1 -2; 1] * conj(dom_map[2; 1]) *
+                                           conj(QC[-3; 2])
+
+                #make a new state using the updated A's
+                new_ALs[loc] = AL1
+                new_ALs[loc + 1] = AL2
+            end
+        end
+        state = InfiniteMPS(new_ALs; tol=alg.tol_gauge)
+    end
+    return state, envs
+end
+
 function changebonds(state::InfiniteMPS, H, alg::VUMPSSvdCut, envs=environments(state, H))
     if (length(state) == 1)
         return changebonds_1(state, H, alg, envs)
     else
-        return changebonds_n(state, H, alg, envs)
+        @static if Defaults.parallelize_sites
+            return changebonds_parallel_n(state, H, alg, envs)
+        else
+            return changebonds_n(state, H, alg, envs)
+        end
     end
 end

src/algorithms/grassmann.jl

@@ -11,6 +11,7 @@
 module GrassmannMPS
 
 using ..MPSKit
+using ..Defaults
 using TensorKit
 import TensorKitManifolds.Grassmann
 
@@ -69,22 +70,38 @@
 
 function ManifoldPoint(state::Union{InfiniteMPS,FiniteMPS}, envs)
     al_d = similar(state.AL)
-    for i in 1:length(state)
-        al_d[i] = MPSKit.∂∂AC(i, state, envs.opp, envs) * state.AC[i]
+    @static if Defaults.parallelize_sites
+        @sync for i in 1:length(state)
+            Threads.@spawn begin
+                al_d[i] = MPSKit.∂∂AC(i, state, envs.opp, envs) * state.AC[i]
+            end
+        end
+    else
+        for i in 1:length(state)
+            al_d[i] = MPSKit.∂∂AC(i, state, envs.opp, envs) * state.AC[i]
+        end
     end
-
     g = Grassmann.project.(al_d, state.AL)
 
     Rhoreg = Vector{eltype(state.CR)}(undef, length(state))
     δmin = sqrt(eps(real(scalartype(state))))
-    for i in 1:length(state)
-        Rhoreg[i] = regularize(state.CR[i], max(norm(g[i]) / 10, δmin))
+    @static if Defaults.parallelize_sites
+        @sync for i in 1:length(state)
+            Threads.@spawn begin
+                Rhoreg[i] = regularize(state.CR[i], max(norm(g[i]) / 10, δmin))
+            end
+        end
+    else
+        for i in 1:length(state)
+            Rhoreg[i] = regularize(state.CR[i], max(norm(g[i]) / 10, δmin))
+        end
     end
 
     return ManifoldPoint(state, envs, g, Rhoreg)
 end
 
 function ManifoldPoint(state::MPSMultiline, envs)
+    #TODO : support parralelize_sites
     # FIXME: add support for unitcells
     @assert length(state.AL) == 1 "GradientGrassmann only supports MPSMultiline without unitcells for now"
 
@@ -115,9 +132,16 @@
 function fg(x::ManifoldPoint{T}) where {T<:Union{InfiniteMPS,FiniteMPS}}
     # the gradient I want to return is the preconditioned gradient!
     g_prec = Vector{PrecGrad{eltype(x.g),eltype(x.Rhoreg)}}(undef, length(x.g))
-
-    for i in 1:length(x.state)
-        g_prec[i] = PrecGrad(rmul!(copy(x.g[i]), x.state.CR[i]'), x.Rhoreg[i])
+    @static if Defaults.parallelize_sites
+        @sync for i in 1:length(x.state)
+            Threads.@spawn begin
+                g_prec[i] = PrecGrad(rmul!(copy(x.g[i]), x.state.CR[i]'), x.Rhoreg[i])
+            end
+        end
+    else
+        for i in 1:length(x.state)
+            g_prec[i] = PrecGrad(rmul!(copy(x.g[i]), x.state.CR[i]'), x.Rhoreg[i])
+        end
     end
 
     # TODO: the operator really should not be part of the environments, and this should
@@ -128,6 +152,7 @@
     return real(f), g_prec
 end
 function fg(x::ManifoldPoint{<:MPSMultiline})
+    #TODO : support parralelize_sites
     @assert length(x.state) == 1 "GradientGrassmann only supports MPSMultiline without unitcells for now"
     # the gradient I want to return is the preconditioned gradient!
     g_prec = map(enumerate(x.g)) do (i, cg)
@@ -147,6 +172,7 @@
 Retract a left-canonical MPSMultiline along Grassmann tangent `g` by distance `alpha`.
 """
 function retract(x::ManifoldPoint{<:MPSMultiline}, tg, alpha)
+    #TODO : support parralelize_sites
     g = reshape(tg, size(x.state))
 
     nal = similar(x.state.AL)
@@ -170,11 +196,19 @@
     envs = x.envs
     nal = similar(state.AL)
     h = similar(g)  # The tangent at the end-point
-    for i in 1:length(g)
-        nal[i], th = Grassmann.retract(state.AL[i], g[i].Pg, alpha)
-        h[i] = PrecGrad(th)
+    @static if Defaults.parallelize_sites
+        @sync for i in 1:length(g)
+            Threads.@spawn begin
+                nal[i], th = Grassmann.retract(state.AL[i], g[i].Pg, alpha)
+                h[i] = PrecGrad(th)
+            end
+        end
+    else
+        for i in 1:length(g)
+            nal[i], th = Grassmann.retract(state.AL[i], g[i].Pg, alpha)
+            h[i] = PrecGrad(th)
+        end
     end
-
     nstate = InfiniteMPS(nal, state.CR[end])
 
     newpoint = ManifoldPoint(nstate, envs)
@@ -186,6 +220,7 @@
 Retract a left-canonical finite MPS along Grassmann tangent `g` by distance `alpha`.
 """
 function retract(x::ManifoldPoint{<:FiniteMPS}, g, alpha)
+    #TODO : support parralelize_sites.
     state = x.state
     envs = x.envs
 
@@ -208,9 +243,18 @@
 `alpha`. `xp` is the end-point of the retraction.
 """
 function transport!(h, x, g, alpha, xp)
-    for i in 1:length(h)
-        h[i] = PrecGrad(Grassmann.transport!(h[i].Pg, x.state.AL[i], g[i].Pg, alpha,
-                                             xp.state.AL[i]))
+    @static if Defaults.parallelize_sites
+        @sync for i in 1:length(h)
+            Threads.@spawn begin
+                h[i] = PrecGrad(Grassmann.transport!(h[i].Pg, x.state.AL[i], g[i].Pg, alpha,
+                                                     xp.state.AL[i]))
+            end
+        end
+    else
+        for i in 1:length(h)
+            h[i] = PrecGrad(Grassmann.transport!(h[i].Pg, x.state.AL[i], g[i].Pg, alpha,
+                                                 xp.state.AL[i]))
+        end
     end
     return h
 end

test/algorithms.jl

@@ -161,6 +161,98 @@ end
     end
 end
 
+@testset "InfiniteMPS 3-site groundstates" verbose = true begin
+    tol = 1e-8
+    g = 4.0
+    D = 6
+    H = repeat(force_planar(transverse_field_ising(; g)), 3)
+
+    @testset "VUMPS" begin
+        ψ₀ = repeat(InfiniteMPS(ℙ^2, ℙ^D), 3)
+        v₀ = variance(ψ₀, H)
+
+        # test logging
+        ψ, envs, δ = find_groundstate(ψ₀, H, VUMPS(; tol, verbosity=5, maxiter=2))
+
+        ψ, envs, δ = find_groundstate(ψ, H, VUMPS(; tol, verbosity=1))
+        v = variance(ψ, H, envs)
+
+        # test using low variance
+        @test sum(δ) ≈ 0 atol = 1e-3
+        @test v < v₀
+        @test v < 1e-2
+    end
+
+    @testset "IDMRG1" begin
+        ψ₀ = repeat(InfiniteMPS(ℙ^2, ℙ^D), 3)
+        v₀ = variance(ψ₀, H)
+
+        # test logging
+        ψ, envs, δ = find_groundstate(ψ₀, H, IDMRG1(; tol, verbosity=5, maxiter=2))
+
+        ψ, envs, δ = find_groundstate(ψ, H, IDMRG1(; tol, verbosity=1))
+        v = variance(ψ, H, envs)
+
+        # test using low variance
+        @test sum(δ) ≈ 0 atol = 1e-3
+        @test v < v₀
+        @test v < 1e-2
+    end
+
+    @testset "IDMRG2" begin
+        ψ₀ = repeat(InfiniteMPS(ℙ^2, ℙ^D), 3)
+        v₀ = variance(ψ₀, H)
+        trscheme = truncbelow(1e-8)
+
+        # test logging
+        ψ, envs, δ = find_groundstate(ψ₀, H,
+                                      IDMRG2(; tol, verbosity=5, maxiter=2,
+                                             trscheme))
+
+        ψ, envs, δ = find_groundstate(ψ, H,
+                                      IDMRG2(; tol, verbosity=1, trscheme))
+        v = variance(ψ, H, envs)
+
+        # test using low variance
+        @test sum(δ) ≈ 0 atol = 1e-3
+        @test v < v₀
+        @test v < 1e-2
+    end
+
+    @testset "GradientGrassmann" begin
+        ψ₀ = repeat(InfiniteMPS(ℙ^2, ℙ^D), 3)
+        v₀ = variance(ψ₀, H)
+
+        # test logging
+        ψ, envs, δ = find_groundstate(ψ₀, H,
+                                      GradientGrassmann(; tol, verbosity=5, maxiter=2))
+
+        ψ, envs, δ = find_groundstate(ψ, H, GradientGrassmann(; tol, verbosity=1))
+        v = variance(ψ, H, envs)
+
+        # test using low variance
+        @test sum(δ) ≈ 0 atol = 1e-3
+        @test v < v₀
+        @test v < 1e-2
+    end
+
+    @testset "Combination" begin
+        ψ₀ = repeat(InfiniteMPS(ℙ^2, ℙ^D), 3)
+        v₀ = variance(ψ₀, H)
+
+        alg = VUMPS(; tol=100 * tol, verbosity=1, maxiter=10) &
+              GradientGrassmann(; tol, verbosity=1, maxiter=50)
+        ψ, envs, δ = find_groundstate(ψ₀, H, alg)
+
+        v = variance(ψ, H, envs)
+
+        # test using low variance
+        @test sum(δ) ≈ 0 atol = 1e-3
+        @test v < v₀
+        @test v < 1e-2
+    end
+end
+
 @testset "LazySum FiniteMPS groundstate" verbose = true begin
     tol = 1e-8
     D = 15
@@ -471,25 +563,32 @@ end
         nn = TensorMap(rand, ComplexF64, pspace * pspace, pspace * pspace)
         nn += nn'
 
-        state = InfiniteMPS([pspace, pspace], [Dspace, Dspace])
+        state1 = InfiniteMPS(pspace, Dspace)
+        state2 = InfiniteMPS([pspace, pspace], [Dspace, Dspace])
+        state3 = InfiniteMPS([pspace, pspace, pspace], [Dspace, Dspace, Dspace])
 
-        state_re = changebonds(state,
+        state_re = changebonds(state2,
                                RandExpand(; trscheme=truncdim(dim(Dspace) * dim(Dspace))))
-        @test dot(state, state_re) ≈ 1 atol = 1e-8
+        @test dot(state2, state_re) ≈ 1 atol = 1e-8
 
-        state_oe, _ = changebonds(state,
+        state_oe, _ = changebonds(state2,
                                   repeat(MPOHamiltonian(nn), 2),
                                   OptimalExpand(;
                                                 trscheme=truncdim(dim(Dspace) *
                                                                   dim(Dspace))))
-        @test dot(state, state_oe) ≈ 1 atol = 1e-8
-
-        state_vs, _ = changebonds(state, repeat(MPOHamiltonian(nn), 2),
-                                  VUMPSSvdCut(; trscheme=notrunc()))
-        @test dim(left_virtualspace(state, 1)) < dim(left_virtualspace(state_vs, 1))
-
-        state_vs_tr = changebonds(state_vs, SvdCut(; trscheme=truncdim(dim(Dspace))))
-        @test dim(right_virtualspace(state_vs_tr, 1)) < dim(right_virtualspace(state_vs, 1))
+        @test dot(state2, state_oe) ≈ 1 atol = 1e-8
+
+        #test vumpssvd with different unit cells 
+        for state in [state1, state2, state3]
+            @show length(state)
+            state_vs, _ = changebonds(state, repeat(MPOHamiltonian(nn), length(state)),
+                                      VUMPSSvdCut(; trscheme=notrunc()))
+            @test dim(left_virtualspace(state, 1)) < dim(left_virtualspace(state_vs, 1))
+
+            state_vs_tr = changebonds(state_vs, SvdCut(; trscheme=truncdim(dim(Dspace))))
+            @test dim(right_virtualspace(state_vs_tr, 1)) <
+                  dim(right_virtualspace(state_vs, 1))
+        end
     end
 
     @testset "finite mps" begin