VMPS++
Loading...
Searching...
No Matches
DmrgPivotMatrix1.h
Go to the documentation of this file.
1#ifndef DMRGPIVOTMATRIX1
2#define DMRGPIVOTMATRIX1
3
4#include "pivot/DmrgPivotVector.h"
5#include "tensors/DmrgIndexGymnastics.h"
6//include "DmrgTypedefs.h"
7//include "tensors/Biped.h"
8//include "tensors/Multipede.h"
9
10//-----------<definitions>-----------
11template<typename Symmetry, typename Scalar, typename MpoScalar=double>
12struct PivotMatrix1Terms
13{
14 Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > L;
15 Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > R;
16 vector<vector<vector<Biped<Symmetry,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE> > > > > W;
17
18 vector<qarray<Symmetry::Nq> > qloc;
19 vector<qarray<Symmetry::Nq> > qOp;
20
21// void save_L (string filename)
22// {
23// lout << termcolor::green << "Saving L to: " << filename << termcolor::reset << std::endl;
24// L.save(filename);
25// }
26//
27// void save_R (string filename)
28// {
29// lout << termcolor::green << "Saving R to: " << filename << termcolor::reset << std::endl;
30// R.save(filename);
31// }
32//
33// void load_L (string filename)
34// {
35// lout << termcolor::green << "Loading L from: " << filename << termcolor::reset << std::endl;
36// L.load(filename);
37// }
38//
39// void load_R (string filename)
40// {
41// lout << termcolor::green << "Loading R from: " << filename << termcolor::reset << std::endl;
42// R.load(filename);
43// }
44};
45
46template<typename Symmetry, typename Scalar, typename MpoScalar=double>
47struct PivotMatrix1
48{
49// Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > L;
50// Tripod<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > R;
51// vector<vector<vector<Biped<Symmetry,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE> > > > > W;
52//
53// vector<qarray<Symmetry::Nq> > qloc;
54// vector<qarray<Symmetry::Nq> > qOp;
55
56 vector<PivotMatrix1Terms<Symmetry,Scalar,MpoScalar> > Terms;
57
58 //---<pre-calculated, if Terms.size() == 1>---
59 vector<std::array<size_t,2> > qlhs;
60 vector<vector<std::array<size_t,6> > > qrhs;
61 vector<vector<Scalar> > factor_cgcs;
62 //--------------------------------
63
64 //---<stuff for excited states>---
65 vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > PL; // PL[n]
66 vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > PR; // PL[n]
67 vector<vector<Biped<Symmetry,Matrix<Scalar,Dynamic,Dynamic> > > > A0proj; // A0proj[n][s]
68 double Epenalty = 0;
69 //--------------------------------
70
71 template<typename OtherScalar>
72 PivotMatrix1<Symmetry,OtherScalar,OtherScalar> cast() const
73 {
74 PivotMatrix1<Symmetry,OtherScalar,OtherScalar> Pout;
75
76 Pout.Terms.resize(Terms.size());
77
78 for (size_t t=0; t<Terms.size(); ++t)
79 {
80 Pout.Terms[t].L = Terms[t].L.template cast<Matrix<OtherScalar,Dynamic,Dynamic> >();
81 Pout.Terms[t].R = Terms[t].R.template cast<Matrix<OtherScalar,Dynamic,Dynamic> >();
82// Pout.W = W;
83
84 Pout.Terms[t].qloc = Terms[t].qloc;
85 Pout.Terms[t].qOp = Terms[t].qOp;
86 }
87
88 Pout.qlhs = qlhs;
89 Pout.qrhs = qrhs;
90 Pout.factor_cgcs.resize(factor_cgcs.size());
91 for (int i=0; i<factor_cgcs.size(); ++i) Pout.factor_cgcs[i].resize(factor_cgcs[i].size());
92
93 for (int i=0; i<factor_cgcs.size(); ++i)
94 for (int j=0; j<factor_cgcs[i].size(); ++j)
95 {
96 Pout.factor_cgcs[i][j] = factor_cgcs[i][j];
97 }
98
99// Pout.qloc = qloc;
100// Pout.qOp = qOp;
101
102 return Pout;
103 }
104
105 double memory (MEMUNIT memunit) const
106 {
107 double res = 0.;
108 for (size_t t=0; t<Terms.size(); ++t)
109 {
110 res += Terms[t].L.memory(memunit);
111 res += Terms[t].R.memory(memunit);
112 for (size_t s1=0; s1<Terms[t].W.size(); ++s1)
113 for (size_t s2=0; s2<Terms[t].W[s1].size(); ++s2)
114 for (size_t k=0; k<Terms[t].W[s1][s2].size(); ++k)
115 {
116 res += Terms[t].W[s1][s2][k].memory(memunit);
117 }
118 }
119 return res;
120 }
121};
122
123template<typename Symmetry, typename Scalar, typename MpoScalar>
124void HxV (const PivotMatrix1<Symmetry,Scalar,MpoScalar> &H, const PivotVector<Symmetry,Scalar> &Vin, PivotVector<Symmetry,Scalar> &Vout)
125{
126// Vout.outerResize(Vin);
127 Vout = Vin;
128 OxV(H,Vin,Vout);
129}
130
131template<typename Symmetry, typename Scalar, typename MpoScalar>
132void OxV (const PivotMatrix1<Symmetry,Scalar,MpoScalar> &H, const PivotVector<Symmetry,Scalar> &Vin, PivotVector<Symmetry,Scalar> &Vout)
133{
134 for (size_t s=0; s<Vout.data.size(); ++s) {Vout.data[s].setZero();}
135 vector<PivotVector<Symmetry,Scalar> > Vt(H.Terms.size());
136 for (size_t t=0; t<H.Terms.size(); ++t) Vt[t] = Vout;
137
138 vector<std::array<size_t,2> > qlhs;
139 vector<vector<std::array<size_t,6> > > qrhs;
140 vector<vector<Scalar> > factor_cgcs;
141
142// precalc_blockStructure<Symmetry,Scalar,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE> >
143// (H.L, Vin.data, H.W, Vin.data, H.R, H.qloc, H.qOp, qlhs, qrhs, factor_cgcs);
144// //cout << "qlhs.size()=" << qlhs.size() << endl;
145//
146// #ifdef DMRG_PIVOT1_PARALLELIZE
147// #pragma omp parallel for schedule(dynamic)
148// #endif
149// for (size_t q=0; q<qlhs.size(); ++q)
150// {
151// size_t s1 = qlhs[q][0];//H.qlhs[q][0];
152// size_t q1 = qlhs[q][1];//H.qlhs[q][1];
153//
154// for (size_t p=0; p<qrhs[q].size(); ++p)
155// {
156// size_t s2 = qrhs[q][p][0];//H.qrhs[q][p][0];
157// size_t q2 = qrhs[q][p][1];//H.qrhs[q][p][1];
158// size_t qL = qrhs[q][p][2];//H.qrhs[q][p][2];
159// size_t qR = qrhs[q][p][3];//H.qrhs[q][p][3];
160// size_t k = qrhs[q][p][4];//H.qrhs[q][p][4];
161// size_t qW = qrhs[q][p][5];//H.qrhs[q][p][5];
162//
163// for (int r=0; r<H.W[s1][s2][k].block[qW].outerSize(); ++r)
164// for (typename Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE>::InnerIterator iW(H.W[s1][s2][k].block[qW],r); iW; ++iW)
165// {
166// if (H.L.block[qL][iW.row()][0].size() != 0 and
167// H.R.block[qR][iW.col()][0].size() != 0 and
168// Vin.data[s2].block[q2].size() !=0)
169// {
174//
175// if (Vout.data[s1].block[q1].rows() != H.L.block[qL][iW.row()][0].rows() or
176// Vout.data[s1].block[q1].cols() != H.R.block[qR][iW.col()][0].cols())
177// {
178// Vout.data[s1].block[q1].noalias() = factor_cgcs[q][p] * iW.value() *
179// (H.L.block[qL][iW.row()][0] *
180// Vin.data[s2].block[q2] *
181// H.R.block[qR][iW.col()][0]);
182// }
183// else
184// {
185// Vout.data[s1].block[q1].noalias() += factor_cgcs[q][p] * iW.value() *
186// (H.L.block[qL][iW.row()][0] *
187// Vin.data[s2].block[q2] *
188// H.R.block[qR][iW.col()][0]);
189// }
190// }
191// }
192// }
193// }
194
195 //Stopwatch<> Timer1;
196
197 #ifdef DMRG_PARALLELIZE_TERMS
198 #pragma omp parallel for schedule(dynamic)
199 #endif
200 for (size_t t=0; t<H.Terms.size(); ++t)
201 {
202 vector<std::array<size_t,2> > qlhs;
203 vector<vector<std::array<size_t,6> > > qrhs;
204 vector<vector<Scalar> > factor_cgcs;
205
206 if (H.Terms.size() == 1)
207 {
208 qlhs = H.qlhs;
209 qrhs = H.qrhs;
210 factor_cgcs = H.factor_cgcs;
211 }
212 else
213 {
214 precalc_blockStructure<Symmetry,Scalar,Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE> >
215 (H.Terms[t].L, Vin.data, H.Terms[t].W, Vin.data, H.Terms[t].R, H.Terms[t].qloc, H.Terms[t].qOp, qlhs, qrhs, factor_cgcs);
216 //cout << "t=" << t << ", qlhs.size()=" << qlhs.size() << endl;
217 }
218
219 #ifdef DMRG_PIVOT1_PARALLELIZE
220 #pragma omp parallel for schedule(dynamic)
221 #endif
222 for (size_t q=0; q<qlhs.size(); ++q)
223 {
224 size_t s1 = qlhs[q][0];//H.Terms[t].qlhs[q][0];
225 size_t q1 = qlhs[q][1];//H.Terms[t].qlhs[q][1];
226
227 for (size_t p=0; p<qrhs[q].size(); ++p)
228 {
229 size_t s2 = qrhs[q][p][0];//H.Terms[t].qrhs[q][p][0];
230 size_t q2 = qrhs[q][p][1];//H.Terms[t].qrhs[q][p][1];
231 size_t qL = qrhs[q][p][2];//H.Terms[t].qrhs[q][p][2];
232 size_t qR = qrhs[q][p][3];//H.Terms[t].qrhs[q][p][3];
233 size_t k = qrhs[q][p][4];//H.Terms[t].qrhs[q][p][4];
234 size_t qW = qrhs[q][p][5];//H.Terms[t].qrhs[q][p][5];
235
236 for (int r=0; r<H.Terms[t].W[s1][s2][k].block[qW].outerSize(); ++r)
237 for (typename Eigen::SparseMatrix<MpoScalar,Eigen::ColMajor,EIGEN_DEFAULT_SPARSE_INDEX_TYPE>::InnerIterator iW(H.Terms[t].W[s1][s2][k].block[qW],r); iW; ++iW)
238 {
239// print_size(H.Terms[t].L.block[qL][iW.row()][0],"H.Terms[t].L.block[qL][iW.row()][0]");
240// print_size(Vin.data[s2].block[q2], "Vin.data[s2].block[q2]");
241// print_size(H.Terms[t].R.block[qR][iW.col()][0], "H.Terms[t].R.block[qR][iW.col()][0]");
242// cout << endl;
243
244 if (H.Terms[t].L.block[qL][iW.row()][0].size() != 0 and
245 H.Terms[t].R.block[qR][iW.col()][0].size() != 0 and
246 Vin.data[s2].block[q2].size() !=0)
247 {
248 if (Vt[t].data[s1].block[q1].rows() != H.Terms[t].L.block[qL][iW.row()][0].rows() or
249 Vt[t].data[s1].block[q1].cols() != H.Terms[t].R.block[qR][iW.col()][0].cols())
250 {
251 Vt[t].data[s1].block[q1].noalias() = factor_cgcs[q][p] * iW.value() *
252 (H.Terms[t].L.block[qL][iW.row()][0] *
253 Vin.data[s2].block[q2] *
254 H.Terms[t].R.block[qR][iW.col()][0]);
255 }
256 else
257 {
258 Vt[t].data[s1].block[q1].noalias() += factor_cgcs[q][p] * iW.value() *
259 (H.Terms[t].L.block[qL][iW.row()][0] *
260 Vin.data[s2].block[q2] *
261 H.Terms[t].R.block[qR][iW.col()][0]);
262 }
263 }
264 }
265 }
266 }
267 }
268
269 //double t1 = Timer1.time();
270 //cout << "multiplication: " << t1 << endl;
271 //Stopwatch<> Timer2;
272
273 for (size_t s=0; s<Vout.size(); ++s)
274 {
275 Vout[s] = Vt[0][s];
276 }
277
278 #ifdef DMRG_PARALLELIZE_TERMS
279 #pragma omp parallel for
280 #endif
281 for (size_t s=0; s<Vout.size(); ++s)
282 for (size_t t=1; t<H.Terms.size(); ++t)
283 {
284 Vout[s].addScale(1.,Vt[t][s]);
285 }
286
287 if (H.Terms.size() > 0) for (size_t s=0; s<Vout.size(); ++s) Vout[s] = Vout[s].cleaned();
288
289 //double t2 = Timer2.time();
290 //cout << "sum: " << t2 << endl;
291 //cout << "t2/t1=" << t2/t1 << endl;
292
293// for (size_t s=0; s<Vin.data.size(); ++s)
294// for (size_t q=0; q<Vin.data[s].dim; ++q)
295// {
296// cout << "Vin inout=" << Vin.data[s].in[q] << ", " << Vin.data[s].out[q] << endl;
297// cout << "Vout inout=" << Vout.data[s].in[q] << ", " << Vout.data[s].out[q] << endl;
298// print_size(Vin.data[s].block[q],"Vin.data[s].block[q]");
299// print_size(Vout.data[s].block[q],"Vout.data[s].block[q]");
300// cout << endl;
301// }
302
303 // project out unwanted states (e.g. to get lower spectrum)
304 for (size_t n=0; n<H.A0proj.size(); ++n)
305 {
306 Scalar overlap = 0;
307 for (size_t s=0; s<H.A0proj[n].size(); ++s)
308 {
309 // Note: Adjoint needed because we need <E0|Psi>, not <Psi|E0>
310 overlap += H.A0proj[n][s].adjoint().contract(Vin.data[s]).trace();
311 }
312// cout << "overlap=" << overlap << endl;
313
314 // explicit variant:
315 for (size_t s=0; s<H.A0proj[n].size(); ++s)
316 for (size_t q=0; q<H.A0proj[n][s].dim; ++q)
317 {
318 qarray2<Symmetry::Nq> cmp = {H.A0proj[n][s].in[q], H.A0proj[n][s].out[q]};
319 auto qA = Vout.data[s].dict.find(cmp);
320// assert(qA != Vout.data[s].dict.end() and "Error in HxV(PivotMatrix1,PivotVector): projected block not found!");
321 if (qA != Vout.data[s].dict.end() and H.A0proj[n][s].block[q].size() != 0)
322 {
323 Vout.data[s].block[qA->second] += H.Epenalty * overlap * H.A0proj[n][s].block[q];
324 }
325 }
326// // using Biped::addScale:
327// for (size_t s=0; s<H.A0proj[n].size(); ++s)
328// {
329// Vout.data[s].addScale(H.Epenalty*overlap, H.A0proj[n][s]);
330// }
331 }
332}
333
334//template<typename Symmetry, typename Scalar, typename MpoScalar>
335//void careful_HxV (const PivotMatrix1<Symmetry,Scalar,MpoScalar> &H, const PivotVector<Symmetry,Scalar> &Vin, PivotVector<Symmetry,Scalar> &Vout, std::array<qarray<Nq>,D> qloc)
336//{
337// Vout = Vin;
338// for (size_t s=0; s<D; ++s) {Vout.data[s].setZero();}
339//
340// // for (size_t s1=0; s1<D; ++s1)
370//
374// for (size_t s1=0; s1<D; ++s1)
375// for (size_t s2=0; s2<D; ++s2)
376// for (size_t qL=0; qL<H.L.dim; ++qL)
377// {
378// tuple<qarray3<Nq>,size_t,size_t,size_t> ix;
379// bool FOUND_MATCH = AWA(H.L.in(qL), H.L.out(qL), H.L.mid(qL), s1, s2, qloc, Vout.data, H.W, Vin.data, ix);
380//
381// if (FOUND_MATCH == true)
382// {
383// size_t q1 = get<1>(ix);
384// size_t qW = get<2>(ix);
385// size_t q2 = get<3>(ix);
386// auto qR = H.R.dict.find(get<0>(ix));
387//
388// if (qR != H.R.dict.end())
389// {
390// for (int k=0; k<H.W[s1][s2].block[qW].outerSize(); ++k)
391// for (SparseMatrixXd::InnerIterator iW(H.W[s1][s2].block[qW],k); iW; ++iW)
392// {
393// if (H.L.block[qL][iW.row()][0].rows() != 0 and
394// H.R.block[qR->second][iW.col()][0].rows() != 0)
395// {
396// if (Vout.data[s1].block[q1].rows() != H.L.block[qL][iW.row()][0].rows() or
397// Vout.dataA[s1].block[q1].cols() != H.R.block[qR->second][iW.col()][0].cols())
398// {
399// Vout.data[s1].block[q1] = iW.value() *
400// (H.L.block[qL][iW.row()][0] *
401// Vin.data[s2].block[q2] *
402// H.R.block[qR->second][iW.col()][0]);
403// }
404// else
405// {
406// Vout.data[s1].block[q1] += iW.value() *
407// (H.L.block[qL][iW.row()][0] *
408// Vin.data[s2].block[q2] *
409// H.R.block[qR->second][iW.col()][0]);
410// }
411// }
412// }
413// }
414// }
415// }
416//}
417
418template<typename Symmetry, typename Scalar, typename MpoScalar>
419void HxV (const PivotMatrix1<Symmetry,Scalar,MpoScalar> &H, PivotVector<Symmetry,Scalar> &Vinout)
420{
421 PivotVector<Symmetry,Scalar> Vtmp;
422 HxV(H,Vinout,Vtmp);
423 Vinout = Vtmp;
424}
425
426template<typename Symmetry, typename Scalar, typename MpoScalar>
427inline size_t dim (const PivotMatrix1<Symmetry,Scalar,MpoScalar> &H)
428{
429 return 0;
430}
431
432// How to calculate the Frobenius norm of this?
433template<typename Symmetry, typename Scalar, typename MpoScalar>
434inline double norm (const PivotMatrix1<Symmetry,Scalar,MpoScalar> &H)
435{
436 return H.dim;
437}
438
439#endif
dim
size_t dim(const PivotMatrix1< Symmetry, Scalar, MpoScalar > &H)
Definition: DmrgPivotMatrix1.h:427
norm
double norm(const PivotMatrix1< Symmetry, Scalar, MpoScalar > &H)
Definition: DmrgPivotMatrix1.h:434
HxV
void HxV(const PivotMatrix1< Symmetry, Scalar, MpoScalar > &H, const PivotVector< Symmetry, Scalar > &Vin, PivotVector< Symmetry, Scalar > &Vout)
Definition: DmrgPivotMatrix1.h:124
OxV
void OxV(const PivotMatrix1< Symmetry, Scalar, MpoScalar > &H, const PivotVector< Symmetry, Scalar > &Vin, PivotVector< Symmetry, Scalar > &Vout)
Definition: DmrgPivotMatrix1.h:132
qarray2
std::array< qarray< Nq >, 2 > qarray2
Definition: qarray.h:51
PivotMatrix1Terms::W
vector< vector< vector< Biped< Symmetry, Eigen::SparseMatrix< MpoScalar, Eigen::ColMajor, EIGEN_DEFAULT_SPARSE_INDEX_TYPE > > > > > W
Definition: DmrgPivotMatrix1.h:16
PivotMatrix1Terms::qloc
vector< qarray< Symmetry::Nq > > qloc
Definition: DmrgPivotMatrix1.h:18
PivotMatrix1Terms::L
Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > L
Definition: DmrgPivotMatrix1.h:14
PivotMatrix1Terms::qOp
vector< qarray< Symmetry::Nq > > qOp
Definition: DmrgPivotMatrix1.h:19
PivotMatrix1Terms::R
Tripod< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > R
Definition: DmrgPivotMatrix1.h:15
PivotMatrix1::qrhs
vector< vector< std::array< size_t, 6 > > > qrhs
Definition: DmrgPivotMatrix1.h:60
PivotMatrix1::qlhs
vector< std::array< size_t, 2 > > qlhs
Definition: DmrgPivotMatrix1.h:59
PivotMatrix1::PL
vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > PL
Definition: DmrgPivotMatrix1.h:65
PivotMatrix1::memory
double memory(MEMUNIT memunit) const
Definition: DmrgPivotMatrix1.h:105
PivotMatrix1::PR
vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > PR
Definition: DmrgPivotMatrix1.h:66
PivotMatrix1::A0proj
vector< vector< Biped< Symmetry, Matrix< Scalar, Dynamic, Dynamic > > > > A0proj
Definition: DmrgPivotMatrix1.h:67
PivotMatrix1::cast
PivotMatrix1< Symmetry, OtherScalar, OtherScalar > cast() const
Definition: DmrgPivotMatrix1.h:72
PivotMatrix1::Terms
vector< PivotMatrix1Terms< Symmetry, Scalar, MpoScalar > > Terms
Definition: DmrgPivotMatrix1.h:56
PivotMatrix1::factor_cgcs
vector< vector< Scalar > > factor_cgcs
Definition: DmrgPivotMatrix1.h:61
PivotVector::size
size_t size() const
Definition: DmrgPivotVector.h:74
PivotVector::data
vector< Biped< Symmetry, Matrix< Scalar_, Dynamic, Dynamic > > > data
Definition: DmrgPivotVector.h:116