VMPS++
Loading...
Searching...
No Matches
DmrgContractions.h
Go to the documentation of this file.
1#ifndef STRAWBERRY_DMRGCONTRACTIONS_WITH_Q
2#define STRAWBERRY_DMRGCONTRACTIONS_WITH_Q
3
5#include <unordered_set>
7
8//include "Biped.h"
9//include "Multipede.h"
10#include "tensors/DmrgIndexGymnastics.h"
11//include "symmetry/functions.h"
12
13enum CONTRACT_LR_MODE {FULL, TRIANGULAR, FIXED_ROWS, FIXED_COLS};
14
34template<typename Symmetry, typename MatrixType, typename MatrixType2, typename MpoMatrixType>
35void contract_L (const Tripod<Symmetry,MatrixType2> &Lold,
36 const vector<Biped<Symmetry,MatrixType> > &Abra,
37 const vector<vector<vector<Biped<Symmetry,MpoMatrixType> > > > &W,
38 const vector<Biped<Symmetry,MatrixType> > &Aket,
39 const vector<qarray<Symmetry::Nq> > &qloc,
40 const vector<qarray<Symmetry::Nq> > &qOp,
41 Tripod<Symmetry,MatrixType2> &Lnew,
42 bool RANDOMIZE = false,
43 tuple<CONTRACT_LR_MODE,size_t> MODE_input = make_pair(FULL,0),
44 const std::unordered_map<pair<qarray<Symmetry::Nq>,size_t>,size_t> &basis_order_map = {})
45{
46 Lnew.clear();
47 Lnew.setZero();
48 auto [MODE,fixed_ab] = MODE_input;
49
50 #ifdef DMRG_CONTRACTLANDR_PARALLELIZE
51 #pragma omp parallel for collapse(3) schedule(dynamic)
52 #endif
53 for (size_t s1=0; s1<qloc.size(); ++s1)
54 for (size_t s2=0; s2<qloc.size(); ++s2)
55 for (size_t k=0; k<qOp.size(); ++k)
56 {
57 typename MpoMatrixType::Scalar factor_cgc;
58 if (!Symmetry::validate(qarray3<Symmetry::Nq>{qloc[s2], qOp[k], qloc[s1]})) {continue;}
59
60 for (size_t qL=0; qL<Lold.dim; ++qL)
61 {
62 vector<tuple<qarray3<Symmetry::Nq>,size_t,size_t,size_t> > ix;
63 bool FOUND_MATCH = LAWA(Lold.in(qL), Lold.out(qL), Lold.mid(qL), qloc[s1], qloc[s2], qOp[k], Abra[s1], Aket[s2], W[s1][s2][k], ix);
64
65 if (FOUND_MATCH)
66 {
67 for(size_t n=0; n<ix.size(); n++ )
68 {
69 qarray3<Symmetry::Nq> quple = get<0>(ix[n]);
70 swap(quple[0], quple[1]);
71 size_t qAbra = get<1>(ix[n]);
72 size_t qAket = get<2>(ix[n]);
73 size_t qW = get<3>(ix[n]);
74
75 if (Aket[s2].block[qAket].size() == 0) {continue;}
76 if (Abra[s1].block[qAbra].size() == 0) {continue;}
77
78 if constexpr ( Symmetry::NON_ABELIAN )
79 {
80 factor_cgc = Symmetry::coeff_buildL(Aket[s2].in[qAket], qloc[s2], Aket[s2].out[qAket],
81 Lold.mid(qL), qOp[k], quple[2],
82 Abra[s1].in[qAbra], qloc[s1], Abra[s1].out[qAbra]);
83 //lout << "cgc=" << factor_cgc << ", " << Aket[s2].in[qAket] << ", " << Aket[s2].out[qAket] << ", " << Lold.mid(qL) << ", " << qOp[k] << ", " << quple[2] << ", " << Abra[s1].in[qAbra] << ", " << Abra[s1].out[qAbra] << endl;
84 }
85 else
86 {
87 factor_cgc = 1.;
88 }
89 if (abs(factor_cgc) < ::mynumeric_limits<double>::epsilon()) {continue;}
90
91 for (int r=0; r<W[s1][s2][k].block[qW].outerSize(); ++r)
92 for (typename MpoMatrixType::InnerIterator iW(W[s1][s2][k].block[qW],r); iW; ++iW)
93 {
94 size_t a = iW.row();
95 size_t b = iW.col();
96
97 // 0 is the Hamiltonian block. Only singlet couplings are neccessary.
98 //if (b == 0 and IS_HAMILTONIAN and quple[2] != Symmetry::qvacuum()) {continue;}
99 // if (MODE == TRIANGULAR) {cout << "before check: fixed_ab=" << fixed_ab << ", a=" << a << ", Lold.mid=" << Lold.mid(qL) << ", basis_number=" << basis_order_map.at(make_pair(Lold.mid(qL),a)) << endl;}
100 if (MODE == FULL or
101 (MODE == TRIANGULAR and basis_order_map.at(make_pair(Lold.mid(qL),a))>fixed_ab) or
102 (MODE == FIXED_COLS and basis_order_map.at(make_pair(quple[2],b))==fixed_ab) or
103 (MODE == FIXED_ROWS and basis_order_map.at(make_pair(Lold.mid(qL),a))==fixed_ab)
104 )
105 {
106// if (MODE == FIXED)
107// {
108// cout << "fixed_b=" << fixed_b << ", a=" << a << "/" << W[s1][s2][k].rows() << ", b=" << b << "/" << W[s1][s2][k].cols() << endl;
109// cout << "Lold.block[qL][a][0].size()=" << Lold.block[qL][a][0].size() << endl;
110// }
111
112 if (Lold.block[qL][a][0].size() != 0)
113 {
114 MatrixType2 Mtmp;
115 if (RANDOMIZE)
116 {
117 Mtmp.resize(Abra[s1].block[qAbra].cols(), Aket[s2].block[qAket].cols());
118 Mtmp.setRandom();
119 }
120 else
121 {
122 optimal_multiply(factor_cgc * iW.value(),
123 Abra[s1].block[qAbra].adjoint(),
124 Lold.block[qL][a][0],
125 Aket[s2].block[qAket],
126 Mtmp);
127 }
128
129 #pragma omp critical
130 {
131 auto it = Lnew.dict.find(quple);
132 if (it != Lnew.dict.end())
133 {
134 if (Lnew.block[it->second][b][0].rows() != Mtmp.rows() or
135 Lnew.block[it->second][b][0].cols() != Mtmp.cols())
136 {
137 Lnew.block[it->second][b][0] = Mtmp;
138 }
139 else
140 {
141 Lnew.block[it->second][b][0] += Mtmp;
142 }
143 }
144 else
145 {
146 boost::multi_array<MatrixType2,LEGLIMIT> Mtmpvec(boost::extents[W[s1][s2][k].block[qW].cols()][1]);
147 Mtmpvec[b][0] = Mtmp;
148 Lnew.push_back(quple, Mtmpvec);
149 }
150 }
151 }
152 }
153 }
154 }
155 }
156 }
157 }
158}
159
179template<typename Symmetry, typename MatrixType, typename MatrixType2, typename MpoMatrixType>
180void contract_R (const Tripod<Symmetry,MatrixType2> &Rold,
181 const vector<Biped<Symmetry,MatrixType> > &Abra,
182 const vector<vector<vector<Biped<Symmetry,MpoMatrixType> > > > &W,
183 const vector<Biped<Symmetry,MatrixType> > &Aket,
184 const vector<qarray<Symmetry::Nq> > &qloc,
185 const vector<qarray<Symmetry::Nq> > &qOp,
186 Tripod<Symmetry,MatrixType2> &Rnew,
187 bool RANDOMIZE = false,
188 tuple<CONTRACT_LR_MODE,size_t> MODE_input = make_pair(FULL,0),
189 const std::unordered_map<pair<qarray<Symmetry::Nq>,size_t>,size_t> &basis_order_map = {})
190{
191 Rnew.clear();
192 Rnew.setZero();
193 auto [MODE,fixed_ab] = MODE_input;
194
195 #ifdef DMRG_CONTRACTLANDR_PARALLELIZE
196 #pragma omp parallel for collapse(3) schedule(dynamic)
197 #endif
198 for (size_t s1=0; s1<qloc.size(); ++s1)
199 for (size_t s2=0; s2<qloc.size(); ++s2)
200 for (size_t k=0; k<qOp.size(); ++k)
201 {
202 typename MpoMatrixType::Scalar factor_cgc;
203 if (!Symmetry::validate(qarray3<Symmetry::Nq>{qloc[s2], qOp[k], qloc[s1]})) {continue;}
204
205 for (size_t qR=0; qR<Rold.dim; ++qR)
206 {
207 vector<tuple<qarray3<Symmetry::Nq>,size_t,size_t,size_t> > ix;
208 //bool FOUND_MATCH = AWAR(Rold.in(qR), Rold.out(qR), Rold.mid(qR), s1, s2, qloc, k, qOp, Abra, Aket, ix, IS_HAMILTONIAN);
209 bool FOUND_MATCH = AWAR(Rold.in(qR), Rold.out(qR), Rold.mid(qR), qloc[s1], qloc[s2], qOp[k], Abra[s1], Aket[s2], W[s1][s2][k], ix);
210
211 if (FOUND_MATCH)
212 {
213 for(size_t n=0; n<ix.size(); n++ )
214 {
215 qarray3<Symmetry::Nq> quple = get<0>(ix[n]);
216 swap(quple[0], quple[1]);
217 size_t qAbra = get<1>(ix[n]);
218 size_t qAket = get<2>(ix[n]);
219 size_t qW = get<3>(ix[n]);
220
221 if (Aket[s2].block[qAket].size() == 0) {continue;}
222 if (Abra[s1].block[qAbra].size() == 0) {continue;}
223
224 if constexpr (Symmetry::NON_ABELIAN)
225 {
226 factor_cgc = Symmetry::coeff_buildR(Aket[s2].in[qAket], qloc[s2], Aket[s2].out[qAket],
227 quple[2] , qOp[k], Rold.mid(qR),
228 Abra[s1].in[qAbra], qloc[s1], Abra[s1].out[qAbra]);
229 }
230 else
231 {
232 factor_cgc = 1.;
233 }
234 if (abs(factor_cgc) < ::mynumeric_limits<double>::epsilon()) {continue;}
235
236 for (int r=0; r<W[s1][s2][k].block[qW].outerSize(); ++r)
237 for (typename MpoMatrixType::InnerIterator iW(W[s1][s2][k].block[qW],r); iW; ++iW)
238 {
239 size_t a = iW.row();
240 size_t b = iW.col();
241
242 // Daux-1 is the Hamiltonian block. Only singlet couplings are neccessary.
243 //if (a == W[s1][s2][k].block[qW].rows()-1 and IS_HAMILTONIAN and quple[2] != Symmetry::qvacuum()) {continue;}
244
245 if (MODE == FULL or
246 (MODE == TRIANGULAR and fixed_ab>basis_order_map.at(make_pair(Rold.mid(qR),b))) or
247 (MODE == FIXED_ROWS and basis_order_map.at(make_pair(quple[2],a))==fixed_ab) or
248 (MODE == FIXED_COLS and basis_order_map.at(make_pair(Rold.mid(qR),b))==fixed_ab)
249 )
250 {
251// if (MODE == FIXED)
252// {
253// cout << "fixed_a=" << fixed_a << ", a=" << a << "/" << W[s1][s2][k].rows() << ", b=" << b << "/" << W[s1][s2][k].cols() << endl;
254// cout << "Rold.block[qR][b][0].rows()=" << Rold.block[qR][b][0].rows() << endl;
255// }
256
257 if (Rold.block[qR][b][0].rows() != 0)
258 {
259 MatrixType2 Mtmp;
260 if (RANDOMIZE)
261 {
262 Mtmp.resize(Aket[s2].block[qAket].rows(), Abra[s1].block[qAbra].rows());
263 Mtmp.setRandom();
264 }
265 else
266 {
267// if (Aket[s2].block[qAket].cols() != Rold.block[qR][b][0].rows() or
268// Rold.block[qR][b][0].cols() != Abra[s1].block[qAbra].adjoint().rows())
269// {
270// print_size(Aket[s2].block[qAket],"Aket[s2].block[qAket]");
271// print_size(Rold.block[qR][b][0],"Rold.block[qR][b][0]");
272// print_size(Abra[s1].block[qAbra].adjoint(),"Abra[s1].block[qAbra].adjoint()");
273// }
274
275 optimal_multiply(factor_cgc * iW.value(),
276 Aket[s2].block[qAket],
277 Rold.block[qR][b][0],
278 Abra[s1].block[qAbra].adjoint(),
279 Mtmp);
280 }
281
282 #pragma omp critical
283 {
284 auto it = Rnew.dict.find(quple);
285 if (it != Rnew.dict.end())
286 {
287 if (Rnew.block[it->second][a][0].rows() != Mtmp.rows() or
288 Rnew.block[it->second][a][0].cols() != Mtmp.cols())
289 {
290 Rnew.block[it->second][a][0] = Mtmp;
291 }
292 else
293 {
294 Rnew.block[it->second][a][0] += Mtmp;
295 }
296 }
297 else
298 {
299 boost::multi_array<MatrixType2,LEGLIMIT> Mtmpvec(boost::extents[W[s1][s2][k].block[qW].rows()][1]);
300 Mtmpvec[a][0] = Mtmp;
301 Rnew.push_back(quple, Mtmpvec);
302 }
303 }
304 }
305 }
306 }
307 }
308 }
309 }
310 }
311}
312
327template<typename Symmetry, typename MatrixType, typename MpoScalar>
328void contract_L (const Tripod<Symmetry,MatrixType> &Lold,
329 const vector<Biped<Symmetry,MatrixType> > &Abra,
330 const unordered_map<tuple<size_t,size_t,size_t,qarray<Symmetry::Nq>,qarray<Symmetry::Nq> >,SparseMatrix<MpoScalar> > &V,
331 const vector<Biped<Symmetry,MatrixType> > &Aket,
332 const vector<qarray<Symmetry::Nq> > &qloc,
333 const vector<qarray<Symmetry::Nq> > &qOp,
334 Tripod<Symmetry,MatrixType> &Lnew)
335{
336 MpoScalar factor_cgc;
337 Lnew.clear();
338 Lnew.setZero();
339
340 for (size_t s1=0; s1<qloc.size(); ++s1)
341 for (size_t s2=0; s2<qloc.size(); ++s2)
342 for (size_t k=0; k<qOp.size(); ++k)
343 {
344 if(!Symmetry::validate(qarray3<Symmetry::Nq>{qloc[s2], qOp[k], qloc[s1]})) {continue;}
345 for (size_t qL=0; qL<Lold.dim; ++qL)
346 {
347 vector<tuple<qarray3<Symmetry::Nq>,size_t,size_t> > ix;
348 bool FOUND_MATCH = LAWA(Lold.in(qL), Lold.out(qL), Lold.mid(qL), s1, s2, qloc, k, qOp, Abra, Aket, ix);
349
350 if (FOUND_MATCH == true)
351 {
352 for(size_t n=0; n<ix.size(); n++ )
353 {
354 qarray3<Symmetry::Nq> quple = get<0>(ix[n]);
355 swap(quple[0], quple[1]);
356 size_t qAbra = get<1>(ix[n]);
357 size_t qAket = get<2>(ix[n]);
358 if constexpr ( Symmetry::NON_ABELIAN )
359 {
360 factor_cgc = Symmetry::coeff_buildL(Aket[s2].in[qAket], qloc[s2], Aket[s2].out[qAket],
361 Lold.mid(qL), qOp[k], quple[2],
362 Abra[s1].in[qAbra], qloc[s1], Abra[s1].out[qAbra]);
363 }
364 else
365 {
366 factor_cgc = 1.;
367 }
368 if (std::abs(factor_cgc) < ::mynumeric_limits<double>::epsilon()) { continue; }
369 auto key = make_tuple(s1,s2,k,Lold.mid(qL),quple[2]);
370 if(auto it=V.find(key); it == V.end()) { continue; }
371 for (int r=0; r<V.at(key).outerSize(); ++r)
372 for (typename SparseMatrix<MpoScalar>::InnerIterator iW(V.at(key),r); iW; ++iW)
373 {
374 size_t a1 = iW.row();
375 size_t a2 = iW.col();
376 if (Lold.block[qL][a1][0].size() != 0)
377 {
378 MatrixType Mtmp;
379 optimal_multiply(factor_cgc*iW.value(),
380 Abra[s1].block[qAbra].adjoint(),
381 Lold.block[qL][a1][0],
382 Aket[s2].block[qAket],
383 Mtmp);
384 // if (Mtmp.norm() < ::mynumeric_limits<MpoScalar>::epsilon()) { continue; }
385
386 auto it = Lnew.dict.find(quple);
387 if (it != Lnew.dict.end())
388 {
389 if (Lnew.block[it->second][a2][0].rows() != Mtmp.rows() or
390 Lnew.block[it->second][a2][0].cols() != Mtmp.cols())
391 {
392 Lnew.block[it->second][a2][0] = Mtmp;
393 }
394 else
395 {
396 Lnew.block[it->second][a2][0] += Mtmp;
397 // if (Lnew.block[it->second][a2][0].norm() < ::mynumeric_limits<MpoScalar>::epsilon())
398 // {
399 // Lnew.block[it->second][a2][0].resize(0,0);
400 // continue;
401 // }
402 }
403 }
404 else
405 {
406 boost::multi_array<MatrixType,LEGLIMIT> Mtmpvec(boost::extents[V.at(key).cols()][1]);
407 Mtmpvec[a2][0] = Mtmp;
408 Lnew.push_back(quple, Mtmpvec);
409 }
410 }
411 }
412 }
413 }
414 }
415 }
416}
417
432template<typename Symmetry, typename MatrixType, typename MpoScalar>
433void contract_R (const Tripod<Symmetry,MatrixType> &Rold,
434 const vector<Biped<Symmetry,MatrixType> > &Abra,
435 const unordered_map<tuple<size_t,size_t,size_t,qarray<Symmetry::Nq>,qarray<Symmetry::Nq> >,SparseMatrix<MpoScalar> > &V,
436 const vector<Biped<Symmetry,MatrixType> > &Aket,
437 const vector<qarray<Symmetry::Nq> > &qloc,
438 const vector<qarray<Symmetry::Nq> > &qOp,
439 Tripod<Symmetry,MatrixType> &Rnew)
440{
441 MpoScalar factor_cgc;
442 Rnew.clear();
443 Rnew.setZero();
444
445 for (size_t s1=0; s1<qloc.size(); ++s1)
446 for (size_t s2=0; s2<qloc.size(); ++s2)
447 for (size_t k=0; k<qOp.size(); ++k)
448 {
449 if(!Symmetry::validate(qarray3<Symmetry::Nq>{qloc[s2],qOp[k],qloc[s1]})) {continue;}
450
451 for (size_t qR=0; qR<Rold.dim; ++qR)
452 {
453 auto qRouts = Symmetry::reduceSilent(Rold.out(qR),Symmetry::flip(qloc[s1]));
454 auto qRins = Symmetry::reduceSilent(Rold.in(qR),Symmetry::flip(qloc[s2]));
455
456 for(const auto& qRout : qRouts)
457 for(const auto& qRin : qRins)
458 {
459 qarray2<Symmetry::Nq> cmp1 = {qRout, Rold.out(qR)};
460 qarray2<Symmetry::Nq> cmp2 = {qRin, Rold.in(qR)};
461
462 auto q1 = Abra[s1].dict.find(cmp1);
463 auto q2 = Aket[s2].dict.find(cmp2);
464
465 if (q1!=Abra[s1].dict.end() and
466 q2!=Aket[s2].dict.end())
467 {
468 qarray<Symmetry::Nq> new_qin = Aket[s2].in[q2->second]; // A.in
469 qarray<Symmetry::Nq> new_qout = Abra[s1].in[q1->second]; // A†.out = A.in
470 auto qRmids = Symmetry::reduceSilent(Rold.mid(qR),Symmetry::flip(qOp[k]));
471
472 for(const auto& new_qmid : qRmids)
473 {
474 qarray3<Symmetry::Nq> quple = {new_qin, new_qout, new_qmid};
475 if constexpr (Symmetry::NON_ABELIAN)
476 {
477 factor_cgc = Symmetry::coeff_buildR(Aket[s2].in[q2->second], qloc[s2], Aket[s2].out[q2->second],
478 quple[2] , qOp[k], Rold.mid(qR),
479 Abra[s1].in[q1->second], qloc[s1], Abra[s1].out[q1->second]);
480 }
481 else
482 {
483 factor_cgc = 1.;
484 }
485 if (std::abs(factor_cgc) < ::mynumeric_limits<double>::epsilon()) { continue; }
486 auto key = make_tuple(s1,s2,k,quple[2],Rold.mid(qR));
487 if(auto it=V.find(key); it == V.end()) { continue; }
488 for (int r=0; r<V.at(key).outerSize(); ++r)
489 for (typename SparseMatrix<MpoScalar>::InnerIterator iW(V.at(key),r); iW; ++iW)
490 {
491 size_t a1 = iW.row();
492 size_t a2 = iW.col();
493
494 if (Rold.block[qR][a2][0].rows() != 0)
495 {
496 MatrixType Mtmp;
497 optimal_multiply(factor_cgc*iW.value(),
498 Aket[s2].block[q2->second],
499 Rold.block[qR][a2][0],
500 Abra[s1].block[q1->second].adjoint(),
501 Mtmp);
502 // if(Mtmp.norm() < ::mynumeric_limits<MpoScalar>::epsilon()) { continue; }
503
504 auto it = Rnew.dict.find(quple);
505 if (it != Rnew.dict.end())
506 {
507 if (Rnew.block[it->second][a1][0].rows() != Mtmp.rows() or
508 Rnew.block[it->second][a1][0].cols() != Mtmp.cols())
509 {
510 Rnew.block[it->second][a1][0] = Mtmp;
511 }
512 else
513 {
514 Rnew.block[it->second][a1][0] += Mtmp;
515 // if(Rnew.block[it->second][a1][0].norm() < ::mynumeric_limits<MpoScalar>::epsilon())
516 // {
517 // Rnew.block[it->second][a1][0].resize(0,0);
518 // continue;
519 // }
520 }
521 }
522 else
523 {
524 boost::multi_array<MatrixType,LEGLIMIT> Mtmpvec(boost::extents[V.at(key).rows()][1]);
525 Mtmpvec[a1][0] = Mtmp;
526 Rnew.push_back(quple, Mtmpvec);
527 }
528 }
529 }
530 }
531 }
532 }
533 }
534 }
535}
536
537template<typename Symmetry, typename MatrixType, typename MatrixType2>
538void contract_L (const Biped<Symmetry,MatrixType2> &Lold,
539 const vector<Biped<Symmetry,MatrixType> > &Abra,
540 const vector<Biped<Symmetry,MatrixType> > &Aket,
541 const vector<qarray<Symmetry::Nq> > &qloc,
542 Biped<Symmetry,MatrixType2> &Lnew,
543 bool RANDOMIZE = false,
544 bool CLEAR = false)
545{
546 if (CLEAR)
547 {
548 Lnew.clear();
549 }
550 else
551 {
552 Lnew.outerResize(Lold);
553 Lnew.setZero();
554 }
555
556 for (size_t s=0; s<qloc.size(); ++s)
557 for (size_t qL=0; qL<Lold.dim; ++qL)
558 {
559 vector<tuple<qarray2<Symmetry::Nq>,size_t,size_t> > ix;
560 bool FOUND_MATCH = LAA(Lold.in[qL], Lold.out[qL], s, qloc, Abra, Aket, ix);
561
562 if (FOUND_MATCH)
563 {
564 for (size_t n=0; n<ix.size(); n++)
565 {
566 qarray2<Symmetry::Nq> quple = get<0>(ix[n]);
567 swap(quple[0], quple[1]);
568 if (!Symmetry::validate(quple)) {continue;}
569 size_t qAbra = get<1>(ix[n]);
570 size_t qAket = get<2>(ix[n]);
571
572 if (Lold.block[qL].rows() != 0)
573 {
574 MatrixType2 Mtmp;
575 if (RANDOMIZE)
576 {
577 Mtmp.resize(Abra[s].block[qAbra].cols(), Aket[s].block[qAket].cols());
578 Mtmp.setRandom();
579 }
580 else
581 {
582// print_size(Abra[s].block[qAbra].adjoint(), "Abra[s].block[qAbra].adjoint()");
583// print_size(Lold.block[qL], "Lold.block[qL]");
584// print_size(Aket[s].block[qAket], "Aket[s].block[qAket]");
585// cout << Abra[s].out[qAbra] << ", " << Abra[s].in[qAbra] << endl;
586// cout << Lold.in[qL] << ", " << Lold.out[qL] << endl;
587// cout << Aket[s].in[qAket] << ", " << Aket[s].out[qAket] << endl;
588// cout << endl;
589
590 optimal_multiply(1.,
591 Abra[s].block[qAbra].adjoint(),
592 Lold.block[qL],
593 Aket[s].block[qAket],
594 Mtmp);
595 }
596
597 auto it = Lnew.dict.find(quple);
598 if (it != Lnew.dict.end())
599 {
600 if (Lnew.block[it->second].rows() != Mtmp.rows() or
601 Lnew.block[it->second].cols() != Mtmp.cols())
602 {
603 Lnew.block[it->second] = Mtmp;
604 }
605 else
606 {
607 Lnew.block[it->second] += Mtmp;
608 }
609 }
610 else
611 {
612 Lnew.push_back(quple, Mtmp);
613 }
614 }
615 }
616 }
617 }
618}
619
620template<typename Symmetry, typename MatrixType, typename MatrixType2>
621void contract_R (const Biped<Symmetry,MatrixType2> &Rold,
622 const vector<Biped<Symmetry,MatrixType> > &Abra,
623 const vector<Biped<Symmetry,MatrixType> > &Aket,
624 const vector<qarray<Symmetry::Nq> > &qloc,
625 Biped<Symmetry,MatrixType2> &Rnew,
626 bool RANDOMIZE = false,
627 bool CLEAR = false)
628{
629 if (CLEAR)
630 {
631 Rnew.clear();
632 }
633 else
634 {
635 Rnew.outerResize(Rold);
636 Rnew.setZero();
637 }
638
639 for (size_t s=0; s<qloc.size(); ++s)
640 for (size_t qR=0; qR<Rold.dim; ++qR)
641 {
642 vector<tuple<qarray2<Symmetry::Nq>,size_t,size_t> > ix;
643 bool FOUND_MATCH = AAR(Rold.in[qR], Rold.out[qR], s, qloc, Abra, Aket, ix);
644
645 if (FOUND_MATCH)
646 {
647 for (size_t n=0; n<ix.size(); n++)
648 {
649 qarray2<Symmetry::Nq> quple = get<0>(ix[n]);
650 swap(quple[0], quple[1]);
651 if (!Symmetry::validate(quple)) {continue;}
652 size_t qAbra = get<1>(ix[n]);
653 size_t qAket = get<2>(ix[n]);
654
655 double factor_cgc = Symmetry::coeff_rightOrtho(Abra[s].out[qAbra],
656 Abra[s].in [qAbra]);
657
658 if (Rold.block[qR].rows() != 0)
659 {
660 MatrixType2 Mtmp;
661 if (RANDOMIZE)
662 {
663 Mtmp.resize(Aket[s].block[qAket].rows(), Abra[s].block[qAbra].rows());
664 Mtmp.setRandom();
665 }
666 else
667 {
668 optimal_multiply(factor_cgc,
669 Aket[s].block[qAket],
670 Rold.block[qR],
671 Abra[s].block[qAbra].adjoint(),
672 Mtmp);
673 }
674
675 auto it = Rnew.dict.find(quple);
676 if (it != Rnew.dict.end())
677 {
678 if (Rnew.block[it->second].rows() != Mtmp.rows() or
679 Rnew.block[it->second].cols() != Mtmp.cols())
680 {
681 Rnew.block[it->second] = Mtmp;
682 }
683 else
684 {
685 Rnew.block[it->second] += Mtmp;
686 }
687 }
688 else
689 {
690 Rnew.push_back(quple, Mtmp);
691 }
692 }
693 }
694 }
695 }
696}
697
698// template<typename Symmetry, typename MatrixType>
699// void contract_TT (const Tripod<Symmetry,MatrixType> &T1,
700// const Tripod<Symmetry,MatrixType> &T2,
701// Biped<Symmetry,MatrixType> &Res)
702// {
703// Res.clear();
704// Res.setZero();
705
706// // for (size_t s=0; s<qloc.size(); ++s)
707// // {
708// for (size_t qT1=0; qT1<T1.dim; ++qT1)
709// {
710// // auto qRouts = Symmetry::reduceSilent(Rold.out[qR],Symmetry::flip(qloc[s]));
711// // auto qRins = Symmetry::reduceSilent(Rold.in[qR],Symmetry::flip(qloc[s]));
712
713// // for(const auto& qRout : qRouts)
714// // for(const auto& qRin : qRins)
715// // {
716// qarray3<Symmetry::Nq> quple = { T1.out(qT1),T1.in(qT1),T1.mid(qT1); }
717// if(auto qT2=T2.dict.find(quple); qT2 != T2.dict.end())
718// {
719// qarray<Symmetry::Nq> new_qin = T1.in(qT1);
720// qarray<Symmetry::Nq> new_qout = T1.out(qT2->second);
721// qarray2<Symmetry::Nq> quple2 = {new_qin, new_qout};
722// if (!Symmetry::validate(quple2)) {continue;}
723
724// double factor_cgc = Symmetry::coeff_rightOrtho(Abra[s].out[q1->second],
725// Abra[s].in[q1->second]);
726
727// if (Rold.block[qR].rows() != 0)
728// {
729// MatrixType Mtmp;
730// optimal_multiply(factor_cgc,
731// Aket[s].block[q2->second],
732// Rold.block[qR],
733// Abra[s].block[q1->second].adjoint(),
734// Mtmp);
735
736// auto it = Rnew.dict.find(quple);
737// if (it != Rnew.dict.end())
738// {
739// if (Rnew.block[it->second].rows() != Mtmp.rows() or
740// Rnew.block[it->second].cols() != Mtmp.cols())
741// {
742// Rnew.block[it->second] = Mtmp;
743// }
744// else
745// {
746// Rnew.block[it->second] += Mtmp;
747// }
748// }
749// else
750// {
751// Rnew.push_back(quple, Mtmp);
752// }
753// }
754// }
755// }
756// // }
757// }
758
759template<typename Symmetry, typename Scalar, typename MpoMatrixType>
760void contract_GRALF (const Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &L,
761 const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Abra,
762 const vector<vector<vector<Biped<Symmetry,MpoMatrixType> > > > &W,
763 const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Aket,
764 const Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &R,
765 const vector<qarray<Symmetry::Nq> > &qloc,
766 const vector<qarray<Symmetry::Nq> > &qOp,
767 Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &Tout,
768 DMRG::DIRECTION::OPTION DIR)
769{
770 #ifdef DMRG_PARALLELIZE_GRALF
771 #pragma omp parallel for collapse(3) schedule(dynamic)
772 #endif
773 for (size_t s1=0; s1<qloc.size(); ++s1)
774 for (size_t s2=0; s2<qloc.size(); ++s2)
775 for (size_t k=0; k<qOp.size(); ++k)
776 {
777 std::array<typename Symmetry::qType,3> qCheck;
778 Scalar factor_cgc, factor_merge;
779
780 qCheck = {qloc[s2],qOp[k],qloc[s1]};
781 if(!Symmetry::validate(qCheck)) {continue;}
782
783 for (size_t qL=0; qL<L.dim; ++qL)
784 {
785 vector<tuple<qarray3<Symmetry::Nq>,size_t,size_t, size_t> > ix;
786 bool FOUND_MATCH = LAWA(L.in(qL), L.out(qL), L.mid(qL), qloc[s1], qloc[s2], qOp[k], Abra[s1], Aket[s2], W[s1][s2][k], ix);
787
788 if (FOUND_MATCH == true)
789 {
790 for(size_t n=0; n<ix.size(); n++ )
791 {
792 qarray3<Symmetry::Nq> quple = get<0>(ix[n]);
793 auto qR = R.dict.find(quple);
794
795 if (qR != R.dict.end())
796 {
797 swap(quple[0], quple[1]);
798 size_t qAbra = get<1>(ix[n]);
799 size_t qAket = get<2>(ix[n]);
800 size_t qW = get<3>(ix[n]);
801 if (Aket[s2].block[qAket].size() == 0) {continue;}
802 if (Abra[s1].block[qAbra].size() == 0) {continue;}
803 if constexpr (Symmetry::NON_ABELIAN)
804 {
805 if (DIR == DMRG::DIRECTION::RIGHT)
806 {
807 factor_cgc = Symmetry::coeff_buildL(Aket[s2].in[qAket], qloc[s2], Aket[s2].out[qAket],
808 L.mid(qL) , qOp[k] , quple[2],
809 Abra[s1].in[qAbra], qloc[s1], Abra[s1].out[qAbra]);
810
811 factor_merge = Symmetry::coeff_rightOrtho(Abra[s1].out[qAbra],
812 Aket[s2].out[qAket]);
813 }
814 else if (DIR == DMRG::DIRECTION::LEFT)
815 {
816 factor_cgc = Symmetry::coeff_buildR(Aket[s2].in[qAket], qloc[s2], Aket[s2].out[qAket],
817 L.mid(qL) , qOp[k] , quple[2],
818 Abra[s1].in[qAbra], qloc[s1], Abra[s1].out[qAbra]);
819 factor_merge = Symmetry::coeff_rightOrtho(L.in(qL),
820 L.out(qL));
821 }
822 }
823 else
824 {
825 factor_cgc = 1.;
826 factor_merge = 1.;
827 }
828 if (std::abs(factor_cgc*factor_merge) < ::mynumeric_limits<Scalar>::epsilon()) {continue;}
829
830 for (int r=0; r<W[s1][s2][k].block[qW].outerSize(); ++r)
831 for (typename MpoMatrixType::InnerIterator iW(W[s1][s2][k].block[qW],r); iW; ++iW)
832 {
833 size_t a1 = iW.row();
834 size_t a2 = iW.col();
835
836 if (L.block[qL][a1][0].size() != 0 and
837 R.block[qR->second][a2][0].size() != 0)
838 {
839 Matrix<Scalar,Dynamic,Dynamic> Mtmp;
840 if (DIR == DMRG::DIRECTION::RIGHT)
841 {
842 // RALF
843 optimal_multiply(iW.value() * factor_cgc * factor_merge,
844 R.block[qR->second][a2][0],
845 Abra[s1].block[qAbra].adjoint(),
846 L.block[qL][a1][0],
847 Aket[s2].block[qAket],
848 Mtmp);
849 }
850 else if (DIR == DMRG::DIRECTION::LEFT)
851 {
852 // GRAL
853 optimal_multiply(iW.value() * factor_cgc * factor_merge,
854 Aket[s2].block[qAket],
855 R.block[qR->second][a2][0],
856 Abra[s1].block[qAbra].adjoint(),
857 L.block[qL][a1][0],
858 Mtmp);
859 }
860
861 qarray2<Symmetry::Nq> qTout;
862 if (DIR == DMRG::DIRECTION::RIGHT)
863 {
864 qTout = {Aket[s2].out[qAket], Aket[s2].out[qAket]};
865 }
866 else if (DIR == DMRG::DIRECTION::LEFT)
867 {
868 qTout = {Aket[s2].in[qAket], Aket[s2].in[qAket]};
869 }
870 #pragma omp critical
871 {
872 auto it = Tout.dict.find(qTout);
873 if (it != Tout.dict.end())
874 {
875 if (Tout.block[it->second].rows() != Mtmp.rows() or
876 Tout.block[it->second].cols() != Mtmp.cols())
877 {
878 Tout.block[it->second] = Mtmp;
879 }
880 else
881 {
882 Tout.block[it->second] += Mtmp;
883 }
884 }
885 else
886 {
887 Tout.push_back(qTout,Mtmp);
888 }
889 }
890 }
891 }
892 }
893 }
894 }
895 }
896 }
897}
898
913template<typename Symmetry, typename Scalar, typename MpoMatrixType>
914Scalar contract_LR (const Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &L,
915 const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Abra,
916 const vector<vector<vector<Biped<Symmetry,MpoMatrixType> > > > &W,
917 const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Aket,
918 const Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &R,
919 const vector<qarray<Symmetry::Nq> > &qloc,
920 const vector<qarray<Symmetry::Nq> > &qOp)
921{
922 Scalar res = 0.;
923 std::array<typename Symmetry::qType,3> qCheck;
924 Scalar factor_cgc;
925
926 for (size_t s1=0; s1<qloc.size(); ++s1)
927 for (size_t s2=0; s2<qloc.size(); ++s2)
928 for (size_t k=0; k<qOp.size(); ++k)
929 {
930 qCheck = {qloc[s2],qOp[k],qloc[s1]};
931 if(!Symmetry::validate(qCheck)) {continue;}
932 for (size_t qL=0; qL<L.dim; ++qL)
933 {
934 vector<tuple<qarray3<Symmetry::Nq>,size_t,size_t,size_t> > ix;
935 bool FOUND_MATCH = LAWA(L.in(qL), L.out(qL), L.mid(qL), qloc[s1], qloc[s2], qOp[k], Abra[s1], Aket[s2], W[s1][s2][k], ix);
936 if (FOUND_MATCH == true)
937 {
938 for(size_t n=0; n<ix.size(); n++ )
939 {
940 qarray3<Symmetry::Nq> quple = get<0>(ix[n]);
941 auto qR = R.dict.find(quple);
942
943 if (qR != R.dict.end())
944 {
945 swap(quple[0], quple[1]);
946 size_t qAbra = get<1>(ix[n]);
947 size_t qAket = get<2>(ix[n]);
948 size_t qW = get<3>(ix[n]);
949 if constexpr ( Symmetry::NON_ABELIAN )
950 {
951 factor_cgc = Symmetry::coeff_buildL(Aket[s2].in[qAket], qloc[s2], Aket[s2].out[qAket],
952 L.mid(qL) , qOp[k] , quple[2],
953 Abra[s1].in[qAbra], qloc[s1], Abra[s1].out[qAbra]);
954 }
955 else
956 {
957 factor_cgc = 1.;
958 }
959 if (std::abs(factor_cgc) < ::mynumeric_limits<double>::epsilon()) { continue; }
960
961 for (int r=0; r<W[s1][s2][k].block[qW].outerSize(); ++r)
962 for (typename MpoMatrixType::InnerIterator iW(W[s1][s2][k].block[qW],r); iW; ++iW)
963 {
964 size_t a1 = iW.row();
965 size_t a2 = iW.col();
966
967 if (L.block[qL][a1][0].rows() != 0 and
968 R.block[qR->second][a2][0].rows() != 0)
969 {
970// Matrix<Scalar,Dynamic,Dynamic> Mtmp = iW.value() *
971// (Abra[s1].block[qAbra].adjoint() *
972// L.block[qL][a1][0] *
973// Aket[s2].block[qAket] *
974// R.block[qR->second][a2][0]);
975// res += Mtmp.trace();
976
977 Matrix<Scalar,Dynamic,Dynamic> Mtmp = L.block[qL][a1][0] *
978 Aket[s2].block[qAket] *
979 R.block[qR->second][a2][0];
980 for (size_t i=0; i<Abra[s1].block[qAbra].cols(); ++i)
981 {
982 res += factor_cgc * iW.value() * Abra[s1].block[qAbra].col(i).dot(Mtmp.col(i));
983 }
984 }
985 }
986 }
987 }
988 }
989 }
990 }
991 return res;
992}
993
994template<typename Symmetry, typename Scalar>
995Scalar contract_LR (pair<qarray<Symmetry::Nq>,size_t> fixed_b,
996 const Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &L,
997 const Biped <Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &R)
998{
999 Scalar res = 0;
1000 assert(fixed_b.first == Symmetry::qvacuum());
1001 for (size_t qL=0; qL<L.dim; ++qL)
1002 {
1003 // Not necessarily vacuum for the structure factor, so this function cannot be used!
1004 assert(L.out(qL) == L.in(qL) and "contract_LR(Tripod,Biped) error!");
1005 if (L.mid(qL) == Symmetry::qvacuum())
1006 {
1007 qarray2<Symmetry::Nq> quple = {L.out(qL), L.in(qL)};
1008 auto qR = R.dict.find(quple);
1009
1010 if (qR != R.dict.end())
1011 {
1012 if (L.block[qL][fixed_b.second][0].size() != 0 and
1013 R.block[qR->second].size() != 0)
1014 {
1015 res += (L.block[qL][fixed_b.second][0] * R.block[qR->second]).trace() * Symmetry::coeff_dot(L.out(qL));
1016 }
1017 }
1018 }
1019 }
1020
1021 return res;
1022}
1023
1024template<typename Symmetry, typename Scalar>
1025Scalar contract_LR (pair<qarray<Symmetry::Nq>,size_t> fixed_a,
1026 const Biped <Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &L,
1027 const Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &R)
1028{
1029 Scalar res = 0;
1030 assert(fixed_a.first == Symmetry::qvacuum());
1031 for (size_t qR=0; qR<R.dim; ++qR)
1032 {
1033 // Not necessarily vacuum for the structure factor, so this function cannot be used!
1034 assert(R.out(qR) == R.in(qR) and "contract_LR(Biped,Tripod) error!");
1035
1036 if (R.mid(qR) == Symmetry::qvacuum())
1037 {
1038 qarray2<Symmetry::Nq> quple = {R.out(qR), R.in(qR)};
1039 auto qL = L.dict.find(quple);
1040
1041 if (qL != L.dict.end())
1042 {
1043 if (R.block[qR][fixed_a.second][0].size() != 0 and
1044 L.block[qL->second].size() != 0)
1045 {
1046// print_size(L.block[qL->second],"L.block[qL->second]");
1047// print_size(R.block[qR][fixed_a.second][0],"R.block[qR][fixed_a.second][0]");
1048 res += (L.block[qL->second] * R.block[qR][fixed_a.second][0]).trace() * Symmetry::coeff_dot(R.out(qR));
1049 }
1050 }
1051 }
1052 }
1053 return res;
1054}
1055
1056template<typename Symmetry, typename Scalar>
1057Scalar contract_LR (const Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &L,
1058 const Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &R)
1059{
1060 Scalar res = 0;
1061
1062 for (size_t qL=0; qL<L.dim; ++qL)
1063 {
1064 qarray3<Symmetry::Nq> quple = {L.out(qL), L.in(qL), L.mid(qL)};
1065 auto qR = R.dict.find(quple);
1066
1067 if (qR != R.dict.end())
1068 {
1069 if (L.block[qL].shape()[0] != R.block[qR->second].shape()[0])
1070 {
1071 cout << "L.block[qL].shape()[0]=" << L.block[qL].shape()[0] << endl;
1072 cout << "R.block[qR->second].shape()[0]=" << R.block[qR->second].shape()[0] << endl;
1073 }
1074 assert(L.block[qL].shape()[0] == R.block[qR->second].shape()[0]);
1075
1076 for (size_t a=0; a<L.block[qL].shape()[0]; ++a)
1077 {
1078 if (L.block[qL][a][0].size() != 0 and
1079 R.block[qR->second][a][0].size() != 0)
1080 {
1081 res += (L.block[qL][a][0] * R.block[qR->second][a][0]).trace() * Symmetry::coeff_dot(L.in(qL));
1082 }
1083 }
1084 }
1085 }
1086
1087 return res;
1088}
1089
1090//template<typename Symmetry, typename MatrixType>
1091//void contract_LR (const Tripod<Symmetry,MatrixType> &L,
1092// const Tripod<Symmetry,MatrixType> &R,
1093// const std::array<qarray<Symmetry::Nq>,D> &qloc,
1094// Tripod<Symmetry,MatrixType> &Bres)
1095//{
1096// Bres.clear();
1097// Bres.setZero();
1098//
1099// for (size_t s1=0; s1<D; ++s1)
1100// for (size_t s2=0; s2<D; ++s2)
1101// for (size_t qL=0; qL<L.dim; ++qL)
1102// {
1103// qarray3<Symmetry::Nq> quple = {L.out(qL), L.in(qL), L.mid(qL)};
1104// auto qR = R.dict.find(quple);
1105//
1106// if (qR != R.dict.end())
1107// {
1108// if (L.block[qL][a1][0].rows() != 0 and
1109// R.block[qR->second][a2][0].rows() != 0)
1110// {
1116//
1117// MatrixType Mtmp = L.block[qL][a1][0] * R.block[qR->second][a2][0]);
1118//
1119// cout << Mtmp.rows() << "\t" << Mtmp.cols() << endl << Mtmp << endl << endl;
1120//
1141// }
1142// }
1143// }
1144//}
1145
1146//template<typename Symmetry, typename MatrixType>
1147//void dryContract_L (const Tripod<Symmetry,MatrixType> &Lold,
1148// const vector<Biped<Symmetry,MatrixType> > &Abra,
1149// const std::array<std::array<Biped<Symmetry,SparseMatrixXd>,D>,D> &W,
1150// const vector<Biped<Symmetry,MatrixType> > &Aket,
1151// const std::array<qarray<Symmetry::Nq>,D> &qloc,
1152// Tripod<Symmetry,MatrixType> &Lnew,
1153// vector<tuple<qarray3<Symmetry::Nq>,std::array<size_t,8> > > &ix)
1154//{
1155// Lnew.setZero();
1156//
1157// MatrixType Mtmp(1,1); Mtmp << 1.;
1158//
1159// for (size_t s1=0; s1<D; ++s1)
1160// for (size_t s2=0; s2<D; ++s2)
1161// for (size_t qL=0; qL<Lold.dim; ++qL)
1162// {
1163// qarray2<Symmetry::Nq> cmp1 = {Lold.in(qL), Lold.in(qL)+qloc[s1]};
1164// qarray2<Symmetry::Nq> cmp2 = {Lold.out(qL), Lold.out(qL)+qloc[s2]};
1165// qarray2<Symmetry::Nq> cmpW = {Lold.mid(qL), Lold.mid(qL)+qloc[s1]-qloc[s2]};
1166//
1167// auto q1 = Abra[s1].dict.find(cmp1);
1168// auto q2 = Aket[s2].dict.find(cmp2);
1169// auto qW = W[s1][s2].dict.find(cmpW);
1170//
1171// if (q1!=Abra[s1].dict.end() and
1172// q2!=Aket[s2].dict.end() and
1173// qW!=W[s1][s2].dict.end())
1174// {
1175// qarray<Symmetry::Nq> new_qin = Abra[s1].out[q1->second]; // A†.in = A.out
1176// qarray<Symmetry::Nq> new_qout = Aket[s2].out[q2->second]; // A.in
1177// qarray<Symmetry::Nq> new_qmid = W[s1][s2].out[qW->second];
1178// qarray3<Symmetry::Nq> quple = {new_qin, new_qout, new_qmid};
1179//
1180// size_t Wcols = W[s1][s2].block[qW->second].cols();
1181//
1182// for (int k=0; k<W[s1][s2].block[qW->second].outerSize(); ++k)
1183// for (SparseMatrixXd::InnerIterator iW(W[s1][s2].block[qW->second],k); iW; ++iW)
1184// {
1185// size_t a1 = iW.row();
1186// size_t a2 = iW.col();
1187//
1188// if (Lold.block[qL][a1][0].rows() != 0)
1189// {
1190// std::array<size_t,9> juple = {s1, s2, q1->second, qW->second, q2->second, qL, a1, a2};
1191// ix.push_back(make_tuple(quple,juple));
1192//
1193// auto it = Lnew.dict.find(quple);
1194// if (it != Lnew.dict.end())
1195// {
1196// if (Lnew.block[it->second][a2][0].rows() == 0)
1197// {
1198// Lnew.block[it->second][a2][0] = Mtmp;
1199// }
1200// }
1201// else
1202// {
1203// boost::multi_array<MatrixType,LEGLIMIT> Mtmpvec(boost::extents[Wcols][1]);
1204// Mtmpvec[a2][0] = Mtmp;
1205// Lnew.push_back(quple, Mtmpvec);
1206// }
1207// }
1208// }
1209// }
1210// }
1211//}
1212
1213//template<typename Symmetry, typename MatrixType>
1214//void contract_L (const Tripod<Symmetry,MatrixType> &Lold,
1215// const vector<tuple<qarray3<Symmetry::Nq>,std::array<size_t,8> > > ix,
1216// const vector<Biped<Symmetry,MatrixType> > &Abra,
1217// const std::array<std::array<Biped<Symmetry,SparseMatrixXd>,D>,D> &W,
1218// const vector<Biped<Symmetry,MatrixType> > &Aket,
1219// const std::array<qarray<Symmetry::Nq>,D> &qloc,
1220// Tripod<Symmetry,MatrixType> &Lnew)
1221//{
1222// Lnew.setZero();
1223//
1224// for (size_t i=0; i<ix.size(); ++i)
1225// {
1226// auto quple = get<0>(ix);
1227// size_t s1 = get<1>(ix)[0];
1228// size_t s2 = get<1>(ix)[1];
1229// size_t q1 = get<1>(ix)[2];
1230// size_t qW = get<1>(ix)[3];
1231// size_t q2 = get<1>(ix)[4];
1232// size_t qL = get<1>(ix)[5];
1233// size_t a1 = get<1>(ix)[6];
1234// size_t a2 = get<1>(ix)[7];
1235//
1236// for (int k=0; k<W[s1][s2].block.block[qW].outerSize(); ++k)
1237// for (SparseMatrixXd::InnerIterator iW(W[s1][s2].block.block[qW],k); iW; ++iW)
1238// {
1239// size_t a1 = iW.row();
1240// size_t a2 = iW.col();
1241//
1242// if (Lold.block[qL][a1][0].rows() != 0)
1243// {
1244// MatrixType Mtmp = iW.value() *
1245// (Abra[s1].block[q1].adjoint() *
1246// Lold.block[qL][a1][0] *
1247// Aket[s2].block[q2]);
1248//
1249// auto it = Lnew.dict.find(quple);
1250// if (it != Lnew.dict.end())
1251// {
1252// if (Lnew.block[it->second][a2][0].rows() != Mtmp.rows() or
1253// Lnew.block[it->second][a2][0].cols() != Mtmp.cols())
1254// {
1255// Lnew.block[it->second][a2][0] = Mtmp;
1256// }
1257// else
1258// {
1259// Lnew.block[it->second][a2][0] += Mtmp;
1260// }
1261// }
1262// else
1263// {
1264// boost::multi_array<MatrixType,LEGLIMIT> Mtmpvec(boost::extents[W[s1][s2].block.block[qW].cols()][1]);
1265// Mtmpvec[a2][0] = Mtmp;
1266// Lnew.push_back(quple, Mtmpvec);
1267// }
1268// }
1269// }
1270// }
1271//}
1272
1278template<typename Symmetry, typename MatrixType, typename MpoMatrixType>
1279void contract_R (const Tripod<Symmetry,MatrixType> &Rold,
1280 const vector<Biped<Symmetry,MatrixType> > &Abra,
1281 const vector<vector<vector<Biped<Symmetry,MpoMatrixType> > > > &Wbot,
1282 const vector<vector<vector<Biped<Symmetry,MpoMatrixType> > > > &Wtop,
1283 const vector<Biped<Symmetry,MatrixType> > &Aket,
1284 const vector<qarray<Symmetry::Nq> > &qloc,
1285 const vector<qarray<Symmetry::Nq> > &qOpBot,
1286 const vector<qarray<Symmetry::Nq> > &qOpTop,
1287 Tripod<Symmetry,MatrixType> &Rnew)
1288{
1289 // cout << baseRightTop << endl << baseLeftTop << endl;
1290 //auto leftTopQs = baseLeftTop.unordered_qs();
1291 //auto leftBotQs = baseLeftBot.unordered_qs();
1292
1293 Qbasis<Symmetry> baseRightBot, baseRightTop, baseLeftBot, baseLeftTop;
1294 baseLeftBot.pullData(Wbot, 0);
1295 baseRightBot.pullData(Wbot, 1);
1296 baseLeftTop.pullData(Wtop, 0);
1297 baseRightTop.pullData(Wtop, 1);
1298
1299 auto TensorBaseRight = baseRightTop.combine(baseRightBot);
1300 auto TensorBaseLeft = baseLeftTop.combine(baseLeftBot);
1301
1302 std::array<typename Symmetry::qType,3> qCheck;
1303
1304 typename MpoMatrixType::Scalar factor_cgc, factor_merge, factor_check;
1305 Rnew.clear();
1306 Rnew.setZero();
1307
1308 for (size_t s1=0; s1<qloc.size(); ++s1)
1309 for (size_t s2=0; s2<qloc.size(); ++s2)
1310 for (size_t s3=0; s3<qloc.size(); ++s3)
1311 for (size_t k1=0; k1<qOpTop.size(); ++k1)
1312 for (size_t k2=0; k2<qOpBot.size(); ++k2)
1313 {
1314 qCheck = {qloc[s3],qOpTop[k1],qloc[s2]};
1315 if(!Symmetry::validate(qCheck)) {continue;}
1316 qCheck = {qloc[s2],qOpBot[k2],qloc[s1]};
1317 if(!Symmetry::validate(qCheck)) {continue;}
1318
1319 auto ks = Symmetry::reduceSilent(qOpTop[k1],qOpBot[k2]);
1320 for(const auto& k : ks)
1321 {
1322 qCheck = {qloc[s3],k,qloc[s1]};
1323 if(!Symmetry::validate(qCheck)) {continue;}
1324
1325 // product in physical space:
1326 factor_check = Symmetry::coeff_MPOprod6(qloc[s1],qOpBot[k2],qloc[s2],
1327 qOpTop[k1],qloc[s3],k);
1328 if (std::abs(factor_check) < ::mynumeric_limits<double>::epsilon()) { continue; }
1329 for (size_t qR=0; qR<Rold.dim; ++qR)
1330 {
1331 auto qRouts = Symmetry::reduceSilent(Rold.out(qR),Symmetry::flip(qloc[s1]));
1332 auto qRins = Symmetry::reduceSilent(Rold.in(qR),Symmetry::flip(qloc[s3]));
1333 auto qrightAuxs = Sym::split<Symmetry>(Rold.mid(qR),baseRightTop.qs(),baseRightBot.qs());
1334
1335 for(const auto& qRout : qRouts)
1336 for(const auto& qRin : qRins)
1337 {
1338 auto q1 = Abra[s1].dict.find({qRout, Rold.out(qR)});
1339 auto q3 = Aket[s3].dict.find({qRin, Rold.in(qR)});
1340 if (q1!=Abra[s1].dict.end() and q3!=Aket[s3].dict.end())
1341 {
1342 auto qRmids = Symmetry::reduceSilent(Rold.mid(qR),Symmetry::flip(k));
1343 for(const auto& new_qmid : qRmids)
1344 {
1345 qarray3<Symmetry::Nq> quple = {Aket[s3].in[q3->second], Abra[s1].in[q1->second], new_qmid};
1346 factor_cgc = Symmetry::coeff_buildR(Aket[s3].in[q3->second], qloc[s3], Aket[s3].out[q3->second],
1347 new_qmid , k , Rold.mid(qR),
1348 Abra[s1].in[q1->second], qloc[s1], Abra[s1].out[q1->second]);
1349
1350 if (std::abs(factor_cgc) < std::abs(::mynumeric_limits<double>::epsilon())) { continue; }
1351 for(const auto& [qrightAux,qrightAuxP] : qrightAuxs)
1352 {
1353 Eigen::Index left2=TensorBaseRight.leftAmount(Rold.mid(qR),{qrightAux, qrightAuxP});
1354 auto qleftAuxs = Symmetry::reduceSilent(qrightAux,Symmetry::flip(qOpTop[k1]));
1355 for(const auto& qleftAux : qleftAuxs)
1356 {
1357 auto qWtop = Wtop[s2][s3][k1].dict.find({qleftAux,qrightAux});
1358 if(qWtop != Wtop[s2][s3][k1].dict.end())
1359 //if(auto it=leftTopQs.find(qleftAux) != leftTopQs.end())
1360 {
1361 auto qleftAuxPs = Symmetry::reduceSilent(qrightAuxP,Symmetry::flip(qOpBot[k2]));
1362 for(const auto& qleftAuxP : qleftAuxPs)
1363 {
1364 auto qWbot = Wbot[s1][s2][k2].dict.find({qleftAuxP,qrightAuxP});
1365 if(qWbot != Wbot[s1][s2][k2].dict.end())
1366 //if(auto it=leftBotQs.find(qleftAuxP) != leftBotQs.end())
1367 {
1368 factor_merge = Symmetry::coeff_MPOprod9(qleftAux,qleftAuxP,new_qmid,
1369 qOpTop[k1],qOpBot[k2],k,
1370 qrightAux,qrightAuxP,Rold.mid(qR));
1371 if (std::abs(factor_merge) < std::abs(::mynumeric_limits<double>::epsilon())) { continue; }
1372 Eigen::Index left1=TensorBaseLeft.leftAmount(new_qmid,{qleftAux, qleftAuxP});
1373 for (int ktop=0; ktop<Wtop[s2][s3][k1].block[qWtop->second].outerSize(); ++ktop)
1374 for (typename MpoMatrixType::InnerIterator iWtop(Wtop[s2][s3][k1].block[qWtop->second],ktop); iWtop; ++iWtop)
1375 for (int kbot=0; kbot<Wbot[s1][s2][k2].block[qWbot->second].outerSize(); ++kbot)
1376 for (typename MpoMatrixType::InnerIterator iWbot(Wbot[s1][s2][k2].block[qWbot->second],kbot); iWbot; ++iWbot)
1377 {
1378 size_t br = iWbot.row();
1379 size_t bc = iWbot.col();
1380 size_t tr = iWtop.row();
1381 size_t tc = iWtop.col();
1382 typename MpoMatrixType::Scalar Wfactor = iWbot.value() * iWtop.value();
1383
1384 // size_t a1 = left1+br*Wtop[s2][s3][k1].block[qWtop->second].rows()+tr;
1385 // size_t a2 = left2+bc*Wtop[s2][s3][k1].block[qWtop->second].cols()+tc;
1386
1387 size_t a1 = left1+tr*Wbot[s1][s2][k2].block[qWbot->second].rows()+br;
1388 size_t a2 = left2+tc*Wbot[s1][s2][k2].block[qWbot->second].cols()+bc;
1389
1390 if (Rold.block[qR][a2][0].rows() != 0)
1391 {
1392 MatrixType Mtmp;
1393 optimal_multiply(factor_check*factor_merge*factor_cgc*Wfactor,
1394 Aket[s3].block[q3->second],
1395 Rold.block[qR][a2][0],
1396 Abra[s1].block[q1->second].adjoint(),
1397 Mtmp);
1398 auto it = Rnew.dict.find(quple);
1399
1400 if (it != Rnew.dict.end())
1401 {
1402 if (Rnew.block[it->second][a1][0].rows() != Mtmp.rows() or
1403 Rnew.block[it->second][a1][0].cols() != Mtmp.cols())
1404 {
1405 Rnew.block[it->second][a1][0] = Mtmp;
1406 }
1407 else
1408 {
1409 Rnew.block[it->second][a1][0] += Mtmp;
1410 }
1411 }
1412 else
1413 {
1414 boost::multi_array<MatrixType,LEGLIMIT> Mtmpvec(boost::extents[TensorBaseLeft.inner_dim(new_qmid)][1]);
1415 Mtmpvec[a1][0] = Mtmp;
1416 Rnew.push_back(quple, Mtmpvec);
1417 }
1418 }
1419 }
1420 }
1421 }
1422 }
1423 }
1424 }
1425 }
1426 }
1427 }
1428 }
1429 }
1430 }
1431}
1432
1438template<typename Symmetry, typename MatrixType, typename MpoMatrixType>
1439void contract_L (const Multipede<4,Symmetry,MatrixType> &Lold,
1440 const vector<Biped<Symmetry,MatrixType> > &Abra,
1441 const vector<vector<vector<Biped<Symmetry,MpoMatrixType> > > > &Wbot,
1442 const vector<vector<vector<Biped<Symmetry,MpoMatrixType> > > > &Wtop,
1443 const vector<Biped<Symmetry,MatrixType> > &Aket,
1444 const vector<qarray<Symmetry::Nq> > &qloc,
1445 const vector<qarray<Symmetry::Nq> > &qOpBot,
1446 const vector<qarray<Symmetry::Nq> > &qOpTop,
1447 Multipede<4,Symmetry,MatrixType> &Lnew)
1448{
1449 std::array<typename Symmetry::qType,3> qCheck;
1450
1451 Lnew.setZero();
1452
1453 for (size_t s1=0; s1<qloc.size(); ++s1)
1454 for (size_t s2=0; s2<qloc.size(); ++s2)
1455 for (size_t s3=0; s3<qloc.size(); ++s3)
1456 for (size_t k1=0; k1<qOpBot.size(); ++k1)
1457 for (size_t k2=0; k2<qOpTop.size(); ++k2)
1458 {
1459 qCheck = {qloc[s2],qOpBot[k1],qloc[s1]};
1460 if(!Symmetry::validate(qCheck)) {continue;}
1461 qCheck = {qloc[s3],qOpTop[k2],qloc[s2]};
1462 if(!Symmetry::validate(qCheck)) {continue;}
1463 for (size_t qL=0; qL<Lold.dim; ++qL)
1464 {
1465 tuple<qarray4<Symmetry::Nq>,size_t,size_t, size_t, size_t> ix;
1466 bool FOUND_MATCH = AWWA(Lold.in(qL), Lold.out(qL), Lold.bot(qL), Lold.top(qL),
1467 qloc[s1], qloc[s2], qloc[s3], qOpBot[k1], qOpTop[k2], Abra[s1], Aket[s3], Wbot[s1][s2][k1], Wtop[s2][s3][k2], ix);
1468 auto quple = get<0>(ix);
1469 swap(quple[0],quple[1]);
1470 size_t qAbra = get<1>(ix);
1471 size_t qAket = get<2>(ix);
1472 size_t qWbot = get<3>(ix);
1473 size_t qWtop = get<4>(ix);
1474
1475 if (FOUND_MATCH == true)
1476 {
1477 for (int kbot=0; kbot<Wbot[s1][s2][k1].block[qWbot].outerSize(); ++kbot)
1478 for (typename MpoMatrixType::InnerIterator iWbot(Wbot[s1][s2][k1].block[qWbot],kbot); iWbot; ++iWbot)
1479 for (int ktop=0; ktop<Wtop[s2][s3][k2].block[qWtop].outerSize(); ++ktop)
1480 for (typename MpoMatrixType::InnerIterator iWtop(Wtop[s2][s3][k2].block[qWtop],ktop); iWtop; ++iWtop)
1481 {
1482 size_t br = iWbot.row();
1483 size_t bc = iWbot.col();
1484 size_t tr = iWtop.row();
1485 size_t tc = iWtop.col();
1486 typename MpoMatrixType::Scalar Wfactor = iWbot.value() * iWtop.value();
1487
1488 if (Lold.block[qL][br][tr].rows() != 0)
1489 {
1490 MatrixType Mtmp;
1491 optimal_multiply(Wfactor,
1492 Abra[s1].block[qAbra].adjoint(),
1493 Lold.block[qL][br][tr],
1494 Aket[s3].block[qAket],
1495 Mtmp);
1496
1497 if (Mtmp.norm() != 0.)
1498 {
1499 auto it = Lnew.dict.find(quple);
1500 if (it != Lnew.dict.end())
1501 {
1502 if (Lnew.block[it->second][bc][tc].rows() != Mtmp.rows() or
1503 Lnew.block[it->second][bc][tc].cols() != Mtmp.cols())
1504 {
1505 Lnew.block[it->second][bc][tc] = Mtmp;
1506 }
1507 else
1508 {
1509 Lnew.block[it->second][bc][tc] += Mtmp;
1510 }
1511 }
1512 else
1513 {
1514 size_t bcols = Wbot[s1][s2][k1].block[qWbot].cols();
1515 size_t tcols = Wtop[s2][s3][k2].block[qWtop].cols();
1516 boost::multi_array<MatrixType,LEGLIMIT> Mtmparray(boost::extents[bcols][tcols]);
1517 Mtmparray[bc][tc] = Mtmp;
1518 Lnew.push_back(quple, Mtmparray);
1519 }
1520 }
1521 }
1522 }
1523 }
1524 }
1525 }
1526}
1527
1529template<typename Symmetry, typename MatrixType, typename MpoMatrixType>
1530void contract_C0 (vector<qarray<Symmetry::Nq> > qloc,
1531 vector<qarray<Symmetry::Nq> > qOp,
1532 const vector<vector<vector<Biped<Symmetry,MpoMatrixType> > > > &W,
1533 const vector<Biped<Symmetry,MatrixType> > &Aket,
1534 vector<Tripod<Symmetry,MatrixType> > &Cnext)
1535{
1536 Cnext.clear();
1537 Cnext.resize(qloc.size());
1538 std::array<typename Symmetry::qType,3> qCheck;
1539
1540 for (size_t s2=0; s2<qloc.size(); ++s2)
1541 {
1542 qarray2<Symmetry::Nq> cmpA = {Symmetry::qvacuum(), Symmetry::qvacuum()+qloc[s2]};
1543 auto qA = Aket[s2].dict.find(cmpA);
1544
1545 if (qA != Aket[s2].dict.end())
1546 {
1547 for (size_t s1=0; s1<qloc.size(); ++s1)
1548 for (size_t k=0; k<qOp.size(); ++k)
1549 {
1550 qCheck = {qloc[s2],qOp[k],qloc[s1]};
1551 if(!Symmetry::validate(qCheck)) {continue;}
1552 for (int r=0; r<W[s1][s2][k].outerSize(); ++r)
1553 for (typename MpoMatrixType::InnerIterator iW(W[s1][s2][k][0],r); iW; ++iW)
1554 {
1555 MatrixType Mtmp = iW.value() * Aket[s2].block[qA->second];
1556
1557 qarray3<Symmetry::Nq> cmpC = {Symmetry::qvacuum(), Aket[s2].out[qA->second], Symmetry::qvacuum()+qloc[s1]-qloc[s2]};
1558 auto qCnext = Cnext[s1].dict.find(cmpC);
1559 if (qCnext != Cnext[s1].dict.end())
1560 {
1561 if (Cnext[s1].block[qCnext->second][iW.col()][0].rows() != Mtmp.rows() or
1562 Cnext[s1].block[qCnext->second][iW.col()][0].cols() != Mtmp.cols())
1563 {
1564 Cnext[s1].block[qCnext->second][iW.col()][0] = Mtmp;
1565 }
1566 else
1567 {
1568 Cnext[s1].block[qCnext->second][iW.col()][0] += Mtmp;
1569 }
1570 }
1571 else
1572 {
1573 boost::multi_array<MatrixType,LEGLIMIT> Mtmpvec(boost::extents[W[s1][s2][k][0].cols()][1]);
1574 Mtmpvec[iW.col()][0] = Mtmp;
1575 Cnext[s1].push_back({Symmetry::qvacuum(), Aket[s2].out[qA->second], Symmetry::qvacuum()+qloc[s1]-qloc[s2]}, Mtmpvec);
1576 }
1577 }
1578 }
1579 }
1580 }
1581}
1582
1585template<typename Symmetry, typename MatrixType, typename MpoMatrixType>
1586void contract_C (vector<qarray<Symmetry::Nq> > qloc,
1587 vector<qarray<Symmetry::Nq> > qOp,
1588 const vector<Biped<Symmetry,MatrixType> > &Abra,
1589 const vector<vector<vector<Biped<Symmetry,MpoMatrixType> > > > &W,
1590 const vector<Biped<Symmetry,MatrixType> > &Aket,
1591 const vector<Tripod<Symmetry,MatrixType> > &C,
1592 vector<Tripod<Symmetry,MatrixType> > &Cnext)
1593{
1594 Cnext.clear();
1595 Cnext.resize(qloc.size());
1596 std::array<typename Symmetry::qType,3> qCheck;
1597
1598 for (size_t s=0; s<qloc.size(); ++s)
1599 for (size_t qC=0; qC<C[s].dim; ++qC)
1600 {
1601 qarray2<Symmetry::Nq> cmpU = {C[s].in(qC), C[s].in(qC)+qloc[s]};
1602 auto qU = Abra[s].dict.find(cmpU);
1603
1604 if (qU != Abra[s].dict.end())
1605 {
1606 for (size_t s1=0; s1<qloc.size(); ++s1)
1607 for (size_t s2=0; s2<qloc.size(); ++s2)
1608 for (size_t k=0; k<qOp.size(); ++k)
1609 {
1610 qCheck = {qloc[s2],qOp[k],qloc[s1]};
1611 if(!Symmetry::validate(qCheck)) {continue;}
1612
1613 qarray2<Symmetry::Nq> cmpA = {C[s].out(qC), C[s].out(qC)+qloc[s2]};
1614 auto qA = Aket[s2].dict.find(cmpA);
1615
1616 if (qA != Aket[s2].dict.end())
1617 {
1618 for (int r=0; r<W[s1][s2][k][0].outerSize(); ++r)
1619 for (typename MpoMatrixType::InnerIterator iW(W[s1][s2][k][0],r); iW; ++iW)
1620 {
1621 if (C[s].block[qC][iW.row()][0].rows() != 0)
1622 {
1623// MatrixType Mtmp = iW.value() * (Abra[s].block[qU->second].adjoint() *
1624// C[s].block[qC][iW.row()][0] *
1625// Aket[s2].block[qA->second]);
1626 MatrixType Mtmp;
1627 optimal_multiply(iW.value(),
1628 Abra[s].block[qU->second].adjoint(),
1629 C[s].block[qC][iW.row()][0],
1630 Aket[s2].block[qA->second],
1631 Mtmp);
1632
1633 qarray3<Symmetry::Nq> cmpC = {Abra[s].out[qU->second], Aket[s2].out[qA->second], C[s].mid(qC)+qloc[s1]-qloc[s2]};
1634 auto qCnext = Cnext[s1].dict.find(cmpC);
1635 if (qCnext != Cnext[s1].dict.end())
1636 {
1637 if (Cnext[s1].block[qCnext->second][iW.col()][0].rows() != Mtmp.rows() or
1638 Cnext[s1].block[qCnext->second][iW.col()][0].cols() != Mtmp.cols())
1639 {
1640 Cnext[s1].block[qCnext->second][iW.col()][0] = Mtmp;
1641 }
1642 else
1643 {
1644 Cnext[s1].block[qCnext->second][iW.col()][0] += Mtmp;
1645 }
1646 }
1647 else
1648 {
1649 boost::multi_array<MatrixType,LEGLIMIT> Mtmpvec(boost::extents[W[s1][s2][0][0].cols()][1]);
1650 Mtmpvec[iW.col()][0] = Mtmp;
1651 Cnext[s1].push_back({Abra[s].out[qU->second], Aket[s2].out[qA->second], C[s].mid(qC)+qloc[s1]-qloc[s2]}, Mtmpvec);
1652 }
1653 }
1654 }
1655 }
1656 }
1657 }
1658 }
1659}
1660
1661//template<typename Symmetry, typename Scalar>
1662//void contract_AA (const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A1,
1663// vector<qarray<Symmetry::Nq> > qloc1,
1664// const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A2,
1665// vector<qarray<Symmetry::Nq> > qloc2,
1666// vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Apair)
1667//{
1668// auto tensor_basis = Symmetry::tensorProd(qloc1,qloc2);
1669// Apair.resize(tensor_basis.size());
1670//
1671// for (size_t s1=0; s1<qloc1.size(); ++s1)
1672// for (size_t s2=0; s2<qloc2.size(); ++s2)
1673// {
1674// auto qmerges = Symmetry::reduceSilent(qloc1[s1], qloc2[s2]);
1675//
1676// for (const auto &qmerge:qmerges)
1677// {
1678// auto qtensor = make_tuple(qloc1[s1], s1, qloc2[s2], s2, qmerge);
1679// auto s1s2 = distance(tensor_basis.begin(), find(tensor_basis.begin(), tensor_basis.end(), qtensor));
1680//
1681// for (size_t q1=0; q1<A1[s1].dim; ++q1)
1682// {
1683// auto qmids = Symmetry::reduceSilent(A1[s1].out[q1], qloc2[s2]);
1684//
1685// for (const auto &qmid:qmids)
1686// {
1687// qarray2<Symmetry::Nq> quple = {A1[s1].out[q1], qmid};
1688// auto q2 = A2[s2].dict.find(quple);
1689//
1690// if (q2 != A2[s2].dict.end())
1691// {
1692// Scalar factor_cgc = Symmetry::coeff_Apair(A1[s1].in[q1], qloc1[s1], A1[s1].out[q1],
1693// qloc2[s2], A2[s2].out[q2->second], qmerge);
1694//
1695// if (abs(factor_cgc) > abs(mynumeric_limits<double>::epsilon()))
1696// {
1697// Matrix<Scalar,Dynamic,Dynamic> Mtmp = factor_cgc * A1[s1].block[q1] * A2[s2].block[q2->second];
1698//
1699// qarray2<Symmetry::Nq> qupleApair = {A1[s1].in[q1], A2[s2].out[q2->second]};
1700//
1701// auto qApair = Apair[s1s2].dict.find(qupleApair);
1702//
1703// if (qApair != Apair[s1s2].dict.end())
1704// {
1705// Apair[s1s2].block[qApair->second] += Mtmp;
1706// }
1707// else
1708// {
1709// Apair[s1s2].push_back(qupleApair, Mtmp);
1710// }
1711// }
1712// }
1713// }
1714// }
1715// }
1716// }
1717//}
1718
1719// FORCE_QTOT to create only one final block;
1720// otherwise crashes when using the result for further calculations (e.g. ground-state sweeping)
1721template<typename Symmetry, typename Scalar, typename MpoMatrixType>
1722void contract_AW (const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Ain,
1723 const vector<qarray<Symmetry::Nq> > &qloc,
1724 const vector<vector<vector<Biped<Symmetry,MpoMatrixType> > > > &W,
1725 const vector<qarray<Symmetry::Nq> > &qOp,
1726 const Qbasis<Symmetry> &qauxAl,
1727 const Qbasis<Symmetry> &qauxWl,
1728 const Qbasis<Symmetry> &qauxAr,
1729 const Qbasis<Symmetry> &qauxWr,
1730 vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Aout,
1731 bool FORCE_QTOT=false,
1732 qarray<Symmetry::Nq> Qtot=Symmetry::qvacuum())
1733{
1734 for (size_t s=0; s<qloc.size(); ++s)
1735 {
1736 Aout[s].clear();
1737 }
1738
1739 double factor_cgc;
1740
1741 auto tensorBasis_l = qauxAl.combine(qauxWl);
1742 auto tensorBasis_r = qauxAr.combine(qauxWr);
1743
1744 for (size_t s1=0; s1<qloc.size(); ++s1)
1745 for (size_t s2=0; s2<qloc.size(); ++s2)
1746 for (size_t k=0; k<qOp.size(); ++k)
1747 {
1748 qarray3<Symmetry::Nq> qCheck_ = {qloc[s2],qOp[k],qloc[s1]};
1749 if (!Symmetry::validate(qCheck_)) { continue; }
1750 // if(W[s1][s2][k].size() == 0) { continue; } //Checks whether QNs s1, s2 and k fit together.
1751
1752 for (size_t q=0; q<Ain[s2].size(); q++)
1753 // for (const auto &[qWl,qWl_dim,qWl_plain] : qauxWl) // cpp is too stupid to call cbegin() and cend() here...
1754 for (auto it=qauxWl.cbegin(); it != qauxWl.cend(); it++)
1755 {
1756 auto [qWl, qWl_dim, qWl_plain] = *it;
1757 auto qWrs = Symmetry::reduceSilent(qWl,qOp[k]);
1758
1759 for (const auto &qWr : qWrs)
1760 {
1761 if (qauxWr.find(qWr) == false) {continue;}
1762
1763 auto qmerge_ls = Symmetry::reduceSilent(Ain[s2].in[q] ,qWl);
1764 auto qmerge_rs = Symmetry::reduceSilent(Ain[s2].out[q],qWr);
1765
1766 for (const auto qmerge_l : qmerge_ls)
1767 for (const auto qmerge_r : qmerge_rs)
1768 {
1769 qarray3<Symmetry::Nq> qCheck = {qmerge_l,qloc[s1],qmerge_r};
1770
1771 if (!Symmetry::validate(qCheck)) { continue; }
1772 if (tensorBasis_l.find(qmerge_l) == false) {continue;}
1773 if (tensorBasis_r.find(qmerge_r) == false) {continue;}
1774
1775 if constexpr (Symmetry::NON_ABELIAN)
1776 {
1777 factor_cgc = Symmetry::coeff_AW(Ain[s2].in[q], qloc[s2], Ain[s2].out[q],
1778 qWl , qOp[k] , qWr,
1779 qmerge_l , qloc[s1], qmerge_r);
1780 }
1781 else { factor_cgc = 1.; }
1782
1783 if (abs(factor_cgc) < abs(mynumeric_limits<double>::epsilon())) { continue; }
1784
1785 Matrix<Scalar,Dynamic,Dynamic> Mtmp(tensorBasis_l.inner_dim(qmerge_l), tensorBasis_r.inner_dim(qmerge_r)); Mtmp.setZero();
1786 int left_l = tensorBasis_l.leftAmount(qmerge_l, { Ain[s2].in[q] , qWl } );
1787 int left_r = tensorBasis_r.leftAmount(qmerge_r, { Ain[s2].out[q], qWr } );
1788
1789 auto dict_entry = W[s1][s2][k].dict.find({qWl,qWr});
1790 if(dict_entry == W[s1][s2][k].dict.end()) continue;
1791 size_t qW = dict_entry->second;
1792
1793 for (int r=0; r<W[s1][s2][k].block[qW].outerSize(); ++r)
1794 for (typename MpoMatrixType::InnerIterator iW(W[s1][s2][k].block[qW],r); iW; ++iW)
1795 {
1796 size_t wr = iW.row() * Ain[s2].block[q].rows();
1797 size_t wc = iW.col() * Ain[s2].block[q].cols();
1798 Mtmp.block(wr+left_l,wc+left_r,Ain[s2].block[q].rows(),Ain[s2].block[q].cols()) += Ain[s2].block[q] * iW.value() * factor_cgc;
1799 }
1800
1801 qarray2<Symmetry::Nq> cmp = qarray2<Symmetry::Nq>{qmerge_l, qmerge_r}; //auxiliary quantum numbers of Aout
1802 auto it = Aout[s1].dict.find(cmp);
1803
1804 if (!FORCE_QTOT)
1805 {
1806 if (it == Aout[s1].dict.end())
1807 {
1808 Aout[s1].push_back(cmp,Mtmp);
1809 }
1810 else
1811 {
1812 Aout[s1].block[it->second] += Mtmp;
1813 }
1814 }
1815 else
1816 {
1817 if (it == Aout[s1].dict.end() and cmp[1] == Qtot)
1818 {
1819 Aout[s1].push_back(cmp,Mtmp);
1820 }
1821 else if (it != Aout[s1].dict.end() and cmp[1] == Qtot)
1822 {
1823 Aout[s1].block[it->second] += Mtmp;
1824 }
1825 }
1826 }
1827 }
1828 }
1829 }
1830}
1831
1832template<typename Symmetry, typename Scalar>
1833vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > extend_A (const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A, const vector<qarray<Symmetry::Nq> > &qloc)
1834{
1835 vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > Aout(A.size());
1836 Qbasis<Symmetry> in; in.pullData(A,0);
1837 Qbasis<Symmetry> out; out.pullData(A,1);
1838 Qbasis<Symmetry> loc; loc.pullData(qloc);
1839 for (size_t s=0; s<qloc.size(); s++)
1840 for (const auto & [qin,num,plain] : in)
1841 {
1842 auto qouts = Symmetry::reduceSilent(qin,qloc[s]);
1843 for (const auto & qout:qouts)
1844 {
1845 if (!out.find(qout)) {continue;}
1846 qarray2<Symmetry::Nq> cmp{qin,qout};
1847 auto it_A = A[s].dict.find(cmp);
1848 if (it_A != A[s].dict.end())
1849 {
1850 Aout[s].push_back(qin,qout,A[s].block[it_A->second]);
1851 }
1852 else
1853 {
1854 Matrix<Scalar,Dynamic,Dynamic> Mtmp(in.inner_dim(qin), out.inner_dim(qout)); Mtmp.setZero();
1855 Aout[s].push_back(qin,qout,Mtmp);
1856 }
1857 }
1858 }
1859 return Aout;
1860}
1861
1862template<typename Symmetry, typename Scalar>
1863void extend_A (vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A, const vector<qarray<Symmetry::Nq> > &qloc)
1864{
1865 const auto copy = A;
1866 A = extend_A(copy,qloc);
1867}
1868
1869template<typename Symmetry, typename Scalar>
1870void contract_AA2 (const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A1,
1871 const vector<qarray<Symmetry::Nq> > &qloc1,
1872 const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A2,
1873 const vector<qarray<Symmetry::Nq> > &qloc2,
1874 vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Apair,
1875 bool DRY = false)
1876{
1877 Qbasis<Symmetry> locBasis1; locBasis1.pullData(qloc1);
1878 Qbasis<Symmetry> locBasis2; locBasis2.pullData(qloc2);
1879 auto locBasis = locBasis1.combine(locBasis2);
1880
1881 Apair.resize(locBasis.size());
1882 for (size_t s1=0; s1<qloc1.size(); s1++)
1883 for (size_t s2=0; s2<qloc2.size(); s2++)
1884 {
1885 auto q_s1s2s = Symmetry::reduceSilent(qloc1[s1],qloc2[s2]);
1886 for (const auto &q_s1s2: q_s1s2s)
1887 {
1888 //get the number of the basis state with QN q_s1s2 which results from state (s1 with QN qloc1[s1] of locBasis1) and (s2 with QN qloc2[s2] of locBasis2)
1889 size_t s1s2 = locBasis.outer_num(q_s1s2) + locBasis.leftOffset(q_s1s2,{qloc1[s1],qloc2[s2]},{locBasis1.inner_num(s1),locBasis2.inner_num(s2)});
1890 // size_t s1s2 = locBasis.outer_num(q_s1s2) + locBasis.leftAmount(q_s1s2,{qloc1[s1],qloc2[s2]}) + locBasis1.inner_num(s1) + locBasis2.inner_num(s2)*locBasis1.inner_dim(qloc1[s1]);
1891 for (size_t q1=0; q1<A1[s1].size(); q1++)
1892 {
1893 typename Symmetry::qType qm = A1[s1].out[q1];
1894 auto q2s = Symmetry::reduceSilent(qm,qloc2[s2]);
1895 for (const auto &q2 : q2s)
1896 {
1897 auto it_q2 = A2[s2].dict.find({qm,q2});
1898 if ( it_q2 == A2[s2].dict.end()) {continue;}
1899 Eigen::Matrix<Scalar,-1,-1> Mtmp(A1[s1].block[q1].rows(),A2[s2].block[it_q2->second].cols());
1900 Mtmp.setZero();
1901 Scalar factor_cgc = Symmetry::coeff_twoSiteGate(A1[s1].in[q1], qloc1[s1], qm,
1902 qloc2[s2] , q2 , q_s1s2);
1903 if (abs(factor_cgc) < ::mynumeric_limits<double>::epsilon()) {continue;}
1904 // cout << "ql=" << A1[s1].in[q1] << ", s1=" << qloc1[s1] << ", qm=" << qm << ", s2=" << qloc2[s2] << ", q2=" << q2 << ", s1s2=" << q_s1s2 << endl;
1905 // cout << "factor_cgc=" << factor_cgc << endl;
1906 // cout << "l=" << l << ", s1p=" << s1p << ", ql=" << ql << ", s2p=" << s2p << "qr=" << it_qr->second << endl;
1907 // print_size(A[l][s1p].block[ql], "A[l][s1p].block[ql]");
1908 // print_size(A[l+1][s2p].block[it_qr->second], "A[l+1][s2p].block[it_qr->second]");
1909 if (!DRY)
1910 {
1911 Mtmp = factor_cgc * A1[s1].block[q1] * A2[s2].block[it_q2->second];
1912 }
1913 // cout << ", qmid[k]=" << q_s1s2 << ", s1s2=" << s1s2 << ", q_s1s2=" << locBasis.find(s1s2) << endl;
1914 auto it_pair = Apair[s1s2].dict.find({A1[s1].in[q1],q2});
1915 if (it_pair == Apair[s1s2].dict.end())
1916 {
1917 Apair[s1s2].push_back(A1[s1].in[q1],q2,Mtmp);
1918 }
1919 else
1920 {
1921 Apair[s1s2].block[it_pair->second] += Mtmp;
1922 }
1923 }
1924 }
1925 }
1926 }
1927}
1928
1929template<typename Symmetry, typename Scalar>
1930void contract_AA (const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A1,
1931 const vector<qarray<Symmetry::Nq> > &qloc1,
1932 const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A2,
1933 const vector<qarray<Symmetry::Nq> > &qloc2,
1934 const qarray<Symmetry::Nq> &Qtop,
1935 const qarray<Symmetry::Nq> &Qbot,
1936 vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Apair,
1937 bool DRY = false)
1938{
1939 auto tensor_basis = Symmetry::tensorProd(qloc1,qloc2);
1940 Apair.resize(tensor_basis.size());
1941
1942 vector<qarray<Symmetry::Nq> > qsplit = calc_qsplit (A1, qloc1, A2, qloc2, Qtop, Qbot);
1943
1944 vector<qarray<Symmetry::Nq> > A1in;
1945 vector<qarray<Symmetry::Nq> > A2out;
1946
1947 // gather all qin at the left:
1948 for (size_t s1=0; s1<qloc1.size(); ++s1)
1949 for (size_t q=0; q<A1[s1].dim; ++q)
1950 {
1951 A1in.push_back(A1[s1].in[q]);
1952 }
1953 // gather all qout at the right:
1954 for (size_t s2=0; s2<qloc2.size(); ++s2)
1955 for (size_t q=0; q<A2[s2].dim; ++q)
1956 {
1957 A2out.push_back(A2[s2].out[q]);
1958 }
1959
1960 for (size_t s1=0; s1<qloc1.size(); ++s1)
1961 for (size_t m=0; m<qsplit.size(); ++m)
1962 {
1963 auto qins = Symmetry::reduceSilent(qsplit[m], Symmetry::flip(qloc1[s1]));
1964
1965 for (const auto &qin:qins)
1966 {
1967 for (size_t s2=0; s2<qloc2.size(); ++s2)
1968 {
1969 auto qmerges = Symmetry::reduceSilent(qloc1[s1], qloc2[s2]);
1970
1971 for (const auto &qmerge:qmerges)
1972 {
1973 auto qtensor = make_tuple(qloc1[s1], s1, qloc2[s2], s2, qmerge);
1974 auto s1s2 = distance(tensor_basis.begin(), find(tensor_basis.begin(), tensor_basis.end(), qtensor));
1975
1976 auto qouts = Symmetry::reduceSilent(qsplit[m], qloc2[s2]);
1977 for (const auto &qout:qouts)
1978 {
1979 auto qA1 = find(A1in.begin(), A1in.end(), qin);
1980 auto qA2 = find(A2out.begin(), A2out.end(), qout);
1981 if (qA1 != A1in.end() and qA2 != A2out.end())
1982 {
1983 Apair[s1s2].try_create_block({qin,qout});
1984 }
1985 }
1986 }
1987 }
1988 }
1989 }
1990
1991 #ifdef DMRG_CONTRACTAA_PARALLELIZE
1992 #pragma omp parallel for collapse(2) schedule(dynamic)
1993 #endif
1994 for (size_t s1=0; s1<qloc1.size(); ++s1)
1995 for (size_t s2=0; s2<qloc2.size(); ++s2)
1996 {
1997 auto qmerges = Symmetry::reduceSilent(qloc1[s1], qloc2[s2]);
1998
1999 for (const auto &qmerge:qmerges)
2000 {
2001 auto qtensor = make_tuple(qloc1[s1], s1, qloc2[s2], s2, qmerge);
2002 auto s1s2 = distance(tensor_basis.begin(), find(tensor_basis.begin(), tensor_basis.end(), qtensor));
2003
2004 for (size_t q1=0; q1<A1[s1].dim; ++q1)
2005 {
2006 auto qouts = Symmetry::reduceSilent(A1[s1].out[q1], qloc2[s2]);
2007
2008 for (const auto &qout:qouts)
2009 {
2010 qarray2<Symmetry::Nq> quple = {A1[s1].out[q1], qout};
2011 auto q2 = A2[s2].dict.find(quple);
2012
2013 if (q2 != A2[s2].dict.end())
2014 {
2015 Scalar factor_cgc = Symmetry::coeff_Apair(A1[s1].in[q1], qloc1[s1], A1[s1].out[q1],
2016 qloc2[s2], A2[s2].out[q2->second], qmerge);
2017 if (abs(factor_cgc) > abs(mynumeric_limits<double>::epsilon()))
2018 {
2019 Matrix<Scalar,Dynamic,Dynamic> Mtmp;
2020 if (!DRY)
2021 {
2022 Mtmp = factor_cgc * A1[s1].block[q1] * A2[s2].block[q2->second];
2023 }
2024
2025 if (Mtmp.size() != 0)
2026 {
2027 #pragma omp critical
2028 {
2029 qarray2<Symmetry::Nq> qupleApair = {A1[s1].in[q1], A2[s2].out[q2->second]};
2030
2031 auto qApair = Apair[s1s2].dict.find(qupleApair);
2032
2033 if (qApair != Apair[s1s2].dict.end() and
2034 Apair[s1s2].block[qApair->second].size() == Mtmp.size())
2035 {
2036 Apair[s1s2].block[qApair->second] += Mtmp;
2037 }
2038 else if (qApair != Apair[s1s2].dict.end() and
2039 Apair[s1s2].block[qApair->second].size() == 0)
2040 {
2041 Apair[s1s2].block[qApair->second] = Mtmp;
2042 }
2043 else
2044 {
2045 Apair[s1s2].push_back(qupleApair, Mtmp);
2046 }
2047 }
2048 }
2049 }
2050 }
2051 }
2052 }
2053 }
2054 }
2055}
2056
2058template<typename Symmetry, typename Scalar>
2059void contract_AAAA (const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A1,
2060 vector<qarray<Symmetry::Nq> > qloc1,
2061 const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A2,
2062 vector<qarray<Symmetry::Nq> > qloc2,
2063 const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A3,
2064 vector<qarray<Symmetry::Nq> > qloc3,
2065 const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &A4,
2066 vector<qarray<Symmetry::Nq> > qloc4,
2067 boost::multi_array<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> >,4> &Aquartett)
2068{
2069 Aquartett.resize(boost::extents[qloc1.size()][qloc2.size()][qloc3.size()][qloc4.size()]);
2070
2071 for (size_t s1=0; s1<qloc1.size(); ++s1)
2072 for (size_t s2=0; s2<qloc2.size(); ++s2)
2073 for (size_t s3=0; s3<qloc3.size(); ++s3)
2074 for (size_t s4=0; s4<qloc4.size(); ++s4)
2075 for (size_t q1=0; q1<A1[s1].dim; ++q1)
2076 {
2077 qarray2<Symmetry::Nq> quple2 = {A1[s1].out[q1], A1[s1].out[q1]+qloc2[s2]};
2078 auto q2 = A2[s2].dict.find(quple2);
2079
2080 if (q2 != A2[s2].dict.end())
2081 {
2082 qarray2<Symmetry::Nq> quple3 = {A2[s2].out[q2->second], A2[s2].out[q2->second]+qloc3[s3]};
2083 auto q3 = A3[s3].dict.find(quple3);
2084
2085 if (q3 != A3[s3].dict.end())
2086 {
2087 qarray2<Symmetry::Nq> quple4 = {A3[s3].out[q3->second], A3[s3].out[q3->second]+qloc4[s4]};
2088 auto q4 = A4[s4].dict.find(quple4);
2089
2090 if (q4 != A4[s4].dict.end())
2091 {
2092 Matrix<Scalar,Dynamic,Dynamic> Mtmp = A1[s1].block[q1] *
2093 A2[s2].block[q2->second] *
2094 A3[s3].block[q3->second] *
2095 A4[s4].block[q4->second];
2096
2097 qarray2<Symmetry::Nq> qupleAquartett = {A1[s1].in[q1], A4[s4].out[q4->second]};
2098 auto qAquartett = Aquartett[s1][s2][s3][s4].dict.find(qupleAquartett);
2099
2100 if (qAquartett != Aquartett[s1][s2][s3][s4].dict.end())
2101 {
2102 Aquartett[s1][s2][s3][s4].block[qAquartett->second] += Mtmp;
2103 }
2104 else
2105 {
2106 Aquartett[s1][s2][s3][s4].push_back(qupleAquartett, Mtmp);
2107 }
2108 }
2109 }
2110 }
2111 }
2112}
2113
2114template<typename Symmetry, typename Scalar>
2115void split_AA (DMRG::DIRECTION::OPTION DIR, const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Apair,
2116 const vector<qarray<Symmetry::Nq> >& qloc_l, vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Al,
2117 const vector<qarray<Symmetry::Nq> >& qloc_r, vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Ar,
2118 const qarray<Symmetry::Nq>& qtop, const qarray<Symmetry::Nq>& qbot, double eps_svd, size_t min_Nsv, size_t max_Nsv)
2119{
2120 Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > Cdump;
2121 double truncDump, Sdump;
2122 split_AA(DIR, Apair, qloc_l, Al, qloc_r, Ar, qtop, qbot,
2123 Cdump, false, truncDump, Sdump, eps_svd,min_Nsv,max_Nsv);
2124}
2125
2126template<typename Symmetry, typename Scalar>
2127void split_AA (DMRG::DIRECTION::OPTION DIR, const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Apair,
2128 const vector<qarray<Symmetry::Nq> >& qloc_l, vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Al,
2129 const vector<qarray<Symmetry::Nq> >& qloc_r, vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Ar,
2130 const qarray<Symmetry::Nq>& qtop, const qarray<Symmetry::Nq>& qbot,
2131 Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &C, bool SEPARATE_SV, double &truncWeight, double &S, double eps_svd, size_t min_Nsv, size_t max_Nsv)
2132{
2133 vector<qarray<Symmetry::Nq> > midset = calc_qsplit(Al, qloc_l, Ar, qloc_r, qtop, qbot);
2134
2135 for (size_t s=0; s<qloc_l.size(); ++s)
2136 {
2137 Al[s].clear();
2138 }
2139 for (size_t s=0; s<qloc_r.size(); ++s)
2140 {
2141 Ar[s].clear();
2142 }
2143
2144 ArrayXd truncWeightSub(midset.size()); truncWeightSub.setZero();
2145 ArrayXd entropySub(midset.size()); entropySub.setZero();
2146
2147 auto tensor_basis = Symmetry::tensorProd(qloc_l, qloc_r);
2148
2149 #ifdef DMRG_SPLITAA_PARALLELIZE
2150 #pragma omp parallel for schedule(dynamic)
2151 #endif
2152 for (size_t qmid=0; qmid<midset.size(); ++qmid)
2153 {
2154 map<pair<size_t,qarray<Symmetry::Nq> >,vector<pair<size_t,qarray<Symmetry::Nq> > > > s13map;
2155 map<tuple<size_t,qarray<Symmetry::Nq>,size_t,qarray<Symmetry::Nq> >,vector<Scalar> > cgcmap;
2156 map<tuple<size_t,qarray<Symmetry::Nq>,size_t,qarray<Symmetry::Nq> >,vector<size_t> > q13map;
2157 map<tuple<size_t,qarray<Symmetry::Nq>,size_t,qarray<Symmetry::Nq> >,vector<size_t> > s1s3map;
2158
2159 for (size_t s1=0; s1<qloc_l.size(); ++s1)
2160 for (size_t s3=0; s3<qloc_r.size(); ++s3)
2161 {
2162 auto qmerges = Symmetry::reduceSilent(qloc_l[s1], qloc_r[s3]);
2163
2164 for (const auto &qmerge:qmerges)
2165 {
2166 auto qtensor = make_tuple(qloc_l[s1], s1, qloc_r[s3], s3, qmerge);
2167 auto s1s3 = distance(tensor_basis.begin(), find(tensor_basis.begin(), tensor_basis.end(), qtensor));
2168
2169 for (size_t q13=0; q13<Apair[s1s3].dim; ++q13)
2170 {
2171 auto qlmids = Symmetry::reduceSilent(Apair[s1s3].in[q13], qloc_l[s1]);
2172 auto qrmids = Symmetry::reduceSilent(Apair[s1s3].out[q13], Symmetry::flip(qloc_r[s3]));
2173
2174 for (const auto &qlmid:qlmids)
2175 for (const auto &qrmid:qrmids)
2176 {
2177 if (qlmid == midset[qmid] and qrmid == midset[qmid])
2178 {
2179 s13map[make_pair(s1,Apair[s1s3].in[q13])].push_back(make_pair(s3,Apair[s1s3].out[q13]));
2180
2181 Scalar factor_cgc = Symmetry::coeff_Apair(Apair[s1s3].in[q13], qloc_l[s1], midset[qmid],
2182 qloc_r[s3], Apair[s1s3].out[q13], qmerge);
2183 if (DIR==DMRG::DIRECTION::LEFT)
2184 {
2185 factor_cgc *= Symmetry::coeff_leftSweep(Apair[s1s3].out[q13],midset[qmid]);
2186 }
2187
2188 cgcmap[make_tuple(s1,Apair[s1s3].in[q13],s3,Apair[s1s3].out[q13])].push_back(factor_cgc);
2189 q13map[make_tuple(s1,Apair[s1s3].in[q13],s3,Apair[s1s3].out[q13])].push_back(q13);
2190 s1s3map[make_tuple(s1,Apair[s1s3].in[q13],s3,Apair[s1s3].out[q13])].push_back(s1s3);
2191 }
2192 }
2193 }
2194 }
2195 }
2196
2197 if (s13map.size() != 0)
2198 {
2199 map<pair<size_t,qarray<Symmetry::Nq> >,Matrix<Scalar,Dynamic,Dynamic> > Aclumpvec;
2200 size_t istitch = 0;
2201 size_t jstitch = 0;
2202 vector<size_t> get_s3;
2203 vector<size_t> get_Ncols;
2204 vector<qarray<Symmetry::Nq> > get_qr;
2205 bool COLS_ARE_KNOWN = false;
2206
2207 for (size_t s1=0; s1<qloc_l.size(); ++s1)
2208 {
2209 auto qls = Symmetry::reduceSilent(midset[qmid], Symmetry::flip(qloc_l[s1]));
2210
2211 for (const auto &ql:qls)
2212 {
2213 for (size_t s3=0; s3<qloc_r.size(); ++s3)
2214 {
2215 auto qrs = Symmetry::reduceSilent(midset[qmid], qloc_r[s3]);
2216
2217 for (const auto &qr:qrs)
2218 {
2219 auto s3block = find(s13map[make_pair(s1,ql)].begin(), s13map[make_pair(s1,ql)].end(), make_pair(s3,qr));
2220
2221 if (s3block != s13map[make_pair(s1,ql)].end())
2222 {
2223 Matrix<Scalar,Dynamic,Dynamic> Mtmp;
2224 for (size_t i=0; i<q13map[make_tuple(s1,ql,s3,qr)].size(); ++i)
2225 {
2226 size_t q13 = q13map[make_tuple(s1,ql,s3,qr)][i];
2227 size_t s1s3 = s1s3map[make_tuple(s1,ql,s3,qr)][i];
2228
2229 if (Mtmp.size() == 0)
2230 {
2231 Mtmp = cgcmap[make_tuple(s1,ql,s3,qr)][i] * Apair[s1s3].block[q13];
2232 }
2233 else if (Mtmp.size() > 0 and Apair[s1s3].block[q13].size() > 0)
2234 {
2235 Mtmp += cgcmap[make_tuple(s1,ql,s3,qr)][i] * Apair[s1s3].block[q13];
2236 }
2237 }
2238 if (Mtmp.size() == 0) {continue;}
2239
2240 addRight(Mtmp, Aclumpvec[make_pair(s1,ql)]);
2241
2242 if (COLS_ARE_KNOWN == false)
2243 {
2244 get_s3.push_back(s3);
2245 get_Ncols.push_back(Mtmp.cols());
2246 get_qr.push_back(qr);
2247 }
2248 }
2249 }
2250 }
2251 if (get_s3.size() != 0) {COLS_ARE_KNOWN = true;}
2252 }
2253 }
2254
2255 vector<size_t> get_s1;
2256 vector<size_t> get_Nrows;
2257 vector<qarray<Symmetry::Nq> > get_ql;
2258 Matrix<Scalar,Dynamic,Dynamic> Aclump;
2259 for (size_t s1=0; s1<qloc_l.size(); ++s1)
2260 {
2261 auto qls = Symmetry::reduceSilent(midset[qmid], Symmetry::flip(qloc_l[s1]));
2262
2263 for (const auto &ql:qls)
2264 {
2265 size_t Aclump_rows_old = Aclump.rows();
2266
2267 // If cols don't match, it means that zeros were cut, restore them
2268 // (happens in MpsCompressor::polyCompress):
2269 if (Aclumpvec[make_pair(s1,ql)].cols() < Aclump.cols())
2270 {
2271 size_t dcols = Aclump.cols() - Aclumpvec[make_pair(s1,ql)].cols();
2272 Aclumpvec[make_pair(s1,ql)].conservativeResize(Aclumpvec[make_pair(s1,ql)].rows(), Aclump.cols());
2273 Aclumpvec[make_pair(s1,ql)].rightCols(dcols).setZero();
2274 }
2275 else if (Aclumpvec[make_pair(s1,ql)].cols() > Aclump.cols())
2276 {
2277 size_t dcols = Aclumpvec[make_pair(s1,ql)].cols() - Aclump.cols();
2278 Aclump.conservativeResize(Aclump.rows(), Aclump.cols()+dcols);
2279 Aclump.rightCols(dcols).setZero();
2280 }
2281
2282 addBottom(Aclumpvec[make_pair(s1,ql)], Aclump);
2283
2284 if (Aclump.rows() > Aclump_rows_old)
2285 {
2286 get_s1.push_back(s1);
2287 get_Nrows.push_back(Aclump.rows()-Aclump_rows_old);
2288 get_ql.push_back(ql);
2289 }
2290 }
2291 }
2292 if (Aclump.size() == 0)
2293 {
2294// if (DIR == DMRG::DIRECTION::RIGHT)
2295// {
2296// this->pivot = (loc==this->N_sites-1)? this->N_sites-1 : loc+1;
2297// }
2298// else
2299// {
2300// this->pivot = (loc==0)? 0 : loc;
2301// }
2302 continue;
2303 }
2304
2305 #ifdef DONT_USE_BDCSVD
2306 JacobiSVD<Matrix<Scalar,Dynamic,Dynamic> > Jack; // standard SVD
2307 #else
2308 BDCSVD<Matrix<Scalar,Dynamic,Dynamic> > Jack; // "Divide and conquer" SVD (only available in Eigen)
2309 #endif
2310 Jack.compute(Aclump,ComputeThinU|ComputeThinV);
2311 VectorXd SV = Jack.singularValues();
2312
2313 // retained states:
2314 size_t Nret = (SV.array().abs() > eps_svd).count();
2315 Nret = max(Nret, min_Nsv);
2316 Nret = min(Nret, max_Nsv);
2317 truncWeightSub(qmid) = Symmetry::degeneracy(midset[qmid]) * SV.tail(SV.rows()-Nret).cwiseAbs2().sum();
2318 size_t Nnz = (Jack.singularValues().array() > 0.).count();
2319 entropySub(qmid) = -Symmetry::degeneracy(midset[qmid]) *
2320 (SV.head(Nnz).array().square() * SV.head(Nnz).array().square().log()).sum();
2321
2322 Matrix<Scalar,Dynamic,Dynamic> Aleft, Aright, Cmatrix;
2323 if (DIR == DMRG::DIRECTION::RIGHT)
2324 {
2325 Aleft = Jack.matrixU().leftCols(Nret);
2326 if (SEPARATE_SV)
2327 {
2328 Aright = Jack.matrixV().adjoint().topRows(Nret);
2329 Cmatrix = Jack.singularValues().head(Nret).asDiagonal();
2330 }
2331 else
2332 {
2333 Aright = Jack.singularValues().head(Nret).asDiagonal() * Jack.matrixV().adjoint().topRows(Nret);
2334 }
2335// this->pivot = (loc==this->N_sites-1)? this->N_sites-1 : loc+1;
2336 }
2337 else
2338 {
2339 Aright = Jack.matrixV().adjoint().topRows(Nret);
2340 if (SEPARATE_SV)
2341 {
2342 Aleft = Jack.matrixU().leftCols(Nret);
2343 Cmatrix = Jack.singularValues().head(Nret).asDiagonal();
2344 }
2345 else
2346 {
2347 Aleft = Jack.matrixU().leftCols(Nret) * Jack.singularValues().head(Nret).asDiagonal();
2348 }
2349// this->pivot = (loc==0)? 0 : loc;
2350 }
2351
2352 // update Al
2353 #pragma omp critical
2354 {
2355 istitch = 0;
2356 for (size_t i=0; i<get_s1.size(); ++i)
2357 {
2358 size_t s1 = get_s1[i];
2359 size_t Nrows = get_Nrows[i];
2360
2361 qarray2<Symmetry::Nq> quple = {get_ql[i], midset[qmid]};
2362 auto q = Al[s1].dict.find(quple);
2363 if (q != Al[s1].dict.end())
2364 {
2365 Al[s1].block[q->second] += Aleft.block(istitch,0, Nrows,Nret);
2366 }
2367 else
2368 {
2369 Al[s1].push_back(get_ql[i], midset[qmid], Aleft.block(istitch,0, Nrows,Nret));
2370 }
2371 istitch += Nrows;
2372 }
2373
2374 // update Ar
2375 jstitch = 0;
2376 for (size_t i=0; i<get_s3.size(); ++i)
2377 {
2378 size_t s3 = get_s3[i];
2379 size_t Ncols = get_Ncols[i];
2380
2381 qarray2<Symmetry::Nq> quple = {midset[qmid], get_qr[i]};
2382 auto q = Ar[s3].dict.find(quple);
2383 Scalar factor_cgc3 = (DIR==DMRG::DIRECTION::LEFT)? Symmetry::coeff_leftSweep(midset[qmid],
2384 get_qr[i]):1.;
2385 if (q != Ar[s3].dict.end())
2386 {
2387 Ar[s3].block[q->second] += factor_cgc3 * Aright.block(0,jstitch, Nret,Ncols);
2388 }
2389 else
2390 {
2391 Ar[s3].push_back(midset[qmid], get_qr[i], factor_cgc3 * Aright.block(0,jstitch, Nret,Ncols));
2392 }
2393 jstitch += Ncols;
2394 }
2395
2396 if (SEPARATE_SV)
2397 {
2398 qarray2<Symmetry::Nq> quple = {midset[qmid], midset[qmid]};
2399 auto q = C.dict.find(quple);
2400 if (q != C.dict.end())
2401 {
2402 C.block[q->second] += Cmatrix;
2403 }
2404 else
2405 {
2406 C.push_back(midset[qmid], midset[qmid], Cmatrix);
2407 }
2408 }
2409 }
2410 }
2411 }
2412
2413 // remove unwanted zero-sized blocks
2414 for (size_t s=0; s<qloc_l.size(); ++s)
2415 {
2416 Al[s] = Al[s].cleaned();
2417 }
2418 for (size_t s=0; s<qloc_r.size(); ++s)
2419 {
2420 Ar[s] = Ar[s].cleaned();
2421 }
2422
2423 truncWeight = truncWeightSub.sum();
2424 S = entropySub.sum();
2425 // if (DIR == DMRG::DIRECTION::RIGHT)
2426 // {
2427 // int bond = (loc==this->N_sites-1)? -1 : loc;
2428 // if (bond != -1)
2429 // {
2430 // S(loc) = entropySub.sum();
2431 // }
2432 // }
2433 // else
2434 // {
2435 // int bond = (loc==0)? -1 : loc;
2436 // if (bond != -1)
2437 // {
2438 // S(loc-1) = entropySub.sum();
2439 // }
2440 // }
2441}
2442
2443template<typename Symmetry, typename Scalar>
2444void split_AA2 (DMRG::DIRECTION::OPTION DIR, const Qbasis<Symmetry>& locBasis, const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Apair,
2445 const vector<qarray<Symmetry::Nq> >& qloc_l, vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Al,
2446 const vector<qarray<Symmetry::Nq> >& qloc_r, vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Ar,
2447 const qarray<Symmetry::Nq>& qtop, const qarray<Symmetry::Nq>& qbot, double eps_svd, size_t min_Nsv, size_t max_Nsv)
2448{
2449 Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > Cdump;
2450 double truncDump, Sdump;
2451 map<qarray<Symmetry::Nq>,ArrayXd> SVspec_dumb;
2452 split_AA2(DIR, locBasis, Apair, qloc_l, Al, qloc_r, Ar, qtop, qbot,
2453 Cdump, false, truncDump, Sdump, SVspec_dumb, eps_svd,min_Nsv,max_Nsv);
2454}
2455
2456// M
2457// --ooooo-- = --o---o--
2458// | | | |
2459//M is the possibly enlarged basis
2460template<typename Symmetry, typename Scalar>
2461void split_AA2 (DMRG::DIRECTION::OPTION DIR, const Qbasis<Symmetry>& locBasis, const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Apair,
2462 const vector<qarray<Symmetry::Nq> >& qloc_l, vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Al,
2463 const vector<qarray<Symmetry::Nq> >& qloc_r, vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Ar,
2464 const qarray<Symmetry::Nq>& qtop, const qarray<Symmetry::Nq>& qbot,
2465 Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &C, bool SEPARATE_SV, double &truncWeight, double &S, map<qarray<Symmetry::Nq>,ArrayXd> &SVspec,
2466 double eps_truncWeight, size_t min_Nsv, size_t max_Nsv)
2467{
2468 truncWeight=0.;
2469 S=0;
2470
2471 Qbasis<Symmetry> leftBasis; leftBasis.pullData(Al,0);
2472 Qbasis<Symmetry> rightBasis; rightBasis.pullData(Ar,1);
2473 // cout << "left x right:" << endl << locBasis << endl;
2474 Qbasis<Symmetry> locBasis_l; locBasis_l.pullData(qloc_l);
2475 // cout << "left:" << endl << locBasis_l << endl;
2476 Qbasis<Symmetry> locBasis_r; locBasis_r.pullData(qloc_r);
2477 // cout << "right" << endl << locBasis_r << endl;
2478 Qbasis<Symmetry> leftTot = leftBasis.combine(locBasis_l);
2479 // cout << "combined Basis left:" << endl << leftTot << endl;
2480 Qbasis<Symmetry> rightTot = rightBasis.combine(locBasis_r, true); //true means flip locBasis_r before the combination
2481 // cout << "combined Basis right:" << endl << rightTot << endl;
2482
2483 for (size_t s=0; s<qloc_l.size(); ++s)
2484 {
2485 Al[s].clear();
2486 }
2487 for (size_t s=0; s<qloc_r.size(); ++s)
2488 {
2489 Ar[s].clear();
2490 }
2491
2492 Biped<Symmetry,Eigen::Matrix<Scalar,-1,-1> > Aclump;
2493 for (size_t sl=0; sl<qloc_l.size(); sl++)
2494 for (size_t sr=0; sr<qloc_r.size(); sr++)
2495 {
2496 auto q_slsrs = Symmetry::reduceSilent(qloc_l[sl],qloc_r[sr]);
2497 for (const auto &q_slsr:q_slsrs)
2498 {
2499 //get the number of the basis state with QN q_slsr which results from state (sl with QN qloc_l[sl] of locBasis_l) and (sr with QN qloc_r[sr] of locBasis_r)
2500 size_t s = locBasis.outer_num(q_slsr) + locBasis.leftOffset(q_slsr,{qloc_l[sl],qloc_r[sr]},{locBasis_l.inner_num(sl),locBasis_r.inner_num(sr)});
2501 // size_t s = locBasis.outer_num(q_slsr) + locBasis.leftAmount(q_slsr,{qloc_l[sl],qloc_r[sr]}) + locBasis_l.inner_num(sl) + locBasis_r.inner_num(sr)*locBasis_l.inner_dim(qloc_l[sl]);
2502 assert(locBasis.find(s) == q_slsr);
2503 // if (!Symmetry::triangle(qarray3<Symmetry::Nq>{qloc_l[sl],qloc_r[sr],qloc[s]})) {continue;}
2504 for (size_t q=0; q<Apair[s].size(); q++)
2505 {
2506 auto Qls = Symmetry::reduceSilent(Apair[s].in[q],qloc_l[sl]);
2507 auto Qrs = Symmetry::reduceSilent(Apair[s].out[q],Symmetry::flip(qloc_r[sr]));
2508 for (const auto &Ql:Qls)
2509 for (const auto &Qr:Qrs)
2510 {
2511 if (Ql != Qr) {continue;} //Aclump is a matrix so Qin=Qout is mandatory. here Qin=Ql and Qout=Qr
2512 // Scalar factor_cgc = Symmetry::coeff_splitAA(qloc_l[sl] , qloc_r[sr] , q_slsr,
2513 // Apair[s].out[q], Apair[s].in[q], Ql);
2514 Scalar factor_cgc = Symmetry::coeff_splitAA(Apair[s].in[q], Apair[s].out[q], q_slsr,
2515 qloc_r[sr] , qloc_l[sl] , Ql);
2516 if (abs(factor_cgc) < ::mynumeric_limits<double>::epsilon()) {continue;}
2517 // cout << "ql=" << Apair[s].in[q] << ", qr=" << Apair[s].out[q] << ", slsr=" << q_slsr << ", sr=" << qloc_r[sr] << ", sl=" << qloc_l[sl] << ", Ql=" << Ql << endl;
2518 // cout << "factor_cgc=" << factor_cgc << endl;
2519 Eigen::Matrix<Scalar,-1,-1> Mtmp(leftTot.inner_dim(Ql),rightTot.inner_dim(Qr));
2520 Mtmp.setZero();
2521 // Index plain_left1 = leftBasis.inner_dim(Apair[s].in[q])*locBasis_l.inner_num(sl);
2522 // Index left1 = leftTot.leftAmount(Ql,{Apair[s].in[q], qloc_l[sl]}) + plain_left1;
2523 Index left1 = leftTot.leftOffset(Ql, {Apair[s].in[q], qloc_l[sl]}, {0, locBasis_l.inner_num(sl)});
2524
2525 array<qarray<Symmetry::Nq>,2> source = {Apair[s].out[q], Symmetry::flip(qloc_r[sr])};
2526 // Index plain_left2 = rightBasis.inner_dim(Apair[s].out[q])*locBasis_r.inner_num(sr);
2527 // Index left2 = rightTot.leftAmount(Qr,source) + plain_left2;
2528 Index left2 = rightTot.leftOffset(Qr, source, {0, locBasis_r.inner_num(sr)});
2529 Mtmp.block(left1, left2, Apair[s].block[q].rows(), Apair[s].block[q].cols()) += factor_cgc*Apair[s].block[q];
2530 auto it_Aclump = Aclump.dict.find({Ql,Qr});
2531 if (it_Aclump == Aclump.dict.end())
2532 {
2533 Aclump.push_back(Ql,Qr,Mtmp);
2534 }
2535 else
2536 {
2537 Aclump.block[it_Aclump->second] += Mtmp;
2538 }
2539 }
2540 }
2541 }
2542 }
2543 // cout << "Aclump:" << endl << Aclump.print(true) << endl;
2544// cout << "begin truncate SVD" << endl;
2545 auto [U,Sigma,Vdag] = Aclump.truncateSVD(min_Nsv,max_Nsv,eps_truncWeight,truncWeight,S,SVspec,false,false); //false: DONT PRESERVE MULTIPLETS
2546// cout << "end truncate SVD" << endl;
2547 // cout << "U,Sigma,Vdag:" << endl << U.print(true) << endl << Sigma.print(true) << endl << Vdag.print(true) << endl;
2548 Biped<Symmetry,Eigen::Matrix<Scalar,-1,-1> > left,right;
2549 if (DIR == DMRG::DIRECTION::RIGHT)
2550 {
2551 left = U;
2552 if (SEPARATE_SV)
2553 {
2554 right = Vdag;
2555 C = Sigma;
2556 }
2557 else
2558 {
2559 right = Sigma.contract(Vdag);
2560 }
2561 }
2562 else
2563 {
2564 right = Vdag;
2565 if (SEPARATE_SV)
2566 {
2567 left = U;
2568 C = Sigma;
2569 }
2570 else
2571 {
2572 left = U.contract(Sigma);
2573 }
2574 }
2575
2576 //reshape left back and write result to Al
2577 for (const auto &[ql, dim_l, plain] : leftBasis)
2578 for (size_t s1=0; s1<qloc_l.size(); s1++)
2579 {
2580 auto Qls = Symmetry::reduceSilent(ql,qloc_l[s1]);
2581 for (const auto &Ql:Qls)
2582 {
2583 auto it_left = left.dict.find({Ql,Ql});
2584 if (it_left == left.dict.end()) {continue;}
2585 if (Ql > qtop or Ql < qbot) {continue;}
2586 Eigen::Matrix<Scalar,-1,-1> Mtmp(leftBasis.inner_dim(ql),left.block[it_left->second].cols());
2587 // Index left1 = leftTot.leftAmount (Ql,{ql, qloc_l[s1]}) + leftBasis.inner_dim(ql)*locBasis_l.inner_num(s1);
2588 Index left1 = leftTot.leftOffset (Ql,{ql, qloc_l[s1]}, {0, locBasis_l.inner_num(s1)});
2589 Mtmp = left.block[it_left->second].block(left1, 0, leftBasis.inner_dim(ql), left.block[it_left->second].cols());
2590 auto it_A = Al[s1].dict.find({ql,Ql});
2591 if (it_A == Al[s1].dict.end())
2592 {
2593 Al[s1].push_back(ql,Ql,Mtmp);
2594 }
2595 else
2596 {
2597 Al[s1].block[it_A->second] += Mtmp;
2598 }
2599 }
2600 }
2601
2602 //reshape right back and write result to Al
2603 for (const auto &[qr, dim_r, plain] : rightBasis)
2604 for (size_t s2=0; s2<qloc_r.size(); s2++)
2605 {
2606 auto Qrs = Symmetry::reduceSilent(qr,Symmetry::flip(qloc_r[s2]));
2607 for (const auto &Qr:Qrs)
2608 {
2609 auto it_right = right.dict.find({Qr,Qr});
2610 if (it_right == right.dict.end()) {continue;}
2611 if (Qr > qtop or Qr < qbot) {continue;}
2612 Eigen::Matrix<Scalar,-1,-1> Mtmp(right.block[it_right->second].rows(),rightBasis.inner_dim(qr));
2613 array<qarray<Symmetry::Nq>,2> source = {qr, Symmetry::flip(qloc_r[s2])};
2614 // Index left2 = rightTot.leftAmount (Qr,source) + rightBasis.inner_dim(qr)*locBasis_r.inner_num(s2);
2615 Index left2 = rightTot.leftOffset (Qr, source, {0, locBasis_r.inner_num(s2)});
2616 Mtmp = Symmetry::coeff_splitAA(Qr,qr,qloc_r[s2])*right.block[it_right->second].block(0, left2, right.block[it_right->second].rows(), rightBasis.inner_dim(qr));
2617 // cout << "Qr=" << Qr << ", qr=" << qr << ", sr=" << qloc_r[s2] << endl;
2618 // cout << "factor_cgc=" << Symmetry::coeff_splitAA(Qr,qr,qloc_r[s2]) << endl;
2619 auto it_A = Ar[s2].dict.find({Qr,qr});
2620 if (it_A == Al[s2].dict.end())
2621 {
2622 Ar[s2].push_back(Qr,qr,Mtmp);
2623 }
2624 else
2625 {
2626 Ar[s2].block[it_A->second] += Mtmp;
2627 }
2628 }
2629 }
2630}
2631#endif
addBottom
void addBottom(const MatrixType1 &Min, MatrixType2 &Mout)
Definition: DmrgConglutinations.h:83
addRight
void addRight(const MatrixType1 &Min, MatrixType2 &Mout)
Definition: DmrgConglutinations.h:29
contract_AAAA
void contract_AAAA(const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &A1, vector< qarray< Symmetry::Nq > > qloc1, const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &A2, vector< qarray< Symmetry::Nq > > qloc2, const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &A3, vector< qarray< Symmetry::Nq > > qloc3, const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &A4, vector< qarray< Symmetry::Nq > > qloc4, boost::multi_array< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > >, 4 > &Aquartett)
Definition: DmrgContractions.h:2059
contract_AA
void contract_AA(const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &A1, const vector< qarray< Symmetry::Nq > > &qloc1, const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &A2, const vector< qarray< Symmetry::Nq > > &qloc2, const qarray< Symmetry::Nq > &Qtop, const qarray< Symmetry::Nq > &Qbot, vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Apair, bool DRY=false)
Definition: DmrgContractions.h:1930
contract_AA2
void contract_AA2(const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &A1, const vector< qarray< Symmetry::Nq > > &qloc1, const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &A2, const vector< qarray< Symmetry::Nq > > &qloc2, vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Apair, bool DRY=false)
Definition: DmrgContractions.h:1870
contract_C0
void contract_C0(vector< qarray< Symmetry::Nq > > qloc, vector< qarray< Symmetry::Nq > > qOp, const vector< vector< vector< Biped< Symmetry, MpoMatrixType > > > > &W, const vector< Biped< Symmetry, MatrixType > > &Aket, vector< Tripod< Symmetry, MatrixType > > &Cnext)
Definition: DmrgContractions.h:1530
extend_A
vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > extend_A(const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &A, const vector< qarray< Symmetry::Nq > > &qloc)
Definition: DmrgContractions.h:1833
contract_GRALF
void contract_GRALF(const Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > &L, const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Abra, const vector< vector< vector< Biped< Symmetry, MpoMatrixType > > > > &W, const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Aket, const Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > &R, const vector< qarray< Symmetry::Nq > > &qloc, const vector< qarray< Symmetry::Nq > > &qOp, Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > &Tout, DMRG::DIRECTION::OPTION DIR)
Definition: DmrgContractions.h:760
contract_AW
void contract_AW(const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Ain, const vector< qarray< Symmetry::Nq > > &qloc, const vector< vector< vector< Biped< Symmetry, MpoMatrixType > > > > &W, const vector< qarray< Symmetry::Nq > > &qOp, const Qbasis< Symmetry > &qauxAl, const Qbasis< Symmetry > &qauxWl, const Qbasis< Symmetry > &qauxAr, const Qbasis< Symmetry > &qauxWr, vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Aout, bool FORCE_QTOT=false, qarray< Symmetry::Nq > Qtot=Symmetry::qvacuum())
Definition: DmrgContractions.h:1722
split_AA2
void split_AA2(DMRG::DIRECTION::OPTION DIR, const Qbasis< Symmetry > &locBasis, const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Apair, const vector< qarray< Symmetry::Nq > > &qloc_l, vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Al, const vector< qarray< Symmetry::Nq > > &qloc_r, vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Ar, const qarray< Symmetry::Nq > &qtop, const qarray< Symmetry::Nq > &qbot, double eps_svd, size_t min_Nsv, size_t max_Nsv)
Definition: DmrgContractions.h:2444
contract_C
void contract_C(vector< qarray< Symmetry::Nq > > qloc, vector< qarray< Symmetry::Nq > > qOp, const vector< Biped< Symmetry, MatrixType > > &Abra, const vector< vector< vector< Biped< Symmetry, MpoMatrixType > > > > &W, const vector< Biped< Symmetry, MatrixType > > &Aket, const vector< Tripod< Symmetry, MatrixType > > &C, vector< Tripod< Symmetry, MatrixType > > &Cnext)
Definition: DmrgContractions.h:1586
split_AA
void split_AA(DMRG::DIRECTION::OPTION DIR, const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Apair, const vector< qarray< Symmetry::Nq > > &qloc_l, vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Al, const vector< qarray< Symmetry::Nq > > &qloc_r, vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Ar, const qarray< Symmetry::Nq > &qtop, const qarray< Symmetry::Nq > &qbot, double eps_svd, size_t min_Nsv, size_t max_Nsv)
Definition: DmrgContractions.h:2115
CONTRACT_LR_MODE
CONTRACT_LR_MODE
Definition: DmrgContractions.h:13
FIXED_ROWS
@ FIXED_ROWS
Definition: DmrgContractions.h:13
FIXED_COLS
@ FIXED_COLS
Definition: DmrgContractions.h:13
TRIANGULAR
@ TRIANGULAR
Definition: DmrgContractions.h:13
FULL
@ FULL
Definition: DmrgContractions.h:13
optimal_multiply
void optimal_multiply(Scalar alpha, const MatrixTypeA &A, const MatrixTypeB &B, const MatrixTypeC &C, MatrixTypeR &result, bool DEBUG=false)
Definition: DmrgExternal.h:191
AWWA
bool AWWA(qarray< Symmetry::Nq > Lin, qarray< Symmetry::Nq > Lout, qarray< Symmetry::Nq > Lbot, qarray< Symmetry::Nq > Ltop, qarray< Symmetry::Nq > qloc1, qarray< Symmetry::Nq > qloc2, qarray< Symmetry::Nq > qloc3, qarray< Symmetry::Nq > qOpBot, qarray< Symmetry::Nq > qOpTop, const Biped< Symmetry, MatrixType > &Abra, const Biped< Symmetry, MatrixType > &Aket, const Biped< Symmetry, MpoMatrixType > &Wbot, const Biped< Symmetry, MpoMatrixType > &Wtop, tuple< qarray4< Symmetry::Nq >, size_t, size_t, size_t, size_t > &result)
Definition: DmrgIndexGymnastics.h:387
LAWA
bool LAWA(const qarray< Symmetry::Nq > &Lin, const qarray< Symmetry::Nq > &Lout, const qarray< Symmetry::Nq > &Lmid, const qarray< Symmetry::Nq > &qloc1, const qarray< Symmetry::Nq > &qloc2, const qarray< Symmetry::Nq > &qOp, const Biped< Symmetry, MatrixType > &Abra, const Biped< Symmetry, MatrixType > &Aket, const Biped< Symmetry, MatrixType2 > &W, vector< tuple< qarray3< Symmetry::Nq >, size_t, size_t, size_t > > &result)
Definition: DmrgIndexGymnastics.h:33
LAA
bool LAA(qarray< Symmetry::Nq > Lin, qarray< Symmetry::Nq > Lout, size_t s, vector< qarray< Symmetry::Nq > > qloc, const vector< Biped< Symmetry, MatrixType > > &Abra, const vector< Biped< Symmetry, MatrixType > > &Aket, vector< tuple< qarray2< Symmetry::Nq >, size_t, size_t > > &result)
Definition: DmrgIndexGymnastics.h:127
calc_qsplit
vector< qarray< Symmetry::Nq > > calc_qsplit(const vector< Biped< Symmetry, Eigen::Matrix< Scalar, Dynamic, Dynamic > > > &A1, const vector< qarray< Symmetry::Nq > > &qloc1, const vector< Biped< Symmetry, Eigen::Matrix< Scalar, Dynamic, Dynamic > > > &A2, vector< qarray< Symmetry::Nq > > qloc2, const qarray< Symmetry::Nq > &Qtop, const qarray< Symmetry::Nq > &Qbot)
Definition: DmrgIndexGymnastics.h:265
AAR
bool AAR(qarray< Symmetry::Nq > Rin, qarray< Symmetry::Nq > Rout, size_t s, vector< qarray< Symmetry::Nq > > qloc, const vector< Biped< Symmetry, MatrixType > > &Abra, const vector< Biped< Symmetry, MatrixType > > &Aket, vector< tuple< qarray2< Symmetry::Nq >, size_t, size_t > > &result)
Definition: DmrgIndexGymnastics.h:163
AWAR
bool AWAR(const qarray< Symmetry::Nq > &Rin, const qarray< Symmetry::Nq > &Rout, const qarray< Symmetry::Nq > &Rmid, const qarray< Symmetry::Nq > &qloc1, const qarray< Symmetry::Nq > &qloc2, const qarray< Symmetry::Nq > &qOp, const Biped< Symmetry, MatrixType > &Abra, const Biped< Symmetry, MatrixType > &Aket, const Biped< Symmetry, MatrixType2 > &W, vector< tuple< qarray3< Symmetry::Nq >, size_t, size_t, size_t > > &result)
Definition: DmrgIndexGymnastics.h:80
sum
Hamiltonian sum(const Hamiltonian &H1, const Hamiltonian &H2, DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: DmrgLinearAlgebra.h:866
swap
void swap(Mps< Symmetry, Scalar > &V1, Mps< Symmetry, Scalar > &V2)
Definition: DmrgLinearAlgebra.h:65
A
@ A
Definition: DmrgTypedefs.h:130
Qbasis
Definition: Qbasis.h:39
Qbasis::qs
const std::vector< qType > qs() const
Definition: Qbasis.h:279
Qbasis::find
qType find(const std::string &ident) const
Definition: Qbasis.h:297
Qbasis::inner_num
Eigen::Index inner_num(const Eigen::Index &outer_num) const
Definition: Qbasis.h:335
Qbasis::combine
Qbasis< Symmetry > combine(const Qbasis< Symmetry > &other, bool FLIP=false) const
Definition: Qbasis.h:686
Qbasis::cbegin
std::vector< std::tuple< qType, Eigen::Index, Basis > >::const_iterator cbegin() const
Definition: Qbasis.h:189
Qbasis::leftOffset
Eigen::Index leftOffset(const qType &qnew, const std::array< qType, 2 > &qold, const std::array< Eigen::Index, 2 > &plain_old) const
Definition: Qbasis.h:440
Qbasis::pullData
void pullData(const vector< Biped< Symmetry, MatrixType > > &A, const Eigen::Index &leg)
Definition: Qbasis.h:486
Qbasis::cend
std::vector< std::tuple< qType, Eigen::Index, Basis > >::const_iterator cend() const
Definition: Qbasis.h:190
Qbasis::outer_num
Eigen::Index outer_num(const qType &q) const
Definition: Qbasis.h:378
Qbasis::inner_dim
Eigen::Index inner_dim(const Eigen::Index &num_in) const
Definition: Qbasis.h:391
contract_R
void contract_R(const Tripod< Symmetry, MatrixType2 > &Rold, const vector< Biped< Symmetry, MatrixType > > &Abra, const vector< vector< vector< Biped< Symmetry, MpoMatrixType > > > > &W, const vector< Biped< Symmetry, MatrixType > > &Aket, const vector< qarray< Symmetry::Nq > > &qloc, const vector< qarray< Symmetry::Nq > > &qOp, Tripod< Symmetry, MatrixType2 > &Rnew, bool RANDOMIZE=false, tuple< CONTRACT_LR_MODE, size_t > MODE_input=make_pair(FULL, 0), const std::unordered_map< pair< qarray< Symmetry::Nq >, size_t >, size_t > &basis_order_map={})
Definition: DmrgContractions.h:180
contract_L
void contract_L(const Tripod< Symmetry, MatrixType2 > &Lold, const vector< Biped< Symmetry, MatrixType > > &Abra, const vector< vector< vector< Biped< Symmetry, MpoMatrixType > > > > &W, const vector< Biped< Symmetry, MatrixType > > &Aket, const vector< qarray< Symmetry::Nq > > &qloc, const vector< qarray< Symmetry::Nq > > &qOp, Tripod< Symmetry, MatrixType2 > &Lnew, bool RANDOMIZE=false, tuple< CONTRACT_LR_MODE, size_t > MODE_input=make_pair(FULL, 0), const std::unordered_map< pair< qarray< Symmetry::Nq >, size_t >, size_t > &basis_order_map={})
Definition: DmrgContractions.h:35
contract_LR
Scalar contract_LR(const Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > &L, const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Abra, const vector< vector< vector< Biped< Symmetry, MpoMatrixType > > > > &W, const vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > &Aket, const Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > &R, const vector< qarray< Symmetry::Nq > > &qloc, const vector< qarray< Symmetry::Nq > > &qOp)
Definition: DmrgContractions.h:914
contract::MODE
MODE
Definition: Biped.h:48
qarray2
std::array< qarray< Nq >, 2 > qarray2
Definition: qarray.h:51
qarray3
std::array< qarray< Nq >, 3 > qarray3
Definition: qarray.h:52
Biped
Definition: Biped.h:64
Biped::dict
std::unordered_map< std::array< qType, 2 >, std::size_t > dict
Definition: Biped.h:104
Biped::block
std::vector< MatrixType_ > block
Definition: Biped.h:96
Biped::outerResize
void outerResize(const Biped< Symmetry, OtherMatrixType > Brhs)
Definition: Biped.h:258
Biped::in
std::vector< qType > in
Definition: Biped.h:87
Biped::contract
Biped< Symmetry, MatrixType_ > contract(const Biped< Symmetry, MatrixType_ > &A, const contract::MODE MODE=contract::MODE::UNITY) const
Definition: Biped.h:976
Biped::push_back
void push_back(qType qin, qType qout, const MatrixType_ &M)
Definition: Biped.h:529
Biped::clear
void clear()
Definition: Biped.h:325
Biped::setZero
void setZero()
Definition: Biped.h:336
Biped::dim
std::size_t dim
Definition: Biped.h:82
Biped::out
std::vector< qType > out
Definition: Biped.h:90
Multipede::clear
void clear()
Definition: Multipede.h:409
Multipede::dict
unordered_map< std::array< qType, Nlegs >, size_t > dict
Definition: Multipede.h:95
Multipede::mid
qType mid(size_t q) const
Definition: Multipede.h:154
Multipede::dim
size_t dim
Definition: Multipede.h:66
Multipede::top
qType top(size_t q) const
Definition: Multipede.h:157
Multipede::setZero
void setZero()
Definition: Multipede.h:398
Multipede::out
qType out(size_t q) const
Definition: Multipede.h:153
Multipede::in
qType in(size_t q) const
Definition: Multipede.h:152
Multipede::push_back
void push_back(std::array< qType, Nlegs > quple, const boost::multi_array< MatrixType, LEGLIMIT > &M)
Definition: Multipede.h:439
Multipede::block
vector< boost::multi_array< MatrixType, LEGLIMIT > > block
Definition: Multipede.h:88
Multipede::bot
qType bot(size_t q) const
Definition: Multipede.h:156
qarray
Definition: qarray.h:26