Add right_canonicalize!

Author: lkdvos

src/operators/ortho.jl

@@ -63,3 +63,36 @@ function left_canonicalize!(H::MPOHamilonian, i::Int; alg = QRPos())
     H[i + 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′)
     return H
 end
+
+function right_canonicalize!(H::MPOHamilonian, i::Int; alg = RQpos())
+    @assert i != length(H) "TBA"
+
+    W = H[i]
+
+    # orthogonalize second row against last
+    WI = removeunit(W[end, 1, 1, end], 4)
+    @tensor t[l; r] := conj(WI[r p; p']) * W.B[l p; p']
+    @plansor B′[l p; p'] := -WI[r p; p'] * t[l; r]
+    add!(B′, W.B)
+
+    # LQ of second row
+    AB = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2)))
+    R, Q = rightorth!(AB; alg)
+    Q′ = transpose(Q, ((1, 4), (2, 3)))
+    Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end])
+    Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 4)
+    V = BlockTensorKit.oplus(oneunit(spacetype(W)), left_virtualspace(Q′), oneunit(spacetype(W)))
+    H[i] = JordanMPOTensor(V ⊗ physicalspace(W) ← domain(W), Q1, Q2, W.C, W.D)
+
+    # absorb into previous site
+    W′ = H[i - 1]
+    @plansor A′[l p; p' r] := W′.A[l p; p' r] * R[r; r] 
+    @plansor B′[l p; p' r'] := W′.A[l p; p' r'] * t[r'; r] R[l; r] * W′.B[r p; p']
+    B′ = add!(removeunit(B′, 4), W′.B)
+    @plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r]
+    @plansor D′[p; p' r] := W′.C[p; p' l] *t[l; r]  
+    D′ = add!(removeunit(D, 3), W′.D)
+
+    H[i - 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′)
+    return H
+end