Add MPSKit.transfer_spectrum method instead, fix correlation_length

src/algorithms/ctmrg.jl

@@ -336,25 +336,32 @@ Compute the PEPS correlation length based on the horizontal and vertical
 transfer matrices. Additionally the (normalized) eigenvalue spectrum is
 returned. Specify the number of computed eigenvalues with `howmany`.
 """
-function MPSKit.correlation_length(peps::InfinitePEPS, env::CTMRGEnv; howmany=2)
-    ξ_h = Vector{Float64}(undef, size(peps, 1))
-    ξ_v = Vector{Float64}(undef, size(peps, 2))
-    λ_h = Matrix{ComplexF64}(undef, size(peps, 1), howmany)
-    λ_v = Matrix{ComplexF64}(undef, size(peps, 2), howmany)
+function MPSKit.correlation_length(peps::InfinitePEPS, env::CTMRGEnv; num_vals=2)
+    T = scalartype(peps)
+    ξ_h = Vector{real(T)}(undef, size(peps, 1))
+    ξ_v = Vector{real(T)}(undef, size(peps, 2))
+    λ_h = Vector{Vector{T}}(undef, size(peps, 1))
+    λ_v = Vector{Vector{T}}(undef, size(peps, 2))
 
     # Horizontal
     above_h = MPSMultiline(map(r -> InfiniteMPS(env.edges[1, r, :]), 1:size(peps, 1)))
-    below_h = MPSMultiline(map(r -> InfiniteMPS(env.edges[3, r, :]), 1:size(peps, 1)))
+    respaced_edges_h = map(zip(space.(env.edges)[1, :, :], env.edges[3, :, :])) do (V1, T3)
+        return TensorMap(T3.data, V1)
+    end
+    below_h = MPSMultiline(map(r -> InfiniteMPS(respaced_edges_h[r, :]), 1:size(peps, 1)))
     transfer_peps_h = TransferPEPSMultiline(peps, NORTH)
-    _, (_, vals_h) = MPSKit.mixed_fixpoints(above_h, transfer_peps_h, below_h; howmany)
+    vals_h = MPSKit.transfer_spectrum(above_h, transfer_peps_h, below_h; num_vals)
     λ_h = map(λ_row -> λ_row / abs(λ_row[1]), vals_h)  # Normalize largest eigenvalue
     ξ_h = map(λ_row -> -1 / log(abs(λ_row[2])), λ_h)
 
     # Vertical
     above_v = MPSMultiline(map(c -> InfiniteMPS(env.edges[2, :, c]), 1:size(peps, 2)))
-    below_v = MPSMultiline(map(c -> InfiniteMPS(env.edges[4, :, c]), 1:size(peps, 2)))
+    respaced_edges_v = map(zip(space.(env.edges)[2, :, :], env.edges[4, :, :])) do (V2, T4)
+        return TensorMap(T4.data, V2)
+    end
+    below_v = MPSMultiline(map(c -> InfiniteMPS(respaced_edges_v[:, c]), 1:size(peps, 2)))
     transfer_peps_v = TransferPEPSMultiline(peps, EAST)
-    _, (_, vals_v) = MPSKit.mixed_fixpoints(above_v, transfer_peps_v, below_v; howmany)
+    vals_v = MPSKit.transfer_spectrum(above_v, transfer_peps_v, below_v; num_vals)
     λ_v = map(λ_row -> λ_row / abs(λ_row[1]), vals_v)  # Normalize largest eigenvalue
     ξ_v = map(λ_row -> -1 / log(abs(λ_row[2])), λ_v)
 

src/mpskit_glue/transferpeps_environments.jl

@@ -20,7 +20,6 @@ function MPSKit.mixed_fixpoints(
     below::MPSMultiline,
     init=gen_init_fps(above, O, below);
     solver=MPSKit.Defaults.eigsolver,
-    howmany=1
 )
     T = eltype(above)
 
@@ -31,8 +30,6 @@ function MPSKit.mixed_fixpoints(
     envtype = eltype(init[1])
     lefties = PeriodicArray{envtype,2}(undef, numrows, numcols)
     righties = PeriodicArray{envtype,2}(undef, numrows, numcols)
-    lvals = Vector{Vector{scalartype(envtype)}}(undef, numrows)
-    rvals = Vector{Vector{scalartype(envtype)}}(undef, numrows)
 
     @threads for cr in 1:numrows
         c_above = above[cr]
@@ -43,20 +40,18 @@ function MPSKit.mixed_fixpoints(
         @sync begin
             Threads.@spawn begin
                 E_LL = TransferMatrix($c_above.AL, $O[cr], $c_below.AL)
-                (Lvals, Ls, convhist) = eigsolve(flip(E_LL), $L0, howmany, :LM, $solver)
-                convhist.converged < howmany &&
+                (_, Ls, convhist) = eigsolve(flip(E_LL), $L0, 1, :LM, $solver)
+                convhist.converged < 1 &&
                     @info "left eigenvalue failed to converge $(convhist.normres)"
                 L0 = first(Ls)
-                lvals[cr] = Lvals[1:howmany]
             end
 
             Threads.@spawn begin
                 E_RR = TransferMatrix($c_above.AR, $O[cr], $c_below.AR)
-                (Rvals, Rs, convhist) = eigsolve(E_RR, $R0, howmany, :LM, $solver)
-                convhist.converged < howmany &&
+                (_, Rs, convhist) = eigsolve(E_RR, $R0, 1, :LM, $solver)
+                convhist.converged < 1 &&
                     @info "right eigenvalue failed to converge $(convhist.normres)"
                 R0 = first(Rs)
-                rvals[cr] = Rvals[1:howmany]
             end
         end
 
@@ -86,7 +81,7 @@ function MPSKit.mixed_fixpoints(
         end
     end
 
-    return (lefties, righties), (lvals, rvals)
+    return (lefties, righties)
 end
 
 function gen_init_fps(above::MPSMultiline, O::TransferPEPSMultiline, below::MPSMultiline)
@@ -112,3 +107,31 @@ function gen_init_fps(above::MPSMultiline, O::TransferPEPSMultiline, below::MPSM
         (L0, R0)
     end
 end
+
+function MPSKit.transfer_spectrum(
+    above::MPSMultiline,
+    O::TransferPEPSMultiline,
+    below::MPSMultiline,
+    init=gen_init_fps(above, O, below);
+    num_vals=2,
+    solver=MPSKit.Defaults.eigsolver,
+)
+    @assert size(above) == size(O)
+    @assert size(below) == size(O)
+
+    numrows = size(above, 2)
+    envtype = eltype(init[1])
+    eigenvals = Vector{Vector{scalartype(envtype)}}(undef, numrows)
+
+    @threads for cr in 1:numrows
+        L0, = init[cr]
+
+        E_LL = TransferMatrix(above[cr].AL, O[cr + 1], below[cr + 2].AL)
+        λ, _, convhist = eigsolve(flip(E_LL), L0, num_vals, :LM, solver)
+        convhist.converged < num_vals &&
+            @warn "correlation length failed to converge: normres = $(convhist.normres)"
+        eigenvals[cr] = λ
+    end
+
+    return eigenvals
+end