VMPS++
Loading...
Searching...
No Matches
MpsCompressor.h
Go to the documentation of this file.
1#ifndef STRAWBERRY_MPSCOMPRESSOR_WITH_Q
2#define STRAWBERRY_MPSCOMPRESSOR_WITH_Q
3
4#ifndef DMRG_POLYCOMPRESS_TOL
5#define DMRG_POLYCOMPRESS_TOL 1e-4
6#endif
7
8#ifndef DMRG_POLYCOMPRESS_MIN
9#define DMRG_POLYCOMPRESS_MIN 1
10#endif
11
12#ifndef DMRG_POLYCOMPRESS_MAX
13#define DMRG_POLYCOMPRESS_MAX 32
14#endif
15
16#define MPSQCOMPRESSOR_DONT_USE_OPENMP
17
19#include "termcolor.hpp" //from https://github.com/ikalnytskyi/termcolor
21
22#include "Stopwatch.h" // from TOOLS
23#include "LanczosSolver.h" // from ALGS
24
25#include "pivot/DmrgPivotMatrix1.h"
26#include "pivot/DmrgPivotMatrix2.h"
27#include "pivot/DmrgPivotOverlap1.h"
28#include "pivot/DmrgPivotOverlap2.h"
29
36template<typename Symmetry, typename Scalar, typename MpoScalar=double>
37class MpsCompressor
38{
39 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
40
41public:
42
43 MpsCompressor (DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::SILENT)
44 :CHOSEN_VERBOSITY(VERBOSITY)
45 {};
46
47 //---info stuff---
49
50 string info() const;
51
52 string t_info() const;
53
55 double memory (MEMUNIT memunit=GB) const;
57
58 //---compression schemes---
60
72 void stateCompress (const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout,
73 size_t Minit, size_t Mincr=100, size_t Mlimit=10000, double tol=1e-4, size_t max_halfsweeps=40, size_t min_halfsweeps=1);
74
93 template<typename MpOperator>
94 void prodCompress (const MpOperator &H, const MpOperator &Hdag, const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout,
95 qarray<Symmetry::Nq> Qtot_input, size_t Minit, size_t Mincr=100, size_t Mlimit=10000, double tol=1e-4,
96 size_t max_halfsweeps=100, size_t min_halfsweeps=1, std::size_t savePeriod = 0,
97 std::string saveName="backup", bool MEASURE_DISTANCE=true, const MpOperator * HdagH = NULL
98 );
99
119 template<typename MpOperator>
120 void polyCompress (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin1, double polyB, const Mps<Symmetry,Scalar> &Vin2, Mps<Symmetry,Scalar> &Vout,
121 size_t Minit, size_t Mincr=100, size_t Mlimit=10000, double tol=DMRG_POLYCOMPRESS_TOL,
122 size_t max_halfsweeps=DMRG_POLYCOMPRESS_MAX, size_t min_halfsweeps=DMRG_POLYCOMPRESS_MIN);
124
125 void lincomboCompress (const vector<Mps<Symmetry,Scalar> > &Vin, const vector<Scalar> &factors, Mps<Symmetry,Scalar> &Vout, const Mps<Symmetry,Scalar> &Vguess,
126 size_t Minit, size_t Mincr=100, size_t Mlimit=10000, double tol=1e-4, size_t max_halfsweeps=40, size_t min_halfsweeps=1);
127
128private:
129
130 DMRG::VERBOSITY::OPTION CHOSEN_VERBOSITY;
131
132 string print_dist() const;
133
134 // for |Vout> ≈ |Vin>
135 vector<PivotOverlap1<Symmetry,Scalar> > Env;
136
137 void prepSweep (const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout);
138
139 void stateOptimize1 (const Mps<Symmetry,Scalar> &Vin, const Mps<Symmetry,Scalar> &Vout, PivotVector<Symmetry,Scalar> &Aout);
140
141 void stateOptimize1 (const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout);
142
143 void stateOptimize2 (const Mps<Symmetry,Scalar> &Vin, const Mps<Symmetry,Scalar> &Vout,
144 PivotVector<Symmetry,Scalar> &ApairOut);
145
146 void stateOptimize2 (const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout);
147
148 void build_L (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const Mps<Symmetry,Scalar> &Vket, bool RANDOMIZE=false);
149
150 void build_R (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const Mps<Symmetry,Scalar> &Vket, bool RANDOMIZE=false);
151
152 // for |Vout> ≈ H*|Vin>
153 vector<PivotMatrix1<Symmetry,Scalar,MpoScalar> > Heff;
154
155 template<typename MpOperator>
156 void prepSweep (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout, bool RANDOMIZE = false);
157
158 template<typename MpOperator>
159 void prodOptimize1 (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin, const Mps<Symmetry,Scalar> &Vout,
160 PivotVector<Symmetry,Scalar> &Aout);
161
162 template<typename MpOperator>
163 void prodOptimize1 (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout);
164
165 template<typename MpOperator>
166 void prodOptimize2 (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout);
167
168 template<typename MpOperator>
169 void prodOptimize2 (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin, const Mps<Symmetry,Scalar> &Vout,
170 PivotVector<Symmetry,Scalar> &ApairOut);
171
172 template<typename MpOperator>
173 void build_LW (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const MpOperator &H, const Mps<Symmetry,Scalar> &Vket, bool RANDOMIZE=false);
174
175 template<typename MpOperator>
176 void build_RW (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const MpOperator &H, const Mps<Symmetry,Scalar> &Vket, bool RANDOMIZE=false);
177
178 // for |Vout> ≈ H*|Vin1> - polyB*|Vin2>
179 template<typename MpOperator>
180 void prepSweep (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin1, const Mps<Symmetry,Scalar> &Vin2, Mps<Symmetry,Scalar> &Vout,
181 bool RANDOMIZE = false);
182
183 template<typename MpOperator>
184 void build_LRW (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin1, const Mps<Symmetry,Scalar> &Vin2, Mps<Symmetry,Scalar> &Vout);
185
186 // for |Vout> ≈ \sum_i alpha_i |Vin>
187 vector<vector<PivotOverlap1<Symmetry,Scalar> > > Envs;
188
189 void prepSweep (const vector<Mps<Symmetry,Scalar> > &Vin, Mps<Symmetry,Scalar> &Vout);
190
191 void stateOptimize1 (const vector<Mps<Symmetry,Scalar> > &Vin, const Mps<Symmetry,Scalar> &Vout, vector<PivotVector<Symmetry,Scalar> > &Aout);
192
193// void stateOptimize1 (const vector<Mps<Symmetry,Scalar> > &Vin, Mps<Symmetry,Scalar> &Vout);
194
195 void stateOptimize2 (const vector<Mps<Symmetry,Scalar> > &Vin, const Mps<Symmetry,Scalar> &Vout,
196 vector<PivotVector<Symmetry,Scalar> > &ApairOut);
197
198// void stateOptimize2 (const vector<Mps<Symmetry,Scalar> > &Vin, Mps<Symmetry,Scalar> &Vout);
199
200 void build_L (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const vector<Mps<Symmetry,Scalar> > &Vket, bool RANDOMIZE=false);
201
202 void build_R (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const vector<Mps<Symmetry,Scalar> > &Vket, bool RANDOMIZE=false);
203
204 void build_LR (const vector<Mps<Symmetry,Scalar> > &Vin, Mps<Symmetry,Scalar> &Vout);
205
206
207 inline size_t loc1() const {return (CURRENT_DIRECTION==DMRG::DIRECTION::RIGHT)? pivot : pivot-1;};
208 inline size_t loc2() const {return (CURRENT_DIRECTION==DMRG::DIRECTION::RIGHT)? pivot+1 : pivot;};
209
210 void sweep_to_edge (const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout, bool BUILD_LR);
211
212 template<typename MpOperator>
213 void sweep_to_edge (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin1, const Mps<Symmetry,Scalar> &Vin2,
214 Mps<Symmetry,Scalar> &Vout, bool BUILD_LR, bool BUILD_LWRW);
215
216 void sweep_to_edge (const vector<Mps<Symmetry,Scalar> > &Vin, Mps<Symmetry,Scalar> &Vout, bool BUILD_LR);
217
218 size_t N_sites;
219 size_t N_sweepsteps, N_halfsweeps;
220 size_t Mcutoff, Mcutoff_new;
221 size_t Mmax, Mmax_new;
222 double sqdist, tol;
223
224 int pivot;
225 DMRG::DIRECTION::OPTION CURRENT_DIRECTION;
226
227 double t_opt = 0; // optimization
228 double t_AA = 0; // contract_AA
229 double t_sweep = 0; // sweepStep, sweepStep2
230 double t_LR = 0; // build L, R, LW, RW
231 double t_ohead = 0; // precalc_blockStructure
232 double t_tot = 0; // full time step
233};
234
235template<typename Symmetry, typename Scalar, typename MpoScalar>
236string MpsCompressor<Symmetry,Scalar,MpoScalar>::
237info() const
238{
239 stringstream ss;
240 ss << "MpsCompressor: ";
241 ss << "Mcutoff=" << Mcutoff;
242 if (Mcutoff != Mcutoff_new)
243 {
244 ss << "→" << Mcutoff_new << ", ";
245 }
246 else
247 {
248 ss << " (not resized), ";
249 }
250 ss << "Mmax=" << Mmax;
251 if (Mmax != Mmax_new)
252 {
253 ss << "→" << Mmax_new << ", ";
254 }
255 else
256 {
257 ss << " (not changed), ";
258 }
259
260 ss << "|Vlhs-Vrhs|^2=";
261 if (sqdist <= tol) {ss << termcolor::green;}
262 else {ss << termcolor::yellow;}
263 ss << termcolor::reset << sqdist << ", ";
264 ss << "halfsweeps=" << N_halfsweeps << ", ";
265 ss << "mem=" << round(memory(GB),3) << "GB";
266 return ss.str();
267}
268
269template<typename Symmetry, typename Scalar, typename MpoScalar>
270string MpsCompressor<Symmetry,Scalar,MpoScalar>::
271t_info() const
272{
273 stringstream ss;
274 ss << "t[s]=" << termcolor::bold << t_tot << termcolor::reset
275 << ", opt=" << round(t_opt/t_tot*100.,0) << "%"
276 << ", sweep=" << round(t_sweep/t_tot*100.) << "%"
277 << ", LR=" << round(t_LR/t_tot*100.) << "%";
278 if (t_ohead > 0.)
279 {
280 ss << ", ohead=" << round(t_ohead/t_tot*100.) << "%";
281 }
282 if (t_AA > 0.)
283 {
284 ss << ", AA=" << round(t_AA/t_tot*100.) << "%";
285 }
286 return ss.str();
287}
288
289template<typename Symmetry, typename Scalar, typename MpoScalar>
290double MpsCompressor<Symmetry,Scalar,MpoScalar>::
291memory (MEMUNIT memunit) const
292{
293 double res = 0.;
294 for (size_t l=0; l<Env.size(); ++l)
295 {
296 res += Env[l].L.memory(memunit);
297 res += Env[l].R.memory(memunit);
298 }
299 for (size_t l=0; l<Heff.size(); ++l)
300 {
301 res += Heff[l].memory(memunit);
302 }
303 return res;
304}
305
306template<typename Symmetry, typename Scalar, typename MpoScalar>
307template<typename MpOperator>
308void MpsCompressor<Symmetry,Scalar,MpoScalar>::
309sweep_to_edge (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin1, const Mps<Symmetry,Scalar> &Vin2,
310 Mps<Symmetry,Scalar> &Vout, bool BUILD_LR, bool BUILD_LWRW)
311{
312 assert(pivot == 0 or pivot==1 or pivot==N_sites-2 or pivot==N_sites-1);
313
314 if (pivot==1)
315 {
316 Vout.sweep(0,DMRG::BROOM::QR);
317 pivot = 0;
318 if (BUILD_LWRW)
319 {
320 build_RW(0,Vout,H,Vin1);
321 }
322 if (BUILD_LR)
323 {
324 build_R (0,Vout, Vin2);
325 }
326 }
327 else if (pivot==N_sites-2)
328 {
329 Vout.sweep(N_sites-1,DMRG::BROOM::QR);
330 pivot = N_sites-1;
331 if (BUILD_LWRW)
332 {
333 build_LW(N_sites-1,Vout,H,Vin1);
334 }
335 if (BUILD_LR)
336 {
337 build_L (N_sites-1,Vout, Vin2);
338 }
339 }
340}
341
342template<typename Symmetry, typename Scalar, typename MpoScalar>
343void MpsCompressor<Symmetry,Scalar,MpoScalar>::
344sweep_to_edge (const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout, bool BUILD_LR)
345{
346 assert(pivot == 0 or pivot==1 or pivot==N_sites-2 or pivot==N_sites-1);
347
348 if (pivot==1)
349 {
350 Vout.sweep(0,DMRG::BROOM::QR);
351 pivot = 0;
352 if (BUILD_LR)
353 {
354 build_R(0,Vout,Vin);
355 }
356 }
357 else if (pivot==N_sites-2)
358 {
359 Vout.sweep(N_sites-1,DMRG::BROOM::QR);
360 pivot = N_sites-1;
361 if (BUILD_LR)
362 {
363 build_L(N_sites-1,Vout,Vin);
364 }
365 }
366}
367
368template<typename Symmetry, typename Scalar, typename MpoScalar>
369void MpsCompressor<Symmetry,Scalar,MpoScalar>::
370sweep_to_edge (const vector<Mps<Symmetry,Scalar> > &Vin, Mps<Symmetry,Scalar> &Vout, bool BUILD_LR)
371{
372 assert(pivot == 0 or pivot==1 or pivot==N_sites-2 or pivot==N_sites-1);
373
374 if (pivot==1)
375 {
376 Vout.sweep(0,DMRG::BROOM::QR);
377 pivot = 0;
378 if (BUILD_LR)
379 {
380 build_R(0,Vout,Vin);
381 }
382 }
383 else if (pivot==N_sites-2)
384 {
385 Vout.sweep(N_sites-1,DMRG::BROOM::QR);
386 pivot = N_sites-1;
387 if (BUILD_LR)
388 {
389 build_L(N_sites-1,Vout,Vin);
390 }
391 }
392}
393
394//---------------------------compression of |Psi>---------------------------
395// |Vout> ≈ |Vin>, M(Vout) < M(Vin)
396// convention in program: <Vout|Vin>
397
398template<typename Symmetry, typename Scalar, typename MpoScalar>
399void MpsCompressor<Symmetry,Scalar,MpoScalar>::
400stateCompress (const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout,
401 size_t Minit, size_t Mincr, size_t Mlimit, double tol_input, size_t max_halfsweeps, size_t min_halfsweeps)
402{
403 Stopwatch<> Chronos;
404 N_sites = Vin.length();
405 tol = tol_input;
406 double sqnormVin = isReal(dot(Vin,Vin));
407 N_halfsweeps = 0;
408 N_sweepsteps = 0;
409 Mcutoff = Mcutoff_new = Minit;
410
411 // set initial guess
412 Vout = Vin;
413 Vout.innerResize(Minit);
414 Vout.max_Nsv = Minit;
415
416 // set L&R edges
417 Env.clear();
418 Env.resize(N_sites);
419// Env[N_sites-1].R.setTarget(Vin.Qmultitarget());
420// Env[0].L.setVacuum();
421 for (size_t l=0; l<N_sites; ++l)
422 {
423 Env[l].qloc = Vin.locBasis(l);
424 }
425 Env[N_sites-1].R.setIdentity(Vout.outBasis(N_sites-1), Vin.outBasis(N_sites-1));
426 Env[0].L.setIdentity(Vout.inBasis(0), Vin.inBasis(0));
427
428 Mmax = Vout.calc_Mmax();
429 prepSweep(Vin,Vout);
430 sqdist = 1.;
431 size_t halfSweepRange = N_sites;
432
433 if (CHOSEN_VERBOSITY >= DMRG::VERBOSITY::HALFSWEEPWISE)
434 {
435 lout << Chronos.info("preparation stateCompress") << endl;
436 }
437
438 // must achieve sqdist > tol or break off after max_halfsweeps, do at least min_halfsweeps
439 while ((sqdist > tol and N_halfsweeps < max_halfsweeps) or N_halfsweeps < min_halfsweeps or N_halfsweeps%2 != 0)
440 {
441 t_opt = 0;
442 t_AA = 0;
443 t_sweep = 0;
444 t_LR = 0;
445 t_ohead = 0;
446 t_tot = 0;
447 Stopwatch<> FullSweepTimer;
448
449 // A 2-site sweep is necessary! Move pivot back to edge.
450 if (N_halfsweeps%4 == 0 and N_halfsweeps > 1)
451 {
452 sweep_to_edge(Vin,Vout,true); // BUILD_LR = true
453 }
454
455 for (size_t j=1; j<=halfSweepRange; ++j)
456 {
457 turnaround(pivot, N_sites, CURRENT_DIRECTION);
458 if (N_halfsweeps%4 == 0 and N_halfsweeps > 1)
459 {
460 stateOptimize2(Vin,Vout);
461 }
462 else
463 {
464 stateOptimize1(Vin,Vout);
465 }
466 ++N_sweepsteps;
467 }
468 halfSweepRange = N_sites-1;
469 ++N_halfsweeps;
470
471// cout << "sqnormVin=" << sqnormVin << ", Vout.squaredNorm()=" << Vout.squaredNorm() << endl;
472 sqdist = abs(sqnormVin-Vout.squaredNorm());
473 assert(!std::isnan(sqdist));
474
475 if (CHOSEN_VERBOSITY >= DMRG::VERBOSITY::HALFSWEEPWISE)
476 {
477 lout << " distance^2=";
478 if (sqdist <= tol) {lout << termcolor::green;}
479 else {lout << termcolor::yellow;}
480 lout << sqdist << termcolor::reset << ", ";
481 t_tot = FullSweepTimer.time();
482 lout << t_info() << endl;
483 }
484
485 if (N_halfsweeps%4 == 0 and
486 N_halfsweeps > 1 and
487 N_halfsweeps != max_halfsweeps and
488 sqdist > tol)
489 {
490 auto max_Nsv_old = Vout.max_Nsv;
491 Vout.max_Nsv += Mincr;
492 Vout.max_Nsv = min(Vout.max_Nsv,Mlimit);
493 Mcutoff_new = Vout.max_Nsv;
494 if (CHOSEN_VERBOSITY >= DMRG::VERBOSITY::HALFSWEEPWISE)
495 {
496 lout << "resize: " << max_Nsv_old << "→" << Vout.max_Nsv << endl;
497 }
498 }
499
500 Mmax_new = Vout.calc_Mmax();
501 }
502
503 // last sweep
504 sweep_to_edge(Vin,Vout,false);
505}
506
507template<typename Symmetry, typename Scalar, typename MpoScalar>
508void MpsCompressor<Symmetry,Scalar,MpoScalar>::
509prepSweep (const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout)
510{
511 assert(Vout.pivot == 0 or Vout.pivot == N_sites-1 or Vout.pivot == -1);
512
513 if (Vout.pivot == N_sites-1 or Vout.pivot == -1)
514 {
515 for (int l=N_sites-1; l>0; --l)
516 {
517 Vout.sweepStep(DMRG::DIRECTION::LEFT, l, DMRG::BROOM::QR, NULL,true);
518 build_R(l-1,Vout,Vin,true); //true: randomize Env[l].R
519 }
520 CURRENT_DIRECTION = DMRG::DIRECTION::RIGHT;
521 }
522 else if (Vout.pivot == 0)
523 {
524 for (size_t l=0; l<N_sites-1; ++l)
525 {
526 Vout.sweepStep(DMRG::DIRECTION::RIGHT, l, DMRG::BROOM::QR, NULL,true);
527 build_L(l+1,Vout,Vin,true); //true: randomize Env[l].L
528 }
529 CURRENT_DIRECTION = DMRG::DIRECTION::LEFT;
530 }
531 pivot = Vout.pivot;
532}
533
534template<typename Symmetry, typename Scalar, typename MpoScalar>
535void MpsCompressor<Symmetry,Scalar,MpoScalar>::
536stateOptimize1 (const Mps<Symmetry,Scalar> &Vin, const Mps<Symmetry,Scalar> &Vout, PivotVector<Symmetry,Scalar> &Aout)
537{
538 PivotVector<Symmetry,Scalar> Ain(Vin.A[pivot]);
539 Stopwatch<> OptTimer;
540 LRxV(Env[pivot], Ain, Aout);
541 t_opt += OptTimer.time();
542
543 for (size_t s=0; s<Aout.size(); ++s)
544 {
545 Aout[s] = Aout[s].cleaned();
546 }
547}
548
549template<typename Symmetry, typename Scalar, typename MpoScalar>
550void MpsCompressor<Symmetry,Scalar,MpoScalar>::
551stateOptimize1 (const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout)
552{
553 PivotVector<Symmetry,Scalar> Aout;
554 stateOptimize1(Vin,Vout,Aout);
555 Vout.A[pivot] = Aout.data;
556
557// // safeguard against sudden norm loss:
558// if (Vout.squaredNorm() < 1e-7)
559// {
560// cout << "Vout.squaredNorm()=" << Vout.squaredNorm() << endl;
561// if (CHOSEN_VERBOSITY > 0)
562// {
563// lout << termcolor::bold << termcolor::red << "WARNING: small norm encountered at pivot=" << pivot << "!" << termcolor::reset << endl;
564// }
565// Vout /= sqrt(Vout.squaredNorm());
566// }
567
568 Stopwatch<> SweepTimer;
569 Vout.sweepStep(CURRENT_DIRECTION, pivot, DMRG::BROOM::SVD);
570 t_sweep += SweepTimer.time();
571 pivot = Vout.get_pivot();
572 (CURRENT_DIRECTION == DMRG::DIRECTION::RIGHT)? build_L(pivot,Vout,Vin) : build_R(pivot,Vout,Vin);
573}
574
575template<typename Symmetry, typename Scalar, typename MpoScalar>
576void MpsCompressor<Symmetry,Scalar,MpoScalar>::
577stateOptimize2 (const Mps<Symmetry,Scalar> &Vin, const Mps<Symmetry,Scalar> &Vout, PivotVector<Symmetry,Scalar> &ApairOut)
578{
579 Stopwatch<> AAtimer;
580 PivotVector<Symmetry,Scalar> ApairIn(Vin.A[loc1()], Vin.locBasis(loc1()),
581 Vin.A[loc2()], Vin.locBasis(loc2()),
582 Vin.QoutTop[loc1()], Vin.QoutBot[loc1()]);
583
584 ApairOut = PivotVector<Symmetry,Scalar>(Vout.A[loc1()], Vout.locBasis(loc1()),
585 Vout.A[loc2()], Vout.locBasis(loc2()),
586 Vout.QoutTop[loc1()], Vout.QoutBot[loc1()],
587 true); // dry run: do not multiply matrices, just set blocks
588 t_AA += AAtimer.time();
589
590 Stopwatch<> OptTimer;
591 PivotOverlap2<Symmetry,Scalar> Env2(Env[loc1()].L, Env[loc2()].R, Vin.locBasis(loc1()), Vin.locBasis(loc2()));
592 LRxV(Env2, ApairIn, ApairOut);
593 t_opt += OptTimer.time();
594
595 for (size_t s=0; s<ApairOut.data.size(); ++s)
596 {
597 ApairOut.data[s] = ApairOut.data[s].cleaned();
598 }
599}
600
601template<typename Symmetry, typename Scalar, typename MpoScalar>
602void MpsCompressor<Symmetry,Scalar,MpoScalar>::
603stateOptimize2 (const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout)
604{
605 PivotVector<Symmetry,Scalar> Apair;
606 stateOptimize2(Vin,Vout,Apair);
607
608 Stopwatch<> SweepTimer;
609 Vout.sweepStep2(CURRENT_DIRECTION, loc1(), Apair.data);
610 t_sweep += SweepTimer.time();
611 pivot = Vout.get_pivot();
612 (CURRENT_DIRECTION == DMRG::DIRECTION::RIGHT)? build_L(pivot,Vout,Vin) : build_R(pivot,Vout,Vin);
613}
614
615template<typename Symmetry, typename Scalar, typename MpoScalar>
616void MpsCompressor<Symmetry,Scalar,MpoScalar>::
617build_L (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const Mps<Symmetry,Scalar> &Vket, bool RANDOMIZE)
618{
619 Stopwatch<> LRtimer;
620 contract_L(Env[loc-1].L, Vbra.A[loc-1], Vket.A[loc-1], Vket.locBasis(loc-1), Env[loc].L, RANDOMIZE, true); // CLEAR=true
621 t_LR += LRtimer.time();
622}
623
624template<typename Symmetry, typename Scalar, typename MpoScalar>
625void MpsCompressor<Symmetry,Scalar,MpoScalar>::
626build_R (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const Mps<Symmetry,Scalar> &Vket, bool RANDOMIZE)
627{
628 Stopwatch<> LRtimer;
629 contract_R(Env[loc+1].R, Vbra.A[loc+1], Vket.A[loc+1], Vket.locBasis(loc+1), Env[loc].R, RANDOMIZE, true); // CLEAR=true
630 t_LR += LRtimer.time();
631}
632
633//---------------------------compression of \sum_n \alpha_n |Psi_n>---------------------------
634// |Vout> ≈ Σₙ αₙ |Ψₙ>
635// convention in program: <Vout|H|Vin>
636
637template<typename Symmetry, typename Scalar, typename MpoScalar>
638void MpsCompressor<Symmetry,Scalar,MpoScalar>::
639lincomboCompress (const vector<Mps<Symmetry,Scalar> > &Vin, const vector<Scalar> &factors, Mps<Symmetry,Scalar> &Vout, const Mps<Symmetry,Scalar> &Vguess,
640 size_t Minit, size_t Mincr, size_t Mlimit, double tol_input, size_t max_halfsweeps, size_t min_halfsweeps)
641{
642 Stopwatch<> Chronos;
643 N_sites = Vin[0].length();
644 tol = tol_input;
645 size_t N_vecs = Vin.size();
646 N_halfsweeps = 0;
647 N_sweepsteps = 0;
648 Mcutoff = Mcutoff_new = Minit;
649
650 // Calculate Hermitian overlap matrix Σₙₘ αₘ* αₙ <Ψₘ|Ψₙ>
651 Matrix<Scalar,Dynamic,Dynamic> overlapsVin(N_vecs,N_vecs);
652 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
653 #pragma omp parallel for
654 #endif
655 for (int i=0; i<N_vecs; ++i)
656 for (int j=0; j<=i; ++j)
657 {
658 overlapsVin(i,j) = conjIfcomplex(factors[i]) * factors[j] * dot(Vin[i],Vin[j]);
659// #pragma omp critical
660// {
661// cout << "i=" << i << ", j=" << j << ", overlapsVin=" << overlapsVin(i,j) << endl;
662// }
663 }
664 overlapsVin.template triangularView<Upper>() = overlapsVin.adjoint();
665 double overlapsVinSum = isReal(overlapsVin.sum());
666
667 // set L&R edges
668 Envs.resize(N_vecs);
669 for (int i=0; i<N_vecs; ++i)
670 {
671 Envs[i].clear();
672 Envs[i].resize(N_sites);
673 Envs[i][N_sites-1].R.setTarget(Vin[i].Qmultitarget());
674 Envs[i][0].L.setVacuum();
675 for (size_t l=0; l<N_sites; ++l)
676 {
677 Envs[i][l].qloc = Vin[i].locBasis(l);
678 }
679 }
680 // set initial guess
681// Vout = Vin[0];
682 Vout = Vguess;
683// Vout.innerResize(Mcutoff);
684// Vout.max_Nsv = Mcutoff;
685
686 Mmax = Vout.calc_Mmax();
687 prepSweep(Vin,Vout);
688 sqdist = 1.;
689 size_t halfSweepRange = N_sites;
690
691 if (CHOSEN_VERBOSITY>=2)
692 {
693 lout << Chronos.info("preparation lincomboCompress") << endl;
694 }
695
696 // must achieve sqdist > tol or break off after max_halfsweeps, do at least min_halfsweeps
697 while ((sqdist > tol and N_halfsweeps < max_halfsweeps) or N_halfsweeps < min_halfsweeps or N_halfsweeps%2 != 0)
698 {
699 t_opt = 0;
700 t_AA = 0;
701 t_sweep = 0;
702 t_LR = 0;
703 t_ohead = 0;
704 t_tot = 0;
705 Stopwatch<> FullSweepTimer;
706
707 // A 2-site sweep is necessary! Move pivot back to edge.
708 if (N_halfsweeps%4 == 0 and N_halfsweeps > 1)
709 {
710// cout << "sweep_to_edge" << endl;
711 sweep_to_edge(Vin,Vout,true); // BUILD_LR = true
712 }
713
714 for (size_t j=1; j<=halfSweepRange; ++j)
715 {
716// cout << "j=" << j << endl;
717 turnaround(pivot, N_sites, CURRENT_DIRECTION);
718 Stopwatch<> Chronos;
719
720 if (N_halfsweeps%4 == 0 and N_halfsweeps > 1)
721 {
722 vector<PivotVector<Symmetry,Scalar> > Apair(N_vecs);
723// cout << "begin stateOptimize2" << endl;
724 stateOptimize2(Vin,Vout,Apair);
725// cout << "end stateOptimize2" << endl;
726
727 PivotVector<Symmetry,Scalar> Asum;
728 Asum.outerResize(Apair[0]);
729// cout << "outer resize done!" << endl;
730
731 for (int i=0; i<N_vecs; ++i)
732 for (size_t s=0; s<Asum.size(); ++s)
733 {
734 Asum[s].addScale(factors[i], Apair[i][s]);
735 }
736// cout << "addScale done!" << endl;
737
738 for (size_t s=0; s<Asum.size(); ++s)
739 {
740 Asum[s] = Asum[s].cleaned();
741 }
742// cout << "cleaned done!" << endl;
743
744 Stopwatch<> SweepTimer;
745// cout << "begin sweepStep2" << endl;
746 Vout.sweepStep2(CURRENT_DIRECTION, loc1(), Asum.data);
747// cout << "end sweepStep2" << endl;
748 t_sweep += SweepTimer.time();
749 pivot = Vout.get_pivot();
750 }
751 else
752 {
753 vector<PivotVector<Symmetry,Scalar> > Ares(N_vecs);
754// cout << "begin stateOptimize1" << endl;
755 stateOptimize1(Vin,Vout,Ares);
756// cout << "end stateOptimize1" << endl;
757
758 Stopwatch<> OheadTimer;
759 double Vsqnorm = Vout.squaredNorm();
760 t_ohead += OheadTimer.time();
761 if (Vsqnorm < 1e-7)
762 {
763 /*if (CHOSEN_VERBOSITY > 0)
764 {
765 lout << termcolor::bold << termcolor::yellow << "WARNING: small norm encountered at pivot=" << pivot << "!" << termcolor::reset << endl;
766 }*/
767 Vout /= sqrt(Vsqnorm);
768 }
769
770 PivotVector<Symmetry,Scalar> Asum;
771 Asum.outerResize(Ares[0]);
772
773 for (int i=0; i<N_vecs; ++i)
774 for (size_t s=0; s<Asum.size(); ++s)
775 {
776 Asum[s].addScale(factors[i], Ares[i][s]);
777 }
778
779 for (size_t s=0; s<Asum.size(); ++s)
780 {
781 Asum[s] = Asum[s].cleaned();
782 }
783
784 Vout.A[pivot] = Asum.data;
785
786 Stopwatch<> SweepTimer;
787 Vout.sweepStep(CURRENT_DIRECTION, pivot, DMRG::BROOM::SVD);
788 t_sweep += SweepTimer.time();
789 pivot = Vout.get_pivot();
790 }
791
792 build_LR(Vin,Vout);
793 ++N_sweepsteps;
794 }
795 halfSweepRange = N_sites-1;
796 ++N_halfsweeps;
797
798 //cout << "sqnormVin=" << overlapsVin.sum() << ", Vout.squaredNorm()=" << Vout.squaredNorm() << endl;
799 sqdist = abs(overlapsVinSum-Vout.squaredNorm());
800// cout << "Vout.squaredNorm()=" << Vout.squaredNorm() << endl;
801 assert(!std::isnan(sqdist));
802
803 if (CHOSEN_VERBOSITY>=2)
804 {
805 lout << " distance^2=";
806 if (sqdist <= tol) {lout << termcolor::green;}
807 else {lout << termcolor::yellow;}
808 lout << sqdist << termcolor::reset << ", ";
809 t_tot = FullSweepTimer.time();
810 lout << t_info() << endl;
811 }
812
813 if (N_halfsweeps%4 == 0 and
814 N_halfsweeps > 1 and
815 N_halfsweeps != max_halfsweeps and
816 sqdist > tol)
817 {
818 auto max_Nsv_old = Vout.max_Nsv;
819 Vout.max_Nsv += Mincr;
820 Vout.max_Nsv = min(Vout.max_Nsv,Mlimit);
821 Mcutoff_new = Vout.max_Nsv;
822 if (CHOSEN_VERBOSITY >= DMRG::VERBOSITY::HALFSWEEPWISE)
823 {
824 lout << "resize: " << max_Nsv_old << "→" << Vout.max_Nsv << endl;
825 }
826 }
827
828 Mmax_new = Vout.calc_Mmax();
829 }
830
831 // last sweep
832 sweep_to_edge(Vin,Vout,false);
833}
834
835template<typename Symmetry, typename Scalar, typename MpoScalar>
836void MpsCompressor<Symmetry,Scalar,MpoScalar>::
837prepSweep (const vector<Mps<Symmetry,Scalar> > &Vin, Mps<Symmetry,Scalar> &Vout)
838{
839 assert(Vout.pivot == 0 or Vout.pivot == N_sites-1 or Vout.pivot == -1);
840// Vout.setRandom();
841
842 if (Vout.pivot == N_sites-1 or
843 Vout.pivot == -1)
844 {
845 for (int l=N_sites-1; l>0; --l)
846 {
847 Vout.sweepStep(DMRG::DIRECTION::LEFT, l, DMRG::BROOM::QR, NULL,true);
848 build_R(l-1,Vout,Vin,true); // true: randomize Env[l].R
849 }
850 CURRENT_DIRECTION = DMRG::DIRECTION::RIGHT;
851 }
852 else if (Vout.pivot == 0)
853 {
854 for (size_t l=0; l<N_sites-1; ++l)
855 {
856 Vout.sweepStep(DMRG::DIRECTION::RIGHT, l, DMRG::BROOM::QR, NULL,true);
857 build_L(l+1,Vout,Vin,true); // true: randomize Env[l].L
858 }
859 CURRENT_DIRECTION = DMRG::DIRECTION::LEFT;
860 }
861 pivot = Vout.pivot;
862}
863
864template<typename Symmetry, typename Scalar, typename MpoScalar>
865void MpsCompressor<Symmetry,Scalar,MpoScalar>::
866stateOptimize1 (const vector<Mps<Symmetry,Scalar> > &Vin, const Mps<Symmetry,Scalar> &Vout, vector<PivotVector<Symmetry,Scalar> > &Aout)
867{
868 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
869 #pragma omp parallel for
870 #endif
871 for (int i=0; i<Vin.size(); ++i)
872 {
873 PivotVector<Symmetry,Scalar> Ain(Vin[i].A[pivot]);
874 Stopwatch<> OptTimer;
875 LRxV(Envs[i][pivot], Ain, Aout[i]);
876 t_opt += OptTimer.time();
877
878 for (size_t s=0; s<Aout[i].size(); ++s)
879 {
880 Aout[i][s] = Aout[i][s].cleaned();
881 }
882 }
883}
884
885template<typename Symmetry, typename Scalar, typename MpoScalar>
886void MpsCompressor<Symmetry,Scalar,MpoScalar>::
887stateOptimize2 (const vector<Mps<Symmetry,Scalar> > &Vin, const Mps<Symmetry,Scalar> &Vout, vector<PivotVector<Symmetry,Scalar> > &ApairOut)
888{
889 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
890 #pragma omp parallel for
891 #endif
892 for (int i=0; i<Vin.size(); ++i)
893 {
894 Stopwatch<> AAtimer;
895 PivotVector<Symmetry,Scalar> ApairIn(Vin[i].A[loc1()], Vin[i].locBasis(loc1()),
896 Vin[i].A[loc2()], Vin[i].locBasis(loc2()),
897 Vin[i].QoutTop[loc1()], Vin[i].QoutBot[loc1()]);
898
899 ApairOut[i] = PivotVector<Symmetry,Scalar>(Vout.A[loc1()], Vout.locBasis(loc1()),
900 Vout.A[loc2()], Vout.locBasis(loc2()),
901 Vout.QoutTop[loc1()], Vout.QoutBot[loc1()],
902 true); // dry run: do not multiply matrices, just set blocks
903 t_AA += AAtimer.time();
904
905 Stopwatch<> OptTimer;
906 PivotOverlap2<Symmetry,Scalar> Env2(Envs[i][loc1()].L, Envs[i][loc2()].R, Vin[i].locBasis(loc1()), Vin[i].locBasis(loc2()));
907 LRxV(Env2, ApairIn, ApairOut[i]);
908 t_opt += OptTimer.time();
909
910 for (size_t s=0; s<ApairOut[i].data.size(); ++s)
911 {
912 ApairOut[i].data[s] = ApairOut[i].data[s].cleaned();
913 }
914 }
915}
916
917template<typename Symmetry, typename Scalar, typename MpoScalar>
918void MpsCompressor<Symmetry,Scalar,MpoScalar>::
919build_L (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const vector<Mps<Symmetry,Scalar> > &Vket, bool RANDOMIZE)
920{
921 Stopwatch<> LRtimer;
922 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
923 #pragma omp parallel for
924 #endif
925 for (int i=0; i<Vket.size(); ++i)
926 {
927 contract_L(Envs[i][loc-1].L, Vbra.A[loc-1], Vket[i].A[loc-1], Vket[i].locBasis(loc-1), Envs[i][loc].L, RANDOMIZE, true); // CLEAR=true
928 }
929 t_LR += LRtimer.time();
930}
931
932template<typename Symmetry, typename Scalar, typename MpoScalar>
933void MpsCompressor<Symmetry,Scalar,MpoScalar>::
934build_R (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const vector<Mps<Symmetry,Scalar> > &Vket, bool RANDOMIZE)
935{
936 Stopwatch<> LRtimer;
937 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
938 #pragma omp parallel for
939 #endif
940 for (int i=0; i<Vket.size(); ++i)
941 {
942 contract_R(Envs[i][loc+1].R, Vbra.A[loc+1], Vket[i].A[loc+1], Vket[i].locBasis(loc+1), Envs[i][loc].R, RANDOMIZE, true); // CLEAR=true
943 }
944 t_LR += LRtimer.time();
945}
946
947template<typename Symmetry, typename Scalar, typename MpoScalar>
948void MpsCompressor<Symmetry,Scalar,MpoScalar>::
949build_LR (const vector<Mps<Symmetry,Scalar> > &Vin, Mps<Symmetry,Scalar> &Vout)
950{
951 (CURRENT_DIRECTION == DMRG::DIRECTION::RIGHT)? build_L (pivot,Vout,Vin) : build_R(pivot,Vout,Vin);
952}
953
954//---------------------------compression of H*|Psi>---------------------------
955// |Vout> ≈ H|Vin>
956// convention in program: <Vout|H|Vin>
957
958template<typename Symmetry, typename Scalar, typename MpoScalar>
959template<typename MpOperator>
960void MpsCompressor<Symmetry,Scalar,MpoScalar>::
961prodCompress (const MpOperator &H, const MpOperator &Hdag, const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout, qarray<Symmetry::Nq> Qtot_input,
962 size_t Minit, size_t Mincr, size_t Mlimit, double tol_input, size_t max_halfsweeps, size_t min_halfsweeps,
963 std::size_t savePeriod, std::string saveName, bool MEASURE_DISTANCE, const MpOperator * HdagH)
964{
965 if (CHOSEN_VERBOSITY>=2)
966 {
967 lout << endl << termcolor::colorize << termcolor::bold
968 << "———————————————————————————————————————————prodCompress: |Φ> = H|Ψ>————————————————————————————————————————————"
969 << termcolor::reset << endl;
970 }
971
972 N_sites = Vin.length();
973 Stopwatch<> Chronos;
974 tol = tol_input;
975 N_halfsweeps = 0;
976 N_sweepsteps = 0;
977 Mcutoff = Mcutoff_new = Minit;
978
979 if (H.Qtarget() == Symmetry::qvacuum())
980 {
981 Vout = Vin;
982 }
983 else
984 {
985 Vout = Mps<Symmetry,Scalar>(H, Mcutoff, Qtot_input, max(static_cast<int>(Vin.calc_Nqmax()), DMRG::CONTROL::DEFAULT::Qinit));
986 }
987
988 // prepare edges of LW & RW
989 Heff.clear();
990 Heff.resize(N_sites);
991 for (int l=0; l<N_sites; ++l) Heff[l].Terms.resize(1);
992 Heff[0].Terms[0].L.setVacuum();
993
994 vector<qarray3<Symmetry::Nq> > Qt;
995 for (size_t i=0; i<Vin.Qmultitarget().size(); ++i)
996 {
997 Qt.push_back(qarray3<Symmetry::Nq>{Vin.Qmultitarget()[i], Vout.Qmultitarget()[i], H.Qtarget()});
998 }
999 Heff[N_sites-1].Terms[0].R.setTarget(Qt);
1000
1001 for (size_t l=0; l<N_sites; ++l)
1002 {
1003 Heff[l].Terms[0].W = H.W[l];
1004 }
1005
1006 Vout.max_Nsv = Mcutoff;
1007 Vout.min_Nsv = Vin.min_Nsv;
1008 Mmax = Vout.calc_Mmax();
1009 double avgHsqVin;
1010
1011 if (MEASURE_DISTANCE)
1012 {
1013 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1014 #pragma omp parallel sections
1015 #endif
1016 {
1017 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1018 #pragma omp section
1019 #endif
1020 {
1021 if (H.IS_UNITARY())
1022 {
1023 avgHsqVin = Vin.squaredNorm();
1024 }
1025 else
1026 {
1027 if (HdagH != NULL)
1028 {
1029 avgHsqVin = isReal(avg(Vin,*HdagH,Vin));
1030 }
1031 else
1032 {
1033 if (H.IS_HERMITIAN())
1034 {
1035 avgHsqVin = (H.maxPower()>=2)? isReal(avg(Vin,H,Vin,2)) : isReal(avg(Vin,H,H,Vin));
1036 }
1037 else
1038 {
1039 avgHsqVin = isReal(avg(Vin,Hdag,H,Vin));
1040 }
1041 }
1042 }
1043 }
1044 }
1045 }
1046 prepSweep(H,Vin,Vout);
1047 sqdist = 1.;
1048 size_t halfSweepRange = N_sites;
1049
1050 if (CHOSEN_VERBOSITY>=2)
1051 {
1052 lout << Chronos.info("• preparation prodCompress") << endl;
1053 size_t standard_precision = cout.precision();
1054 lout << "• initial state : " << Vout.info() << endl;
1055 lout << "• bond dim. increase by ";
1056 cout << termcolor::underline;
1057 lout << Mincr;
1058 cout << termcolor::reset;
1059 lout << " every ";
1060 cout << termcolor::underline;
1061 lout << "4";
1062 cout << termcolor::reset;
1063 lout << " half-sweeps" << endl;
1064 lout << "• make between ";
1065 cout << termcolor::underline;
1066 lout << min_halfsweeps;
1067 cout << termcolor::reset;
1068 lout << " and ";
1069 cout << termcolor::underline;
1070 lout << max_halfsweeps;
1071 cout << termcolor::reset;
1072 lout << " half-sweep iterations" << endl;
1073 lout << "• convergence tolerance: ";
1074 cout << termcolor::underline;
1075 lout << tol_input;
1076 cout << termcolor::reset;
1077 lout << endl << endl;
1078 }
1079
1080 // must achieve sqdist > tol or break off after max_halfsweeps, do at least min_halfsweeps
1081 while ((sqdist > tol and N_halfsweeps < max_halfsweeps) or N_halfsweeps < min_halfsweeps)
1082 {
1083 t_opt = 0;
1084 t_AA = 0;
1085 t_sweep = 0;
1086 t_LR = 0;
1087 t_ohead = 0;
1088 t_tot = 0;
1089 Stopwatch<> FullSweepTimer;
1090
1091 // A 2-site sweep is necessary! Move pivot back to edge.
1092 if (N_halfsweeps%4 == 0 and N_halfsweeps > 1)
1093 {
1094 sweep_to_edge(H,Vin,Vin,Vout,false,true); // build_LWRW = true
1095 }
1096
1097 double norm_old = Vout.squaredNorm();
1098
1099 // optimization
1100 for (size_t j=1; j<=halfSweepRange; ++j)
1101 {
1102 turnaround(pivot, N_sites, CURRENT_DIRECTION);
1103 if (N_halfsweeps%4 == 0 and N_halfsweeps > 1)
1104 {
1105 prodOptimize2(H,Vin,Vout);
1106 }
1107 else
1108 {
1109 prodOptimize1(H,Vin,Vout);
1110 }
1111 ++N_sweepsteps;
1112 }
1113 halfSweepRange = N_sites-1;
1114 ++N_halfsweeps;
1115
1116 if (MEASURE_DISTANCE)
1117 {
1118 cout << "avgHsqVin=" << avgHsqVin << endl;
1119 cout << "Vout.squaredNorm()=" << Vout.squaredNorm() << endl;
1120 sqdist = abs(avgHsqVin - Vout.squaredNorm());
1121 assert(!std::isnan(sqdist));
1122 }
1123 else
1124 {
1125 double norm_new = Vout.squaredNorm();
1126 sqdist = abs(norm_new-norm_old);
1127 }
1128
1129 if (CHOSEN_VERBOSITY>=2)
1130 {
1131 cout << termcolor::underline;
1132 lout << "half-sweeps=" << N_halfsweeps;
1133 cout << termcolor::reset;
1134 size_t standard_precision = cout.precision();
1135 lout << ", distance^2=";
1136 if (sqdist <= tol) {cout << termcolor::green;}
1137 else {cout << termcolor::yellow;}
1138 lout << sqdist;
1139 cout << termcolor::reset;
1140 lout << endl;
1141 t_tot = FullSweepTimer.time();
1142 lout << t_info() << endl;
1143 lout << Vout.info() << endl;
1144 }
1145
1146 bool RESIZED = false;
1147 if (N_halfsweeps%4 == 0 and
1148 N_halfsweeps > 1 and
1149 N_halfsweeps != max_halfsweeps and
1150 sqdist > tol)
1151 {
1152 //size_t Delta_Nsv = (sqdist>10.*tol)? 40:20;
1153// Vout.max_Nsv += Delta_Nsv;
1154 auto max_Nsv_old = Vout.max_Nsv;
1155 Vout.max_Nsv += Mincr;
1156 Vout.max_Nsv = min(Vout.max_Nsv,Mlimit);
1157 Mcutoff_new = Vout.max_Nsv;
1158 if (CHOSEN_VERBOSITY >= DMRG::VERBOSITY::HALFSWEEPWISE)
1159 {
1160 lout << "resize: " << max_Nsv_old << "→" << Vout.max_Nsv << endl;
1161 }
1162 }
1163
1164 Mmax_new = Vout.calc_Mmax();
1165
1166 #ifdef USE_HDF5_STORAGE
1167 if (savePeriod != 0 and N_halfsweeps%savePeriod == 0)
1168 {
1169 Vout.save(saveName,H.info());
1170 cout << termcolor::green;
1171 lout << "Saved state to: " << saveName;
1172 cout << termcolor::reset;
1173 lout << endl;
1174 }
1175 #endif
1176
1177 #ifdef COMPRESSOR_RESTART_FROM_RANDOM
1178 if (N_halfsweeps == max_halfsweeps/2 and sqdist > tol)
1179 {
1180 cout << termcolor::red;
1181 lout << "Warning: Could not reach tolerance, restarting from random!";
1182 cout << termcolor::reset;
1183 lout << endl;
1184 prepSweep(H,Vin,Vout,true);
1185 }
1186 #endif
1187 if (CHOSEN_VERBOSITY>=2)
1188 {
1189 lout << endl;
1190 }
1191 }
1192
1193 // move pivot to edge at the end
1194// if (pivot==1) {Vout.sweep(0,DMRG::BROOM::QR);}
1195// else if (pivot==N_sites-2) {Vout.sweep(N_sites-1,DMRG::BROOM::QR);}
1196 Scalar factor_cgc = 1.;//Symmetry::degeneracy(H.Qtarget());
1197 Vout *= factor_cgc;
1198 sweep_to_edge(H,Vin,Vin,Vout,false,false);
1199}
1200
1201template<typename Symmetry, typename Scalar, typename MpoScalar>
1202template<typename MpOperator>
1203void MpsCompressor<Symmetry,Scalar,MpoScalar>::
1204prepSweep (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout, bool RANDOMIZE)
1205{
1206 assert(Vout.pivot == 0 or Vout.pivot == N_sites-1 or Vout.pivot == -1);
1207// Vout.setRandom();
1208
1209// bool RANDOMIZE;
1210// if (H.IS_HERMITIAN())
1211// {
1212// Vout = Vin;
1213// RANDOMIZE = false;
1214// }
1215// else
1216// {
1217// RANDOMIZE = true;
1218// }
1219
1220 if (Vout.pivot == N_sites-1 or Vout.pivot == -1)
1221 {
1222 for (int l=N_sites-1; l>0; --l)
1223 {
1224 Vout.sweepStep(DMRG::DIRECTION::LEFT, l, DMRG::BROOM::QR, NULL,true);
1225 build_RW(l-1,Vout,H,Vin,RANDOMIZE);
1226 }
1227 CURRENT_DIRECTION = DMRG::DIRECTION::RIGHT;
1228 }
1229 else if (Vout.pivot == 0)
1230 {
1231 for (size_t l=0; l<N_sites-1; ++l)
1232 {
1233 Vout.sweepStep(DMRG::DIRECTION::RIGHT, l, DMRG::BROOM::QR, NULL,true);
1234 build_LW(l+1,Vout,H,Vin,RANDOMIZE);
1235 }
1236 CURRENT_DIRECTION = DMRG::DIRECTION::LEFT;
1237 }
1238
1239 pivot = Vout.pivot;
1240}
1241
1242template<typename Symmetry, typename Scalar, typename MpoScalar>
1243template<typename MpOperator>
1244void MpsCompressor<Symmetry,Scalar,MpoScalar>::
1245prodOptimize1 (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin, const Mps<Symmetry,Scalar> &Vout, PivotVector<Symmetry,Scalar> &Aout)
1246{
1247 Stopwatch<> OheadTimer;
1248 precalc_blockStructure (Heff[pivot].Terms[0].L, Vout.A[pivot], Heff[pivot].Terms[0].W, Vin.A[pivot], Heff[pivot].Terms[0].R,
1249 H.locBasis(pivot), H.opBasis(pivot),
1250 Heff[pivot].qlhs, Heff[pivot].qrhs, Heff[pivot].factor_cgcs);
1251
1252 PivotVector<Symmetry,Scalar> Ain(Vin.A[pivot]);
1253 Aout = PivotVector<Symmetry,Scalar>(Vout.A[pivot]);
1254 t_ohead += OheadTimer.time();
1255
1256 Stopwatch<> OptTimer;
1257 OxV(Heff[pivot], Ain, Aout);
1258 t_opt += OptTimer.time();
1259}
1260
1261template<typename Symmetry, typename Scalar, typename MpoScalar>
1262template<typename MpOperator>
1263void MpsCompressor<Symmetry,Scalar,MpoScalar>::
1264prodOptimize1 (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout)
1265{
1266 PivotVector<Symmetry,Scalar> Aout;
1267 prodOptimize1(H,Vin,Vout,Aout);
1268 Vout.A[pivot] = Aout.data;
1269
1270 // safeguard against sudden norm loss:
1271 Stopwatch<> OheadTimer;
1272 double Vsqnorm = Vout.squaredNorm();
1273 t_ohead += OheadTimer.time();
1274 if (Vsqnorm < 1e-7)
1275 {
1276 /*if (CHOSEN_VERBOSITY > 0)
1277 {
1278 lout << termcolor::bold << termcolor::yellow << "WARNING: small norm encountered at pivot=" << pivot << "!" << termcolor::reset << endl;
1279 }*/
1280 Vout /= sqrt(Vsqnorm);
1281 }
1282
1283 Stopwatch<> SweepTimer;
1284 Vout.sweepStep(CURRENT_DIRECTION, pivot, DMRG::BROOM::SVD);
1285 t_sweep += SweepTimer.time();
1286
1287 pivot = Vout.get_pivot();
1288 (CURRENT_DIRECTION==DMRG::DIRECTION::RIGHT)? build_LW(pivot,Vout,H,Vin) : build_RW(pivot,Vout,H,Vin);
1289}
1290
1291template<typename Symmetry, typename Scalar, typename MpoScalar>
1292template<typename MpOperator>
1293void MpsCompressor<Symmetry,Scalar,MpoScalar>::
1294prodOptimize2 (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin, const Mps<Symmetry,Scalar> &Vout, PivotVector<Symmetry,Scalar> &ApairOut)
1295{
1296 Stopwatch<> AAtimer;
1297 PivotVector<Symmetry,Scalar> ApairIn(Vin.A[loc1()], Vin.locBasis(loc1()),
1298 Vin.A[loc2()], Vin.locBasis(loc2()),
1299 Vin.QoutTop[loc1()], Vin.QoutBot[loc1()]);
1300
1301 ApairOut = PivotVector<Symmetry,Scalar>(Vout.A[loc1()], Vout.locBasis(loc1()),
1302 Vout.A[loc2()], Vout.locBasis(loc2()),
1303 Vout.QoutTop[loc1()], Vout.QoutBot[loc1()],
1304 true); // dry run: do not multiply matrices, just set blocks
1305 t_AA += AAtimer.time();
1306
1307 PivotMatrix2<Symmetry,Scalar,MpoScalar> Heff2(Heff[loc1()].Terms[0].L, Heff[loc2()].Terms[0].R,
1308 H.W[loc1()], H.W[loc2()],
1309 H.locBasis(loc1()), H.locBasis(loc2()),
1310 H.opBasis (loc1()), H.opBasis (loc2()));
1311
1312 Stopwatch<> OheadTimer;
1313 precalc_blockStructure (Heff[loc1()].Terms[0].L, ApairOut.data, Heff2.Terms[0].W12, Heff2.Terms[0].W34, ApairIn.data, Heff[loc2()].Terms[0].R,
1314 H.locBasis(loc1()), H.locBasis(loc2()), H.opBasis(loc1()), H.opBasis(loc2()),
1315 Heff2.qlhs, Heff2.qrhs, Heff2.factor_cgcs);
1316 t_ohead += OheadTimer.time();
1317
1318 Stopwatch<> OptTimer;
1319 OxV(Heff2, ApairIn, ApairOut);
1320 t_opt += OptTimer.time();
1321
1322 for (size_t s=0; s<ApairOut.data.size(); ++s)
1323 {
1324 ApairOut.data[s] = ApairOut.data[s].cleaned();
1325 }
1326}
1327
1328template<typename Symmetry, typename Scalar, typename MpoScalar>
1329template<typename MpOperator>
1330void MpsCompressor<Symmetry,Scalar,MpoScalar>::
1331prodOptimize2 (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin, Mps<Symmetry,Scalar> &Vout)
1332{
1333 Stopwatch<> Chronos;
1334
1335 Stopwatch<> OptTimer;
1336 PivotVector<Symmetry,Scalar> Apair;
1337 prodOptimize2(H,Vin,Vout,Apair);
1338 t_opt += OptTimer.time();
1339
1340 Stopwatch<> SweepTimer;
1341 Vout.sweepStep2(CURRENT_DIRECTION, loc1(), Apair.data);
1342 t_sweep += SweepTimer.time();
1343
1344 pivot = Vout.get_pivot();
1345
1346 (CURRENT_DIRECTION == DMRG::DIRECTION::RIGHT)? build_LW(pivot,Vout,H,Vin) : build_RW(pivot,Vout,H,Vin);
1347}
1348
1349template<typename Symmetry, typename Scalar, typename MpoScalar>
1350template<typename MpOperator>
1351void MpsCompressor<Symmetry,Scalar,MpoScalar>::
1352build_LW (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const MpOperator &H, const Mps<Symmetry,Scalar> &Vket, bool RANDOMIZE)
1353{
1354 Stopwatch<> LRtimer;
1355 contract_L(Heff[loc-1].Terms[0].L, Vbra.A[loc-1], H.W[loc-1], Vket.A[loc-1], H.locBasis(loc-1), H.opBasis(loc-1), Heff[loc].Terms[0].L, RANDOMIZE);
1356 t_LR += LRtimer.time();
1357}
1358
1359template<typename Symmetry, typename Scalar, typename MpoScalar>
1360template<typename MpOperator>
1361void MpsCompressor<Symmetry,Scalar,MpoScalar>::
1362build_RW (size_t loc, const Mps<Symmetry,Scalar> &Vbra, const MpOperator &H, const Mps<Symmetry,Scalar> &Vket, bool RANDOMIZE)
1363{
1364 Stopwatch<> LRtimer;
1365 contract_R(Heff[loc+1].Terms[0].R, Vbra.A[loc+1], H.W[loc+1], Vket.A[loc+1], H.locBasis(loc+1), H.opBasis(loc+1), Heff[loc].Terms[0].R, RANDOMIZE);
1366 t_LR += LRtimer.time();
1367}
1368
1369template<typename Symmetry, typename Scalar, typename MpoScalar>
1370template<typename MpOperator>
1371void MpsCompressor<Symmetry,Scalar,MpoScalar>::
1372polyCompress (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin1, double polyB, const Mps<Symmetry,Scalar> &Vin2, Mps<Symmetry,Scalar> &Vout,
1373 size_t Minit, size_t Mincr, size_t Mlimit, double tol_input, size_t max_halfsweeps, size_t min_halfsweeps)
1374{
1375 N_sites = Vin1.length();
1376 tol = tol_input;
1377 Stopwatch<> Chronos;
1378 N_halfsweeps = 0;
1379 N_sweepsteps = 0;
1380 Mcutoff = Mcutoff_new = Minit;
1381
1382 Vout = Vin1;
1383// Vout = Mps<Symmetry,Scalar>(H, Mcutoff, Vin1.Qtarget(), max(Vin1.calc_Nqmax(), DMRG::CONTROL::DEFAULT::Qinit));
1384// Vout.setRandom();
1385 if (CHOSEN_VERBOSITY>=2)
1386 {
1387 lout << "Vin1: " << Vin1.info() << endl;
1388 lout << "Vin2: " << Vin2.info() << endl;
1389 }
1390
1391// // prepare edges of LW & RW
1392// Heff.clear();
1393// Heff.resize(N_sites);
1394// Heff[0].L.setVacuum();
1395//
1397// vector<qarray3<Symmetry::Nq> > Qt;
1398// for (size_t i=0; i<Vin1.Qmultitarget().size(); ++i)
1399// {
1400// Qt.push_back(qarray3<Symmetry::Nq>{Vin1.Qmultitarget()[i], Vout.Qmultitarget()[i], H.Qtarget()});
1401// }
1402// Heff[N_sites-1].R.setTarget(Qt);
1403
1404 Heff.clear();
1405 Heff.resize(N_sites);
1406 Heff[0].L = Vin1.get_boundaryTensor(DMRG::DIRECTION::LEFT);
1407 Heff[N_sites-1].R = Vin1.get_boundaryTensor(DMRG::DIRECTION::RIGHT);
1408
1409 for (size_t l=0; l<N_sites; ++l)
1410 {
1411 Heff[l].W = H.W[l];
1412 }
1413
1414// // set L&R edges
1415// Env.clear();
1416// Env.resize(N_sites);
1417// Env[N_sites-1].R.setTarget(Vin2.Qmultitarget());
1418// Env[0].L.setVacuum();
1419// for (size_t l=0; l<N_sites; ++l)
1420// {
1421// Env[l].qloc = H.locBasis(l);
1422// }
1423
1424 Env.clear();
1425 Env.resize(N_sites);
1426 for (size_t l=0; l<N_sites; ++l)
1427 {
1428 Env[l].qloc = Vin2.locBasis(l);
1429 }
1430 Env[N_sites-1].R.setIdentity(Vin2.outBasis(N_sites-1), Vout.outBasis(N_sites-1));
1431 Env[0].L.setIdentity(Vin2.inBasis(0), Vout.inBasis(0));
1432
1433 double avgHsqV1, sqnormV2, overlapV12;
1434 sqnormV2 = Vin2.squaredNorm();
1435
1436 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1437 #pragma omp parallel sections
1438 #endif
1439 {
1440 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1441 #pragma omp section
1442 #endif
1443 {
1444 if (Vin1.Boundaries.IS_TRIVIAL())
1445 {
1446 avgHsqV1 = (H.maxPower()>=2)? isReal(avg(Vin1,H,Vin1,2)) : isReal(avg(Vin1,H,H,Vin1));
1447 }
1448 else
1449 {
1450 avgHsqV1 = avg_hetero(Vin1,H,Vin1,true,true); // USE_BOUNDARY=true, USE_SQUARE=true
1451 }
1452 }
1453 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1454 #pragma omp section
1455 #endif
1456 {
1457 if (Vin1.Boundaries.IS_TRIVIAL())
1458 {
1459 overlapV12 = isReal(avg(Vin2,H,Vin1));
1460 }
1461 else
1462 {
1463 overlapV12 = isReal(avg_hetero(Vin2,H,Vin1,true)); // USE_BOUNDARY=true
1464 }
1465 }
1466 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1467 #pragma omp section
1468 #endif
1469 {
1470 prepSweep(H,Vin1,Vin2,Vout);
1471 }
1472 }
1473 sqdist = 1.;
1474 size_t halfSweepRange = N_sites;
1475
1476 if (CHOSEN_VERBOSITY>=2)
1477 {
1478 lout << Chronos.info("preparation polyCompress") << endl;
1479 }
1480
1481 Mmax = Vout.calc_Mmax();
1482 Vout.max_Nsv = Mcutoff;
1483 Vout.min_Nsv = Vin1.min_Nsv;
1484 // In order to avoid block loss for small Hilbert spaces:
1485 if (Vout.calc_Nqmax() <= 4)
1486 {
1487 Vout.min_Nsv = 1;
1488 }
1489
1490 // must achieve sqdist > tol or break off after max_halfsweeps, do at least min_halfsweeps
1491 while ((sqdist > tol and N_halfsweeps < max_halfsweeps) or N_halfsweeps < min_halfsweeps)
1492 {
1493 t_opt = 0;
1494 t_AA = 0;
1495 t_sweep = 0;
1496 t_LR = 0;
1497 t_ohead = 0;
1498 t_tot = 0;
1499 Stopwatch<> FullSweepTimer;
1500
1501 // A 2-site sweep is necessary! Move pivot back to edge.
1502 if (N_halfsweeps%4 == 0 and N_halfsweeps > 1)
1503 {
1504 sweep_to_edge(H,Vin1,Vin2,Vout,true,true); // build_LR = true, build LWRW = true
1505 }
1506
1507 for (size_t j=1; j<=halfSweepRange; ++j)
1508 {
1509 turnaround(pivot, N_sites, CURRENT_DIRECTION);
1510 Stopwatch<> Chronos;
1511
1512 if (N_halfsweeps%4 == 0 and N_halfsweeps > 1)
1513 {
1514 PivotVector<Symmetry,Scalar> Apair1;
1515 prodOptimize2(H,Vin1,Vout,Apair1);
1516
1517 PivotVector<Symmetry,Scalar> Apair2;
1518 stateOptimize2(Vin2,Vout,Apair2);
1519
1520 for (size_t s=0; s<Apair1.size(); ++s)
1521 {
1522 Apair1[s].addScale(-polyB, Apair2[s]);
1523 Apair1[s] = Apair1[s].cleaned();
1524 }
1525
1526 Stopwatch<> SweepTimer;
1527 Vout.sweepStep2(CURRENT_DIRECTION, loc1(), Apair1.data);
1528 t_sweep += SweepTimer.time();
1529 pivot = Vout.get_pivot();
1530 }
1531 else
1532 {
1533 PivotVector<Symmetry,Scalar> A1;
1534 prodOptimize1(H,Vin1,Vout,A1);
1535
1536 PivotVector<Symmetry,Scalar> A2;
1537 stateOptimize1(Vin2,Vout,A2);
1538
1539 for (size_t s=0; s<A1.size(); ++s)
1540 {
1541 A1[s].addScale(-polyB, A2[s]);
1542 A1[s] = A1[s].cleaned();
1543 }
1544 Vout.A[pivot] = A1.data;
1545
1546 Stopwatch<> SweepTimer;
1547 Vout.sweepStep(CURRENT_DIRECTION, pivot, DMRG::BROOM::SVD);
1548 t_sweep += SweepTimer.time();
1549 pivot = Vout.get_pivot();
1550 }
1551
1552 build_LRW(H,Vin1,Vin2,Vout);
1553 ++N_sweepsteps;
1554 }
1555 halfSweepRange = N_sites-1;
1556 ++N_halfsweeps;
1557
1558// cout << "avgHsqV1=" << avgHsqV1
1559// << ", Vout.squaredNorm()=" << Vout.squaredNorm()
1560// << ", polyB*polyB*sqnormV2=" << polyB*polyB*sqnormV2
1561// << ", 2.*polyB*overlapV12=" << 2.*polyB*overlapV12
1562// << endl;
1563 double sqdist_ = sqdist;
1564 sqdist = abs(avgHsqV1 - Vout.squaredNorm() + polyB*polyB*sqnormV2 - 2.*polyB*overlapV12);
1565// lout << "diff=" << abs(sqdist_-sqdist) << endl;
1566 assert(!std::isnan(sqdist));
1567
1568 if (CHOSEN_VERBOSITY>=2)
1569 {
1570 lout << " distance^2=";
1571 if (sqdist <= tol) {lout << termcolor::green;}
1572 else {lout << termcolor::yellow;}
1573 lout << sqdist << termcolor::reset << ", ";
1574 t_tot = FullSweepTimer.time();
1575 lout << t_info() << endl;
1576 }
1577
1578 bool RESIZED = false;
1579 if (N_halfsweeps%4 == 0 and
1580 N_halfsweeps > 1 and
1581 N_halfsweeps != max_halfsweeps and
1582 sqdist > tol)
1583 {
1584 auto max_Nsv_old = Vout.max_Nsv;
1585 Vout.max_Nsv += Mincr;
1586 Vout.max_Nsv = min(Vout.max_Nsv,Mlimit);
1587 Mcutoff_new = Vout.max_Nsv;
1588 if (CHOSEN_VERBOSITY >= DMRG::VERBOSITY::HALFSWEEPWISE)
1589 {
1590 lout << "resize: " << max_Nsv_old << "→" << Vout.max_Nsv << endl;
1591 }
1592 }
1593
1594 if (N_halfsweeps == max_halfsweeps/2 and sqdist > tol)
1595 {
1596 lout << termcolor::red << "Warning: Could not reach tolerance, restarting from random!" << termcolor::reset << endl;
1597 prepSweep(H,Vin1,Vin2,Vout,true);
1598 }
1599
1600 Mmax_new = Vout.calc_Mmax();
1601 }
1602
1603 // last sweep
1604 sweep_to_edge(H,Vin1,Vin2,Vout,false,false);
1605}
1606
1607template<typename Symmetry, typename Scalar, typename MpoScalar>
1608template<typename MpOperator>
1609void MpsCompressor<Symmetry,Scalar,MpoScalar>::
1610prepSweep (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin1, const Mps<Symmetry,Scalar> &Vin2, Mps<Symmetry,Scalar> &Vout, bool RANDOMIZE)
1611{
1612 assert(Vout.pivot == 0 or Vout.pivot == N_sites-1 or Vout.pivot == -1);
1613// if (RANDOMIZE) {Vout.setRandom();}
1614
1615 if (Vout.pivot == N_sites-1)
1616 {
1617 for (int l=N_sites-1; l>0; --l)
1618 {
1619 Vout.sweepStep(DMRG::DIRECTION::LEFT, l, DMRG::BROOM::QR, NULL,RANDOMIZE);
1620 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1621 #pragma omp parallel sections
1622 #endif
1623 {
1624 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1625 #pragma omp section
1626 #endif
1627 {
1628 build_RW(l-1,Vout,H,Vin1,RANDOMIZE);
1629 }
1630 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1631 #pragma omp section
1632 #endif
1633 {
1634 build_R(l-1,Vout,Vin2,RANDOMIZE);
1635 }
1636 }
1637 }
1638 CURRENT_DIRECTION = DMRG::DIRECTION::RIGHT;
1639 }
1640 else if (Vout.pivot == 0 or Vout.pivot == -1)
1641 {
1642 for (size_t l=0; l<N_sites-1; ++l)
1643 {
1644 Vout.sweepStep(DMRG::DIRECTION::RIGHT, l, DMRG::BROOM::QR, NULL,RANDOMIZE);
1645 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1646 #pragma omp parallel sections
1647 #endif
1648 {
1649 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1650 #pragma omp section
1651 #endif
1652 {
1653 build_LW(l+1,Vout,H,Vin1,RANDOMIZE);
1654 }
1655 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1656 #pragma omp section
1657 #endif
1658 {
1659 build_L(l+1,Vout,Vin2,RANDOMIZE);
1660 }
1661 }
1662 }
1663 CURRENT_DIRECTION = DMRG::DIRECTION::LEFT;
1664 }
1665 pivot = Vout.pivot;
1666}
1667
1668template<typename Symmetry, typename Scalar, typename MpoScalar>
1669template<typename MpOperator>
1670void MpsCompressor<Symmetry,Scalar,MpoScalar>::
1671build_LRW (const MpOperator &H, const Mps<Symmetry,Scalar> &Vin1, const Mps<Symmetry,Scalar> &Vin2, Mps<Symmetry,Scalar> &Vout)
1672{
1673 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1674 #pragma omp parallel sections
1675 #endif
1676 {
1677 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1678 #pragma omp section
1679 #endif
1680 {
1681 (CURRENT_DIRECTION == DMRG::DIRECTION::RIGHT)? build_LW(pivot,Vout,H,Vin1) : build_RW(pivot,Vout,H,Vin1);
1682 }
1683 #ifndef MPSQCOMPRESSOR_DONT_USE_OPENMP
1684 #pragma omp section
1685 #endif
1686 {
1687 (CURRENT_DIRECTION == DMRG::DIRECTION::RIGHT)? build_L (pivot,Vout,Vin2) : build_R(pivot,Vout,Vin2);
1688 }
1689 }
1690}
1691
1692#endif
precalc_blockStructure
void precalc_blockStructure(const Tripod< Symmetry, Eigen::Matrix< Scalar, Dynamic, Dynamic > > &L, const vector< Biped< Symmetry, Eigen::Matrix< Scalar, Dynamic, Dynamic > > > &Abra, const vector< vector< vector< Biped< Symmetry, MpoMatrixType > > > > &W, const vector< Biped< Symmetry, Eigen::Matrix< Scalar, Dynamic, Dynamic > > > &Aket, const Tripod< Symmetry, Eigen::Matrix< Scalar, Dynamic, Dynamic > > &R, const vector< qarray< Symmetry::Nq > > &qloc, const vector< qarray< Symmetry::Nq > > &qOp, vector< std::array< size_t, 2 > > &qlhs, vector< vector< std::array< size_t, 6 > > > &qrhs, vector< vector< Scalar > > &factor_cgcs)
Definition: DmrgIndexGymnastics.h:488
turnaround
void turnaround(int pivot, size_t L, DMRG::DIRECTION::OPTION &DIR)
Flips the sweep direction when the edge is reached.
Definition: DmrgJanitor.h:7
dot
Scalar dot(const Mps< Symmetry, Scalar > &Vbra, const Mps< Symmetry, Scalar > &Vket)
Definition: DmrgLinearAlgebra.h:24
OxV
void OxV(const Mpo< Symmetry, MpoScalar > &H, const Mpo< Symmetry, MpoScalar > &Hdag, const Mps< Symmetry, Scalar > &Vin, Mps< Symmetry, Scalar > &Vout, bool VERBOSE=true, double tol_compr=1e-4, int Mincr=100, int Mlimit=10000, int max_halfsweeps=100, int min_halfsweeps=1, bool MEASURE_DISTANCE=true)
Definition: DmrgLinearAlgebra.h:593
avg_hetero
Scalar avg_hetero(const Mps< Symmetry, Scalar > &Vbra, const Mpo< Symmetry, MpoScalar > &O, const Mps< Symmetry, Scalar > &Vket, bool USE_BOUNDARY=false, size_t usePower=1ul, const qarray< Symmetry::Nq > &Qmid=Symmetry::qvacuum())
Definition: DmrgLinearAlgebra.h:287
avg
Scalar avg(const Mps< Symmetry, Scalar > &Vbra, const Mpo< Symmetry, MpoScalar > &O, const Mps< Symmetry, Scalar > &Vket, size_t power_of_O=1ul, DMRG::DIRECTION::OPTION DIR=DMRG::DIRECTION::RIGHT, DMRG::VERBOSITY::OPTION CHOSEN_VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: DmrgLinearAlgebra.h:149
LRxV
void LRxV(const PivotOverlap1< Symmetry, Scalar > &H, const PivotVector< Symmetry, Scalar > &Vin, PivotVector< Symmetry, Scalar > &Vout)
Definition: DmrgPivotOverlap1.h:26
conjIfcomplex
double conjIfcomplex(double x)
Definition: DmrgTypedefs.h:27
isReal
double isReal(double x)
Definition: DmrgTypedefs.h:21
A
@ A
Definition: DmrgTypedefs.h:130
DMRG_POLYCOMPRESS_TOL
#define DMRG_POLYCOMPRESS_TOL
Definition: MpsCompressor.h:5
DMRG_POLYCOMPRESS_MIN
#define DMRG_POLYCOMPRESS_MIN
Definition: MpsCompressor.h:9
DMRG_POLYCOMPRESS_MAX
#define DMRG_POLYCOMPRESS_MAX
Definition: MpsCompressor.h:13
DmrgJanitor::sweepStep
void sweepStep(DMRG::DIRECTION::OPTION DIR, size_t loc, DMRG::BROOM::OPTION TOOL, PivotMatrixType *H=NULL, bool DISCARD_U_or_V=false)
Definition: DmrgJanitor.h:218
DmrgJanitor::sweep
void sweep(size_t loc, DMRG::BROOM::OPTION TOOL, PivotMatrixType *H=NULL)
Definition: DmrgJanitor.h:198
DmrgJanitor::max_Nsv
size_t max_Nsv
Definition: DmrgJanitor.h:98
DmrgJanitor::length
size_t length() const
Definition: DmrgJanitor.h:92
DmrgJanitor::min_Nsv
size_t min_Nsv
Definition: DmrgJanitor.h:98
MpsCompressor::loc2
size_t loc2() const
Definition: MpsCompressor.h:208
MpsCompressor::stateOptimize2
void stateOptimize2(const Mps< Symmetry, Scalar > &Vin, const Mps< Symmetry, Scalar > &Vout, PivotVector< Symmetry, Scalar > &ApairOut)
Definition: MpsCompressor.h:577
MpsCompressor::prodOptimize2
void prodOptimize2(const MpOperator &H, const Mps< Symmetry, Scalar > &Vin, Mps< Symmetry, Scalar > &Vout)
Definition: MpsCompressor.h:1331
MpsCompressor::MatrixType
Matrix< Scalar, Dynamic, Dynamic > MatrixType
Definition: MpsCompressor.h:39
MpsCompressor::info
string info() const
Definition: MpsCompressor.h:237
MpsCompressor::stateOptimize1
void stateOptimize1(const Mps< Symmetry, Scalar > &Vin, const Mps< Symmetry, Scalar > &Vout, PivotVector< Symmetry, Scalar > &Aout)
Definition: MpsCompressor.h:536
MpsCompressor::prodCompress
void prodCompress(const MpOperator &H, const MpOperator &Hdag, const Mps< Symmetry, Scalar > &Vin, Mps< Symmetry, Scalar > &Vout, qarray< Symmetry::Nq > Qtot_input, size_t Minit, size_t Mincr=100, size_t Mlimit=10000, double tol=1e-4, size_t max_halfsweeps=100, size_t min_halfsweeps=1, std::size_t savePeriod=0, std::string saveName="backup", bool MEASURE_DISTANCE=true, const MpOperator *HdagH=NULL)
Definition: MpsCompressor.h:961
MpsCompressor::Mcutoff_new
size_t Mcutoff_new
Definition: MpsCompressor.h:220
MpsCompressor::CHOSEN_VERBOSITY
DMRG::VERBOSITY::OPTION CHOSEN_VERBOSITY
Definition: MpsCompressor.h:130
MpsCompressor::Envs
vector< vector< PivotOverlap1< Symmetry, Scalar > > > Envs
Definition: MpsCompressor.h:187
MpsCompressor::lincomboCompress
void lincomboCompress(const vector< Mps< Symmetry, Scalar > > &Vin, const vector< Scalar > &factors, Mps< Symmetry, Scalar > &Vout, const Mps< Symmetry, Scalar > &Vguess, size_t Minit, size_t Mincr=100, size_t Mlimit=10000, double tol=1e-4, size_t max_halfsweeps=40, size_t min_halfsweeps=1)
Definition: MpsCompressor.h:639
MpsCompressor::polyCompress
void polyCompress(const MpOperator &H, const Mps< Symmetry, Scalar > &Vin1, double polyB, const Mps< Symmetry, Scalar > &Vin2, Mps< Symmetry, Scalar > &Vout, size_t Minit, size_t Mincr=100, size_t Mlimit=10000, double tol=DMRG_POLYCOMPRESS_TOL, size_t max_halfsweeps=DMRG_POLYCOMPRESS_MAX, size_t min_halfsweeps=DMRG_POLYCOMPRESS_MIN)
Definition: MpsCompressor.h:1372
MpsCompressor::build_R
void build_R(size_t loc, const Mps< Symmetry, Scalar > &Vbra, const Mps< Symmetry, Scalar > &Vket, bool RANDOMIZE=false)
Definition: MpsCompressor.h:626
MpsCompressor::memory
double memory(MEMUNIT memunit=GB) const
Definition: MpsCompressor.h:291
MpsCompressor::Heff
vector< PivotMatrix1< Symmetry, Scalar, MpoScalar > > Heff
Definition: MpsCompressor.h:153
MpsCompressor::stateCompress
void stateCompress(const Mps< Symmetry, Scalar > &Vin, Mps< Symmetry, Scalar > &Vout, size_t Minit, size_t Mincr=100, size_t Mlimit=10000, double tol=1e-4, size_t max_halfsweeps=40, size_t min_halfsweeps=1)
Definition: MpsCompressor.h:400
MpsCompressor::build_L
void build_L(size_t loc, const Mps< Symmetry, Scalar > &Vbra, const Mps< Symmetry, Scalar > &Vket, bool RANDOMIZE=false)
Definition: MpsCompressor.h:617
MpsCompressor::print_dist
string print_dist() const
MpsCompressor::N_halfsweeps
size_t N_halfsweeps
Definition: MpsCompressor.h:219
MpsCompressor::build_LRW
void build_LRW(const MpOperator &H, const Mps< Symmetry, Scalar > &Vin1, const Mps< Symmetry, Scalar > &Vin2, Mps< Symmetry, Scalar > &Vout)
Definition: MpsCompressor.h:1671
MpsCompressor::N_sweepsteps
size_t N_sweepsteps
Definition: MpsCompressor.h:219
MpsCompressor::sweep_to_edge
void sweep_to_edge(const Mps< Symmetry, Scalar > &Vin, Mps< Symmetry, Scalar > &Vout, bool BUILD_LR)
Definition: MpsCompressor.h:344
MpsCompressor::CURRENT_DIRECTION
DMRG::DIRECTION::OPTION CURRENT_DIRECTION
Definition: MpsCompressor.h:225
MpsCompressor::loc1
size_t loc1() const
Definition: MpsCompressor.h:207
MpsCompressor::build_RW
void build_RW(size_t loc, const Mps< Symmetry, Scalar > &Vbra, const MpOperator &H, const Mps< Symmetry, Scalar > &Vket, bool RANDOMIZE=false)
Definition: MpsCompressor.h:1362
MpsCompressor::MpsCompressor
MpsCompressor(DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: MpsCompressor.h:43
MpsCompressor::build_LW
void build_LW(size_t loc, const Mps< Symmetry, Scalar > &Vbra, const MpOperator &H, const Mps< Symmetry, Scalar > &Vket, bool RANDOMIZE=false)
Definition: MpsCompressor.h:1352
MpsCompressor::prodOptimize1
void prodOptimize1(const MpOperator &H, const Mps< Symmetry, Scalar > &Vin, const Mps< Symmetry, Scalar > &Vout, PivotVector< Symmetry, Scalar > &Aout)
Definition: MpsCompressor.h:1245
MpsCompressor::Env
vector< PivotOverlap1< Symmetry, Scalar > > Env
Definition: MpsCompressor.h:135
MpsCompressor::prepSweep
void prepSweep(const Mps< Symmetry, Scalar > &Vin, Mps< Symmetry, Scalar > &Vout)
Definition: MpsCompressor.h:509
MpsCompressor::t_info
string t_info() const
Definition: MpsCompressor.h:271
MpsCompressor::build_LR
void build_LR(const vector< Mps< Symmetry, Scalar > > &Vin, Mps< Symmetry, Scalar > &Vout)
Definition: MpsCompressor.h:949
Mps
Definition: Mps.h:39
Mps::locBasis
vector< qarray< Nq > > locBasis(size_t loc) const
Definition: Mps.h:444
Mps::get_boundaryTensor
Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > get_boundaryTensor(DMRG::DIRECTION::OPTION DIR, size_t usePower=1ul) const
Definition: Mps.h:488
Mps::calc_Nqmax
size_t calc_Nqmax() const
Mps::Boundaries
MpsBoundaries< Symmetry, Scalar > Boundaries
Definition: Mps.h:481
Mps::QoutBot
vector< qarray< Nq > > QoutBot
Definition: Mps.h:642
Mps::innerResize
void innerResize(size_t Mmax)
Mps::QoutTop
vector< qarray< Nq > > QoutTop
Definition: Mps.h:641
Mps::info
string info() const
Mps::squaredNorm
double squaredNorm() const
Mps::inBasis
Qbasis< Symmetry > inBasis(size_t loc) const
Definition: Mps.h:448
Mps::calc_Mmax
size_t calc_Mmax() const
Mps::get_pivot
int get_pivot() const
Definition: Mps.h:462
Mps::Qmultitarget
vector< qarray< Nq > > Qmultitarget() const
Definition: Mps.h:441
Mps::sweepStep2
void sweepStep2(DMRG::DIRECTION::OPTION DIR, size_t loc, const vector< Biped< Symmetry, MatrixType > > &Apair, bool SEPARATE_SV=false)
Mps::outBasis
Qbasis< Symmetry > outBasis(size_t loc) const
Definition: Mps.h:452
Mps::A
vector< vector< Biped< Symmetry, MatrixType > > > A
Definition: Mps.h:624
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
qarray3
std::array< qarray< Nq >, 3 > qarray3
Definition: qarray.h:52
DMRG::CONTROL::DEFAULT::Qinit
static constexpr int Qinit
Definition: DmrgTypedefs.h:342
PivotMatrix2::qrhs
vector< vector< std::array< size_t, 12 > > > qrhs
Definition: DmrgPivotMatrix2.h:68
PivotMatrix2::qlhs
vector< std::array< size_t, 2 > > qlhs
Definition: DmrgPivotMatrix2.h:67
PivotMatrix2::factor_cgcs
vector< vector< Scalar > > factor_cgcs
Definition: DmrgPivotMatrix2.h:69
PivotMatrix2::Terms
vector< PivotMatrix2Terms< Symmetry, Scalar, MpoScalar > > Terms
Definition: DmrgPivotMatrix2.h:64
PivotVector::size
size_t size() const
Definition: DmrgPivotVector.h:74
PivotVector::outerResize
void outerResize(const PivotVector &Vrhs)
Definition: DmrgPivotVector.h:51
PivotVector::data
vector< Biped< Symmetry, Matrix< Scalar_, Dynamic, Dynamic > > > data
Definition: DmrgPivotVector.h:116
qarray
Definition: qarray.h:26