23 using BRanks =
typename Tb::Ranks;
25 auto const toleranceLocal =
30 [&](
auto const... seq)
33 MHNop.template core<seq>()...);
41 decltype(x0.template core<k>())>::size(0),
44 decltype(x0.template core<k>())>::size(2)>(
45 x.template core<k>(), 0, 0, 0) = x0.template core<k>();
54 auto z = round<ZRanks>(
74 reshape<1, 1>(std::get<0>(psiphiA)) =
76 reshape<1, 1>(std::get<D>(psiphiA)) =
91 reshape<1, 1>(std::get<0>(psiphib)) =
93 reshape<1, 1>(std::get<D>(psiphib)) =
111 reshape<1, 1>(std::get<0>(psiphihatA)) =
113 reshape<1, 1>(std::get<D>(psiphihatA)) =
118 [](
auto const... seq)
128 reshape<1, 1>(std::get<0>(psiphihatb)) =
130 reshape<1, 1>(std::get<D>(psiphihatb)) =
139 orthogonalize<pRC::Position::Right>(
x.template core<k>());
141 x.template core<k>() = q;
142 x.template core<k - 1>() =
143 contract<2, 0>(
x.template core<k - 1>(), l);
146 auto const [l, q] = orthogonalize<pRC::Position::Right>(
147 z.template core<k>());
148 z.template core<k>() = q;
149 z.template core<k - 1>() =
150 contract<2, 0>(
z.template core<k - 1>(), l);
154 std::get<k>(psiphiA) =
155 contract<1, 2, 1, 3>(conj(
x.template core<k>()),
156 contract<2, 3, 1, 3>(op.template core<k>(),
157 eval(contract<2, 2>(
x.template core<k>(),
158 std::get<k + 1>(psiphiA)))));
160 std::get<k>(psiphib) =
161 contract<1, 2, 1, 2>(conj(
x.template core<k>()),
162 eval(contract<2, 1>(b.template core<k>(),
163 std::get<k + 1>(psiphib))));
165 std::get<k>(psiphihatA) =
166 contract<1, 2, 1, 3>(conj(
z.template core<k>()),
167 contract<2, 3, 1, 3>(op.template core<k>(),
168 eval(contract<2, 2>(
x.template core<k>(),
169 get<k + 1>(psiphihatA)))));
171 std::get<k>(psiphihatb) =
172 contract<1, 2, 1, 2>(conj(
z.template core<k>()),
173 eval(contract<2, 1>(b.template core<k>(),
174 get<k + 1>(psiphihatb))));
187 decltype(op.template core<k>())>::size(0),
190 decltype(op.template core<k + 1>())>::size(3)>(
191 permute<0, 1, 3, 2, 4, 5>(contract<3, 0>(
192 op.template core<k>(), op.template core<k + 1>())));
196 decltype(b.template core<k>())>::size(0),
199 decltype(b.template core<k + 1>())>::size(2)>(
200 contract<2, 0>(b.template core<k>(),
201 b.template core<k + 1>()));
205 decltype(
x.template core<k>())>::size(0),
208 decltype(
x.template core<k + 1>())>::size(2)>(
209 contract<2, 0>(
x.template core<k>(),
210 x.template core<k + 1>()));
212 auto const Ak = matricize(permute<0, 2, 4, 1, 3, 5>(
213 contract<1, 0>(std::get<k>(psiphiA),
214 eval(contract<3, 1>(sA,
215 std::get<k + 2>(psiphiA))))));
217 auto const bk = linearize(contract<1, 0>(
218 std::get<k>(psiphib),
219 eval(contract<2, 1>(sB, std::get<k + 2>(psiphib)))));
221 auto const hatAk = matricize(permute<0, 2, 4, 1, 3, 5>(
222 contract<1, 0>(std::get<k>(psiphihatA),
223 eval(contract<3, 1>(sA,
224 std::get<k + 2>(psiphihatA))))));
226 auto const hatbk = linearize(contract<1, 0>(
227 std::get<k>(psiphihatb),
228 eval(contract<2, 1>(sB, std::get<k + 2>(psiphihatb)))));
230 auto const rAk = matricize(permute<0, 2, 4, 1, 3, 5>(
231 contract<1, 0>(std::get<k>(psiphiA),
232 eval(contract<3, 1>(sA,
233 std::get<k + 2>(psiphihatA))))));
235 auto const rbk = linearize(contract<1, 0>(
236 std::get<k>(psiphib),
237 eval(contract<2, 1>(sB, std::get<k + 2>(psiphihatb)))));
242 decltype(
x.template core<k>())>::size(0),
245 decltype(
x.template core<k + 1>())>::size(2)>
248 linearize(sX), toleranceLocal);
251 auto const [u, s, v] = svd<Ranks::size(k)>(matricize(sol));
255 decltype(
x.template core<k>())>::size(0),
260 decltype(
x.template core<k>())>::size(0),
261 2, Ranks::size(k)>(ex, 0, 0, 0)) = linearize(u);
263 if constexpr(ERanks::size(k) - Ranks::size(k) > 0)
266 decltype(
x.template core<k>())>::size(0);
268 decltype(
z.template core<k + 1>())>::size(2);
270 auto const tt = eval(rbk - rAk * linearize(sol));
273 auto const [zu, zs, zv] =
274 svd<ERanks::size(k) - Ranks::size(k)>(
279 decltype(
x.template core<k>())>::size(0),
280 2, ERanks::size(k) - Ranks::size(k)>(ex, 0, 0,
282 linearize(zu * fromDiagonal(zs));
287 decltype(
z.template core<k>())>::size(0),
290 decltype(
z.template core<k + 1>())>::size(2)>
292 linearize(zz) = hatbk - hatAk * linearize(sol);
294 auto const [zu, zs, zv] =
295 svd<ZRanks::size(k)>(matricize(zz));
297 linearize(
z.template core<k>()) = linearize(zu);
301 decltype(
x.template core<k + 1>())>::size(2)>
303 folding<pRC::Position::Right>(tNext) =
304 fromDiagonal(s) * adjoint(v);
307 decltype(
x.template core<k + 1>())>::size(2)>
309 slice<Ranks::size(k), 2,
311 decltype(
x.template core<k + 1>())>::size(2)>(eNext,
313 slice<ERanks::size(k) - Ranks::size(k), 2,
315 decltype(
x.template core<k + 1>())>::size(2)>(eNext,
318 auto const [qq, rr] =
319 orthogonalize<pRC::Position::Left>(ex);
320 x.template core<k>() = qq;
321 x.template core<k + 1>() = contract<1, 0>(rr, eNext);
324 std::get<k + 1>(psiphiA) =
325 contract<1, 0, 0, 3>(conj(
x.template core<k>()),
326 contract<2, 0, 0, 3>(op.template core<k>(),
327 eval(contract<0, 2>(
x.template core<k>(),
328 std::get<k>(psiphiA)))));
331 std::get<k + 1>(psiphib) =
332 contract<0, 1, 2, 0>(conj(
x.template core<k>()),
333 eval(contract<0, 1>(b.template core<k>(),
334 std::get<k>(psiphib))));
337 std::get<k + 1>(psiphihatA) =
338 contract<1, 0, 0, 3>(conj(
z.template core<k>()),
339 contract<2, 0, 0, 3>(op.template core<k>(),
340 eval(contract<0, 2>(
x.template core<k>(),
341 std::get<k>(psiphihatA)))));
344 std::get<k + 1>(psiphihatb) =
345 contract<0, 1, 2, 0>(conj(
z.template core<k>()),
346 eval(contract<0, 1>(b.template core<k>(),
347 std::get<k>(psiphihatb))));
358 decltype(op.template core<k - 1>())>::size(0),
361 decltype(op.template core<k>())>::size(3)>(
362 permute<0, 1, 3, 2, 4, 5>(contract<3, 0>(
363 op.template core<k - 1>(), op.template core<k>())));
367 decltype(b.template core<k - 1>())>::size(0),
370 decltype(b.template core<k>())>::size(2)>(
371 contract<2, 0>(b.template core<k - 1>(),
372 b.template core<k>()));
376 decltype(
x.template core<k - 1>())>::size(0),
379 decltype(
x.template core<k>())>::size(2)>(
380 contract<2, 0>(
x.template core<k - 1>(),
381 x.template core<k>()));
383 auto const Ak = matricize(permute<0, 2, 4, 1, 3, 5>(
384 contract<1, 0>(std::get<k - 1>(psiphiA),
385 eval(contract<3, 1>(sA,
386 std::get<k + 1>(psiphiA))))));
388 auto const bk = linearize(contract<1, 0>(
389 std::get<k - 1>(psiphib),
390 eval(contract<2, 1>(sB, std::get<k + 1>(psiphib)))));
392 auto const hatAk = matricize(permute<0, 2, 4, 1, 3, 5>(
393 contract<1, 0>(std::get<k - 1>(psiphihatA),
394 eval(contract<3, 1>(sA,
395 std::get<k + 1>(psiphihatA))))));
397 auto const hatbk = linearize(contract<1, 0>(
398 std::get<k - 1>(psiphihatb),
399 eval(contract<2, 1>(sB, std::get<k + 1>(psiphihatb)))));
401 auto const rAk = matricize(permute<0, 2, 4, 1, 3, 5>(
402 contract<1, 0>(std::get<k - 1>(psiphihatA),
403 eval(contract<3, 1>(sA,
404 std::get<k + 1>(psiphiA))))));
406 auto const rbk = linearize(contract<1, 0>(
407 std::get<k - 1>(psiphihatb),
408 eval(contract<2, 1>(sB, std::get<k + 1>(psiphib)))));
413 decltype(
x.template core<k - 1>())>::size(0),
416 decltype(
x.template core<k>())>::size(2)>
419 linearize(sX), tolerance);
422 auto const [u, s, v] =
423 svd<Ranks::size(k - 1)>(matricize(sol));
427 decltype(
x.template core<k>())>::size(2)>
429 linearize(slice<Ranks::size(k - 1), 2,
431 decltype(
x.template core<k>())>::size(2)>(ex, 0, 0,
432 0)) = linearize(adjoint(v));
434 if constexpr(ERanks::size(k - 1) - Ranks::size(k - 1) > 0)
437 decltype(
z.template core<k - 1>())>::size(0);
439 decltype(
x.template core<k>())>::size(2);
441 auto const tt = eval(rbk - rAk * linearize(sol));
444 auto const [zu, zs, zv] =
445 svd<ERanks::size(k - 1) - Ranks::size(k - 1)>(
449 slice<ERanks::size(k - 1) - Ranks::size(k - 1), 2,
451 decltype(
x.template core<k>())>::size(2)>(
452 ex, Ranks::size(k - 1), 0, 0)) =
453 linearize(fromDiagonal(zs) * adjoint(zv));
458 decltype(
z.template core<k - 1>())>::size(0),
461 decltype(
z.template core<k>())>::size(2)>
463 linearize(zz) = hatbk - hatAk * linearize(sol);
465 auto const [zu, zs, zv] =
466 svd<ZRanks::size(k - 1)>(matricize(zz));
468 linearize(
z.template core<k>()) = linearize(adjoint(zv));
472 decltype(
x.template core<k - 1>())>::size(0),
473 2, Ranks::size(k - 1)>
475 folding<pRC::Position::Left>(tNext) = u * fromDiagonal(s);
478 decltype(
x.template core<k - 1>())>::size(0),
479 2, ERanks::size(k - 1)>
482 decltype(
x.template core<k - 1>())>::size(0),
483 2, Ranks::size(k - 1)>(eNext, 0, 0, 0) = tNext;
486 decltype(
x.template core<k - 1>())>::size(0),
487 2, ERanks::size(k - 1) - Ranks::size(k - 1)>(eNext,
488 0, 0, Ranks::size(k - 1))) =
pRC::zero();
490 auto const [ll, qq] =
491 orthogonalize<pRC::Position::Right>(ex);
492 x.template core<k>() = qq;
493 x.template core<k - 1>() = contract<2, 0>(eNext, ll);
497 std::get<k>(psiphiA) =
498 contract<1, 2, 1, 3>(conj(
x.template core<k>()),
499 contract<2, 3, 1, 3>(op.template core<k>(),
500 eval(contract<2, 2>(
x.template core<k>(),
501 std::get<k + 1>(psiphiA)))));
504 std::get<k>(psiphib) =
505 contract<1, 2, 1, 2>(conj(
x.template core<k>()),
506 eval(contract<2, 1>(b.template core<k>(),
507 std::get<k + 1>(psiphib))));
510 std::get<k>(psiphihatA) =
511 contract<1, 2, 1, 3>(conj(
z.template core<k>()),
512 contract<2, 3, 1, 3>(op.template core<k>(),
513 eval(contract<2, 2>(
x.template core<k>(),
514 std::get<k + 1>(psiphihatA)))));
517 std::get<k>(psiphihatb) =
518 contract<1, 2, 1, 2>(conj(
z.template core<k>()),
519 eval(contract<2, 1>(b.template core<k>(),
520 std::get<k + 1>(psiphihatb))));
523 return round<Ranks>(
x);