NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_utilities_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_UTILITIES_IMPL_H
2 #define CNTR_UTILITIES_IMPL_H
3 
5 //#include "cntr_exception.hpp"
6 #include "cntr_elements.hpp"
12 #include "cntr_function_decl.hpp"
14 
15 namespace cntr {
16 
18 // K-th order polynomila extrapolate of realtime functions (G^ret, G^vt,
19 // G^les)
20 // to time t=(n+1), using information at times t=(n-j) [j=0...k]
21 template <typename T, class GG, int SIZE1>
22 void extrapolate_timestep_dispatch(int n, GG &G,
24  typedef std::complex<T> cplx;
25  T *p1;
26  cplx z1, *gtemp, *gtemp1;
27  int m, l, j, ntau, nt, n1, k, sg, size1 = G.size1();
28  ntau = G.ntau();
29  nt = G.nt();
30  n1 = n + 1;
31  k = I.k();
32  sg = G.element_size();
33  gtemp = new cplx[sg];
34  gtemp1 = new cplx[sg];
35  if (n1 > nt || n < k) {
36  std::cerr << " k= " << k << " n= " << n << " nt= " << nt << std::endl;
37  std::cerr << "extrapolate_tstep: n out of range" << std::endl;
38  abort();
39  }
40  if (k == 0) {
41  for (m = 0; m <= ntau; m++)
42  element_set<T, SIZE1>(size1, G.tvptr(n1, m), G.tvptr(n, m));
43  for (j = 0; j <= n; j++)
44  element_set<T, SIZE1>(size1, G.retptr(n1, n1 - j),
45  G.retptr(n, n - j));
46  element_set<T, SIZE1>(size1, G.retptr(n1, 0), G.retptr(n, 0));
47  for (j = 1; j <= n1; j++)
48  element_set<T, SIZE1>(size1, G.lesptr(j, n1), G.lesptr(j - 1, n));
49  element_set<T, SIZE1>(size1, G.lesptr(0, n1), G.lesptr(0, n));
50  } else {
51  p1 = new T[k + 1];
52  for (m = 0; m <= k; m++) {
53  p1[m] = 0.0;
54  for (l = 0; l <= k; l++)
55  p1[m] += (1 - 2 * (l % 2)) * I.poly_interpolation(l, m);
56  }
57  // vt-component
58  for (m = 0; m <= ntau; m++) {
59  element_set_zero<T, SIZE1>(size1, gtemp);
60  for (l = 0; l <= k; l++)
61  element_incr<T, SIZE1>(size1, gtemp, p1[l],
62  G.tvptr(n - l, m));
63  element_set<T, SIZE1>(size1, G.tvptr(n1, m), gtemp);
64  }
65  // ret-component
66  for (j = 0; j <= n - k; j++) {
67  element_set_zero<T, SIZE1>(size1, gtemp);
68  for (l = 0; l <= k; l++)
69  element_incr<T, SIZE1>(size1, gtemp, p1[l],
70  G.retptr(n - l, n - l - j));
71  element_set<T, SIZE1>(size1, G.retptr(n1, n1 - j), gtemp);
72  }
73  for (j = 0; j <= k; j++) {
74  element_set_zero<T, SIZE1>(size1, gtemp);
75  for (l = 0; l <= k; l++) {
76  if (n - l <= j) {
77  element_set<T, SIZE1>(size1, gtemp1, G.retptr(j, n - l));
78  element_conj<T, SIZE1>(size1, gtemp1);
79  element_smul<T, SIZE1>(size1, gtemp1, -1);
80  } else {
81  element_set<T, SIZE1>(size1, gtemp1, G.retptr(n - l, j));
82  }
83  element_incr<T, SIZE1>(size1, gtemp, p1[l], gtemp1);
84  }
85  element_set<T, SIZE1>(size1, G.retptr(n1, j), gtemp);
86  }
87  // les-component
88  for (j = 0; j <= k; j++) {
89  element_set_zero<T, SIZE1>(size1, gtemp);
90  for (l = 0; l <= k; l++) {
91  if (n - l < j) {
92  element_set<T, SIZE1>(size1, gtemp1, G.lesptr(n - l, j));
93  element_conj<T, SIZE1>(size1, gtemp1);
94  element_smul<T, SIZE1>(size1, gtemp1, -1);
95  } else {
96  element_set<T, SIZE1>(size1, gtemp1, G.lesptr(j, n - l));
97  }
98  element_incr<T, SIZE1>(size1, gtemp, p1[l], gtemp1);
99  }
100  element_set<T, SIZE1>(size1, G.lesptr(j, n1), gtemp);
101  }
102  for (j = k + 1; j <= n1; j++) {
103  element_set_zero<T, SIZE1>(size1, gtemp);
104  for (l = 0; l <= k; l++)
105  element_incr<T, SIZE1>(size1, gtemp, p1[l],
106  G.lesptr(j - l - 1, n - l));
107  element_set<T, SIZE1>(size1, G.lesptr(j, n1), gtemp);
108  }
109  delete[] p1;
110  }
111  delete[] gtemp;
112  delete[] gtemp1;
113 }
115 
134 template <typename T>
135 void extrapolate_timestep(int n, herm_matrix<T> &G,
137  if (G.size1() == 1)
138  extrapolate_timestep_dispatch<T, herm_matrix<T>, 1>(n, G, I);
139  else
140  extrapolate_timestep_dispatch<T, herm_matrix<T>, LARGESIZE>(n, G, I);
141 }
142 
143 
163 template <typename T>
164 void extrapolate_timestep(int n, herm_matrix<T> &G, int SolveOrder) {
165  if (G.size1() == 1)
166  extrapolate_timestep_dispatch<T, herm_matrix<T>, 1>(n, G, integration::I<T>(SolveOrder));
167  else
168  extrapolate_timestep_dispatch<T, herm_matrix<T>, LARGESIZE>(n, G, integration::I<T>(SolveOrder));
169 }
170 
173 template <typename T>
174 void extrapolate_timestep(int n, herm_pseudo<T> &G,
176  if (G.size1() == 1)
177  extrapolate_timestep_dispatch<T, herm_pseudo<T>, 1>(n, G, I);
178  else
179  extrapolate_timestep_dispatch<T, herm_pseudo<T>, LARGESIZE>(n, G, I);
180 }
181 
182 
184 // k-th order polynomial extrapolate of contour function
185 // to time t=(n+1), using information at times t=(n-j) [j=0...k]
186 template <typename T, int SIZE1>
187 void extrapolate_timestep_dispatch(int n, cntr::function<T> &f,
189 
190  typedef std::complex<T> cplx;
191  int size1 = f.size1_;
192  int size2 = f.size2_;
193  int nt = f.nt_;
194  int kt = I.k();
195  assert(n+1<=nt);
196  assert(n>=kt);
197  cdmatrix value(size1,size2),output(size1,size2);
198  value.setZero();
199  output.setZero();
200  if (kt == 0) {
201  f.get_value(n,value);
202  output=value;
203  f.set_value(n+1,output);
204  } else {
205  std::vector<double> weight(kt+1);
206  // Set up weights for extrapolation from the interpolated ones
207  for(int m=0; m<=kt; m++){
208  weight[m]=0.0;
209  for(int l=0; l <=kt;l++){
210  weight[m] += (1 - 2 * (l % 2))*I.poly_interpolation(l, m);
211  }
212  }
213  for(int l=0;l<= kt;l++){
214  f.get_value(n-l,value);
215  // std::cout << "extra2d " << value << " " << weight[l] << std::endl;
216  output+=weight[l]*value;
217  // std::cout << "extra2e " << output << std::endl;
218  }
219  // std::cout << "extra2f " << output << std::endl;
220  f.set_value(n+1,output);
221  // std::cout << "extra2g " << output << std::endl;
222  }
223 }
225 
242 template <typename T>
243 void extrapolate_timestep(int n, function<T> &f,
245  if (f.size1() == 1)
246  extrapolate_timestep_dispatch<T, 1>(n, f, I);
247  else
248  extrapolate_timestep_dispatch<T, LARGESIZE>(n, f, I);
249 }
250 
251 
269 template <typename T>
270 void extrapolate_timestep(int n, function<T> &f, int ExtrapolationOrder) {
271  if (f.size1() == 1)
272  extrapolate_timestep_dispatch<T, 1>(n, f, integration::I<T>(ExtrapolationOrder));
273  else
274  extrapolate_timestep_dispatch<T, LARGESIZE>(n, f, integration::I<T>(ExtrapolationOrder));
275 }
276 
277 
279 // k-th order polynomial interpolation of contour function
280 // to arbitrary point on time interval [tstp-kt,tstp],
281 // using information at times t=(n-j) [j=0...k]
282 // TODO: test this more throughly
283 // TODO: can we do something smart for edge points ?
284 template <typename T, int SIZE1>
285 cdmatrix interpolation_dispatch(int tstp, double tinter, cntr::function<T> &f,
287  typedef std::complex<T> cplx;
288  int kt=I.k();
289  int size=f.size1();
290  cdmatrix output(size,size),tmp(size,size);
291  output.setZero();
292  double t0=float(tstp-kt);
293  double t1;
294 
295  // The k-th order approximation polynomial through function
296  // values {(x_l,f(x_l)),l=0...k} is given by
297  // P(x) = sum_{p,l=0}^{k} x^p*I.interpolation(p,l)*f_l
298 
299  for(int l=0;l<=kt;l++){
300  t1=1.0;
301  double weight = I.poly_interpolation(0,l);
302  f.get_value(tstp-kt+l,tmp);
303  for(int p=1;p<=kt;p++){
304  t1*=tinter-t0;
305  weight +=t1 * I.poly_interpolation(p,l);
306  }
307 
308  output+=tmp*weight;
309 
310  }
311  return output;
312 }
314 
336 template <typename T>
337 cdmatrix interpolation(int tstp,double tinter,function<T> &f,
339 
340  int kt=I.k();
341  assert(tstp>=kt);
342  assert(tinter<=tstp);
343  assert(tinter>=(tstp-kt));
344 
345  if (f.size1() == 1)
346  return interpolation_dispatch<T, 1>(tstp,tinter, f, I);
347  else
348  return interpolation_dispatch<T, LARGESIZE>(tstp,tinter, f, I);
349 }
350 
373 template <typename T>
374 cdmatrix interpolation(int tstp,double tinter,function<T> &f,
375  int InterpolationOrder) {
376 
377  assert(tstp>=InterpolationOrder);
378  assert(tinter<=tstp);
379  assert(tinter>=(tstp-InterpolationOrder));
380 
381  if (f.size1() == 1)
382  return interpolation_dispatch<T, 1>(tstp,tinter, f, integration::I<T>(InterpolationOrder));
383  else
384  return interpolation_dispatch<T, LARGESIZE>(tstp,tinter, f, integration::I<T>(InterpolationOrder));
385 }
386 
399 template <typename T>
401  int m, ntau = G.ntau(), size1 = G.size1();
402  for (m = 0; m <= ntau; m++) {
403  element_set<T, LARGESIZE>(size1, G.tvptr(0, m), G.matptr(ntau - m));
404  element_smul<T, LARGESIZE>(size1, G.tvptr(0, m),
405  std::complex<T>(0, G.sig()));
406  }
407  element_set<T, LARGESIZE>(size1, G.lesptr(0, 0), G.tvptr(0, 0));
408  element_set<T, LARGESIZE>(size1, G.retptr(0, 0), G.matptr(0));
409  element_smul<T, LARGESIZE>(size1, G.retptr(0, 0), std::complex<T>(0, 1));
410  element_incr<T, LARGESIZE>(size1, G.retptr(0, 0), -1, G.lesptr(0, 0));
411 }
412 
413 
414 // Fix first kt points from mat
432 template <typename T>
434  int m,n,ntau=G.ntau(),size1=G.size1();
435  for(n=0;n<=kt;n++){
436  for(m=0;m<=ntau;m++){
437  element_set<T,LARGESIZE>(size1,G.tvptr(n,m),G.matptr(ntau-m));
438  element_smul<T,LARGESIZE>(size1,G.tvptr(n,m),std::complex<T>(0,G.sig()));
439  }
440  element_set<T,LARGESIZE>(size1,G.lesptr(n,n),G.tvptr(0,0));
441  element_set<T,LARGESIZE>(size1,G.retptr(n,n),G.matptr(0));
442  element_smul<T,LARGESIZE>(size1,G.retptr(n,n),std::complex<T>(0,1));
443  element_incr<T,LARGESIZE>(size1,G.retptr(n,n),-1,G.lesptr(0,0));
444  }
445 }
446 
448 template <typename T>
449 void set_t0_from_mat(herm_pseudo<T> &G) {
450  int m, ntau = G.ntau(), size1 = G.size1();
451  for (m = 0; m <= ntau; m++) {
452  element_set<T, LARGESIZE>(size1, G.tvptr(0, m), G.matptr(ntau - m));
453  element_smul<T, LARGESIZE>(size1, G.tvptr(0, m),
454  std::complex<T>(0, G.sig()));
455  }
456  element_set<T, LARGESIZE>(size1, G.lesptr(0, 0), G.tvptr(0, 0));
457  element_set<T, LARGESIZE>(size1, G.retptr(0, 0), G.matptr(0));
458  element_smul<T, LARGESIZE>(size1, G.retptr(0, 0), std::complex<T>(0, 1));
459  // element_incr<T,LARGESIZE>(size1,G.retptr(0,0),-1,G.lesptr(0,0));
460 }
461 
489 template <typename T>
491  integration::Integrator<T> &I, T beta, T h){
492  int size1=G.size1();
493  std::complex<T> trGxSGM;
494  std::complex<T> *GxSGM = new std::complex<T>[size1*size1];
495 
496  convolution_density_matrix<T,herm_matrix<T>>(tstp, GxSGM, Sigma, G, I, beta, h);
497  trGxSGM = (0.0,0.0);
498  for(int i=0; i< size1; i++){
499  trGxSGM += GxSGM[i*size1 + i];
500  }
501  delete GxSGM;
502  return 0.5*trGxSGM.real();
503 }
504 
532 template <typename T>
534  T beta, T h, int SolveOrder){
535  int size1=G.size1();
536  std::complex<T> trGxSGM;
537  std::complex<T> *GxSGM = new std::complex<T>[size1*size1];
538 
539  convolution_density_matrix<T,herm_matrix<T>>(tstp, GxSGM, Sigma, G, integration::I<T>(SolveOrder), beta, h);
540  trGxSGM = (0.0,0.0);
541  for(int i=0; i< size1; i++){
542  trGxSGM += GxSGM[i*size1 + i];
543  }
544  return 0.5*trGxSGM.real();
545 
546  delete GxSGM;
547 }
548 
549 
551 template <typename T, class GG, int SIZE1>
552 T distance_norm2_ret_dispatch(int tstp, GG &g1, GG &g2) {
553  int size1 = g1.size1(), i;
554  T err = 0.0;
555  std::complex<T> *temp = new std::complex<T>[size1 * size1];
556  assert(tstp>=-1 && g1.nt()>=tstp && g2.nt()>=tstp && g1.size1() == g2.size1());
557 
558  if (tstp == -1)
559  return 0.0;
560  for (i = 0; i <= tstp; i++) {
561  element_set<T, SIZE1>(size1, temp, g1.retptr(tstp, i));
562  element_incr<T, SIZE1>(size1, temp, -1.0, g2.retptr(tstp, i));
563  err += element_norm2<T, SIZE1>(size1, temp);
564  }
565  delete[] temp;
566  return err;
567 }
569 template <typename T, class GG, int SIZE1>
570 T distance_norm2_tv_dispatch(int tstp, GG &g1, GG &g2) {
571  int size1 = g1.size1(), ntau = g1.ntau(), i;
572  T err = 0.0;
573  std::complex<T> *temp = new std::complex<T>[size1 * size1];
574  assert(tstp>=-1 && g1.nt()>=tstp && g2.nt() >= tstp);
575  assert(g1.ntau()==g2.ntau() && g1.size1()==g2.size1());
576  if (tstp == -1)
577  return 0.0;
578  for (i = 0; i <= ntau; i++) {
579  element_set<T, SIZE1>(size1, temp, g1.tvptr(tstp, i));
580  element_incr<T, SIZE1>(size1, temp, -1.0, g2.tvptr(tstp, i));
581  err += element_norm2<T, SIZE1>(size1, temp);
582  }
583  delete[] temp;
584  return err;
585 }
587 template <typename T, class GG, int SIZE1>
588 T distance_norm2_les_dispatch(int tstp, GG &g1, GG &g2) {
589  int size1 = g1.size1(), i;
590  T err = 0.0;
591  std::complex<T> *temp = new std::complex<T>[size1 * size1];
592 
593  assert(tstp>=-1 && g1.nt() >= tstp && g2.nt() >= tstp && g1.size1() == g2.size1());
594  if (tstp == -1)
595  return 0.0;
596  for (i = 0; i <= tstp; i++) {
597  element_set<T, SIZE1>(size1, temp, g1.lesptr(i, tstp));
598  element_incr<T, SIZE1>(size1, temp, -1.0, g2.lesptr(i, tstp));
599  err += element_norm2<T, SIZE1>(size1, temp);
600  }
601  delete[] temp;
602  return err;
603 }
604 
606 template <typename T, class GG, int SIZE1>
607 T distance_norm2_dispatch(int tstp, GG &g1, GG &g2) {
608  int size1 = g1.size1(), ntau = g1.ntau(), i;
609  T err = 0.0;
610  std::complex<T> *temp = new std::complex<T>[size1 * size1];
611 
612  assert(tstp >= -1 && g1.nt() >= tstp && g2.nt() >= tstp);
613  assert(g1.ntau() == g2.ntau());
614  assert(g1.size1() == g2.size1());
615 
616  if (tstp == -1) {
617  for (i = 0; i <= ntau; i++) {
618  element_set<T, SIZE1>(size1, temp, g1.matptr(i));
619  element_incr<T, SIZE1>(size1, temp, -1.0, g2.matptr(i));
620  err += element_norm2<T, SIZE1>(size1, temp);
621  }
622  } else {
623  for (i = 0; i <= tstp; i++) {
624  element_set<T, SIZE1>(size1, temp, g1.retptr(tstp, i));
625  element_incr<T, SIZE1>(size1, temp, -1.0, g2.retptr(tstp, i));
626  err += element_norm2<T, SIZE1>(size1, temp);
627  }
628  for (i = 0; i <= ntau; i++) {
629  element_set<T, SIZE1>(size1, temp, g1.tvptr(tstp, i));
630  element_incr<T, SIZE1>(size1, temp, -1.0, g2.tvptr(tstp, i));
631  err += element_norm2<T, SIZE1>(size1, temp);
632  }
633  for (i = 0; i <= tstp; i++) {
634  element_set<T, SIZE1>(size1, temp, g1.lesptr(i, tstp));
635  element_incr<T, SIZE1>(size1, temp, -1.0, g2.lesptr(i, tstp));
636  err += element_norm2<T, SIZE1>(size1, temp);
637  }
638  }
639  delete[] temp;
640  return err;
641 }
642 
643 
644 
645 
646 
647 
648 
649 
651 template <typename T, int SIZE1>
652 T distance_norm2_ret_dispatch(int tstp, herm_matrix_timestep_view<T> &g1, herm_matrix_timestep_view<T> &g2) {
653  int size1 = g1.size1();
654  int s1 = size1*size1;
655  T err = 0.0;
656  std::complex<T> *temp = new std::complex<T>[size1 * size1];
657 
658  if (tstp == -1)
659  return 0.0;
660  for (int i = 0; i <= tstp; i++) {
661  element_set<T, SIZE1>(size1, temp, g1.ret_ + i * s1);
662  element_incr<T, SIZE1>(size1, temp, -1.0, g2.ret_ + i * s1);
663  err += element_norm2<T, SIZE1>(size1, temp);
664  }
665  delete[] temp;
666  return err;
667 }
669 template <typename T, int SIZE1>
670 T distance_norm2_tv_dispatch(int tstp, herm_matrix_timestep_view<T> &g1, herm_matrix_timestep_view<T> &g2) {
671  int size1 = g1.size1(), ntau = g1.ntau();
672  int s1 = size1*size1;
673  T err = 0.0;
674  std::complex<T> *temp = new std::complex<T>[size1 * size1];
675 
676  if (tstp == -1)
677  return 0.0;
678  for (int i = 0; i <= ntau; i++) {
679  element_set<T, SIZE1>(size1, temp, g1.tv_ + i * s1);
680  element_incr<T, SIZE1>(size1, temp, -1.0, g2.tv_ + i * s1);
681  err += element_norm2<T, SIZE1>(size1, temp);
682  }
683  delete[] temp;
684  return err;
685 }
687 template <typename T, int SIZE1>
688 T distance_norm2_les_dispatch(int tstp, herm_matrix_timestep_view<T> &g1, herm_matrix_timestep_view<T> &g2) {
689  int size1 = g1.size1();
690  int s1 = size1*size1;
691  T err = 0.0;
692  std::complex<T> *temp = new std::complex<T>[size1 * size1];
693 
694  if (tstp == -1)
695  return 0.0;
696  for (int i = 0; i <= tstp; i++) {
697  element_set<T, SIZE1>(size1, temp, g1.les_ + i * s1);
698  element_incr<T, SIZE1>(size1, temp, -1.0, g2.les_ + i * s1);
699  err += element_norm2<T, SIZE1>(size1, temp);
700  }
701  delete[] temp;
702  return err;
703 }
704 
706 template <typename T, int SIZE1>
707 T distance_norm2_dispatch(int tstp, herm_matrix_timestep_view<T> &g1, herm_matrix_timestep_view<T> &g2) {
708  int size1 = g1.size1(), ntau = g1.ntau();
709  int s1 = size1*size1;
710  T err = 0.0;
711  std::complex<T> *temp = new std::complex<T>[size1 * size1];
712 
713  if (tstp == -1) {
714  for (int i = 0; i <= ntau; i++) {
715  element_set<T, SIZE1>(size1, temp, g1.mat_ + i * s1);
716  element_incr<T, SIZE1>(size1, temp, -1.0, g2.mat_ + i * s1);
717  err += element_norm2<T, SIZE1>(size1, temp);
718  }
719  } else {
720  for (int i = 0; i <= tstp; i++) {
721  element_set<T, SIZE1>(size1, temp, g1.ret_ + i * s1);
722  element_incr<T, SIZE1>(size1, temp, -1.0, g2.ret_ + i * s1);
723  err += element_norm2<T, SIZE1>(size1, temp);
724  }
725  for (int i = 0; i <= ntau; i++) {
726  element_set<T, SIZE1>(size1, temp, g1.tv_ + i * s1);
727  element_incr<T, SIZE1>(size1, temp, -1.0, g2.tv_ + i * s1);
728  err += element_norm2<T, SIZE1>(size1, temp);
729  }
730  for (int i = 0; i <= tstp; i++) {
731  element_set<T, SIZE1>(size1, temp, g1.les_ + i * s1);
732  element_incr<T, SIZE1>(size1, temp, -1.0, g2.les_ + i * s1);
733  err += element_norm2<T, SIZE1>(size1, temp);
734  }
735  }
736  delete[] temp;
737  return err;
738 }
739 
740 
759 template <typename T>
761  assert(g1.size1() == g2.size1());
762  assert(g1.ntau() == g2.ntau());
763  assert(g1.nt() >= tstp);
764  assert(g2.nt() >= tstp);
765  if (g1.size1() == 1)
766  return distance_norm2_ret_dispatch<T, herm_matrix<T>, 1>(tstp, g1,
767  g2);
768  else
769  return distance_norm2_ret_dispatch<T, herm_matrix<T>, LARGESIZE>(
770  tstp, g1, g2);
771 }
772 
791 template <typename T>
793  assert(g1.size1() == g2.size1());
794  assert(g1.ntau() == g2.ntau());
795  assert(g1.tstp() == tstp);
796  assert(g2.nt() >= tstp);
797  herm_matrix_timestep_view<T> g1_view(tstp, g1);
798  herm_matrix_timestep_view<T> g2_view(tstp, g2);
799  if (g1.size1() == 1)
800  return distance_norm2_ret_dispatch<T, 1>(tstp, g1_view, g2_view);
801  else
802  return distance_norm2_ret_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
803 }
804 
823 template <typename T>
825  assert(g1.size1() == g2.size1());
826  assert(g1.ntau() == g2.ntau());
827  assert(g1.nt() >= tstp);
828  assert(g2.tstp() == tstp);
829  herm_matrix_timestep_view<T> g1_view(tstp, g1);
830  herm_matrix_timestep_view<T> g2_view(tstp, g2);
831  if (g1.size1() == 1)
832  return distance_norm2_ret_dispatch<T, 1>(tstp, g1_view, g2_view);
833  else
834  return distance_norm2_ret_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
835 }
836 
837 
856 template <typename T>
858  assert(g1.size1() == g2.size1());
859  assert(g1.ntau() == g2.ntau());
860  assert(g1.tstp() == tstp);
861  assert(g2.tstp() == tstp);
862  herm_matrix_timestep_view<T> g1_view(tstp, g1);
863  herm_matrix_timestep_view<T> g2_view(tstp, g2);
864  if (g1.size1() == 1)
865  return distance_norm2_ret_dispatch<T, 1>(tstp, g1_view, g2_view);
866  else
867  return distance_norm2_ret_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
868 }
869 
888 template <typename T>
890  assert(g1.size1() == g2.size1());
891  assert(g1.ntau() == g2.ntau());
892  assert(g1.tstp() == tstp);
893  assert(g2.nt() >= tstp);
894  herm_matrix_timestep_view<T> g2_view(tstp, g2);
895  if (g1.size1() == 1)
896  return distance_norm2_ret_dispatch<T, 1>(tstp, g1, g2_view);
897  else
898  return distance_norm2_ret_dispatch<T, LARGESIZE>(tstp, g1, g2_view);
899 }
900 
919 template <typename T>
921  assert(g1.size1() == g2.size1());
922  assert(g1.ntau() == g2.ntau());
923  assert(g1.nt() >= tstp);
924  assert(g2.tstp() == tstp);
925  herm_matrix_timestep_view<T> g1_view(tstp, g1);
926  if (g1.size1() == 1)
927  return distance_norm2_ret_dispatch<T, 1>(tstp, g1_view, g2);
928  else
929  return distance_norm2_ret_dispatch<T, LARGESIZE>(tstp, g1_view, g2);
930 }
931 
932 
951 template <typename T>
953  assert(g1.size1() == g2.size1());
954  assert(g1.ntau() == g2.ntau());
955  assert(g1.tstp() == tstp);
956  assert(g2.tstp() == tstp);
957  herm_matrix_timestep_view<T> g1_view(tstp, g1);
958  if (g1.size1() == 1)
959  return distance_norm2_ret_dispatch<T, 1>(tstp, g1_view, g2);
960  else
961  return distance_norm2_ret_dispatch<T, LARGESIZE>(tstp, g1_view, g2);
962 }
963 
982 template <typename T>
984  assert(g1.size1() == g2.size1());
985  assert(g1.ntau() == g2.ntau());
986  assert(g1.tstp() == tstp);
987  assert(g2.tstp() == tstp);
988  herm_matrix_timestep_view<T> g2_view(tstp, g2);
989  if (g1.size1() == 1)
990  return distance_norm2_ret_dispatch<T, 1>(tstp, g1, g2_view);
991  else
992  return distance_norm2_ret_dispatch<T, LARGESIZE>(tstp, g1, g2_view);
993 }
994 
995 
1014 template <typename T>
1016  assert(g1.size1() == g2.size1());
1017  assert(g1.ntau() == g2.ntau());
1018  assert(g1.tstp() == tstp);
1019  assert(g2.tstp() == tstp);
1020  if (g1.size1() == 1)
1021  return distance_norm2_ret_dispatch<T, 1>(tstp, g1, g2);
1022  else
1023  return distance_norm2_ret_dispatch<T, LARGESIZE>(tstp, g1, g2);
1024 }
1025 
1026 
1027 
1046 template <typename T>
1048  assert(g1.size1() == g2.size1());
1049  assert(g1.ntau() == g2.ntau());
1050  assert(g1.nt() >= tstp);
1051  assert(g2.nt() >= tstp);
1052  if (g1.size1() == 1)
1053  return distance_norm2_tv_dispatch<T, herm_matrix<T>, 1>(tstp, g1, g2);
1054  else
1055  return distance_norm2_tv_dispatch<T, herm_matrix<T>, LARGESIZE>(
1056  tstp, g1, g2);
1057 }
1058 
1077 template <typename T>
1079  assert(g1.size1() == g2.size1());
1080  assert(g1.ntau() == g2.ntau());
1081  assert(g1.tstp() == tstp);
1082  assert(g2.nt() >= tstp);
1083  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1084  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1085  if (g1.size1() == 1)
1086  return distance_norm2_tv_dispatch<T, 1>(tstp, g1_view, g2_view);
1087  else
1088  return distance_norm2_tv_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
1089 }
1090 
1109 template <typename T>
1111  assert(g1.size1() == g2.size1());
1112  assert(g1.ntau() == g2.ntau());
1113  assert(g1.nt() >= tstp);
1114  assert(g2.tstp() == tstp);
1115  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1116  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1117  if (g1.size1() == 1)
1118  return distance_norm2_tv_dispatch<T, 1>(tstp, g1_view, g2_view);
1119  else
1120  return distance_norm2_tv_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
1121 }
1122 
1123 
1142 template <typename T>
1144  assert(g1.size1() == g2.size1());
1145  assert(g1.ntau() == g2.ntau());
1146  assert(g1.tstp() == tstp);
1147  assert(g2.tstp() == tstp);
1148  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1149  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1150  if (g1.size1() == 1)
1151  return distance_norm2_tv_dispatch<T, 1>(tstp, g1_view, g2_view);
1152  else
1153  return distance_norm2_tv_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
1154 }
1155 
1174 template <typename T>
1176  assert(g1.size1() == g2.size1());
1177  assert(g1.ntau() == g2.ntau());
1178  assert(g1.tstp() == tstp);
1179  assert(g2.nt() >= tstp);
1180  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1181  if (g1.size1() == 1)
1182  return distance_norm2_tv_dispatch<T, 1>(tstp, g1, g2_view);
1183  else
1184  return distance_norm2_tv_dispatch<T, LARGESIZE>(tstp, g1, g2_view);
1185 }
1186 
1205 template <typename T>
1207  assert(g1.size1() == g2.size1());
1208  assert(g1.ntau() == g2.ntau());
1209  assert(g1.nt() >= tstp);
1210  assert(g2.tstp() == tstp);
1211  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1212  if (g1.size1() == 1)
1213  return distance_norm2_tv_dispatch<T, 1>(tstp, g1_view, g2);
1214  else
1215  return distance_norm2_tv_dispatch<T, LARGESIZE>(tstp, g1_view, g2);
1216 }
1217 
1218 
1237 template <typename T>
1239  assert(g1.size1() == g2.size1());
1240  assert(g1.ntau() == g2.ntau());
1241  assert(g1.tstp() == tstp);
1242  assert(g2.tstp() == tstp);
1243  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1244  if (g1.size1() == 1)
1245  return distance_norm2_tv_dispatch<T, 1>(tstp, g1_view, g2);
1246  else
1247  return distance_norm2_tv_dispatch<T, LARGESIZE>(tstp, g1_view, g2);
1248 }
1249 
1268 template <typename T>
1270  assert(g1.size1() == g2.size1());
1271  assert(g1.ntau() == g2.ntau());
1272  assert(g1.tstp() == tstp);
1273  assert(g2.tstp() == tstp);
1274  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1275  if (g1.size1() == 1)
1276  return distance_norm2_tv_dispatch<T, 1>(tstp, g1, g2_view);
1277  else
1278  return distance_norm2_tv_dispatch<T, LARGESIZE>(tstp, g1, g2_view);
1279 }
1280 
1281 
1300 template <typename T>
1302  assert(g1.size1() == g2.size1());
1303  assert(g1.ntau() == g2.ntau());
1304  assert(g1.tstp() == tstp);
1305  assert(g2.tstp() == tstp);
1306  if (g1.size1() == 1)
1307  return distance_norm2_tv_dispatch<T, 1>(tstp, g1, g2);
1308  else
1309  return distance_norm2_tv_dispatch<T, LARGESIZE>(tstp, g1, g2);
1310 }
1311 
1312 
1331 template <typename T>
1333  assert(g1.size1() == g2.size1());
1334  assert(g1.ntau() == g2.ntau());
1335  assert(g1.nt() >= tstp);
1336  assert(g2.nt() >= tstp);
1337  if (g1.size1() == 1)
1338  return distance_norm2_les_dispatch<T, herm_matrix<T>, 1>(tstp, g1,
1339  g2);
1340  else
1341  return distance_norm2_les_dispatch<T, herm_matrix<T>, LARGESIZE>(
1342  tstp, g1, g2);
1343 }
1344 
1345 
1346 
1365 template <typename T>
1367  assert(g1.size1() == g2.size1());
1368  assert(g1.ntau() == g2.ntau());
1369  assert(g1.tstp() == tstp);
1370  assert(g2.nt() >= tstp);
1371  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1372  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1373  if (g1.size1() == 1)
1374  return distance_norm2_les_dispatch<T, 1>(tstp, g1_view, g2_view);
1375  else
1376  return distance_norm2_les_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
1377 }
1378 
1397 template <typename T>
1399  assert(g1.size1() == g2.size1());
1400  assert(g1.ntau() == g2.ntau());
1401  assert(g1.nt() >= tstp);
1402  assert(g2.tstp() == tstp);
1403  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1404  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1405  if (g1.size1() == 1)
1406  return distance_norm2_les_dispatch<T, 1>(tstp, g1_view, g2_view);
1407  else
1408  return distance_norm2_les_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
1409 }
1410 
1411 
1430 template <typename T>
1432  assert(g1.size1() == g2.size1());
1433  assert(g1.ntau() == g2.ntau());
1434  assert(g1.tstp() == tstp);
1435  assert(g2.tstp() == tstp);
1436  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1437  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1438  if (g1.size1() == 1)
1439  return distance_norm2_les_dispatch<T, 1>(tstp, g1_view, g2_view);
1440  else
1441  return distance_norm2_les_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
1442 }
1443 
1462 template <typename T>
1464  assert(g1.size1() == g2.size1());
1465  assert(g1.ntau() == g2.ntau());
1466  assert(g1.tstp() == tstp);
1467  assert(g2.nt() >= tstp);
1468  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1469  if (g1.size1() == 1)
1470  return distance_norm2_les_dispatch<T, 1>(tstp, g1, g2_view);
1471  else
1472  return distance_norm2_les_dispatch<T, LARGESIZE>(tstp, g1, g2_view);
1473 }
1474 
1493 template <typename T>
1495  assert(g1.size1() == g2.size1());
1496  assert(g1.ntau() == g2.ntau());
1497  assert(g1.nt() >= tstp);
1498  assert(g2.tstp() == tstp);
1499  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1500  if (g1.size1() == 1)
1501  return distance_norm2_les_dispatch<T, 1>(tstp, g1_view, g2);
1502  else
1503  return distance_norm2_les_dispatch<T, LARGESIZE>(tstp, g1_view, g2);
1504 }
1505 
1506 
1525 template <typename T>
1527  assert(g1.size1() == g2.size1());
1528  assert(g1.ntau() == g2.ntau());
1529  assert(g1.tstp() == tstp);
1530  assert(g2.tstp() == tstp);
1531  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1532  if (g1.size1() == 1)
1533  return distance_norm2_les_dispatch<T, 1>(tstp, g1_view, g2);
1534  else
1535  return distance_norm2_les_dispatch<T, LARGESIZE>(tstp, g1_view, g2);
1536 }
1537 
1556 template <typename T>
1558  assert(g1.size1() == g2.size1());
1559  assert(g1.ntau() == g2.ntau());
1560  assert(g1.tstp() == tstp);
1561  assert(g2.tstp() == tstp);
1562  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1563  if (g1.size1() == 1)
1564  return distance_norm2_les_dispatch<T, 1>(tstp, g1, g2_view);
1565  else
1566  return distance_norm2_les_dispatch<T, LARGESIZE>(tstp, g1, g2_view);
1567 }
1568 
1587 template <typename T>
1589  assert(g1.size1() == g2.size1());
1590  assert(g1.ntau() == g2.ntau());
1591  assert(g1.tstp() == tstp);
1592  assert(g2.tstp() == tstp);
1593  if (g1.size1() == 1)
1594  return distance_norm2_les_dispatch<T, 1>(tstp, g1, g2);
1595  else
1596  return distance_norm2_les_dispatch<T, LARGESIZE>(tstp, g1, g2);
1597 }
1598 
1599 
1600 
1620 template <typename T>
1622  assert(g1.size1() == g2.size1());
1623  assert(g1.ntau() == g2.ntau());
1624  assert(g1.nt() >= tstp);
1625  assert(g2.nt() >= tstp);
1626  if (g1.size1() == 1)
1627  return distance_norm2_dispatch<T, herm_matrix<T>, 1>(tstp, g1, g2);
1628  else
1629  return distance_norm2_dispatch<T, herm_matrix<T>, LARGESIZE>(tstp, g1,
1630  g2);
1631 }
1632 
1633 
1652 template <typename T>
1654  assert(g1.size1() == g2.size1());
1655  assert(g1.ntau() == g2.ntau());
1656  assert(g1.tstp() == tstp);
1657  assert(g2.nt() >= tstp);
1658  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1659  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1660  if (g1.size1() == 1)
1661  return distance_norm2_dispatch<T, 1>(tstp, g1_view, g2_view);
1662  else
1663  return distance_norm2_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
1664 }
1665 
1684 template <typename T>
1686  assert(g1.size1() == g2.size1());
1687  assert(g1.ntau() == g2.ntau());
1688  assert(g1.nt() >= tstp);
1689  assert(g2.tstp() == tstp);
1690  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1691  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1692  if (g1.size1() == 1)
1693  return distance_norm2_dispatch<T, 1>(tstp, g1_view, g2_view);
1694  else
1695  return distance_norm2_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
1696 }
1697 
1698 
1717 template <typename T>
1719  assert(g1.size1() == g2.size1());
1720  assert(g1.ntau() == g2.ntau());
1721  assert(g1.tstp() == tstp);
1722  assert(g2.tstp() == tstp);
1723  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1724  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1725  if (g1.size1() == 1)
1726  return distance_norm2_dispatch<T, 1>(tstp, g1_view, g2_view);
1727  else
1728  return distance_norm2_dispatch<T, LARGESIZE>(tstp, g1_view, g2_view);
1729 }
1730 
1731 
1733 template <typename T>
1735  assert(g1.size1() == g2.size1());
1736  assert(g1.ntau() == g2.ntau());
1737  assert(g1.tstp() == g2.tstp());
1738  herm_matrix_timestep_view<T> g1_view(g1.tstp(), g1);
1739  herm_matrix_timestep_view<T> g2_view(g1.tstp(), g2);
1740  if (g1.size1() == 1)
1741  return distance_norm2_dispatch<T, 1>(g1.tstp(), g1_view, g2_view);
1742  else
1743  return distance_norm2_dispatch<T, LARGESIZE>(g1.tstp(), g1_view, g2_view);
1744 }
1745 
1747 template <typename T>
1749  assert(g1.size1() == g2.size1());
1750  assert(g1.ntau() == g2.ntau());
1751  assert(g1.tstp() == g2.tstp());
1752  herm_matrix_timestep_view<T> g2_view(g1.tstp(), g2);
1753  if (g1.size1() == 1)
1754  return distance_norm2_dispatch<T, 1>(g1.tstp(), g1, g2_view);
1755  else
1756  return distance_norm2_dispatch<T, LARGESIZE>(g1.tstp(), g1, g2_view);
1757 }
1758 
1777 template <typename T>
1779  assert(g1.size1() == g2.size1());
1780  assert(g1.ntau() == g2.ntau());
1781  assert(g1.tstp() == tstp);
1782  assert(g2.nt() >= tstp);
1783  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1784  if (g1.size1() == 1)
1785  return distance_norm2_dispatch<T, 1>(tstp, g1, g2_view);
1786  else
1787  return distance_norm2_dispatch<T, LARGESIZE>(tstp, g1, g2_view);
1788 }
1789 
1808 template <typename T>
1810  assert(g1.size1() == g2.size1());
1811  assert(g1.ntau() == g2.ntau());
1812  assert(g1.nt() >= tstp);
1813  assert(g2.tstp() == tstp);
1814  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1815  if (g1.size1() == 1)
1816  return distance_norm2_dispatch<T, 1>(tstp, g1_view, g2);
1817  else
1818  return distance_norm2_dispatch<T, LARGESIZE>(tstp, g1_view, g2);
1819 }
1820 
1821 
1840 template <typename T>
1842  assert(g1.size1() == g2.size1());
1843  assert(g1.ntau() == g2.ntau());
1844  assert(g1.tstp() == tstp);
1845  assert(g2.tstp() == tstp);
1846  herm_matrix_timestep_view<T> g1_view(tstp, g1);
1847  if (g1.size1() == 1)
1848  return distance_norm2_dispatch<T, 1>(tstp, g1_view, g2);
1849  else
1850  return distance_norm2_dispatch<T, LARGESIZE>(tstp, g1_view, g2);
1851 }
1852 
1871 template <typename T>
1873  assert(g1.size1() == g2.size1());
1874  assert(g1.ntau() == g2.ntau());
1875  assert(g1.tstp() == tstp);
1876  assert(g2.tstp() == tstp);
1877  herm_matrix_timestep_view<T> g2_view(tstp, g2);
1878  if (g1.size1() == 1)
1879  return distance_norm2_dispatch<T, 1>(tstp, g1, g2_view);
1880  else
1881  return distance_norm2_dispatch<T, LARGESIZE>(tstp, g1, g2_view);
1882 }
1883 
1884 
1903 template <typename T>
1905  assert(g1.size1() == g2.size1());
1906  assert(g1.ntau() == g2.ntau());
1907  assert(g1.tstp() == tstp);
1908  assert(g2.tstp() == tstp);
1909  if (g1.size1() == 1)
1910  return distance_norm2_dispatch<T, 1>(tstp, g1, g2);
1911  else
1912  return distance_norm2_dispatch<T, LARGESIZE>(tstp, g1, g2);
1913 }
1914 
1915 
1916 // template <typename T>
1917 // T distance_norm2(int tstp, herm_matrix_timestep_view<T> &g1, herm_matrix<T> &g2) {
1918 // herm_matrix_timestep<T> tmp(tstp,g2.ntau(),g2.size1(),g2.size2(),g2.sig());
1919 // g1.get_data(tmp);
1920 // return distance_norm2(tstp,tmp,g2);
1921 // }
1922 
1923 // template <typename T>
1924 // T distance_norm2(int tstp, herm_matrix_timestep_view<T> &g1, herm_matrix_timestep<T> &g2) {
1925 // herm_matrix_timestep<T> tmp(tstp,g2.ntau(),g2.size1(),g2.size2(),g2.sig());
1926 // g1.get_data(tmp);
1927 // return distance_norm2(tmp,g2);
1928 // }
1929 
1930 template <typename T>
1931 T distance_norm2(int tstp, herm_pseudo<T> &g1, herm_pseudo<T> &g2) {
1932  if (g1.size1() == 1)
1933  return distance_norm2_dispatch<T, herm_pseudo<T>, 1>(tstp, g1, g2);
1934  else
1935  return distance_norm2_dispatch<T, herm_pseudo<T>, LARGESIZE>(tstp, g1,
1936  g2);
1937 }
1938 
1959 template <typename T>
1961  int size1=g2.size1(),size2=g2.size2(),ntau=g2.ntau();
1962  T err=0.0;
1963 
1964  assert(g1.tstp_==tstp && g1.tstp_<=g2.nt());
1965  assert(g1.ntau_==g2.ntau() && g1.size1_==g2.size1() && g1.size2_==g2.size2() );
1966  cdmatrix matg1(size1,size2);
1967  cdmatrix matg2(size1,size2);
1968 
1969  if(tstp==-1){
1970  for(int i=0;i<=ntau;i++){
1971  g1.get_mat(i,matg1);
1972  g2.get_mat(i,matg2);
1973  err+=(matg1-matg2).norm();
1974  }
1975  }else{
1976  // Ret
1977  for(int i=0;i<=tstp;i++){
1978  g1.get_ret_tstp_t(i,matg1);
1979  g2.get_ret(tstp,i,matg2);
1980  err+=(matg1-matg2).norm();
1981  }
1982  // tv
1983  for(int i=0;i<=ntau;i++){
1984  g1.get_tv(i,matg1);
1985  g2.get_tv(tstp,i,matg2);
1986  err+=(matg1-matg2).norm();
1987  }
1988  // Les
1989  for(int i=0;i<=tstp;i++){
1990  g1.get_les_t_tstp(i,matg1);
1991  g2.get_les(i,tstp,matg2);
1992  err+=(matg1-matg2).norm();
1993  }
1994  }
1995 }
1996 
2014 template <typename T>
2015 size_t mem_herm_matrix(int nt, int ntau, int size) {
2016  size_t mem = ((nt + 1) * (nt + 2) + (nt + 2) * (ntau + 1)) * size * size *
2017  sizeof(std::complex<T>);
2018  return mem;
2019 }
2020 
2037 template <typename T>
2038 size_t mem_herm_matrix_timestep(int tstp, int ntau, int size) {
2039  size_t mem = ((tstp + 2) + (ntau + 1)) * size * size *
2040  sizeof(std::complex<T>);
2041  return mem;
2042 }
2043 
2058 template <typename T>
2059 size_t mem_function(int nt, int size) {
2060  size_t mem = (nt + 2) * size * size * sizeof(std::complex<T>);
2061  return mem;
2062 }
2063 
2076 template <typename T>
2078  int ntau = G.ntau(), m, s1 = G.size1(), p1, p2;
2079  cdmatrix Gmat,Gmat_herm;
2080  for (m = 0; m <= ntau; m++) {
2081  G.get_mat(m,Gmat);
2082  Gmat_herm = 0.5*(Gmat + Gmat.adjoint());
2083  G.set_mat(m,Gmat_herm);
2084  }
2085 }
2086 
2099 template <class GG>
2101  int ntau = G.ntau(), m, s1 = G.size1(), p1, p2;
2102  cdmatrix Gmat,Gmat_herm;
2103  for (m = 0; m <= ntau; m++) {
2104  G.get_mat(m,Gmat);
2105  Gmat_herm = 0.5*(Gmat + Gmat.adjoint());
2106  G.set_mat(m,Gmat_herm);
2107  }
2108 }
2109 
2110 } // namespace cntr
2111 
2112 #endif // CNTR_UTILITIES_IMPL_H
Class herm_matrix_timestep_view serves for interfacing with class herm_matrix_timestep without copyi...
int size2_
Number of the rows in the Matrix form.
Class Integrator contains all kinds of weights for integration and differentiation of a function at ...
T correlation_energy(int tstp, herm_matrix< T > &G, herm_matrix< T > &Sigma, T beta, T h, int SolveOrder=MAX_SOLVE_ORDER)
Returns the result of the correlation energy at a given time-step from the time-diagonal convolution...
int size1(void) const
void force_matsubara_hermitian(herm_matrix< T > &G)
Force the Matsubara component of herm_matrix to be a hermitian matrix at each . ...
Class herm_matrix_timestep deals with contour objects at a particular timestep .
size_t mem_herm_matrix(int nt, int ntau, int size)
Evaluate the memory necessary for a herm_matrix.
std::complex< double > cplx
Definition: fourier.cpp:11
Integrator< T > & I(int k)
T distance_norm2_ret(int tstp, herm_matrix< T > &g1, herm_matrix< T > &g2)
Evaluate the Euclidean norm between the retarded components of two herm_matrix at a given time step...
int size1_
Number of the colums in the Matrix form.
Class function for objects with time on real axis.
void get_value(int tstp, EigenMatrix &M) const
Get matrix value of this function object at a specific time point
T distance_norm2_eigen(int tstp, herm_matrix_timestep< T > &g1, herm_matrix< T > &g2)
Evaluate the Euclidean norm between &#39;herm_matrix_timestep&#39; and &#39;herm_matrix&#39; at a given time step...
T distance_norm2(int tstp, herm_matrix< T > &g1, herm_matrix< T > &g2)
Evaluate the Euclidean norm between two herm_matrix at a given time step.
T distance_norm2_tv(int tstp, herm_matrix< T > &g1, herm_matrix< T > &g2)
Evaluate the Euclidean norm between the left-mixing components of two herm_matrces at a given time s...
void set_tk_from_mat(herm_matrix< T > &G, int SolveOrder)
Set components of herm_matrix at t=0,1,..kt from the Matsubara component.
void set_t0_from_mat(herm_matrix< T > &G)
Set t=0 components of the two-time contour object from the Matsubara component. ...
Class herm_matrix for two-time contour objects with hermitian symmetry.
T distance_norm2_les(int tstp, herm_matrix< T > &g1, herm_matrix< T > &g2)
Evaluate the Euclidean norm between the lesser components of two herm_matrces at a given time step...
size_t mem_function(int nt, int size)
Evaluate the memory necessary for a contour function.
int nt_
Maximum number of the time steps.
void set_value(int tstp, EigenMatrix &M)
Set matrix value at a specific time point
size_t mem_herm_matrix_timestep(int tstp, int ntau, int size)
Evaluate the memory necessary for a herm_matrix_timestep.