Skip to content

Commit dcb2587

Browse files
committed
Handle more edge cases
1 parent b058626 commit dcb2587

File tree

1 file changed

+87
-42
lines changed

1 file changed

+87
-42
lines changed

src/operators/ortho.jl

Lines changed: 87 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -15,58 +15,83 @@ function left_canonicalize!(
1515
@plansor C′[p; p' r] := -WI[p; p' l] * t[l; r]
1616
add!(C′, W.C)
1717

18-
# QR of second colum
19-
if i == 1
18+
# QR of second column
19+
if size(W, 1) == 1
2020
tmp = transpose(C′, ((2, 1), (3,)))
21+
2122
if trscheme == notrunc()
2223
Q, R = leftorth!(tmp; alg)
2324
else
2425
@assert alg == SVD() || alg == SDD()
25-
# TODO: truncation might not be using normalized tensors?
2626
Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme)
2727
R = Σ * Vᴴ
2828
end
29-
scale!(Q, d)
30-
scale!(R, inv(d))
31-
32-
V = BlockTensorKit.oplus(oneunit(S), space(R, 1), oneunit(S))
3329

34-
Q1 = typeof(W.A)(undef, SumSpace{S}() physicalspace(W) physicalspace(W) space(R, 1))
35-
Q2 = transpose(Q, ((2,), (1, 3)))
30+
if dim(R) == 0 # fully truncated
31+
V = BlockTensorKit.oplus(oneunit(S), oneunit(S))
32+
Q1 = typeof(W.A)(undef, SumSpace{S}() physicalspace(W) physicalspace(W) SumSpace{S}())
33+
Q2 = typeof(W.C)(undef, physicalspace(W) physicalspace(W) SumSpace{S}())
34+
else
35+
V = BlockTensorKit.oplus(oneunit(S), space(R, 1), oneunit(S))
36+
scale!(Q, d)
37+
scale!(R, inv(d))
38+
Q1 = typeof(W.A)(undef, SumSpace{S}() physicalspace(W) physicalspace(W) space(R, 1))
39+
Q2 = transpose(Q, ((2,), (1, 3)))
40+
end
41+
H[i] = JordanMPOTensor(codomain(W) physicalspace(W) V, Q1, W.B, Q2, W.D)
3642
else
3743
tmp = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,)))
3844
if trscheme == notrunc()
3945
Q, R = leftorth!(tmp; alg)
4046
else
4147
@assert alg == SVD() || alg == SDD()
42-
# TODO: truncation might not be using normalized tensors?
4348
Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme)
4449
R = Σ * Vᴴ
4550
end
46-
scale!(Q, d)
47-
scale!(R, inv(d))
48-
Q′ = transpose(Q, ((2, 3), (1, 4)))
49-
Q1 = Q′[2:end, 1, 1, 1]
50-
Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1)
51-
V = BlockTensorKit.oplus(oneunit(S), right_virtualspace(Q′), oneunit(S))
51+
if dim(R) == 0 # fully truncated
52+
V = BlockTensorKit.oplus(oneunit(S), oneunit(S))
53+
Q1 = typeof(W.A)(undef, SumSpace{S}() physicalspace(W) physicalspace(W) SumSpace{S}())
54+
Q2 = typeof(W.C)(undef, physicalspace(W) physicalspace(W) SumSpace{S}())
55+
else
56+
scale!(Q, d)
57+
scale!(R, inv(d))
58+
Q′ = transpose(Q, ((2, 3), (1, 4)))
59+
Q1 = Q′[2:end, 1, 1, 1]
60+
Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1)
61+
V = BlockTensorKit.oplus(oneunit(S), right_virtualspace(Q′), oneunit(S))
62+
end
63+
H[i] = JordanMPOTensor(codomain(W) physicalspace(W) V, Q1, W.B, Q2, W.D)
5264
end
53-
H[i] = JordanMPOTensor(codomain(W) physicalspace(W) V, Q1, W.B, Q2, W.D)
5465

5566
# absorb into next site
5667
W′ = H[i + 1]
57-
if i != length(H) - 1
68+
69+
if size(W′, 4) > 1 && dim(R) != 0
5870
@plansor A′[l p; p' r] := R[l; r'] * W′.A[r' p; p' r]
71+
else
72+
A′ = similar(W′.A, right_virtualspace(H[i].A) physicalspace(W′) domain(W′.A))
73+
end
74+
75+
if size(W′, 4) > 1
5976
@plansor C′[l p; p' r] := t[l; r'] * W′.A[r' p; p' r]
6077
C′ = add!(removeunit(C′, 1), W′.C)
6178
else
62-
A′ = similar(W′.A, space(R, 1) physicalspace(W′) domain(W′.A))
63-
C′ = W′.C
79+
C′ = W′.C # empty
6480
end
65-
@plansor B′[l p; p'] := R[l; r] * W′.B[r p; p']
81+
82+
if dim(R) != 0
83+
@plansor B′[l p; p'] := R[l; r] * W′.B[r p; p']
84+
else
85+
B′ = similar(W′.B, right_virtualspace(H[i].A) physicalspace(W′) domain(W′.B))
86+
end
87+
6688
@plansor D′[l p; p'] := t[l; r] * W′.B[r p; p']
6789
D′ = add!(removeunit(D′, 1), W′.D)
6890

69-
H[i + 1] = JordanMPOTensor(V physicalspace(W′) domain(W′), A′, B′, C′, D′)
91+
H[i + 1] = JordanMPOTensor(
92+
right_virtualspace(H[i]) physicalspace(W′) domain(W′),
93+
A′, B′, C′, D′
94+
)
7095
return H
7196
end
7297

@@ -87,56 +112,76 @@ function right_canonicalize!(
87112
add!(B′, W.B)
88113

89114
# LQ of second row
90-
if i == length(H)
115+
if size(W, 4) == 1
91116
tmp = transpose(B′, ((1,), (3, 2)))
92117
if trscheme == notrunc()
93118
R, Q = rightorth!(tmp; alg)
94119
else
95120
@assert alg == SVD() || alg == SDD()
96-
# TODO: truncation might not be using normalized tensors?
97121
U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme)
98122
R = U * Σ
99123
end
100-
scale!(Q, d)
101-
scale!(R, inv(d))
102124

103-
R, Q = rightorth!(transpose(B′, ((1,), (3, 2))); alg)
104-
V = BlockTensorKit.oplus(oneunit(S), space(Q, 1), oneunit(S))
105-
Q1 = typeof(W.A)(undef, space(Q, 1) physicalspace(W) physicalspace(W) SumSpace{S}())
106-
Q2 = transpose(Q, ((1, 3), (2,)))
125+
if dim(R) == 0
126+
V = BlockTensorKit.oplus(oneunit(S), oneunit(S))
127+
Q1 = typeof(W.A)(undef, SumSpace{S}() physicalspace(W) physicalspace(W) SumSpace{S}())
128+
Q2 = typeof(W.B)(undef, SumSpace{S}() physicalspace(W) physicalspace(W))
129+
else
130+
V = BlockTensorKit.oplus(oneunit(S), space(Q, 1), oneunit(S))
131+
scale!(Q, d)
132+
scale!(R, inv(d))
133+
Q1 = typeof(W.A)(undef, space(Q, 1) physicalspace(W) physicalspace(W) SumSpace{S}())
134+
Q2 = transpose(Q, ((1, 3), (2,)))
135+
end
136+
H[i] = JordanMPOTensor(V physicalspace(W) domain(W), Q1, Q2, W.C, W.D)
107137
else
108138
tmp = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2)))
109139
if trscheme == notrunc()
110140
R, Q = rightorth!(tmp; alg)
111141
else
112142
@assert alg == SVD() || alg == SDD()
113-
# TODO: truncation might not be using normalized tensors?
114143
U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme)
115144
R = U * Σ
116145
end
117-
scale!(Q, d)
118-
scale!(R, inv(d))
119-
Q′ = transpose(Q, ((1, 4), (2, 3)))
120-
Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end])
121-
Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4)
122-
V = BlockTensorKit.oplus(oneunit(S), left_virtualspace(Q′), oneunit(S))
146+
if dim(R) == 0
147+
V = BlockTensorKit.oplus(oneunit(S), oneunit(S))
148+
Q1 = typeof(W.A)(undef, SumSpace{S}() physicalspace(W) physicalspace(W) SumSpace{S}())
149+
Q2 = typeof(W.B)(undef, SumSpace{S}() physicalspace(W) physicalspace(W))
150+
else
151+
scale!(Q, d)
152+
scale!(R, inv(d))
153+
Q′ = transpose(Q, ((1, 4), (2, 3)))
154+
Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end])
155+
Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4)
156+
V = BlockTensorKit.oplus(oneunit(S), left_virtualspace(Q′), oneunit(S))
157+
end
158+
H[i] = JordanMPOTensor(V physicalspace(W) domain(W), Q1, Q2, W.C, W.D)
123159
end
124-
H[i] = JordanMPOTensor(V physicalspace(W) domain(W), Q1, Q2, W.C, W.D)
125160

126161
# absorb into previous site
127162
W′ = H[i - 1]
128-
if i != 2
163+
164+
if size(W′, 1) > 1 && dim(R) != 0
129165
@plansor A′[l p; p' r] := W′.A[l p; p' r'] * R[r'; r]
166+
else
167+
A′ = similar(W′.A, codomain(W′.A) physicalspace(W′.A) left_virtualspace(H[i].A))
168+
end
169+
170+
if size(W′, 1) > 1
130171
@plansor B′[l p; p' r] := W′.A[l p; p' r'] * t[r'; r]
131172
B′ = add!(removeunit(B′, 4), W′.B)
132173
else
133-
A′ = similar(W′.A, codomain(W′.A) physicalspace(W′.A) space(R, 2)')
134174
B′ = W′.B
135175
end
136-
@plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r]
176+
177+
if dim(R) != 0
178+
@plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r]
179+
else
180+
C′ = similar(W′.C, codomain(W′.C) physicalspace(W′) left_virtualspace(H[i].A))
181+
end
182+
137183
@plansor D′[p; p' r] := W′.C[p; p' r'] * t[r'; r]
138184
D′ = add!(removeunit(D′, 3), W′.D)
139-
140185
H[i - 1] = JordanMPOTensor(codomain(W′) physicalspace(W′) V, A′, B′, C′, D′)
141186
return H
142187
end

0 commit comments

Comments
 (0)