NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_vie2_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_VIE2_IMPL_H
2 #define CNTR_VIE2_IMPL_H
3 
4 #include "cntr_vie2_decl.hpp"
5 #include "fourier.hpp"
6 #include "cntr_elements.hpp"
7 #include "cntr_herm_matrix_timestep_decl.hpp"
8 #include "cntr_herm_matrix_decl.hpp"
9 #include "cntr_matsubara_impl.hpp"
10 #include "cntr_convolution_impl.hpp"
11 
12 namespace cntr {
13 
14 #define CPLX std::complex<T>
15 /* #######################################################################################
16 #
17 # [1+F]G=Q, G,Q hermitian
18 #
19 ###########################################################################################*/
20 
21 /*###########################################################################################
22 #
23 # RETARDED FUNCTION:
24 #
25 ###########################################################################################*/
27 
59 template <typename T, class GG, int SIZE1>
60 void vie2_timestep_ret(int n, GG &G, GG &F, GG &Fcc, GG &Q, integration::Integrator<T> &I,
61  T h) {
62  typedef std::complex<T> cplx;
63  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1;
64  int ss, sg, n1, l, j, i, p, q, m, j1, j2, size1 = G.size1();
65  cplx *gret, *gtemp, *mm, *qq, *qqj, *stemp, *one, cplx_i, cweight, *diffw;
66  cplx *sret, w0;
67  T weight;
68  cplx_i = cplx(0, 1);
69 
70  ss = F.element_size();
71  sg = G.element_size();
72  gtemp = new cplx[k * sg];
73  diffw = new cplx[k1 + 1];
74  qq = new cplx[(n + 1) * sg];
75  one = new cplx[sg];
76  mm = new cplx[k * k * sg];
77  stemp = new cplx[sg]; // sic
78  element_set<T, SIZE1>(size1, one, 1);
79  // check consistency:
80  assert(n > k);
81  assert(F.nt() >= n);
82  assert(G.nt() >= n);
83  assert(G.sig() == F.sig());
84  // SET ENTRIES IN TIMESTEP TO 0
85  gret = G.retptr(n, 0);
86  n1 = (n + 1) * sg;
87  for (i = 0; i < n1; i++)
88  gret[i] = 0;
89  // INITIAL VALUE t' = n
90  element_set<T, SIZE1>(size1, G.retptr(n, n), Q.retptr(n, n));
91  // START VALUES t' = n-j, j = 1...k: solve a kxk problem
92  for (i = 0; i < k * k * sg; i++)
93  mm[i] = 0;
94  for (i = 0; i < k * sg; i++)
95  qq[i] = 0;
96  for (j = 1; j <= k; j++) {
97  p = j - 1;
98  element_incr<T, SIZE1>(size1, qq + p * sg, Q.retptr(n, n - j));
99  element_incr<T, SIZE1>(size1, mm + sg * (p + k * p), one);
100  // integral
101  for (l = 0; l <= k; l++) {
102  weight = -h * I.gregory_weights(j, l);
103  if (n - l >= n - j) {
104  element_set<T, SIZE1>(
105  size1, stemp, F.retptr(n - l, n - j)); // stemp is element of type G!!
106  } else {
107  element_set<T, SIZE1>(size1, stemp, Fcc.retptr(n - j, n - l));
108  element_conj<T, SIZE1>(size1, stemp);
109  weight *= -1;
110  }
111  if (l == 0) {
112  element_incr<T, SIZE1>(size1, qq + p * sg, weight, G.retptr(n, n), stemp);
113  } else {
114  q = l - 1;
115  element_incr<T, SIZE1>(size1, mm + sg * (p * k + q), -weight, stemp);
116  }
117  }
118  }
119  element_linsolve_left<T, SIZE1>(size1, k, gtemp, mm, qq); // gtemp * mm = qq
120  for (j = 1; j <= k; j++)
121  element_set<T, SIZE1>(size1, G.retptr(n, n - j), gtemp + (j - 1) * sg);
122  // Compute the contribution to the convolution int dm G(n,m)Sigma(m,j)
123  // from G(n,m=n-k..n) to m=0...n-k, store into qq(j), without factor h!!
124  for (i = 0; i < sg * (n + 1); i++)
125  qq[i] = 0;
126  for (m = n - k; m <= n; m++) {
127  weight = -I.gregory_omega(n - m);
128  gret = G.retptr(n, m);
129  sret = F.retptr(m, 0);
130  for (j = 0; j <= n - k2; j++) {
131  element_incr<T, SIZE1>(size1, qq + j * sg, -I.gregory_weights(n - j, n - m),
132  gret, F.retptr(m, j));
133  sret += ss;
134  }
135  j1 = (n - k2 + 1 < 0 ? 0 : n - k2 + 1);
136  for (j = j1; j < n - k; j++) { // start weights
137  weight = I.gregory_weights(n - j, n - m);
138  element_incr<T, SIZE1>(size1, qq + j * sg, -I.gregory_weights(n - j, n - m),
139  gret, F.retptr(m, j));
140  sret += ss;
141  }
142  }
143  // DER REST VOM SCHUETZENFEST: t' = n-l, l = k+1 ... n:
144  w0 = -h * I.gregory_omega(0);
145  for (l = k + 1; l <= n; l++) {
146  j = n - l;
147  // set up mm and qqj for 1x1 problem:
148  qqj = qq + j * sg;
149  for (i = 0; i < sg; i++) {
150  qqj[i] *= h;
151  }
152  element_set<T, SIZE1>(size1, stemp, F.retptr(j, j));
153  for (i = 0; i < sg; i++)
154  mm[i] = one[i] - w0 * stemp[i];
155  element_incr<T, SIZE1>(size1, qqj, Q.retptr(n, j));
156  element_linsolve_left<T, SIZE1>(size1, G.retptr(n, j), mm, qqj);
157  // compute the contribution of Gret(n,j) to
158  // int dm G(n,j)Sigma(j,j2) for j2<j, store into qq(j2), without h
159  gret = G.retptr(n, j);
160  sret = F.retptr(j, 0);
161  for (j2 = 0; j2 < j - k; j2++) {
162  element_incr<T, SIZE1>(size1, qq + j2 * sg, -I.gregory_weights(n - j2, n - j),
163  gret, sret);
164  sret += ss;
165  }
166  j1 = (j - k < 0 ? 0 : j - k);
167  for (j2 = j1; j2 < j; j2++) { // start weights
168  weight = -I.gregory_omega(j - j2);
169  element_incr<T, SIZE1>(size1, qq + j2 * sg, -I.gregory_weights(n - j2, n - j),
170  gret, sret);
171  sret += ss;
172  }
173  }
174  delete[] diffw;
175  delete[] stemp;
176  delete[] qq;
177  delete[] mm;
178  delete[] gtemp;
179  delete[] one;
180  return;
181 }
183 
214 template <typename T, class GG, int SIZE1>
215 void vie2_start_ret(GG &G, GG &F, GG &Fcc, GG &Q, integration::Integrator<T> &I, T h) {
216  typedef std::complex<T> cplx;
217  int k = I.get_k(), k1 = k + 1;
218  int sg, l, n, j, p, q, i, dim, size1 = G.size1();
219  cplx ih, *gtemp, *mm, *qq, *stemp, minusi, *one;
220  T weight;
221 
222  sg = G.element_size();
223  assert(F.nt() >= k);
224  assert(G.nt() >= k);
225  assert(G.sig() == F.sig());
226 
227  // temporary storage
228  qq = new cplx[k1 * sg];
229  mm = new cplx[k1 * k1 * sg];
230  one = new cplx[sg];
231  gtemp = new cplx[k1 * sg];
232  stemp = new cplx[sg]; // sic
233 
234  element_set<T, SIZE1>(size1, one, 1.0);
235  minusi = cplx(0, -1);
236 
237  // set initial values:
238  for (n = 0; n <= k; n++)
239  element_set<T, SIZE1>(size1, G.retptr(n, n), Q.retptr(n, n));
240 
241  for (j = 0; j < k; j++) { // determine G(n,j), n=j+1..k
242  dim = k - j;
243  for (i = 0; i < k * sg; i++)
244  qq[i] = 0;
245  for (i = 0; i < k * k * sg; i++)
246  mm[i] = 0;
247  // fill the matrix mm, indices p,q
248  for (n = j + 1; n <= k; n++) {
249  p = n - (j + 1);
250  element_incr<T, SIZE1>(size1, qq + p * sg, Q.retptr(n, j));
251  for (l = 0; l <= k; l++) {
252  q = l - (j + 1);
253  if (l <= j) { // this involves G which is known and goes into qq(p)
254  element_conj<T, SIZE1>(size1, gtemp, G.retptr(j, l));
255  element_smul<T, SIZE1>(size1, gtemp, -1);
256  element_set_zero<T, SIZE1>(size1, stemp);
257  element_incr<T, SIZE1>(size1, stemp, F.retptr(n, l), gtemp);
258  for (i = 0; i < sg; i++)
259  qq[p * sg + i] += -h * I.poly_integration(j, n, l) * stemp[i];
260  } else { // this goes into mm(p,q)
261  for (i = 0; i < sg; i++)
262  mm[sg * (p * dim + q) + i] = 0.0;
263  if (l == n) {
264  for (i = 0; i < sg; i++) {
265  mm[sg * (p * dim + q) + i] += one[i];
266  }
267  }
268  weight = -h * I.poly_integration(j, n, l);
269  if (n >= l) {
270  element_set<T, SIZE1>(size1, stemp, F.retptr(n, l));
271  } else {
272  element_set<T, SIZE1>(size1, stemp, Fcc.retptr(l, n));
273  element_conj<T, SIZE1>(size1, stemp);
274  weight *= -1;
275  }
276  for (i = 0; i < sg; i++)
277  mm[sg * (p * dim + q) + i] -= weight * stemp[i];
278  }
279  }
280  }
281  element_linsolve_right<T, SIZE1>(size1, dim, gtemp, mm, qq);
282  for (n = j + 1; n <= k; n++) {
283  p = n - (j + 1);
284  element_set<T, SIZE1>(size1, G.retptr(n, j), gtemp + p * sg);
285  }
286  }
287  delete[] stemp;
288  delete[] qq;
289  delete[] mm;
290  delete[] gtemp;
291  delete[] one;
292  return;
293 }
294 /*####################################################################################
295 #
296 # MATSUBARA FUNCTION g(tau) = 1/beta sum_n e^{-iomn tau} (iomn-H-Sigma(iomn))^{-1}
297 #
298 ####################################################################################*/
300 // TODO:
301 // For Sigma=0 the boundary values seem to be a very badly approximated, although
302 // the jump=1 is accounted far correctly. Sigma!=0 somehow makes the performance
303 // better
304 template <typename T, class GG, int SIZE1>
305 void vie2_mat_fourier_dispatch(GG &G, GG &F, GG &Fcc, GG &Q, T beta, int pcf = 20, int order = 3) {
306  typedef std::complex<T> cplx;
307  cplx *fmdft, *fiomn, *qiomn, *qiomn1, *qmdft, *qmasy, *z1, *z2, *z3, *zinv, *one;
308  cplx *expfac, *xmat;
309  int ntau, m, r, p, m2, sig, size1 = G.size1(), l, sg, ss, matsub_one;
310  T dtau, omn;
311  // T arg;
312 
313  assert(G.ntau() == F.ntau());
314  sig = G.sig();
315  assert(sig == -1 || sig == 1); // Tested for both Bosons and Fermions (StrandH)
316  assert(G.sig() == F.sig());
317  sg = G.element_size();
318  ss = F.element_size();
319  ntau = G.ntau();
320  dtau = beta / ntau;
321 
322  matsub_one = (sig == -1 ? 1.0 : 0.0); // 1 for Fermions, 0 for Bosons
323 
324  // pcf=10;
325  if (ntau % 2 == 1) {
326  std::cerr << "matsubara_inverse: ntau odd" << std::endl;
327  abort();
328  }
329  m2 = ntau / 2;
330  xmat = new cplx[(ntau + 1) * sg];
331  fmdft = new cplx[(ntau + 1) * sg];
332  qmdft = new cplx[(ntau + 1) * sg];
333  expfac = new cplx[ntau + 1];
334  qmasy = new cplx[sg];
335  z1 = new cplx[sg];
336  z2 = new cplx[sg];
337  z3 = new cplx[sg];
338  zinv = new cplx[sg];
339  fiomn = new cplx[sg];
340  qiomn = new cplx[sg];
341  qiomn1 = new cplx[sg];
342  one = new cplx[sg];
343  element_set<T, SIZE1>(size1, one, 1.0);
344 
345 
346  // First order tail correction factor
347  // qmasy = -1 * ( Q(0) + Q(\beta) ) NB! Only valid for Fermions
348  // The general expression reads: qmasy = sig * Q(\beta) - Q(0)
349 
350  element_set<T, SIZE1>(size1, qmasy, Q.matptr(0));
351  element_smul<T, SIZE1>(size1, qmasy, -1.0);
352  element_set<T, SIZE1>(size1, z1, Q.matptr(ntau));
353  element_smul<T, SIZE1>(size1, z1, sig);
354  element_incr<T, SIZE1>(size1, qmasy, z1);
355 
356  // Fermion specific code
357  // element_set<T,SIZE1>(size1,qmasy,Q.matptr(0));
358  // element_incr<T,SIZE1>(size1,qmasy,Q.matptr(ntau));
359  // element_smul<T,SIZE1>(size1,qmasy,-1.0);
360 
361  set_first_order_tail<T, SIZE1>(xmat, qmasy, beta, sg, ntau, sig, size1);
362 
363  // Raw Discrete Fourier Transform of F(\tau) & G(\tau)
364  // Sig determines whether the Matsubara frequencies are Fermionic or Bosonic
365  // F(i\omega_n) computed for \omega_n \ge 0 (i.e. including zero for Bosons)
366  matsubara_dft<T, GG, SIZE1>(fmdft, F, sig);
367  matsubara_dft<T, GG, SIZE1>(qmdft, Q, sig);
368 
369  // Oversampling loop over Matsubara freqencies
370  // nomega = (2*pcf + 1) * ntau
371  for (p = -pcf; p <= pcf; p++) {
372 
373  // For Bosons we sample the middle interval including the zero frequency.
374  // (For Fermions there is no need to fiddle with the frequency interval)
375  int m_low = 0, m_high = 0;
376  if (sig == 1) {
377  m_high = (p < 0 ? 0 : 1);
378  m_low = (p > 0 ? -1 : 0);
379  }
380 
381  // Loop over matsubara frequencies
382  for (m = -m2 - m_low; m <= m2 - 1 + m_high; m++) {
383 
384  // omn=(2*(m+p*ntau)+1)*PI/(ntau*dtau); // Fermionic Matsubara freq.
385  omn = (2 * (m + p * ntau) + matsub_one) * PI /
386  (ntau * dtau); // General Matsubara freq.
387 
388  matsubara_ft<T, GG, SIZE1>(fiomn, m + p * ntau, F, fmdft, sig, beta, order);
389  matsubara_ft<T, GG, SIZE1>(qiomn, m + p * ntau, Q, qmdft, sig, beta, order);
390 
391  // This is our tail correction in frequency space
392  // z3 = -i/omn * qmasy = 1/(i*omn) * qmasy
393  element_set<T, SIZE1>(size1, z3, qmasy);
394 
395  // For Bosons we need to skip the zeroth frequency
396  if (sig == 1 && m + p * ntau == 0)
397  element_smul<T, SIZE1>(size1, z3, cplx(0.0, 0.0));
398  else
399  element_smul<T, SIZE1>(size1, z3, cplx(0.0, -1.0 / omn));
400 
401  // qiomn1 = qiomn - 1/(i*omn) * qmasy, remove tail correction from q
402  for (l = 0; l < sg; l++)
403  qiomn1[l] = qiomn[l] - z3[l];
404 
405  // z1 = 1 + f
406  for (l = 0; l < sg; l++)
407  z1[l] = fiomn[l] + one[l];
408  // z2 = f * z3 = f * qmasy/(i*omn) ?? multiply f with tail correction
409  element_mult<T, SIZE1>(size1, z2, fiomn, z3);
410 
411  // add correction to z2 = qi - z2 = q - qmasy/(i*omn) - f * qmasy/(i*omn)
412  // so z2 = q - (1 + f) * qmasy/(i*omn)
413  for (l = 0; l < sg; l++)
414  z2[l] = qiomn1[l] - z2[l];
415 
416  // zinv = 1./z1 = (1 + f)^{-1}
417  element_inverse<T, SIZE1>(size1, zinv, z1);
418  element_mult<T, SIZE1>(size1, z1, zinv, z2); // REIHENFOLGE ?? Oroa dig inte ;)
419  element_smul<T, SIZE1>(size1, z1, 1.0 / beta);
420 
421  // z1 = 1/beta * zinv * z2 = 1/beta (1+f)^{-1} * z2
422  // => z1 = 1/beta (1+f)^{-1} * [q - (1+f)*qmasy/(i*omn)]
423 
424  // So actually we could simplify this to
425  // z1 = 1/beta * [ (1+f)^{-1} * q - qmasy/(i*omn) ] ?
426 
427  for (r = 0; r <= ntau; r++) { // Loop over tau
428 
429  double arg = omn * r * dtau;
430  std::complex<double> expfactor = cplx(cos(arg), -sin(arg));
431  element_set<T, SIZE1>(size1, z2, z1);
432  element_smul<T, SIZE1>(size1, z2, expfactor);
433  element_incr<T, SIZE1>(size1, xmat + r * sg, z2);
434 
435  } // End tau loop
436  } // End matsubara freq loop
437  } // End oversampling loop
438 
439  for (r = 0; r <= ntau; r++)
440  element_set<T, SIZE1>(size1, G.matptr(r), xmat + r * sg);
441 
442  delete[] xmat;
443  delete[] fmdft;
444  delete[] qmdft;
445  delete[] expfac;
446  delete[] qmasy;
447  delete[] z1;
448  delete[] z2;
449  delete[] z3;
450  delete[] zinv;
451  delete[] qiomn;
452  delete[] qiomn1;
453  delete[] fiomn;
454  delete[] one;
455  return;
456 }
457 
459 template <typename T, class GG, int SIZE1>
460 void vie2_mat_fixpoint_dispatch(GG &G, GG &F, GG &Fcc, GG &Q, T beta,
461  integration::Integrator<T> &I, int nfixpoint, int pcf = 5,
462  int order = 3) {
463  int fixpoint, ntau = G.ntau(), size1 = G.size1(), r;
464 
465  vie2_mat_fourier_dispatch<T, GG, SIZE1>(G, F, Fcc, Q, beta, pcf, order);
466 
467  if (nfixpoint > 0) {
468  GG Qn(-1, ntau, size1, Q.sig()), dG(-1, ntau, size1, G.sig());
469  for (fixpoint = 0; fixpoint < nfixpoint; fixpoint++) {
470  convolution_matsubara(Qn, F, G, I, beta);
471  for (r = 0; r <= ntau; r++) {
472  element_incr<T, SIZE1>(size1, Qn.matptr(r), G.matptr(r));
473  element_incr<T, SIZE1>(size1, Qn.matptr(r), -1.0, Q.matptr(r));
474  }
475  vie2_mat_fourier_dispatch<T, GG, SIZE1>(dG, F, Fcc, Qn, beta, pcf, order);
476  for (r = 0; r <= ntau; r++) {
477  element_incr<T, SIZE1>(size1, G.matptr(r), -1.0, dG.matptr(r));
478  }
479  }
480  }
481  return;
482 }
483 
485 template <typename T, class GG, int SIZE1>
486 void vie2_mat_steep_dispatch(GG &G, GG &F, GG &Fcc, GG &Q, T beta,
487  integration::Integrator<T> &I, int maxiter, T maxerr, int pcf = 3,
488  int order = 3) {
489  int iter, ntau = G.ntau(), size1 = G.size1(), r;
490  std::complex<T> *temp = new std::complex<T>[size1 * size1];
491  std::complex<T> *Rcc = new std::complex<T>[size1 * size1];
492  std::complex<T> *Pcc = new std::complex<T>[size1 * size1];
493 
494  vie2_mat_fourier_dispatch<T, GG, SIZE1>(G, F, Fcc, Q, beta, pcf, order);
495 
496  if (maxiter > 0) {
497  double alpha,c1,c2,err,rsold,rsnew,weight;
498  GG R(-1, ntau, size1, Q.sig()), AR(-1, ntau, size1, Q.sig());
499  GG P(-1, ntau, size1, Q.sig()), AP(-1, ntau, size1, Q.sig());
500  GG C1(-1, ntau, size1, Q.sig());
501  GG C2(-1, ntau, size1, Q.sig());
502 
503  //initial guees for residue
504  convolution_matsubara(C1, F, G, I, beta);
505  convolution_matsubara(C2, G, Fcc, I, beta);
506  for (r = 0; r <= ntau; r++) {
507  element_set<T, SIZE1>(size1, R.matptr(r), C1.matptr(r));
508  element_incr<T, SIZE1>(size1, R.matptr(r), 1.0, C2.matptr(r));
509  element_smul<T, SIZE1>(size1, R.matptr(r), 0.5);
510  element_incr<T, SIZE1>(size1, R.matptr(r), 1.0, G.matptr(r));
511  element_incr<T, SIZE1>(size1, R.matptr(r), -1.0, Q.matptr(r));
512  }
513  R.smul(-1,-1.0);
514  for (r=0; r <= ntau; r++)
515  element_set<T, SIZE1>(size1, P.matptr(r), R.matptr(r));
516 
517 
518  // norm of residue
519  rsold = 0.0;
520  for (r = 0; r <= ntau; r++) {
521  if(r <= I.get_k()) {
522  weight = I.gregory_omega(r);
523  } else if(r >= ntau - I.get_k()) {
524  weight = I.gregory_omega(ntau-r);
525  } else {
526  weight = 1.0;
527  }
528  element_conj<T, SIZE1>(size1, Rcc, R.matptr(r));
529  element_mult<T, SIZE1>(size1, temp, Rcc, R.matptr(r));
530  for(int s=0; s<size1*size1; s++)
531  rsold += weight * (temp[s]).real();
532  }
533 
534  // if norm of residue < threshold, no iteration
535  if(rsold > maxerr) {
536  for (iter = 0; iter < maxiter; iter++) {
537 
538  // compute A.P
539  convolution_matsubara(C1, F, P, I, beta);
540  convolution_matsubara(C2, P, Fcc, I, beta);
541  for (r = 0; r <= ntau; r++) {
542  element_set<T, SIZE1>(size1, AP.matptr(r), C1.matptr(r));
543  element_incr<T, SIZE1>(size1, AP.matptr(r), 1.0, C2.matptr(r));
544  element_smul<T, SIZE1>(size1, AP.matptr(r), 0.5);
545  element_incr<T, SIZE1>(size1, AP.matptr(r), 1.0, P.matptr(r));
546  }
547 
548  // compute P.A.P
549  c2 = 0.0;
550  for (r = 0; r <= ntau; r++) {
551  if(r <= I.get_k()) {
552  weight = I.gregory_omega(r);
553  } else if(r >= ntau - I.get_k()) {
554  weight = I.gregory_omega(ntau-r);
555  } else {
556  weight = 1.0;
557  }
558  element_conj<T, SIZE1>(size1, Pcc, P.matptr(r));
559  element_mult<T, SIZE1>(size1, temp, Pcc, AP.matptr(r));
560  for(int s=0; s<size1*size1; s++)
561  c2 += weight * (temp[s]).real();
562  }
563 
564  alpha = rsold / c2;
565 
566  for (r = 0; r <= ntau; r++) {
567  element_incr<T, SIZE1>(size1, G.matptr(r), alpha, P.matptr(r));
568  element_incr<T, SIZE1>(size1, R.matptr(r), -alpha, AP.matptr(r));
569  }
570 
571  rsnew = 0.0;
572  for (r = 0; r <= ntau; r++) {
573  if(r <= I.get_k()) {
574  weight = I.gregory_omega(r);
575  } else if(r >= ntau - I.get_k()) {
576  weight = I.gregory_omega(ntau-r);
577  } else {
578  weight = 1.0;
579  }
580  element_conj<T, SIZE1>(size1, Rcc, R.matptr(r));
581  element_mult<T, SIZE1>(size1, temp, Rcc, R.matptr(r));
582  for(int s=0; s<size1*size1; s++)
583  rsnew += weight * (temp[s]).real();
584  }
585 
586  for (r = 0; r <= ntau; r++) {
587  element_set<T, SIZE1>(size1, temp, R.matptr(r));
588  element_incr<T, SIZE1>(size1, temp, rsnew/rsold, P.matptr(r));
589  element_set<T, SIZE1>(size1, P.matptr(r), temp);
590  }
591 
592  rsold = rsnew;
593 
594  if (sqrt(rsold) < maxerr) break;
595  }
596  }
597  delete[] temp;
598  delete[] Rcc;
599  delete[] Pcc;
600  return;
601  }
602 }
603 
604 /*###########################################################################################
605 #
606 # tv-FUNCTION:
607 #
608 ###########################################################################################*/
610 
644 template <typename T, class GG, int SIZE1>
645 void vie2_timestep_tv(int n, GG &G, GG &F, GG &Fcc, GG &Q, integration::Integrator<T> &I,
646  T beta, T h) {
647  typedef std::complex<T> cplx;
648  int k = I.get_k();
649  int sg, n1, l, j, ntau, size1 = G.size1();
650  cplx *gtv, cweight, ih, *mm, *qq, *stemp, minusi, *one;
651  T weight;
652 
653  sg = G.element_size();
654  ntau = G.ntau();
655  // check consistency: (more assertations follow in convolution)
656  assert(n > k);
657  assert(F.nt() >= n);
658  assert(G.nt() >= n);
659  assert(G.sig() == F.sig());
660 
661  one = new cplx[sg];
662  qq = new cplx[sg];
663  mm = new cplx[sg];
664  stemp = new cplx[sg]; // sic
665  element_set<T, SIZE1>(size1, one, 1.0);
666  // SET ENTRIES IN TIMESTEP(TV) TO 0
667  gtv = G.tvptr(n, 0);
668  n1 = (ntau + 1) * sg;
669  for (l = 0; l < n1; l++)
670  gtv[l] = 0;
671  // CONVOLUTION SIGMA*G: ---> Gtv(n,m)
672  convolution_timestep_tv<T, herm_matrix<T>, SIZE1>(n, G, F, Fcc, G, G, I, beta,
673  h); // note: this sets only tv
674  // Now solve
675  // [ 1 - h w(n,0) Sigma(n,n) ] G(n,m) = Q(m),
676  // where Q is initially stored in G(n,m)
677  element_set<T, SIZE1>(size1, stemp, F.retptr(n, n));
678  weight = h * I.gregory_weights(n, 0);
679  for (l = 0; l < sg; l++)
680  mm[l] = one[l] + weight * stemp[l];
681  for (j = 0; j <= ntau; j++) {
682  for (l = 0; l < sg; l++)
683  qq[l] = -G.tvptr(n, j)[l] + Q.tvptr(n, j)[l];
684  element_linsolve_right<T, SIZE1>(size1, G.tvptr(n, j), mm, qq);
685  }
686  delete[] stemp;
687  delete[] qq;
688  delete[] mm;
689  delete[] one;
690  return;
691 }
693 
726 template <typename T, class GG, int SIZE1>
727 void vie2_start_tv(GG &G, GG &F, GG &Fcc, GG &Q, integration::Integrator<T> &I, T beta,
728  T h) {
729  typedef std::complex<T> cplx;
730  int k = I.get_k(), k1 = k + 1;
731  int sg, l, m, j, ntau, p, q, n, sig, size1 = G.size1();
732  cplx cweight, *mm, *qq, *stemp, *one, *gtemp;
733  cplx cplx_i = cplx(0.0, 1.0);
734  T dtau;
735  sg = G.element_size();
736  ntau = G.ntau();
737  dtau = beta / ntau;
738  sig = G.sig();
739  // check consistency: (more assertations follow in convolution)
740  assert(F.nt() >= k);
741  assert(G.nt() >= k);
742  assert(F.ntau() == ntau);
743  assert(G.sig() == F.sig());
744 
745  one = new cplx[sg];
746  qq = new cplx[k1 * sg];
747  mm = new cplx[k1 * k1 * sg];
748  stemp = new cplx[sg]; // sic
749  gtemp = new cplx[k1 * sg];
750  element_set<T, SIZE1>(size1, one, 1.0);
751 
752  // CONVOLUTION -i int dtau Sigma^tv(n,tau)G^mat(tau,m) ---> Gtv(n,m)
753  for (n = 0; n <= k; n++) {
754  for (m = 0; m <= ntau; m++) {
755  matsubara_integral_2<T, SIZE1>(size1, m, ntau, gtemp, F.tvptr(n, 0), G.matptr(0),
756  I, G.sig());
757  for (l = 0; l < sg; l++)
758  G.tvptr(n, m)[l] = -dtau * gtemp[l];
759  }
760  }
761  // loop over m --> determine G^tv(n,m) for n=0...k
762  for (m = 0; m <= ntau; m++) {
763  for (l = 0; l < k1 * k1 * sg; l++)
764  mm[l] = 0;
765  for (l = 0; l < k1 * sg; l++)
766  qq[l] = 0;
767  // derive linear equations
768  // mm(p,q)*G(q)=Q(p) for p=n-1=0...k-1, q=n-1=0...k
769  // G(p)=G(p+1,m)
770  for (n = 0; n <= k; n++) {
771  // 1 -- goes into m(p,p)
772  element_incr<T, SIZE1>(size1, mm + sg * (n * k1 + n), one);
773  // integral 0..n
774  for (j = 0; j <= k; j++) {
775  cweight = h * I.gregory_weights(n, j);
776  if (n >= j) {
777  element_set<T, SIZE1>(size1, stemp, F.retptr(n, j));
778  } else {
779  element_set<T, SIZE1>(size1, stemp, Fcc.retptr(j, n));
780  element_conj<T, SIZE1>(size1, stemp);
781  element_smul<T, SIZE1>(size1, stemp, -1);
782  }
783  for (l = 0; l < sg; l++)
784  mm[sg * (n * k1 + j) + l] += cweight * stemp[l];
785  }
786  // integral Sigmatv*Gmat --> take from Gtv(n,m), write into qq
787  element_incr<T, SIZE1>(size1, qq + n * sg, G.tvptr(n, m));
788  element_incr<T, SIZE1>(size1, qq + n * sg, Q.tvptr(n, m));
789  }
790  element_linsolve_right<T, SIZE1>(size1, k1, gtemp, mm,
791  qq); // solve kXk problem mm*gtemp=qq
792  // write elements into Gtv
793  for (n = 0; n <= k; n++)
794  element_set<T, SIZE1>(size1, G.tvptr(n, m), gtemp + n * sg);
795  }
796  delete[] stemp;
797  delete[] gtemp;
798  delete[] qq;
799  delete[] mm;
800  delete[] one;
801  return;
802 }
803 /*###########################################################################################
804 #
805 # les-FUNCTION: use id_t G(t,nh)= xxx , t=0...n
806 # note: G(t,n) is not continuous in memory!!!
807 # for n<k: start routine
808 #
809 ###########################################################################################*/
811 
846 template <typename T, class GG, int SIZE1>
847 void vie2_timestep_les(int n, GG &G, GG &F, GG &Fcc, GG &Q, integration::Integrator<T> &I,
848  T beta, T h) {
849  typedef std::complex<T> cplx;
850  int k = I.get_k(), k1 = k + 1;
851  int sg, n1, l, m, j, ntau, p, q, sig, size1 = G.size1();
852  cplx *gles, cweight, ih, *mm, *qq, *stemp, cplx_i = cplx(0, 1), *one, *gtemp, *qtemp;
853 
854  n1 = (n > k ? n : k);
855  sg = G.element_size();
856  ntau = G.ntau();
857  sig = G.sig();
858  // check consistency: (more assertations follow in convolution)
859  assert(F.nt() >= n1);
860  assert(G.nt() >= n1);
861  assert(G.sig() == F.sig());
862 
863  one = new cplx[sg];
864  qq = new cplx[k1 * sg];
865  mm = new cplx[k1 * k1 * sg];
866  gtemp = new cplx[sg];
867  qtemp = new cplx[sg];
868  stemp = new cplx[sg]; // sic
869  gles = new cplx[(n1 + 1) * sg];
870  for (j = 0; j <= n1; j++)
871  element_set_zero<T, SIZE1>(size1, gles + j * sg);
872  element_set<T, SIZE1>(size1, one, 1.0);
873 
874  // CONVOLUTION SIGMA*G: ---> G^les(j,n)
875  // Note: this is only the tv*vt + les*adv part, Gles is not adressed
876  // Note: gles is determioned for j=0...max(k,n)!!!
877  convolution_timestep_les_tvvt<T, GG, SIZE1>(n, gles, G, F, Fcc, G, G, I, beta,
878  h); // sig needed!!!
879  convolution_timestep_les_lesadv<T, GG, SIZE1>(n, gles, G, F, Fcc, G, G, I, beta, h);
880  for (j = 0; j <= n1; j++)
881  element_smul<T, SIZE1>(size1, gles + j * sg, -1.0);
882 
883  // Start for integrodifferential equation: j=1...k
884  // .... the usual mess:
885  // derive linear equations
886  // mm(p,q)*G(q)=Q(p) for p=j-1=0...k-1, q=j-1=0...k
887  // G(p)=G(p,m)
888  for (l = 0; l < k1 * k1 * sg; l++)
889  mm[l] = 0;
890  for (l = 0; l < k1 * sg; l++)
891  qq[l] = 0;
892  for (j = 0; j <= k; j++) {
893  p = j;
894  if (j <= n) {
895  element_incr<T, SIZE1>(size1, qq + p * sg, Q.lesptr(j, n));
896  } else {
897  element_conj<T, SIZE1>(size1, qtemp, Q.lesptr(n, j));
898  element_smul<T, SIZE1>(size1, qtemp, -1.0);
899  element_incr<T, SIZE1>(size1, qq + p * sg, qtemp);
900  }
901  // integral Sigmatv*Gvt + Sigma^les*G^adv --> take from gles(j)
902  // after that, gles can be overwritten
903  element_incr<T, SIZE1>(size1, qq + p * sg, gles + j * sg);
904  // -- goes into m(p,p)
905  element_incr<T, SIZE1>(size1, mm + sg * (p * k1 + p), one);
906  // integral Sigma^ret(j,m)G^les(m,n)
907  for (m = 0; m <= k; m++) {
908  cweight = -h * I.gregory_weights(j, m);
909  if (0 && m == 0) { // goes into qq(p)
910  element_incr<T, SIZE1>(size1, qq + p * sg, cweight, F.retptr(j, 0),
911  gles + 0 * sg);
912  } else { // goes into mm(p,q)
913  // q=m-1;
914  q = m;
915  if (j >= m) {
916  element_set<T, SIZE1>(size1, stemp, F.retptr(j, m));
917  } else {
918  element_set<T, SIZE1>(size1, stemp, Fcc.retptr(m, j));
919  element_conj<T, SIZE1>(size1, stemp);
920  element_smul<T, SIZE1>(size1, stemp, -1);
921  }
922  for (l = 0; l < sg; l++)
923  mm[sg * (p * k1 + q) + l] += -cweight * stemp[l];
924  }
925  }
926  }
927  element_linsolve_right<T, SIZE1>(size1, k1, gles, mm,
928  qq); // solve k1Xk1 problem mm*gtemp=qq
929  // integrodifferential equation k+1...n
930  for (j = k + 1; j <= n; j++) {
931  element_set_zero<T, SIZE1>(size1, mm);
932  // CONTRIBUTION FROM INTEGRAL tv*vt+les*adv
933  element_set<T, SIZE1>(size1, qq, gles + j * sg); // now gles(j) may be overwritten
934  element_incr<T, SIZE1>(size1, qq, Q.lesptr(j, n));
935  element_incr<T, SIZE1>(size1, mm, one);
936  // CONTRIBUTION FROM INTEGRAL Sigma^ret(j,m)*G^les(m,n)
937  element_set<T, SIZE1>(size1, stemp, F.retptr(j, j));
938  for (l = 0; l < sg; l++)
939  mm[l] += h * stemp[l] * I.gregory_weights(j, j);
940  for (m = 0; m < j; m++) {
941  cweight = -h * I.gregory_weights(j, m);
942  element_incr<T, SIZE1>(size1, qq, cweight, F.retptr(j, m), gles + m * sg);
943  }
944  element_linsolve_right<T, SIZE1>(size1, gles + j * sg, mm, qq);
945  }
946  // write elements into Gles
947  for (j = 0; j <= n; j++)
948  element_set<T, SIZE1>(size1, G.lesptr(j, n), gles + j * sg);
949  delete[] stemp;
950  delete[] gtemp;
951  delete[] qtemp;
952  delete[] gles;
953  delete[] qq;
954  delete[] mm;
955  delete[] one;
956  return;
957 }
959 
992 template <typename T, class GG, int SIZE1>
993 void vie2_start_les(GG &G, GG &F, GG &Fcc, GG &Q, integration::Integrator<T> &I, T beta,
994  T h) {
995  int k = I.get_k(), n;
996  for (n = 0; n <= k; n++)
997  vie2_timestep_les<T, GG, SIZE1>(n, G, F, Fcc, Q, I, beta, h);
998  return;
999 }
1001 
1029 template <typename T>
1030 void vie2_mat_fourier(herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc, herm_matrix<T> &Q,
1031  T beta, int order) {
1032  int size1 = G.size1(), pcf = 20;
1033  assert(G.size1() == F.size1());
1034  assert(G.ntau() == F.ntau());
1035  switch (size1) {
1036  case 1:
1037  vie2_mat_fourier_dispatch<T, herm_matrix<T>, 1>(G, F, Fcc, Q, beta, pcf, order);
1038  break;
1039  case 2:
1040  vie2_mat_fourier_dispatch<T, herm_matrix<T>, 2>(G, F, Fcc, Q, beta, pcf, order);
1041  break;
1042  case 3:
1043  vie2_mat_fourier_dispatch<T, herm_matrix<T>, 3>(G, F, Fcc, Q, beta, pcf, order);
1044  break;
1045  case 4:
1046  vie2_mat_fourier_dispatch<T, herm_matrix<T>, 4>(G, F, Fcc, Q, beta, pcf, order);
1047  break;
1048  case 5:
1049  vie2_mat_fourier_dispatch<T, herm_matrix<T>, 5>(G, F, Fcc, Q, beta, pcf, order);
1050  break;
1051  case 8:
1052  vie2_mat_fourier_dispatch<T, herm_matrix<T>, 8>(G, F, Fcc, Q, beta, pcf, order);
1053  break;
1054  default:
1055  vie2_mat_fourier_dispatch<T, herm_matrix<T>, LARGESIZE>(G, F, Fcc, Q, beta, pcf, order);
1056  break;
1057  }
1058 }
1060 
1090 template <typename T>
1091 void vie2_mat_fixpoint(herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc,
1092  herm_matrix<T> &Q, T beta, integration::Integrator<T> &I, int nfixpoint) {
1093  int size1 = G.size1(), pcf = 5, order = 3;
1094  assert(G.size1() == F.size1());
1095  assert(G.ntau() == F.ntau());
1096  switch (size1) {
1097  case 1:
1098  vie2_mat_fixpoint_dispatch<T, herm_matrix<T>, 1>(G, F, Fcc, Q, beta, I, nfixpoint, pcf,
1099  order);
1100  break;
1101  case 2:
1102  vie2_mat_fixpoint_dispatch<T, herm_matrix<T>, 2>(G, F, Fcc, Q, beta, I, nfixpoint, pcf,
1103  order);
1104  break;
1105  case 3:
1106  vie2_mat_fixpoint_dispatch<T, herm_matrix<T>, 3>(G, F, Fcc, Q, beta, I, nfixpoint, pcf,
1107  order);
1108  break;
1109  case 4:
1110  vie2_mat_fixpoint_dispatch<T, herm_matrix<T>, 4>(G, F, Fcc, Q, beta, I, nfixpoint, pcf,
1111  order);
1112  break;
1113  case 5:
1114  vie2_mat_fixpoint_dispatch<T, herm_matrix<T>, 5>(G, F, Fcc, Q, beta, I, nfixpoint, pcf,
1115  order);
1116  break;
1117  case 8:
1118  vie2_mat_fixpoint_dispatch<T, herm_matrix<T>, 8>(G, F, Fcc, Q, beta, I, nfixpoint, pcf,
1119  order);
1120  break;
1121  default:
1122  vie2_mat_fixpoint_dispatch<T, herm_matrix<T>, LARGESIZE>(G, F, Fcc, Q, beta, I, nfixpoint,
1123  pcf, order);
1124  break;
1125  }
1126 }
1128 
1160 template <typename T>
1161 void vie2_mat_steep(herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc,
1162  herm_matrix<T> &Q, T beta, integration::Integrator<T> &I, int maxiter, T tol) {
1163  int size1 = G.size1(), pcf = 3, order = 3;
1164  assert(G.size1() == F.size1());
1165  assert(G.ntau() == F.ntau());
1166  switch (size1) {
1167  case 1:
1168  vie2_mat_steep_dispatch<T, herm_matrix<T>, 1>(G, F, Fcc, Q, beta, I, maxiter, tol, pcf,
1169  order);
1170  break;
1171  case 2:
1172  vie2_mat_steep_dispatch<T, herm_matrix<T>, 2>(G, F, Fcc, Q, beta, I, maxiter, tol, pcf,
1173  order);
1174  break;
1175  case 3:
1176  vie2_mat_steep_dispatch<T, herm_matrix<T>, 3>(G, F, Fcc, Q, beta, I, maxiter, tol, pcf,
1177  order);
1178  break;
1179  case 4:
1180  vie2_mat_steep_dispatch<T, herm_matrix<T>, 4>(G, F, Fcc, Q, beta, I, maxiter, tol, pcf,
1181  order);
1182  break;
1183  case 5:
1184  vie2_mat_steep_dispatch<T, herm_matrix<T>, 5>(G, F, Fcc, Q, beta, I, maxiter, tol, pcf,
1185  order);
1186  break;
1187  case 8:
1188  vie2_mat_steep_dispatch<T, herm_matrix<T>, 8>(G, F, Fcc, Q, beta, I, maxiter, tol, pcf,
1189  order);
1190  break;
1191  default:
1192  vie2_mat_steep_dispatch<T, herm_matrix<T>, LARGESIZE>(G, F, Fcc, Q, beta, I, maxiter, tol,
1193  pcf, order);
1194  break;
1195  }
1196 }
1198 
1228 template <typename T>
1229 void vie2_mat(herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc,
1230  herm_matrix<T> &Q, T beta, integration::Integrator<T> &I, const int method){
1231 
1232  const int fourier_order=3;
1233  int maxiter;
1234  T tol = 1.0e-12;
1235 
1236  switch(method) {
1237  case CNTR_MAT_FOURIER:
1238  vie2_mat_fourier(G, F, Fcc, Q, beta, fourier_order);
1239  break;
1240  case CNTR_MAT_CG:
1241  maxiter = 40;
1242  vie2_mat_steep(G, F, Fcc, Q, beta, I, maxiter, tol);
1243  break;
1244  default:
1245  maxiter = 6;
1246  vie2_mat_fixpoint(G, F, Fcc, Q, beta, I, maxiter);
1247  break;
1248  }
1249 
1250 }
1251 
1282 template <typename T>
1283 void vie2_mat(herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc,
1284  herm_matrix<T> &Q, T beta, const int SolveOrder, const int method){
1285 
1286  const int fourier_order=3;
1287  int maxiter;
1288  T tol = 1.0e-12;
1289 
1290  switch(method) {
1291  case CNTR_MAT_FOURIER:
1292  vie2_mat_fourier(G, F, Fcc, Q, beta, fourier_order);
1293  break;
1294  case CNTR_MAT_CG:
1295  maxiter = 40;
1296  vie2_mat_steep(G, F, Fcc, Q, beta, integration::I<T>(SolveOrder), maxiter, tol);
1297  break;
1298  default:
1299  maxiter = 6;
1300  vie2_mat_fixpoint(G, F, Fcc, Q, beta, integration::I<T>(SolveOrder), maxiter);
1301  break;
1302  }
1303 
1304 }
1306 
1337 template <typename T>
1338 void vie2_start(herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc, herm_matrix<T> &Q,
1339  integration::Integrator<T> &I, T beta, T h) {
1340  int size1 = G.size1(), k = I.k();
1341  assert(G.size1() == F.size1());
1342  assert(G.ntau() == F.ntau());
1343  assert(G.nt() >= k);
1344  assert(F.nt() >= k);
1345  assert(G.size1() == Fcc.size1());
1346  assert(G.ntau() == Fcc.ntau());
1347  assert(Fcc.nt() >= k);
1348 
1349  switch (size1) {
1350  case 1:
1351  vie2_start_ret<T, herm_matrix<T>, 1>(G, F, Fcc, Q, I, h);
1352  vie2_start_tv<T, herm_matrix<T>, 1>(G, F, Fcc, Q, I, beta, h);
1353  vie2_start_les<T, herm_matrix<T>, 1>(G, F, Fcc, Q, I, beta, h);
1354  break;
1355  case 2:
1356  vie2_start_ret<T, herm_matrix<T>, 2>(G, F, Fcc, Q, I, h);
1357  vie2_start_tv<T, herm_matrix<T>, 2>(G, F, Fcc, Q, I, beta, h);
1358  vie2_start_les<T, herm_matrix<T>, 2>(G, F, Fcc, Q, I, beta, h);
1359  break;
1360  case 3:
1361  vie2_start_ret<T, herm_matrix<T>, 3>(G, F, Fcc, Q, I, h);
1362  vie2_start_tv<T, herm_matrix<T>, 3>(G, F, Fcc, Q, I, beta, h);
1363  vie2_start_les<T, herm_matrix<T>, 3>(G, F, Fcc, Q, I, beta, h);
1364  break;
1365  case 4:
1366  vie2_start_ret<T, herm_matrix<T>, 4>(G, F, Fcc, Q, I, h);
1367  vie2_start_tv<T, herm_matrix<T>, 4>(G, F, Fcc, Q, I, beta, h);
1368  vie2_start_les<T, herm_matrix<T>, 4>(G, F, Fcc, Q, I, beta, h);
1369  break;
1370  case 5:
1371  vie2_start_ret<T, herm_matrix<T>, 5>(G, F, Fcc, Q, I, h);
1372  vie2_start_tv<T, herm_matrix<T>, 5>(G, F, Fcc, Q, I, beta, h);
1373  vie2_start_les<T, herm_matrix<T>, 5>(G, F, Fcc, Q, I, beta, h);
1374  break;
1375  case 8:
1376  vie2_start_ret<T, herm_matrix<T>, 8>(G, F, Fcc, Q, I, h);
1377  vie2_start_tv<T, herm_matrix<T>, 8>(G, F, Fcc, Q, I, beta, h);
1378  vie2_start_les<T, herm_matrix<T>, 8>(G, F, Fcc, Q, I, beta, h);
1379  break;
1380  default:
1381  vie2_start_ret<T, herm_matrix<T>, LARGESIZE>(G, F, Fcc, Q, I, h);
1382  vie2_start_tv<T, herm_matrix<T>, LARGESIZE>(G, F, Fcc, Q, I, beta, h);
1383  vie2_start_les<T, herm_matrix<T>, LARGESIZE>(G, F, Fcc, Q, I, beta, h);
1384  break;
1385  }
1386 }
1387 
1419 template <typename T>
1420 void vie2_start(herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc, herm_matrix<T> &Q,
1421  T beta, T h, const int SolveOrder) {
1422  int size1 = G.size1();
1423  assert(G.size1() == F.size1());
1424  assert(G.ntau() == F.ntau());
1425  assert(G.nt() >= SolveOrder);
1426  assert(F.nt() >= SolveOrder);
1427  assert(G.size1() == Fcc.size1());
1428  assert(G.ntau() == Fcc.ntau());
1429  assert(Fcc.nt() >= SolveOrder);
1430 
1431  switch (size1) {
1432  case 1:
1433  vie2_start_ret<T, herm_matrix<T>, 1>(G, F, Fcc, Q, integration::I<T>(SolveOrder), h);
1434  vie2_start_tv<T, herm_matrix<T>, 1>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1435  vie2_start_les<T, herm_matrix<T>, 1>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1436  break;
1437  case 2:
1438  vie2_start_ret<T, herm_matrix<T>, 2>(G, F, Fcc, Q, integration::I<T>(SolveOrder), h);
1439  vie2_start_tv<T, herm_matrix<T>, 2>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1440  vie2_start_les<T, herm_matrix<T>, 2>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1441  break;
1442  case 3:
1443  vie2_start_ret<T, herm_matrix<T>, 3>(G, F, Fcc, Q, integration::I<T>(SolveOrder), h);
1444  vie2_start_tv<T, herm_matrix<T>, 3>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1445  vie2_start_les<T, herm_matrix<T>, 3>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1446  break;
1447  case 4:
1448  vie2_start_ret<T, herm_matrix<T>, 4>(G, F, Fcc, Q, integration::I<T>(SolveOrder), h);
1449  vie2_start_tv<T, herm_matrix<T>, 4>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1450  vie2_start_les<T, herm_matrix<T>, 4>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1451  break;
1452  case 5:
1453  vie2_start_ret<T, herm_matrix<T>, 5>(G, F, Fcc, Q, integration::I<T>(SolveOrder), h);
1454  vie2_start_tv<T, herm_matrix<T>, 5>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1455  vie2_start_les<T, herm_matrix<T>, 5>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1456  break;
1457  case 8:
1458  vie2_start_ret<T, herm_matrix<T>, 8>(G, F, Fcc, Q, integration::I<T>(SolveOrder), h);
1459  vie2_start_tv<T, herm_matrix<T>, 8>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1460  vie2_start_les<T, herm_matrix<T>, 8>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1461  break;
1462  default:
1463  vie2_start_ret<T, herm_matrix<T>, LARGESIZE>(G, F, Fcc, Q, integration::I<T>(SolveOrder), h);
1464  vie2_start_tv<T, herm_matrix<T>, LARGESIZE>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1465  vie2_start_les<T, herm_matrix<T>, LARGESIZE>(G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1466  break;
1467  }
1468 }
1470 
1504 template <typename T>
1505 void vie2_timestep(int n, herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc,
1506  herm_matrix<T> &Q, integration::Integrator<T> &I, T beta, T h, const int matsubara_method) {
1507  int size1 = G.size1(), k = I.k();
1508  assert(G.size1() == F.size1());
1509  assert(G.size1() == Fcc.size1());
1510  assert(G.ntau() == F.ntau());
1511  assert(G.ntau() == Fcc.ntau());
1512  assert(G.nt() >= n);
1513  assert(F.nt() >= n);
1514  assert(Fcc.nt() >= n);
1515 
1516  if (n==-1){
1517  vie2_mat(G, F, Fcc, Q, beta, I, matsubara_method);
1518  }else if(n<=k){
1519  vie2_start(G,F,Fcc,Q,integration::I<T>(n),beta,h);
1520  }else{
1521  switch (size1) {
1522  case 1:
1523  vie2_timestep_ret<T, herm_matrix<T>, 1>(n, G, Fcc, F, Q, I, h);
1524  vie2_timestep_tv<T, herm_matrix<T>, 1>(n, G, F, Fcc, Q, I, beta, h);
1525  vie2_timestep_les<T, herm_matrix<T>, 1>(n, G, F, Fcc, Q, I, beta, h);
1526  break;
1527  case 2:
1528  vie2_timestep_ret<T, herm_matrix<T>, 2>(n, G, Fcc, F, Q, I, h);
1529  vie2_timestep_tv<T, herm_matrix<T>, 2>(n, G, F, Fcc, Q, I, beta, h);
1530  vie2_timestep_les<T, herm_matrix<T>, 2>(n, G, F, Fcc, Q, I, beta, h);
1531  break;
1532  case 3:
1533  vie2_timestep_ret<T, herm_matrix<T>, 3>(n, G, Fcc, F, Q, I, h);
1534  vie2_timestep_tv<T, herm_matrix<T>, 3>(n, G, F, Fcc, Q, I, beta, h);
1535  vie2_timestep_les<T, herm_matrix<T>, 3>(n, G, F, Fcc, Q, I, beta, h);
1536  break;
1537  case 4:
1538  vie2_timestep_ret<T, herm_matrix<T>, 4>(n, G, Fcc, F, Q, I, h);
1539  vie2_timestep_tv<T, herm_matrix<T>, 4>(n, G, F, Fcc, Q, I, beta, h);
1540  vie2_timestep_les<T, herm_matrix<T>, 4>(n, G, F, Fcc, Q, I, beta, h);
1541  break;
1542  case 5:
1543  vie2_timestep_ret<T, herm_matrix<T>, 5>(n, G, Fcc, F, Q, I, h);
1544  vie2_timestep_tv<T, herm_matrix<T>, 5>(n, G, F, Fcc, Q, I, beta, h);
1545  vie2_timestep_les<T, herm_matrix<T>, 5>(n, G, F, Fcc, Q, I, beta, h);
1546  break;
1547  case 8:
1548  vie2_timestep_ret<T, herm_matrix<T>, 8>(n, G, Fcc, F, Q, I, h);
1549  vie2_timestep_tv<T, herm_matrix<T>, 8>(n, G, F, Fcc, Q, I, beta, h);
1550  vie2_timestep_les<T, herm_matrix<T>, 8>(n, G, F, Fcc, Q, I, beta, h);
1551  break;
1552  default:
1553  vie2_timestep_ret<T, herm_matrix<T>, LARGESIZE>(n, G, Fcc, F, Q, I, h);
1554  vie2_timestep_tv<T, herm_matrix<T>, LARGESIZE>(n, G, F, Fcc, Q, I, beta, h);
1555  vie2_timestep_les<T, herm_matrix<T>, LARGESIZE>(n, G, F, Fcc, Q, I, beta, h);
1556  break;
1557  }
1558  }
1559 }
1560 
1561 
1596 template <typename T>
1597 void vie2_timestep(int n, herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc,
1598  herm_matrix<T> &Q, T beta, T h, const int SolveOrder, const int matsubara_method) {
1599  int size1 = G.size1();
1600  assert(G.size1() == F.size1());
1601  assert(G.size1() == Fcc.size1());
1602  assert(G.ntau() == F.ntau());
1603  assert(G.ntau() == Fcc.ntau());
1604  assert(G.nt() >= n);
1605  assert(F.nt() >= n);
1606  assert(Fcc.nt() >= n);
1607 
1608  if (n==-1){
1609  vie2_mat(G, F, Fcc, Q, beta, integration::I<T>(SolveOrder), matsubara_method);
1610  }else if(n<=SolveOrder){
1611  vie2_start(G,F,Fcc,Q,integration::I<T>(n),beta,h);
1612  }else{
1613  switch (size1) {
1614  case 1:
1615  vie2_timestep_ret<T, herm_matrix<T>, 1>(n, G, Fcc, F, Q, integration::I<T>(SolveOrder), h);
1616  vie2_timestep_tv<T, herm_matrix<T>, 1>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1617  vie2_timestep_les<T, herm_matrix<T>, 1>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1618  break;
1619  case 2:
1620  vie2_timestep_ret<T, herm_matrix<T>, 2>(n, G, Fcc, F, Q, integration::I<T>(SolveOrder), h);
1621  vie2_timestep_tv<T, herm_matrix<T>, 2>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1622  vie2_timestep_les<T, herm_matrix<T>, 2>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1623  break;
1624  case 3:
1625  vie2_timestep_ret<T, herm_matrix<T>, 3>(n, G, Fcc, F, Q, integration::I<T>(SolveOrder), h);
1626  vie2_timestep_tv<T, herm_matrix<T>, 3>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1627  vie2_timestep_les<T, herm_matrix<T>, 3>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1628  break;
1629  case 4:
1630  vie2_timestep_ret<T, herm_matrix<T>, 4>(n, G, Fcc, F, Q, integration::I<T>(SolveOrder), h);
1631  vie2_timestep_tv<T, herm_matrix<T>, 4>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1632  vie2_timestep_les<T, herm_matrix<T>, 4>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1633  break;
1634  case 5:
1635  vie2_timestep_ret<T, herm_matrix<T>, 5>(n, G, Fcc, F, Q, integration::I<T>(SolveOrder), h);
1636  vie2_timestep_tv<T, herm_matrix<T>, 5>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1637  vie2_timestep_les<T, herm_matrix<T>, 5>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1638  break;
1639  case 8:
1640  vie2_timestep_ret<T, herm_matrix<T>, 8>(n, G, Fcc, F, Q, integration::I<T>(SolveOrder), h);
1641  vie2_timestep_tv<T, herm_matrix<T>, 8>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1642  vie2_timestep_les<T, herm_matrix<T>, 8>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1643  break;
1644  default:
1645  vie2_timestep_ret<T, herm_matrix<T>, LARGESIZE>(n, G, Fcc, F, Q, integration::I<T>(SolveOrder), h);
1646  vie2_timestep_tv<T, herm_matrix<T>, LARGESIZE>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1647  vie2_timestep_les<T, herm_matrix<T>, LARGESIZE>(n, G, F, Fcc, Q, integration::I<T>(SolveOrder), beta, h);
1648  break;
1649  }
1650  }
1651 }
1653 
1685 template <typename T>
1686 void vie2(herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc, herm_matrix<T> &Q,
1687  integration::Integrator<T> &I, T beta, T h, const int matsubara_method) {
1688  int tstp, k = I.k();
1689  vie2_mat(G, F, Fcc, Q, beta, I, matsubara_method);
1690  if (G.nt() >= 0)
1691  vie2_start(G, F, Fcc, Q, I, beta, h);
1692  for (tstp = k + 1; tstp <= G.nt(); tstp++)
1693  vie2_timestep(tstp, G, F, Fcc, Q, I, beta, h);
1694 }
1695 
1728 template <typename T>
1729 void vie2(herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc, herm_matrix<T> &Q,
1730  T beta, T h, const int SolveOrder, const int matsubara_method) {
1731  int tstp;
1732  vie2_mat(G, F, Fcc, Q, beta, SolveOrder, matsubara_method);
1733  if (G.nt() >= 0)
1734  vie2_start(G, F, Fcc, Q, beta, h, SolveOrder);
1735  for (tstp = SolveOrder + 1; tstp <= G.nt(); tstp++)
1736  vie2_timestep(tstp, G, F, Fcc, Q, beta, h, SolveOrder);
1737 }
1738 
1778 template < typename T>
1779 void vie2_timestep_sin(int n,herm_matrix<T> &G,function<T> &Gsin,herm_matrix<T> &F,herm_matrix<T> &Fcc, function<T> &Fsin ,herm_matrix<T> &Q,function<T> &Qsin,T beta,T h,int SolveOrder){
1780 
1781  int n1=(n<=SolveOrder && n>=0 ? 0 : n);
1782  int n2=(n<=SolveOrder && n>=0 ? SolveOrder : n);
1783  int nt=F.nt(),ntau=F.ntau(),size1=F.size1(),sig=F.sig();
1784 
1785  assert(n >= -1);
1786  assert(ntau > 0);
1787  assert(SolveOrder > 0 && SolveOrder <= 5);
1788  assert(SolveOrder <= 2 * ntau + 2);
1789  assert(ntau == Fcc.ntau());
1790  assert(ntau == F.ntau());
1791  assert(ntau == Q.ntau());
1792  assert(F.sig()== G.sig());
1793  assert(Fcc.sig()== G.sig());
1794  assert(Q.sig()== G.sig());
1795  assert(F.size1()== size1);
1796  assert(F.size2()== size1);
1797  assert(Fcc.size1()== size1);
1798  assert(Fcc.size2()== size1);
1799  assert(Q.size1()== size1);
1800  assert(Q.size2()== size1);
1801  assert(n1 <= F.nt());
1802  assert(n1 <= Fcc.nt());
1803  assert(n1 <= G.nt());
1804  assert(n1 <= Q.nt());
1805 
1806  function<T> funFinv(n2,size1);
1807 
1808  cntr::herm_matrix<T> tmpF(nt,ntau,size1,sig),tmpFcc(nt,ntau,size1,sig),tmpF1(nt,ntau,size1,sig),tmpQ(nt,ntau,size1,sig);
1809  cdmatrix cdF,cdFinv,cdQ,cdG;
1810  //Check consistency
1811  assert(G.sig()==F.sig());
1812  assert(G.nt()==F.nt());
1813  assert(G.nt()==Q.nt());
1814 
1815  for(int n=-1;n<=n2;n++){
1816  Fsin.get_value(n,cdF);
1817  Qsin.get_value(n,cdQ);
1818  cdF=cdF+cdmatrix::Identity(size1,size1);
1819  cdFinv=cdF.inverse(); //Expected that these matrices are small
1820  cdG=cdFinv*cdQ;
1821  funFinv.set_value(n,cdFinv);
1822  Gsin.set_value(n,cdG);
1823  }
1824 
1825  tmpF=F;
1826  tmpFcc=Fcc;
1827  tmpQ=Q;
1828 
1829  //Set new F and Q
1830  for(int n=-1;n<=n2;n++){
1831  tmpF.left_multiply(n,funFinv,1.0);
1832  tmpFcc.right_multiply(n,funFinv,1.0);
1833  tmpQ.left_multiply(n,funFinv,1.0);
1834  }
1835  tmpF1=tmpF;
1836 
1837  for(int n=-1;n<=n2;n++){
1838  tmpF1.right_multiply(n,Gsin,1.0);
1839  }
1840  for(int i=-1;i<=n2;i++){
1841  tmpQ.incr_timestep(i,tmpF1,-1.0);
1842  }
1843  vie2_timestep(n,G,tmpF,tmpFcc,tmpQ,integration::I<T>(SolveOrder),beta,h);
1844  }
1845 
1846 
1848 // function calls with alpha and possibly function
1849 #if CNTR_USE_OMP == 1
1850 template <typename T, class GG, int SIZE1>
1851 void vie2_timestep_omp_dispatch(int omp_num_threads, int tstp, GG &B, CPLX alpha, GG &A,
1852  GG &Acc, CPLX *f0, CPLX *ft, GG &Q,
1853  integration::Integrator<T> &I, T beta, T h) {
1854  int kt = I.get_k();
1855  int ntau = A.ntau();
1856  int size1 = A.size1();
1857  int sc = size1 * size1;
1858  bool func = (ft == NULL ? false : true);
1859  int n1 = (tstp == -1 || tstp > kt ? tstp : kt);
1860  CPLX *ftcc, *f0cc;
1861  herm_matrix_timestep<T> Qtstp(tstp, ntau, size1);
1862  // save Q
1863  Q.get_timestep(tstp, Qtstp);
1864  B.set_timestep_zero(tstp);
1865  if (func) {
1866  ftcc = new CPLX[sc * n1];
1867  f0cc = new CPLX[sc];
1868  element_conj<T, SIZE1>(size1, f0cc, f0);
1869  for (int n = 0; n <= n1; n++)
1870  element_conj<T, SIZE1>(size1, ftcc + n * sc, ft + n * sc);
1871  } else {
1872  ftcc = NULL;
1873  f0cc = NULL;
1874  }
1875  {
1876  CPLX *mtemp = new CPLX[sc];
1877  CPLX *one = new CPLX[sc];
1878  // mtemp= A.retptr(tstp,tstp)*ft(tstp)*alpha*h [NB1 x NB1]
1879  if (func) {
1880  element_mult<T, SIZE1>(size1, mtemp, A.retptr(tstp, tstp), ft + sc * tstp);
1881  element_smul<T, SIZE1>(size1, mtemp, h * alpha);
1882  } else {
1883  element_set<T, SIZE1>(size1, mtemp, A.retptr(tstp, tstp));
1884  element_smul<T, SIZE1>(size1, mtemp, h * alpha);
1885  }
1886  // one = diag(1) [NB1 x NB1]
1887  element_set<T, SIZE1>(size1, one, CPLX(1, 0));
1888 #pragma omp parallel num_threads(omp_num_threads)
1889  {
1890  int nomp = omp_get_num_threads();
1891  int tid = omp_get_thread_num();
1892  int i, n, m;
1893  CPLX *mm = new CPLX[sc];
1894  T wt;
1895  std::vector<bool> mask_tv, mask_les, mask_ret;
1896  mask_tv = std::vector<bool>(ntau + 1, false);
1897  mask_ret = std::vector<bool>(tstp + 1, false);
1898  mask_les = std::vector<bool>(tstp + 1, false);
1899  for (i = 0; i <= ntau; i++)
1900  if (i % nomp == tid)
1901  mask_tv[i] = true;
1902  for (i = 0; i <= tstp; i++)
1903  if (i % nomp == tid)
1904  mask_ret[i] = true;
1905  for (int i = 0; i <= tstp; i++)
1906  if (mask_ret[i])
1907  mask_les[tstp - i] = true;
1908  mask_les[tstp] = false; // computed separately at the end
1909  // retarded part:
1910  // Q1 = Q - alpha*A*ft*B, with Bret(tstp,n)=0
1911  incr_convolution_ret<T, GG, SIZE1>(tstp, mask_ret, -alpha, Q, A, Acc, ft, B, B,
1912  I, h);
1913  for (n = 0; n <= tstp; n++) {
1914  if (mask_ret[n]) {
1915  wt = (tstp < kt ? I.poly_integration(n, tstp, tstp)
1916  : I.gregory_weights(tstp - n, 0));
1917  element_set<T, SIZE1>(size1, mm, one);
1918  element_incr<T, SIZE1>(size1, mm, wt, mtemp);
1919  element_linsolve_right<T, SIZE1>(size1, 1, B.retptr(tstp, n), mm,
1920  Q.retptr(tstp, n));
1921  }
1922  }
1923  // tv part:
1924  // Q -> Q - alpha*A*ft*B, with Btv(tstp,n)=0
1925  incr_convolution_tv<T, GG, SIZE1>(tstp, mask_tv, -alpha, Q, A, Acc, f0, ft, B, B,
1926  I, beta, h);
1927  for (m = 0; m <= ntau; m++) {
1928  if (mask_tv[m]) {
1929  wt = I.gregory_weights(tstp, tstp);
1930  element_set<T, SIZE1>(size1, mm, one);
1931  element_incr<T, SIZE1>(size1, mm, wt, mtemp);
1932  element_linsolve_right<T, SIZE1>(size1, 1, B.tvptr(tstp, m), mm,
1933  Q.tvptr(tstp, m));
1934  }
1935  }
1936  // les part:
1937  // Q -> Q - conj(alpha)*B*ftcc*Acc, with Bles(n,tstp)=0, n=0...tstp-1 (Bret
1938  // already known!)
1939  mask_les[tstp] = false;
1940  incr_convolution_les<T, GG, SIZE1>(tstp, mask_les, -conj(alpha), Q, B, B, f0cc,
1941  ftcc, Acc, A, I, beta, h);
1942  for (n = 0; n < tstp; n++) {
1943  if (mask_les[n]) {
1944  wt = I.gregory_weights(tstp, tstp);
1945  element_set<T, SIZE1>(size1, mm, one);
1946  element_incr<T, SIZE1>(size1, mm, wt, mtemp);
1947  element_conj<T, SIZE1>(size1, mm);
1948  element_linsolve_left<T, SIZE1>(size1, 1, B.lesptr(n, tstp), mm,
1949  Q.lesptr(n, tstp));
1950  }
1951  }
1952  }
1953  // Bles(tstp,tstp) (Bles(n<tstp,tstp) etc. enters the convolution!)
1954  {
1955  int n;
1956  CPLX *mm = new CPLX[sc];
1957  std::vector<bool> mask_les;
1958  T wt;
1959  n = tstp;
1960  mask_les = std::vector<bool>(tstp + 1, false);
1961  mask_les[n] = true;
1962  incr_convolution_les<T, GG, SIZE1>(tstp, mask_les, -conj(alpha), Q, B, B, f0cc,
1963  ftcc, Acc, A, I, beta, h);
1964  wt = I.gregory_weights(tstp, tstp);
1965  element_set<T, SIZE1>(size1, mm, one);
1966  element_incr<T, SIZE1>(size1, mm, wt, mtemp);
1967  element_conj<T, SIZE1>(size1, mm);
1968  element_linsolve_left<T, SIZE1>(size1, 1, B.lesptr(n, tstp), mm,
1969  Q.lesptr(n, tstp));
1970  delete[] mm;
1971  }
1972  delete[] one;
1973  delete[] mtemp;
1974  if (func) {
1975  delete[] f0cc;
1976  delete[] ftcc;
1977  }
1978  // restore Q
1979  Q.set_timestep(tstp, Qtstp);
1980  }
1981 }
1982 // same function call as above: anyway, only for timestep
1984 
2021 template <typename T>
2022 void vie2_timestep_omp(int omp_num_threads, int tstp, herm_matrix<T> &G, herm_matrix<T> &F,
2023  herm_matrix<T> &Fcc, herm_matrix<T> &Q, integration::Integrator<T> &I,
2024  T beta, T h, const int matsubara_method) {
2025  int ntau = G.ntau();
2026  int size1 = G.size1(), kt = I.k();
2027  int n1 = (tstp >= kt ? tstp : kt);
2028 
2029  assert(tstp >= 0);
2030  assert(ntau > 0);
2031  assert(kt > 0 && kt <= 5);
2032  assert(kt <= 2 * ntau + 2);
2033  assert(ntau == Fcc.ntau());
2034  assert(ntau == F.ntau());
2035  assert(ntau == Q.ntau());
2036  assert(F.sig()== G.sig());
2037  assert(Fcc.sig()== G.sig());
2038  assert(Q.sig()== G.sig());
2039  assert(F.size1()== size1);
2040  assert(F.size2()== size1);
2041  assert(Fcc.size1()== size1);
2042  assert(Fcc.size2()== size1);
2043  assert(Q.size1()== size1);
2044  assert(Q.size2()== size1);
2045  assert(n1 <= F.nt());
2046  assert(n1 <= Fcc.nt());
2047  assert(n1 <= G.nt());
2048  assert(n1 <= Q.nt());
2049 
2050  if (tstp==-1){
2051  vie2_mat(G, F, Fcc, Q, beta, I, matsubara_method);
2052  }else if(tstp<=kt){
2053  cntr::vie2_start(G,F,Fcc,Q,integration::I<double>(tstp),beta,h);
2054  }else{
2055  switch (size1){
2056  case 1:
2057  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 1>(
2058  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, I, beta, h);
2059  break;
2060  case 2:
2061  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 2>(
2062  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, I, beta, h);
2063  break;
2064  case 3:
2065  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 3>(
2066  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, I, beta, h);
2067  break;
2068  case 4:
2069  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 4>(
2070  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, I, beta, h);
2071  break;
2072  case 5:
2073  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 5>(
2074  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, I, beta, h);
2075  break;
2076  case 8:
2077  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 8>(
2078  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, I, beta, h);
2079  break;
2080  }
2081 
2082  }
2083 }
2084 
2085 
2086 
2087 
2125 template <typename T>
2126 void vie2_timestep_omp(int omp_num_threads, int tstp, herm_matrix<T> &G, herm_matrix<T> &F,
2127  herm_matrix<T> &Fcc, herm_matrix<T> &Q,
2128  T beta, T h, const int SolveOrder, const int matsubara_method) {
2129  int ntau = G.ntau();
2130  int size1 = G.size1();
2131  int n1 = (tstp >= SolveOrder ? tstp : SolveOrder);
2132 
2133  assert(tstp >= 0);
2134  assert(ntau > 0);
2135  assert(SolveOrder > 0 && SolveOrder <= 5);
2136  assert(SolveOrder <= 2 * ntau + 2);
2137  assert(ntau == Fcc.ntau());
2138  assert(ntau == F.ntau());
2139  assert(ntau == Q.ntau());
2140  assert(F.sig()== G.sig());
2141  assert(Fcc.sig()== G.sig());
2142  assert(Q.sig()== G.sig());
2143  assert(F.size1()== size1);
2144  assert(F.size2()== size1);
2145  assert(Fcc.size1()== size1);
2146  assert(Fcc.size2()== size1);
2147  assert(Q.size1()== size1);
2148  assert(Q.size2()== size1);
2149  assert(n1 <= F.nt());
2150  assert(n1 <= Fcc.nt());
2151  assert(n1 <= G.nt());
2152  assert(n1 <= Q.nt());
2153 
2154  if (tstp==-1){
2155  vie2_mat(G, F, Fcc, Q, beta, integration::I<T>(SolveOrder), matsubara_method);
2156  }else if(tstp<=SolveOrder){
2157  cntr::vie2_start(G,F,Fcc,Q,integration::I<T>(tstp),beta,h);
2158  }else{
2159  switch (size1){
2160  case 1:
2161  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 1>(
2162  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, integration::I<T>(SolveOrder), beta, h);
2163  break;
2164  case 2:
2165  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 2>(
2166  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, integration::I<T>(SolveOrder), beta, h);
2167  break;
2168  case 3:
2169  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 3>(
2170  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, integration::I<T>(SolveOrder), beta, h);
2171  break;
2172  case 4:
2173  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 4>(
2174  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, integration::I<T>(SolveOrder), beta, h);
2175  break;
2176  case 5:
2177  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 5>(
2178  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, integration::I<T>(SolveOrder), beta, h);
2179  break;
2180  case 8:
2181  vie2_timestep_omp_dispatch<T, herm_matrix<T>, 8>(
2182  omp_num_threads, tstp, G, CPLX(1, 0), F, Fcc, NULL, NULL, Q, integration::I<T>(SolveOrder), beta, h);
2183  break;
2184  }
2185 
2186  }
2187 }
2188 
2189 
2191 
2234 template < typename T>
2235 void vie2_timestep_sin_omp(int omp_num_threads, int tstp,herm_matrix<T> &G,function<T> &Gsin,
2236  herm_matrix<T> &F,herm_matrix<T> &Fcc, function<T> &Fsin ,
2237  herm_matrix<T> &Q,function<T> &Qsin,integration::Integrator<T> &I,T beta,T h){
2238 
2239  int kt=I.k();
2240  int n1=(tstp<=kt && tstp>=0 ? 0 : tstp);
2241  int n2=(tstp<=kt && tstp>=0 ? kt : tstp);
2242  int nt=F.nt(),ntau=F.ntau(),size1=F.size1(),sig=F.sig();
2243 
2244  assert(tstp >= 0);
2245  assert(ntau > 0);
2246  assert(kt > 0 && kt <= 5);
2247  assert(kt <= 2 * ntau + 2);
2248  assert(ntau == Fcc.ntau());
2249  assert(ntau == F.ntau());
2250  assert(ntau == Q.ntau());
2251  assert(F.sig()== G.sig());
2252  assert(Fcc.sig()== G.sig());
2253  assert(Q.sig()== G.sig());
2254  assert(F.size1()== size1);
2255  assert(F.size2()== size1);
2256  assert(Fcc.size1()== size1);
2257  assert(Fcc.size2()== size1);
2258  assert(Q.size1()== size1);
2259  assert(Q.size2()== size1);
2260  assert(n1 <= F.nt());
2261  assert(n1 <= Fcc.nt());
2262  assert(n1 <= G.nt());
2263  assert(n1 <= Q.nt());
2264 
2265  function<T> funFinv(n2,size1);
2266 
2267  cntr::herm_matrix<T> tmpF(nt,ntau,size1,sig),tmpFcc(nt,ntau,size1,sig),tmpF1(nt,ntau,size1,sig),tmpQ(nt,ntau,size1,sig);
2268  cdmatrix cdF,cdFinv,cdQ,cdG;
2269  //Check consistency
2270  assert(G.sig()==F.sig());
2271  assert(G.nt()==F.nt());
2272  assert(G.nt()==Q.nt());
2273 
2274  for(int n=-1;n<=n2;n++){
2275  Fsin.get_value(n,cdF);
2276  Qsin.get_value(n,cdQ);
2277  cdF=cdF+cdmatrix::Identity(size1,size1);
2278  cdFinv=cdF.inverse(); //Expected that these matrices are small
2279  cdG=cdFinv*cdQ;
2280  funFinv.set_value(n,cdFinv);
2281  Gsin.set_value(n,cdG);
2282  }
2283 
2284  tmpF=F;
2285  tmpFcc=Fcc;
2286  tmpQ=Q;
2287 
2288  //Set new F and Q
2289  for(int n=-1;n<=n2;n++){
2290  tmpF.left_multiply(n,funFinv,1.0);
2291  tmpFcc.right_multiply(n,funFinv,1.0);
2292  tmpQ.left_multiply(n,funFinv,1.0);
2293  }
2294  tmpF1=tmpF;
2295 
2296  for(int n=-1;n<=n2;n++){
2297  tmpF1.right_multiply(n,Gsin,1.0);
2298  }
2299  for(int i=-1;i<=n2;i++){
2300  tmpQ.incr_timestep(i,tmpF1,-1.0);
2301  }
2302 
2303  vie2_timestep_omp(omp_num_threads,tstp,G,tmpF,tmpFcc,tmpQ,I,beta,h);
2304  }
2305 
2306 
2307 
2351 template < typename T>
2352 void vie2_timestep_sin_omp(int omp_num_threads, int tstp,herm_matrix<T> &G,function<T> &Gsin,
2353  herm_matrix<T> &F,herm_matrix<T> &Fcc, function<T> &Fsin ,
2354  herm_matrix<T> &Q,function<T> &Qsin,T beta,T h, const int SolveOrder){
2355 
2356  int n1=(tstp<=SolveOrder && tstp>=0 ? 0 : tstp);
2357  int n2=(tstp<=SolveOrder && tstp>=0 ? SolveOrder : tstp);
2358  int nt=F.nt(),ntau=F.ntau(),size1=F.size1(),sig=F.sig();
2359 
2360  assert(tstp >= 0);
2361  assert(ntau > 0);
2362  assert(SolveOrder > 0 && SolveOrder <= 5);
2363  assert(SolveOrder <= 2 * ntau + 2);
2364  assert(ntau == Fcc.ntau());
2365  assert(ntau == F.ntau());
2366  assert(ntau == Q.ntau());
2367  assert(F.sig()== G.sig());
2368  assert(Fcc.sig()== G.sig());
2369  assert(Q.sig()== G.sig());
2370  assert(F.size1()== size1);
2371  assert(F.size2()== size1);
2372  assert(Fcc.size1()== size1);
2373  assert(Fcc.size2()== size1);
2374  assert(Q.size1()== size1);
2375  assert(Q.size2()== size1);
2376  assert(n1 <= F.nt());
2377  assert(n1 <= Fcc.nt());
2378  assert(n1 <= G.nt());
2379  assert(n1 <= Q.nt());
2380 
2381  function<T> funFinv(n2,size1);
2382 
2383  cntr::herm_matrix<T> tmpF(nt,ntau,size1,sig),tmpFcc(nt,ntau,size1,sig),tmpF1(nt,ntau,size1,sig),tmpQ(nt,ntau,size1,sig);
2384  cdmatrix cdF,cdFinv,cdQ,cdG;
2385  //Check consistency
2386  assert(G.sig()==F.sig());
2387  assert(G.nt()==F.nt());
2388  assert(G.nt()==Q.nt());
2389 
2390  for(int n=-1;n<=n2;n++){
2391  Fsin.get_value(n,cdF);
2392  Qsin.get_value(n,cdQ);
2393  cdF=cdF+cdmatrix::Identity(size1,size1);
2394  cdFinv=cdF.inverse(); //Expected that these matrices are small
2395  cdG=cdFinv*cdQ;
2396  funFinv.set_value(n,cdFinv);
2397  Gsin.set_value(n,cdG);
2398  }
2399 
2400  tmpF=F;
2401  tmpFcc=Fcc;
2402  tmpQ=Q;
2403 
2404  //Set new F and Q
2405  for(int n=-1;n<=n2;n++){
2406  tmpF.left_multiply(n,funFinv,1.0);
2407  tmpFcc.right_multiply(n,funFinv,1.0);
2408  tmpQ.left_multiply(n,funFinv,1.0);
2409  }
2410  tmpF1=tmpF;
2411 
2412  for(int n=-1;n<=n2;n++){
2413  tmpF1.right_multiply(n,Gsin,1.0);
2414  }
2415  for(int i=-1;i<=n2;i++){
2416  tmpQ.incr_timestep(i,tmpF1,-1.0);
2417  }
2418 
2419  vie2_timestep_omp(omp_num_threads,tstp,G,tmpF,tmpFcc,tmpQ,beta,h,SolveOrder);
2420  }
2421 
2422 
2424 
2458 template <typename T>
2459 void vie2_omp(int omp_num_threads, herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc,
2460  herm_matrix<T> &Q, integration::Integrator<T> &I, T beta, T h, const int matsubara_method) {
2461  int tstp, k = I.k();
2462  vie2_mat(G, F, Fcc, Q, beta, I, matsubara_method);
2463  // vie2_mat(G,F,Fcc,Q,beta,3);
2464  if (G.nt() >= 0)
2465  vie2_start(G, F, Fcc, Q, I, beta, h);
2466  for (tstp = k + 1; tstp <= G.nt(); tstp++)
2467  vie2_timestep_omp(omp_num_threads, tstp, G, F, Fcc, Q, I, beta, h);
2468 }
2469 
2470 
2471 
2472 
2507 template <typename T>
2508 void vie2_omp(int omp_num_threads, herm_matrix<T> &G, herm_matrix<T> &F, herm_matrix<T> &Fcc,
2509  herm_matrix<T> &Q, T beta, T h, const int SolveOrder, const int matsubara_method) {
2510  int tstp;
2511  vie2_mat(G, F, Fcc, Q, beta, SolveOrder, matsubara_method);
2512  if (G.nt() >= 0)
2513  vie2_start(G, F, Fcc, Q, beta, h, SolveOrder);
2514  for (tstp = SolveOrder + 1; tstp <= G.nt(); tstp++)
2515  vie2_timestep_omp(omp_num_threads, tstp, G, F, Fcc, Q, beta, h, SolveOrder);
2516 }
2517 
2518 #endif // CNTR_USE_OMP
2519 
2520 #undef CPLX
2521 
2522 } // namespace cntr
2523 
2524 #endif // CNTR_VIE2_IMPL_H
integration::Integrator
Class Integrator contains all kinds of weights for integration and differentiation of a function at ...
Definition: integration.hpp:128
cntr_herm_matrix_timestep_decl.hpp
cntr_elements.hpp
cntr::vie2_timestep_sin
void vie2_timestep_sin(int n, herm_matrix< T > &G, function< T > &Gsin, herm_matrix< T > &F, herm_matrix< T > &Fcc, function< T > &Fsin, herm_matrix< T > &Q, function< T > &Qsin, T beta, T h, int SolveOrder=MAX_SOLVE_ORDER)
One step VIE solver for Green&#39;s function with instantaneous contributions for given integration ord...
Definition: cntr_vie2_impl.hpp:1779
cntr_matsubara_impl.hpp
fourier::cplx
std::complex< double > cplx
Definition: fourier.cpp:11
cntr::herm_matrix::nt
int nt(void) const
Definition: cntr_herm_matrix_decl.hpp:71
cntr_convolution_impl.hpp
cntr::herm_matrix::incr_timestep
void incr_timestep(int tstp, herm_matrix_timestep< T > &timestep, cplx alpha)
Adds a herm_matrix_timestep with given weight to the herm_matrix.
Definition: cntr_herm_matrix_impl.hpp:2567
integration::I
Integrator< T > & I(int k)
Definition: integration.hpp:506
cntr::herm_matrix::size2
int size2(void) const
Definition: cntr_herm_matrix_decl.hpp:69
cntr::herm_matrix::sig
int sig(void) const
Definition: cntr_herm_matrix_decl.hpp:72
cntr::function
Class function for objects with time on real axis.
Definition: cntr_convolution_decl.hpp:9
cntr::function::get_value
void get_value(int tstp, EigenMatrix &M) const
Get matrix value of this function object at a specific time point
Definition: cntr_function_impl.hpp:414
cntr_herm_matrix_decl.hpp
cntr::herm_matrix
Class herm_matrix for two-time contour objects with hermitian symmetry.
Definition: cntr_convolution_decl.hpp:10
cntr::vie2_timestep_sin_omp
void vie2_timestep_sin_omp(int omp_num_threads, int tstp, herm_matrix< T > &G, function< T > &Gsin, herm_matrix< T > &F, herm_matrix< T > &Fcc, function< T > &Fsin, herm_matrix< T > &Q, function< T > &Qsin, T beta, T h, const int SolveOrder=MAX_SOLVE_ORDER)
One step VIE solver for Green&#39;s function with instantaneous contributions for given integration ord...
Definition: cntr_vie2_impl.hpp:2352
cntr::herm_matrix::ntau
int ntau(void) const
Definition: cntr_herm_matrix_decl.hpp:70
cntr_vie2_decl.hpp
cntr::vie2_start
void vie2_start(herm_matrix< T > &G, herm_matrix< T > &F, herm_matrix< T > &Fcc, herm_matrix< T > &Q, T beta, T h, const int SolveOrder=MAX_SOLVE_ORDER)
VIE solver for a Green&#39;s function for the first k timesteps
Definition: cntr_vie2_impl.hpp:1420
cntr::function::set_value
void set_value(int tstp, EigenMatrix &M)
Set matrix value at a specific time point
Definition: cntr_function_impl.hpp:383
cntr::herm_matrix::size1
int size1(void) const
Definition: cntr_herm_matrix_decl.hpp:68
integration::I< double >
template Integrator< double > & I< double >(int k)