Add (fast) full update

src/PEPSKit.jl

@@ -77,6 +77,7 @@ 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")
+include("algorithms/time_evolution/fullupdate.jl")
 
 include("algorithms/transfermatrix.jl")
 include("algorithms/toolbox.jl")
@@ -105,7 +106,7 @@ export fixedpoint
 
 export absorb_weight
 export ALSTruncation, FullEnvTruncation
-export SimpleUpdate
+export SimpleUpdate, FullUpdate
 export TimeEvolver, timestep, time_evolve
 
 export InfiniteSquareNetwork

src/algorithms/contractions/bondenv/als_solve.jl

@@ -206,3 +206,13 @@ function _solve_ab(
     x1, info = linsolve(f, Sx, x0, 0, 1)
     return x1, info
 end
+
+function _solve_ab_pinv!(
+        Rx::AbstractTensorMap{T, S, 2, 2}, Sx::AbstractTensorMap{T, S, 2, 1}; kwargs...
+    ) where {T <: Number, S <: ElementarySpace}
+    Rx_inv, ϵ = _pinv!(copy(Rx); kwargs...)
+    is = filter(i -> isdual(codomain(Rx_inv, i)), 1:numout(Rx_inv))
+    x = Rx_inv * Sx
+    twist!(x, is)
+    return x, ϵ
+end

src/algorithms/contractions/bondenv/benv_ctm.jl

@@ -8,7 +8,8 @@ Construct the environment (norm) tensor
     C4---T3---------T3---C3   r+1
     c-1  c         c+1  c+2
 ```
-where `XX = X' X` and `YY = Y' Y` (stacked together).
+where `X, Y` are unitary tensors produced when finding the reduced site tensors 
+with `_qr_bond`; and `XX = X' X` and `YY = Y' Y` (stacked together).
 
 Axis order: `[DX1 DY1; DX0 DY0]`, as in
 ```
@@ -21,7 +22,9 @@ Axis order: `[DX1 DY1; DX0 DY0]`, as in
     └---------------------┘
 ```
 """
-function bondenv_fu(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv)
+function bondenv_fu(
+        row::Int, col::Int, X::PO, Y::PO, env::CTMRGEnv
+    ) where {PO <: Union{PEPSOrth, PEPOOrth}}
     Nr, Nc = size(env.corners)[[2, 3]]
     cm1 = _prev(col, Nc)
     cp1 = _next(col, Nc)
@@ -49,12 +52,23 @@ function bondenv_fu(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv)
     C4--χ3--T3X---------χ5---------T3Y--χ7---C3     r+1
     c-1      c                      c+1     c+2
     =#
-    @autoopt @tensor benv[DX1 DY1; DX0 DY0] :=
-        c4[χ3 χ1] * t4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * t3X[χ5 DSX0 DSX1 χ3] *
-        X[DNX0 DX0 DSX0 DWX0] * conj(X[DNX1 DX1 DSX1 DWX1]) * t1X[χ4 DNX0 DNX1 χ6] *
-        c3[χ9 χ7] * t2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * t3Y[χ7 DSY0 DSY1 χ5] *
-        Y[DNY0 DEY0 DSY0 DY0] * conj(Y[DNY1 DEY1 DSY1 DY1]) * t1Y[χ6 DNY0 DNY1 χ8]
-
+    benv = nothing
+    if PO <: PEPSOrth
+        @autoopt @tensor benv[DX1 DY1; DX0 DY0] :=
+            c4[χ3 χ1] * t4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * t3X[χ5 DSX0 DSX1 χ3] *
+            X[DNX0 DX0 DSX0 DWX0] * conj(X[DNX1 DX1 DSX1 DWX1]) * t1X[χ4 DNX0 DNX1 χ6] *
+            c3[χ9 χ7] * t2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * t3Y[χ7 DSY0 DSY1 χ5] *
+            Y[DNY0 DEY0 DSY0 DY0] * conj(Y[DNY1 DEY1 DSY1 DY1]) * t1Y[χ6 DNY0 DNY1 χ8]
+    else
+        # eliminate fermion sign when contracting the remaining physical axis in X, Y
+        X2 = isdual(space(X, 1)) ? twist(X, 1) : X
+        Y2 = isdual(space(Y, 1)) ? twist(Y, 1) : Y
+        @autoopt @tensor benv[DX1 DY1; DX0 DY0] :=
+            c4[χ3 χ1] * t4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * t3X[χ5 DSX0 DSX1 χ3] *
+            X2[pX DNX0 DX0 DSX0 DWX0] * conj(X[pX DNX1 DX1 DSX1 DWX1]) * t1X[χ4 DNX0 DNX1 χ6] *
+            c3[χ9 χ7] * t2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * t3Y[χ7 DSY0 DSY1 χ5] *
+            Y2[pY DNY0 DEY0 DSY0 DY0] * conj(Y[pY DNY1 DEY1 DSY1 DY1]) * t1Y[χ6 DNY0 DNY1 χ8]
+    end
     normalize!(benv, Inf)
     return benv
 end

src/algorithms/contractions/bondenv/gaugefix.jl

@@ -92,12 +92,13 @@ function fixgauge_benv(
 end
 
 """
-When the (half) bond environment `Z` consists of two `PEPSOrth` tensors `X`, `Y` as
+When the (half) bond environment `Z` consists of 
+two `PEPSOrth` or `PEPOOrth` tensors `X`, `Y` as
 ```
-    ┌---------------┐   ┌-------------------┐
-    |               | = |                   | ,
-    └---Z--       --┘   └--Z0---X--    --Y--┘
-        ↓                  ↓
+    ┌-----------------------┐
+    |                       |
+    └---Z---(X)--   --(Y)---┘
+        ↓
 ```
 apply the gauge transformation `Linv`, `Rinv` for `Z` to `X`, `Y`:
 ```
@@ -106,15 +107,25 @@ apply the gauge transformation `Linv`, `Rinv` for `Z` to `X`, `Y`:
     -4 - X - 1 - Rinv - -2      -4 - Linv - 1 - Y - -2
          |                                      |
         -3                                     -3
+
+        -2                                     -2
+         |                                      |
+    -5 - X - 1 - Rinv - -3      -5 - Linv - 1 - Y - -3
+         | ╲                                    | ╲
+        -4  -1                                 -4  -1
 ```
 """
 function _fixgauge_benvXY(
-        X::PEPSOrth{T, S},
-        Y::PEPSOrth{T, S},
-        Linv::AbstractTensorMap{T, S, 1, 1},
-        Rinv::AbstractTensorMap{T, S, 1, 1},
-    ) where {T <: Number, S <: ElementarySpace}
+        X::PEPSOrth, Y::PEPSOrth, Linv::MPSBondTensor, Rinv::MPSBondTensor,
+    )
     @plansor X[-1 -2 -3 -4] := X[-1 1 -3 -4] * Rinv[1; -2]
     @plansor Y[-1 -2 -3 -4] := Y[-1 -2 -3 1] * Linv[1; -4]
     return X, Y
 end
+function _fixgauge_benvXY(
+        X::PEPOOrth, Y::PEPOOrth, Linv::MPSBondTensor, Rinv::MPSBondTensor,
+    )
+    @plansor X[-1 -2 -3 -4 -5] := X[-1 -2 1 -4 -5] * Rinv[1; -3]
+    @plansor Y[-1 -2 -3 -4 -5] := Y[-1 -2 -3 -4 1] * Linv[1; -5]
+    return X, Y
+end

src/algorithms/ctmrg/sequential.jl

@@ -57,6 +57,65 @@ function ctmrg_leftmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG)
     return env, info
 end
 
+"""
+    ctmrg_rightmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG)
+
+Perform sequential CTMRG right move on the `col`-th column.
+"""
+function ctmrg_rightmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG)
+    #= 
+        right move <---
+        ←-- T1 ← C2     r-1
+            ‖    ↑
+        === M' = T2     r
+            ‖    ↑
+        --→ T3 → C3     r+1
+            c   c+1
+    =#
+    Nc = size(network)[2]
+    @assert 1 <= col <= Nc
+    env, info = ctmrg_leftmove(Nc + 1 - col, rot180(network), rot180(env), alg)
+    return rot180(env), info
+end
+
+"""
+    ctmrg_upmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG)
+
+Perform sequential CTMRG up move on the `row`-th row.
+"""
+function ctmrg_upmove(row::Int, network, env::CTMRGEnv, alg::SequentialCTMRG)
+    #= 
+            c-1  c   c+1
+        r-1 C1 ← T1 ← C2
+            ↓    ‖    ↑     |  up
+        r   T4 = M == T2    ↓  move
+            ↓    ‖    ↑
+    =#
+    Nr = size(network)[1]
+    @assert 1 <= row <= Nr
+    env, info = ctmrg_leftmove(Nr + 1 - row, rotl90(network), rotl90(env), alg)
+    return rotr90(env), info
+end
+
+"""
+    ctmrg_downmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG)
+
+Perform sequential CTMRG down move on the `row`-th row.
+"""
+function ctmrg_downmove(row::Int, network, env::CTMRGEnv, alg::SequentialCTMRG)
+    #= 
+            ↓    ‖    ↑
+        r   T4 = M == T2    ↑  down
+            ↓    ‖    ↑     |  move
+        r+1 C4 → T3 → C3
+            c-1  c   c+1
+    =#
+    Nr = size(network)[1]
+    @assert 1 <= row <= Nr
+    env, info = ctmrg_leftmove(row, rotr90(network), rotr90(env), alg)
+    return rotl90(env), info
+end
+
 function ctmrg_iteration(network, env::CTMRGEnv, alg::SequentialCTMRG)
     truncation_error = zero(real(scalartype(network)))
     condition_number = zero(real(scalartype(network)))