VMPS++
Loading...
Searching...
No Matches
DmrgPivotMatrix2.h
Go to the documentation of this file.
1#ifndef STRAWBERRY_DMRG_HEFF_STUFF_2SITE_WITH_Q
2#define STRAWBERRY_DMRG_HEFF_STUFF_2SITE_WITH_Q
3
5#include <vector>
7
8//include "DmrgTypedefs.h"
9//include "tensors/Biped.h"
10//include "tensors/Multipede.h"
11//include "Mps.h"
12#include "pivot/DmrgPivotVector.h"
13#include "pivot/DmrgPivotOverlap2.h"
14
15template<typename Symmetry, typename Scalar, typename MpoScalar=double>
16struct PivotMatrix2Terms
17{
18 Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > L;
19 Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > R;
20 vector<vector<vector<Biped<Symmetry,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE>> > > > W12;
21 vector<vector<vector<Biped<Symmetry,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE>> > > > W34;
22
23 vector<qarray<Symmetry::Nq> > qloc12;
24 vector<qarray<Symmetry::Nq> > qloc34;
25 vector<qarray<Symmetry::Nq> > qOp12;
26 vector<qarray<Symmetry::Nq> > qOp34;
27};
28
29template<typename Symmetry, typename Scalar, typename MpoScalar=double>
30struct PivotMatrix2
31{
32 PivotMatrix2(){};
33
34 PivotMatrix2 (const Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &L_input,
35 const Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > &R_input,
36 const vector<vector<vector<Biped<Symmetry,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE>> > > >& W12_input,
37 const vector<vector<vector<Biped<Symmetry,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE>> > > >& W34_input,
38 const vector<qarray<Symmetry::Nq> > &qloc12_input,
39 const vector<qarray<Symmetry::Nq> > &qloc34_input,
40 const vector<qarray<Symmetry::Nq> > &qOp12_input,
41 const vector<qarray<Symmetry::Nq> > &qOp34_input)
42 {
43 Terms.resize(1);
44 Terms[0].L = L_input;
45 Terms[0].R = R_input;
46 Terms[0].W12 = W12_input;
47 Terms[0].W34 = W34_input;
48 Terms[0].qloc12 = qloc12_input;
49 Terms[0].qloc34 = qloc34_input;
50 Terms[0].qOp12 = qOp12_input;
51 Terms[0].qOp34 = qOp34_input;
52 }
53
54// Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > L;
55// Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > R;
56// vector<vector<vector<Biped<Symmetry,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE>> > > > W12;
57// vector<vector<vector<Biped<Symmetry,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE>> > > > W34;
58//
59// vector<qarray<Symmetry::Nq> > qloc12;
60// vector<qarray<Symmetry::Nq> > qloc34;
61// vector<qarray<Symmetry::Nq> > qOp12;
62// vector<qarray<Symmetry::Nq> > qOp34;
63
64 vector<PivotMatrix2Terms<Symmetry,Scalar,MpoScalar> > Terms;
65
66 //---<pre-calculated, if Terms.size() == 1>---
67 vector<std::array<size_t,2> > qlhs;
68 vector<vector<std::array<size_t,12> > > qrhs;
69 vector<vector<Scalar> > factor_cgcs;
70 //--------------------------------
71
72 //---<stuff for excited states>---
73 vector<vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > > A0proj; // A0proj[n][s]
74 double Epenalty = 0;
75 //--------------------------------
76};
77
78template<typename Symmetry, typename Scalar, typename MpoScalar>
79void HxV (const PivotMatrix2<Symmetry,Scalar,MpoScalar> &H, const PivotVector<Symmetry,Scalar> &Vin, PivotVector<Symmetry,Scalar> &Vout)
80{
81 Vout = Vin;
82 OxV(H,Vin,Vout);
83}
84
85template<typename Symmetry, typename Scalar, typename MpoScalar>
86void OxV (const PivotMatrix2<Symmetry,Scalar,MpoScalar> &H, const PivotVector<Symmetry,Scalar> &Vin, PivotVector<Symmetry,Scalar> &Vout)
87{
88 for (size_t i=0; i<Vout.data.size(); ++i) {Vout.data[i].setZero();}
89 vector<PivotVector<Symmetry,Scalar> > Vt(H.Terms.size());
90 for (size_t t=0; t<H.Terms.size(); ++t) Vt[t] = Vout;
91
92// vector<std::array<size_t,2> > qlhs;
93// vector<vector<std::array<size_t,12> > > qrhs;
94// vector<vector<Scalar> > factor_cgcs;
95// precalc_blockStructure<Symmetry,Scalar,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE> >
96// (H.L, Vin.data, H.W12, H.W34, Vin.data, H.R, H.qloc12, H.qloc34, H.qOp12, H.qOp34, qlhs, qrhs, factor_cgcs);
97//
98// #ifdef DMRG_PIVOT2_PARALLELIZE
99// #pragma omp parallel for schedule(dynamic)
100// #endif
101// for (size_t q=0; q<qlhs.size(); ++q)
102// {
103// size_t s1s3 = qlhs[q][0];//H.qlhs[q][0];
104// size_t q13 = qlhs[q][1];//H.qlhs[q][1];
105//
106// for (size_t p=0; p<qrhs[q].size(); ++p)
107// {
108// size_t s2s4 = qrhs[q][p][0];//H.qrhs[q][p][0];
109// size_t q24 = qrhs[q][p][1];//H.qrhs[q][p][1];
110// size_t qL = qrhs[q][p][2];//H.qrhs[q][p][2];
111// size_t qR = qrhs[q][p][3];//H.qrhs[q][p][3];
112// size_t s1 = qrhs[q][p][4];//H.qrhs[q][p][4];
113// size_t s2 = qrhs[q][p][5];//H.qrhs[q][p][5];
114// size_t k12 = qrhs[q][p][6];//H.qrhs[q][p][6];
115// size_t qW12 = qrhs[q][p][7];//H.qrhs[q][p][7];
116// size_t s3 = qrhs[q][p][8];//H.qrhs[q][p][8];
117// size_t s4 = qrhs[q][p][9];//H.qrhs[q][p][9];
118// size_t k34 = qrhs[q][p][10];//H.qrhs[q][p][10];
119// size_t qW34 = qrhs[q][p][11];//H.qrhs[q][p][11];
120//
121// for (int r12=0; r12<H.W12[s1][s2][k12].block[qW12].outerSize(); ++r12)
122// for (typename SparseMatrix<MpoScalar>::InnerIterator iW12(H.W12[s1][s2][k12].block[qW12],r12); iW12; ++iW12)
123// for (int r34=0; r34<H.W34[s3][s4][k34].block[qW34].outerSize(); ++r34)
124// for (typename SparseMatrix<MpoScalar>::InnerIterator iW34(H.W34[s3][s4][k34].block[qW34],r34); iW34; ++iW34)
125// {
126// if (H.L.block[qL][iW12.row()][0].size() != 0 and
127// H.R.block[qR][iW34.col()][0].size() != 0 and
128// Vin.data[s2s4].block[q24].size() != 0 and
129// iW12.col() == iW34.row())
130// {
137//
138// if (Vout.data[s1s3].block[q13].rows() != H.L.block[qL][iW12.row()][0].rows() or
139// Vout.data[s1s3].block[q13].cols() != H.R.block[qR][iW34.col()][0].cols())
140// {
141// Vout.data[s1s3].block[q13].noalias() = factor_cgcs[q][p] * iW12.value() * iW34.value() *
142// (H.L.block[qL][iW12.row()][0] *
143// Vin.data[s2s4].block[q24] *
144// H.R.block[qR][iW34.col()][0]);
145// }
146// else
147// {
148// Vout.data[s1s3].block[q13].noalias() += factor_cgcs[q][p] * iW12.value() * iW34.value() *
149// (H.L.block[qL][iW12.row()][0] *
150// Vin.data[s2s4].block[q24] *
151// H.R.block[qR][iW34.col()][0]);
152// }
153// }
154// }
155// }
156// }
157
158 #ifdef DMRG_PARALLELIZE_TERMS
159 #pragma omp parallel for schedule(dynamic)
160 #endif
161 for (size_t t=0; t<H.Terms.size(); ++t)
162 {
163 vector<std::array<size_t,2> > qlhs;
164 vector<vector<std::array<size_t,12> > > qrhs;
165 vector<vector<Scalar> > factor_cgcs;
166
167 if (H.Terms.size() == 1)
168 {
169 qlhs = H.qlhs;
170 qrhs = H.qrhs;
171 factor_cgcs = H.factor_cgcs;
172 }
173 else
174 {
175 precalc_blockStructure<Symmetry,Scalar,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE> >
176 (H.Terms[t].L, Vin.data, H.Terms[t].W12, H.Terms[t].W34, Vin.data, H.Terms[t].R, H.Terms[t].qloc12, H.Terms[t].qloc34, H.Terms[t].qOp12, H.Terms[t].qOp34, qlhs, qrhs, factor_cgcs);
177 }
178
179 #ifdef DMRG_PIVOT2_PARALLELIZE
180 #pragma omp parallel for schedule(dynamic)
181 #endif
182 for (size_t q=0; q<qlhs.size(); ++q)
183 {
184 size_t s1s3 = qlhs[q][0];
185 size_t q13 = qlhs[q][1];
186
187 for (size_t p=0; p<qrhs[q].size(); ++p)
188 {
189 size_t s2s4 = qrhs[q][p][0];
190 size_t q24 = qrhs[q][p][1];
191 size_t qL = qrhs[q][p][2];
192 size_t qR = qrhs[q][p][3];
193 size_t s1 = qrhs[q][p][4];
194 size_t s2 = qrhs[q][p][5];
195 size_t k12 = qrhs[q][p][6];
196 size_t qW12 = qrhs[q][p][7];
197 size_t s3 = qrhs[q][p][8];
198 size_t s4 = qrhs[q][p][9];
199 size_t k34 = qrhs[q][p][10];
200 size_t qW34 = qrhs[q][p][11];
201
202 for (int r12=0; r12<H.Terms[t].W12[s1][s2][k12].block[qW12].outerSize(); ++r12)
203 for (typename SparseMatrix<MpoScalar>::InnerIterator iW12(H.Terms[t].W12[s1][s2][k12].block[qW12],r12); iW12; ++iW12)
204 for (int r34=0; r34<H.Terms[t].W34[s3][s4][k34].block[qW34].outerSize(); ++r34)
205 for (typename SparseMatrix<MpoScalar>::InnerIterator iW34(H.Terms[t].W34[s3][s4][k34].block[qW34],r34); iW34; ++iW34)
206 {
207 if (H.Terms[t].L.block[qL][iW12.row()][0].size() != 0 and
208 H.Terms[t].R.block[qR][iW34.col()][0].size() != 0 and
209 Vin.data[s2s4].block[q24].size() != 0 and
210 iW12.col() == iW34.row())
211 {
212 if (Vt[t].data[s1s3].block[q13].rows() != H.Terms[t].L.block[qL][iW12.row()][0].rows() or
213 Vt[t].data[s1s3].block[q13].cols() != H.Terms[t].R.block[qR][iW34.col()][0].cols())
214 {
215 Vt[t].data[s1s3].block[q13].noalias() = factor_cgcs[q][p] * iW12.value() * iW34.value() *
216 (H.Terms[t].L.block[qL][iW12.row()][0] *
217 Vin.data[s2s4].block[q24] *
218 H.Terms[t].R.block[qR][iW34.col()][0]);
219 }
220 else
221 {
222 Vt[t].data[s1s3].block[q13].noalias() += factor_cgcs[q][p] * iW12.value() * iW34.value() *
223 (H.Terms[t].L.block[qL][iW12.row()][0] *
224 Vin.data[s2s4].block[q24] *
225 H.Terms[t].R.block[qR][iW34.col()][0]);
226 }
227 }
228 }
229 }
230 }
231 }
232
233 for (size_t s=0; s<Vout.size(); ++s)
234 {
235 Vout[s] = Vt[0][s];
236 }
237
238 #ifdef DMRG_PARALLELIZE_TERMS
239 #pragma omp parallel for
240 #endif
241 for (size_t s=0; s<Vout.size(); ++s)
242 for (size_t t=1; t<H.Terms.size(); ++t)
243 {
244 Vout[s].addScale(1.,Vt[t][s]);
245 }
246
247 if (H.Terms.size() > 0) for (size_t s=0; s<Vout.size(); ++s) Vout[s] = Vout[s].cleaned();
248
249 // project out unwanted states (e.g. to get lower spectrum)
250 for (size_t n=0; n<H.A0proj.size(); ++n)
251 {
252 Scalar overlap = 0;
253 for (size_t s=0; s<H.A0proj[n].size(); ++s)
254 {
255 // Note: Adjoint needed because we need <E0|Psi>, not <Psi|E0>
256 overlap += H.A0proj[n][s].adjoint().contract(Vin.data[s]).trace();
257 }
258// cout << "overlap=" << overlap << endl;
259
260 for (size_t s=0; s<H.A0proj[n].size(); ++s)
261 for (size_t q=0; q<H.A0proj[n][s].dim; ++q)
262 {
263 qarray2<Symmetry::Nq> cmp = {H.A0proj[n][s].in[q], H.A0proj[n][s].out[q]};
264 auto qA = Vout.data[s].dict.find(cmp);
265// assert(qA != Vout.data[s].dict.end() and "Error in HxV(PivotMatrix1,PivotVector): projected block not found!");
266
267 if (qA != Vout.data[s].dict.end() and H.A0proj[n][s].block[q].size() != 0)
268 {
269 Vout.data[s].block[qA->second] += H.Epenalty * overlap * H.A0proj[n][s].block[q];
270 }
271 }
272 }
273}
274
275//template<typename Symmetry, typename Scalar, typename MpoScalar>
276//void OxV (const PivotMatrix2<Symmetry,Scalar,MpoScalar> &H, const PivotVector<Symmetry,Scalar> &Vin, PivotVector<Symmetry,Scalar> &Vout)
277//{
278// auto tensor_basis = Symmetry::tensorProd(H.qloc12, H.qloc34);
279// for (size_t i=0; i<Vout.data.size(); ++i)
280// {
281// Vout.data[i].setZero();
282// }
283//
284// Stopwatch Wtot;
285// double ttot=0;
286// double tmult=0;
287// double tcgc=0;
288// double tmatch=0;
289// int N_cgc_calc=0;
290//
291// for (size_t s1=0; s1<H.qloc12.size(); ++s1)
292// for (size_t s2=0; s2<H.qloc12.size(); ++s2)
293// for (size_t k12=0; k12<H.qOp12.size(); ++k12)
294// {
295// if (!Symmetry::validate(qarray3<Symmetry::Nq>{H.qloc12[s2], H.qOp12[k12], H.qloc12[s1]})) {continue;}
296//
297// for (size_t s3=0; s3<H.qloc34.size(); ++s3)
298// for (size_t s4=0; s4<H.qloc34.size(); ++s4)
299// for (size_t k34=0; k34<H.qOp34.size(); ++k34)
300// {
301// if (!Symmetry::validate(qarray3<Symmetry::Nq>{H.qloc34[s4], H.qOp34[k34], H.qloc34[s3]})) {continue;}
302//
303// auto qOps = Symmetry::reduceSilent(H.qOp12[k12], H.qOp34[k34]);
304//
305// for (const auto &qOp:qOps)
306// {
307// auto qmerges13 = Symmetry::reduceSilent(H.qloc12[s1], H.qloc34[s3]);
308// auto qmerges24 = Symmetry::reduceSilent(H.qloc12[s2], H.qloc34[s4]);
309//
310// for (const auto &qmerge13:qmerges13)
311// for (const auto &qmerge24:qmerges24)
312// {
313// auto qtensor13 = make_tuple(H.qloc12[s1], s1, H.qloc34[s3], s3, qmerge13);
314// auto s1s3 = distance(tensor_basis.begin(), find(tensor_basis.begin(), tensor_basis.end(), qtensor13));
315//
316// auto qtensor24 = make_tuple(H.qloc12[s2], s2, H.qloc34[s4], s4, qmerge24);
317// auto s2s4 = distance(tensor_basis.begin(), find(tensor_basis.begin(), tensor_basis.end(), qtensor24));
318//
319// // tensor product of the MPO operators in the physical space
320// Stopwatch<> Wcgc9;
321// Scalar factor_cgc9 = (Symmetry::NON_ABELIAN)?
322// Symmetry::coeff_buildR(H.qloc12[s2], H.qloc34[s4], qmerge24,
323// H.qOp12[k12], H.qOp34[k34], qOp,
324// H.qloc12[s1], H.qloc34[s3], qmerge13)
325// :1.;
326// tcgc += Wcgc9.time();
327// ++N_cgc_calc;
328// if (abs(factor_cgc9) < abs(mynumeric_limits<Scalar>::epsilon())) {continue;}
329//
330// for (size_t qL=0; qL<H.L.dim; ++qL)
331// {
332// vector<tuple<qarray3<Symmetry::Nq>,qarray<Symmetry::Nq>,size_t,size_t> > ixs;
333// Stopwatch Wmatch;
334// bool FOUND_MATCH = AAWWAA(H.L.in(qL), H.L.out(qL), H.L.mid(qL),
335// k12, H.qOp12, k34, H.qOp34,
336// s1s3, qmerge13, s2s4, qmerge24,
337// Vout.data, Vin.data, ixs);
338// tmatch += Wmatch.time();
339//
340// if (FOUND_MATCH)
341// {
342// for (const auto &ix:ixs)
343// {
344// auto qR = H.R.dict.find(get<0>(ix));
345// auto qW = get<1>(ix);
346// size_t qA13 = get<2>(ix);
347// size_t qA24 = get<3>(ix);
348//
349// // multiplication of Op12, Op34 in the auxiliary space
350// Stopwatch<> Wcgc6;
351// Scalar factor_cgc6 = (Symmetry::NON_ABELIAN)?
352// Symmetry::coeff_Apair(H.L.mid(qL), H.qOp12[k12], qW,
353// H.qOp34[k34], get<0>(ix)[2], qOp)
354// :1.;
355// tcgc += Wcgc6.time();
356// ++N_cgc_calc;
357// if (abs(factor_cgc6) < abs(mynumeric_limits<Scalar>::epsilon())) {continue;}
358//
359// if (qR != H.R.dict.end())
360// {
361// // standard coefficient for H*Psi with environments
362// Stopwatch<> WcgcHPsi;
363// Scalar factor_cgcHPsi = (Symmetry::NON_ABELIAN)?
364// Symmetry::coeff_HPsi(Vin.data[s2s4].out[qA24], qmerge24, Vin.data[s2s4].in[qA24],
365// H.R.mid(qR->second), qOp, H.L.mid(qL),
366// Vout.data[s1s3].out[qA13], qmerge13, Vout.data[s1s3].in[qA13])
367// :1.;
368// ++N_cgc_calc;
369// tcgc += WcgcHPsi.time();
370//
371// for (int r12=0; r12<H.W12[s1][s2][k12].outerSize(); ++r12)
372// for (typename SparseMatrix<MpoScalar>::InnerIterator iW12(H.W12[s1][s2][k12],r12); iW12; ++iW12)
373// for (int r34=0; r34<H.W34[s3][s4][k34].outerSize(); ++r34)
374// for (typename SparseMatrix<MpoScalar>::InnerIterator iW34(H.W34[s3][s4][k34],r34); iW34; ++iW34)
375// {
376// Matrix<Scalar,Dynamic,Dynamic> Mtmp;
377// auto Wfactor = iW12.value() * iW34.value() * factor_cgc6 * factor_cgc9 * factor_cgcHPsi;
378//
379// if (H.L.block[qL][iW12.row()][0].size() != 0 and
380// H.R.block[qR->second][iW34.col()][0].size() !=0 and
381// iW12.col() == iW34.row())
382// {
383// Stopwatch Wmult;
384// optimal_multiply(Wfactor,
385// H.L.block[qL][iW12.row()][0],
386// Vin.data[s2s4].block[qA24],
387// H.R.block[qR->second][iW34.col()][0],
388// Mtmp);
389// tmult += Wmult.time();
390// }
391//
392// if (Mtmp.size() != 0)
393// {
394// if (Vout.data[s1s3].block[qA13].rows() == Mtmp.rows() and
395// Vout.data[s1s3].block[qA13].cols() == Mtmp.cols())
396// {
397// Vout.data[s1s3].block[qA13] += Mtmp;
398// }
399// else
400// {
401// Vout.data[s1s3].block[qA13] = Mtmp;
402// }
403// }
404// }
405// }
406// }
407// }
408// }
409// }
410// }
411// }
412// }
413//
414// ttot = Wtot.time();
415//
417//}
418
419template<typename Symmetry, typename Scalar, typename MpoScalar>
420void HxV (const PivotMatrix2<Symmetry,Scalar,MpoScalar> &H, PivotVector<Symmetry,Scalar> &Vinout)
421{
422 PivotVector<Symmetry,Scalar> Vtmp;
423 HxV(H,Vinout,Vtmp);
424 Vinout = Vtmp;
425}
426
427template<typename Symmetry, typename Scalar, typename MpoScalar>
428inline size_t dim (const PivotMatrix2<Symmetry,Scalar,MpoScalar> &H)
429{
430 return 0;
431}
432
433// How to calculate the Frobenius norm of this?
434template<typename Symmetry, typename Scalar, typename MpoScalar>
435inline double norm (const PivotMatrix2<Symmetry,Scalar,MpoScalar> &H)
436{
437 return 0;
438}
439
440//template<typename Symmetry, typename Scalar, typename MpoScalar>
441//void HxV (const PivotMatrix1<Symmetry,Scalar,MpoScalar> &H1,
442// const PivotMatrix1<Symmetry,Scalar,MpoScalar> &H2,
443// const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Aket1,
444// const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Aket2,
445// const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Abra1,
446// const vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > &Abra2,
447// const vector<qarray<Symmetry::Nq> > &qloc1, const vector<qarray<Symmetry::Nq> > &qloc2,
448// vector<vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > > &Apair)
449//{
450// Apair.resize(qloc1.size());
451// for (size_t s1=0; s1<qloc1.size(); ++s1)
452// {
453// Apair[s1].resize(qloc2.size());
454// }
455//
456// for (size_t s1=0; s1<qloc1.size(); ++s1)
457// for (size_t s2=0; s2<qloc1.size(); ++s2)
458// for (size_t qL=0; qL<H1.L.dim; ++qL)
459// {
460// tuple<qarray3<Symmetry::Nq>,size_t,size_t> ix12;
461// bool FOUND_MATCH12 = AWA(H1.L.in(qL), H1.L.out(qL), H1.L.mid(qL), s1, s2, qloc1, Abra1, Aket1, ix12);
462//
463// if (FOUND_MATCH12)
464// {
465// qarray3<Symmetry::Nq> quple12 = get<0>(ix12);
466// swap(quple12[0], quple12[1]);
467// size_t qA12 = get<2>(ix12);
468//
469// for (size_t s3=0; s3<qloc2.size(); ++s3)
470// for (size_t s4=0; s4<qloc2.size(); ++s4)
471// {
472// tuple<qarray3<Symmetry::Nq>,size_t,size_t> ix34;
473// bool FOUND_MATCH34 = AWA(quple12[0], quple12[1], quple12[2], s3, s4, qloc2, Abra2, Aket2, ix34);
474//
475// if (FOUND_MATCH34)
476// {
477// qarray3<Symmetry::Nq> quple34 = get<0>(ix34);
478// size_t qA34 = get<2>(ix34);
479// auto qR = H2.R.dict.find(quple34);
480//
481// if (qR != H2.R.dict.end())
482// {
483// if (H1.L.mid(qL) + qloc1[s1] - qloc1[s2] ==
484// H2.R.mid(qR->second) - qloc2[s3] + qloc2[s4])
485// {
486// for (int k12=0; k12<H1.W[s1][s2][0].outerSize(); ++k12)
487// for (typename SparseMatrix<MpoScalar>::InnerIterator iW12(H1.W[s1][s2][0],k12); iW12; ++iW12)
488// for (int k34=0; k34<H2.W[s3][s4][0].outerSize(); ++k34)
489// for (typename SparseMatrix<MpoScalar>::InnerIterator iW34(H2.W[s3][s4][0],k34); iW34; ++iW34)
490// {
491// Matrix<Scalar,Dynamic,Dynamic> Mtmp;
492// MpoScalar Wfactor = iW12.value() * iW34.value();
493//
494// if (H1.L.block[qL][iW12.row()][0].rows() != 0 and
495// H2.R.block[qR->second][iW34.col()][0].rows() !=0 and
496// iW12.col() == iW34.row())
497// {
503// optimal_multiply(Wfactor,
504// H1.L.block[qL][iW12.row()][0],
505// Aket1[s2].block[qA12],
506// Aket2[s4].block[qA34],
507// H2.R.block[qR->second][iW34.col()][0],
508// Mtmp);
509// }
510//
511// if (Mtmp.rows() != 0)
512// {
513// qarray2<Symmetry::Nq> qupleApair = {H1.L.in(qL), H2.R.out(qR->second)};
514// auto qApair = Apair[s1][s3].dict.find(qupleApair);
515//
516// if (qApair != Apair[s1][s3].dict.end())
517// {
518// Apair[s1][s3].block[qApair->second] += Mtmp;
519// }
520// else
521// {
522// Apair[s1][s3].push_back(qupleApair, Mtmp);
523// }
524// }
525// }
526// }
527// }
528// }
529// }
530// }
531// }
532//}
533
534#endif
dim
size_t dim(const PivotMatrix2< Symmetry, Scalar, MpoScalar > &H)
Definition: DmrgPivotMatrix2.h:428
norm
double norm(const PivotMatrix2< Symmetry, Scalar, MpoScalar > &H)
Definition: DmrgPivotMatrix2.h:435
OxV
void OxV(const PivotMatrix2< Symmetry, Scalar, MpoScalar > &H, const PivotVector< Symmetry, Scalar > &Vin, PivotVector< Symmetry, Scalar > &Vout)
Definition: DmrgPivotMatrix2.h:86
HxV
void HxV(const PivotMatrix2< Symmetry, Scalar, MpoScalar > &H, const PivotVector< Symmetry, Scalar > &Vin, PivotVector< Symmetry, Scalar > &Vout)
Definition: DmrgPivotMatrix2.h:79
qarray2
std::array< qarray< Nq >, 2 > qarray2
Definition: qarray.h:51
Biped
Definition: Biped.h:64
PivotMatrix2Terms::W12
vector< vector< vector< Biped< Symmetry, Eigen::SparseMatrix< MpoScalar, Eigen::ColMajor, EIGEN_DEFAULT_SPARSE_INDEX_TYPE > > > > > W12
Definition: DmrgPivotMatrix2.h:20
PivotMatrix2Terms::R
Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > R
Definition: DmrgPivotMatrix2.h:19
PivotMatrix2Terms::L
Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > L
Definition: DmrgPivotMatrix2.h:18
PivotMatrix2Terms::W34
vector< vector< vector< Biped< Symmetry, Eigen::SparseMatrix< MpoScalar, Eigen::ColMajor, EIGEN_DEFAULT_SPARSE_INDEX_TYPE > > > > > W34
Definition: DmrgPivotMatrix2.h:21
PivotMatrix2Terms::qloc12
vector< qarray< Symmetry::Nq > > qloc12
Definition: DmrgPivotMatrix2.h:23
PivotMatrix2Terms::qOp34
vector< qarray< Symmetry::Nq > > qOp34
Definition: DmrgPivotMatrix2.h:26
PivotMatrix2Terms::qloc34
vector< qarray< Symmetry::Nq > > qloc34
Definition: DmrgPivotMatrix2.h:24
PivotMatrix2Terms::qOp12
vector< qarray< Symmetry::Nq > > qOp12
Definition: DmrgPivotMatrix2.h:25
PivotMatrix2::A0proj
vector< vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > > A0proj
Definition: DmrgPivotMatrix2.h:73
PivotMatrix2::qrhs
vector< vector< std::array< size_t, 12 > > > qrhs
Definition: DmrgPivotMatrix2.h:68
PivotMatrix2::qlhs
vector< std::array< size_t, 2 > > qlhs
Definition: DmrgPivotMatrix2.h:67
PivotMatrix2::factor_cgcs
vector< vector< Scalar > > factor_cgcs
Definition: DmrgPivotMatrix2.h:69
PivotMatrix2::PivotMatrix2
PivotMatrix2(const Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > &L_input, const Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > &R_input, const vector< vector< vector< Biped< Symmetry, Eigen::SparseMatrix< MpoScalar, Eigen::ColMajor, EIGEN_DEFAULT_SPARSE_INDEX_TYPE > > > > > &W12_input, const vector< vector< vector< Biped< Symmetry, Eigen::SparseMatrix< MpoScalar, Eigen::ColMajor, EIGEN_DEFAULT_SPARSE_INDEX_TYPE > > > > > &W34_input, const vector< qarray< Symmetry::Nq > > &qloc12_input, const vector< qarray< Symmetry::Nq > > &qloc34_input, const vector< qarray< Symmetry::Nq > > &qOp12_input, const vector< qarray< Symmetry::Nq > > &qOp34_input)
Definition: DmrgPivotMatrix2.h:34
PivotMatrix2::Terms
vector< PivotMatrix2Terms< Symmetry, Scalar, MpoScalar > > Terms
Definition: DmrgPivotMatrix2.h:64
PivotVector::size
size_t size() const
Definition: DmrgPivotVector.h:74
PivotVector::data
vector< Biped< Symmetry, Matrix< Scalar_, Dynamic, Dynamic > > > data
Definition: DmrgPivotVector.h:116
qarray
Definition: qarray.h:26