VMPS++
Loading...
Searching...
No Matches
Umps.h
Go to the documentation of this file.
1#ifndef VANILLA_Umps
2#define VANILLA_Umps
3
5#include <set>
6#include <numeric>
7#include <algorithm>
8#include <ctime>
9#include <type_traits>
10#include <iostream>
11#include <fstream>
13
14#include "RandomVector.h" // from ALGS
15
16#include "VUMPS/VumpsTypedefs.h"
17#include "VumpsContractions.h"
18#include "Blocker.h"
19#include "VumpsTransferMatrixSF.h"
20#include "VumpsTransferMatrixQ.h"
21#include "Mps.h"
22
23#ifdef USE_HDF5_STORAGE
24 #include <HDF5Interface.h>
25#endif
26//include "PolychromaticConsole.h" // from TOOLS
27//include "tensors/Biped.h"
28//include "LanczosSolver.h" // from ALGS
29//include "ArnoldiSolver.h" // from ALGS
30//include "VUMPS/VumpsTransferMatrix.h"
31//include "Mpo.h"
32//include "tensors/DmrgConglutinations.h"
33
40template<typename Symmetry, typename Scalar=double>
41class Umps
42{
43 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
44 static constexpr size_t Nq = Symmetry::Nq;
45
46 template<typename Symmetry_, typename MpHamiltonian, typename Scalar_> friend class VumpsSolver;
47 template<typename Symmetry_, typename S1, typename S2> friend class MpsCompressor;
48
49public:
50
52 Umps(){};
53
55 template<typename Hamiltonian> Umps (const Hamiltonian &H, qarray<Nq> Qtot_input, size_t L_input, size_t Mmax, size_t Nqmax, bool INIT_TO_HALF_INTEGER_SPIN);
56
58 Umps (const vector<qarray<Symmetry::Nq> > &qloc_input, qarray<Nq> Qtot_input, size_t L_input, size_t Mmax, size_t Nqmax, bool INIT_TO_HALF_INTEGER_SPIN);
59
61 Umps (const vector<vector<qarray<Symmetry::Nq> > > &qloc_input, qarray<Nq> Qtot_input, size_t L_input, size_t Mmax, size_t Nqmax, bool INIT_TO_HALF_INTEGER_SPIN);
62
63 #ifdef USE_HDF5_STORAGE
69 Umps (string filename) {load(filename);}
70 #endif //USE_HDF5_STORAGE
71
73 string info() const;
74
76 void graph (string filename) const;
77
79 string test_ortho (double tol=1e-6) const;
80
82 void setRandom();
83
85 void normalize_C();
86
88 void resize (size_t Mmax_input, size_t Nqmax_input, bool INIT_TO_HALF_INTEGER_SPIN);
89
95 template<typename OtherMatrixType> void resize (const Umps<Symmetry,OtherMatrixType> &V);
96
98 void resize_arrays();
99
102 void svdDecompose (size_t loc, GAUGE::OPTION gauge = GAUGE::C);
103
106 void polarDecompose (size_t loc, GAUGE::OPTION gauge = GAUGE::C);
107
109 VectorXd entropy() const {return S;};
110
112 inline vector<map<qarray<Nq>,tuple<ArrayXd,int> > > entanglementSpectrum() const {return SVspec;};
113
115 std::pair<vector<qarray<Symmetry::Nq> >, ArrayXd> entanglementSpectrumLoc (size_t loc) const;
116
120 template<typename OtherScalar> Umps<Symmetry,OtherScalar> cast() const;
121
126 Umps<Symmetry,double> real() const;
127
129 inline vector<qarray<Symmetry::Nq> > locBasis (size_t loc) const {return qloc[loc];}
130 inline vector<vector<qarray<Symmetry::Nq> > > locBasis() const {return qloc;}
131
133 inline Qbasis<Symmetry> inBasis (size_t loc) const {return inbase[loc];}
134 inline vector<Qbasis<Symmetry> > inBasis() const {return inbase;}
135
137 inline Qbasis<Symmetry> outBasis (size_t loc) const {return outbase[loc];}
138 inline vector<Qbasis<Symmetry> > outBasis() const {return outbase;}
139
141 size_t get_frst_rows() const {return A[GAUGE::C][0][0].block[0].rows();}
142
144 size_t get_last_cols() const {return A[GAUGE::C][N_sites-1][0].block[0].cols();}
145
147 size_t length() const {return N_sites;}
148
150 Scalar calc_epsLRsq (GAUGE::OPTION gauge, size_t loc) const;
151
152 #ifdef USE_HDF5_STORAGE
154
161 void save (string filename, string info="none", double energy=std::nan("1"), double err_var=std::nan("1"), double err_state=std::nan("1"));
162
168 void load (string filename, double &energy=dump_Mps, double &err_var=dump_Mps, double &err_state=dump_Mps);
169 #endif //USE_HDF5_STORAGE
170
174 size_t calc_Mmax() const;
175
179 size_t calc_fullMmax() const;
180
184 size_t calc_Dmax() const;
185
189 size_t calc_Nqmax() const;
190
192 double memory (MEMUNIT memunit) const;
193
198 double dot (const Umps<Symmetry,Scalar> &Vket) const;
199
201 const vector<Biped<Symmetry,MatrixType> > &A_at (GAUGE::OPTION g, size_t loc) const {return A[g][loc];};
202
204 qarray<Symmetry::Nq> Qtop (size_t loc) const;
205 qarray<Symmetry::Nq> Qbot (size_t loc) const;
206
208 inline size_t minus1modL (size_t l) const {return (l==0)? N_sites-1 : (l-1);}
209
211 inline qarray<Nq> Qtarget() const {return Qtot;};
212
213 void calc_N (DMRG::DIRECTION::OPTION DIR, size_t loc, vector<Biped<Symmetry,MatrixType> > &N) const;
214
222 void truncate(bool SET_AC_RANDOM=true);
223
230 void orthogonalize_left (GAUGE::OPTION g, vector<Biped<Symmetry,MatrixType> > &G_L);
237 void orthogonalize_right(GAUGE::OPTION g, vector<Biped<Symmetry,MatrixType> > &G_R);
238
245 std::pair<complex<double>, Biped<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > >
246 calc_dominant_1symm (GAUGE::OPTION g, DMRG::DIRECTION::OPTION DIR, const Mpo<Symmetry,complex<double>> &R, bool TRANSPOSE, bool CONJUGATE) const;
247
248 std::pair<complex<double>, Biped<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > >
249 calc_dominant_2symm (GAUGE::OPTION g, DMRG::DIRECTION::OPTION DIR, const Mpo<Symmetry,complex<double>> &R1, const Mpo<Symmetry,complex<double>> &R2) const;
250
251 vector<std::pair<complex<double>,Biped<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > > >
252 calc_dominant (GAUGE::OPTION g=GAUGE::R, DMRG::DIRECTION::OPTION DIR=DMRG::DIRECTION::RIGHT, int N=2, double tol=1e-15, int dimK=-1, qarray<Symmetry::Nq> Qtot=Symmetry::qvacuum(), string label="") const;
253
254 template<typename MpoScalar>
255 vector<std::pair<complex<double>,Tripod<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > > >
256 calc_dominant_Q (const Mpo<Symmetry,MpoScalar> &O, GAUGE::OPTION g=GAUGE::R, DMRG::DIRECTION::OPTION DIR=DMRG::DIRECTION::RIGHT, int N=2, double tol=1e-15, int dimK=-1, string label="") const;
257
265 void adjustQN (const size_t number_cells);
266
270 void sort_A (size_t loc, GAUGE::OPTION g, bool SORT_ALL_GAUGES=false);
271
276 void updateC (size_t loc);
283 void updateAC (size_t loc, GAUGE::OPTION g);
284
292 void enrich (size_t loc, GAUGE::OPTION g, const vector<Biped<Symmetry,MatrixType> > &P);
293
306 template<typename MpoScalar>
307 ArrayXXcd intercellSF (const Mpo<Symmetry,MpoScalar> &Oalfa, const Mpo<Symmetry,MpoScalar> &Obeta, int Lx,
308 double kmin=0., double kmax=2.*M_PI, int kpoints=51,
309 DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::ON_EXIT, double tol=1e-12);
310
314 template<typename MpoScalar>
315 complex<Scalar> intercellSFpoint (const Mpo<Symmetry,MpoScalar> &Oalfa, const Mpo<Symmetry,MpoScalar> &Obeta, int Lx,
316 double kval,
317 DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::ON_EXIT);
318
330 template<typename MpoScalar>
331 ArrayXXcd SF (const ArrayXXcd &cellAvg, const vector<Mpo<Symmetry,MpoScalar> > &Oalfa, const vector<Mpo<Symmetry,MpoScalar> > &Obeta, int Lx,
332 double kmin, double kmax, int kpoints,
333 DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::ON_EXIT, double tol=1e-12);
334
338 template<typename MpoScalar>
339 complex<Scalar> SFpoint (const ArrayXXcd &cellAvg, const vector<Mpo<Symmetry,MpoScalar> > &Oalfa, const vector<Mpo<Symmetry,MpoScalar> > &Obeta, int Lx,
340 double kval,
341 DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::ON_EXIT);
342
343//private:
344
346 size_t N_sites;
347 size_t Mmax, Nqmax;
348 double eps_svd = 1e-13;
349 double eps_truncWeight = 1e-14;
350 size_t max_Nsv=100000ul, min_Nsv=1ul;
351 int max_Nrich;
352
353 qarray<Nq> Qtot;
354
356 void calc_entropy (size_t loc, bool PRINT=false);
357
359 void calc_entropy (bool PRINT=false) {for (size_t l=0; l<N_sites; ++l) calc_entropy(l,PRINT);};
360
362 ArrayXd truncWeight;
363
365 vector<vector<qarray<Symmetry::Nq> > > qloc;
366
368 std::array<vector<vector<Biped<Symmetry,MatrixType> > >,3> A; // A[L/R/C][l][s].block[q]
369
371 vector<Biped<Symmetry,MatrixType> > C; // zero-site part C[l]
372
374 // std::array<vector<vector<Biped<Symmetry,MatrixType> > >,3> N; // N[L/R/C][l][s].block[q]
375
376 VectorXd S;
377
378 vector<map<qarray<Nq>,tuple<ArrayXd,int> > > SVspec;
379
381 vector<Qbasis<Symmetry> > inbase;
382 vector<Qbasis<Symmetry> > outbase;
383
385 void update_inbase (size_t loc, GAUGE::OPTION g = GAUGE::C);
386 void update_outbase (size_t loc, GAUGE::OPTION g = GAUGE::C);
387 void update_inbase (GAUGE::OPTION g = GAUGE::C) {for (size_t l=0; l<this->N_sites; l++) {update_inbase (l,g);}}
388 void update_outbase (GAUGE::OPTION g = GAUGE::C) {for (size_t l=0; l<this->N_sites; l++) {update_outbase(l,g);}}
389};
390
391template<typename Symmetry, typename Scalar>
392string Umps<Symmetry,Scalar>::
393info() const
394{
395 stringstream ss;
396 ss << "Umps: ";
397 ss << Symmetry::name() << ", ";
398 //ss << ", " << Symmetry::name() << ", ";
399 if (Nq != 0)
400 {
401 ss << "(";
402 for (size_t q=0; q<Nq; ++q)
403 {
404 ss << Symmetry::kind()[q];
405 if (q!=Nq-1) {ss << ",";}
406 }
407 ss << ")=(" << Sym::format<Symmetry>(Qtot) << "), ";
408 }
409 ss << "Lcell=" << N_sites << ", ";
410 ss << "Mmax=" << calc_Mmax() << " (";
411 if (Symmetry::NON_ABELIAN)
412 {
413 ss << "full=" << calc_fullMmax() << ", ";
414 }
415 ss << "Dmax=" << calc_Dmax() << "), ";
416 ss << "Nqmax=" << calc_Nqmax() << ", ";
417 ss << "S=(" << S.transpose() << "), ";
418 ss << "mem=" << round(memory(GB),3) << "GB";
419
420 return ss.str();
421}
422
423template<typename Symmetry, typename Scalar>
424template<typename Hamiltonian>
425Umps<Symmetry,Scalar>::
426Umps (const Hamiltonian &H, qarray<Nq> Qtot_input, size_t L_input, size_t Mmax, size_t Nqmax, bool INIT_TO_HALF_INTEGER_SPIN)
427:N_sites(L_input), Qtot(Qtot_input)
428{
429 qloc = H.locBasis();
430 resize(Mmax,Nqmax,INIT_TO_HALF_INTEGER_SPIN);
431}
432
433template<typename Symmetry, typename Scalar>
434Umps<Symmetry,Scalar>::
435Umps (const vector<qarray<Symmetry::Nq> > &qloc_input, qarray<Nq> Qtot_input, size_t L_input, size_t Mmax, size_t Nqmax, bool INIT_TO_HALF_INTEGER_SPIN)
436:N_sites(L_input), Qtot(Qtot_input)
437{
438 qloc.resize(N_sites);
439 for (size_t l=0; l<N_sites; ++l) {qloc[l] = qloc_input;}
440 resize(Mmax,Nqmax,INIT_TO_HALF_INTEGER_SPIN);
441 ::transform_base<Symmetry>(qloc,Qtot); // from symmetry/functions.h
442}
443
444template<typename Symmetry, typename Scalar>
445Umps<Symmetry,Scalar>::
446Umps (const vector<vector<qarray<Symmetry::Nq> > > &qloc_input, qarray<Nq> Qtot_input, size_t L_input, size_t Mmax, size_t Nqmax, bool INIT_TO_HALF_INTEGER_SPIN)
447:N_sites(L_input), Qtot(Qtot_input)
448{
449 qloc.resize(N_sites);
450 for (size_t l=0; l<N_sites; ++l) {qloc[l] = qloc_input[l];}
451 resize(Mmax,Nqmax,INIT_TO_HALF_INTEGER_SPIN);
452 ::transform_base<Symmetry>(qloc,Qtot); // from symmetry/functions.h
453}
454
455template<typename Symmetry, typename Scalar>
456double Umps<Symmetry,Scalar>::
457memory (MEMUNIT memunit) const
458{
459 double res = 0.;
460 for (size_t l=0; l<N_sites; ++l)
461 {
462 res += C[l].memory(memunit);
463 for (size_t g=0; g<3; ++g)
464 for (size_t s=0; s<qloc[l].size(); ++s)
465 {
466 res += A[g][l][s].memory(memunit);
467 }
468 }
469 return res;
470}
471
472template<typename Symmetry, typename Scalar>
473size_t Umps<Symmetry,Scalar>::
474calc_Nqmax() const
475{
476 size_t res = 0;
477 for (size_t l=0; l<this->N_sites; ++l)
478 {
479 if (inbase[l].Nq() > res) {res = inbase[l].Nq();}
480 if (outbase[l].Nq() > res) {res = outbase[l].Nq();}
481 }
482 return res;
483}
484
485template<typename Symmetry, typename Scalar>
486size_t Umps<Symmetry,Scalar>::
487calc_Dmax() const
488{
489 size_t res = 0;
490 for (size_t l=0; l<this->N_sites; ++l)
491 {
492 if (inbase[l].Dmax() > res) {res = inbase[l].Dmax();}
493 if (outbase[l].Dmax() > res) {res = outbase[l].Dmax();}
494 }
495 return res;
496}
497
498template<typename Symmetry, typename Scalar>
499size_t Umps<Symmetry,Scalar>::
500calc_Mmax() const
501{
502 size_t res = 0;
503 for (size_t l=0; l<this->N_sites; ++l)
504 {
505 if (inbase[l].M() > res) {res = inbase[l].M();}
506 if (outbase[l].M() > res) {res = outbase[l].M();}
507 }
508 return res;
509}
510
511template<typename Symmetry, typename Scalar>
512std::size_t Umps<Symmetry,Scalar>::
513calc_fullMmax () const
514{
515 size_t res = 0;
516 for (size_t l=0; l<this->N_sites; ++l)
517 {
518 if (inbase[l].fullM() > res) {res = inbase[l].fullM();}
519 if (outbase[l].fullM() > res) {res = outbase[l].fullM();}
520 }
521 return res;
522}
523
524template<typename Symmetry, typename Scalar>
525void Umps<Symmetry,Scalar>::
526update_inbase (size_t loc, GAUGE::OPTION g)
527{
528 inbase[loc].clear();
529 inbase[loc].pullData(A[g][loc],0);
530}
531
532template<typename Symmetry, typename Scalar>
533void Umps<Symmetry,Scalar>::
534update_outbase (size_t loc, GAUGE::OPTION g)
535{
536 outbase[loc].clear();
537 outbase[loc].pullData(A[g][loc],1);
538}
539
540template<typename Symmetry, typename Scalar>
541qarray<Symmetry::Nq> Umps<Symmetry,Scalar>::
542Qtop (size_t loc) const
543{
544 return qplusinf<Symmetry::Nq>();
545}
546
547template<typename Symmetry, typename Scalar>
548qarray<Symmetry::Nq> Umps<Symmetry,Scalar>::
549Qbot (size_t loc) const
550{
551 return qminusinf<Symmetry::Nq>();
552}
553
554template<typename Symmetry, typename Scalar>
555void Umps<Symmetry,Scalar>::
556resize_arrays ()
557{
558 truncWeight.resize(N_sites);
559 for (size_t g=0; g<3; ++g)
560 {
561 A[g].resize(N_sites);
562 for (size_t l=0; l<N_sites; ++l)
563 {
564 A[g][l].resize(qloc[l].size());
565 }
566
567 // N[g].resize(N_sites);
568 // for (size_t l=0; l<N_sites; ++l)
569 // {
570 // N[g][l].resize(qloc[l].size());
571 // }
572 }
573 C.resize(N_sites);
574 inbase.resize(N_sites);
575 outbase.resize(N_sites);
576 S.resize(N_sites);
577 SVspec.resize(N_sites);
578}
579
580template<typename Symmetry, typename Scalar>
581template<typename OtherMatrixType>
582void Umps<Symmetry,Scalar>::
583resize (const Umps<Symmetry,OtherMatrixType> &V)
584{
585 N_sites = V.N_sites;
586 // N_phys = V.N_phys;
587 qloc = V.qloc;
588 Qtot = V.Qtot;
589
590 inbase = V.inbase;
591 outbase = V.outbase;
592
593 for (size_t g=0; g<3; g++) {A[g].resize(this->N_sites);}
594 C.resize(this->N_sites);
595
596 truncWeight.resize(this->N_sites); truncWeight.setZero();
597 S.resize(this->N_sites-1); S.setConstant(numeric_limits<double>::quiet_NaN());
598 SVspec.resize(N_sites);
599
600 for (size_t g=0; g<3; g++)
601 for (size_t l=0; l<V.N_sites; ++l)
602 {
603 A[g][l].resize(qloc[l].size());
604
605 for (size_t s=0; s<qloc[l].size(); ++s)
606 {
607 A[g][l][s].in = V.A[g][l][s].in;
608 A[g][l][s].out = V.A[g][l][s].out;
609 A[g][l][s].block.resize(V.A[g][l][s].dim);
610 A[g][l][s].dict = V.A[g][l][s].dict;
611 A[g][l][s].dim = V.A[g][l][s].dim;
612 }
613 C[l].in = V.C[l].in;
614 C[l].out = V.C[l].out;
615 C[l].block.resize(V.C[l].dim);
616 C[l].dict = V.C[l].dict;
617 C[l].dim = V.C[l].dim;
618 }
619}
620
621template<typename Symmetry, typename Scalar>
622void Umps<Symmetry,Scalar>::
623resize (size_t Mmax_input, size_t Nqmax_input, bool INIT_TO_HALF_INTEGER_SPIN)
624{
625// if (!Symmetry::NON_ABELIAN)
626// {
627 Mmax = Mmax_input;
628 Nqmax = Nqmax_input;
629 if (Symmetry::IS_TRIVIAL) {Nqmax = 1;}
630 else if (Symmetry::IS_MODULAR) {Nqmax = min(static_cast<size_t>(Symmetry::MOD_N),Nqmax_input);}
631
632 resize_arrays();
633
634 auto take_first_elems = [this] (const vector<qarray<Nq> > &qs) -> vector<qarray<Nq> >
635 {
636 vector<qarray<Nq> > out = qs;
637 // sort the vector first according to the distance to qvacuum
638 sort(out.begin(),out.end(),[this] (qarray<Nq> q1, qarray<Nq> q2)
639 {
640 VectorXd dist_q1(Nq);
641 VectorXd dist_q2(Nq);
642 for (size_t q=0; q<Nq; q++)
643 {
644 double Delta = 1.; // QinTop[loc][q] - QinBot[loc][q];
645 dist_q1(q) = q1[q] / Delta;
646 dist_q2(q) = q2[q] / Delta;
647 }
648 return (dist_q1.norm() < dist_q2.norm())? true:false;
649 });
650 return out;
651 };
652
653 vector<set<qarray<Symmetry::Nq> > > qinset(N_sites);
654 vector<set<qarray<Symmetry::Nq> > > qoutset(N_sites);
655 if (INIT_TO_HALF_INTEGER_SPIN)
656 {
657 for (const auto & q:Symmetry::lowest_qs()) { qoutset[N_sites-1].insert(q); }
658 }
659 else
660 {
661 qoutset[N_sites-1].insert(Symmetry::qvacuum());
662 }
663 ArrayXi inSize(N_sites); inSize = 0;
664 ArrayXi outSize(N_sites); outSize = 0;
665
666 while (not((inSize == Nqmax).all() and (outSize == Nqmax).all()))
667 {
668 for (size_t l=0; l<N_sites; ++l)
669 {
670 size_t index = (l==0)? N_sites-1 : (l-1)%N_sites;
671 for (const auto &t:qoutset[index])
672 {
673 if (qinset[l].size() < Nqmax)
674 {
675 qinset[l].insert(t);
676 }
677 }
678 inSize[l] = qinset[l].size();
679
680 vector<qarray<Symmetry::Nq> > qinvec(qinset[l].size());
681 copy(qinset[l].begin(), qinset[l].end(), qinvec.begin());
682
683 auto tmp = Symmetry::reduceSilent(qinvec, qloc[l], true);
684 tmp = take_first_elems(tmp);
685 for (const auto &t:tmp)
686 {
687 if (qoutset[l].size() < Nqmax)
688 {
689 qoutset[l].insert(t);
690 }
691 }
692 outSize[l] = qoutset[l].size();
693 }
694 }
695
696 // symmetrization
697 if (Qtot == Symmetry::qvacuum())
698 {
699 for (size_t l=0; l<N_sites; ++l)
700 {
701 auto qinset_tmp = qinset[l];
702 for (const auto &q:qinset_tmp)
703 {
704 if (auto it=qinset_tmp.find(Symmetry::flip(q)); it==qinset_tmp.end())
705 {
706 qinset[l].insert(Symmetry::flip(q));
707 }
708 }
709
710 auto qoutset_tmp = qoutset[l];
711 for (const auto &q:qoutset_tmp)
712 {
713 if (auto it=qoutset_tmp.find(Symmetry::flip(q)); it==qoutset_tmp.end())
714 {
715 qoutset[l].insert(Symmetry::flip(q));
716 }
717 }
718 }
719 }
720
721 for (size_t l=0; l<N_sites; ++l)
722 {
723 vector<qarray<Symmetry::Nq> > qins(qinset[l].size());
724 copy(qinset[l].begin(), qinset[l].end(), qins.begin());
725
726 vector<qarray<Symmetry::Nq> > qouts(qoutset[l].size());
727 copy(qoutset[l].begin(), qoutset[l].end(), qouts.begin());
728 assert(Mmax >= qins.size() and "Choose a greater Minit to have at least one state per QN block.");
729 size_t Dmax_in = Mmax / qins.size();
730 size_t Dmax_in_remainder = Mmax%qins.size();
731 assert(Mmax >= qouts.size() and "Choose a greater Minit to have at least one state per QN block.");
732 size_t Dmax_out = Mmax / qouts.size();
733 size_t Dmax_out_remainder = Mmax%qouts.size();
734
735 assert(Dmax_in*qins.size()+Dmax_in_remainder == Mmax and "Strange thing in Umps::resize");
736 assert(Dmax_out*qouts.size()+Dmax_out_remainder == Mmax and "Strange thing in Umps::resize");
737
738 MatrixXd Mtmpqq(Dmax_in,Dmax_out); Mtmpqq.setZero();
739 MatrixXd Mtmp0q(Dmax_in+Dmax_in_remainder,Dmax_out); Mtmp0q.setZero();
740 MatrixXd Mtmpq0(Dmax_in,Dmax_out+Dmax_out_remainder); Mtmpq0.setZero();
741 MatrixXd Mtmp00(Dmax_in+Dmax_in_remainder,Dmax_out+Dmax_out_remainder); Mtmp00.setZero();
742
743 for (size_t g=0; g<3; ++g)
744 for (size_t s=0; s<qloc[l].size(); ++s)
745 {
746 for (const auto &qin:qins)
747 {
748 auto qouts = Symmetry::reduceSilent(qloc[l][s], qin);
749 for (const auto &qout:qouts)
750 {
751 if (auto it=qoutset[l].find(qout); it!=qoutset[l].end())
752 {
753 qarray2<Symmetry::Nq> qinout = {qin,qout};
754 if (qin != Symmetry::qvacuum() and qout != Symmetry::qvacuum() ) {A[g][l][s].try_push_back(qinout, Mtmpqq);}
755 else if (qin != Symmetry::qvacuum() and qout == Symmetry::qvacuum() ) {A[g][l][s].try_push_back(qinout, Mtmpq0);}
756 else if (qin == Symmetry::qvacuum() and qout != Symmetry::qvacuum() ) {A[g][l][s].try_push_back(qinout, Mtmp0q);}
757 else if (qin == Symmetry::qvacuum() and qout == Symmetry::qvacuum() ) {A[g][l][s].try_push_back(qinout, Mtmp00);}
758 }
759 }
760 }
761 }
762 }
763
764 // for (size_t g=0; g<3; ++g)
765 // for (size_t l=0; l<N_sites; ++l)
766 // for (size_t s=0; s<qloc[l].size(); ++s)
767 // for (size_t q=0; q<A[g][l][s].dim; ++q)
768 // {
769 // A[g][l][s].block[q].resize(Dmax,Dmax);
770 // }
771
772 update_inbase();
773 update_outbase();
774
775 for (size_t l=0; l<N_sites; ++l)
776 {
777 size_t Dmax = Mmax / outbase[l].Nq();
778 size_t Dmax_remainder = Mmax%outbase[l].Nq();
779 assert(Dmax*outbase[l].Nq()+Dmax_remainder == Mmax and "Strange thing in Umps::resize");
780
781 MatrixXd Mtmpqq(Dmax,Dmax); Mtmpqq.setZero();
782 MatrixXd Mtmp00(Dmax+Dmax_remainder,Dmax+Dmax_remainder); Mtmp00.setZero();
783 for (size_t qout=0; qout<outbase[l].Nq(); ++qout)
784 {
785 if (outbase[l][qout] != Symmetry::qvacuum()) {C[l].try_push_back(qarray2<Symmetry::Nq>{outbase[l][qout], outbase[l][qout]}, Mtmpqq);}
786 else if (outbase[l][qout] == Symmetry::qvacuum()) {C[l].try_push_back(qarray2<Symmetry::Nq>{outbase[l][qout], outbase[l][qout]}, Mtmp00);}
787 }
788 }
789
790 for (size_t l=0; l<N_sites; ++l)
791 {
792 C[l] = C[l].sorted();
793 }
794 // graph("init");
795
796 setRandom();
797 // for (size_t l=0; l<N_sites; ++l) svdDecompose(l);
798// }
799// else
800// {
801// Dmax = Dmax_input;
802// Nqmax = Nqmax_input;
803// resize_arrays();
804//
805// Mps<Symmetry,Scalar> Tmp(N_sites, qloc, Symmetry::qvacuum(), N_sites, Nqmax, true);
806// Tmp.innerResize(Dmax);
807// Tmp.setRandom();
808//
809// A[GAUGE::C] = Tmp.A;
810// A[GAUGE::L] = Tmp.A;
811// A[GAUGE::R] = Tmp.A;
812//
813// for (size_t l=0; l<N_sites; ++l) Tmp.rightSplitStep(l,C[l]);
814//
815// normalize_C();
816//
817// update_inbase(GAUGE::C);
818// update_outbase(GAUGE::C);
819//
820// for (size_t l=0; l<N_sites; ++l) svdDecompose(l,GAUGE::C);
821// }
822}
823
824//template<typename Symmetry, typename Scalar>
825//void Umps<Symmetry,Scalar>::
826//resize (size_t Dmax_input, size_t Nqmax_input)
827//{
828// Dmax = Dmax_input;
829// Nqmax = Nqmax_input;
830// resize_arrays();
831//
832// Mps<Symmetry,Scalar> Tmp(N_sites, qloc, Symmetry::qvacuum(), N_sites, Nqmax, true);
833// Tmp.innerResize(Dmax);
834// Tmp.setRandom();
835//
836// A[GAUGE::C] = Tmp.A;
837// A[GAUGE::L] = Tmp.A;
838// A[GAUGE::R] = Tmp.A;
839//
840// for (size_t l=0; l<N_sites; ++l) Tmp.rightSplitStep(l,C[l]);
841//
842// normalize_C();
843//
844// update_inbase(GAUGE::C);
845// update_outbase(GAUGE::C);
846//
847// for (size_t l=0; l<N_sites; ++l) svdDecompose(l,GAUGE::C);
848//}
849
850template<typename Symmetry, typename Scalar>
851void Umps<Symmetry,Scalar>::
852normalize_C()
853{
854 // normalize the centre matrices for proper wavefunction norm: Tr(C*C†)=1
855 for (size_t l=0; l<N_sites; ++l)
856 {
857 C[l] = 1./sqrt((C[l].contract(C[l].adjoint())).trace()) * C[l];
858 }
859}
860
861template<typename Symmetry, typename Scalar>
862void Umps<Symmetry,Scalar>::
863setRandom()
864{
865 for (size_t l=0; l<N_sites; ++l)
866 for (size_t q=0; q<C[l].dim; ++q)
867 {
868 for (size_t a1=0; a1<C[l].block[q].rows(); ++a1)
869// for (size_t a2=0; a2<C[l].block[0].cols(); ++a2)
870 for (size_t a2=0; a2<=a1; ++a2)
871 {
872 C[l].block[q](a1,a2) = threadSafeRandUniform<Scalar,double>(-1.,1.);
873 C[l].block[q](a2,a1) = C[l].block[q](a1,a2);
874 }
875 }
876
877 normalize_C();
878
879 for (size_t l=0; l<N_sites; ++l)
880 for (size_t s=0; s<qloc[l].size(); ++s)
881 for (size_t q=0; q<A[GAUGE::C][l][s].dim; ++q)
882 for (size_t a1=0; a1<A[GAUGE::C][l][s].block[q].rows(); ++a1)
883 for (size_t a2=0; a2<A[GAUGE::C][l][s].block[q].cols(); ++a2)
884 // for (size_t a2=0; a2<=a1; ++a2)
885 {
886 A[GAUGE::C][l][s].block[q](a1,a2) = threadSafeRandUniform<Scalar,double>(-1.,1.);
887 // A[GAUGE::C][l][s].block[q](a2,a1) = A[GAUGE::C][l][s].block[q](a1,a2);
888 }
889
890 for (size_t l=0; l<N_sites; ++l)
891 {
892 svdDecompose(l);
893 }
894
895 calc_entropy();
896}
897
898template<typename Symmetry, typename Scalar>
899void Umps<Symmetry,Scalar>::
900graph (string filename) const
901{
902 stringstream ss;
903
904 ss << "#!/usr/bin/dot dot -Tpdf -o " << filename << ".pdf\n\n";
905 ss << "digraph G\n{\n";
906 ss << "rankdir = LR;\n";
907 ss << "labelloc=\"t\";\n";
908 ss << "label=\"Umps: cell size=" << N_sites << ", Q=(";
909 for (size_t q=0; q<Nq; ++q)
910 {
911 ss << Symmetry::kind()[q];
912 if (q!=Nq-1) {ss << ",";}
913 }
914 ss << ")" << "\";\n";
915
916 // first node
917 ss << "\"l=" << 0 << ", " << Sym::format<Symmetry>(Symmetry::qvacuum()) << "\"";
918 ss << "[label=" << "\"" << Sym::format<Symmetry>(Symmetry::qvacuum()) << "\"" << "];\n";
919
920 // site nodes
921 for (size_t l=0; l<N_sites; ++l)
922 {
923 ss << "subgraph" << " cluster_" << l << "\n{\n";
924 for (size_t s=0; s<qloc[l].size(); ++s)
925 for (size_t q=0; q<A[GAUGE::C][l][s].dim; ++q)
926 {
927 string qin = Sym::format<Symmetry>(A[GAUGE::C][l][s].in[q]);
928 ss << "\"l=" << l << ", " << qin << "\"";
929 ss << "[label=" << "\"" << qin << "\"" << "];\n";
930 }
931 ss << "label=\"l=" << l << "\"\n";
932 ss << "}\n";
933 }
934
935 // last node
936 ss << "subgraph" << " cluster_" << N_sites << "\n{\n";
937 for (size_t s=0; s<qloc[N_sites-1].size(); ++s)
938 for (size_t q=0; q<A[GAUGE::C][N_sites-1][s].dim; ++q)
939 {
940 string qout = Sym::format<Symmetry>(A[GAUGE::C][N_sites-1][s].out[q]);
941 ss << "\"l=" << N_sites << ", " << qout << "\"";
942 ss << "[label=" << "\"" << qout << "\"" << "];\n";
943 }
944 ss << "label=\"l=" << N_sites << "=0" << "\"\n";
945 ss << "}\n";
946
947 // edges
948 for (size_t l=0; l<N_sites; ++l)
949 {
950 for (size_t s=0; s<qloc[l].size(); ++s)
951 for (size_t q=0; q<A[GAUGE::C][l][s].dim; ++q)
952 {
953 string qin = Sym::format<Symmetry>(A[GAUGE::C][l][s].in[q]);
954 string qout = Sym::format<Symmetry>(A[GAUGE::C][l][s].out[q]);
955 ss << "\"l=" << l << ", " << qin << "\"";
956 ss << "->";
957 ss << "\"l=" << l+1 << ", " << qout << "\"";
958 ss << " [label=\"" << A[GAUGE::C][l][s].block[q].rows() << "x" << A[GAUGE::C][l][s].block[q].cols() << "\"";
959 ss << "];\n";
960 }
961 }
962
963 ss << "\n}";
964
965 ofstream f(filename+".dot");
966 f << ss.str();
967 f.close();
968}
969
970template<typename Symmetry, typename Scalar>
971string Umps<Symmetry,Scalar>::
972test_ortho (double tol) const
973{
974 stringstream sout;
975 std::array<string,4> normal_token = {"A","B","M","X"};
976 std::array<string,4> special_token = {"\e[4mA\e[0m","\e[4mB\e[0m","\e[4mM\e[0m","\e[4mX\e[0m"};
977 Array<Scalar,Dynamic,1> norm(N_sites);
978
979 for (int l=0; l<N_sites; ++l)
980 {
981 // check for A
982 Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > Test = A[GAUGE::L][l][0].adjoint().contract(A[GAUGE::L][l][0]);
983 for (size_t s=1; s<qloc[l].size(); ++s)
984 {
985 Test += A[GAUGE::L][l][s].adjoint().contract(A[GAUGE::L][l][s]);
986 }
987 // cout << "AL test:" << endl << Test.print(true) << endl;
988 vector<bool> A_CHECK(Test.dim);
989 vector<double> A_infnorm(Test.dim);
990 for (size_t q=0; q<Test.dim; ++q)
991 {
992 Test.block[q] -= MatrixType::Identity(Test.block[q].rows(), Test.block[q].cols());
993 A_CHECK[q] = Test.block[q].norm()<tol ? true : false;
994 A_infnorm[q] = Test.block[q].norm();
995 // cout << "q=" << Test.in[q] << ", A_infnorm[q]=" << A_infnorm[q] << endl;
996 }
997
998 // check for B
999 Test.clear();
1000 Test = A[GAUGE::R][l][0].contract(A[GAUGE::R][l][0].adjoint(), contract::MODE::OORR);
1001 for (size_t s=1; s<qloc[l].size(); ++s)
1002 {
1003 Test += A[GAUGE::R][l][s].contract(A[GAUGE::R][l][s].adjoint(), contract::MODE::OORR);
1004 }
1005 // cout << Test.print(true) << endl;
1006 vector<bool> B_CHECK(Test.dim);
1007 vector<double> B_infnorm(Test.dim);
1008 for (size_t q=0; q<Test.dim; ++q)
1009 {
1010 Test.block[q] -= MatrixType::Identity(Test.block[q].rows(), Test.block[q].cols());
1011 B_CHECK[q] = Test.block[q].template lpNorm<Infinity>()<tol ? true : false;
1012 B_infnorm[q] = Test.block[q].template lpNorm<Infinity>();
1013// cout << "q=" << Test.in[q] << ", B_infnorm[q]=" << B_infnorm[q] << endl;
1014 }
1015
1016 // check for AL*C = AC
1017 for (size_t s=0; s<qloc[l].size(); ++s)
1018 {
1019 Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > Test = A[GAUGE::L][l][s] * C[l];
1020 for (size_t q=0; q<Test.dim; ++q)
1021 {
1022 qarray2<Symmetry::Nq> quple = {Test.in[q], Test.out[q]};
1023 auto it = A[GAUGE::C][l][s].dict.find(quple);
1024 if (it != A[GAUGE::C][l][s].dict.end())
1025 {
1026 Test.block[q] -= A[GAUGE::C][l][s].block[it->second];
1027 }
1028 }
1029 vector<double> T_CHECK(Test.dim);
1030// cout << "Test.dim=" << Test.dim << endl;
1031 double normsum = 0;
1032 for (size_t q=0; q<Test.dim; ++q)
1033 {
1034// cout << "g=L, " << "s=" << s << ", q=" << Test.in[q] << ", " << Test.out[q]
1035// << ", norm=" << Test.block[q].template lpNorm<Infinity>() << endl;
1036 normsum += Test.block[q].norm();
1037 T_CHECK[q] = Test.block[q].norm()<tol ? true : false;
1038 }
1039 if (all_of(T_CHECK.begin(),T_CHECK.end(),[](bool x){return x;}))
1040 {
1041 cout << "l=" << l << ", s=" << s << ", AL[" << l << "]*C[" << l << "]=AC[" << l << "]="
1042 << termcolor::green << "true" << termcolor::reset << ", normsum=" << normsum << endl;
1043 }
1044 else
1045 {
1046 cout << "l=" << l << ", s=" << s << ", AL[" << l << "]*C[" << l << "]=AC[" << l << "]="
1047 << termcolor::red << "false" << termcolor::reset << ", normsum=" << normsum << endl;
1048 }
1049 }
1050
1051 // check for C*AR = AC
1052 size_t locC = minus1modL(l);
1053 for (size_t s=0; s<qloc[l].size(); ++s)
1054 {
1055 Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > Test = C[locC] * A[GAUGE::R][l][s];
1056 for (size_t q=0; q<Test.dim; ++q)
1057 {
1058 qarray2<Symmetry::Nq> quple = {Test.in[q], Test.out[q]};
1059 auto it = A[GAUGE::C][l][s].dict.find(quple);
1060 if (it != A[GAUGE::C][l][s].dict.end())
1061 {
1062 Test.block[q] -= A[GAUGE::C][l][s].block[it->second];
1063 }
1064 }
1065 vector<double> T_CHECK(Test.dim);
1066// cout << "Test.dim=" << Test.dim << endl;
1067 double normsum = 0;
1068 for (size_t q=0; q<Test.dim; ++q)
1069 {
1070// cout << "g=R, " << "s=" << s << ", q=" << Test.in[q] << ", " << Test.out[q]
1071// << ", norm=" << Test.block[q].template lpNorm<Infinity>() << endl;
1072 normsum += Test.block[q].template lpNorm<Infinity>();
1073 T_CHECK[q] = Test.block[q].template lpNorm<Infinity>()<tol ? true : false;
1074 }
1075 if (all_of(T_CHECK.begin(),T_CHECK.end(),[](bool x){return x;}))
1076 {
1077 cout << "l=" << l << ", s=" << s << ", C[" << locC << "]*AR[" << l << "]=AC[" << l << "]="
1078 << termcolor::green << "true" << termcolor::reset << ", normsum=" << normsum << endl;
1079 }
1080 else
1081 {
1082 cout << "l=" << l << ", s=" << s << ", C[" << locC << "]*AR[" << l << "]=AC[" << l << "]="
1083 << termcolor::red << "false" << termcolor::reset << ", normsum=" << normsum << endl;
1084 }
1085 }
1086
1087 norm(l) = (C[l].contract(C[l].adjoint())).trace();
1088
1089 // interpret result
1090 if (all_of(A_CHECK.begin(),A_CHECK.end(),[](bool x){return x;}))
1091 {
1092 sout << termcolor::red;
1093 sout << normal_token[0]; // A
1094 }
1095 else
1096 {
1097 // assert(1!=1 and "AL is wrong");
1098 sout << termcolor::green;
1099 sout << normal_token[2]; // M
1100 }
1101
1102 if (all_of(B_CHECK.begin(),B_CHECK.end(),[](bool x){return x;}))
1103 {
1104 sout << termcolor::blue;
1105 sout << normal_token[1]; // B
1106 }
1107 else
1108 {
1109 // assert(1!=1 and "AR is wrong");
1110 sout << termcolor::green;
1111 sout << normal_token[2]; // M
1112 }
1113 }
1114
1115 sout << termcolor::reset;
1116 sout << ", norm=" << norm.transpose();
1117 return sout.str();
1118}
1119
1120template<typename Symmetry, typename Scalar>
1121double Umps<Symmetry,Scalar>::
1122dot (const Umps<Symmetry,Scalar> &Vket) const
1123{
1124// double outL = calc_LReigen(VMPS::DIRECTION::LEFT, A[GAUGE::L], Vket.A[GAUGE::L], inBasis(0), Vket.inBasis(0), qloc).energy;
1125// // for testing:
1126// double outR = calc_LReigen(VMPS::DIRECTION::RIGHT, A[GAUGE::R], Vket.A[GAUGE::R], outBasis(N_sites-1), Vket.outBasis(N_sites-1), qloc).energy;
1127// lout << "dot consistency check: using AL: " << outL << ", using AR: " << outR << endl;
1128// return outL;
1129
1130 Eigenstate<Biped<Symmetry,Matrix<complex<Scalar>,Dynamic,Dynamic> > > Leigen, Reigen;
1131
1132 #pragma omp parallel sections
1133 {
1134 #pragma omp section
1135 {
1136 Leigen = calc_LReigen(VMPS::DIRECTION::LEFT, A[GAUGE::L], Vket.A[GAUGE::L], inBasis(0), Vket.inBasis(0), qloc);
1137 }
1138 #pragma omp section
1139 {
1140 Reigen = calc_LReigen(VMPS::DIRECTION::RIGHT, A[GAUGE::R], Vket.A[GAUGE::R], outBasis(N_sites-1), Vket.outBasis(N_sites-1), qloc);
1141 }
1142 }
1143
1144 auto LxCket = Leigen.state.contract(Vket.C[N_sites-1].template cast<MatrixXcd>());
1145 auto RxCdag = Reigen.state.contract( C[N_sites-1].adjoint().template cast<MatrixXcd>());
1146 auto mixed = RxCdag.contract(LxCket).trace();
1147
1148 lout << "dot: L gauge: " << Leigen.energy
1149 << ", R gauge: " << Reigen.energy
1150 << ", diff=" << abs(Leigen.energy-Reigen.energy)
1151 << ", mixed gauge (?): " << mixed << endl;
1152 return Leigen.energy;
1153}
1154
1155template<typename Symmetry, typename Scalar>
1156void Umps<Symmetry,Scalar>::
1157calc_entropy (size_t loc, bool PRINT)
1158{
1159 S(loc) = 0;
1160 SVspec[loc].clear();
1161
1162 if (PRINT)
1163 {
1164 lout << termcolor::magenta << "loc=" << loc << termcolor::reset << endl;
1165 }
1166 for (size_t q=0; q<C[loc].dim; ++q)
1167 {
1168 #ifdef DONT_USE_BDCSVD
1169 JacobiSVD<MatrixType> Jack; // standard SVD
1170 #else
1171 BDCSVD<MatrixType> Jack; // "Divide and conquer" SVD (only available in Eigen)
1172 #endif
1173
1174 Jack.compute(C[loc].block[q], ComputeThinU|ComputeThinV);
1175// Csingular[loc] += Jack.singularValues();
1176
1177 size_t Nnz = (Jack.singularValues().array() > 0.).count();
1178 double Scontrib = -Symmetry::degeneracy(C[loc].in[q]) *
1179 (Jack.singularValues().head(Nnz).array().square() *
1180 Jack.singularValues().head(Nnz).array().square().log()
1181 ).sum();
1182
1183// lout << "loc=" << loc << ", q=" << q
1184// << ", C[loc].in[q]=" << Sym::format<Symmetry>(C[loc].in[q])
1185// << ", Scontrib=" << Scontrib
1186// << ", deg=" << Symmetry::degeneracy(C[loc].in[q])
1187// << endl;
1188//
1189// for (int i=0; i<Nnz; ++i)
1190// {
1191// lout << "i=" << i << ", "
1192// << pow(Jack.singularValues()(i),2) << ", "
1193// << log(pow(Jack.singularValues()(i),2)) << ", "
1194// << -pow(Jack.singularValues()(i),2) * log(pow(Jack.singularValues()(i),2))
1195// << endl;
1196// }
1197// lout << endl;
1198 S(loc) += Scontrib;
1199
1200 SVspec[loc].insert(pair<qarray<Symmetry::Nq>,tuple<ArrayXd,int> >(C[loc].in[q], make_tuple(Jack.singularValues(), Symmetry::degeneracy(C[loc].in[q]))));
1201
1202 if (PRINT)
1203 {
1204 lout << termcolor::magenta
1205 << "S(" << C[loc].in[q] << ","
1206 << C[loc].out[q] << ")=" << Scontrib
1207 << ", size=" << C[loc].block[q].rows() << "x" << C[loc].block[q].cols()
1208 << ", deg=" << Symmetry::degeneracy(C[loc].in[q])
1209 << ", #sv=" << Jack.singularValues().rows()
1210 << ", svs=" << Jack.singularValues().head(min(20,int(Jack.singularValues().rows()))).transpose()
1211 << termcolor::reset << endl;
1212 }
1213 }
1214 if (PRINT)
1215 {
1216 lout << endl;
1217 }
1218
1219// lout << termcolor::blue << "S=" << S << termcolor::reset << endl;
1220}
1221
1222template<typename Symmetry, typename Scalar>
1223Scalar Umps<Symmetry,Scalar>::
1224calc_epsLRsq (GAUGE::OPTION gauge, size_t loc) const
1225{
1226 for (size_t s=0; s<qloc[loc].size(); ++s)
1227 for (size_t q=0; q<A[GAUGE::C][loc][s].dim; ++q)
1228 {
1229 std::array<qarray2<Symmetry::Nq>,3> quple;
1230 for (size_t g=0; g<3; ++g)
1231 {
1232 quple[g] = {A[g][loc][s].in[q], A[g][loc][s].out[q]};
1233 }
1234
1235 for (size_t g=1; g<3; ++g)
1236 {
1237 assert(quple[0] == quple[g]);
1238 }
1239 }
1240
1241 Scalar res = 0;
1242
1243 if (gauge == GAUGE::L)
1244 {
1245 for (size_t qout=0; qout<outbase[loc].Nq(); ++qout)
1246 {
1247 qarray2<Symmetry::Nq> quple = {outbase[loc][qout], outbase[loc][qout]};
1248 auto it = C[loc].dict.find(quple);
1249 assert(it != C[loc].dict.end());
1250 size_t qC = it->second;
1251
1252 // Determine how many A's to glue together
1253 vector<size_t> svec, qvec, Nrowsvec;
1254 for (size_t s=0; s<qloc[loc].size(); ++s)
1255 for (size_t q=0; q<A[GAUGE::C][loc][s].dim; ++q)
1256 {
1257 if (A[GAUGE::C][loc][s].out[q] == outbase[loc][qout])
1258 {
1259 svec.push_back(s);
1260 qvec.push_back(q);
1261 Nrowsvec.push_back(A[GAUGE::C][loc][s].block[q].rows());
1262 }
1263 }
1264
1265 // Do the glue
1266 size_t Ncols = A[GAUGE::C][loc][svec[0]].block[qvec[0]].cols();
1267 for (size_t i=1; i<svec.size(); ++i) {assert(A[GAUGE::C][loc][svec[i]].block[qvec[i]].cols() == Ncols);}
1268 size_t Nrows = accumulate(Nrowsvec.begin(),Nrowsvec.end(),0);
1269
1270 MatrixType Aclump(Nrows,Ncols);
1271 MatrixType Acmp(Nrows,Ncols);
1272 Aclump.setZero();
1273 Acmp.setZero();
1274 size_t stitch = 0;
1275 for (size_t i=0; i<svec.size(); ++i)
1276 {
1277 Aclump.block(stitch,0, Nrowsvec[i],Ncols) = A[GAUGE::C][loc][svec[i]].block[qvec[i]];
1278 Acmp.block (stitch,0, Nrowsvec[i],Ncols) = A[GAUGE::L][loc][svec[i]].block[qvec[i]];
1279 stitch += Nrowsvec[i];
1280 }
1281
1282 double diff = (Aclump-Acmp*C[loc].block[qC]).squaredNorm() * Symmetry::coeff_dot(C[loc].in[qC]);
1283 double summ = (Aclump+Acmp*C[loc].block[qC]).squaredNorm() * Symmetry::coeff_dot(C[loc].in[qC]);
1284// res += (Aclump-Acmp*C[loc].block[qC]).squaredNorm() * Symmetry::coeff_dot(C[loc].in[qC]);
1285 res += min(diff,summ);
1286// cout << "L, loc=" << loc
1287// << ", qout=" << qout
1288// << ", diff=" << (Aclump-Acmp*C[loc].block[qC]).squaredNorm() * Symmetry::coeff_dot(C[loc].in[qC])
1289// << ", summ=" << (Aclump+Acmp*C[loc].block[qC]).squaredNorm() * Symmetry::coeff_dot(C[loc].in[qC])
1290// << endl;
1291 }
1292 }
1293 else if (gauge == GAUGE::R)
1294 {
1295 size_t locC = minus1modL(loc);
1296
1297 for (size_t qin=0; qin<inbase[loc].Nq(); ++qin)
1298 {
1299 qarray2<Symmetry::Nq> quple = {inbase[loc][qin], inbase[loc][qin]};
1300 auto it = C[locC].dict.find(quple);
1301 assert(it != C[locC].dict.end());
1302 size_t qC = it->second;
1303
1304 // Determine how many A's to glue together
1305 vector<size_t> svec, qvec, Ncolsvec;
1306 for (size_t s=0; s<qloc[loc].size(); ++s)
1307 for (size_t q=0; q<A[GAUGE::C][loc][s].dim; ++q)
1308 {
1309 if (A[GAUGE::C][loc][s].in[q] == inbase[loc][qin])
1310 {
1311 svec.push_back(s);
1312 qvec.push_back(q);
1313 Ncolsvec.push_back(A[GAUGE::C][loc][s].block[q].cols());
1314 }
1315 }
1316
1317 // Do the glue
1318 size_t Nrows = A[GAUGE::C][loc][svec[0]].block[qvec[0]].rows();
1319 for (size_t i=1; i<svec.size(); ++i) {assert(A[GAUGE::C][loc][svec[i]].block[qvec[i]].rows() == Nrows);}
1320 size_t Ncols = accumulate(Ncolsvec.begin(), Ncolsvec.end(), 0);
1321
1322 MatrixType Aclump(Nrows,Ncols);
1323 MatrixType Acmp(Nrows,Ncols);
1324 Aclump.setZero();
1325 Acmp.setZero();
1326 size_t stitch = 0;
1327 for (size_t i=0; i<svec.size(); ++i)
1328 {
1329 Aclump.block(0,stitch, Nrows,Ncolsvec[i]) = A[GAUGE::C][loc][svec[i]].block[qvec[i]]*
1330 Symmetry::coeff_leftSweep(
1331 A[GAUGE::C][loc][svec[i]].in[qvec[i]],
1332 A[GAUGE::C][loc][svec[i]].out[qvec[i]]);
1333
1334 Acmp.block (0,stitch, Nrows,Ncolsvec[i]) = A[GAUGE::R][loc][svec[i]].block[qvec[i]]*
1335 Symmetry::coeff_leftSweep(
1336 A[GAUGE::R][loc][svec[i]].in[qvec[i]],
1337 A[GAUGE::R][loc][svec[i]].out[qvec[i]]);
1338 stitch += Ncolsvec[i];
1339 }
1340
1341 double diff = (Aclump-C[locC].block[qC]*Acmp).squaredNorm() * Symmetry::coeff_dot(C[locC].in[qC]);
1342 double summ = (Aclump+C[locC].block[qC]*Acmp).squaredNorm() * Symmetry::coeff_dot(C[locC].in[qC]);
1343// res += (Aclump-C[locC].block[qC]*Acmp).squaredNorm() * Symmetry::coeff_dot(C[locC].in[qC]);
1344 res += min(diff,summ);
1345// cout << "R, loc=" << loc
1346// << ", qin=" << qin
1347// << ", diff=" << (Aclump-C[locC].block[qC]*Acmp).squaredNorm() * Symmetry::coeff_dot(C[locC].in[qC])
1348// << ", summ=" << (Aclump+C[locC].block[qC]*Acmp).squaredNorm() * Symmetry::coeff_dot(C[locC].in[qC])
1349// << endl;
1350 }
1351 }
1352
1353 return res;
1354}
1355
1356template<typename Symmetry, typename Scalar>
1357void Umps<Symmetry,Scalar>::
1358polarDecompose (size_t loc, GAUGE::OPTION gauge)
1359{
1360 // check that blocks are the same for all gauges
1361 for (size_t s=0; s<qloc[loc].size(); ++s)
1362 for (size_t q=0; q<A[GAUGE::C][loc][s].dim; ++q)
1363 {
1364 std::array<qarray2<Symmetry::Nq>,3> quple;
1365 for (size_t g=0; g<3; ++g)
1366 {
1367 quple[g] = {A[g][loc][s].in[q], A[g][loc][s].out[q]};
1368 }
1369
1370 for (size_t g=1; g<3; ++g)
1371 {
1372 if (quple[0] != quple[g])
1373 {
1374 cout << "g=" << g << ", quple[0]=(" << quple[0][0] << "," << quple[0][1] << ")" << ", quple[g]=(" << quple[g][0] << "," << quple[g][1] << ")" << endl;
1375 }
1376 assert(quple[0] == quple[g]);
1377 }
1378 }
1379
1380 #ifdef DONT_USE_BDCSVD
1381 JacobiSVD<MatrixType> Jack; // standard SVD
1382 #else
1383 BDCSVD<MatrixType> Jack; // "Divide and conquer" SVD (only available in Eigen)
1384 #endif
1385
1386 if (gauge == GAUGE::L or gauge == GAUGE::C)
1387 {
1388// S(loc) = 0;
1389 vector<MatrixType> UC(C[loc].dim);
1390 for (size_t q=0; q<C[loc].dim; ++q)
1391 {
1392 Jack.compute(C[loc].block[q], ComputeThinU|ComputeThinV);
1393 UC[q] = Jack.matrixU() * Jack.matrixV().adjoint();
1394 }
1395
1396 for (size_t qout=0; qout<outbase[loc].Nq(); ++qout)
1397 {
1398 qarray2<Symmetry::Nq> quple = {outbase[loc][qout], outbase[loc][qout]};
1399
1400 // Determine how many A's to glue together
1401 vector<size_t> svec, qvec, Nrowsvec;
1402 for (size_t s=0; s<qloc[loc].size(); ++s)
1403 for (size_t q=0; q<A[GAUGE::C][loc][s].dim; ++q)
1404 {
1405 if (A[GAUGE::C][loc][s].out[q] == outbase[loc][qout])
1406 {
1407 svec.push_back(s);
1408 qvec.push_back(q);
1409 Nrowsvec.push_back(A[GAUGE::C][loc][s].block[q].rows());
1410 }
1411 }
1412
1413 // Do the glue
1414 size_t Ncols = A[GAUGE::C][loc][svec[0]].block[qvec[0]].cols();
1415 for (size_t i=1; i<svec.size(); ++i) {assert(A[GAUGE::C][loc][svec[i]].block[qvec[i]].cols() == Ncols);}
1416 size_t Nrows = accumulate(Nrowsvec.begin(),Nrowsvec.end(),0);
1417
1418 MatrixType Aclump(Nrows,Ncols);
1419 Aclump.setZero();
1420 size_t stitch = 0;
1421 for (size_t i=0; i<svec.size(); ++i)
1422 {
1423 Aclump.block(stitch,0, Nrowsvec[i],Ncols) = A[GAUGE::C][loc][svec[i]].block[qvec[i]];
1424 stitch += Nrowsvec[i];
1425 }
1426
1427 Jack.compute(Aclump,ComputeThinU|ComputeThinV);
1428 MatrixType UL = Jack.matrixU() * Jack.matrixV().adjoint();
1429
1430 auto it = C[loc].dict.find(quple);
1431 assert(it != C[loc].dict.end());
1432 size_t qC = it->second;
1433
1434 // Update AL
1435 stitch = 0;
1436 for (size_t i=0; i<svec.size(); ++i)
1437 {
1438 A[GAUGE::L][loc][svec[i]].block[qvec[i]] = UL.block(stitch,0, Nrowsvec[i],Ncols) * UC[qC].adjoint();
1439 stitch += Nrowsvec[i];
1440 }
1441 }
1442 }
1443
1444 size_t locC = minus1modL(loc);
1445
1446 if (gauge == GAUGE::R or gauge == GAUGE::C)
1447 {
1448 vector<MatrixType> UC(C[locC].dim);
1449// cout << "polarDecompose AR from C at " << locC << endl;
1450
1451 for (size_t q=0; q<C[locC].dim; ++q)
1452 {
1453 Jack.compute(C[locC].block[q], ComputeThinU|ComputeThinV);
1454 UC[q] = Jack.matrixU() * Jack.matrixV().adjoint();
1455 }
1456
1457 for (size_t qin=0; qin<inbase[loc].Nq(); ++qin)
1458 {
1459 qarray2<Symmetry::Nq> quple = {inbase[loc][qin], inbase[loc][qin]};
1460
1461 // Determine how many A's to glue together
1462 vector<size_t> svec, qvec, Ncolsvec;
1463 for (size_t s=0; s<qloc[loc].size(); ++s)
1464 for (size_t q=0; q<A[GAUGE::C][loc][s].dim; ++q)
1465 {
1466 if (A[GAUGE::C][loc][s].in[q] == inbase[loc][qin])
1467 {
1468 svec.push_back(s);
1469 qvec.push_back(q);
1470 Ncolsvec.push_back(A[GAUGE::C][loc][s].block[q].cols());
1471 }
1472 }
1473
1474 // Do the glue
1475 size_t Nrows = A[GAUGE::C][loc][svec[0]].block[qvec[0]].rows();
1476 for (size_t i=1; i<svec.size(); ++i) {assert(A[GAUGE::C][loc][svec[i]].block[qvec[i]].rows() == Nrows);}
1477 size_t Ncols = accumulate(Ncolsvec.begin(), Ncolsvec.end(), 0);
1478
1479 MatrixType Aclump(Nrows,Ncols);
1480 size_t stitch = 0;
1481 for (size_t i=0; i<svec.size(); ++i)
1482 {
1483 Aclump.block(0,stitch, Nrows,Ncolsvec[i]) = A[GAUGE::C][loc][svec[i]].block[qvec[i]]*
1484 Symmetry::coeff_leftSweep(
1485 A[GAUGE::C][loc][svec[i]].out[qvec[i]],
1486 A[GAUGE::C][loc][svec[i]].in[qvec[i]]);
1487 stitch += Ncolsvec[i];
1488 }
1489
1490 Jack.compute(Aclump,ComputeThinU|ComputeThinV);
1491 MatrixType UR = Jack.matrixU() * Jack.matrixV().adjoint();
1492
1493 auto it = C[locC].dict.find(quple);
1494 assert(it != C[locC].dict.end());
1495 size_t qC = it->second;
1496
1497 // Update AR
1498 stitch = 0;
1499 for (size_t i=0; i<svec.size(); ++i)
1500 {
1501 A[GAUGE::R][loc][svec[i]].block[qvec[i]] = UC[qC].adjoint() * UR.block(0,stitch, Nrows,Ncolsvec[i])*
1502 Symmetry::coeff_leftSweep(
1503 A[GAUGE::C][loc][svec[i]].in[qvec[i]],
1504 A[GAUGE::C][loc][svec[i]].out[qvec[i]]);
1505 stitch += Ncolsvec[i];
1506 }
1507 }
1508 }
1509}
1510
1511template<typename Symmetry, typename Scalar>
1512void Umps<Symmetry,Scalar>::
1513svdDecompose (size_t loc, GAUGE::OPTION gauge)
1514{
1515 // check that blocks are the same for all gauges
1516 for (size_t s=0; s<qloc[loc].size(); ++s)
1517 for (size_t q=0; q<A[GAUGE::C][loc][s].dim; ++q)
1518 {
1519 std::array<qarray2<Symmetry::Nq>,3> quple;
1520 for (size_t g=0; g<3; ++g)
1521 {
1522 quple[g] = {A[g][loc][s].in[q], A[g][loc][s].out[q]};
1523 }
1524
1525 for (size_t g=1; g<3; ++g)
1526 {
1527 assert(quple[0] == quple[g]);
1528 }
1529 }
1530
1531 if (gauge == GAUGE::L or gauge == GAUGE::C)
1532 {
1533 ArrayXd truncWeightSub(outbase[loc].Nq()); truncWeightSub.setZero();
1534
1535 vector<Biped<Symmetry,MatrixType> > Atmp(qloc[loc].size());
1536 for (size_t s=0; s<qloc[loc].size(); ++s)
1537 {
1538 Atmp[s] = A[GAUGE::C][loc][s] * C[loc].adjoint();
1539 }
1540
1541// cout << "svdDecompose AL from C at " << loc << endl;
1542 for (size_t qout=0; qout<outbase[loc].Nq(); ++qout)
1543 {
1544 // qarray2<Symmetry::Nq> quple = {outbase[loc][qout], outbase[loc][qout]};
1545 // auto it = C[loc].dict.find(quple);
1546 // assert(it != C[loc].dict.end());
1547 // size_t qC = it->second;
1548
1549 // Determine how many A's to glue together
1550 vector<size_t> svec, qvec, Nrowsvec;
1551 for (size_t s=0; s<qloc[loc].size(); ++s)
1552 for (size_t q=0; q<Atmp[s].dim; ++q)
1553 {
1554 if (Atmp[s].out[q] == outbase[loc][qout])
1555 {
1556 svec.push_back(s);
1557 qvec.push_back(q);
1558 Nrowsvec.push_back(Atmp[s].block[q].rows());
1559 }
1560 }
1561
1562 // Do the glue
1563 size_t Ncols = Atmp[svec[0]].block[qvec[0]].cols();
1564 for (size_t i=1; i<svec.size(); ++i) {assert(Atmp[svec[i]].block[qvec[i]].cols() == Ncols);}
1565 size_t Nrows = accumulate(Nrowsvec.begin(),Nrowsvec.end(),0);
1566
1567 MatrixType Aclump(Nrows,Ncols);
1568 Aclump.setZero();
1569 size_t stitch = 0;
1570 for (size_t i=0; i<svec.size(); ++i)
1571 {
1572 Aclump.block(stitch,0, Nrowsvec[i],Ncols) = Atmp[svec[i]].block[qvec[i]];
1573 stitch += Nrowsvec[i];
1574 }
1575
1576// Aclump *= C[loc].block[qC].adjoint();
1577
1578 #ifdef DONT_USE_BDCSVD
1579 JacobiSVD<MatrixType> Jack; // standard SVD
1580 #else
1581 BDCSVD<MatrixType> Jack; // "Divide and conquer" SVD (only available in Eigen)
1582 #endif
1583 Jack.compute(Aclump,ComputeThinU|ComputeThinV);
1584 VectorXd SV = Jack.singularValues();
1585
1586 //Here is probably the place for truncations of the Mps by taking Nret dependent on singluarValues() < eps_svd
1587 size_t Nret = Jack.singularValues().rows();
1588 // size_t Nret = (SV.array() > this->eps_svd).count();
1589 // Nret = max(Nret, this->min_Nsv);
1590 // Nret = min(Nret, this->max_Nsv);
1591 // cout << "L: Nret=" << Nret << ", full=" << Jack.singularValues().rows() << endl;
1592 // truncWeightSub(qout) = Symmetry::degeneracy(outbase[loc][qout]) * SV.tail(SV.rows()-Nret).cwiseAbs2().sum();
1593
1594 // Update AL
1595 stitch = 0;
1596 for (size_t i=0; i<svec.size(); ++i)
1597 {
1598 A[GAUGE::L][loc][svec[i]].block[qvec[i]] = Jack.matrixU().block(stitch,0, Nrowsvec[i],Nret) *
1599 Jack.matrixV().adjoint().topRows(Nret);
1600 stitch += Nrowsvec[i];
1601 }
1602 }
1603 truncWeight(loc) = truncWeightSub.sum();
1604 }
1605
1606 size_t locC = minus1modL(loc);
1607
1608 if (gauge == GAUGE::R or gauge == GAUGE::C)
1609 {
1610 ArrayXd truncWeightSub(inbase[loc].Nq()); truncWeightSub.setZero();
1611 vector<Biped<Symmetry,MatrixType> > Atmp(qloc[loc].size());
1612 for (size_t s=0; s<qloc[loc].size(); ++s)
1613 {
1614 // Atmp[s] = C[locC].adjoint().contract(A[GAUGE::C][loc][s]);
1615 Atmp[s] = C[locC].adjoint() * A[GAUGE::C][loc][s];
1616 }
1617
1618// cout << "svdDecompose AR from C at " << locC << endl;
1619
1620 for (size_t qin=0; qin<inbase[loc].Nq(); ++qin)
1621 {
1622// qarray2<Symmetry::Nq> quple = {inbase[loc][qin], inbase[loc][qin]};
1623// auto it = C[locC].dict.find(quple);
1624// assert(it != C[locC].dict.end());
1625// size_t qC = it->second;
1626
1627 // Determine how many A's to glue together
1628 vector<size_t> svec, qvec, Ncolsvec;
1629 for (size_t s=0; s<qloc[loc].size(); ++s)
1630 for (size_t q=0; q<Atmp[s].dim; ++q)
1631 {
1632 if (Atmp[s].in[q] == inbase[loc][qin])
1633 {
1634 svec.push_back(s);
1635 qvec.push_back(q);
1636 Ncolsvec.push_back(Atmp[s].block[q].cols());
1637 }
1638 }
1639
1640 // Do the glue
1641 size_t Nrows = Atmp[svec[0]].block[qvec[0]].rows();
1642 for (size_t i=1; i<svec.size(); ++i) {assert(Atmp[svec[i]].block[qvec[i]].rows() == Nrows);}
1643 size_t Ncols = accumulate(Ncolsvec.begin(), Ncolsvec.end(), 0);
1644
1645 MatrixType Aclump(Nrows,Ncols);
1646 Aclump.setZero();
1647 size_t stitch = 0;
1648 for (size_t i=0; i<svec.size(); ++i)
1649 {
1650 Aclump.block(0,stitch, Nrows,Ncolsvec[i]) = Atmp[svec[i]].block[qvec[i]]*
1651 Symmetry::coeff_leftSweep(
1652 Atmp[svec[i]].out[qvec[i]],
1653 Atmp[svec[i]].in[qvec[i]]);
1654 stitch += Ncolsvec[i];
1655 }
1656
1657// Aclump = C[locC].block[qC].adjoint() * Aclump;
1658
1659 #ifdef DONT_USE_BDCSVD
1660 JacobiSVD<MatrixType> Jack; // standard SVD
1661 #else
1662 BDCSVD<MatrixType> Jack; // "Divide and conquer" SVD (only available in Eigen)
1663 #endif
1664 Jack.compute(Aclump,ComputeThinU|ComputeThinV);
1665 VectorXd SV = Jack.singularValues();
1666
1667 //Here is probably the place for truncations of the Mps by taking Nret dependent on singluarValues() < eps_svd
1668 size_t Nret = Jack.singularValues().rows();
1669 // size_t Nret = (SV.array() > this->eps_svd).count();
1670 // Nret = max(Nret, this->min_Nsv);
1671 // Nret = min(Nret, this->max_Nsv);
1672 // cout << "R: Nret=" << Nret << ", full=" << Jack.singularValues().rows() << endl;
1673 // truncWeightSub(qin) = Symmetry::degeneracy(inbase[loc][qin]) * SV.tail(SV.rows()-Nret).cwiseAbs2().sum();
1674
1675 // Update AR
1676 stitch = 0;
1677 for (size_t i=0; i<svec.size(); ++i)
1678 {
1679 A[GAUGE::R][loc][svec[i]].block[qvec[i]] = Jack.matrixU().leftCols(Nret) *
1680 Jack.matrixV().adjoint().block(0,stitch, Nret,Ncolsvec[i])
1681 *
1682 Symmetry::coeff_leftSweep(
1683 Atmp[svec[i]].in[qvec[i]],
1684 Atmp[svec[i]].out[qvec[i]]);
1685 stitch += Ncolsvec[i];
1686 }
1687 }
1688 truncWeight(loc) = truncWeightSub.sum();
1689 }
1690}
1691
1692template<typename Symmetry, typename Scalar>
1693vector<vector<Biped<Symmetry,Matrix<complex<Scalar>,Dynamic,Dynamic>>>>
1694apply_symm (const vector<vector<Biped<Symmetry,Matrix<complex<Scalar>,Dynamic,Dynamic>>>> &A,
1695 const Mpo<Symmetry,complex<Scalar>> &R,
1696 const vector<vector<qarray<Symmetry::Nq> > > &qloc,
1697 const vector<Qbasis<Symmetry> > &qauxAl,
1698 const vector<Qbasis<Symmetry> > &qauxAr,
1699 bool TRANSPOSE=false,
1700 bool CONJUGATE=false)
1701{
1702 vector<vector<Biped<Symmetry,Matrix<complex<Scalar>,Dynamic,Dynamic>>>> Ares(A.size());
1703 for (int l=0; l<A.size(); ++l) Ares[l].resize(A[l].size());
1704
1705 auto Aprep = A;
1706 auto qloc_cpy = qloc;
1707 auto qauxAl_cpy = qauxAl;
1708 auto qauxAr_cpy = qauxAr;
1709
1710 if (CONJUGATE and not TRANSPOSE)
1711 {
1712 for (int l=0; l<A.size(); ++l)
1713 for (int s=0; s<A[l].size(); ++s)
1714 {
1715 Aprep[l][s] = A[l][s].conjugate();
1716 }
1717 }
1718 else if (TRANSPOSE)
1719 {
1720 int linv = A.size()-1;
1721 for (int l=0; l<A.size(); ++l)
1722 {
1723 for (int s=0; s<A[l].size(); ++s)
1724 {
1725// cout << "l=" << l << ", s=" << s << endl;
1726 if (CONJUGATE)
1727 {
1728 Aprep[l][s] = A[linv][s].adjoint();
1729 }
1730 else
1731 {
1732 Aprep[l][s] = A[linv][s].transpose();
1733// cout << Aprep[l][s] << endl << endl;
1734 }
1735 }
1736
1737 qloc_cpy[l] = qloc[linv];
1738 qauxAl_cpy[l] = qauxAr[linv];
1739 qauxAr_cpy[l] = qauxAl[linv];
1740
1741 --linv;
1742 }
1743 }
1744
1745 for (size_t l=0; l<A.size(); ++l)
1746 {
1747 contract_AW(Aprep[l], qloc_cpy[l], R.W_at(0), R.opBasis(0),
1748 qauxAl_cpy[l], R.inBasis(0),
1749 qauxAr_cpy[l], R.outBasis(0),
1750 Ares[l],
1751 false);
1752 }
1753
1754 return Ares;
1755}
1756
1757//template<typename Symmetry, typename Scalar>
1758//vector<vector<Biped<Symmetry,Matrix<complex<Scalar>,Dynamic,Dynamic>>>>
1759//spin_rotation (const vector<vector<Biped<Symmetry,Matrix<complex<Scalar>,Dynamic,Dynamic>>>> &A,
1760// const Mpo<Symmetry> &R,
1761// const vector<vector<qarray<Symmetry::Nq> > > &qloc,
1762// const vector<Qbasis<Symmetry> > &qauxAl,
1763// const vector<Qbasis<Symmetry> > &qauxAr)
1764//{
1765// vector<vector<Biped<Symmetry,Matrix<complex<Scalar>,Dynamic,Dynamic>>>> Ares(A.size());
1766// for (int l=0; l<A.size(); ++l) Ares[l].resize(A[l].size());
1767//
1768// for (size_t l=0; l<A.size(); ++l)
1769// {
1770// contract_AW(A[l], qloc[l], R.W_at(0), R.opBasis(0),
1771// qauxAl[l], R.inBasis(0),
1772// qauxAr[l], R.outBasis(0),
1773// Ares[l],
1774// false, {});
1775// }
1776//
1777// return Ares;
1778//}
1779
1780template<typename Symmetry, typename Scalar>
1781void Umps<Symmetry,Scalar>::
1782calc_N (DMRG::DIRECTION::OPTION DIR, size_t loc, vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &N) const
1783{
1784 N.clear();
1785 N.resize(qloc[loc].size());
1786
1787 if (DIR == DMRG::DIRECTION::LEFT)
1788 {
1789 for (size_t qin=0; qin<inbase[loc].Nq(); ++qin)
1790 {
1791 // determine how many A's to glue together
1792 vector<size_t> svec, qvec, Ncolsvec;
1793 for (size_t s=0; s<qloc[loc].size(); ++s)
1794 for (size_t q=0; q<A[GAUGE::R][loc][s].dim; ++q)
1795 {
1796 if (A[GAUGE::R][loc][s].in[q] == inbase[loc][qin])
1797 {
1798 svec.push_back(s);
1799 qvec.push_back(q);
1800 Ncolsvec.push_back(A[GAUGE::R][loc][s].block[q].cols());
1801 }
1802 }
1803
1804 if (Ncolsvec.size() > 0)
1805 {
1806 // do the glue
1807 size_t Nrows = A[GAUGE::R][loc][svec[0]].block[qvec[0]].rows();
1808 for (size_t i=1; i<svec.size(); ++i) {assert(A[GAUGE::R][loc][svec[i]].block[qvec[i]].rows() == Nrows);}
1809 size_t Ncols = accumulate(Ncolsvec.begin(), Ncolsvec.end(), 0);
1810
1811 MatrixType Aclump(Nrows,Ncols);
1812 size_t stitch = 0;
1813 for (size_t i=0; i<svec.size(); ++i)
1814 {
1815 Aclump.block(0,stitch, Nrows,Ncolsvec[i]) = A[GAUGE::R][loc][svec[i]].block[qvec[i]]*
1816 Symmetry::coeff_leftSweep(
1817 A[GAUGE::R][loc][svec[i]].out[qvec[i]],
1818 A[GAUGE::R][loc][svec[i]].in[qvec[i]]);
1819 stitch += Ncolsvec[i];
1820 }
1821
1822 HouseholderQR<MatrixType> Quirinus(Aclump.adjoint());
1823 MatrixType Qmatrix = Quirinus.householderQ().adjoint();
1824 size_t Nret = Nrows; // retained states
1825
1826 // fill N
1827 stitch = 0;
1828 for (size_t i=0; i<svec.size(); ++i)
1829 {
1830 if (Qmatrix.rows() > Nret)
1831 {
1832 size_t Nnull = Qmatrix.rows()-Nret;
1833 MatrixType Mtmp = Qmatrix.block(Nret,stitch, Nnull,Ncolsvec[i])*
1834 Symmetry::coeff_leftSweep(
1835 A[GAUGE::R][loc][svec[i]].in[qvec[i]],
1836 A[GAUGE::R][loc][svec[i]].out[qvec[i]]);
1837 N[svec[i]].try_push_back(A[GAUGE::R][loc][svec[i]].in[qvec[i]], A[GAUGE::R][loc][svec[i]].out[qvec[i]], Mtmp);
1838 }
1839 stitch += Ncolsvec[i];
1840 }
1841 }
1842 }
1843
1844 Qbasis<Symmetry> qloc_(qloc[loc]);
1845 Qbasis<Symmetry> qcomb = outbase[loc].combine(qloc_,true);
1846
1847 for (size_t qout=0; qout<outbase[loc].Nq(); ++qout)
1848 for (size_t s=0; s<qloc[loc].size(); ++s)
1849 {
1850 auto qfulls = Symmetry::reduceSilent(outbase[loc][qout], Symmetry::flip(qloc[loc][s]));
1851 for (const auto &qfull:qfulls)
1852 {
1853 qarray2<Symmetry::Nq> quple = {qfull,outbase[loc][qout]};
1854 auto it = A[GAUGE::R][loc][s].dict.find(quple);
1855 if (it == A[GAUGE::R][loc][s].dict.end())
1856 {
1857 MatrixType Mtmp(qcomb.inner_dim(qfull), outbase[loc].inner_dim(outbase[loc][qout]));
1858 Mtmp.setZero();
1859 Index down=qcomb.leftAmount(qfull,{outbase[loc][qout], Symmetry::flip(qloc[loc][s])});
1860
1861 size_t source_dim;
1862 auto it = qcomb.history.find(qfull);
1863 for (size_t i=0; i<(it->second).size(); i++)
1864 {
1865 if ((it->second)[i].source == qarray2<Nq>{outbase[loc][qout], Symmetry::flip(qloc[loc][s])})
1866 {
1867 source_dim = (it->second)[i].dim;
1868 }
1869 }
1870 Mtmp.block(down,0,source_dim,outbase[loc].inner_dim(outbase[loc][qout])).setIdentity();
1871 Mtmp.block(down,0,source_dim,outbase[loc].inner_dim(outbase[loc][qout])) *= Symmetry::coeff_leftSweep( qfull,
1872 outbase[loc][qout]);
1873 // Mtmp.setIdentity();
1874 // cout << "push_back with q=" << quple[0] << "," << quple[1] << endl;
1875 N[s].push_back(quple, Mtmp);
1876
1877 // Mtmp.setIdentity();
1878 // Mtmp *= Symmetry::coeff_sign( outbase[loc][qout], qfull, qloc[loc][s] );
1879 // N[s].push_back(quple, Mtmp);
1880 }
1881 }
1882 }
1883 }
1884 else if (DIR == DMRG::DIRECTION::RIGHT)
1885 {
1886 for (size_t qout=0; qout<outbase[loc].Nq(); ++qout)
1887 {
1888 // determine how many A's to glue together
1889 vector<size_t> svec, qvec, Nrowsvec;
1890 for (size_t s=0; s<qloc[loc].size(); ++s)
1891 for (size_t q=0; q<A[GAUGE::L][loc][s].dim; ++q)
1892 {
1893 if (A[GAUGE::L][loc][s].out[q] == outbase[loc][qout])
1894 {
1895 svec.push_back(s);
1896 qvec.push_back(q);
1897 Nrowsvec.push_back(A[GAUGE::L][loc][s].block[q].rows());
1898 }
1899 }
1900
1901 if (Nrowsvec.size() > 0)
1902 {
1903 // do the glue
1904 size_t Ncols = A[GAUGE::L][loc][svec[0]].block[qvec[0]].cols();
1905 for (size_t i=1; i<svec.size(); ++i) {assert(A[GAUGE::L][loc][svec[i]].block[qvec[i]].cols() == Ncols);}
1906 size_t Nrows = accumulate(Nrowsvec.begin(),Nrowsvec.end(),0);
1907
1908 MatrixType Aclump(Nrows,Ncols);
1909 Aclump.setZero();
1910 size_t stitch = 0;
1911 for (size_t i=0; i<svec.size(); ++i)
1912 {
1913 Aclump.block(stitch,0, Nrowsvec[i],Ncols) = A[GAUGE::L][loc][svec[i]].block[qvec[i]];
1914 stitch += Nrowsvec[i];
1915 }
1916 HouseholderQR<MatrixType> Quirinus(Aclump);
1917 MatrixType Qmatrix = Quirinus.householderQ();
1918 size_t Nret = Ncols; // retained states
1919
1920 // fill N
1921 stitch = 0;
1922 for (size_t i=0; i<svec.size(); ++i)
1923 {
1924 if (Qmatrix.cols() > Nret)
1925 {
1926 size_t Nnull = Qmatrix.cols()-Nret;
1927 MatrixType Mtmp = Qmatrix.block(stitch,Nret, Nrowsvec[i],Nnull);
1928 N[svec[i]].try_push_back(A[GAUGE::L][loc][svec[i]].in[qvec[i]], A[GAUGE::L][loc][svec[i]].out[qvec[i]], Mtmp);
1929 }
1930 stitch += Nrowsvec[i];
1931 }
1932 }
1933 }
1934
1935 Qbasis<Symmetry> qloc_(qloc[loc]);
1936 Qbasis<Symmetry> qcomb = inbase[loc].combine(qloc_);
1937
1938 for (size_t qin=0; qin<inbase[loc].Nq(); ++qin)
1939 for (size_t s=0; s<qloc[loc].size(); ++s)
1940 {
1941 auto qfulls = Symmetry::reduceSilent(inbase[loc][qin], qloc[loc][s]);
1942 for (const auto &qfull:qfulls)
1943 {
1944 qarray2<Symmetry::Nq> quple = {inbase[loc][qin], qfull};
1945 auto it = A[GAUGE::L][loc][s].dict.find(quple);
1946 if (it == A[GAUGE::L][loc][s].dict.end())
1947 {
1948 MatrixType Mtmp(inbase[loc].inner_dim(inbase[loc][qin]), qcomb.inner_dim(qfull));
1949 Mtmp.setZero();
1950 Index left=qcomb.leftAmount(qfull,{inbase[loc][qin], qloc[loc][s]});
1951
1952 size_t source_dim;
1953 auto it = qcomb.history.find(qfull);
1954 for (size_t i=0; i<(it->second).size(); i++)
1955 {
1956 if ((it->second)[i].source == qarray2<Nq>{inbase[loc][qin], qloc[loc][s]})
1957 {
1958 source_dim = (it->second)[i].dim;
1959 }
1960 }
1961 Mtmp.block(0,left,inbase[loc].inner_dim(inbase[loc][qin]),source_dim).setIdentity();
1962 // Mtmp.setIdentity();
1963 N[s].push_back(quple, Mtmp);
1964 }
1965 }
1966 }
1967 }
1968}
1969
1970template<typename Symmetry, typename Scalar>
1971template<typename OtherScalar>
1972Umps<Symmetry,OtherScalar> Umps<Symmetry,Scalar>::
1973cast() const
1974{
1975 Umps<Symmetry,OtherScalar> Vout;
1976 Vout.resize(*this);
1977
1978 for (size_t g=0; g<3; ++g)
1979 for (size_t l=0; l<this->N_sites; ++l)
1980 for (size_t s=0; s<qloc[l].size(); ++s)
1981 for (size_t q=0; q<A[g][l][s].dim; ++q)
1982 {
1983 Vout.A[g][l][s].block[q] = A[g][l][s].block[q].template cast<OtherScalar>();
1984 }
1985 for (size_t l=0; l<this->N_sites; ++l)
1986 for (size_t q=0; q<C[l].dim; ++q)
1987 {
1988 Vout.C[l].block[q] = C[l].block[q].template cast<OtherScalar>();
1989 }
1990 Vout.eps_svd = this->eps_svd;
1991 Vout.eps_truncWeight = this->eps_truncWeight;
1992 Vout.truncWeight = truncWeight;
1993
1994 return Vout;
1995}
1996
1997template<typename Symmetry, typename Scalar>
1998Umps<Symmetry,double> Umps<Symmetry,Scalar>::
1999real() const
2000{
2001 Umps<Symmetry,double> Vout;
2002 Vout.resize(*this);
2003
2004 for (size_t g=0; g<3; ++g)
2005 for (size_t l=0; l<this->N_sites; ++l)
2006 for (size_t s=0; s<qloc[l].size(); ++s)
2007 for (size_t q=0; q<A[g][l][s].dim; ++q)
2008 {
2009 Vout.A[g][l][s].block[q] = A[g][l][s].block[q].real();
2010 }
2011 for (size_t l=0; l<this->N_sites; ++l)
2012 for (size_t q=0; q<C[l].dim; ++q)
2013 {
2014 Vout.C[l].block[q] = C[l].block[q].real();
2015 }
2016 Vout.eps_svd = this->eps_svd;
2017 Vout.eps_truncWeight = this->eps_truncWeight;
2018 Vout.truncWeight = truncWeight;
2019
2020 return Vout;
2021}
2022
2023template<typename Symmetry, typename Scalar>
2024void Umps<Symmetry,Scalar>::
2025adjustQN (const size_t number_cells)
2026{
2027 //transform quantum number in all Bipeds
2028 for (size_t g=0; g<3; ++g)
2029 for (size_t l=0; l<N_sites; ++l)
2030 for (size_t s=0; s<qloc[l].size(); ++s)
2031 {
2032 A[g][l][s] = A[g][l][s].adjustQN(number_cells);
2033 }
2034
2035 //transform physical quantum numbers
2036 for (size_t l=0; l<N_sites; ++l)
2037 for (size_t s=0; s<qloc[l].size(); ++s)
2038 {
2039 qloc[l][s] = ::adjustQN<Symmetry>(qloc[l][s],number_cells);
2040 }
2041
2042 update_inbase();
2043 update_outbase();
2044};
2045
2046#ifdef USE_HDF5_STORAGE
2047template<typename Symmetry, typename Scalar>
2048void Umps<Symmetry,Scalar>::
2049save (string filename, string info, double energy, double err_var, double err_state)
2050{
2051 filename+=".h5";
2052 HDF5Interface target(filename, WRITE);
2053 target.create_group("As");
2054 target.create_group("Cs");
2055 target.create_group("qloc");
2056 target.create_group("Qtot");
2057
2058 string add_infoLabel = "add_info";
2059
2060 //save scalar values
2061 if (!isnan(energy))
2062 {
2063 target.save_scalar(energy,"energy");
2064 }
2065 if (!isnan(err_var))
2066 {
2067 target.save_scalar(err_var,"err_var");
2068 }
2069 if (!isnan(err_state))
2070 {
2071 target.save_scalar(err_state,"err_state");
2072 }
2073 target.save_scalar(this->N_sites,"L");
2074 for (size_t q=0; q<Nq; q++)
2075 {
2076 stringstream ss; ss << "q=" << q;
2077 target.save_scalar(this->Qtot[q],ss.str(),"Qtot");
2078 }
2079 target.save_scalar(this->calc_Dmax(),"Dmax");
2080 target.save_scalar(this->calc_Nqmax(),"Nqmax");
2081 target.save_scalar(this->min_Nsv,"min_Nsv");
2082 target.save_scalar(this->max_Nsv,"max_Nsv");
2083 target.save_scalar(this->eps_svd,"eps_svd");
2084 target.save_scalar(this->eps_truncWeight,"eps_truncWeight");
2085 target.save_char(info,add_infoLabel.c_str());
2086
2087 //save qloc
2088 for (size_t l=0; l<this->N_sites; ++l)
2089 {
2090 stringstream ss; ss << "l=" << l;
2091 target.save_scalar(qloc[l].size(),ss.str(),"qloc");
2092 for (size_t s=0; s<qloc[l].size(); ++s)
2093 for (size_t q=0; q<Nq; q++)
2094 {
2095 stringstream tt; tt << "l=" << l << ",s=" << s << ",q=" << q;
2096 target.save_scalar((qloc[l][s])[q],tt.str(),"qloc");
2097 }
2098 }
2099
2100 //save the A-matrices
2101 string label;
2102 for (size_t g=0; g<3; ++g)
2103 for (size_t l=0; l<this->N_sites; ++l)
2104 for (size_t s=0; s<qloc[l].size(); ++s)
2105 {
2106 stringstream tt; tt << "g=" << g << ",l=" << l << ",s=" << s;
2107 target.save_scalar(A[g][l][s].dim,tt.str());
2108 for (size_t q=0; q<A[g][l][s].dim; ++q)
2109 {
2110 for (size_t p=0; p<Nq; p++)
2111 {
2112 stringstream in; in << "in,g=" << g << ",l=" << l << ",s=" << s << ",q=" << q << ",p=" << p;
2113 stringstream out; out << "out,g=" << g << ",l=" << l << ",s=" << s << ",q=" << q << ",p=" << p;
2114 target.save_scalar((A[g][l][s].in[q])[p],in.str(),"As");
2115 target.save_scalar((A[g][l][s].out[q])[p],out.str(),"As");
2116 }
2117 stringstream ss;
2118 ss << g << "_" << l << "_" << s << "_" << "(" << A[g][l][s].in[q] << "," << A[g][l][s].out[q] << ")";
2119 label = ss.str();
2120 if constexpr (std::is_same<Scalar,complex<double>>::value)
2121 {
2122 MatrixXd Re = A[g][l][s].block[q].real();
2123 MatrixXd Im = A[g][l][s].block[q].imag();
2124 target.save_matrix(Re, label+"Re", "As");
2125 target.save_matrix(Im, label+"Im", "As");
2126 }
2127 else
2128 {
2129 target.save_matrix(A[g][l][s].block[q], label, "As");
2130 }
2131 }
2132 }
2133
2134 //save the C-matrices
2135 for (size_t l=0; l<this->N_sites; ++l)
2136 {
2137 stringstream tt; tt << "l=" << l;
2138 target.save_scalar(C[l].dim,tt.str());
2139 for (size_t q=0; q<C[l].dim; ++q)
2140 {
2141 for (size_t p=0; p<Nq; p++)
2142 {
2143 stringstream in; in << "l=" << l << ",q=" << q << ",p=" << p;
2144 target.save_scalar((C[l].in[q])[p],in.str(),"Cs");
2145 }
2146 stringstream ss;
2147 ss << l << "_" << "(" << C[l].in[q] << ")";
2148 label = ss.str();
2149 if constexpr (std::is_same<Scalar,complex<double>>::value)
2150 {
2151 MatrixXd Re = C[l].block[q].real();
2152 MatrixXd Im = C[l].block[q].imag();
2153 target.save_matrix(Re, label+"Re", "Cs");
2154 target.save_matrix(Im, label+"Im", "Cs");
2155 }
2156 else
2157 {
2158 target.save_matrix(C[l].block[q],label,"Cs");
2159 }
2160 }
2161 }
2162 target.close();
2163}
2164
2165template<typename Symmetry, typename Scalar>
2166void Umps<Symmetry,Scalar>::
2167load (string filename, double &energy, double &err_var, double &err_state)
2168{
2169 filename+=".h5";
2170 HDF5Interface source(filename, READ);
2171
2172 //load the scalars
2173 if (source.CHECK("energy"))
2174 {
2175 source.load_scalar(energy,"energy");
2176 }
2177 if (source.CHECK("err_var"))
2178 {
2179 source.load_scalar(energy,"err_var");
2180 }
2181 if (source.CHECK("err_state"))
2182 {
2183 source.load_scalar(energy,"err_state");
2184 }
2185 source.load_scalar(this->N_sites,"L");
2186 for (size_t q=0; q<Nq; q++)
2187 {
2188 stringstream ss; ss << "q=" << q;
2189 source.load_scalar(this->Qtot[q],ss.str(),"Qtot");
2190 }
2191 source.load_scalar(this->eps_svd,"eps_svd");
2192 // To ensure older files can be loaded, make check here
2193 // HAS_GROUP is the same for groups and single objects
2194 if (source.HAS_GROUP("eps_truncWeight")) source.load_scalar(this->eps_truncWeight,"eps_truncWeight");
2195 source.load_scalar(this->min_Nsv,"min_Nsv");
2196 source.load_scalar(this->max_Nsv,"max_Nsv");
2197
2198 //load qloc
2199 qloc.resize(this->N_sites);
2200 for (size_t l=0; l<this->N_sites; ++l)
2201 {
2202 stringstream ss; ss << "l=" << l;
2203 size_t qloc_size;
2204 source.load_scalar(qloc_size,ss.str(),"qloc");
2205 qloc[l].resize(qloc_size);
2206 for (size_t s=0; s<qloc[l].size(); ++s)
2207 for (size_t q=0; q<Nq; q++)
2208 {
2209 stringstream tt; tt << "l=" << l << ",s=" << s << ",q=" << q;
2210 int Q;
2211 source.load_scalar(Q,tt.str(),"qloc");
2212 (qloc[l][s])[q] = Q;
2213 }
2214 }
2215 this->resize_arrays();
2216
2217 //load the A-matrices
2218 string label;
2219 for (size_t g=0; g<3; ++g)
2220 for (size_t l=0; l<this->N_sites; ++l)
2221 for (size_t s=0; s<qloc[l].size(); ++s)
2222 {
2223 size_t Asize;
2224 stringstream tt; tt << "g=" << g << ",l=" << l << ",s=" << s;
2225 source.load_scalar(Asize,tt.str());
2226 for (size_t q=0; q<Asize; ++q)
2227 {
2228 qarray<Nq> qin,qout;
2229 for (size_t p=0; p<Nq; p++)
2230 {
2231 stringstream in; in << "in,g=" << g << ",l=" << l << ",s=" << s << ",q=" << q << ",p=" << p;
2232 stringstream out; out << "out,g=" << g << ",l=" << l << ",s=" << s << ",q=" << q << ",p=" << p;
2233 source.load_scalar(qin[p],in.str(),"As");
2234 source.load_scalar(qout[p],out.str(),"As");
2235 }
2236 stringstream ss;
2237 ss << g << "_" << l << "_" << s << "_" << "(" << qin << "," << qout << ")";
2238 label = ss.str();
2239 MatrixType mat;
2240 if constexpr (std::is_same<Scalar,complex<double>>::value)
2241 {
2242 MatrixXd Re, Im;
2243 source.load_matrix(Re, label+"Re", "As");
2244 source.load_matrix(Im, label+"Im", "As");
2245 mat = Re+1.i*Im;
2246 }
2247 else
2248 {
2249 source.load_matrix(mat, label, "As");
2250 }
2251 A[g][l][s].push_back(qin,qout,mat);
2252 }
2253 }
2254
2255 //load the C-matrices
2256 label.clear();
2257 for (size_t l=0; l<this->N_sites; ++l)
2258 {
2259 size_t Asize;
2260 stringstream tt; tt << "l=" << l;
2261 source.load_scalar(Asize,tt.str());
2262 for (size_t q=0; q<Asize; ++q)
2263 {
2264 qarray<Nq> qVal;
2265 for (size_t p=0; p<Nq; p++)
2266 {
2267 stringstream qq; qq << "l=" << l << ",q=" << q << ",p=" << p;
2268 source.load_scalar(qVal[p],qq.str(),"Cs");
2269 }
2270 stringstream ss;
2271 ss << l << "_" << "(" << qVal << ")";
2272 label = ss.str();
2273 MatrixType mat;
2274 if constexpr (std::is_same<Scalar,complex<double>>::value)
2275 {
2276 MatrixXd Re, Im;
2277 source.load_matrix(Re, label+"Re", "Cs");
2278 source.load_matrix(Im, label+"Im", "Cs");
2279 mat = Re+1.i*Im;
2280 }
2281 else
2282 {
2283 source.load_matrix(mat, label, "Cs");
2284 }
2285 C[l].push_back(qVal,qVal,mat);
2286 }
2287 }
2288 source.close();
2289 update_inbase();
2290 update_outbase();
2291}
2292#endif //USE_HDF5_STORAGE
2293
2294template<typename Symmetry, typename Scalar>
2295void Umps<Symmetry,Scalar>::
2296sort_A(size_t loc, GAUGE::OPTION g, bool SORT_ALL_GAUGES)
2297{
2298 if (SORT_ALL_GAUGES)
2299 {
2300 for (size_t gP=0; gP<3; ++gP)
2301 for (size_t s=0; s<locBasis(loc).size(); ++s)
2302 {
2303 A[gP][loc][s] = A[gP][loc][s].sorted();
2304 }
2305 }
2306 else
2307 {
2308 for (size_t s=0; s<locBasis(loc).size(); ++s)
2309 {
2310 A[g][loc][s] = A[g][loc][s].sorted();
2311 }
2312 }
2313}
2314
2315template<typename Symmetry, typename Scalar>
2316void Umps<Symmetry,Scalar>::
2317updateC(size_t loc)
2318{
2319 for (size_t q=0; q<outBasis(loc).Nq(); ++q)
2320 {
2321 qarray2<Symmetry::Nq> quple = {outBasis(loc)[q], outBasis(loc)[q]};
2322 auto qC = C[loc].dict.find(quple);
2323 size_t r = outBasis(loc).inner_dim(outBasis(loc)[q]);
2324 size_t c = r;
2325 if (qC != C[loc].dict.end())
2326 {
2327 int dr = r-C[loc].block[qC->second].rows();
2328 int dc = c-C[loc].block[qC->second].cols();
2329
2330 C[loc].block[qC->second].conservativeResize(r,c);
2331
2332 C[loc].block[qC->second].bottomRows(dr).setZero();
2333 C[loc].block[qC->second].rightCols(dc).setZero();
2334 }
2335 else
2336 {
2337 MatrixType Mtmp(r,c);
2338 Mtmp.setZero();
2339 C[loc].push_back(quple, Mtmp);
2340 }
2341 }
2342}
2343
2344template<typename Symmetry, typename Scalar>
2345void Umps<Symmetry,Scalar>::
2346updateAC(size_t loc, GAUGE::OPTION g)
2347{
2348 assert(g != GAUGE::C and "Oouuhh.. you tried to update AC with itself, but we have no bootstrap implemented ;). Use AL or AR.");
2349 for (size_t s=0; s<locBasis(loc).size(); ++s)
2350 for (size_t q=0; q<A[g][loc][s].size(); ++q)
2351 {
2352 qarray2<Symmetry::Nq> quple = {A[g][loc][s].in[q], A[g][loc][s].out[q]};
2353 auto it = A[GAUGE::C][loc][s].dict.find(quple);
2354 if (it != A[GAUGE::C][loc][s].dict.end())
2355 {
2356 int dr = A[g][loc][s].block[q].rows() - A[GAUGE::C][loc][s].block[it->second].rows();
2357 int dc = A[g][loc][s].block[q].cols() - A[GAUGE::C][loc][s].block[it->second].cols();
2358 assert(dr >= 0 and dc >= 0 and "Something went wrong in expand_basis during the VUMPS Algorithm.");
2359 MatrixType Mtmp(dr,A[GAUGE::C][loc][s].block[it->second].cols()); Mtmp.setZero();
2360 addBottom(Mtmp, A[GAUGE::C][loc][s].block[it->second]);
2361 Mtmp.resize(A[GAUGE::C][loc][s].block[it->second].rows(),dc); Mtmp.setZero();
2362 addRight(Mtmp, A[GAUGE::C][loc][s].block[it->second]);
2363 }
2364 else
2365 {
2366 MatrixType Mtmp(A[g][loc][s].block[q].rows(), A[g][loc][s].block[q].cols());
2367 Mtmp.setZero();
2368 A[GAUGE::C][loc][s].push_back(quple,Mtmp);
2369 }
2370 }
2371}
2372
2373template<typename Symmetry, typename Scalar>
2374void Umps<Symmetry,Scalar>::
2375enrich (size_t loc, GAUGE::OPTION g, const vector<Biped<Symmetry,MatrixType> > &P)
2376{
2377 auto Pcopy = P;
2378 size_t loc1,loc2;
2379 if (g == GAUGE::L) {loc1 = loc; loc2 = (loc+1)%N_sites;}
2380 else if (g == GAUGE::R) {loc1 = (loc+1)%N_sites; loc2 = loc;}
2381
2382 Qbasis<Symmetry> base_P; if (g == GAUGE::L) {base_P.pullData(P,1);} else if (g == GAUGE::R) {base_P.pullData(P,0);}
2383
2384 Qbasis<Symmetry> expanded_base; if (g == GAUGE::L) {expanded_base=outBasis(loc1).add(base_P);} else if (g == GAUGE::R) {expanded_base=inBasis(loc1).add(base_P);}
2385 // cout << "new basis states for the expansion" << endl << expanded_base << endl;
2386
2387 if (g == GAUGE::L and loc1 != loc2)
2388 {
2389 for (const auto & [qin, num_in, plain_in] : inBasis(loc1))
2390 for (size_t s=0; s<locBasis(loc1).size(); ++s)
2391 {
2392 bool QIN_IS_IN_P=false;
2393 for (size_t qP=0; qP<Pcopy[s].size(); ++qP)
2394 {
2395 if (Pcopy[s].in[qP] == qin)
2396 {
2397 if (Pcopy[s].block[qP].rows() != plain_in.size()) //A[g][loc1][s].block[qA->second].rows()
2398 {
2399 addBottom(Eigen::Matrix<Scalar,-1,-1>::Zero(plain_in.size() - Pcopy[s].block[qP].rows(),Pcopy[s].block[qP].cols()), Pcopy[s].block[qP]);
2400 }
2401 QIN_IS_IN_P=true;
2402 }
2403 }
2404 if (QIN_IS_IN_P == false)
2405 {
2406 auto qouts = Symmetry::reduceSilent(qin, locBasis(loc1)[s]);
2407 for (const auto &qout : qouts)
2408 {
2409 if (!base_P.find(qout)) {continue;}
2410 Pcopy[s].push_back(qin,qout,Eigen::Matrix<Scalar,-1,-1>::Zero(plain_in.size(),base_P.inner_dim(qout)));
2411 }
2412 }
2413 }
2414 }
2415 else if (g == GAUGE::R and loc1 != loc2)
2416 {
2417 for (const auto & [qout, num_out, plain_out] : outBasis(loc1))
2418 for (size_t s=0; s<locBasis(loc1).size(); ++s)
2419 {
2420 bool QOUT_IS_IN_P=false;
2421 for (size_t qP=0; qP<Pcopy[s].size(); ++qP)
2422 {
2423 if (Pcopy[s].out[qP] == qout)
2424 {
2425 if (Pcopy[s].block[qP].cols() != plain_out.size())
2426 {
2427 addRight(Eigen::Matrix<Scalar,-1,-1>::Zero(Pcopy[s].block[qP].rows(), plain_out.size() - Pcopy[s].block[qP].cols()), Pcopy[s].block[qP]);
2428 }
2429 QOUT_IS_IN_P=true;
2430 }
2431 }
2432 if (QOUT_IS_IN_P == false)
2433 {
2434 auto qins = Symmetry::reduceSilent(qout, Symmetry::flip(locBasis(loc1)[s]));
2435 for (const auto &qin : qins)
2436 {
2437 if (!base_P.find(qin)) {continue;}
2438 Pcopy[s].push_back(qin,qout,Eigen::Matrix<Scalar,-1,-1>::Zero(base_P.inner_dim(qin),plain_out.size()));
2439 }
2440 }
2441 }
2442 }
2443
2444 for (size_t s=0; s<locBasis(loc1).size(); ++s)
2445 for (size_t qP=0; qP<Pcopy[s].size(); ++qP)
2446 {
2447 qarray2<Symmetry::Nq> quple = {Pcopy[s].in[qP], Pcopy[s].out[qP]};
2448 auto qA = A[g][loc1][s].dict.find(quple);
2449
2450 if (qA != A[g][loc1][s].dict.end())
2451 {
2452 if (g == GAUGE::L)
2453 {
2454 addRight(Pcopy[s].block[qP], A[g][loc1][s].block[qA->second]);
2455 }
2456 else if (g == GAUGE::R)
2457 {
2458 addBottom(Pcopy[s].block[qP], A[g][loc1][s].block[qA->second]);
2459 }
2460 }
2461 else
2462 {
2463 if (g == GAUGE::L)
2464 {
2465 if (Pcopy[s].block[qP].rows() != inBasis(loc1).inner_dim(quple[0]))
2466 {
2467 addBottom(Eigen::Matrix<Scalar,-1,-1>::Zero(inBasis(loc1).inner_dim(quple[0])-Pcopy[s].block[qP].rows(),Pcopy[s].block[qP].cols()),Pcopy[s].block[qP]);
2468 }
2469 }
2470 if (g == GAUGE::R)
2471 {
2472 if (Pcopy[s].block[qP].cols() != outBasis(loc1).inner_dim(quple[1]))
2473 {
2474 addRight(Eigen::Matrix<Scalar,-1,-1>::Zero(Pcopy[s].block[qP].rows(),outBasis(loc1).inner_dim(quple[1])-Pcopy[s].block[qP].cols()),Pcopy[s].block[qP]);
2475 }
2476 }
2477 A[g][loc1][s].push_back(quple, Pcopy[s].block[qP]);
2478 }
2479 }
2480
2481 //resize blocks which was not present in P with zeros.
2482 for (size_t s=0; s<A[g][loc1].size(); s++)
2483 for (size_t qA=0; qA<A[g][loc1][s].size(); qA++)
2484 {
2485 if (g == GAUGE::L)
2486 {
2487 if (A[g][loc1][s].block[qA].cols() != expanded_base.inner_dim(A[g][loc1][s].out[qA]))
2488 {
2489 addRight(Eigen::Matrix<Scalar,-1,-1>::Zero(A[g][loc1][s].block[qA].rows(), expanded_base.inner_dim(A[g][loc1][s].out[qA]) - A[g][loc1][s].block[qA].cols()), A[g][loc1][s].block[qA]);
2490 }
2491 }
2492 else if (g == GAUGE::R)
2493 {
2494 if (A[g][loc1][s].block[qA].rows() != expanded_base.inner_dim(A[g][loc1][s].in[qA]))
2495 {
2496 addBottom(Eigen::Matrix<Scalar,-1,-1>::Zero(expanded_base.inner_dim(A[g][loc1][s].in[qA]) - A[g][loc1][s].block[qA].rows(), A[g][loc1][s].block[qA].cols()), A[g][loc1][s].block[qA]);
2497 }
2498 }
2499 }
2500
2501 // update the inleg from AL (outleg from AR) at site loc2 with zeros
2502 update_inbase(loc1,g);
2503 update_outbase(loc1,g);
2504
2505 for (const auto &[qval,qdim,plain]:base_P)
2506 for (size_t s=0; s<locBasis(loc2).size(); ++s)
2507 {
2508 std::vector<qarray<Symmetry::Nq> > qins_outs;
2509 if (g == GAUGE::L) {qins_outs = Symmetry::reduceSilent(qval, locBasis(loc2)[s]);}
2510 else if (g == GAUGE::R) {qins_outs = Symmetry::reduceSilent(qval, Symmetry::flip(locBasis(loc2)[s]));}
2511 for (const auto &qin_out:qins_outs)
2512 {
2513 qarray2<Symmetry::Nq> quple;
2514 if (g == GAUGE::L)
2515 {
2516 if (outBasis(loc2).find(qin_out) == false) {continue;}
2517 quple = {qval, qin_out};
2518 }
2519 else if (g == GAUGE::R)
2520 {
2521 if (inBasis(loc2).find(qin_out) == false) {continue;}
2522 quple = {qin_out, qval};
2523 }
2524 auto it = A[g][loc2][s].dict.find(quple);
2525 if (it != A[g][loc2][s].dict.end())
2526 {
2527 if (g == GAUGE::L)
2528 {
2529 if (A[g][loc2][s].block[it->second].rows() != expanded_base.inner_dim(quple[0]))
2530 {
2531 addBottom(Eigen::Matrix<Scalar,-1,-1>::Zero(base_P.inner_dim(qval),A[g][loc2][s].block[it->second].cols()), A[g][loc2][s].block[it->second]);
2532 }
2533 }
2534 else if (g == GAUGE::R)
2535 {
2536 if (A[g][loc2][s].block[it->second].cols() != expanded_base.inner_dim(quple[1]))
2537 {
2538 addRight(Eigen::Matrix<Scalar,-1,-1>::Zero(A[g][loc2][s].block[it->second].rows(),base_P.inner_dim(qval)), A[g][loc2][s].block[it->second]);
2539 }
2540 }
2541 }
2542 else
2543 {
2544 MatrixType Mtmp;
2545 if (g == GAUGE::L) {Mtmp.resize(base_P.inner_dim(qval), outBasis(loc2).inner_dim(qin_out));}
2546 else if (g == GAUGE::R) {Mtmp.resize(inBasis(loc2).inner_dim(qin_out), base_P.inner_dim(qval));}
2547 Mtmp.setZero();
2548 A[g][loc2][s].push_back(quple, Mtmp);
2549 }
2550 }
2551 }
2552 //resize blocks which was not present in P with zeros.
2553 // for (size_t s=0; s<A[g][loc2].size(); s++)
2554 // for (size_t qA=0; qA<A[g][loc2][s].size(); qA++)
2555 // {
2556 // if (g == GAUGE::R)
2557 // {
2558 // if (A[g][loc2][s].block[qA].cols() != expanded_base.inner_dim(A[g][loc2][s].out[qA]))
2559 // {
2560 // addRight(Eigen::Matrix<Scalar,-1,-1>::Zero(A[g][loc2][s].block[qA].rows(), expanded_base.inner_dim(A[g][loc2][s].out[qA]) - A[g][loc2][s].block[qA].cols()), A[g][loc2][s].block[qA]);
2561 // }
2562 // }
2563 // else if (g == GAUGE::L)
2564 // {
2565 // if (A[g][loc2][s].block[qA].rows() != expanded_base.inner_dim(A[g][loc2][s].in[qA]))
2566 // {
2567 // addBottom(Eigen::Matrix<Scalar,-1,-1>::Zero(expanded_base.inner_dim(A[g][loc2][s].in[qA]) - A[g][loc2][s].block[qA].rows(), A[g][loc2][s].block[qA].cols()), A[g][loc2][s].block[qA]);
2568 // }
2569 // }
2570 // }
2571}
2572
2573template<typename Symmetry, typename Scalar>
2574void Umps<Symmetry,Scalar>::
2575orthogonalize_right (GAUGE::OPTION g, vector<Biped<Symmetry,MatrixType> > &G_R)
2576{
2577 vector<vector<Biped<Symmetry,MatrixType> > > A_ortho(N_sites);
2578 for (size_t l=0; l<N_sites; l++)
2579 {
2580 A_ortho[l].resize(qloc[l].size());
2581 }
2582
2583 vector<Biped<Symmetry,MatrixType> > Rprev(N_sites);
2584 Biped<Symmetry,MatrixType> Xnext;
2585 Rprev[N_sites-1].setRandom(outbase[N_sites-1],outbase[N_sites-1]);
2586 Rprev[0] = 1./sqrt(Rprev[0].squaredNorm().sum()) * Rprev[0];
2587
2588 vector<vector<Biped<Symmetry,MatrixType> > > AxR(N_sites);
2589 for (size_t loc=0; loc<N_sites; loc++)
2590 {
2591 AxR[loc].resize(qloc[loc].size());
2592 }
2593 double tol = 1.e-10, err=1.;
2594 size_t i=0, imax = 10001;
2595 while (err > tol and i < imax)
2596 {
2597 for (size_t loc=N_sites-1; loc!=-1; --loc)
2598 {
2599 for (size_t s=0; s<qloc[loc].size(); s++)
2600 {
2601 AxR[loc][s] = A[g][loc][s] * Rprev[loc];
2602 }
2603 Blocker<Symmetry,Scalar> Jim(AxR[loc], qloc[loc], inbase[loc], outbase[loc]);
2604 auto A_blocked = Jim.Aclump(DMRG::DIRECTION::RIGHT);
2605 Biped<Symmetry,MatrixType> Q,R;
2606 for (size_t q=0; q<A_blocked.dim; ++q)
2607 {
2608 HouseholderQR<MatrixType> Quirinus;
2609 Quirinus.compute(A_blocked.block[q].adjoint());
2610
2611 MatrixType Qmat = Quirinus.householderQ() * MatrixType::Identity(A_blocked.block[q].cols(),A_blocked.block[q].rows());
2612 MatrixType Rmat = MatrixType::Identity(A_blocked.block[q].rows(),A_blocked.block[q].cols()) * Quirinus.matrixQR().template triangularView<Upper>();
2613 //make the QR decomposition unique by enforcing the diagonal of R to be positive.
2614 DiagonalMatrix<Scalar,Dynamic> Sign = Rmat.diagonal().cwiseSign().matrix().asDiagonal();
2615
2616 Rmat = Sign*Rmat;
2617 Qmat = Qmat*Sign;
2618
2619 Q.push_back(A_blocked.in[q], A_blocked.out[q], (Qmat.adjoint()));
2620 R.push_back(A_blocked.in[q], A_blocked.out[q], Rmat.adjoint());
2621 }
2622 R = 1./R.operatorNorm(false) * R;
2623 if (loc>0) {Rprev[loc-1] = R;}
2624 else {Xnext = R;}
2625 A_ortho[loc] = Jim.reblock(Q, DMRG::DIRECTION::RIGHT);
2626 }
2627 err = (Xnext - Rprev[N_sites-1]).norm();
2628 Rprev[N_sites-1] = Xnext;
2629 i++;
2630 // cout << "iteration number=" << i << ", err=" << err << endl << endl;
2631 }
2632 G_R = Rprev;
2633 lout << "Orhtogonalize right: iteration number=" << i << ", err=" << err << endl << endl;
2634
2635 A[g] = A_ortho;
2636}
2637
2638template<typename Symmetry, typename Scalar>
2639void Umps<Symmetry,Scalar>::
2640orthogonalize_left(GAUGE::OPTION g, vector<Biped<Symmetry,MatrixType> > &G_L)
2641{
2642 vector<vector<Biped<Symmetry,MatrixType> > > A_ortho(N_sites);
2643 for (size_t l=0; l<N_sites; l++)
2644 {
2645 A_ortho[l].resize(qloc[l].size());
2646 }
2647
2648 vector<Biped<Symmetry,MatrixType> > Lprev(N_sites);
2649 Biped<Symmetry,MatrixType> Lnext;
2650 Lprev[0].setRandom(inbase[0],inbase[0]);
2651 Lprev[0] = 1./sqrt(Lprev[0].squaredNorm().sum()) * Lprev[0];
2652 vector<vector<Biped<Symmetry,MatrixType> > > LxA(N_sites);
2653 for (size_t loc=0; loc<N_sites; loc++)
2654 {
2655 LxA[loc].resize(qloc[loc].size());
2656 }
2657 double tol = 1.e-10, err=1.;
2658 size_t i=0, imax = 10001;
2659 while (err > tol and i < imax)
2660 {
2661 for (size_t loc=0; loc<N_sites; loc++)
2662 {
2663 for (size_t s=0; s<qloc[loc].size(); s++)
2664 {
2665 LxA[loc][s] = Lprev[loc] * A[g][loc][s];
2666 }
2667 Blocker<Symmetry,Scalar> Jim(LxA[loc], qloc[loc], inbase[loc], outbase[loc]);
2668 auto A_blocked = Jim.Aclump(DMRG::DIRECTION::LEFT);
2669 Biped<Symmetry,MatrixType> Q,R;
2670 for (size_t q=0; q<A_blocked.dim; ++q)
2671 {
2672 HouseholderQR<MatrixType> Quirinus;
2673 Quirinus.compute(A_blocked.block[q]);
2674
2675 MatrixType Qmat = Quirinus.householderQ() * MatrixType::Identity(A_blocked.block[q].rows(),A_blocked.block[q].cols());
2676 MatrixType Rmat = MatrixType::Identity(A_blocked.block[q].cols(),A_blocked.block[q].rows()) * Quirinus.matrixQR().template triangularView<Upper>();
2677 //make the QR decomposition unique by enforcing the diagonal of R to be positive.
2678 DiagonalMatrix<Scalar,Dynamic> Sign = Rmat.diagonal().cwiseSign().matrix().asDiagonal();
2679
2680 Rmat = Sign*Rmat;
2681 Qmat = Qmat*Sign;
2682
2683 Q.push_back(A_blocked.in[q], A_blocked.out[q], Qmat);
2684 R.push_back(A_blocked.in[q], A_blocked.out[q], Rmat);
2685 }
2686 R = 1./R.operatorNorm(false) * R;
2687 if (loc<N_sites-1) {Lprev[loc+1] = R;}
2688 else {Lnext = R;}
2689 A_ortho[loc] = Jim.reblock(Q, DMRG::DIRECTION::LEFT);
2690 }
2691 err = (Lnext - Lprev[0]).norm();
2692 Lprev[0] = Lnext;
2693 i++;
2694 // cout << "iteration number=" << i << ", err=" << err << endl << endl;
2695 }
2696 G_L = Lprev;
2697 A[g] = A_ortho;
2698 lout << "Orthogonalize left: iteration number=" << i << ", err=" << err << endl << endl;
2699}
2700
2701template<typename Symmetry, typename Scalar>
2702void Umps<Symmetry,Scalar>::
2703truncate (bool SET_AC_RANDOM)
2704{
2705 //isometries from the truncated SVD from the center-matrix C
2706
2707 vector<Biped<Symmetry,MatrixType> > U(N_sites);
2708 vector<Biped<Symmetry,MatrixType> > Vdag(N_sites);
2709
2710 //decompose C by SVD and write isometries to U and Vdag and the singular (Schmidt) values into C.
2711 for (size_t l=0; l<N_sites; ++l)
2712 {
2713 // cout << "**********************l=" << l << "*************************" << endl;
2714 // cout << C[l].print(false) << endl;
2715 double trunc=0.;
2716 auto [trunc_U,Sigma,trunc_Vdag] = C[l].truncateSVD(min_Nsv,max_Nsv,this->eps_truncWeight,trunc,true); //true: PRESERVE_MULTIPLETS
2717 U[l] = trunc_U;
2718 Vdag[l] = trunc_Vdag;
2719 C[l] = Sigma;
2720 }
2721
2722 //update AL and AR
2723 for (size_t l=0; l<N_sites; ++l)
2724 {
2725 for (size_t s=0; s<qloc[l].size(); ++s)
2726 {
2727 A[GAUGE::L][l][s] = U[minus1modL(l)].adjoint() * A[GAUGE::L][l][s] * U[l];
2728 A[GAUGE::R][l][s] = Vdag[minus1modL(l)] * A[GAUGE::R][l][s] * Vdag[l].adjoint();
2729 }
2730 }
2731 update_outbase(GAUGE::L);
2732 update_inbase(GAUGE::L);
2733
2734 //Orthogonalize AL and AR again and safe gauge transformation into L and R. (AL -> L*AL*Linv, AR -> Rinv*AR*R)
2735 //L and R need to be multiplied into the center matrix afterwards
2736 vector<Biped<Symmetry,MatrixType> > L(N_sites),R(N_sites);
2737 orthogonalize_right(GAUGE::R,R);
2738 orthogonalize_left(GAUGE::L,L);
2739 for (size_t l=0; l<N_sites; ++l)
2740 {
2741 C[l] = L[(l+1)%N_sites] * C[l] * R[l];
2742 C[l] = C[l].sorted();
2743 }
2744
2745 //normalize the state to get rid off small norm changes due to the truncation
2746 normalize_C();
2747
2748 //Update AC so that it has the truncated sizes and sort the A tensors.
2749 for (size_t l=0; l<N_sites; l++)
2750 {
2751 for (size_t s=0; s<qloc[l].size(); ++s)
2752 {
2753 A[GAUGE::C][l][s] = C[minus1modL(l)] * A[GAUGE::R][l][s];
2754 if (SET_AC_RANDOM) { A[GAUGE::C][l][s].setRandom(); }
2755 }
2756 sort_A(l,GAUGE::L,true); //true means sort all GAUGES. First parameter has no consequences
2757 }
2758 calc_entropy();
2759 // cout << test_ortho() << endl;
2760}
2761
2762template<typename Symmetry, typename Scalar>
2763vector<std::pair<complex<double>,Biped<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > > > Umps<Symmetry,Scalar>::
2764calc_dominant (GAUGE::OPTION g, DMRG::DIRECTION::OPTION DIR, int N, double tol, int dimK, qarray<Symmetry::Nq> Qtot, string label) const
2765{
2766 vector<std::pair<complex<double>,Biped<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > > > res(N);
2767
2768 Umps<Symmetry,complex<double> > Compl = this->template cast<complex<double> > ();
2769 complex<double> lambda1;
2770
2771 TransferMatrix<Symmetry,complex<double> > T;
2772 if (DIR == DMRG::DIRECTION::LEFT)
2773 {
2774 T = TransferMatrix<Symmetry,complex<double> >
2775 (VMPS::DIRECTION::RIGHT, Compl.A[g], Compl.A[g], locBasis(), false, Qtot);
2776 }
2777 else
2778 {
2779 T = TransferMatrix<Symmetry,complex<double> >
2780 (VMPS::DIRECTION::LEFT, Compl.A[g], Compl.A[g], locBasis(), false, Qtot);
2781 }
2782
2783 Biped<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > RandBiped;
2784 if (DIR == DMRG::DIRECTION::LEFT)
2785 {
2786 RandBiped.setRandom(inBasis(0), inBasis(0));
2787 }
2788 else
2789 {
2790 RandBiped.setRandom(outBasis(N_sites-1), outBasis(N_sites-1));
2791 }
2792 RandBiped = 1./RandBiped.norm() * RandBiped;
2793 TransferVector<Symmetry,complex<double> > x(RandBiped);
2794
2795 ArnoldiSolver<TransferMatrix<Symmetry,complex<double> >,TransferVector<Symmetry,complex<double> > > John(N,tol);
2796 if (dimK != -1)
2797 {
2798 John.set_dimK(dimK);
2799 }
2800 John.calc_dominant(T,x);
2801 lout << "Fixed point(" << label << "): GAUGE=" << g << ", DIR=" << DIR << ": " << John.info() << endl;
2802 //Normalize the Fixed point and try to make it real.
2803// x.data = exp(-1.i*arg(x.data.block[0](0,0))) * (1./x.data.norm()) * x.data;
2804
2805 res[0].first = John.get_lambda(0);
2806 res[0].second = x.data;
2807
2808 for (int n=0; n<N-1; ++n)
2809 {
2810 res[n+1].first = John.get_lambda(n+1);
2811 res[n+1].second = John.get_excited(n).data;
2812 }
2813
2814 if (abs(lambda1.imag()) > 1e1*tol)
2815 {
2816 lout << John.info() << endl;
2817 lout << termcolor::red << "Non-zero imaginary part of dominant eigenvalue λ=" << lambda1 << ", |λ|=" << abs(lambda1) << termcolor::reset << endl;
2818 }
2819
2820 lout << "norm test=" << x.data.norm() << endl;
2821
2822// LanczosSolver<TransferMatrix<Symmetry,double>,TransferVector<Symmetry,double>,double> Lutz(LANCZOS::REORTHO::FULL);
2823// Eigenstate<TransferVector<Symmetry,double>> z;
2824// z.state = TransferVector<Symmetry,double>(RandBipedr);;
2825// Lutz.edgeState(Tr, z, LANCZOS::EDGE::ROOF, 1e-7,1e-4, false);
2826//
2827// cout << Lutz.info() << endl;
2828// cout << "z.energy=" << z.energy << endl;
2829// lout << "z.norm test=" << z.state.data.norm() << endl;
2830
2831 for (int n=1; n<N; ++n)
2832 {
2833 res[n].first = John.get_lambda(n);
2834 }
2835
2836 lout << endl;
2837 // Note: corr.length ξ=-L/ln(|lambda[1]|")
2838
2839 return res;
2840}
2841
2842template<typename Symmetry, typename Scalar>
2843template<typename MpoScalar>
2844vector<std::pair<complex<double>,Tripod<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > > > Umps<Symmetry,Scalar>::
2845calc_dominant_Q (const Mpo<Symmetry,MpoScalar> &O, GAUGE::OPTION g, DMRG::DIRECTION::OPTION DIR, int N, double tol, int dimK, string label) const
2846{
2847 vector<std::pair<complex<double>,Tripod<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > > > res(N);
2848
2849 Umps<Symmetry,complex<double> > Compl = this->template cast<complex<double> > ();
2850 complex<double> lambda1;
2851
2852 auto Obs = O;
2853 Obs.transform_base(Qtarget(),false);
2854
2855 TransferMatrixQ<Symmetry,complex<double> > T;
2856 if (DIR == DMRG::DIRECTION::LEFT)
2857 {
2858 T = TransferMatrixQ<Symmetry,complex<double> >
2859 (VMPS::DIRECTION::RIGHT, Compl.A[g], Compl.A[g], locBasis(), Obs.Qtarget()); // contract from right to left
2860 }
2861 else
2862 {
2863 T = TransferMatrixQ<Symmetry,complex<double> >
2864 (VMPS::DIRECTION::LEFT, Compl.A[g], Compl.A[g], locBasis(), Obs.Qtarget()); // contract from left to right
2865 }
2866
2867 Tripod<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > Lid;
2868 Lid.setIdentity(Obs.inBasis(0).inner_dim(Symmetry::qvacuum()), 1, inBasis(0));
2869
2870 Tripod<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > Rid;
2871 Rid.setIdentity(Obs.outBasis(Obs.length()-1).inner_dim(Symmetry::qvacuum()), 1, outBasis(Obs.length()-1));
2872
2873 Tripod<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > TripodInit;
2874 Tripod<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > TripodInitTmp;
2875
2876 if (DIR == DMRG::DIRECTION::LEFT) // contract from right to left
2877 {
2878 contract_R(Rid, A[GAUGE::R][N_sites-1], Obs.W_at(N_sites-1), A[GAUGE::R][N_sites-1], Obs.locBasis(N_sites-1), Obs.opBasis(N_sites-1), TripodInitTmp);
2879 //lout << TripodInitTmp.print() << endl;
2880 // shift backward in cell
2881 for (int l=N_sites-2; l>=0; --l)
2882 {
2883 contract_R(TripodInitTmp, A[GAUGE::R][l], Obs.W_at(l), A[GAUGE::R][l], Obs.locBasis(l), Obs.opBasis(l), TripodInit);
2884 TripodInitTmp = TripodInit;
2885 //lout << TripodInitTmp.print() << endl;
2886 }
2887 }
2888 else // contract from left to right
2889 {
2890 contract_L(Lid, A[GAUGE::L][0], Obs.W_at(0), A[GAUGE::L][0], Obs.locBasis(0), Obs.opBasis(0), TripodInitTmp);
2891 // shift forward in cell
2892 for (size_t l=1; l<N_sites; ++l)
2893 {
2894 contract_L(TripodInitTmp, A[GAUGE::L][l], Obs.W_at(l), A[GAUGE::L][l], Obs.locBasis(l), Obs.opBasis(l), TripodInit);
2895 TripodInitTmp = TripodInit;
2896 }
2897 }
2898
2899 MpoTransferVector<Symmetry,complex<double> > x(TripodInit, make_pair(Obs.Qtarget(),0));
2900
2901 ArnoldiSolver<TransferMatrixQ<Symmetry,complex<double> >,MpoTransferVector<Symmetry,complex<double> > > Arnie(N,tol);
2902 if (dimK != -1) Arnie.set_dimK(dimK);
2903 Arnie.calc_dominant(T,x);
2904 lout << "Fixed point(" << label << "):" << " GAUGE=" << g << ", DIR=" << DIR << ": " << Arnie.info() << endl;
2905
2906 res[0].first = Arnie.get_lambda(0);
2907 res[0].second = x.data;
2908
2909 for (int n=0; n<N-1; ++n)
2910 {
2911 res[n+1].first = Arnie.get_lambda(n+1);
2912 res[n+1].second = Arnie.get_excited(n).data;
2913 }
2914
2915 if (abs(lambda1.imag()) > 1e1*tol)
2916 {
2917 lout << Arnie.info() << endl;
2918 lout << termcolor::red << "Non-zero imaginary part of dominant eigenvalue λ=" << lambda1 << ", |λ|=" << abs(lambda1) << termcolor::reset << endl;
2919 }
2920
2921 for (int n=1; n<N; ++n)
2922 {
2923 res[n].first = Arnie.get_lambda(n);
2924 }
2925
2926 lout << endl;
2927 // Note: corr.length ξ=-L/ln(|lambda[1]|")
2928
2929 return res;
2930}
2931
2932template<typename Symmetry, typename Scalar>
2933std::pair<complex<double>,Biped<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > > Umps<Symmetry,Scalar>::
2934calc_dominant_1symm (GAUGE::OPTION g, DMRG::DIRECTION::OPTION DIR, const Mpo<Symmetry,complex<double>> &R, bool TRANSPOSE, bool CONJUGATE) const
2935{
2936 Umps<Symmetry,complex<double> > Compl = this->template cast<complex<double> > ();
2937 complex<double> lambda;
2938
2939 TransferMatrix<Symmetry,complex<double> > T;
2940 if (DIR == DMRG::DIRECTION::LEFT)
2941 {
2942 // time_reverse(Compl.A[g],R,locBasis(),inBasis(),outBasis())
2943 T = TransferMatrix<Symmetry,complex<double> >
2944 (VMPS::DIRECTION::RIGHT, Compl.A[g], apply_symm(Compl.A[g],R,locBasis(),inBasis(),outBasis(),TRANSPOSE,CONJUGATE), locBasis());
2945 }
2946 else
2947 {
2948 //Compl.A[g]
2949 T = TransferMatrix<Symmetry,complex<double> >
2950 (VMPS::DIRECTION::LEFT, Compl.A[g], apply_symm(Compl.A[g],R,locBasis(),inBasis(),outBasis(),TRANSPOSE,CONJUGATE), locBasis());
2951 }
2952
2953 Biped<Symmetry,Matrix<complex<double>, Dynamic,Dynamic> > RandBiped;
2954 if (DIR == DMRG::DIRECTION::LEFT)
2955 {
2956 RandBiped.setRandom(inBasis(0), inBasis(0));
2957 }
2958 else
2959 {
2960 RandBiped.setRandom(outBasis(N_sites-1), outBasis(N_sites-1));
2961 }
2962 RandBiped = 1./RandBiped.norm() * RandBiped;
2963 TransferVector<Symmetry,complex<double> > x(RandBiped);
2964
2965 ArnoldiSolver<TransferMatrix<Symmetry,complex<double> >,TransferVector<Symmetry,complex<double> > > John(T,x,lambda);
2966
2967 lout << "fixed point, gauge=" << g << ", DIR=" << DIR << ": " << John.info() << endl;
2968 //Normalize the Fixed point and try to make it real.
2969// x.data = exp(-1.i*arg(x.data.block[0](0,0))) * (1./x.data.norm()) * x.data;
2970
2971 auto U = x.data.adjoint();
2972
2973 lout << boolalpha << "TRANSPOSE=" << TRANSPOSE << ", CONJUGATE=" << CONJUGATE << endl;
2974 lout << R.info() << endl;
2975// lout << "U.norm()=" << U.norm() << "\t" << U.adjoint().contract(U).trace() << endl;
2976 complex<double> O = (U.contract(U.conjugate())).trace();
2977 lout << "O raw result: " << O << endl;
2978
2979 if (abs(abs(lambda)-1.)>1e-2) O=0.;
2980 lout << termcolor::blue << "O=" << O.real() << termcolor::reset << endl << endl;
2981 return std::make_pair(O,x.data);
2982}
2983
2984template<typename Symmetry, typename Scalar>
2985std::pair<complex<double>,Biped<Symmetry,Matrix<complex<double>,Dynamic,Dynamic> > > Umps<Symmetry,Scalar>::
2986calc_dominant_2symm (GAUGE::OPTION g, DMRG::DIRECTION::OPTION DIR, const Mpo<Symmetry,complex<double>> &R1, const Mpo<Symmetry,complex<double>> &R2) const
2987{
2988 Umps<Symmetry,complex<double> > Compl = this->template cast<complex<double> > ();
2989 complex<double> lambda1, lambda2;
2990
2991 TransferMatrix<Symmetry,complex<double> > T1, T2;
2992 if (DIR == DMRG::DIRECTION::LEFT)
2993 {
2994 T1 = TransferMatrix<Symmetry,complex<double> >
2995 (VMPS::DIRECTION::RIGHT, Compl.A[g], apply_symm(Compl.A[g],R1,locBasis(),inBasis(),outBasis()), locBasis());
2996 T2 = TransferMatrix<Symmetry,complex<double> >
2997 (VMPS::DIRECTION::RIGHT, Compl.A[g], apply_symm(Compl.A[g],R2,locBasis(),inBasis(),outBasis()), locBasis());
2998 }
2999 else
3000 {
3001 //Compl.A[g]
3002 T1 = TransferMatrix<Symmetry,complex<double> >
3003 (VMPS::DIRECTION::LEFT, Compl.A[g], apply_symm(Compl.A[g],R1,locBasis(),inBasis(),outBasis()), locBasis());
3004 T2 = TransferMatrix<Symmetry,complex<double> >
3005 (VMPS::DIRECTION::LEFT, Compl.A[g], apply_symm(Compl.A[g],R2,locBasis(),inBasis(),outBasis()), locBasis());
3006 }
3007
3008 Biped<Symmetry,Matrix<complex<double>, Dynamic,Dynamic> > RandBiped;
3009 if (DIR == DMRG::DIRECTION::LEFT)
3010 {
3011 RandBiped.setRandom(inBasis(0), inBasis(0));
3012 }
3013 else
3014 {
3015 RandBiped.setRandom(outBasis(N_sites-1), outBasis(N_sites-1));
3016 }
3017 RandBiped = 1./RandBiped.norm() * RandBiped;
3018
3019 TransferVector<Symmetry,complex<double> > x1(RandBiped);
3020 TransferVector<Symmetry,complex<double> > x2(RandBiped);
3021
3022 ArnoldiSolver<TransferMatrix<Symmetry,complex<double> >,TransferVector<Symmetry,complex<double> > > John, Jane;
3023 #pragma omp sections
3024 {
3025 #pragma omp section
3026 {
3027 John.calc_dominant(T1,x1,lambda1);
3028 }
3029 #pragma omp section
3030 {
3031 Jane.calc_dominant(T2,x2,lambda2);
3032 }
3033 }
3034
3035 lout << "fixed point, gauge=" << g << ", DIR=" << DIR << ": " << John.info() << endl;
3036 lout << "fixed point, gauge=" << g << ", DIR=" << DIR << ": " << Jane.info() << endl;
3037
3038 auto U1 = x1.data.adjoint();
3039 auto U2 = x2.data.adjoint();
3040
3041 lout << R1.info() << endl;
3042 lout << R2.info() << endl;
3043// lout << "U1.norm()=" << U1.norm() << "\t" << U1.adjoint().contract(U1).trace() << "\t" << U1.contract(U1.adjoint()).trace() << endl;
3044// lout << "U2.norm()=" << U2.norm() << "\t" << U2.adjoint().contract(U2).trace() << "\t" << U1.contract(U1.adjoint()).trace() << endl;
3045
3046 complex<double> O12 = (U1.contract(U2.contract(U1.adjoint().contract(U2.adjoint())))).trace();
3047 O12 *= double(U1.block[0].rows());
3048 // Note: Pollmann normalizes U*Udag=Id, tr(U*Udag)=Chi; we normalize U*Udag=1/Chi, tr(U*Udag)=1
3049 lout << "O12 raw result=" << O12 << endl;
3050
3051 lout << "commut=" << U1.block[0].rows()*(U1.block[0]*U2.block[0]-U2.block[0]*U1.block[0]).norm() << endl;
3052 lout << "anticommut=" << U1.block[0].rows()*(U1.block[0]*U2.block[0]+U2.block[0]*U1.block[0]).norm() << endl;
3053
3054 if (abs(abs(lambda1)-1.)>1e-2 or abs(abs(lambda2)-1.)>1e-2) O12=0.;
3055 lout << termcolor::blue << "O12=" << O12.real() << termcolor::reset << endl << endl;
3056 return std::make_pair(O12,x1.data);
3057}
3058
3059template<typename Symmetry, typename Scalar>
3060std::pair<vector<qarray<Symmetry::Nq> >, ArrayXd> Umps<Symmetry,Scalar>::
3061entanglementSpectrumLoc (size_t loc) const
3062{
3063 vector<pair<qarray<Nq>, double> > Svals;
3064 for (const auto &x : SVspec[loc])
3065 for (int i=0; i<std::get<0>(x.second).size(); ++i)
3066 {
3067 Svals.push_back(std::make_pair(x.first,std::get<0>(x.second)(i)));
3068 }
3069 sort(Svals.begin(), Svals.end(), [] (const pair<qarray<Nq>, double> &p1, const pair<qarray<Nq>, double> &p2) { return p2.second < p1.second;});
3070
3071 ArrayXd Sout(Svals.size());
3072 vector<qarray<Nq> > Qout(Svals.size());
3073 for (int i=0; i<Svals.size(); ++i)
3074 {
3075 Sout(i) = Svals[i].second;
3076 Qout[i] = Svals[i].first;
3077 }
3078 return std::make_pair(Qout,Sout);
3079}
3080
3081template<typename Symmetry, typename Scalar>
3082template<typename MpoScalar>
3083ArrayXXcd Umps<Symmetry,Scalar>::
3084intercellSF (const Mpo<Symmetry,MpoScalar> &Oalfa, const Mpo<Symmetry,MpoScalar> &Obeta, int Lx, double kmin, double kmax, int kpoints, DMRG::VERBOSITY::OPTION VERB, double tol)
3085{
3086 double t_tot=0.;
3087 double t_LReigen=0.;
3088 double t_GMRES=0.;
3089 double t_contraction=0.;
3090
3091 Stopwatch<> TotTimer;
3092
3093 Biped<Symmetry,Matrix<complex<Scalar>,Dynamic,Dynamic> > Reigen_LR, Leigen_LR, Reigen_RL, Leigen_RL;
3094
3095 Stopwatch<> LReigenTimer;
3096
3097 // T_L^R, right eigenvector
3098 Reigen_LR = calc_LReigen(VMPS::DIRECTION::RIGHT, A[GAUGE::L], A[GAUGE::R], outBasis(N_sites-1), outBasis(N_sites-1), qloc, 100ul, tol).state;
3099 // T_L^R, left eigenvector
3100 Leigen_LR = calc_LReigen(VMPS::DIRECTION::LEFT, A[GAUGE::L], A[GAUGE::R], inBasis(0), inBasis(0), qloc, 100ul, tol).state;
3101 // T_R^L, right eigenvector
3102 Reigen_RL = calc_LReigen(VMPS::DIRECTION::RIGHT, A[GAUGE::R], A[GAUGE::L], outBasis(N_sites-1), outBasis(N_sites-1), qloc, 100ul, tol).state;
3103 // T_R^L, left eigenvector
3104 Leigen_RL = calc_LReigen(VMPS::DIRECTION::LEFT, A[GAUGE::R], A[GAUGE::L], inBasis(0), inBasis(0), qloc, 100ul, tol).state;
3105
3106 t_LReigen += LReigenTimer.time();
3107
3108 // b (edge tensor of contraction) for alfa, beta and exp(-i*Lcell*k), exp(+i*Lcell*k)
3109 // Note: AC is set to the locality of beta in the bra state and to alfa in the ket state
3110
3111 Stopwatch<> ContractionTimer;
3112
3113 lout << Oalfa.info() << endl;
3114 lout << Obeta.info() << endl;
3115 //lout << Oalfa.print(true) << endl;
3116
3117 Tripod<Symmetry,Matrix<MpoScalar,Dynamic,Dynamic> > Lid; Lid.setIdentity(1,1,inBasis(0));
3118 Tripod<Symmetry,Matrix<MpoScalar,Dynamic,Dynamic> > Rid; Rid.setIdentity(1,1,outBasis(N_sites-1));
3119
3120 // term exp(-i*Lcell*k), alfa
3121 vector<Tripod<Symmetry,Matrix<MpoScalar,Dynamic,Dynamic> > > bmalfaTripod(N_sites);
3122 contract_L(Lid, A[GAUGE::L][0], Oalfa.W_at(0), A[GAUGE::C][0],
3123 Oalfa.locBasis(0), Oalfa.opBasis(0), bmalfaTripod[0]);
3124 // shift forward in cell
3125 for (size_t l=1; l<N_sites; ++l)
3126 {
3127 contract_L(bmalfaTripod[l-1], A[GAUGE::L][l], Oalfa.W_at(l), A[GAUGE::R][l],
3128 Oalfa.locBasis(l), Oalfa.opBasis(l), bmalfaTripod[l]);
3129 }
3130
3131 // term exp(+i*Lcell*k), alfa
3132 vector<Tripod<Symmetry,Matrix<MpoScalar,Dynamic,Dynamic> > > bpalfaTripod(N_sites);
3133 contract_R(Rid, A[GAUGE::R][N_sites-1], Oalfa.reversed.W[N_sites-1], A[GAUGE::C][N_sites-1],
3134 Oalfa.locBasis(N_sites-1), Oalfa.opBasis(N_sites-1), bpalfaTripod[N_sites-1]);
3135 // shift backward in cell
3136 for (int l=N_sites-2; l>=0; --l)
3137 {
3138 contract_R(bpalfaTripod[l+1], A[GAUGE::R][l], Oalfa.reversed.W[l], A[GAUGE::L][l],
3139 Oalfa.locBasis(l), Oalfa.opBasis(l), bpalfaTripod[l]);
3140 }
3141 assert(bpalfaTripod[0].size() > 0);
3142
3143 // term exp(-i*Lcell*k), beta
3144 vector<Tripod<Symmetry,Matrix<MpoScalar,Dynamic,Dynamic> > > bmbetaTripod(N_sites);
3145 contract_R(Rid, A[GAUGE::C][N_sites-1], Obeta.reversed.W[N_sites-1], A[GAUGE::R][N_sites-1],
3146 Obeta.locBasis(N_sites-1), Obeta.opBasis(N_sites-1), bmbetaTripod[N_sites-1]);
3147 // shift backward in cell
3148 for (int l=N_sites-2; l>=0; --l)
3149 {
3150 contract_R(bmbetaTripod[l+1], A[GAUGE::L][l], Obeta.reversed.W[l], A[GAUGE::R][l],
3151 Obeta.locBasis(l), Obeta.opBasis(l), bmbetaTripod[l]);
3152 }
3153 assert(bmbetaTripod[0].size() > 0);
3154
3155 // term exp(+i*Lcell*k), beta
3156 vector<Tripod<Symmetry,Matrix<MpoScalar,Dynamic,Dynamic> > > bpbetaTripod(N_sites);
3157 contract_L(Lid, A[GAUGE::C][0], Obeta.W_at(0), A[GAUGE::L][0],
3158 Obeta.locBasis(0), Obeta.opBasis(0), bpbetaTripod[0]);
3159 // shift forward in cell
3160 for (size_t l=1; l<N_sites; ++l)
3161 {
3162 contract_L(bpbetaTripod[l-1], A[GAUGE::R][l], Obeta.W_at(l), A[GAUGE::L][l],
3163 Obeta.locBasis(l), Obeta.opBasis(l), bpbetaTripod[l]);
3164 }
3165
3166 // wrap bmalfa, bpalfa by MpoTransferVector for GMRES
3167 // Note: the Tripods has only a single quantum number with inner dimension 1 on their mid leg. We need to pass this information to MpoTransferVector.
3168 MpoTransferVector<Symmetry,complex<Scalar> > bmalfa(bmalfaTripod[N_sites-1].template cast<MatrixXcd >(), make_pair(bmalfaTripod[N_sites-1].mid(0),0));
3169 MpoTransferVector<Symmetry,complex<Scalar> > bpalfa(bpalfaTripod[0].template cast<MatrixXcd>(), make_pair(bpalfaTripod[0].mid(0),0));
3170
3171 // cast bmbeta, bpbeta to complex Tripod for final contraction
3172 Tripod<Symmetry,Matrix<complex<Scalar>,Dynamic,Dynamic> > bmbeta = bmbetaTripod[0].template cast<MatrixXcd >();
3173 Tripod<Symmetry,Matrix<complex<Scalar>,Dynamic,Dynamic> > bpbeta = bpbetaTripod[N_sites-1].template cast<MatrixXcd>();
3174
3175 t_contraction += ContractionTimer.time();
3176
3177 ArrayXXcd out(kpoints,2);
3178 if (kmin==kmax) {out.resize(1,2);} // only one k-point needed in case of kmin=kmax
3179
3180 // solve linear systems
3181 Stopwatch<> GMRES_Timer;
3182 #pragma omp parallel for
3183 for (int ik=0; ik<out.rows(); ++ik)
3184 {
3185 // the last k-point is repeated, therefore kpoints-1 independent points:
3186 double kval = (kmin==kmax)? kmin : kmin + ik*(kmax-kmin)/(kpoints-1);
3187
3188 GMResSolver<TransferMatrixSF<Symmetry,Scalar>,MpoTransferVector<Symmetry,complex<Scalar> > > Gimli;
3189
3190 // term exp(-i*Lcell*k)
3191 TransferMatrixSF<Symmetry,Scalar> Tm(VMPS::DIRECTION::LEFT, A[GAUGE::L], A[GAUGE::R], Leigen_LR, Reigen_LR, qloc, Lx*kval, bmalfaTripod[N_sites-1].mid(0));
3192 Gimli.set_dimK(min(100ul,dim(bmalfa)));
3193 assert(dim(bmalfa) > 0);
3194 MpoTransferVector<Symmetry,complex<Scalar> > Fmalfa;
3195 Gimli.solve_linear(Tm, bmalfa, Fmalfa, tol, true);
3196 if (VERB >= DMRG::VERBOSITY::STEPWISE)
3197 {
3198 lout << ik << ", k/Ï€=" << kval/M_PI << ", term exp(-i*Lcell*k), " << Gimli.info() << "; dim(bmalfa)=" << dim(bmalfa) << endl;
3199 }
3200
3201 // term exp(+i*Lcell*k)
3202 TransferMatrixSF<Symmetry,Scalar> Tp(VMPS::DIRECTION::RIGHT, A[GAUGE::R], A[GAUGE::L], Leigen_RL, Reigen_RL, qloc, Lx*kval, bpalfaTripod[0].mid(0));
3203 Gimli.set_dimK(min(100ul,dim(bpalfa)));
3204 assert(dim(bpalfa) > 0);
3205 MpoTransferVector<Symmetry,complex<Scalar> > Fpalfa;
3206 Gimli.solve_linear(Tp, bpalfa, Fpalfa, tol, true);
3207 if (VERB >= DMRG::VERBOSITY::STEPWISE)
3208 {
3209 lout << ik << ", k/Ï€=" << kval/M_PI << ", term exp(+i*Lcell*k), " << Gimli.info() << "; dim(bpalfa)=" << dim(bpalfa) << endl;
3210 }
3211
3212 complex<double> resm = contract_LR(Fmalfa.data, bmbeta);
3213 complex<double> resp = contract_LR(bpbeta, Fpalfa.data);
3214 // cout << "resm=" << resm << ", resp=" << resp << endl;
3215
3216 // result
3217 out(ik,0) = kval;
3218 out(ik,1) = exp(-1.i*static_cast<double>(Lx)*kval) * resm + exp(+1.i*static_cast<double>(Lx)*kval) * resp;
3219 }
3220
3221 t_GMRES += GMRES_Timer.time();
3222
3223 t_tot = TotTimer.time();
3224
3225 if (VERB >= DMRG::VERBOSITY::ON_EXIT)
3226 {
3227 lout << TotTimer.info("StructureFactor")
3228 << " (LReigen=" << round(t_LReigen/t_tot*100.,0) << "%, "
3229 << "GMRES=" << round(t_GMRES/t_tot*100.,0) << "%, "
3230 << "contractions=" << round(t_contraction/t_tot*100.,0) << "%)"
3231 << ", kmin/Ï€=" << kmin/M_PI << ", kmax/Ï€=" << kmax/M_PI << ", kpoints=" << out.rows() << endl;
3232 lout << "\t" << Oalfa.info() << endl;
3233 lout << "\t" << Obeta.info() << endl;
3234 }
3235
3236 return out;
3237}
3238
3239template<typename Symmetry, typename Scalar>
3240template<typename MpoScalar>
3241complex<Scalar> Umps<Symmetry,Scalar>::
3242intercellSFpoint (const Mpo<Symmetry,MpoScalar> &Oalfa, const Mpo<Symmetry,MpoScalar> &Obeta, int Lx, double kval, DMRG::VERBOSITY::OPTION VERB)
3243{
3244 ArrayXXcd res = intercellSF(Oalfa, Obeta, Lx, kval, kval, 1, VERB);
3245 return res(0,1);
3246}
3247
3248template<typename Symmetry, typename Scalar>
3249template<typename MpoScalar>
3250complex<Scalar> Umps<Symmetry,Scalar>::
3251SFpoint (const ArrayXXcd &cellAvg, const vector<Mpo<Symmetry,MpoScalar> > &Oalfa, const vector<Mpo<Symmetry,MpoScalar> > &Obeta,
3252 int Lx, double kval, DMRG::VERBOSITY::OPTION VERB)
3253{
3254 assert(Oalfa.size() == Lx and Obeta.size() == Lx and cellAvg.rows() == Lx and cellAvg.cols() == Lx);
3255
3256 complex<double> res = 0;
3257
3258 ArrayXXcd Sijk = cellAvg;
3259
3260 #ifndef UMPS_DONT_PARALLELIZE_SF_LOOPS
3261 #pragma omp parallel for collapse(2)
3262 #endif
3263 for (size_t i0=0; i0<Lx; ++i0)
3264 for (size_t j0=0; j0<Lx; ++j0)
3265 {
3266 Sijk(i0,j0) += intercellSFpoint(Oalfa[i0],Obeta[j0], Lx, kval, VERB);
3267 }
3268
3269 for (size_t i0=0; i0<Lx; ++i0)
3270 for (size_t j0=0; j0<Lx; ++j0)
3271 {
3272 // Careful: Must first convert to double and then subtract, since the difference can become negative!
3273 res += 1./static_cast<double>(Lx) * exp(-1.i*kval*(static_cast<double>(i0)-static_cast<double>(j0))) * Sijk(j0,i0);
3274 // Attention: order (j0,i0) in argument is correct!
3275 }
3276
3277 return res;
3278}
3279
3280template<typename Symmetry, typename Scalar>
3281template<typename MpoScalar>
3282ArrayXXcd Umps<Symmetry,Scalar>::
3283SF (const ArrayXXcd &cellAvg, const vector<Mpo<Symmetry,MpoScalar> > &Oalfa, const vector<Mpo<Symmetry,MpoScalar> > &Obeta,
3284 int Lx, double kmin, double kmax, int kpoints, DMRG::VERBOSITY::OPTION VERB, double tol)
3285{
3286 assert(Oalfa.size() == Lx and Obeta.size() == Lx and cellAvg.rows() == Lx and cellAvg.cols() == Lx);
3287
3288 vector<vector<ArrayXXcd> > Sijk(Lx);
3289 for (size_t i0=0; i0<Lx; ++i0)
3290 {
3291 Sijk[i0].resize(Lx);
3292 for (size_t j0=0; j0<Lx; ++j0)
3293 {
3294 Sijk[i0][j0].resize(kpoints,2);
3295 Sijk[i0][j0] = 0;
3296 }
3297 }
3298
3299 #ifndef UMPS_DONT_PARALLELIZE_SF_LOOPS
3300 #pragma omp parallel for collapse(2)
3301 #endif
3302 for (size_t i0=0; i0<Lx; ++i0)
3303 for (size_t j0=0; j0<Lx; ++j0)
3304 {
3305 Sijk[i0][j0] = intercellSF(Oalfa[i0],Obeta[j0],Lx,kmin,kmax,kpoints,VERB,tol);
3306 Sijk[i0][j0].col(1) += cellAvg(i0,j0);
3307 }
3308
3309 ArrayXXcd res(kpoints,2); res=0;
3310
3311 for (size_t ik=0; ik<kpoints; ++ik)
3312 for (size_t i0=0; i0<Lx; ++i0)
3313 for (size_t j0=0; j0<Lx; ++j0)
3314 {
3315 double kval = Sijk[i0][j0](ik,0).real();
3316 res(ik,0) = kval;
3317 // Careful: Must first convert to double and then subtract, since the difference can become negative!
3318 //res(ik,1) += 1./static_cast<double>(Lx) * exp(-1.i*kval*(static_cast<double>(i0)-static_cast<double>(j0))) * Sijk[i0][j0](ik,1);
3319 res(ik,1) += 1./static_cast<double>(Lx) * exp(-1.i*kval*(static_cast<double>(i0)-static_cast<double>(j0))) * Sijk[j0][i0](ik,1);
3320 // Attention: order [j0][i0] in argument is correct!
3321 // Or maybe not?...
3322 }
3323
3324 return res;
3325}
3326
3327#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_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
sum
Hamiltonian sum(const Hamiltonian &H1, const Hamiltonian &H2, DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: DmrgLinearAlgebra.h:866
diff
Mpo< Symmetry, Scalar > diff(const Mpo< Symmetry, Scalar > &H1, const Mpo< Symmetry, Scalar > &H2, DMRG::VERBOSITY::OPTION VERBOSITY=DMRG::VERBOSITY::SILENT)
Definition: DmrgLinearAlgebra.h:914
norm
double norm(const PivotMatrix0< Symmetry, Scalar, MpoScalar > &H)
Definition: DmrgPivotMatrix0.h:114
dim
size_t dim(const PivotMatrix0< Symmetry, Scalar, MpoScalar > &H)
Definition: DmrgPivotMatrix0.h:107
squaredNorm
double squaredNorm(const PivotVector< Symmetry, Scalar_ > &V)
Definition: DmrgPivotVector.h:342
A
@ A
Definition: DmrgTypedefs.h:130
apply_symm
vector< vector< Biped< Symmetry, Matrix< complex< Scalar >, Dynamic, Dynamic > > > > apply_symm(const vector< vector< Biped< Symmetry, Matrix< complex< Scalar >, Dynamic, Dynamic > > > > &A, const Mpo< Symmetry, complex< Scalar > > &R, const vector< vector< qarray< Symmetry::Nq > > > &qloc, const vector< Qbasis< Symmetry > > &qauxAl, const vector< Qbasis< Symmetry > > &qauxAr, bool TRANSPOSE=false, bool CONJUGATE=false)
Definition: Umps.h:1694
calc_LReigen
Eigenstate< Biped< Symmetry, Matrix< complex< Scalar >, Dynamic, Dynamic > > > calc_LReigen(VMPS::DIRECTION::OPTION DIR, const vector< vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > > &Abra, const vector< vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > > &Aket, const Qbasis< Symmetry > &basisBra, const Qbasis< Symmetry > &basisKet, const vector< vector< qarray< Symmetry::Nq > > > &qloc, size_t dimK=100ul, double tol_input=1e-12, Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > *LReigenTop=NULL)
Definition: VumpsContractions.h:17
Blocker::Aclump
Biped< Symmetry, MatrixType > Aclump(DMRG::DIRECTION::OPTION DIR)
Definition: Blocker.h:26
Blocker::reblock
vector< Biped< Symmetry, MatrixType > > reblock(const Biped< Symmetry, MatrixType > &B, DMRG::DIRECTION::OPTION DIR)
Definition: Blocker.h:58
MpoTerms::W_at
const std::vector< std::vector< std::vector< Biped< Symmetry, MatrixType > > > > & W_at(const std::size_t loc) const
Definition: MpoTerms.h:699
MpoTerms::locBasis
const std::vector< std::vector< qType > > & locBasis() const
Definition: MpoTerms.h:702
MpoTerms::transform_base
void transform_base(const qType &qShift, const bool PRINT=false, const int factor=-1, const std::size_t powre=0ul)
Definition: MpoTerms.h:3113
MpoTerms::reversed
reversedData reversed
Definition: MpoTerms.h:409
MpoTerms::opBasis
const std::vector< std::vector< qType > > & opBasis() const
Definition: MpoTerms.h:710
Mpo
Definition: Mpo.h:40
Mpo::info
std::string info(bool REDUCED=false) const
MpsCompressor::loc2
size_t loc2() const
Definition: MpsCompressor.h:208
MpsCompressor::MatrixType
Matrix< Scalar, Dynamic, Dynamic > MatrixType
Definition: MpsCompressor.h:39
MpsCompressor::info
string info() const
Definition: MpsCompressor.h:237
MpsCompressor::memory
double memory(MEMUNIT memunit=GB) const
Definition: MpsCompressor.h:291
MpsCompressor::loc1
size_t loc1() const
Definition: MpsCompressor.h:207
Qbasis
Definition: Qbasis.h:39
Qbasis::leftAmount
Eigen::Index leftAmount(const qType &qnew, const std::array< qType, 2 > &qold) const
Definition: Qbasis.h:421
Qbasis::find
qType find(const std::string &ident) const
Definition: Qbasis.h:297
Qbasis::combine
Qbasis< Symmetry > combine(const Qbasis< Symmetry > &other, bool FLIP=false) const
Definition: Qbasis.h:686
Qbasis::add
Qbasis< Symmetry > add(const Qbasis< Symmetry > &other) const
Definition: Qbasis.h:640
Qbasis::pullData
void pullData(const vector< Biped< Symmetry, MatrixType > > &A, const Eigen::Index &leg)
Definition: Qbasis.h:486
Qbasis::history
std::unordered_map< qType, std::vector< fuseData > > history
Definition: Qbasis.h:209
Qbasis::inner_dim
Eigen::Index inner_dim(const Eigen::Index &num_in) const
Definition: Qbasis.h:391
Umps
Definition: Umps.h:42
Umps::S
VectorXd S
Definition: Umps.h:376
Umps::Qtop
qarray< Symmetry::Nq > Qtop(size_t loc) const
Definition: Umps.h:542
Umps::Qbot
qarray< Symmetry::Nq > Qbot(size_t loc) const
Definition: Umps.h:549
Umps::dot
double dot(const Umps< Symmetry, Scalar > &Vket) const
Definition: Umps.h:1122
Umps::cast
Umps< Symmetry, OtherScalar > cast() const
Definition: Umps.h:1973
Umps::locBasis
vector< qarray< Symmetry::Nq > > locBasis(size_t loc) const
Definition: Umps.h:129
Umps::qloc
vector< vector< qarray< Symmetry::Nq > > > qloc
Definition: Umps.h:365
Umps::A_at
const vector< Biped< Symmetry, MatrixType > > & A_at(GAUGE::OPTION g, size_t loc) const
Definition: Umps.h:201
Umps::Umps
Umps()
Definition: Umps.h:52
Umps::graph
void graph(string filename) const
Definition: Umps.h:900
Umps::N_sites
size_t N_sites
Definition: Umps.h:346
Umps::get_frst_rows
size_t get_frst_rows() const
Definition: Umps.h:141
Umps::Mmax
size_t Mmax
Definition: Umps.h:347
Umps::resize_arrays
void resize_arrays()
Definition: Umps.h:556
Umps::length
size_t length() const
Definition: Umps.h:147
Umps::calc_N
void calc_N(DMRG::DIRECTION::OPTION DIR, size_t loc, vector< Biped< Symmetry, MatrixType > > &N) const
Definition: Umps.h:1782
Umps::entanglementSpectrum
vector< map< qarray< Nq >, tuple< ArrayXd, int > > > entanglementSpectrum() const
Definition: Umps.h:112
Umps::calc_entropy
void calc_entropy(bool PRINT=false)
Definition: Umps.h:359
Umps::real
Umps< Symmetry, double > real() const
Definition: Umps.h:1999
Umps::inBasis
Qbasis< Symmetry > inBasis(size_t loc) const
Definition: Umps.h:133
Umps::updateC
void updateC(size_t loc)
Definition: Umps.h:2317
Umps::inbase
vector< Qbasis< Symmetry > > inbase
Definition: Umps.h:381
Umps::svdDecompose
void svdDecompose(size_t loc, GAUGE::OPTION gauge=GAUGE::C)
Definition: Umps.h:1513
Umps::intercellSF
ArrayXXcd intercellSF(const Mpo< Symmetry, MpoScalar > &Oalfa, const Mpo< Symmetry, MpoScalar > &Obeta, int Lx, double kmin=0., double kmax=2.*M_PI, int kpoints=51, DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::ON_EXIT, double tol=1e-12)
Definition: Umps.h:3084
Umps::truncate
void truncate(bool SET_AC_RANDOM=true)
Definition: Umps.h:2703
Umps::outBasis
vector< Qbasis< Symmetry > > outBasis() const
Definition: Umps.h:138
Umps::polarDecompose
void polarDecompose(size_t loc, GAUGE::OPTION gauge=GAUGE::C)
Definition: Umps.h:1358
Umps::min_Nsv
size_t min_Nsv
Definition: Umps.h:350
Umps::Nq
static constexpr size_t Nq
Definition: Umps.h:44
Umps::max_Nsv
size_t max_Nsv
Definition: Umps.h:350
Umps::Nqmax
size_t Nqmax
Definition: Umps.h:347
Umps::memory
double memory(MEMUNIT memunit) const
Definition: Umps.h:457
Umps::max_Nrich
int max_Nrich
Definition: Umps.h:351
Umps::calc_Nqmax
size_t calc_Nqmax() const
Definition: Umps.h:474
Umps::calc_entropy
void calc_entropy(size_t loc, bool PRINT=false)
Definition: Umps.h:1157
Umps::Qtarget
qarray< Nq > Qtarget() const
Definition: Umps.h:211
Umps::calc_fullMmax
size_t calc_fullMmax() const
Definition: Umps.h:513
Umps::Qtot
qarray< Nq > Qtot
Definition: Umps.h:353
Umps::outBasis
Qbasis< Symmetry > outBasis(size_t loc) const
Definition: Umps.h:137
Umps::update_inbase
void update_inbase(GAUGE::OPTION g=GAUGE::C)
Definition: Umps.h:387
Umps::SVspec
vector< map< qarray< Nq >, tuple< ArrayXd, int > > > SVspec
Definition: Umps.h:378
Umps::enrich
void enrich(size_t loc, GAUGE::OPTION g, const vector< Biped< Symmetry, MatrixType > > &P)
Definition: Umps.h:2375
Umps::adjustQN
void adjustQN(const size_t number_cells)
Definition: Umps.h:2025
Umps::calc_epsLRsq
Scalar calc_epsLRsq(GAUGE::OPTION gauge, size_t loc) const
Definition: Umps.h:1224
Umps::get_last_cols
size_t get_last_cols() const
Definition: Umps.h:144
Umps::calc_dominant_Q
vector< std::pair< complex< double >, Tripod< Symmetry, Matrix< complex< double >, Dynamic, Dynamic > > > > calc_dominant_Q(const Mpo< Symmetry, MpoScalar > &O, GAUGE::OPTION g=GAUGE::R, DMRG::DIRECTION::OPTION DIR=DMRG::DIRECTION::RIGHT, int N=2, double tol=1e-15, int dimK=-1, string label="") const
Definition: Umps.h:2845
Umps::minus1modL
size_t minus1modL(size_t l) const
Definition: Umps.h:208
Umps::truncWeight
ArrayXd truncWeight
Definition: Umps.h:362
Umps::SF
ArrayXXcd SF(const ArrayXXcd &cellAvg, const vector< Mpo< Symmetry, MpoScalar > > &Oalfa, const vector< Mpo< Symmetry, MpoScalar > > &Obeta, int Lx, double kmin, double kmax, int kpoints, DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::ON_EXIT, double tol=1e-12)
Definition: Umps.h:3283
Umps::test_ortho
string test_ortho(double tol=1e-6) const
Definition: Umps.h:972
Umps::setRandom
void setRandom()
Definition: Umps.h:863
Umps::updateAC
void updateAC(size_t loc, GAUGE::OPTION g)
Definition: Umps.h:2346
Umps::update_inbase
void update_inbase(size_t loc, GAUGE::OPTION g=GAUGE::C)
Definition: Umps.h:526
Umps::resize
void resize(size_t Mmax_input, size_t Nqmax_input, bool INIT_TO_HALF_INTEGER_SPIN)
Definition: Umps.h:623
Umps::orthogonalize_right
void orthogonalize_right(GAUGE::OPTION g, vector< Biped< Symmetry, MatrixType > > &G_R)
Definition: Umps.h:2575
Umps::eps_truncWeight
double eps_truncWeight
Definition: Umps.h:349
Umps::A
std::array< vector< vector< Biped< Symmetry, MatrixType > > >, 3 > A
Definition: Umps.h:368
Umps::outbase
vector< Qbasis< Symmetry > > outbase
Definition: Umps.h:382
Umps::intercellSFpoint
complex< Scalar > intercellSFpoint(const Mpo< Symmetry, MpoScalar > &Oalfa, const Mpo< Symmetry, MpoScalar > &Obeta, int Lx, double kval, DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::ON_EXIT)
Definition: Umps.h:3242
Umps::calc_Mmax
size_t calc_Mmax() const
Definition: Umps.h:500
Umps::calc_Dmax
size_t calc_Dmax() const
Definition: Umps.h:487
Umps::eps_svd
double eps_svd
Definition: Umps.h:348
Umps::calc_dominant
vector< std::pair< complex< double >, Biped< Symmetry, Matrix< complex< double >, Dynamic, Dynamic > > > > calc_dominant(GAUGE::OPTION g=GAUGE::R, DMRG::DIRECTION::OPTION DIR=DMRG::DIRECTION::RIGHT, int N=2, double tol=1e-15, int dimK=-1, qarray< Symmetry::Nq > Qtot=Symmetry::qvacuum(), string label="") const
Definition: Umps.h:2764
Umps::SFpoint
complex< Scalar > SFpoint(const ArrayXXcd &cellAvg, const vector< Mpo< Symmetry, MpoScalar > > &Oalfa, const vector< Mpo< Symmetry, MpoScalar > > &Obeta, int Lx, double kval, DMRG::VERBOSITY::OPTION VERB=DMRG::VERBOSITY::ON_EXIT)
Definition: Umps.h:3251
Umps::update_outbase
void update_outbase(GAUGE::OPTION g=GAUGE::C)
Definition: Umps.h:388
Umps::normalize_C
void normalize_C()
Definition: Umps.h:852
Umps::inBasis
vector< Qbasis< Symmetry > > inBasis() const
Definition: Umps.h:134
Umps::entanglementSpectrumLoc
std::pair< vector< qarray< Symmetry::Nq > >, ArrayXd > entanglementSpectrumLoc(size_t loc) const
Definition: Umps.h:3061
Umps::entropy
VectorXd entropy() const
Definition: Umps.h:109
Umps::update_outbase
void update_outbase(size_t loc, GAUGE::OPTION g=GAUGE::C)
Definition: Umps.h:534
Umps::C
vector< Biped< Symmetry, MatrixType > > C
Definition: Umps.h:371
Umps::calc_dominant_2symm
std::pair< complex< double >, Biped< Symmetry, Matrix< complex< double >, Dynamic, Dynamic > > > calc_dominant_2symm(GAUGE::OPTION g, DMRG::DIRECTION::OPTION DIR, const Mpo< Symmetry, complex< double > > &R1, const Mpo< Symmetry, complex< double > > &R2) const
Definition: Umps.h:2986
Umps::sort_A
void sort_A(size_t loc, GAUGE::OPTION g, bool SORT_ALL_GAUGES=false)
Definition: Umps.h:2296
Umps::locBasis
vector< vector< qarray< Symmetry::Nq > > > locBasis() const
Definition: Umps.h:130
Umps::orthogonalize_left
void orthogonalize_left(GAUGE::OPTION g, vector< Biped< Symmetry, MatrixType > > &G_L)
Definition: Umps.h:2640
Umps::info
string info() const
Definition: Umps.h:393
Umps::MatrixType
Matrix< Scalar, Dynamic, Dynamic > MatrixType
Definition: Umps.h:43
Umps::calc_dominant_1symm
std::pair< complex< double >, Biped< Symmetry, Matrix< complex< double >, Dynamic, Dynamic > > > calc_dominant_1symm(GAUGE::OPTION g, DMRG::DIRECTION::OPTION DIR, const Mpo< Symmetry, complex< double > > &R, bool TRANSPOSE, bool CONJUGATE) const
Definition: Umps.h:2934
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
Definition: Biped.h:47
contract::OORR
@ OORR
Definition: Biped.h:48
qarray2
std::array< qarray< Nq >, 2 > qarray2
Definition: qarray.h:51
Biped
Definition: Biped.h:64
Biped::adjoint
Biped< Symmetry, MatrixType_ > adjoint() const
Definition: Biped.h:630
Biped::block
std::vector< MatrixType_ > block
Definition: Biped.h:96
Biped::in
std::vector< qType > in
Definition: Biped.h:87
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::operatorNorm
double operatorNorm(bool COLWISE=true) const
Definition: Biped.h:494
Biped::setRandom
void setRandom()
Definition: Biped.h:343
Biped::norm
double norm() const
Definition: Biped.h:510
Biped::dim
std::size_t dim
Definition: Biped.h:82
Biped::out
std::vector< qType > out
Definition: Biped.h:90
MpoTerms::reversedData::W
std::vector< std::vector< std::vector< std::vector< Biped< Symmetry, MatrixType > > > > > W
Definition: MpoTerms.h:404
MpoTransferVector::data
Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > data
Definition: VumpsMpoTransferMatrix.h:84
Multipede::setIdentity
void setIdentity(size_t Drows, size_t Dcols, size_t amax=1, size_t bmax=1)
Definition: Multipede.h:509
TransferVector::data
Biped< Symmetry, Matrix< Scalar_, Dynamic, Dynamic > > data
Definition: VumpsTransferMatrix.h:76
qarray
Definition: qarray.h:26