NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_dyson_omp_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_DYSON_OMP_IMPL_H
2 #define CNTR_DYSON_OMP_IMPL_H
3 
4 #include "cntr_dyson_omp_decl.hpp"
5 //#include "cntr_exception.hpp"
6 #include "cntr_elements.hpp"
7 #include "cntr_dyson_decl.hpp"
8 #include "cntr_dyson_impl.hpp"
9 #include "cntr_function_decl.hpp"
10 #include "cntr_herm_matrix_decl.hpp"
11 #include "cntr_herm_pseudo_decl.hpp"
12 #include "cntr_pseudo_convolution_decl.hpp"
13 #include "cntr_pseudodyson_decl.hpp"
14 #include "cntr_convolution_impl.hpp"
15 
16 namespace cntr {
17 
18 #if CNTR_USE_OMP == 1
19 
20 /* #######################################################################################
21 # DYSON EQUATION:
22 #
23 # [ id/dt + mu - H(t) ] G(t,t') - [Sigma*G](t,t') = 1(t,t')
24 # G(t,t')[-id/dt' + mu - H(t')] - [G*Sigma](t,t') = 1(t,t')
25 #
26 # explanation see cntr_dyson
27 # for omp paralellizaion I use conjugate equation, but this does not work for
28 #
29 #
30 ###########################################################################################*/
31 
32 /*###########################################################################################
33 # RETARDED FUNCTION: (GG = herm_matrix or herm_pseudo)
34 ###########################################################################################*/
36 template <typename T, class GG, int SIZE1>
37 void dyson_timestep_ret_omp(int omp_num_threads, int n, GG &G, T mu, std::complex<T> *H,
38  GG &Sigma, integration::Integrator<T> &I, T h) {
39  typedef std::complex<T> cplx;
40  int k = I.get_k(), k1 = k + 1;
41  cplx cplx_i = cplx(0, 1);
42  int size1 = G.size1();
43  int sg = G.element_size();
44  // check consistency:
45  assert(k + 1<= n);
46  assert(n<= Sigma.nt());
47  assert(n<= G.nt());
48  assert(G.sig()== Sigma.sig());
49 
51  // SET ENTRIES IN TIMESTEP TO 0
52  {
53  cplx *gret = G.retptr(n, 0);
54  int n1 = (n + 1) * sg, i;
55  for (i = 0; i < n1; i++)
56  gret[i] = 0;
57  }
59  // INITIAL VALUE t' = n
60  element_set<T, SIZE1>(size1, G.retptr(n, n), -cplx_i);
62  // START VALUES t' = n-j, j = 1...k: solve a kxk problem
63  // this is the same as before, i.e., use conjugate equation ii*d/dt1 G(t,t1) = ...
64  // using the other equation ii*d/dt G(t,t1) = ... seems to be unstable
65  // (may be interesting to investigate this instability in general)
66  {
67  int i, j, p, l, q;
68  cplx w0, cweight;
69  T weight;
70  cplx *gtemp = new cplx[k * sg];
71  cplx *diffw = new cplx[k1 + 1];
72  cplx *qq = new cplx[(n + 1) * sg];
73  cplx *one = new cplx[sg];
74  cplx *mm = new cplx[k * k * sg];
75  cplx *hj = new cplx[sg];
76  cplx *stemp = new cplx[sg]; // sic
77  element_set<T, SIZE1>(size1, one, 1);
78  for (i = 0; i < k * k * sg; i++)
79  mm[i] = 0;
80  for (i = 0; i < k * sg; i++)
81  qq[i] = 0;
82  for (j = 1; j <= k; j++) {
83  p = j - 1;
84  // derivatives:
85  for (l = 0; l <= k; l++) {
86  cweight = cplx_i / h * I.poly_differentiation(j, l);
87  if (l == 0) {
88  element_incr<T, SIZE1>(size1, qq + p * sg, -cweight, G.retptr(n, n));
89  } else {
90  q = l - 1;
91  element_incr<T, SIZE1>(size1, mm + sg * (p * k + q), cweight);
92  }
93  }
94  // H
95  element_set<T, SIZE1>(size1, gtemp, H + (n - j) * sg);
96  element_smul<T, SIZE1>(size1, gtemp, -1);
97  for (i = 0; i < sg; i++)
98  gtemp[i] += mu * one[i];
99  element_incr<T, SIZE1>(size1, mm + sg * (p + k * p), gtemp);
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,
106  Sigma.retptr(n - l, n - j)); // stemp is element of type G!!
107  } else {
108  element_set<T, SIZE1>(size1, stemp, Sigma.retptr(n - j, n - l));
109  element_conj<T, SIZE1>(size1, stemp);
110  weight *= -1;
111  }
112  if (l == 0) {
113  element_incr<T, SIZE1>(size1, qq + p * sg, weight, G.retptr(n, n),
114  stemp);
115  } else {
116  q = l - 1;
117  element_incr<T, SIZE1>(size1, mm + sg * (p * k + q), -weight, stemp);
118  }
119  }
120  }
121  element_linsolve_left<T, SIZE1>(size1, k, gtemp, mm, qq); // gtemp * mm = qq
122  for (j = 1; j <= k; j++)
123  element_set<T, SIZE1>(size1, G.retptr(n, n - j), gtemp + (j - 1) * sg);
124  delete[] gtemp;
125  delete[] diffw;
126  delete[] qq;
127  delete[] one;
128  delete[] mm;
129  delete[] hj;
130  delete[] stemp;
131  }
133 // now use equation ii*d/dt G(t,t1) = ... to compute G(n*h,j*h),j=0 ... n-k-1
134 // OMP parallelization over j
135 
136 #pragma omp parallel num_threads(omp_num_threads)
137  {
138  // convolution Sigma*G ->> written to G, on
139  int j, p, i;
140  int nomp = omp_get_num_threads();
141  int tid = omp_get_thread_num();
142  std::vector<bool> mask_ret(n + 1, false);
143  cplx w0 = h * I.gregory_omega(0);
144  cplx *diffw = new cplx[k1 + 1];
145  cplx *qq = new cplx[sg];
146  cplx *mm = new cplx[sg];
147  for (i = 0; i < n - k; i++)
148  if (i % nomp == tid)
149  mask_ret[i] = true;
150  incr_convolution_ret<T, GG, SIZE1>(n, mask_ret, cplx(1.0, 0.0), G, Sigma, Sigma,
151  NULL, G, G, I, h);
152  for (p = 0; p <= k1; p++)
153  diffw[p] = I.bd_weights(p) * cplx_i / h; // use BD(k+1!!)
154  for (j = 0; j < n - k; j++) {
155  if (mask_ret[j]) {
156  element_set<T, SIZE1>(size1, qq, G.retptr(n, j)); // << Sigma*G(n,j)
157  // set up mm and qqj for 1x1 problem:
158  for (p = 1; p <= k1; p++)
159  element_incr<T, SIZE1>(size1, qq, -diffw[p], G.retptr(n - p, j));
160  element_set<T, SIZE1>(size1, mm, diffw[0] + mu);
161  element_incr<T, SIZE1>(size1, mm, -w0, Sigma.retptr(j, j));
162  element_incr<T, SIZE1>(size1, mm, cplx(-1.0, 0.0), H + n * sg);
163  element_linsolve_right<T, SIZE1>(size1, G.retptr(n, j), mm, qq); // mm*G=qq
164  }
165  }
166  delete[] diffw;
167  delete[] qq;
168  delete[] mm;
169  }
170  return;
171 }
172 // retarded start is not omp-paralellized
173 /*###########################################################################################
174 # TV FUNCTION:
175 ###########################################################################################*/
176 // GG = herm_matrix
178 template <typename T, class GG, int SIZE1>
179 void dyson_timestep_tv_omp(int omp_num_threads, int n, GG &G, T mu, std::complex<T> *Hn,
180  GG &Sigma, integration::Integrator<T> &I, T beta, T h) {
181  typedef std::complex<T> cplx;
182  int size1 = G.size1();
183  int k = I.get_k(), k1 = k + 1;
184  int sg = G.element_size();
185  int ntau = G.ntau();
186  cplx ih = cplx(0, 1.0 / h);
187  assert(k + 1<= n);
188  assert(n<= Sigma.nt());
189  assert(n<= G.nt());
190  assert(G.sig()== Sigma.sig());
191 
192  for (int j = 0; j <= ntau; j++)
193  element_set_zero<T, SIZE1>(size1, G.tvptr(n, j));
194 #pragma omp parallel num_threads(omp_num_threads)
195  {
196  // convolution Sigma*G ->> written to G.tv
197  int j, p, i;
198  int nomp = omp_get_num_threads();
199  int tid = omp_get_thread_num();
200  std::vector<bool> mask(ntau + 1, false);
201  cplx cweight;
202  cplx *diffw = new cplx[k1 + 1];
203  cplx *qq = new cplx[sg];
204  cplx *mm = new cplx[sg];
205  for (i = 0; i <= ntau; i++)
206  if (i % nomp == tid)
207  mask[i] = true;
208  incr_convolution_tv<T, GG, SIZE1>(n, mask, cplx(1.0, 0.0), G, Sigma, Sigma, NULL,
209  NULL, G, G, I, beta, h);
210  // Now solve
211  // [ i/h bd(0) - H - h w(n,0) Sigma(n,n) ] G(n,m) = Q(m),
212  // where Q is initially stored in G(n,m)
213  element_set<T, SIZE1>(size1, mm, ih * I.bd_weights(0) + mu);
214  element_incr<T, SIZE1>(size1, mm, cplx(-1.0, 0.0), Hn);
215  cweight = -h * I.gregory_weights(n, 0);
216  element_incr<T, SIZE1>(size1, mm, cweight, Sigma.retptr(n, n));
217  // ACCUMULATE CONTRIBUTION TO id/dt G(t,t') FROM t=mh, m=n-k..n-1
218  for (p = 0; p <= k1; p++)
219  diffw[p] = ih * I.bd_weights(p); // use BD(k+1!!)
220  for (j = 0; j <= ntau; j++) {
221  if (mask[j]) {
222  element_set<T, SIZE1>(size1, qq, G.tvptr(n, j));
223  for (p = 1; p <= k1; p++)
224  element_incr<T, SIZE1>(size1, qq, -diffw[p], G.tvptr(n - p, j));
225  element_linsolve_right<T, SIZE1>(size1, G.tvptr(n, j), mm, qq);
226  }
227  }
228  delete[] qq;
229  delete[] mm;
230  delete[] diffw;
231  }
232  return;
233 }
234 // GG = pseudo_matrix the only differebnce is the convolution!
236 template <typename T, class GG, int SIZE1>
237 void pseudodyson_timestep_tv_omp(int omp_num_threads, int n, GG &G, T mu,
238  std::complex<T> *Hn, GG &Sigma,
239  integration::Integrator<T> &I, T beta, T h) {
240  typedef std::complex<T> cplx;
241  int size1 = G.size1();
242  int k = I.get_k(), k1 = k + 1;
243  int sg = G.element_size();
244  int ntau = G.ntau();
245  cplx ih = cplx(0, 1.0 / h);
246  assert(k + 1<= n);
247  assert(n<= Sigma.nt());
248  assert(n<= G.nt());
249  assert(G.sig()== Sigma.sig());
250 
251  for (int j = 0; j <= ntau; j++)
252  element_set_zero<T, SIZE1>(size1, G.tvptr(n, j));
253 #pragma omp parallel num_threads(omp_num_threads)
254  {
255  // convolution Sigma*G ->> written to G.tv
256  int j, p, i;
257  int nomp = omp_get_num_threads();
258  int tid = omp_get_thread_num();
259  std::vector<bool> mask(ntau + 1, false);
260  cplx cweight;
261  cplx *diffw = new cplx[k1 + 1];
262  cplx *qq = new cplx[sg];
263  cplx *mm = new cplx[sg];
264  for (i = 0; i <= ntau; i++)
265  if (i % nomp == tid)
266  mask[i] = true;
267  incr_pseudo_convolution_tv<T, GG, SIZE1>(n, mask, cplx(1.0, 0.0), G, Sigma, Sigma,
268  NULL, NULL, G, G, I, beta, h);
269  // Now solve
270  // [ i/h bd(0) - H - h w(n,0) Sigma(n,n) ] G(n,m) = Q(m),
271  // where Q is initially stored in G(n,m)
272  element_set<T, SIZE1>(size1, mm, ih * I.bd_weights(0) + mu);
273  element_incr<T, SIZE1>(size1, mm, cplx(-1.0, 0.0), Hn);
274  cweight = -h * I.gregory_weights(n, 0);
275  element_incr<T, SIZE1>(size1, mm, cweight, Sigma.retptr(n, n));
276  // ACCUMULATE CONTRIBUTION TO id/dt G(t,t') FROM t=mh, m=n-k..n-1
277  for (p = 0; p <= k1; p++)
278  diffw[p] = ih * I.bd_weights(p); // use BD(k+1!!)
279  for (j = 0; j <= ntau; j++) {
280  if (mask[j]) {
281  element_set<T, SIZE1>(size1, qq, G.tvptr(n, j));
282  for (p = 1; p <= k1; p++)
283  element_incr<T, SIZE1>(size1, qq, -diffw[p], G.tvptr(n - p, j));
284  element_linsolve_right<T, SIZE1>(size1, G.tvptr(n, j), mm, qq);
285  }
286  }
287  delete[] qq;
288  delete[] mm;
289  delete[] diffw;
290  }
291  return;
292 }
293 /*###########################################################################################
294 # LESSER FUNCTION: (GG = herm_matrix or herm_pseudo)
295 ###########################################################################################*/
297 template <typename T, class GG, int SIZE1>
298 void dyson_timestep_les_omp(int omp_num_threads, int n, GG &G, T mu, std::complex<T> *H,
299  GG &Sigma, integration::Integrator<T> &I, T beta, T h) {
300  typedef std::complex<T> cplx;
301  cplx cplx_i = cplx(0, 1);
302  int k = I.get_k(), k1 = k + 1;
303  int size1 = G.size1();
304  int sg = G.element_size();
305  int n1 = (n > k ? n : k);
307  // check consistency: (more assertations follow in convolution)
308  assert(k + 1<= n);
309  assert(G.ntau()<= Sigma.ntau());
310  assert(n1<= Sigma.nt());
311  assert(n1<= G.nt());
312  assert(G.sig()== Sigma.sig());
313  // OMP PARALELLIZARION STARTS ONLY FOR n>=2*k+1
314  if (n < 2 * k + 1) {
315  return dyson_timestep_les<T, GG, SIZE1>(n, G, mu, H, Sigma, I, beta, h);
316  }
317  for (int j = 0; j <= n; j++)
318  element_set_zero<T, SIZE1>(size1, G.lesptr(j, n));
320 // get G(j,n), j=0...n-k-1 from d/dt' G(t,t') equation
321 #pragma omp parallel num_threads(omp_num_threads)
322  {
323  int j, p, i;
324  int nomp = omp_get_num_threads();
325  int tid = omp_get_thread_num();
326  std::vector<bool> mask_les(n + 1, false);
327  cplx w0 = h * I.gregory_omega(0);
328  cplx *diffw = new cplx[k1 + 1];
329  cplx *qq = new cplx[sg];
330  cplx *mm = new cplx[sg];
331  cplx *stemp = new cplx[sg];
332  // convolution Sigma*G ->> written to G
333  for (i = 0; i < n - k; i++)
334  if (i % nomp == tid)
335  mask_les[i] = true;
336  incr_convolution_les<T, GG, SIZE1>(n, mask_les, cplx(1.0, 0.0), G, G, G, NULL, NULL,
337  Sigma, Sigma, I, beta, h);
338  for (p = 0; p <= k1; p++)
339  diffw[p] = I.bd_weights(p) * cplx_i / h; // use BD(k+1!!)
340  for (j = 0; j < n - k; j++) {
341  if (mask_les[j]) {
342  element_set<T, SIZE1>(size1, qq, G.lesptr(j, n)); // << G*Sigma(j,n)
343  // set up mm and qqj for 1x1 problem:
344  for (p = 1; p <= k1; p++)
345  element_incr<T, SIZE1>(size1, qq, diffw[p], G.lesptr(j, n - p));
346  element_set<T, SIZE1>(size1, mm, -diffw[0] + mu);
347  element_conj<T, SIZE1>(size1, stemp, Sigma.retptr(j, j));
348  element_incr<T, SIZE1>(size1, mm, -w0, stemp);
349  element_incr<T, SIZE1>(size1, mm, cplx(-1.0, 0.0), H + n * sg);
350  element_linsolve_left<T, SIZE1>(size1, G.lesptr(j, n), mm, qq);
351  }
352  }
353  delete[] diffw;
354  delete[] qq;
355  delete[] mm;
356  delete[] stemp;
357  }
359  // get G(j,n), j=n-k...n from d/dt G(t,t') equation (old implementation)
360  // currently not paralellized
361  {
362  int j, p, m;
363  cplx *gles = new cplx[(n + 1) * sg];
364  cplx *qq = new cplx[k * sg];
365  cplx *mm = new cplx[k * k * sg];
366  cplx cweight;
367 // CONVOLUTION SIGMA*G: ---> G^les(j,n) j=n-k...n
368 // Note: this is only the tv*vt + les*adv part, Gles is not adressed
369 // coupld parallelize this:
370 #pragma omp parallel num_threads(omp_num_threads)
371  {
372  int j1;
373  int nomp = omp_get_num_threads();
374  int tid = omp_get_thread_num();
375  for (j1 = n - k; j1 <= n; j1++) {
376  if ((n - j1) % nomp == tid) {
377  element_set_zero<T, SIZE1>(size1, gles + j1 * sg);
378  convolution_timestep_les_tvvt<T, GG, SIZE1>(n, j1, j1, gles, G, Sigma,
379  Sigma, G, G, I, beta, h);
380  convolution_timestep_les_lesadv<T, GG, SIZE1>(n, j1, j1, gles, G, Sigma,
381  Sigma, G, G, I, beta, h);
382  }
383  }
384  }
385  for (j = n - k; j <= n; j++) {
386  // CONTRIBUTION FROM INTEGRAL tv*vt+les*adv
387  element_set<T, SIZE1>(size1, qq, gles + j * sg);
388  // ACCUMULATE CONTRIBUTION TO id/dt G(j-p,n) p=1...k1 into qq
389  for (p = 1; p <= k1; p++) { // use BD(k+1) !!!
390  cweight = -cplx_i / h * I.bd_weights(p);
391  element_incr<T, SIZE1>(size1, qq, cweight, G.lesptr(j - p, n));
392  }
393  element_set<T, SIZE1>(size1, mm, cplx_i / h * I.bd_weights(0) + mu);
394  element_incr<T, SIZE1>(size1, mm, cplx(-1.0, 0.0), H + sg * j);
395  cweight = -h * I.gregory_weights(j, j);
396  element_incr<T, SIZE1>(size1, mm, cweight, Sigma.retptr(j, j));
397  for (m = 0; m < j; m++) {
398  cweight = h * I.gregory_weights(j, m);
399  element_incr<T, SIZE1>(size1, qq, cweight, Sigma.retptr(j, m),
400  G.lesptr(m, n));
401  }
402  element_linsolve_right<T, SIZE1>(size1, G.lesptr(j, n), mm, qq);
403  }
404  delete[] gles;
405  delete[] qq;
406  delete[] mm;
407  }
408  return;
409 }
411 
413 // main implementation:
414 // with function object:
416 template <typename T>
417 void pseudodyson_timestep_omp(int omp_num_threads, int n, herm_pseudo<T> &G, T lam0,
418  function<T> &H, herm_pseudo<T> &Sigma,
419  integration::Integrator<T> &I, T beta, T h) {
420  int size1 = G.size1(), k = I.k();
421  int omp_num_threads1 = (omp_num_threads == -1 ? omp_get_max_threads() : omp_num_threads);
422  assert(k + 1<= n);
423  assert(n<= Sigma.nt());
424  assert(n<= G.nt());
425  assert(n<= H.nt());
426  assert(G.sig()== Sigma.sig());
427  assert(G.size1()== Sigma.size1());
428  assert(G.size1()== H.size1());
429  assert(G.ntau()== Sigma.ntau());
430  if (size1 == 1) {
431  dyson_timestep_ret_omp<T, herm_pseudo<T>, 1>(omp_num_threads1, n, G, lam0, H.ptr(0),
432  Sigma, I, h);
433  pseudodyson_timestep_tv_omp<T, herm_pseudo<T>, 1>(omp_num_threads1, n, G, lam0,
434  H.ptr(n), Sigma, I, beta, h);
435  dyson_timestep_les_omp<T, herm_pseudo<T>, 1>(omp_num_threads1, n, G, lam0, H.ptr(0),
436  Sigma, I, beta, h);
437  } else {
438  dyson_timestep_ret_omp<T, herm_pseudo<T>, LARGESIZE>(omp_num_threads1, n, G, lam0,
439  H.ptr(0), Sigma, I, h);
440  pseudodyson_timestep_tv_omp<T, herm_pseudo<T>, LARGESIZE>(
441  omp_num_threads1, n, G, lam0, H.ptr(n), Sigma, I, beta, h);
442  dyson_timestep_les_omp<T, herm_pseudo<T>, LARGESIZE>(omp_num_threads1, n, G, lam0,
443  H.ptr(0), Sigma, I, beta, h);
444  }
445 }
446 // with raw pointer:
448 template <typename T>
449 void pseudodyson_timestep_omp(int omp_num_threads, int n, herm_pseudo<T> &G, T lam0,
450  std::complex<T> *Ht, herm_pseudo<T> &Sigma,
451  integration::Integrator<T> &I, T beta, T h) {
452  int size1 = G.size1(), k = I.k();
453  int omp_num_threads1 = (omp_num_threads == -1 ? omp_get_max_threads() : omp_num_threads);
454  assert(k + 1<= n);
455  assert(n<= Sigma.nt());
456  assert(n<= G.nt());
457  assert(G.sig()== Sigma.sig());
458  assert(G.size1()== Sigma.size1());
459  assert(G.ntau()== Sigma.ntau());
460  if (size1 == 1) {
461  dyson_timestep_ret_omp<T, herm_pseudo<T>, 1>(omp_num_threads1, n, G, lam0, Ht, Sigma,
462  I, h);
463  pseudodyson_timestep_tv_omp<T, herm_pseudo<T>, 1>(
464  omp_num_threads1, n, G, lam0, Ht + n * size1 * size1, Sigma, I, beta, h);
465  dyson_timestep_les_omp<T, herm_pseudo<T>, 1>(omp_num_threads1, n, G, lam0, Ht, Sigma,
466  I, beta, h);
467  } else {
468  dyson_timestep_ret_omp<T, herm_pseudo<T>, LARGESIZE>(omp_num_threads1, n, G, lam0,
469  Ht, Sigma, I, h);
470  pseudodyson_timestep_tv_omp<T, herm_pseudo<T>, LARGESIZE>(
471  omp_num_threads1, n, G, lam0, Ht + n * size1 * size1, Sigma, I, beta, h);
472  dyson_timestep_les_omp<T, herm_pseudo<T>, LARGESIZE>(omp_num_threads1, n, G, lam0,
473  Ht, Sigma, I, beta, h);
474  }
475 }
476 // only for compatibility, works for size1 only
478 template <typename T>
479 void pseudodyson_timestep_omp(int omp_num_threads, int n, herm_pseudo<T> &G, T lam0,
480  std::vector<std::complex<T>> &Ht, herm_pseudo<T> &Sigma,
481  integration::Integrator<T> &I, T beta, T h) {
482  int size1 = G.size1(), k = I.k();
483  std::complex<T> *hh = new std::complex<T>[n + 1];
484  int omp_num_threads1 = (omp_num_threads == -1 ? omp_get_max_threads() : omp_num_threads);
485  assert(k + 1<= n);
486  assert(size1== 1);
487  assert(n<= Sigma.nt());
488  assert(n<= G.nt());
489  assert(G.sig()== Sigma.sig());
490  assert(G.size1()== Sigma.size1());
491  assert(G.ntau()== Sigma.ntau());
492  assert(n + 1<= (int)Ht.size());
493  for (int j = 0; j <= n; j++)
494  hh[j] = Ht[j];
495  dyson_timestep_ret_omp<T, herm_pseudo<T>, 1>(omp_num_threads1, n, G, lam0, hh, Sigma, I,
496  h);
497  pseudodyson_timestep_tv_omp<T, herm_pseudo<T>, 1>(
498  omp_num_threads1, n, G, lam0, hh + n * size1 * size1, Sigma, I, beta, h);
499  dyson_timestep_les_omp<T, herm_pseudo<T>, 1>(omp_num_threads1, n, G, lam0, hh, Sigma, I,
500  beta, h);
501  delete[] hh;
502 }
503 
505 template <typename T>
506 void dyson_timestep_omp(int omp_num_threads, int n, herm_matrix<T> &G, T lam0,
507  function<T> &H, herm_matrix<T> &Sigma, integration::Integrator<T> &I,
508  T beta, T h) {
509  int size1 = G.size1(), k = I.k();
510  int omp_num_threads1 = (omp_num_threads == -1 ? omp_get_max_threads() : omp_num_threads);
511  assert(k + 1<= n);
512  assert(n<= Sigma.nt());
513  assert(n<= G.nt());
514  assert(n<= H.nt());
515  assert(G.sig()== Sigma.sig());
516  assert(G.size1()== Sigma.size1());
517  assert(G.size1()== H.size1());
518  assert(G.ntau()== Sigma.ntau());
519  if (size1 == 1) {
520  dyson_timestep_ret_omp<T, herm_matrix<T>, 1>(omp_num_threads1, n, G, lam0, H.ptr(0),
521  Sigma, I, h);
522  dyson_timestep_tv_omp<T, herm_matrix<T>, 1>(omp_num_threads1, n, G, lam0, H.ptr(n),
523  Sigma, I, beta, h);
524  dyson_timestep_les_omp<T, herm_matrix<T>, 1>(omp_num_threads1, n, G, lam0, H.ptr(0),
525  Sigma, I, beta, h);
526  } else {
527  dyson_timestep_ret_omp<T, herm_matrix<T>, LARGESIZE>(omp_num_threads1, n, G, lam0,
528  H.ptr(0), Sigma, I, h);
529  dyson_timestep_tv_omp<T, herm_matrix<T>, LARGESIZE>(omp_num_threads1, n, G, lam0,
530  H.ptr(n), Sigma, I, beta, h);
531  dyson_timestep_les_omp<T, herm_matrix<T>, LARGESIZE>(omp_num_threads1, n, G, lam0,
532  H.ptr(0), Sigma, I, beta, h);
533  }
534 }
535 
576 template <typename T>
577 void dyson_timestep_omp(int omp_num_threads, int n, herm_matrix<T> &G, T lam0,
578  function<T> &H, herm_matrix<T> &Sigma,
579  T beta, T h, int SolveOrder) {
580  int size1 = G.size1();
581  int omp_num_threads1 = (omp_num_threads == -1 ? omp_get_max_threads() : omp_num_threads);
582  assert(SolveOrder + 1<= n);
583  assert(n<= Sigma.nt());
584  assert(n<= G.nt());
585  assert(n<= H.nt());
586  assert(G.sig()== Sigma.sig());
587  assert(G.size1()== Sigma.size1());
588  assert(G.size1()== H.size1());
589  assert(G.ntau()== Sigma.ntau());
590  if (size1 == 1) {
591  dyson_timestep_ret_omp<T, herm_matrix<T>, 1>(omp_num_threads1, n, G, lam0, H.ptr(0),
592  Sigma, integration::I<T>(SolveOrder), h);
593  dyson_timestep_tv_omp<T, herm_matrix<T>, 1>(omp_num_threads1, n, G, lam0, H.ptr(n),
594  Sigma, integration::I<T>(SolveOrder), beta, h);
595  dyson_timestep_les_omp<T, herm_matrix<T>, 1>(omp_num_threads1, n, G, lam0, H.ptr(0),
596  Sigma, integration::I<T>(SolveOrder), beta, h);
597  } else {
598  dyson_timestep_ret_omp<T, herm_matrix<T>, LARGESIZE>(omp_num_threads1, n, G, lam0,
599  H.ptr(0), Sigma, integration::I<T>(SolveOrder), h);
600  dyson_timestep_tv_omp<T, herm_matrix<T>, LARGESIZE>(omp_num_threads1, n, G, lam0,
601  H.ptr(n), Sigma, integration::I<T>(SolveOrder), beta, h);
602  dyson_timestep_les_omp<T, herm_matrix<T>, LARGESIZE>(omp_num_threads1, n, G, lam0,
603  H.ptr(0), Sigma, integration::I<T>(SolveOrder), beta, h);
604  }
605 }
606 
607 #endif // CNTR_USE_OMP
608 
609 } // namespace cntr
610 
611 #endif // CNTR_DYSON_OMP_IMPL_H
integration::Integrator
Class Integrator contains all kinds of weights for integration and differentiation of a function at ...
Definition: integration.hpp:128
cntr::function::size1
int size1(void) const
Definition: cntr_function_decl.hpp:42
cntr_herm_pseudo_decl.hpp
cntr_elements.hpp
cntr_function_decl.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_dyson_impl.hpp
cntr::herm_matrix::sig
int sig(void) const
Definition: cntr_herm_matrix_decl.hpp:72
cntr_pseudodyson_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::function::nt
int nt(void) const
Definition: cntr_function_decl.hpp:44
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_dyson_omp_decl.hpp
cntr_pseudo_convolution_decl.hpp
cntr::herm_matrix::size1
int size1(void) const
Definition: cntr_herm_matrix_decl.hpp:68