NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_dyson_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_DYSON_IMPL_H
2 #define CNTR_DYSON_IMPL_H
3 
4 #include "cntr_dyson_decl.hpp"
5 //#include "cntr_exception.hpp"
6 #include "cntr_elements.hpp"
7 #include "cntr_function_decl.hpp"
8 #include "cntr_herm_matrix_decl.hpp"
9 #include "cntr_matsubara_impl.hpp"
10 #include "cntr_convolution_impl.hpp"
11 #include "cntr_equilibrium_decl.hpp"
12 #include "cntr_vie2_decl.hpp"
13 #include "cntr_utilities_decl.hpp"
14 
15 namespace cntr {
16 
17 /* #######################################################################################
18 # DYSON EQUATION:
19 #
20 # [ id/dt + mu - H(t) ] G(t,t') - [Sigma*G](t,t') = 1(t,t')
21 # G(t,t')[-id/dt' + mu - H(t')] - [G*Sigma](t,t') = 1(t,t')
22 #
23 # (*) "_timestep" computes G(t,t') at timestep n, i.e., G^ret(nh,t'<=nh),
24 # G^les(t<=nh,nh), G^tv(nt,tau=0..beta). Timestep must be >k.
25 #
26 # (*) For the computation of timestep n, G(t,t') and Sigma(t,t') are adressed at
27 # times t,t' <= max(n,k), where k is the Integration order (see Integrator).
28 # I.e., the timesteps n=0..k must be computed seperately, using the routine
29 # "_start", which assumes that the Matsubara component of G and Sigma(t,t')
30 # for t,t'<=k are given.
31 #
32 # (*) The following types are allowed for G and Sigma:
33 # - G=scalar, Sigma=scalar
34 # - G=matrix, Sigma=scalar,matrix, (sparse matrix)
35 # Htype is some type that can be converted into G by G.element_set
36 # The Dyson equation assumes H,G, and Sigma to be hermitian.
37 #
38 ###########################################################################################*/
39 
40 /*###########################################################################################
41 #
42 # RETARDED FUNCTION:
43 #
44 ###########################################################################################*/
45 
47 
79 template <typename T, class GG, int SIZE1>
80 void dyson_timestep_ret(int n, GG &G, T mu, std::complex<T> *H, GG &Sigma,
81  integration::Integrator<T> &I, T h) {
82  typedef std::complex<T> cplx;
83  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1;
84  int ss, sg, n1, l, j, i, p, q, m, j1, j2, size1 = G.size1();
85  cplx *gret, *gtemp, *mm, *qq, *qqj, *stemp, *one, cplx_i, cweight, *diffw;
86  cplx *sret, w0, *hj;
87  T weight;
88  cplx_i = cplx(0, 1);
89 
90  ss = Sigma.element_size();
91  sg = G.element_size();
92  gtemp = new cplx[k * sg];
93  diffw = new cplx[k1 + 1];
94  qq = new cplx[(n + 1) * sg];
95  one = new cplx[sg];
96  mm = new cplx[k * k * sg];
97  hj = new cplx[sg];
98  stemp = new cplx[sg]; // sic
99  element_set<T, SIZE1>(size1, one, 1);
100  // check consistency:
101  assert(n > k);
102  assert(Sigma.nt() >= n);
103  assert(G.nt() >= n);
104  assert(G.sig() == Sigma.sig());
105  // SET ENTRIES IN TIMESTEP TO 0
106  gret = G.retptr(n, 0);
107  n1 = (n + 1) * sg;
108  for (i = 0; i < n1; i++)
109  gret[i] = 0;
110  // INITIAL VALUE t' = n
111  element_set<T, SIZE1>(size1, G.retptr(n, n), -cplx_i);
112  // START VALUES t' = n-j, j = 1...k: solve a kxk problem
113  for (i = 0; i < k * k * sg; i++)
114  mm[i] = 0;
115  for (i = 0; i < k * sg; i++)
116  qq[i] = 0;
117  for (j = 1; j <= k; j++) {
118  p = j - 1;
119  // derivatives:
120  for (l = 0; l <= k; l++) {
121  cweight = cplx_i / h * I.poly_differentiation(j, l);
122  if (l == 0) {
123  element_incr<T, SIZE1>(size1, qq + p * sg, -cweight, G.retptr(n, n));
124  } else {
125  q = l - 1;
126  element_incr<T, SIZE1>(size1, mm + sg * (p * k + q), cweight);
127  }
128  }
129  // H
130  element_set<T, SIZE1>(size1, gtemp, H + (n - j) * sg);
131  element_smul<T, SIZE1>(size1, gtemp, -1);
132  for (i = 0; i < sg; i++)
133  gtemp[i] += mu * one[i];
134  element_incr<T, SIZE1>(size1, mm + sg * (p + k * p), gtemp);
135  // integral
136  for (l = 0; l <= k; l++) {
137  weight = h * I.gregory_weights(j, l);
138  if (n - l >= n - j) {
139  element_set<T, SIZE1>(
140  size1, stemp,
141  Sigma.retptr(n - l, n - j)); // stemp is element of type G!!
142  } else {
143  element_set<T, SIZE1>(size1, stemp, Sigma.retptr(n - j, n - l));
144  element_conj<T, SIZE1>(size1, stemp);
145  weight *= -1;
146  }
147  if (l == 0) {
148  element_incr<T, SIZE1>(size1, qq + p * sg, weight, G.retptr(n, n), stemp);
149  } else {
150  q = l - 1;
151  element_incr<T, SIZE1>(size1, mm + sg * (p * k + q), -weight, stemp);
152  }
153  }
154  }
155  element_linsolve_left<T, SIZE1>(size1, k, gtemp, mm, qq); // gtemp * mm = qq
156  for (j = 1; j <= k; j++)
157  element_set<T, SIZE1>(size1, G.retptr(n, n - j), gtemp + (j - 1) * sg);
158  // Compute the contribution to the convolution int dm G(n,m)Sigma(m,j)
159  // from G(n,m=n-k..n) to m=0...n-k, store into qq(j), without factor h!!
160  for (i = 0; i < sg * (n + 1); i++)
161  qq[i] = 0;
162  for (m = n - k; m <= n; m++) {
163  weight = I.gregory_omega(n - m);
164  gret = G.retptr(n, m);
165  sret = Sigma.retptr(m, 0);
166  for (j = 0; j <= n - k2; j++) {
167  element_incr<T, SIZE1>(size1, qq + j * sg, I.gregory_weights(n - j, n - m), gret,
168  Sigma.retptr(m, j));
169  sret += ss;
170  }
171  j1 = (n - k2 + 1 < 0 ? 0 : n - k2 + 1);
172  for (j = j1; j < n - k; j++) { // start weights
173  weight = I.gregory_weights(n - j, n - m);
174  element_incr<T, SIZE1>(size1, qq + j * sg, I.gregory_weights(n - j, n - m), gret,
175  Sigma.retptr(m, j));
176  sret += ss;
177  }
178  }
179  // DER REST VOM SCHUETZENFEST: t' = n-l, l = k+1 ... n:
180  for (p = 0; p <= k1; p++)
181  diffw[p] = I.bd_weights(p) * cplx_i / h; // use BD(k+1!!)
182  w0 = h * I.gregory_omega(0);
183  for (l = k + 1; l <= n; l++) {
184  j = n - l;
185  element_set<T, SIZE1>(size1, hj, H + j * sg);
186  element_conj<T, SIZE1>(size1, hj);
187  // set up mm and qqj for 1x1 problem:
188  qqj = qq + j * sg;
189  for (i = 0; i < sg; i++) {
190  qqj[i] *= h;
191  for (p = 1; p <= k1; p++)
192  qqj[i] += -diffw[p] * G.retptr(n, n - l + p)[i];
193  }
194  element_set<T, SIZE1>(size1, stemp, Sigma.retptr(j, j));
195  for (i = 0; i < sg; i++)
196  mm[i] = diffw[0] * one[i] + mu * one[i] - hj[i] - w0 * stemp[i];
197  element_linsolve_left<T, SIZE1>(size1, G.retptr(n, j), mm, qqj);
198  // compute the contribution of Gret(n,j) to
199  // int dm G(n,j)Sigma(j,j2) for j2<j, store into qq(j2), without h
200  gret = G.retptr(n, j);
201  sret = Sigma.retptr(j, 0);
202  for (j2 = 0; j2 < j - k; j2++) {
203  element_incr<T, SIZE1>(size1, qq + j2 * sg, I.gregory_weights(n - j2, n - j),
204  gret, sret);
205  sret += ss;
206  }
207  j1 = (j - k < 0 ? 0 : j - k);
208  for (j2 = j1; j2 < j; j2++) { // start weights
209  weight = I.gregory_omega(j - j2);
210  element_incr<T, SIZE1>(size1, qq + j2 * sg, I.gregory_weights(n - j2, n - j),
211  gret, sret);
212  sret += ss;
213  }
214  }
215  delete[] diffw;
216  delete[] stemp;
217  delete[] qq;
218  delete[] mm;
219  delete[] gtemp;
220  delete[] one;
221  delete[] hj;
222  return;
223 }
224 
226 
257 template <typename T, class GG, int SIZE1>
258 void dyson_start_ret(GG &G, T mu, std::complex<T> *H, GG &Sigma,
259  integration::Integrator<T> &I, T h) {
260  typedef std::complex<T> cplx;
261  int k = I.get_k();
262  int sg, l, n, j, p, q, i, dim, size1 = G.size1();
263  cplx ih, *gtemp, *mm, *qq, *stemp, minusi, *one, *hj;
264  T weight;
265 
266  sg = G.element_size();
267  assert(Sigma.nt() >= k);
268  assert(G.nt() >= k);
269  assert(G.sig() == Sigma.sig());
270 
271  // temporary storage
272  qq = new cplx[k * sg];
273  mm = new cplx[k * k * sg];
274  one = new cplx[sg];
275  hj = new cplx[sg];
276  gtemp = new cplx[k * sg];
277  stemp = new cplx[sg]; // sic
278 
279  element_set<T, SIZE1>(size1, one, 1.0);
280  minusi = cplx(0, -1);
281 
282  // set initial values:
283  for (n = 0; n <= k; n++)
284  element_set<T, SIZE1>(size1, G.retptr(n, n), cplx(0, -1.0));
285 
286  for (j = 0; j < k; j++) { // determine G(n,j), n=j+1..k
287  dim = k - j;
288  for (i = 0; i < k * sg; i++)
289  qq[i] = 0;
290  for (i = 0; i < k * k * sg; i++)
291  mm[i] = 0;
292  // fill the matrix mm, indices p,q
293  for (n = j + 1; n <= k; n++) {
294  p = n - (j + 1);
295  for (l = 0; l <= k; l++) {
296  q = l - (j + 1);
297  if (l <= j) { // this involves G which is known and goes into qq(p)
298  element_conj<T, SIZE1>(size1, gtemp, G.retptr(j, l));
299  element_smul<T, SIZE1>(size1, gtemp, -1);
300  element_set_zero<T, SIZE1>(size1, stemp);
301  element_incr<T, SIZE1>(size1, stemp, Sigma.retptr(n, l), gtemp);
302  for (i = 0; i < sg; i++)
303  qq[p * sg + i] +=
304  minusi / h * I.poly_differentiation(n, l) * gtemp[i] +
305  h * I.poly_integration(j, n, l) * stemp[i];
306  } else { // this goes into mm(p,q)
307  for (i = 0; i < sg; i++)
308  mm[sg * (p * dim + q) + i] =
309  -minusi / h * I.poly_differentiation(n, l) * one[i];
310  if (l == n) {
311  element_set<T, SIZE1>(size1, hj, H + l * sg);
312  for (i = 0; i < sg; i++) {
313  mm[sg * (p * dim + q) + i] += mu * one[i] - hj[i];
314  }
315  }
316  weight = h * I.poly_integration(j, n, l);
317  if (n >= l) {
318  element_set<T, SIZE1>(size1, stemp, Sigma.retptr(n, l));
319  } else {
320  element_set<T, SIZE1>(size1, stemp, Sigma.retptr(l, n));
321  element_conj<T, SIZE1>(size1, stemp);
322  weight *= -1;
323  }
324  for (i = 0; i < sg; i++)
325  mm[sg * (p * dim + q) + i] -= weight * stemp[i];
326  }
327  }
328  }
329  // solve mm*gtemp=qq
330  element_linsolve_right<T, SIZE1>(size1, dim, gtemp, mm, qq);
331  for (n = j + 1; n <= k; n++) {
332  p = n - (j + 1);
333  // std::cerr << n << " " << j << " " << gtemp[p] << std:: endl;
334  element_set<T, SIZE1>(size1, G.retptr(n, j), gtemp + p * sg);
335  }
336  }
337  delete[] stemp;
338  delete[] qq;
339  delete[] mm;
340  delete[] gtemp;
341  delete[] one;
342  delete[] hj;
343  return;
344 }
345 /*####################################################################################
346 #
347 # MATSUBARA FUNCTION g(tau) = 1/beta sum_n e^{-iomn tau} (iomn-H-Sigma(iomn))^{-1}
348 #
349 ####################################################################################*/
351 // TODO:
352 // For Sigma=0 the boundary values seem to be a very badly approximated, although
353 // the jump=1 is accounted far correctly. Sigma!=0 somehow makes the performance
354 // better
355 template <typename T, class GG, int SIZE1>
356 void dyson_mat_fourier_dispatch(GG &G, GG &Sigma, T mu, std::complex<T> *H0, T beta, int order = 3) {
357  typedef std::complex<double> cplx;
358  cplx *sigmadft, *sigmaiomn, *z1, *z2, *one;
359  cplx *expfac, *gmat, *hj, iomn, *zinv;
360  int ntau, m, r, pcf, p, m2, sg, ss, l, sig, size1 = G.size1();
361  double dtau;
362 
363  assert(G.ntau() == Sigma.ntau());
364  sig = G.sig();
365  assert(G.sig() == Sigma.sig());
366  sg = G.element_size();
367  ss = Sigma.element_size();
368  ntau = G.ntau();
369  dtau = beta / ntau;
370  if (ntau % 2 == 1) {
371  std::cerr << "matsubara_inverse: ntau odd" << std::endl;
372  abort();
373  }
374  sigmadft = new cplx[(ntau + 1) * ss];
375  sigmaiomn = new cplx[ss];
376  expfac = new cplx[ntau + 1];
377  gmat = new cplx[(ntau + 1) * sg];
378  z1 = new cplx[sg];
379  z2 = new cplx[sg];
380  hj = new cplx[sg];
381  one = new cplx[sg];
382  zinv = new cplx[sg];
383  element_set<T, SIZE1>(size1, one, 1.0);
384  element_set<T, SIZE1>(size1, hj, H0);
385  for (l = 0; l < sg; l++)
386  hj[l] -= mu * one[l];
387  pcf = 10;
388  m2 = ntau / 2;
389  matsubara_dft<T, GG, SIZE1>(sigmadft, Sigma, sig);
390  set_first_order_tail<T, SIZE1>(gmat, one, beta, sg, ntau, sig, size1);
391 
392  for (m = -m2; m <= m2 - 1; m++) {
393  for (p = -pcf; p <= pcf; p++) {
394 
395  iomn = cplx(0, get_omega(m + p * ntau, beta, sig));
396  matsubara_ft<T, GG, SIZE1>(sigmaiomn, m + p * ntau, Sigma, sigmadft, sig, beta,
397  order);
398 
399  element_set<T, SIZE1>(size1, z1, sigmaiomn); // convert
400  element_incr<T, SIZE1>(size1, z1, hj); // z1=H+Sigma(iomn)
401 
402  if (sig == 1 && m + p * ntau == 0) {
403  // For Bosons we need to special treat the zeroth frequency
404  // as 1/iwn diverges.
405 
406  // z2 = 1./(-H - S)
407  element_inverse<T, SIZE1>(size1, zinv, z1);
408  element_set<T, SIZE1>(size1, z2, zinv);
409  element_smul<T, SIZE1>(size1, z2, -1.0);
410 
411  } else {
412 
413  // z1 = H + S
414  // zinv = 1/( iwn * ( iwn - H - S ))
415  // z2 = 1/(iwn - H - S) - 1/iwn = (H + S) / (iwn * (iwn - H - S))
416 
417  for (l = 0; l < sg; l++)
418  z2[l] = iomn * one[l] - z1[l];
419  element_inverse<T, SIZE1>(size1, zinv, z2);
420 
421  element_smul<T, SIZE1>(size1, zinv, 1. / iomn);
422  element_mult<T, SIZE1>(size1, z2, zinv, z1);
423  }
424 
425  element_smul<T, SIZE1>(size1, z2, 1 / beta);
426  for (r = 0; r <= ntau; r++) {
427  cplx expfac(std::exp(-get_tau(r, beta, ntau) * iomn));
428  for (l = 0; l < sg; l++)
429  gmat[r * sg + l] += z2[l] * expfac;
430  }
431  }
432  }
433  for (r = 0; r <= ntau; r++) {
434  element_set<T, SIZE1>(size1, G.matptr(r), gmat + r * sg);
435  }
436 
437  delete[] sigmadft;
438  delete[] sigmaiomn;
439  delete[] expfac;
440  delete[] gmat;
441  delete[] z1;
442  delete[] z2;
443  delete[] hj;
444  delete[] one;
445  delete[] zinv;
446  return;
447 }
448 
450 template <typename T, class GG, int SIZE1>
451 void dyson_mat_fixpoint_dispatch(GG &G, GG &Sigma, T mu, cdmatrix &H0,
452  integration::Integrator<T> &I, T beta, int fixpiter){
453 
454  int k = I.get_k();
455  int ntau = G.ntau(), size1=G.size1();
456  int iter;
457  T hdummy=1.0;
458  herm_matrix<T> G0(-1,ntau,size1,G.sig());
459  herm_matrix<T> G0xSGM(-1,ntau,size1,G.sig());
460 
461  green_from_H(G0, mu, H0, beta, hdummy);
462 
463  convolution_matsubara(G0xSGM, G0, Sigma, I, beta);
464  G0xSGM.smul(-1,-1);
465 
466  vie2_mat_fixpoint(G, G0xSGM, G0xSGM, G0, beta, I, fixpiter);
467 
468 }
469 
471 template <typename T, class GG, int SIZE1>
472 void dyson_mat_fixpoint_dispatch(GG &G, GG &Sigma, T mu, cdmatrix &H0, function<T> &SigmaMF,
473  integration::Integrator<T> &I, T beta, int fixpiter){
474 
475  int k = I.get_k();
476  int ntau = G.ntau(), size1=G.size1();
477  T hdummy=1.0;
478  herm_matrix<T> G0(-1,ntau,size1,G.sig());
479  herm_matrix<T> G0xSGM(-1,ntau,size1,G.sig());
480  herm_matrix<T> G0xMF(-1,ntau,size1,G.sig());
481 
482  green_from_H(G0, mu, H0, beta, hdummy);
483 
484  convolution_matsubara(G0xSGM, G0, Sigma, I, beta);
485  G0xMF.set_timestep(-1,G0);
486  G0xMF.right_multiply(-1,SigmaMF);
487  G0xSGM.incr_timestep(-1,G0xMF);
488  G0xSGM.smul(-1,-1);
489 
490  vie2_mat_fixpoint(G, G0xSGM, G0xSGM, G0, beta, I, fixpiter);
491 
492 }
493 
495 template <typename T, class GG, int SIZE1>
496 void dyson_mat_fixpoint_dispatch(GG &G, GG &Sigma, T mu, cdmatrix &H0, cdmatrix &SigmaMF,
497  integration::Integrator<T> &I, T beta, int fixpiter){
498 
499  int k = I.get_k();
500  int ntau = G.ntau(), size1=G.size1();
501  T hdummy=1.0;
502  herm_matrix<T> G0(-1,ntau,size1,G.sig());
503  herm_matrix<T> G0xSGM(-1,ntau,size1,G.sig());
504  herm_matrix<T> G0xMF(-1,ntau,size1,G.sig());
505  function<T> sgmf(-1,size1);
506 
507  green_from_H(G0, mu, H0, beta, hdummy);
508 
509  convolution_matsubara(G0xSGM, G0, Sigma, I, beta);
510  sgmf.set_value(-1,SigmaMF);
511  G0xMF.set_timestep(-1,G0);
512  G0xMF.right_multiply(-1,sgmf);
513  G0xSGM.incr_timestep(-1,G0xMF);
514  G0xSGM.smul(-1,-1);
515 
516  vie2_mat_fixpoint(G, G0xSGM, G0xSGM, G0, beta, I, fixpiter);
517 
518 }
519 
521 template <typename T, class GG, int SIZE1>
522 void dyson_mat_steep_dispatch(GG &G, GG &Sigma, T mu, cdmatrix &H0,
523  integration::Integrator<T> &I, T beta, int maxiter, T tol){
524 
525  int k = I.get_k();
526  int ntau = G.ntau(), size1=G.size1();
527  T hdummy=1.0;
528  herm_matrix<T> G0(-1,ntau,size1,G.sig());
529  herm_matrix<T> G0xSGM(-1,ntau,size1,G.sig());
530  herm_matrix<T> SGMxG0(-1,ntau,size1,G.sig());
531 
532  green_from_H(G0, mu, H0, beta, hdummy);
533 
534  convolution_matsubara(G0xSGM, G0, Sigma, I, beta);
535  convolution_matsubara(SGMxG0, Sigma, G0, I, beta);
536  G0xSGM.smul(-1,-1);
537  SGMxG0.smul(-1,-1);
538 
539  vie2_mat_steep(G, G0xSGM, SGMxG0, G0, beta, I, maxiter, tol);
540 
541 }
542 
544 template <typename T, class GG, int SIZE1>
545 void dyson_mat_steep_dispatch(GG &G, GG &Sigma, T mu, cdmatrix &H0, function<T> &SigmaMF,
546  integration::Integrator<T> &I, T beta, int maxiter, T tol){
547 
548  int k = I.get_k();
549  int ntau = G.ntau(), size1=G.size1();
550  T hdummy=1.0;
551  herm_matrix<T> G0(-1,ntau,size1,G.sig());
552  herm_matrix<T> G0xSGM(-1,ntau,size1,G.sig());
553  herm_matrix<T> SGMxG0(-1,ntau,size1,G.sig());
554  herm_matrix<T> G0xMF(-1,ntau,size1,G.sig());
555  herm_matrix<T> MFxG0(-1,ntau,size1,G.sig());
556 
557  green_from_H(G0, mu, H0, beta, hdummy);
558 
559  G0xMF.set_timestep(-1,G0);
560  G0xMF.right_multiply(-1,SigmaMF);
561  MFxG0.set_timestep(-1,G0);
562  MFxG0.left_multiply(-1,SigmaMF);
563 
564  convolution_matsubara(G0xSGM, G0, Sigma, I, beta);
565  convolution_matsubara(SGMxG0, Sigma, G0, I, beta);
566  G0xSGM.incr_timestep(-1,G0xMF);
567  G0xSGM.smul(-1,-1);
568  SGMxG0.incr_timestep(-1,MFxG0);
569  SGMxG0.smul(-1,-1);
570 
571  vie2_mat_steep(G, G0xSGM, SGMxG0, G0, beta, I, maxiter, tol);
572 
573 }
574 
576 template <typename T, class GG, int SIZE1>
577 void dyson_mat_steep_dispatch(GG &G, GG &Sigma, T mu, cdmatrix &H0, cdmatrix &SigmaMF,
578  integration::Integrator<T> &I, T beta, int maxiter, T tol){
579 
580  int k = I.get_k();
581  int ntau = G.ntau(), size1=G.size1();
582  T hdummy=1.0;
583  herm_matrix<T> G0(-1,ntau,size1,G.sig());
584  herm_matrix<T> G0xSGM(-1,ntau,size1,G.sig());
585  herm_matrix<T> SGMxG0(-1,ntau,size1,G.sig());
586  herm_matrix<T> G0xMF(-1,ntau,size1,G.sig());
587  herm_matrix<T> MFxG0(-1,ntau,size1,G.sig());
588  function<T> sgmf(-1,size1);
589 
590  green_from_H(G0, mu, H0, beta, hdummy);
591 
592  sgmf.set_value(-1,SigmaMF);
593  G0xMF.set_timestep(-1,G0);
594  G0xMF.right_multiply(-1,sgmf);
595  MFxG0.set_timestep(-1,G0);
596  MFxG0.left_multiply(-1,sgmf);
597 
598  convolution_matsubara(G0xSGM, G0, Sigma, I, beta);
599  convolution_matsubara(SGMxG0, G0, Sigma, I, beta);
600  G0xSGM.incr_timestep(-1,G0xMF);
601  G0xSGM.smul(-1,-1);
602  SGMxG0.incr_timestep(-1,MFxG0);
603  SGMxG0.smul(-1,-1);
604 
605  vie2_mat_steep(G, G0xSGM, SGMxG0, G0, beta, I, maxiter, tol);
606 
607 }
608 
609 
610 /*###########################################################################################
611 #
612 # tv-FUNCTION:
613 #
614 ###########################################################################################*/
616 
650 template <typename T, class GG, int SIZE1>
651 void dyson_timestep_tv(int n, GG &G, T mu, std::complex<T> *Hn, GG &Sigma,
652  integration::Integrator<T> &I, T beta, T h) {
653  typedef std::complex<T> cplx;
654  int k = I.get_k();
655  int sg, n1, l, m, j, ntau, size1 = G.size1();
656  cplx *gtv, *gtv1, cweight, ih, *mm, *qq, *stemp, minusi, *one, *htemp;
657  T weight;
658 
659  sg = G.element_size();
660  ntau = G.ntau();
661  // check consistency: (more assertations follow in convolution)
662  assert(n > k);
663  assert(Sigma.nt() >= n);
664  assert(G.nt() >= n);
665  assert(G.sig() == Sigma.sig());
666 
667  one = new cplx[sg];
668  qq = new cplx[sg];
669  mm = new cplx[sg];
670  stemp = new cplx[sg]; // sic
671  htemp = new cplx[sg];
672  element_set<T, SIZE1>(size1, one, 1.0);
673  // SET ENTRIES IN TIMESTEP(TV) TO 0
674  gtv = G.tvptr(n, 0);
675  n1 = (ntau + 1) * sg;
676  for (l = 0; l < n1; l++)
677  gtv[l] = 0;
678  // CONVOLUTION SIGMA*G: ---> Gtv(n,m)
679  convolution_timestep_tv<T, herm_matrix<T>, SIZE1>(n, G, Sigma, Sigma, G, G, I, beta,
680  h); // note: this sets only tv
681  // ACCUMULATE CONTRIBUTION TO id/dt G(t,t') FROM t=mh, m=n-k..n-1
682  ih = cplx(0, 1 / h);
683  for (m = n - k - 1; m < n; m++) {
684  cweight = ih * I.bd_weights(n - m); // use BD(k+1!!)
685  // G(n,j) -= cweight*G(m,j), for j=0...m
686  n1 = (ntau + 1) * sg;
687  gtv = G.tvptr(m, 0);
688  gtv1 = G.tvptr(n, 0);
689  for (l = 0; l < n1; l++)
690  gtv1[l] -= cweight * gtv[l];
691  }
692  // Now solve
693  // [ i/h bd(0) - H - h w(n,0) Sigma(n,n) ] G(n,m) = Q(m),
694  // where Q is initially stored in G(n,m)
695  element_set<T, SIZE1>(size1, stemp, Sigma.retptr(n, n));
696  weight = -h * I.gregory_weights(n, 0);
697  element_set<T, SIZE1>(size1, htemp, Hn);
698  for (l = 0; l < sg; l++)
699  mm[l] = ih * I.bd_weights(0) * one[l] + weight * stemp[l] + mu * one[l] - htemp[l];
700  for (j = 0; j <= ntau; j++) {
701  for (l = 0; l < sg; l++)
702  qq[l] = G.tvptr(n, j)[l];
703  element_linsolve_right<T, SIZE1>(size1, G.tvptr(n, j), mm, qq);
704  }
705  delete[] stemp;
706  delete[] htemp;
707  delete[] qq;
708  delete[] mm;
709  delete[] one;
710  return;
711 }
713 
746 template <typename T, class GG, int SIZE1>
747 void dyson_start_tv(GG &G, T mu, std::complex<T> *H, GG &Sigma,
748  integration::Integrator<T> &I, T beta, T h) {
749  typedef std::complex<T> cplx;
750  int k = I.get_k();
751  int sg, l, m, j, ntau, p, q, n, sig, size1 = G.size1();
752  cplx cweight, *mm, *qq, *stemp, *one, *gtemp;
753  cplx cplx_i = cplx(0.0, 1.0);
754  T dtau;
755  sg = G.element_size();
756  ntau = G.ntau();
757  dtau = beta / ntau;
758  sig = G.sig();
759  // check consistency: (more assertations follow in convolution)
760  assert(Sigma.nt() >= k);
761  assert(G.nt() >= k);
762  assert(Sigma.ntau() == ntau);
763  assert(G.sig() == Sigma.sig());
764 
765  one = new cplx[sg];
766  qq = new cplx[k * sg];
767  mm = new cplx[k * k * sg];
768  stemp = new cplx[sg]; // sic
769  gtemp = new cplx[k * sg];
770  element_set<T, SIZE1>(size1, one, 1.0);
771  // INITIAL VALUE: G^tv(0,tau) = i sgn G^mat(beta-tau) (sgn=Bose/Fermi)
772  for (m = 0; m <= ntau; m++)
773  for (l = 0; l < sg; l++)
774  G.tvptr(0, m)[l] = ((T)sig) * cplx_i * G.matptr(ntau - m)[l];
775  // CONVOLUTION -i int dtau Sigma^tv(n,tau)G^mat(tau,m) ---> Gtv(n,m)
776  for (n = 1; n <= k; n++) {
777  for (m = 0; m <= ntau; m++) {
778  matsubara_integral_2<T, SIZE1>(size1, m, ntau, gtemp, Sigma.tvptr(n, 0),
779  G.matptr(0), I, G.sig());
780  for (l = 0; l < sg; l++)
781  G.tvptr(n, m)[l] = dtau * gtemp[l];
782  }
783  }
784  // loop over m --> determine G^tv(n,m) for n=1...k
785  for (m = 0; m <= ntau; m++) {
786  for (l = 0; l < k * k * sg; l++)
787  mm[l] = 0;
788  for (l = 0; l < k * sg; l++)
789  qq[l] = 0;
790  // derive linear equations
791  // mm(p,q)*G(q)=Q(p) for p=n-1=0...k-1, q=n-1=0...k
792  // G(p)=G(p+1,m)
793  for (n = 1; n <= k; n++) {
794  p = n - 1;
795  // derivative id/dt Gtv(n,m)
796  for (j = 0; j <= k; j++) {
797  cweight = cplx_i / h * I.poly_differentiation(n, j);
798  if (j == 0) {
799  for (l = 0; l < sg; l++)
800  qq[p * sg + l] -= cweight * G.tvptr(0, m)[l];
801  } else { // goes into mm(p,q)
802  q = j - 1;
803  for (l = 0; l < sg; l++)
804  mm[sg * (p * k + q) + l] += cweight * one[l];
805  }
806  }
807  // H -- goes into m(p,p)
808  element_set<T, SIZE1>(size1, gtemp, H + sg * n);
809  element_smul<T, SIZE1>(size1, gtemp, -1.0);
810  for (l = 0; l < sg; l++)
811  gtemp[l] += mu * one[l];
812  element_incr<T, SIZE1>(size1, mm + sg * (p * k + p), gtemp);
813  // integral 0..n
814  for (j = 0; j <= k; j++) {
815  cweight = h * I.gregory_weights(n, j);
816  if (j == 0) { // goes into qq(p)
817  element_incr<T, SIZE1>(size1, qq + p * sg, cweight, Sigma.retptr(n, 0),
818  G.tvptr(0, m));
819  } else { // goes into mm(p,q)
820  q = j - 1;
821  if (n >= j) {
822  element_set<T, SIZE1>(size1, stemp, Sigma.retptr(n, j));
823  } else {
824  element_set<T, SIZE1>(size1, stemp, Sigma.retptr(j, n));
825  element_conj<T, SIZE1>(size1, stemp);
826  element_smul<T, SIZE1>(size1, stemp, -1);
827  }
828  for (l = 0; l < sg; l++)
829  mm[sg * (p * k + q) + l] += -cweight * stemp[l];
830  }
831  }
832  // integral Sigmatv*Gmat --> take from Gtv(n,m), write into qq
833  element_incr<T, SIZE1>(size1, qq + p * sg, G.tvptr(n, m));
834  }
835  element_linsolve_right<T, SIZE1>(size1, k, gtemp, mm,
836  qq); // solve kXk problem mm*gtemp=qq
837  // write elements into Gtv
838  for (n = 1; n <= k; n++)
839  element_set<T, SIZE1>(size1, G.tvptr(n, m), gtemp + (n - 1) * sg);
840  }
841  delete[] stemp;
842  delete[] gtemp;
843  delete[] qq;
844  delete[] mm;
845  delete[] one;
846  return;
847 }
848 /*###########################################################################################
849 #
850 # les-FUNCTION: use id_t G(t,nh)= xxx , t=0...n
851 # note: G(t,n) is not continuous in memory!!!
852 # for n<k: start routine
853 #
854 ###########################################################################################*/
856 
892 template <typename T, class GG, int SIZE1>
893 void dyson_timestep_les(int n, GG &G, T mu, std::complex<T> *H, GG &Sigma,
894  integration::Integrator<T> &I, T beta, T h) {
895  typedef std::complex<T> cplx;
896  int k = I.get_k(), k1 = k + 1;
897  int sg, n1, l, m, j, ntau, p, q, sig, size1 = G.size1();
898  cplx *gles, cweight, ih, *mm, *qq, *stemp, cplx_i = cplx(0, 1), *one, *gtemp;
899 
900  n1 = (n > k ? n : k);
901  sg = G.element_size();
902  ntau = G.ntau();
903  sig = G.sig();
904  // check consistency: (more assertations follow in convolution)
905  assert(Sigma.nt() >= n1);
906  assert(G.nt() >= n1);
907  assert(G.sig() == Sigma.sig());
908 
909  one = new cplx[sg];
910  qq = new cplx[k * sg];
911  mm = new cplx[k * k * sg];
912  gtemp = new cplx[sg];
913  stemp = new cplx[sg]; // sic
914  gles = new cplx[(n1 + 1) * sg];
915  for (j = 0; j <= n1; j++)
916  element_set_zero<T, SIZE1>(size1, gles + j * sg);
917  element_set<T, SIZE1>(size1, one, 1.0);
918 
919  // CONVOLUTION SIGMA*G: ---> G^les(j,n)
920  // Note: this is only the tv*vt + les*adv part, Gles is not adressed
921  // Note: gles is determioned for j=0...max(k,n)!!!
922  convolution_timestep_les_tvvt<T, GG, SIZE1>(n, gles, G, Sigma, Sigma, G, G, I, beta,
923  h); // sig needed!!!
924  convolution_timestep_les_lesadv<T, GG, SIZE1>(n, gles, G, Sigma, Sigma, G, G, I, beta,
925  h);
926  // INITIAL VALUE G^les(0,n) = -G^tv(n,0)^*
927  element_conj<T, SIZE1>(size1, gles, G.tvptr(n, 0));
928  element_smul<T, SIZE1>(size1, gles, -1);
929  // Start for integrodifferential equation: j=1...k
930  // .... the usual mess:
931  // derive linear equations
932  // mm(p,q)*G(q)=Q(p) for p=j-1=0...k-1, q=j-1=0...k
933  // G(p)=G(p+1,m)
934  for (l = 0; l < k * k * sg; l++)
935  mm[l] = 0;
936  for (l = 0; l < k * sg; l++)
937  qq[l] = 0;
938  for (j = 1; j <= k; j++) {
939  p = j - 1;
940  // integral Sigmatv*Gvt + Sigma^les*G^adv --> take from gles(j)
941  // after that, gles can be overwritten
942  element_incr<T, SIZE1>(size1, qq + p * sg, gles + j * sg);
943  // derivative id/dt Gtv(n,m)
944  for (m = 0; m <= k; m++) {
945  cweight = cplx_i / h * I.poly_differentiation(j, m);
946  if (m == 0) {
947  for (l = 0; l < sg; l++)
948  qq[p * sg + l] -= cweight * gles[0 * sg + l];
949  } else {
950  q = m - 1;
951  for (l = 0; l < sg; l++)
952  mm[sg * (p * k + q) + l] += cweight * one[l];
953  }
954  }
955  // H -- goes into m(p,p)
956  element_set<T, SIZE1>(size1, gtemp, H + sg * j);
957  element_smul<T, SIZE1>(size1, gtemp, -1);
958  for (l = 0; l < sg; l++)
959  gtemp[l] += mu * one[l];
960  element_incr<T, SIZE1>(size1, mm + sg * (p * k + p), gtemp);
961  // integral Sigma^ret(j,m)G^les(m,n)
962  for (m = 0; m <= k; m++) {
963  cweight = h * I.gregory_weights(j, m);
964  if (m == 0) { // goes into qq(p)
965  element_incr<T, SIZE1>(size1, qq + p * sg, cweight, Sigma.retptr(j, 0),
966  gles + 0 * sg);
967  } else { // goes into mm(p,q)
968  q = m - 1;
969  if (j >= m) {
970  element_set<T, SIZE1>(size1, stemp, Sigma.retptr(j, m));
971  } else {
972  element_set<T, SIZE1>(size1, stemp, Sigma.retptr(m, j));
973  element_conj<T, SIZE1>(size1, stemp);
974  element_smul<T, SIZE1>(size1, stemp, -1);
975  }
976  for (l = 0; l < sg; l++)
977  mm[sg * (p * k + q) + l] += -cweight * stemp[l];
978  }
979  }
980  }
981  element_linsolve_right<T, SIZE1>(size1, k, gles + 1 * sg, mm,
982  qq); // solve kXk problem mm*gtemp=qq
983  // integrodifferential equation k+1...n
984  for (j = k + 1; j <= n; j++) {
985  element_set_zero<T, SIZE1>(size1, mm);
986  // CONTRIBUTION FROM INTEGRAL tv*vt+les*adv
987  element_set<T, SIZE1>(size1, qq, gles + j * sg); // now gles(j) may be overwritten
988  // ACCUMULATE CONTRIBUTION TO id/dt G(j-p,n) p=1...k1 into qq
989  for (p = 1; p <= k1; p++) { // use BD(k+1) !!!
990  cweight = -cplx_i / h * I.bd_weights(p);
991  for (l = 0; l < sg; l++)
992  qq[l] += cweight * gles[sg * (j - p) + l];
993  }
994  for (l = 0; l < sg; l++)
995  mm[l] += cplx_i / h * I.bd_weights(0) * one[l];
996  // CONTRIBUTION FROM H
997  element_set<T, SIZE1>(size1, gtemp, H + sg * j);
998  element_smul<T, SIZE1>(size1, gtemp, -1);
999  for (l = 0; l < sg; l++)
1000  gtemp[l] += mu * one[l];
1001  element_incr<T, SIZE1>(size1, mm, gtemp);
1002  // CONTRIBUTION FROM INTEGRAL Sigma^ret(j,m)*G^les(m,j)
1003  element_set<T, SIZE1>(size1, stemp, Sigma.retptr(j, j));
1004  for (l = 0; l < sg; l++)
1005  mm[l] += -h * stemp[l] * I.gregory_weights(j, j);
1006  for (m = 0; m < j; m++) {
1007  cweight = h * I.gregory_weights(j, m);
1008  element_incr<T, SIZE1>(size1, qq, cweight, Sigma.retptr(j, m), gles + m * sg);
1009  }
1010  element_linsolve_right<T, SIZE1>(size1, gles + j * sg, mm, qq);
1011  }
1012  // write elements into Gles
1013  for (j = 0; j <= n; j++)
1014  element_set<T, SIZE1>(size1, G.lesptr(j, n), gles + j * sg);
1015  delete[] stemp;
1016  delete[] gtemp;
1017  delete[] gles;
1018  delete[] qq;
1019  delete[] mm;
1020  delete[] one;
1021  return;
1022 }
1023 
1025 
1058 template <typename T, class GG, int SIZE1>
1059 void dyson_start_les(GG &G, T mu, std::complex<T> *H, GG &Sigma,
1060  integration::Integrator<T> &I, T beta, T h) {
1061  int k = I.get_k(), n;
1062  for (n = 0; n <= k; n++)
1063  dyson_timestep_les<T, GG, SIZE1>(n, G, mu, H, Sigma, I, beta, h);
1064  return;
1065 }
1067 // preferred: H passed as a cntr::function object:
1068 template <typename T>
1069 void dyson_mat_fourier(herm_matrix<T> &G, herm_matrix<T> &Sigma, T mu, function<T> &H, T beta,
1070  int order) {
1071  int size1 = G.size1();
1072  assert(G.size1() == Sigma.size1());
1073  assert(G.ntau() == Sigma.ntau());
1074  if (size1 == 1)
1075  dyson_mat_fourier_dispatch<T, herm_matrix<T>, 1>(G, Sigma, mu, H.ptr(-1), beta, order);
1076  else
1077  dyson_mat_fourier_dispatch<T, herm_matrix<T>, LARGESIZE>(G, Sigma, mu, H.ptr(-1), beta,
1078  order);
1079 }
1081 template <typename T>
1082 void dyson_mat_fourier(herm_matrix<T> &G, herm_matrix<T> &Sigma, T mu, function<T> &H,
1083  function<T> &SigmaMF, T beta, int order) {
1084  int size1 = G.size1();
1085  assert(G.size1() == Sigma.size1());
1086  assert(G.ntau() == Sigma.ntau());
1087  std::complex<T> *hmf;
1088  hmf = new std::complex<T>[size1*size1];
1089 
1090  if (size1 == 1){
1091  element_set<T, 1>(size1, hmf, H.ptr(-1));
1092  element_incr<T, 1>(size1, hmf, 1.0, SigmaMF.ptr(-1));
1093  dyson_mat_fourier_dispatch<T, herm_matrix<T>, 1>(G, Sigma, mu, hmf, beta, order);
1094  } else {
1095  element_set<T, LARGESIZE>(size1, hmf, H.ptr(-1));
1096  element_incr<T, LARGESIZE>(size1, hmf, 1.0, SigmaMF.ptr(-1));
1097  dyson_mat_fourier_dispatch<T, herm_matrix<T>, LARGESIZE>(G, Sigma, mu, hmf, beta,
1098  order);
1099  }
1100  delete hmf;
1101 }
1103 template <typename T>
1104 void dyson_mat_fixpoint(herm_matrix<T> &G, herm_matrix<T> &Sigma, T mu, function<T> &H,
1105  integration::Integrator<T> &I,T beta, int fixpiter) {
1106  assert(G.size1() == Sigma.size1());
1107  assert(G.ntau() == Sigma.ntau());
1108  int size1 = G.size1();
1109  cdmatrix h0(size1,size1);
1110 
1111  H.get_value(-1,h0);
1112  if (size1 == 1)
1113  dyson_mat_fixpoint_dispatch<T, herm_matrix<T>, 1>(G, Sigma, mu, h0, I, beta, fixpiter);
1114  else
1115  dyson_mat_fixpoint_dispatch<T, herm_matrix<T>, LARGESIZE>(G, Sigma, mu, h0, I, beta,
1116  fixpiter);
1117 }
1119 template <typename T>
1120 void dyson_mat_fixpoint(herm_matrix<T> &G, herm_matrix<T> &Sigma, T mu, function<T> &H,
1121  function<T> &SigmaMF, integration::Integrator<T> &I,T beta, int fixpiter) {
1122  assert(G.size1() == Sigma.size1());
1123  assert(G.ntau() == Sigma.ntau());
1124  int size1 = G.size1();
1125  cdmatrix h0(size1,size1);
1126 
1127  H.get_value(-1,h0);
1128  if (size1 == 1)
1129  dyson_mat_fixpoint_dispatch<T, herm_matrix<T>, 1>(G, Sigma, mu, h0, SigmaMF, I, beta, fixpiter);
1130  else
1131  dyson_mat_fixpoint_dispatch<T, herm_matrix<T>, LARGESIZE>(G, Sigma, mu, h0, SigmaMF, I, beta,
1132  fixpiter);
1133 }
1134 
1135 
1137 template <typename T>
1138 void dyson_mat_steep(herm_matrix<T> &G, herm_matrix<T> &Sigma, T mu, function<T> &H,
1139  integration::Integrator<T> &I,T beta, int maxiter, T tol) {
1140  assert(G.size1() == Sigma.size1());
1141  assert(G.ntau() == Sigma.ntau());
1142  assert(H.size1() == G.size1() );
1143  int size1 = G.size1();
1144  cdmatrix h0(size1,size1);
1145 
1146  H.get_value(-1,h0);
1147  if (size1 == 1)
1148  dyson_mat_steep_dispatch<T, herm_matrix<T>, 1>(G, Sigma, mu, h0, I, beta, maxiter, tol);
1149  else
1150  dyson_mat_steep_dispatch<T, herm_matrix<T>, LARGESIZE>(G, Sigma, mu, h0, I, beta,
1151  maxiter, tol);
1152 }
1153 
1155 template <typename T>
1156 void dyson_mat_steep(herm_matrix<T> &G, herm_matrix<T> &Sigma, T mu, function<T> &H,
1157  function<T> &SigmaMF, integration::Integrator<T> &I,T beta, int maxiter, T tol) {
1158  assert(G.size1() == Sigma.size1());
1159  assert(G.ntau() == Sigma.ntau());
1160  assert(H.size1() == G.size1() );
1161  int size1 = G.size1();
1162  cdmatrix h0(size1,size1);
1163 
1164  H.get_value(-1,h0);
1165  if (size1 == 1)
1166  dyson_mat_steep_dispatch<T, herm_matrix<T>, 1>(G, Sigma, mu, h0, SigmaMF, I, beta, maxiter, tol);
1167  else
1168  dyson_mat_steep_dispatch<T, herm_matrix<T>, LARGESIZE>(G, Sigma, mu, h0, SigmaMF, I, beta,
1169  maxiter, tol);
1170 }
1171 
1173 
1208 template <typename T>
1209 void dyson_mat(herm_matrix<T> &G, herm_matrix<T> &Sigma, T mu, function<T> &H,
1210  function<T> &SigmaMF, integration::Integrator<T> &I,T beta, const int method,
1211  const bool force_hermitian){
1212  assert(method <= 2 && "UNKNOWN CNTR_MAT_METHOD");
1213 
1214 
1215  const int fourier_order = 3;
1216  const double tol=1.0e-12;
1217  int maxiter;
1218  switch(method){
1219  case CNTR_MAT_FOURIER:
1220  dyson_mat_fourier(G, Sigma, mu, H, SigmaMF, beta, fourier_order);
1221  break;
1222  case CNTR_MAT_CG:
1223  maxiter = 40;
1224  dyson_mat_steep(G, Sigma, mu, H, SigmaMF, I, beta, maxiter, tol);
1225  break;
1226  default:
1227  maxiter = 6;
1228  dyson_mat_fixpoint(G, Sigma, mu, H, SigmaMF, I, beta, maxiter);
1229  break;
1230  }
1231  if(force_hermitian){
1232  force_matsubara_hermitian(G);
1233  }
1234 }
1235 
1236 
1237 // global interface
1273 template <typename T>
1274 void dyson_mat(herm_matrix<T> &G, T mu, function<T> &H,
1275  function<T> &SigmaMF, herm_matrix<T> &Sigma, T beta, const int SolveOrder,const int method,
1276  const bool force_hermitian){
1277  assert(method <= 2 && "UNKNOWN CNTR_MAT_METHOD");
1278  assert(SolveOrder <= MAX_SOLVE_ORDER);
1279 
1280  const int fourier_order = 3;
1281  const double tol=1.0e-12;
1282  int maxiter;
1283  switch(method){
1284  case CNTR_MAT_FOURIER:
1285  dyson_mat_fourier(G, Sigma, mu, H, SigmaMF, beta, fourier_order);
1286  break;
1287  case CNTR_MAT_CG:
1288  maxiter = 40;
1289  dyson_mat_steep(G, Sigma, mu, H, SigmaMF, integration::I<T>(SolveOrder), beta, maxiter, tol);
1290  break;
1291  default:
1292  maxiter = 6;
1293  dyson_mat_fixpoint(G, Sigma, mu, H, SigmaMF, integration::I<T>(SolveOrder), beta, maxiter);
1294  break;
1295  }
1296  if(force_hermitian){
1297  force_matsubara_hermitian(G);
1298  }
1299 }
1301 
1333 template <typename T>
1334 void dyson_mat(herm_matrix<T> &G, herm_matrix<T> &Sigma, T mu, function<T> &H,
1335  integration::Integrator<T> &I,T beta, const int method,const bool force_hermitian){
1336  assert(method <= 2 && "UNKNOWN CNTR_MAT_METHOD");
1337 
1338  const int fourier_order = 3;
1339  const double tol=1.0e-12;
1340  int maxiter;
1341 
1342  switch(method){
1343  case 0:
1344  dyson_mat_fourier(G, Sigma, mu, H, beta, fourier_order);
1345  break;
1346  case 1:
1347  maxiter = 40;
1348  dyson_mat_steep(G, Sigma, mu, H, I, beta, maxiter, tol);
1349  break;
1350  case 2:
1351  maxiter=6;
1352  dyson_mat_fixpoint(G, Sigma, mu, H, I, beta, maxiter);
1353  break;
1354  }
1355  if(force_hermitian){
1356  force_matsubara_hermitian(G);
1357  }
1358 }
1359 
1392 template <typename T>
1393 void dyson_mat(herm_matrix<T> &G, T mu, function<T> &H, herm_matrix<T> &Sigma,
1394  T beta, const int SolveOrder, const int method,const bool force_hermitian){
1395  assert(method <= 2 && "UNKNOWN CNTR_MAT_METHOD");
1396  assert(SolveOrder <= MAX_SOLVE_ORDER);
1397 
1398  const int fourier_order = 3;
1399  const double tol=1.0e-12;
1400  int maxiter;
1401 
1402  switch(method){
1403  case 0:
1404  dyson_mat_fourier(G, Sigma, mu, H, beta, fourier_order);
1405  break;
1406  case 1:
1407  maxiter = 40;
1408  dyson_mat_steep(G, Sigma, mu, H, integration::I<T>(SolveOrder), beta, maxiter, tol);
1409  break;
1410  case 2:
1411  maxiter=6;
1412  dyson_mat_fixpoint(G, Sigma, mu, H, integration::I<T>(SolveOrder), beta, maxiter);
1413  break;
1414  }
1415  if(force_hermitian){
1416  force_matsubara_hermitian(G);
1417  }
1418 }
1420 
1451 template <typename T>
1452 void dyson_start(herm_matrix<T> &G, T mu, function<T> &H, herm_matrix<T> &Sigma,
1453  integration::Integrator<T> &I, T beta, T h) {
1454  int size1 = G.size1(), k = I.k();
1455  assert(G.size1() == Sigma.size1());
1456  assert(G.ntau() == Sigma.ntau());
1457  assert(G.nt() >= k);
1458  assert(Sigma.nt() >= k);
1459  if (size1 == 1) {
1460  dyson_start_ret<T, herm_matrix<T>, 1>(G, mu, H.ptr(0), Sigma, I, h);
1461  dyson_start_tv<T, herm_matrix<T>, 1>(G, mu, H.ptr(0), Sigma, I, beta, h);
1462  dyson_start_les<T, herm_matrix<T>, 1>(G, mu, H.ptr(0), Sigma, I, beta, h);
1463  } else {
1464  dyson_start_ret<T, herm_matrix<T>, LARGESIZE>(G, mu, H.ptr(0), Sigma, I, h);
1465  dyson_start_tv<T, herm_matrix<T>, LARGESIZE>(G, mu, H.ptr(0), Sigma, I, beta, h);
1466  dyson_start_les<T, herm_matrix<T>, LARGESIZE>(G, mu, H.ptr(0), Sigma, I, beta, h);
1467  }
1468 }
1469 
1501 template <typename T>
1502 void dyson_start(herm_matrix<T> &G, T mu, function<T> &H, herm_matrix<T> &Sigma,
1503  T beta, T h, const int SolveOrder) {
1504  int size1 = G.size1();
1505  assert(G.size1() == Sigma.size1());
1506  assert(G.ntau() == Sigma.ntau());
1507  assert(G.nt() >= SolveOrder);
1508  assert(Sigma.nt() >= SolveOrder);
1509  if (size1 == 1) {
1510  dyson_start_ret<T, herm_matrix<T>, 1>(G, mu, H.ptr(0), Sigma, integration::I<T>(SolveOrder), h);
1511  dyson_start_tv<T, herm_matrix<T>, 1>(G, mu, H.ptr(0), Sigma, integration::I<T>(SolveOrder), beta, h);
1512  dyson_start_les<T, herm_matrix<T>, 1>(G, mu, H.ptr(0), Sigma, integration::I<T>(SolveOrder), beta, h);
1513  } else {
1514  dyson_start_ret<T, herm_matrix<T>, LARGESIZE>(G, mu, H.ptr(0), Sigma, integration::I<T>(SolveOrder), h);
1515  dyson_start_tv<T, herm_matrix<T>, LARGESIZE>(G, mu, H.ptr(0), Sigma, integration::I<T>(SolveOrder), beta, h);
1516  dyson_start_les<T, herm_matrix<T>, LARGESIZE>(G, mu, H.ptr(0), Sigma, integration::I<T>(SolveOrder), beta, h);
1517  }
1518 }
1520 
1557 template <typename T>
1558 void dyson_timestep(int n, herm_matrix<T> &G, T mu, function<T> &H, herm_matrix<T> &Sigma,
1559  integration::Integrator<T> &I, T beta, T h) {
1560  int size1 = G.size1(), k = I.k();
1561  assert(G.size1() == Sigma.size1());
1562  assert(G.ntau() == Sigma.ntau());
1563  assert(G.nt() >= n);
1564  assert(Sigma.nt() >= n);
1565  assert(n > k);
1566  if (size1 == 1) {
1567  dyson_timestep_ret<T, herm_matrix<T>, 1>(n, G, mu, H.ptr(0), Sigma, I, h);
1568  dyson_timestep_tv<T, herm_matrix<T>, 1>(n, G, mu, H.ptr(n), Sigma, I, beta, h);
1569  dyson_timestep_les<T, herm_matrix<T>, 1>(n, G, mu, H.ptr(0), Sigma, I, beta, h);
1570  } else {
1571  dyson_timestep_ret<T, herm_matrix<T>, LARGESIZE>(n, G, mu, H.ptr(0), Sigma, I, h);
1572  dyson_timestep_tv<T, herm_matrix<T>, LARGESIZE>(n, G, mu, H.ptr(n), Sigma, I, beta,
1573  h);
1574  dyson_timestep_les<T, herm_matrix<T>, LARGESIZE>(n, G, mu, H.ptr(0), Sigma, I, beta,
1575  h);
1576  }
1577 }
1578 
1579 
1617 template <typename T>
1618 void dyson_timestep(int n, herm_matrix<T> &G, T mu, function<T> &H, herm_matrix<T> &Sigma,
1619  T beta, T h, const int SolveOrder) {
1620  int size1 = G.size1();
1621  assert(G.size1() == Sigma.size1());
1622  assert(G.ntau() == Sigma.ntau());
1623  assert(G.nt() >= n);
1624  assert(Sigma.nt() >= n);
1625  assert(n > SolveOrder);
1626  if (size1 == 1) {
1627  dyson_timestep_ret<T, herm_matrix<T>, 1>(n, G, mu, H.ptr(0), Sigma, integration::I<T>(SolveOrder), h);
1628  dyson_timestep_tv<T, herm_matrix<T>, 1>(n, G, mu, H.ptr(n), Sigma, integration::I<T>(SolveOrder), beta, h);
1629  dyson_timestep_les<T, herm_matrix<T>, 1>(n, G, mu, H.ptr(0), Sigma, integration::I<T>(SolveOrder), beta, h);
1630  } else {
1631  dyson_timestep_ret<T, herm_matrix<T>, LARGESIZE>(n, G, mu, H.ptr(0), Sigma, integration::I<T>(SolveOrder), h);
1632  dyson_timestep_tv<T, herm_matrix<T>, LARGESIZE>(n, G, mu, H.ptr(n), Sigma, integration::I<T>(SolveOrder), beta,
1633  h);
1634  dyson_timestep_les<T, herm_matrix<T>, LARGESIZE>(n, G, mu, H.ptr(0), Sigma, integration::I<T>(SolveOrder), beta,
1635  h);
1636  }
1637 }
1639 
1675 template <typename T>
1676 void dyson(herm_matrix<T> &G, T mu, function<T> &H, herm_matrix<T> &Sigma,
1677  integration::Integrator<T> &I, T beta, T h, const int matsubara_method,
1678  const bool force_hermitian) {
1679  int n, k = I.k(), nt = G.nt();
1680  dyson_mat(G, Sigma, mu, H, I, beta, matsubara_method, force_hermitian);
1681  if (nt >= 0)
1682  dyson_start(G, mu, H, Sigma, I, beta, h);
1683  for (n = k + 1; n <= nt; n++)
1684  dyson_timestep(n, G, mu, H, Sigma, I, beta, h);
1685 }
1686 
1687 
1724 template <typename T>
1725 void dyson(herm_matrix<T> &G, T mu, function<T> &H, herm_matrix<T> &Sigma,
1726  T beta, T h, const int SolveOrder, const int matsubara_method,
1727  const bool force_hermitian) {
1728  int n, nt = G.nt();
1729  dyson_mat(G, mu, H, Sigma, beta, SolveOrder, matsubara_method, force_hermitian);
1730  if (nt >= 0)
1731  dyson_start(G, mu, H, Sigma, beta, h, SolveOrder);
1732  for (n = SolveOrder + 1; n <= nt; n++)
1733  dyson_timestep(n, G, mu, H, Sigma, beta, h, SolveOrder);
1734 }
1735 
1736 }
1737 #endif // CNTR_DYSON_IMPL_H
cntr::green_from_H
void green_from_H(herm_matrix< T > &G, T mu, cdmatrix &eps, T beta, T h)
Propagator for time-independent free Hamiltonian
Definition: cntr_equilibrium_impl.hpp:1142
integration::Integrator
Class Integrator contains all kinds of weights for integration and differentiation of a function at ...
Definition: integration.hpp:128
cntr::force_matsubara_hermitian
void force_matsubara_hermitian(herm_matrix< T > &G)
Force the Matsubara component of herm_matrix to be a hermitian matrix at each . ...
Definition: cntr_utilities_impl.hpp:2077
cntr_equilibrium_decl.hpp
cntr_elements.hpp
cntr_function_decl.hpp
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
integration::I
Integrator< T > & I(int k)
Definition: integration.hpp:506
cntr_dyson_decl.hpp
cntr_utilities_decl.hpp
cntr::function
Class function for objects with time on real axis.
Definition: cntr_convolution_decl.hpp:9
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::herm_matrix::ntau
int ntau(void) const
Definition: cntr_herm_matrix_decl.hpp:70
cntr_vie2_decl.hpp
cntr::herm_matrix::size1
int size1(void) const
Definition: cntr_herm_matrix_decl.hpp:68