NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_convolution_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_CONVOLUTION_IMPL_H
2 #define CNTR_CONVOLUTION_IMPL_H
3 
4 #include "eigen_map.hpp"
5 #include "cntr_convolution_decl.hpp"
6 #include "cntr_elements.hpp"
7 #include "cntr_function_decl.hpp"
8 #include "cntr_herm_matrix_decl.hpp"
9 
10 namespace cntr {
11 
12 /* #######################################################################################
13 #
14 # CONTOUR CONVOLUTION C = A*B
15 #
16 # (0) The routine computes C(t,t') = int_CC dt1 A(t,t1)B(t1,t') at timestep n,
17 # or matsubara convolution. i.e.,
18 # C^ret(nh,t'<=nh), C^les(t<=nh,nh), C^tv(nt,tau=0..beta). If the timestep n==-1,
19 # the Matsubara convolution is performed, and C^ret,C^tv, C^les are untouched
20 #
21 # (1) A and B are contour Greenfunctions of type GA and GB, which may be scalar,
22 # matrix etc. The behavior of the elements for the following operations must
23 # be specified:
24 # - multiplication: routine element_incr<T,GC,GA,GB>
25 # - conjugation: routine element_conjugate<T,G>
26 #
27 # (2) Acc and Bcc are the conjugate functions to A and B. If A or B are hermitian,
28 # just call convolution with A=Acc or B=Bcc, respectively.
29 #
30 # (3) For the computation of timestep n, A(t,t') and B(t,t') are adressed at times
31 # t,t' <= max(n,k), where k is the Integartion order (see Integrator). I.e.,
32 # the timesteps n=0..k can be computed only if A and B are given for t,t'<=k.
33 #
34 ###########################################################################################*/
35 
36 // MATSUBARA INTEGRAL:
37 //
38 // the integral c = int_0^beta dx a(tau-x)b(x) for matsubara_integral_1
39 // the integral c = int_0^beta dx a(x)b(x-tau) for matsubara_integral_2
40 //
41 // for beta=ntau, tau=x,
42 // a and b are antiperiodic, i.e., a(-x+beta) =-a(-x) for 0 <= x <= beta
43 //
44 // a,b,c are square-matrix-vauled, dimension size1, element size
45 // sa=size1*size1
46 // a is a pointer to an array of ntau elements of size sa*cplx: a + i*sa ->
47 // a(i)
48 // b is a pointer to an array of ntau elements of size sb*cplx: b + i*sa ->
49 // b(i)
50 // c is a pointer to an array of one element of size sa*cplx
51 
52 
54 
85 template <typename T, int SIZE1>
86 void matsubara_integral_1(int size1, int m, int ntau, std::complex<T> *C,
87  std::complex<T> *A, std::complex<T> *B,
88  integration::Integrator<T> &I, int sig) {
89  typedef std::complex<T> cplx;
90  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1, j, l, sa1, sb1, sc1;
91  cplx *amat, *bmat, *ctemp1, *ctemp2;
92  T weight;
93 
94  sa1 = size1 * size1;
95  sb1 = sa1;
96  sc1 = sa1;
97  // std::cout << "conv " << m << " " << ntau << std::endl;
98  ctemp1 = new cplx[sc1];
99  ctemp2 = new cplx[sc1];
100  // CONTRIBUTION FROM 0...TAU
101  for (l = 0; l < sc1; l++) {
102  ctemp1[l] = 0;
103  }
104  if (m >= k2 - 1) {
105  amat = A + m * sa1;
106  bmat = B;
107  for (j = 0; j <= k; j++) {
108  weight = I.gregory_omega(j);
109  element_incr<T, SIZE1>(size1, ctemp1, weight, amat, bmat);
110  amat -= sa1;
111  bmat += sb1;
112  }
113  for (j = k1; j < m - k; j++) {
114  element_incr<T, SIZE1>(size1, ctemp1, amat, bmat);
115  amat -= sa1;
116  bmat += sb1;
117  }
118  for (j = m - k; j <= m; j++) {
119  weight = I.gregory_omega(m - j);
120  element_incr<T, SIZE1>(size1, ctemp1, weight, amat, bmat);
121  amat -= sa1;
122  bmat += sb1;
123  }
124  } else if (m >= k) {
125  amat = A + m * sa1;
126  bmat = B;
127  for (j = 0; j <= m; j++) {
128  weight = I.gregory_weights(m, j);
129  element_incr<T, SIZE1>(size1, ctemp1, weight, amat, bmat);
130  amat -= sa1;
131  bmat += sb1;
132  }
133  } else if (m > 0) { // here we need another strage boundary correction
134  for (l = 0; l <= k; l++) {
135  for (j = 0; j <= k; j++) {
136  weight = I.rcorr(m, l, j);
137  element_incr<T, SIZE1>(size1, ctemp1, weight, A + l * sa1,
138  B + j * sb1);
139  }
140  }
141  }
142  // CONTRIBUTION FROM TAU...BETA: mind the minus sign below
143  for (l = 0; l < sc1; l++) {
144  ctemp2[l] = 0;
145  }
146  if (ntau - m >= k2 - 1) { // usual gregory integration
147  amat = A + sa1 * ntau;
148  bmat = B + m * sb1;
149  for (j = m; j <= m + k; j++) {
150  weight = I.gregory_omega(j - m);
151  element_incr<T, SIZE1>(size1, ctemp2, weight, amat, bmat);
152  amat -= sa1;
153  bmat += sb1;
154  }
155  for (j = m + k1; j < ntau - k; j++) {
156  element_incr<T, SIZE1>(size1, ctemp2, amat, bmat);
157  amat -= sa1;
158  bmat += sb1;
159  }
160  for (j = ntau - k; j <= ntau; j++) {
161  weight = I.gregory_omega(ntau - j);
162  element_incr<T, SIZE1>(size1, ctemp2, weight, amat, bmat);
163  amat -= sa1;
164  bmat += sb1;
165  }
166  } else if (ntau - m >= k) {
167  amat = A + sa1 * ntau;
168  bmat = B + sb1 * m;
169  for (j = m; j <= ntau; j++) {
170  weight = I.gregory_weights(ntau - m, ntau - j);
171  element_incr<T, SIZE1>(size1, ctemp2, weight, amat, bmat);
172  amat -= sa1;
173  bmat += sb1;
174  }
175  } else if (ntau - m >
176  0) { // here we need another strange boundary correction
177  // std::cout << "i am here " << std::endl;
178  for (l = 0; l <= k; l++) {
179  for (j = 0; j <= k; j++) {
180  weight = I.rcorr(ntau - m, l, j);
181  element_incr<T, SIZE1>(size1, ctemp2, weight,
182  A + sa1 * (ntau - l),
183  B + sb1 * (ntau - j));
184  }
185  }
186  }
187  // cmat += ctemp1 + sig * ctemp2:
188  for (l = 0; l < sc1; l++)
189  C[l] = ctemp1[l] + std::complex<T>(sig, 0.0) * ctemp2[l];
190  delete[] ctemp1;
191  delete[] ctemp2;
192  return;
193 }
195 
225 template <typename T, int SIZE1>
226 void matsubara_integral_1_1(int size1, int m, int ntau, std::complex<T> *C,
227  std::complex<T> *A, std::complex<T> *B,
228  integration::Integrator<T> &I) {
229  typedef std::complex<T> cplx;
230  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1, j, l, sa1, sb1, sc1;
231  cplx *amat, *bmat, *ctemp1;
232  T weight;
233 
234  sa1 = size1 * size1;
235  sb1 = sa1;
236  sc1 = sa1;
237  // std::cout << "conv " << m << " " << ntau << std::endl;
238  ctemp1 = new cplx[sc1];
239  // CONTRIBUTION FROM 0...TAU
240  for (l = 0; l < sc1; l++) {
241  ctemp1[l] = 0;
242  }
243  if (m >= k2 - 1) {
244  amat = A + m * sa1;
245  bmat = B;
246  for (j = 0; j <= k; j++) {
247  weight = I.gregory_omega(j);
248  element_incr<T, SIZE1>(size1, ctemp1, weight, amat, bmat);
249  amat -= sa1;
250  bmat += sb1;
251  }
252  for (j = k1; j < m - k; j++) {
253  element_incr<T, SIZE1>(size1, ctemp1, amat, bmat);
254  amat -= sa1;
255  bmat += sb1;
256  }
257  for (j = m - k; j <= m; j++) {
258  weight = I.gregory_omega(m - j);
259  element_incr<T, SIZE1>(size1, ctemp1, weight, amat, bmat);
260  amat -= sa1;
261  bmat += sb1;
262  }
263  } else if (m >= k) {
264  amat = A + m * sa1;
265  bmat = B;
266  for (j = 0; j <= m; j++) {
267  weight = I.gregory_weights(m, j);
268  element_incr<T, SIZE1>(size1, ctemp1, weight, amat, bmat);
269  amat -= sa1;
270  bmat += sb1;
271  }
272  } else if (m > 0) { // here we need another strage boundary correction
273  for (l = 0; l <= k; l++) {
274  for (j = 0; j <= k; j++) {
275  weight = I.rcorr(m, l, j);
276  element_incr<T, SIZE1>(size1, ctemp1, weight, A + l * sa1,
277  B + j * sb1);
278  }
279  }
280  }
281  // cmat += ctemp1:
282  for (l = 0; l < sc1; l++)
283  C[l] = ctemp1[l];
284  delete[] ctemp1;
285  return;
286 }
288 
318 template <typename T, int SIZE1>
319 void matsubara_integral_2_2(int size1, int m, int ntau, std::complex<T> *C,
320  std::complex<T> *A, std::complex<T> *B,
321  integration::Integrator<T> &I) {
322  typedef std::complex<T> cplx;
323  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1, j, l, sa1, sb1, sc1;
324  cplx *amat, *bmat, *ctemp2;
325  T weight;
326  sa1 = size1 * size1;
327  sb1 = sa1;
328  sc1 = sa1;
329  ctemp2 = new cplx[sc1];
330  // CONTRIBUTION FROM TAU ... BETA
331  for (l = 0; l < sc1; l++)
332  ctemp2[l] = 0;
333  if (m == ntau) { // donothing
334  } else if (m > ntau - k) {
335  for (l = 0; l <= k; l++) {
336  for (j = 0; j <= k; j++) {
337  weight = I.rcorr(ntau - m, l, j);
338  element_incr<T, SIZE1>(size1, ctemp2, weight,
339  A + sa1 * (ntau - l), B + sb1 * j);
340  }
341  }
342  } else if (m > ntau - k2 + 1) {
343  amat = A + sa1 * m;
344  bmat = B;
345  for (l = 0; l <= ntau - m; l++) {
346  weight = I.gregory_weights(ntau - m, l);
347  element_incr<T, SIZE1>(size1, ctemp2, weight, amat, bmat);
348  amat += sa1;
349  bmat += sb1;
350  }
351  } else {
352  amat = A + sa1 * m;
353  bmat = B;
354  for (l = m; l <= m + k; l++) {
355  weight = I.gregory_omega(l - m);
356  element_incr<T, SIZE1>(size1, ctemp2, weight, amat, bmat);
357  amat += sa1;
358  bmat += sb1;
359  }
360  for (l = m + k1; l < ntau - k; l++) {
361  element_incr<T, SIZE1>(size1, ctemp2, amat, bmat);
362  amat += sa1;
363  bmat += sb1;
364  }
365  for (l = ntau - k; l <= ntau; l++) {
366  weight = I.gregory_omega(ntau - l);
367  element_incr<T, SIZE1>(size1, ctemp2, weight, amat, bmat);
368  amat += sa1;
369  bmat += sb1;
370  }
371  }
372  // ctv += ctemp2 :
373  for (l = 0; l < sc1; l++)
374  C[l] = ctemp2[l];
375  delete[] ctemp2;
376 }
378 
409 template <typename T, int SIZE1>
410 void matsubara_integral_2(int size1, int m, int ntau, std::complex<T> *C,
411  std::complex<T> *A, std::complex<T> *B,
412  integration::Integrator<T> &I, int sig) {
413  typedef std::complex<T> cplx;
414  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1, j, l, sa1, sb1, sc1;
415  cplx *amat, *bmat, *ctemp1, *ctemp2;
416  T weight;
417 
418  sa1 = size1 * size1;
419  sb1 = sa1;
420  sc1 = sa1;
421 
422  ctemp1 = new cplx[sc1];
423  ctemp2 = new cplx[sc1];
424  // CONTRIBUTION FROM 0 ... TAU
425  for (l = 0; l < sc1; l++)
426  ctemp1[l] = 0;
427  if (m == 0) {
428  // donothing
429  } else if (m < k) {
430  for (j = 0; j <= k; j++) {
431  for (l = 0; l <= k; l++) {
432  weight = I.rcorr(m, l, j);
433  element_incr<T, SIZE1>(size1, ctemp1, weight, A + j * sa1,
434  B + sb1 * (ntau - l));
435  }
436  }
437  } else if (m < k2 - 1) {
438  amat = A + sa1 * m;
439  bmat = B + sb1 * ntau;
440  for (l = 0; l <= m; l++) {
441  weight = I.gregory_weights(m, l);
442  element_incr<T, SIZE1>(size1, ctemp1, weight, amat, bmat);
443  amat -= sa1;
444  bmat -= sb1;
445  }
446  } else {
447  amat = A + sa1 * m;
448  bmat = B + sb1 * ntau;
449  for (l = 0; l <= k; l++) {
450  weight = I.gregory_omega(l);
451  element_incr<T, SIZE1>(size1, ctemp1, weight, amat, bmat);
452  amat -= sa1;
453  bmat -= sb1;
454  }
455  for (l = k1; l < m - k; l++) {
456  element_incr<T, SIZE1>(size1, ctemp1, amat, bmat);
457  amat -= sa1;
458  bmat -= sb1;
459  }
460  for (l = m - k; l <= m; l++) {
461  weight = I.gregory_omega(m - l);
462  element_incr<T, SIZE1>(size1, ctemp1, weight, amat, bmat);
463  amat -= sa1;
464  bmat -= sb1;
465  }
466  }
467  // CONTRIBUTION FROM TAU ... BETA
468  for (l = 0; l < sc1; l++)
469  ctemp2[l] = 0;
470  if (m == ntau) { // donothing
471  } else if (m > ntau - k) {
472  for (l = 0; l <= k; l++) {
473  for (j = 0; j <= k; j++) {
474  weight = I.rcorr(ntau - m, l, j);
475  element_incr<T, SIZE1>(size1, ctemp2, weight,
476  A + sa1 * (ntau - l), B + sb1 * j);
477  }
478  }
479  } else if (m > ntau - k2 + 1) {
480  amat = A + sa1 * m;
481  bmat = B;
482  for (l = 0; l <= ntau - m; l++) {
483  weight = I.gregory_weights(ntau - m, l);
484  element_incr<T, SIZE1>(size1, ctemp2, weight, amat, bmat);
485  amat += sa1;
486  bmat += sb1;
487  }
488  } else {
489  amat = A + sa1 * m;
490  bmat = B;
491  for (l = m; l <= m + k; l++) {
492  weight = I.gregory_omega(l - m);
493  element_incr<T, SIZE1>(size1, ctemp2, weight, amat, bmat);
494  amat += sa1;
495  bmat += sb1;
496  }
497  for (l = m + k1; l < ntau - k; l++) {
498  element_incr<T, SIZE1>(size1, ctemp2, amat, bmat);
499  amat += sa1;
500  bmat += sb1;
501  }
502  for (l = ntau - k; l <= ntau; l++) {
503  weight = I.gregory_omega(ntau - l);
504  element_incr<T, SIZE1>(size1, ctemp2, weight, amat, bmat);
505  amat += sa1;
506  bmat += sb1;
507  }
508  }
509  // ctv += ctemp2 + sig * ctemp1:
510  for (l = 0; l < sc1; l++)
511  C[l] = ctemp2[l] + std::complex<T>(sig, 0.0) * ctemp1[l];
512 
513  delete[] ctemp1;
514  delete[] ctemp2;
515 }
516 
518 
544 template <typename T, class GG, int SIZE1>
545 void convolution_matsubara_dispatch(GG &C, GG &A, GG &B,
546  integration::Integrator<T> &I, T beta) {
547  int ntau, l, m, size1 = C.size1();
548  std::complex<T> *cmat;
549  T dtau;
550  ntau = A.ntau();
551  for (m = 0; m <= ntau;
552  m++) { // compute cmat(m*dtau) = int_0^beta dx amat(tau-x) b(x)
553  matsubara_integral_1<T, SIZE1>(size1, m, ntau, C.matptr(m),
554  A.matptr(0), B.matptr(0), I, A.sig());
555  }
556  // multiply by dtau:
557  dtau = beta / ntau;
558  cmat = C.matptr(0);
559  m = (ntau + 1) * C.element_size();
560  for (l = 0; l < m; l++)
561  cmat[l] *= dtau;
562  return;
563 }
565 
593 template <typename T, class GG>
594 void convolution_matsubara(GG &C, GG &A, GG &B, integration::Integrator<T> &I,
595  T beta) {
596  int size1 = A.size1();
597  assert(B.ntau() == A.ntau());
598  assert(C.ntau() == A.ntau());
599  assert(B.size1() == size1);
600  assert(C.size1() == size1);
601  if (size1 == 1)
602  convolution_matsubara_dispatch<T, GG, 1>(C, A, B, I, beta);
603  else
604  convolution_matsubara_dispatch<T, GG, LARGESIZE>(C, A, B, I, beta);
605 }
606 
607 #if CNTR_USE_OMP == 1
608 
610 
638 template <typename T, class GG, int SIZE1>
639 void convolution_matsubara_omp_dispatch(int nomp, GG &C, GG &A, GG &B,
640  integration::Integrator<T> &I, T beta) {
641  int ntau, l, m, size1 = C.size1();
642  std::complex<T> *cmat;
643  T dtau;
644  ntau = A.ntau();
645 #pragma omp parallel for num_threads(nomp)
646  for (m = 0; m <= ntau;
647  m++) { // compute cmat(m*dtau) = int_0^beta dx amat(tau-x) b(x)
648  matsubara_integral_1<T, SIZE1>(size1, m, ntau, C.matptr(m),
649  A.matptr(0), B.matptr(0), I, A.sig());
650  }
651  // multiply by dtau:
652  dtau = beta / ntau;
653  cmat = C.matptr(0);
654  m = (ntau + 1) * C.element_size();
655  for (l = 0; l < m; l++)
656  cmat[l] *= dtau;
657  return;
658 }
660 
688 template <typename T, class GG>
689 void convolution_matsubara_nomp(int nomp, GG &C, GG &A, GG &B, integration::Integrator<T> &I,
690  T beta) {
691  int size1 = A.size1();
692  assert(B.ntau() == A.ntau());
693  assert(C.ntau() == A.ntau());
694  assert(B.size1() == size1);
695  assert(C.size1() == size1);
696  if (size1 == 1)
697  convolution_matsubara_omp_dispatch<T, GG, 1>(nomp, C, A, B, I, beta);
698  else
699  convolution_matsubara_omp_dispatch<T, GG, LARGESIZE>(nomp, C, A, B, I, beta);
700 }
701 
702 #endif // CNTR_USE_OMP
703 
705 
740 template <typename T, class GG, int SIZE1>
741 void convolution_timestep_ret(int n, GG &C, GG &A, GG &Acc, GG &B, GG &Bcc,
742  integration::Integrator<T> &I, T h) {
743  typedef std::complex<T> cplx;
744  int k = I.get_k();
745  int sa, sb, sc, j, m, j1, n1, l, size1 = C.size1();
746  cplx *aret, *cret, *bret, *btemp, *atemp, *result;
747  T weight;
748 
749  // duplicated arguments
750  sa = A.element_size();
751  sb = B.element_size();
752  sc = C.element_size();
753  n1 = (n < k ? k : n);
754  atemp = new cplx[sa];
755  btemp = new cplx[sb];
756  // such that G=Sigma*G can be called without creating a mess,
757  // data are first written in a temporary variable and then written to C at
758  // the end
759  result = new cplx[(n + 1) * sc];
760  for (l = 0; l < (n + 1) * sc; l++)
761  result[l] = 0;
762  // check consistency:
763  assert(sa * sb * sc != 0);
764  assert(Acc.element_size() == sa);
765  assert(Bcc.element_size() == sb);
766  assert(A.nt() >= n1);
767  assert(Acc.nt() >= n1);
768  assert(B.nt() >= n1);
769  assert(Bcc.nt() >= n1);
770  assert(C.nt() >= n);
771 
772  if (n >= k) {
773  // CONTRIBUTION FROM BRET: loop over lines of Bret
774  for (m = 0; m <= n; m++) { // contribution to integral from Bret(m,j)
775  aret = A.retptr(n, m);
776  for (l = 0; l < sa; l++)
777  atemp[l] = aret[l] * h; // here enters h
778  // the triangle j <= m
779  bret = B.retptr(m, 0);
780  cret = result;
781  // in the following sector the weights are 1
782  if (m < n - k) {
783  for (j = 0; j < m - k; j++) {
784  element_incr<T, SIZE1>(size1, cret, atemp,
785  bret); // cret += aret*bret
786  bret += sb;
787  cret += sc;
788  }
789  } else { // m>=n-k
790  weight = I.gregory_omega(n - m);
791  for (j = 0; j < m - k; j++) {
792  element_incr<T, SIZE1>(size1, cret, weight, atemp,
793  bret); // cret += aret*bret
794  bret += sb;
795  cret += sc;
796  }
797  }
798  // contribution from the stripe m-j <= k, with different weight
799  j1 = m - k;
800  if (j1 < 0)
801  j1 = 0;
802  for (j = j1; j <= m; j++) {
803  weight = I.gregory_weights(n - j, n - m);
804  element_incr<T, SIZE1>(size1, cret, weight, atemp,
805  bret); // cret += aret*bret
806  bret += sb;
807  cret += sc;
808  }
809  }
810  // CONTRIBUTION FROM BRET^CONJ:
811  for (m = n - k; m < n; m++) {
812  aret = A.retptr(n, m);
813  for (l = 0; l < sa; l++)
814  atemp[l] = aret[l] * h; // here enters h
815  for (j = m + 1; j <= n; j++) {
816  weight = I.gregory_weights(n - j, n - m);
817  element_conj<T, SIZE1>(size1, btemp, Bcc.retptr(j, m));
818  // mind the minus sign: B(m,j) continued to -Bcc(j,m)*
819  element_incr<T, SIZE1>(size1, result + j * sc, -weight, atemp,
820  btemp);
821  }
822  }
823  } else { // n < k
824  for (j = 0; j <= n; j++) {
825  cret = result + j * sc;
826  for (m = 0; m <= k; m++) {
827  weight = I.poly_integration(j, n, m) * h;
828  if (m >= j) {
829  for (l = 0; l < sb; l++)
830  btemp[l] = B.retptr(m, j)[l];
831  } else {
832  element_conj<T, SIZE1>(size1, btemp, Bcc.retptr(j, m));
833  weight *= -1;
834  }
835  if (n >= m) {
836  for (l = 0; l < sa; l++)
837  atemp[l] = A.retptr(n, m)[l];
838  } else {
839  element_conj<T, SIZE1>(size1, atemp, Acc.retptr(m, n));
840  weight *= -1;
841  }
842  // std::cout << n << " " << m << " " << j << std::endl;
843  element_incr<T, SIZE1>(size1, cret, weight, atemp, btemp);
844  }
845  }
846  }
847 
848  cret = C.retptr(n, 0);
849  for (l = 0; l < (n + 1) * sc; l++)
850  cret[l] = result[l];
851  delete[] result;
852  delete[] atemp;
853  delete[] btemp;
854  return;
855 }
857 
900 template <typename T, class GG, int SIZE1>
901 void convolution_timestep_tv(int n, std::complex<T> *ctv, GG &C, GG &A,
902  GG &Acc, GG &B, GG &Bcc,
903  integration::Integrator<T> &I, T beta, T h) {
904  typedef std::complex<T> cplx;
905 
906  int k = I.get_k(); // order
907  int sa = A.element_size();
908  int sb = B.element_size();
909  int sc = C.element_size();
910  int ntau = A.ntau();
911  int n1 = (n > k ? n : k);
912  int size1 = C.size1();
913  T dtau = beta / ntau;
914  T weight;
915  cplx *ctemp1, *ctv1, *btv, *atemp;
916  int j, m, l;
917 
918  // check consistency:
919  assert(sa * sb * sc != 0);
920  assert(Acc.element_size() == sa);
921  assert(Bcc.element_size() == sb);
922  assert(Acc.ntau() == ntau);
923  assert(B.ntau() == ntau);
924  assert(Bcc.ntau() == ntau);
925  assert(C.ntau() == ntau);
926  assert(A.nt() >= n1);
927  assert(Acc.nt() >= n1);
928  assert(B.nt() >= n1);
929  assert(Bcc.nt() >= n1);
930 
931  // CONTRIBUTION FROM Atv * Bmat:
932  // very similar to computing the matsubara convolution
933  ctemp1 = new cplx[sc];
934  for (m = 0; m <= ntau; m++) {
935  matsubara_integral_2<T, SIZE1>(size1, m, ntau, ctemp1, A.tvptr(n, 0),
936  B.matptr(0), I, B.sig());
937  for (l = 0; l < sc; l++)
938  ctv[m * sc + l] = dtau * ctemp1[l];
939  }
940  delete[] ctemp1;
941 
942  // CONTRIBUTION FROM Aret * Btv:
943  // loop over lines j
944  atemp = new cplx[sa];
945  for (j = 0; j <= n1; j++) { // j <= n1, n1 = max(n, k)
946  weight = I.gregory_weights(n, j);
947  if (n < j) { // j > n ==> n < k ; j <= max(n, k)
948  element_conj<T, SIZE1>(size1, atemp, Acc.retptr(j, n));
949  element_smul<T, SIZE1>(size1, atemp, -1);
950  // atemp = dt (-Acc(j,n)*)
951  } else { // j <= n
952  element_set<T, SIZE1>(size1, atemp, A.retptr(n, j));
953  // atemp = dt A(n,j)
954  }
955  element_smul<T, SIZE1>(size1, atemp, h); // here enters h
956  btv = B.tvptr(j, 0);
957  ctv1 = ctv;
958  // ctv1 = ctv1 + weight atemp . btv
959  if (weight != 1) {
960  for (m = 0; m <= ntau; m++) {
961  element_incr<T, SIZE1>(size1, ctv1, weight, atemp, btv);
962  btv += sb;
963  ctv1 += sc;
964  }
965  } else {
966  for (m = 0; m <= ntau; m++) {
967  element_incr<T, SIZE1>(size1, ctv1, atemp, btv);
968  btv += sb;
969  ctv1 += sc;
970  }
971  }
972  }
973  delete[] atemp;
974 }
976 
1015 template <typename T, class GG, int SIZE1>
1016 void convolution_timestep_tv(int n, GG &C, GG &A, GG &Acc, GG &B, GG &Bcc,
1017  integration::Integrator<T> &I, T beta, T h) {
1018  typedef std::complex<T> cplx;
1019  int ntau, m, sc, size1 = C.size1();
1020  cplx *ctv;
1021  ntau = C.ntau();
1022  sc = C.element_size();
1023 
1024  assert(sc > 0 && ntau > 0);
1025  ctv = new cplx[(ntau + 1) * sc];
1026  convolution_timestep_tv<T, GG, SIZE1>(n, ctv, C, A, Acc, B, Bcc, I, beta,
1027  h);
1028  for (m = 0; m <= ntau; m++)
1029  element_set<T, SIZE1>(size1, C.tvptr(n, m), ctv + m * sc);
1030  delete[] ctv;
1031 }
1033 
1074 template <typename T, class GG, int SIZE1>
1075 void
1076 convolution_timestep_les_tvvt(int n, int j1, int j2, std::complex<T> *cles,
1077  GG &C, GG &A, GG &Acc, GG &B, GG &Bcc,
1078  integration::Integrator<T> &I, T beta, T h) {
1079  typedef std::complex<T> cplx;
1080  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1;
1081  int sa, sb, sc, ntau, j, m, l, n1, sig, size1 = C.size1();
1082  T weight, dtau;
1083  cplx *atv, *btv, *btemp, *cles1, idtau;
1084 
1085  ntau = A.ntau();
1086  sa = A.element_size();
1087  sb = B.element_size();
1088  sc = C.element_size();
1089  dtau = beta / ntau;
1090  sig = A.sig();
1091  n1 = (n < k ? k : n);
1092  // check consistency:
1093  assert(sa * sb * sc != 0);
1094  assert(ntau >= k);
1095  assert(Acc.element_size() == sa);
1096  assert(Bcc.element_size() == sb);
1097  assert(Acc.ntau() == ntau);
1098  assert(B.ntau() == ntau);
1099  assert(Bcc.ntau() == ntau);
1100  assert(C.ntau() == ntau);
1101  assert(A.nt() >= n1);
1102  assert(Acc.nt() >= n1);
1103  assert(B.nt() >= n1);
1104  assert(Bcc.nt() >= n1);
1105  assert(sig == Acc.sig());
1106  assert(sig == B.sig());
1107  assert(sig == Bcc.sig());
1108  assert(0 <= j1 && j1 <= n1);
1109  assert(j1 <= j2 && j2 <= n1);
1110 
1111  // contribution from Atv*Bvt = Atv(jh,tau) * Bcc^tv(nh,beta-tau)^* *
1112  // (-Bose/Fermi)
1113  btemp = new cplx[(ntau + 1) * sb];
1114  idtau = cplx(0, -dtau);
1115  for (m = 0; m <= ntau; m++)
1116  element_conj<T, SIZE1>(size1, btemp + m * sb, Bcc.tvptr(n, ntau - m));
1117  for (l = 0; l < (ntau + 1) * sb; l++)
1118  btemp[l] *= idtau * (-(T)sig);
1119  for (j = j1; j <= j2; j++) {
1120  btv = btemp;
1121  atv = A.tvptr(j, 0);
1122  cles1 = cles + j * sc;
1123  if (ntau < k2 - 1) {
1124  for (m = 0; m <= ntau; m++) {
1125  weight = I.gregory_weights(ntau, m);
1126  element_incr<T, SIZE1>(size1, cles1, weight, atv, btv);
1127  atv += sa;
1128  btv += sb;
1129  }
1130  } else {
1131  for (m = 0; m <= k; m++) {
1132  weight = I.gregory_omega(m);
1133  element_incr<T, SIZE1>(size1, cles1, weight, atv, btv);
1134  atv += sa;
1135  btv += sb;
1136  }
1137  for (m = k1; m < ntau - k; m++) {
1138  element_incr<T, SIZE1>(size1, cles1, atv, btv);
1139  atv += sa;
1140  btv += sb;
1141  }
1142  for (m = ntau - k; m <= ntau; m++) {
1143  weight = I.gregory_omega(ntau - m);
1144  element_incr<T, SIZE1>(size1, cles1, weight, atv, btv);
1145  atv += sa;
1146  btv += sb;
1147  }
1148  }
1149  }
1150  delete[] btemp;
1151  return;
1152 }
1154 
1191 template <typename T, class GG, int SIZE1>
1192 void convolution_timestep_les_tvvt(int n, std::complex<T> *cles, GG &C, GG &A,
1193  GG &Acc, GG &B, GG &Bcc,
1194  integration::Integrator<T> &I, T beta,
1195  T h) {
1196  int k = I.get_k();
1197  int n1 = (n < k ? k : n);
1198  return convolution_timestep_les_tvvt<T, GG, SIZE1>(
1199  n, 0, n1, cles, C, A, Acc, B, Bcc, I, beta, h);
1200 }
1202 
1243 template <typename T, class GG, int SIZE1>
1244 void
1245 convolution_timestep_les_lesadv(int n, int j1, int j2, std::complex<T> *cles,
1246  GG &C, GG &A, GG &Acc, GG &B, GG &Bcc,
1247  integration::Integrator<T> &I, T beta, T h) {
1248  typedef std::complex<T> cplx;
1249  int k = I.get_k();
1250  int sa, sb, sc, ntau, j, m, l, n1, size1 = C.size1();
1251  T weight, dtau;
1252  cplx idtau, *ales, *badv;
1253 
1254  ntau = A.ntau();
1255  sa = A.element_size();
1256  sb = B.element_size();
1257  sc = C.element_size();
1258  dtau = beta / ntau;
1259  n1 = (n < k ? k : n); // j1 = 0, j2 = n1.
1260  // check consistency:
1261  assert(sa * sb * sc != 0);
1262  assert(ntau >= k);
1263  assert(Acc.element_size() == sa);
1264  assert(Bcc.element_size() == sb);
1265  assert(Acc.ntau() == ntau);
1266  assert(B.ntau() == ntau);
1267  assert(Bcc.ntau() == ntau);
1268  assert(C.ntau() == ntau);
1269  assert(A.nt() >= n1);
1270  assert(Acc.nt() >= n1);
1271  assert(B.nt() >= n1);
1272  assert(Bcc.nt() >= n1);
1273  assert(0 <= j1 && j1 <= n1);
1274  assert(j1 <= j2 && j2 <= n1);
1275 
1276  // contribution from Ales(j,m)*Badv(m,n)
1277  ales = new cplx[sa];
1278  badv = new cplx[(n1 + 1) * sb];
1279  for (m = 0; m <= n1; m++) {
1280  weight = h * I.gregory_weights(n, m);
1281  if (m <= n) {
1282  element_conj<T, SIZE1>(size1, badv + m * sb, Bcc.retptr(n, m));
1283  for (l = 0; l < sb; l++)
1284  badv[m * sb + l] *= weight;
1285  } else {
1286  for (l = 0; l < sb; l++)
1287  badv[m * sb + l] = -weight * B.retptr(m, n)[l];
1288  }
1289  }
1290  // Cles(t',t) += \int_0^{t'} ds -Ales*(s,t') * Bret*(t,s)
1291  for (j = j1; j <= j2; j++) {
1292  for (m = 0; m < j; m++) { // inner loop over cache-optimal `s`.
1293  element_minusconj<T, SIZE1>(size1, ales, Acc.lesptr(m, j));
1294  element_incr<T, SIZE1>(size1, cles + j * sc, ales, badv + m * sb);
1295  }
1296  }
1297  // Cles(t',t) += \int_{t'}^t ds Ales(t',s) * Bret*(t,s)
1298  for (m = 0; m <= n1; ++m) {
1299  int jmax = std::min(j2, m);
1300  for (j = j1; j <= jmax; ++j) { // inner loop over cache-optimal `t'`.
1301  element_incr<T, SIZE1>(size1, cles + j * sc, A.lesptr(j, m),
1302  badv + m * sb);
1303  }
1304  }
1305  delete[] ales;
1306  delete[] badv;
1307  return;
1308 }
1310 
1347 template <typename T, class GG, int SIZE1>
1348 void convolution_timestep_les_lesadv(int n, std::complex<T> *cles, GG &C,
1349  GG &A, GG &Acc, GG &B, GG &Bcc,
1350  integration::Integrator<T> &I, T beta,
1351  T h) {
1352  int k = I.get_k();
1353  int n1 = (n < k ? k : n);
1354  return convolution_timestep_les_lesadv<T, GG, SIZE1>(
1355  n, 0, n1, cles, C, A, Acc, B, Bcc, I, beta, h);
1356 }
1358 
1395 template <typename T, class GG, int SIZE1>
1396 void convolution_timestep_les_retles(int n, std::complex<T> *cles, GG &C,
1397  GG &A, GG &Acc, GG &B, GG &Bcc,
1398  integration::Integrator<T> &I, T beta,
1399  T h) {
1400  typedef std::complex<T> cplx;
1401  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1, size1 = C.size1();
1402  int sa, sb, sc, ntau, j, m, l, n1;
1403  T weight, dtau;
1404  cplx *aret, *atemp, *btemp, *cles1, *bles, idtau;
1405 
1406  ntau = A.ntau();
1407  sa = A.element_size();
1408  sb = B.element_size();
1409  sc = C.element_size();
1410  dtau = beta / ntau;
1411  n1 = (n < k ? k : n);
1412  // check consistency:
1413  assert(sa * sb * sc != 0);
1414  assert(ntau >= k);
1415  assert(Acc.element_size() == sa);
1416  assert(Bcc.element_size() == sb);
1417  assert(Acc.ntau() == ntau);
1418  assert(B.ntau() == ntau);
1419  assert(Bcc.ntau() == ntau);
1420  assert(C.ntau() == ntau);
1421  assert(A.nt() >= n1);
1422  assert(Acc.nt() >= n1);
1423  assert(B.nt() >= n1);
1424  assert(Bcc.nt() >= n1);
1425 
1426  // contribution from Aret*Bles
1427  btemp = new cplx[(n1 + 1) * sb];
1428  atemp = new cplx[sa];
1429  for (m = 0; m <= n1; m++) { // btemp(m) --> B^<(m,n)
1430  if (m <= n) {
1431  for (l = 0; l < sb; l++)
1432  btemp[m * sb + l] = h * B.lesptr(m, n)[l];
1433  } else {
1434  element_conj<T, SIZE1>(size1, btemp + m * sb, Bcc.lesptr(n, m));
1435  for (l = 0; l < sb; l++)
1436  btemp[m * sb + l] *= -h;
1437  }
1438  }
1439  for (j = 0; j <= n; j++) { // compute -> cles1(j,n)
1440  cles1 = cles + j * sc;
1441  // CONTRINBUTION FROM A_RET
1442  if (j >= k2 - 1) {
1443  aret = A.retptr(j, 0);
1444  bles = btemp;
1445  for (m = 0; m <= k; m++) {
1446  weight = I.gregory_omega(m);
1447  element_incr<T, SIZE1>(size1, cles1, weight, aret, bles);
1448  bles += sb;
1449  aret += sa;
1450  }
1451  for (m = k1; m < j - k; m++) {
1452  element_incr<T, SIZE1>(size1, cles1, aret, bles);
1453  bles += sb;
1454  aret += sa;
1455  }
1456  for (m = j - k; m <= j; m++) {
1457  weight = I.gregory_omega(j - m);
1458  element_incr<T, SIZE1>(size1, cles1, weight, aret, bles);
1459  bles += sb;
1460  aret += sa;
1461  }
1462  } else {
1463  aret = A.retptr(j, 0);
1464  bles = btemp;
1465  for (m = 0; m <= j; m++) {
1466  weight = I.gregory_weights(j, m);
1467  element_incr<T, SIZE1>(size1, cles1, weight, aret, bles);
1468  bles += sb;
1469  aret += sa;
1470  }
1471  }
1472  // CONTRIBUTION FROM ACC_RET
1473  if (j < k) {
1474  for (m = j + 1; m <= k; m++) {
1475  element_conj<T, SIZE1>(size1, atemp, Acc.retptr(m, j));
1476  weight = -I.gregory_weights(j, m);
1477  bles = btemp + m * sb;
1478  element_incr<T, SIZE1>(size1, cles1, weight, atemp, bles);
1479  }
1480  }
1481  }
1482  delete[] btemp;
1483  delete[] atemp;
1484  return;
1485 }
1487 
1524 template <typename T, class GG, int SIZE1>
1525 void convolution_timestep_les(int n, GG &C, GG &A, GG &Acc, GG &B, GG &Bcc,
1526  integration::Integrator<T> &I, T beta, T h) {
1527  typedef std::complex<T> cplx;
1528  cplx *cles;
1529  int m, sc, n1, k = I.get_k(), size1 = C.size1();
1530  sc = C.element_size();
1531  assert(sc > 0);
1532  n1 = (k > n ? k : n);
1533  cles = new cplx[(n1 + 1) * sc];
1534  for (m = 0; m < sc * (n1 + 1); m++)
1535  cles[m] = 0;
1536  convolution_timestep_les_tvvt<T, GG, SIZE1>(n, cles, C, A, Acc, B, Bcc, I,
1537  beta, h);
1538  convolution_timestep_les_lesadv<T, GG, SIZE1>(n, cles, C, A, Acc, B, Bcc,
1539  I, beta, h);
1540  convolution_timestep_les_retles<T, GG, SIZE1>(n, cles, C, A, Acc, B, Bcc,
1541  I, beta, h);
1542  for (m = 0; m <= n; m++)
1543  element_set<T, SIZE1>(size1, C.lesptr(m, n), cles + m * sc);
1544  delete[] cles;
1545  return;
1546 }
1548 
1581 template <typename T>
1582 void convolution_timestep(int n, herm_matrix<T> &C, herm_matrix<T> &A,
1583  herm_matrix<T> &Acc, herm_matrix<T> &B,
1584  herm_matrix<T> &Bcc, integration::Integrator<T> &I,
1585  T beta, T h) {
1586  int size1 = C.size1(), ntau = C.ntau(), k = I.k(), n1 = (n < k ? k : n);
1587  if (n == -1) {
1588  convolution_matsubara(C, A, B, I, beta);
1589  return;
1590  }
1591  assert(n >= 0);
1592  assert(A.size1() == size1);
1593  assert(Acc.size1() == size1);
1594  assert(B.size1() == size1);
1595  assert(Bcc.size1() == size1);
1596  assert(A.ntau() == ntau);
1597  assert(Acc.ntau() == ntau);
1598  assert(B.ntau() == ntau);
1599  assert(Bcc.ntau() == ntau);
1600  assert(A.nt() >= n1);
1601  assert(Acc.nt() >= n1);
1602  assert(B.nt() >= n1);
1603  assert(Bcc.nt() >= n1);
1604  assert(C.nt() >= n);
1605  if (size1 == 1) {
1606  convolution_timestep_ret<T, herm_matrix<T>, 1>(n, C, A, Acc, B, Bcc,
1607  I, h);
1608  convolution_timestep_tv<T, herm_matrix<T>, 1>(n, C, A, Acc, B, Bcc, I,
1609  beta, h);
1610  convolution_timestep_les<T, herm_matrix<T>, 1>(n, C, A, Acc, B, Bcc,
1611  I, beta, h);
1612  } else {
1613  convolution_timestep_ret<T, herm_matrix<T>, LARGESIZE>(n, C, A, Acc,
1614  B, Bcc, I, h);
1615  convolution_timestep_tv<T, herm_matrix<T>, LARGESIZE>(
1616  n, C, A, Acc, B, Bcc, I, beta, h);
1617  convolution_timestep_les<T, herm_matrix<T>, LARGESIZE>(
1618  n, C, A, Acc, B, Bcc, I, beta, h);
1619  }
1620 }
1621 
1622 
1656 template <typename T>
1657 void convolution_timestep(int n, herm_matrix<T> &C, herm_matrix<T> &A,
1658  herm_matrix<T> &Acc, herm_matrix<T> &B,
1659  herm_matrix<T> &Bcc,
1660  T beta, T h, int SolveOrder) {
1661  int size1 = C.size1(), ntau = C.ntau(), n1 = (n < SolveOrder ? SolveOrder : n);
1662  if (n == -1) {
1663  convolution_matsubara(C, A, B, integration::I<T>(SolveOrder), beta);
1664  return;
1665  }
1666  assert(n >= 0);
1667  assert(A.size1() == size1);
1668  assert(Acc.size1() == size1);
1669  assert(B.size1() == size1);
1670  assert(Bcc.size1() == size1);
1671  assert(A.ntau() == ntau);
1672  assert(Acc.ntau() == ntau);
1673  assert(B.ntau() == ntau);
1674  assert(Bcc.ntau() == ntau);
1675  assert(A.nt() >= n1);
1676  assert(Acc.nt() >= n1);
1677  assert(B.nt() >= n1);
1678  assert(Bcc.nt() >= n1);
1679  assert(C.nt() >= n);
1680  if (size1 == 1) {
1681  convolution_timestep_ret<T, herm_matrix<T>, 1>(n, C, A, Acc, B, Bcc,
1682  integration::I<T>(SolveOrder), h);
1683  convolution_timestep_tv<T, herm_matrix<T>, 1>(n, C, A, Acc, B, Bcc, integration::I<T>(SolveOrder),
1684  beta, h);
1685  convolution_timestep_les<T, herm_matrix<T>, 1>(n, C, A, Acc, B, Bcc,
1686  integration::I<T>(SolveOrder), beta, h);
1687  } else {
1688  convolution_timestep_ret<T, herm_matrix<T>, LARGESIZE>(n, C, A, Acc,
1689  B, Bcc, integration::I<T>(SolveOrder), h);
1690  convolution_timestep_tv<T, herm_matrix<T>, LARGESIZE>(
1691  n, C, A, Acc, B, Bcc, integration::I<T>(SolveOrder), beta, h);
1692  convolution_timestep_les<T, herm_matrix<T>, LARGESIZE>(
1693  n, C, A, Acc, B, Bcc, integration::I<T>(SolveOrder), beta, h);
1694  }
1695 }
1696 
1698 
1728 template <typename T>
1729 void convolution_timestep(int n, herm_matrix<T> &C, herm_matrix<T> &A,
1730  herm_matrix<T> &B, integration::Integrator<T> &I,
1731  T beta, T h) {
1732  convolution_timestep<T>(n, C, A, A, B, B, I, beta, h);
1733 }
1734 
1735 
1766 template <typename T>
1767 void convolution_timestep(int n, herm_matrix<T> &C, herm_matrix<T> &A,
1768  herm_matrix<T> &B,
1769  T beta, T h, int SolveOrder) {
1770  convolution_timestep<T>(n, C, A, A, B, B, beta, h, SolveOrder);
1771 }
1773 
1804 template <typename T>
1805 void convolution(herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &Acc,
1806  herm_matrix<T> &B, herm_matrix<T> &Bcc,
1807  integration::Integrator<T> &I, T beta, T h) {
1808  int tstp;
1809  convolution_matsubara(C, A, B, I, beta);
1810  for (tstp = 0; tstp <= C.nt(); tstp++)
1811  convolution_timestep<T>(tstp, C, A, Acc, B, Bcc, I, beta, h);
1812 }
1813 
1815 
1844 template <typename T, class GG, int SIZE1>
1845 void convolution_matsubara_dispatch(GG &C, GG &A, std::complex<T> *f0, GG &B,
1846  integration::Integrator<T> &I, T beta) {
1847  int ntau, l, m, size1 = C.size1(), sb = B.element_size();
1848  std::complex<T> *cmat, *bmat;
1849  T dtau;
1850  ntau = A.ntau();
1851  bmat = new std::complex<T>[sb * (ntau + 1)];
1852  for (m = 0; m <= ntau; m++)
1853  element_mult<T, SIZE1>(size1, bmat + m * sb, f0, B.matptr(m));
1854  for (m = 0; m <= ntau; m++) { // compute cmat(m*dtau) = int_0^beta dx amat(tau-x) b(x)
1855  matsubara_integral_1<T, SIZE1>(size1, m, ntau, C.matptr(m), A.matptr(0), bmat, I,
1856  A.sig());
1857  }
1858  delete[] bmat;
1859  // multiply by dtau:
1860  dtau = beta / ntau;
1861  cmat = C.matptr(0);
1862  m = (ntau + 1) * C.element_size();
1863  for (l = 0; l < m; l++)
1864  cmat[l] *= dtau;
1865  return;
1866 }
1868 
1896 template <typename T, class GG>
1897 void convolution_matsubara(GG &C, GG &A, std::complex<T> *f0, GG &B,
1898  integration::Integrator<T> &I, T beta) {
1899  int size1 = A.size1();
1900  assert(B.ntau() == A.ntau());
1901  assert(C.ntau() == A.ntau());
1902  assert(B.size1() == size1);
1903  assert(C.size1() == size1);
1904  if (size1 == 1)
1905  convolution_matsubara_dispatch<T, GG, 1>(C, A, f0, B, I, beta);
1906  else
1907  convolution_matsubara_dispatch<T, GG, LARGESIZE>(C, A, f0, B, I, beta);
1908 }
1910 
1947 template <typename T, class GG, int SIZE1>
1948 void convolution_timestep_ret(int n, GG &C, GG &A, GG &Acc, std::complex<T> *ft, GG &B,
1949  GG &Bcc, integration::Integrator<T> &I, T h) {
1950  typedef std::complex<T> cplx;
1951  int k = I.get_k();
1952  int sf, sa, sb, sc, j, m, j1, n1, l, size1 = C.size1();
1953  cplx *aret, *cret, *bret, *btemp, *atemp, *result;
1954  T weight;
1955 
1956  // duplicated arguments
1957  sa = A.element_size();
1958  sb = B.element_size();
1959  sc = C.element_size();
1960  sf = size1 * size1;
1961  n1 = (n < k ? k : n);
1962  atemp = new cplx[sa];
1963  aret = new cplx[sa];
1964  btemp = new cplx[sb];
1965  // such that G=Sigma*G can be called without creating a mess,
1966  // data are first written in a temporary variable and then written to C at the end
1967  result = new cplx[(n + 1) * sc];
1968  for (l = 0; l < (n + 1) * sc; l++)
1969  result[l] = 0;
1970  // check consistency:
1971  assert(sa * sb * sc != 0);
1972  assert(Acc.element_size() == sa);
1973  assert(Bcc.element_size() == sb);
1974  assert(A.nt() >= n1);
1975  assert(Acc.nt() >= n1);
1976  assert(B.nt() >= n1);
1977  assert(Bcc.nt() >= n1);
1978  assert(C.nt() >= n);
1979 
1980  if (n >= k) {
1981  // CONTRIBUTION FROM BRET: loop over lines of Bret
1982  for (m = 0; m <= n; m++) { // contribution to integral from Bret(m,j)
1983  element_mult<T, SIZE1>(size1, aret, A.retptr(n, m), ft + m * sf);
1984  // aret = A.retptr(n,m); // without f
1985  for (l = 0; l < sa; l++)
1986  atemp[l] = aret[l] * h; // here enters h
1987  // the triangle j <= m
1988  bret = B.retptr(m, 0);
1989  cret = result;
1990  // in the following sector the weights are 1
1991  if (m < n - k) {
1992  for (j = 0; j < m - k; j++) {
1993  element_incr<T, SIZE1>(size1, cret, atemp, bret); // cret += aret*bret
1994  bret += sb;
1995  cret += sc;
1996  }
1997  } else { // m>=n-k
1998  weight = I.gregory_omega(n - m);
1999  for (j = 0; j < m - k; j++) {
2000  element_incr<T, SIZE1>(size1, cret, weight, atemp,
2001  bret); // cret += aret*bret
2002  bret += sb;
2003  cret += sc;
2004  }
2005  }
2006  // contribution from the stripe m-j <= k, with different weight
2007  j1 = m - k;
2008  if (j1 < 0)
2009  j1 = 0;
2010  for (j = j1; j <= m; j++) {
2011  weight = I.gregory_weights(n - j, n - m);
2012  element_incr<T, SIZE1>(size1, cret, weight, atemp,
2013  bret); // cret += aret*bret
2014  bret += sb;
2015  cret += sc;
2016  }
2017  }
2018  // CONTRIBUTION FROM BRET^CONJ:
2019  for (m = n - k; m < n; m++) {
2020  element_mult<T, SIZE1>(size1, aret, A.retptr(n, m), ft + m * sf);
2021  // without f ... aret = A.retptr(n,m);
2022  for (l = 0; l < sa; l++)
2023  atemp[l] = aret[l] * h; // here enters h
2024  for (j = m + 1; j <= n; j++) {
2025  weight = I.gregory_weights(n - j, n - m);
2026  element_conj<T, SIZE1>(size1, btemp, Bcc.retptr(j, m));
2027  // mind the minus sign: B(m,j) continued to -Bcc(j,m)*
2028  element_incr<T, SIZE1>(size1, result + j * sc, -weight, atemp, btemp);
2029  }
2030  }
2031  } else { // n < k
2032  for (j = 0; j <= n; j++) {
2033  cret = result + j * sc;
2034  for (m = 0; m <= k; m++) {
2035  weight = I.poly_integration(j, n, m) * h;
2036  if (m >= j) {
2037  for (l = 0; l < sb; l++)
2038  btemp[l] = B.retptr(m, j)[l];
2039  } else {
2040  element_conj<T, SIZE1>(size1, btemp, Bcc.retptr(j, m));
2041  weight *= -1;
2042  }
2043  if (n >= m) {
2044  for (l = 0; l < sa; l++)
2045  atemp[l] = A.retptr(n, m)[l];
2046  } else {
2047  element_conj<T, SIZE1>(size1, atemp, Acc.retptr(m, n));
2048  weight *= -1;
2049  }
2050  element_mult<T, SIZE1>(size1, aret, atemp, ft + m * sf);
2051  element_incr<T, SIZE1>(size1, cret, weight, aret, btemp);
2052  // without f element_incr<T,SIZE1>(size1,cret,weight,atemp,btemp);
2053  }
2054  }
2055  }
2056 
2057  cret = C.retptr(n, 0);
2058  for (l = 0; l < (n + 1) * sc; l++)
2059  cret[l] = result[l];
2060  delete[] result;
2061  delete[] atemp;
2062  delete[] aret;
2063  delete[] btemp;
2064  return;
2065 }
2067 
2115 template <typename T, class GG, int SIZE1>
2116 void convolution_timestep_tv(int n, std::complex<T> *ctv, GG &C, GG &A, GG &Acc,
2117  std::complex<T> *f0, std::complex<T> *ft, GG &B, GG &Bcc,
2118  integration::Integrator<T> &I, T beta, T h) {
2119  typedef std::complex<T> cplx;
2120  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1;
2121  int sf, sa, sb, sc, ntau, j, m, l, n1, size1 = C.size1();
2122  T weight, dtau;
2123  cplx *ctemp1, *ctv1, *btv, *atemp, *bmat, *atemp1;
2124 
2125  ntau = A.ntau();
2126  sa = A.element_size();
2127  sb = B.element_size();
2128  sc = C.element_size();
2129  sf = size1 * size1;
2130  dtau = beta / ntau;
2131  n1 = (n > k ? n : k);
2132 
2133  // check consistency:
2134  assert(sa * sb * sc != 0);
2135  assert(ntau >= k2);
2136  assert(Acc.element_size() == sa);
2137  assert(Bcc.element_size() == sb);
2138  assert(Acc.ntau() == ntau);
2139  assert(B.ntau() == ntau);
2140  assert(Bcc.ntau() == ntau);
2141  assert(C.ntau() == ntau);
2142  assert(A.nt() >= n1);
2143  assert(Acc.nt() >= n1);
2144  assert(B.nt() >= n1);
2145  assert(Bcc.nt() >= n1);
2146 
2147  // CONTRIBUTION FROM Atv * Bmat : very similar to computing the matsubara convolution
2148  ctemp1 = new cplx[sc];
2149  bmat = new cplx[(ntau + 1) * sb];
2150  for (m = 0; m <= ntau; m++)
2151  element_mult<T, SIZE1>(size1, bmat + m * sb, f0, B.matptr(m));
2152  for (m = 0; m <= ntau; m++) {
2153  matsubara_integral_2<T, SIZE1>(size1, m, ntau, ctemp1, A.tvptr(n, 0), bmat, I,
2154  B.sig());
2155  for (l = 0; l < sc; l++)
2156  ctv[m * sc + l] = dtau * ctemp1[l];
2157  }
2158  delete[] bmat;
2159  delete[] ctemp1;
2160  // CONTRIBUTION FROM Aret * Btv: loop over lines j
2161  atemp = new cplx[sa];
2162  atemp1 = new cplx[sa];
2163  n1 = (n > k ? n : k);
2164  for (j = 0; j <= n1; j++) {
2165  weight = I.gregory_weights(n, j);
2166  if (n < j) {
2167  element_conj<T, SIZE1>(size1, atemp, Acc.retptr(j, n));
2168  element_smul<T, SIZE1>(size1, atemp, -1);
2169  } else {
2170  element_set<T, SIZE1>(size1, atemp, A.retptr(n, j));
2171  }
2172  element_smul<T, SIZE1>(size1, atemp, h); // here enters h
2173  element_mult<T, SIZE1>(size1, atemp1, atemp, ft + sf * j);
2174  btv = B.tvptr(j, 0);
2175  ctv1 = ctv;
2176  if (weight != 1) {
2177  for (m = 0; m <= ntau; m++) {
2178  element_incr<T, SIZE1>(size1, ctv1, weight, atemp1, btv);
2179  btv += sb;
2180  ctv1 += sc;
2181  }
2182  } else {
2183  for (m = 0; m <= ntau; m++) {
2184  element_incr<T, SIZE1>(size1, ctv1, atemp1, btv);
2185  btv += sb;
2186  ctv1 += sc;
2187  }
2188  }
2189  }
2190  delete[] atemp;
2191  delete[] atemp1;
2192  return;
2193 }
2195 
2238 template <typename T, class GG, int SIZE1>
2239 void convolution_timestep_tv(int n, GG &C, GG &A, GG &Acc, std::complex<T> *f0,
2240  std::complex<T> *ft, GG &B, GG &Bcc,
2241  integration::Integrator<T> &I, T beta, T h) {
2242  typedef std::complex<T> cplx;
2243  int ntau, m, sc, size1 = C.size1();
2244  cplx *ctv;
2245  ntau = C.ntau();
2246  sc = C.element_size();
2247 
2248  assert(sc > 0 && ntau > 0);
2249  ctv = new cplx[(ntau + 1) * sc];
2250  convolution_timestep_tv<T, GG, SIZE1>(n, ctv, C, A, Acc, f0, ft, B, Bcc, I, beta, h);
2251  for (m = 0; m <= ntau; m++)
2252  element_set<T, SIZE1>(size1, C.tvptr(n, m), ctv + m * sc);
2253  delete[] ctv;
2254 }
2256 
2295 template <typename T, class GG, int SIZE1>
2296 void convolution_timestep_les_tvvt(int n, std::complex<T> *cles, GG &C, GG &A, GG &Acc,
2297  std::complex<T> *f0, GG &B, GG &Bcc,
2298  integration::Integrator<T> &I, T beta, T h) {
2299  typedef std::complex<T> cplx;
2300  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1;
2301  int sa, sb, sc, ntau, j, m, l, n1, sig, size1 = C.size1();
2302  T weight, dtau;
2303  cplx *atv, *btv, *btemp, *btemp1, *cles1, idtau;
2304 
2305  ntau = A.ntau();
2306  sa = A.element_size();
2307  sb = B.element_size();
2308  sc = C.element_size();
2309  dtau = beta / ntau;
2310  sig = A.sig();
2311  n1 = (n < k ? k : n);
2312  // check consistency:
2313  assert(sa * sb * sc != 0);
2314  assert(ntau >= k);
2315  assert(Acc.element_size() == sa);
2316  assert(Bcc.element_size() == sb);
2317  assert(Acc.ntau() == ntau);
2318  assert(B.ntau() == ntau);
2319  assert(Bcc.ntau() == ntau);
2320  assert(C.ntau() == ntau);
2321  assert(A.nt() >= n1);
2322  assert(Acc.nt() >= n1);
2323  assert(B.nt() >= n1);
2324  assert(Bcc.nt() >= n1);
2325  assert(sig == Acc.sig());
2326  assert(sig == B.sig());
2327  assert(sig == Bcc.sig());
2328 
2329  // contribution from Atv*Bvt = Atv(jh,tau) * Bcc^tv(nh,beta-tau)^* * (-Bose/Fermi)
2330  btemp = new cplx[(ntau + 1) * sb];
2331  btemp1 = new cplx[sb];
2332  idtau = cplx(0, -dtau);
2333  for (m = 0; m <= ntau; m++) {
2334  element_conj<T, SIZE1>(size1, btemp1, Bcc.tvptr(n, ntau - m));
2335  element_mult<T, SIZE1>(size1, btemp + m * sb, f0, btemp1);
2336  }
2337  delete[] btemp1;
2338  for (l = 0; l < (ntau + 1) * sb; l++)
2339  btemp[l] *= idtau * (-(T)sig);
2340  for (j = 0; j <= n1; j++) {
2341  btv = btemp;
2342  atv = A.tvptr(j, 0);
2343  cles1 = cles + j * sc;
2344  if (ntau < k2 - 1) {
2345  for (m = 0; m <= ntau; m++) {
2346  weight = I.gregory_weights(ntau, m);
2347  element_incr<T, SIZE1>(size1, cles1, weight, atv, btv);
2348  atv += sa;
2349  btv += sb;
2350  }
2351  } else {
2352  for (m = 0; m <= k; m++) {
2353  weight = I.gregory_omega(m);
2354  element_incr<T, SIZE1>(size1, cles1, weight, atv, btv);
2355  atv += sa;
2356  btv += sb;
2357  }
2358  for (m = k1; m < ntau - k; m++) {
2359  element_incr<T, SIZE1>(size1, cles1, atv, btv);
2360  atv += sa;
2361  btv += sb;
2362  }
2363  for (m = ntau - k; m <= ntau; m++) {
2364  weight = I.gregory_omega(ntau - m);
2365  element_incr<T, SIZE1>(size1, cles1, weight, atv, btv);
2366  atv += sa;
2367  btv += sb;
2368  }
2369  }
2370  }
2371  delete[] btemp;
2372  return;
2373 }
2375 
2414 template <typename T, class GG, int SIZE1>
2415 void convolution_timestep_les_lesadv(int n, std::complex<T> *cles, GG &C, GG &A, GG &Acc,
2416  std::complex<T> *ft, GG &B, GG &Bcc,
2417  integration::Integrator<T> &I, T beta, T h) {
2418  typedef std::complex<T> cplx;
2419  int k = I.get_k();
2420  int sf, sa, sb, sc, ntau, j, m, l, n1, size1 = C.size1();
2421  T weight, dtau;
2422  cplx *atemp, *btemp, idtau, *ales, *badv, *btemp1;
2423 
2424  ntau = A.ntau();
2425  sa = A.element_size();
2426  sb = B.element_size();
2427  sc = C.element_size();
2428  sf = size1 * size1;
2429  dtau = beta / ntau;
2430  n1 = (n < k ? k : n);
2431  // check consistency:
2432  assert(sa * sb * sc != 0);
2433  assert(ntau >= k);
2434  assert(Acc.element_size() == sa);
2435  assert(Bcc.element_size() == sb);
2436  assert(Acc.ntau() == ntau);
2437  assert(B.ntau() == ntau);
2438  assert(Bcc.ntau() == ntau);
2439  assert(C.ntau() == ntau);
2440  assert(A.nt() >= n1);
2441  assert(Acc.nt() >= n1);
2442  assert(B.nt() >= n1);
2443  assert(Bcc.nt() >= n1);
2444 
2445  // contribution from Ales(j,m)*Badv(m,n)
2446  atemp = new cplx[(n1 + 1) * sa];
2447  btemp = new cplx[(n1 + 1) * sb];
2448  btemp1 = new cplx[sb];
2449  for (m = 0; m <= n1; m++) {
2450  weight = h * I.gregory_weights(n, m);
2451  if (m <= n) {
2452  element_conj<T, SIZE1>(size1, btemp1, Bcc.retptr(n, m));
2453  for (l = 0; l < sb; l++)
2454  btemp1[l] *= weight;
2455  } else {
2456  for (l = 0; l < sb; l++)
2457  btemp1[l] = -weight * B.retptr(m, n)[l];
2458  }
2459  element_mult<T, SIZE1>(size1, btemp + m * sb, ft + m * sf, btemp1);
2460  }
2461  delete[] btemp1;
2462  for (j = 0; j <= n1; j++) {
2463  for (m = 0; m <= n1; m++) {
2464  if (m < j) {
2465  element_conj<T, SIZE1>(size1, atemp + m * sa, Acc.lesptr(m, j));
2466  for (l = 0; l < sa; l++)
2467  atemp[m * sa + l] *= -1;
2468  } else {
2469  for (l = 0; l < sa; l++)
2470  atemp[m * sa + l] = A.lesptr(j, m)[l];
2471  }
2472  }
2473  ales = atemp;
2474  badv = btemp;
2475  for (m = 0; m <= n1; m++) {
2476  element_incr<T, SIZE1>(size1, cles + j * sc, ales, badv);
2477  ales += sa;
2478  badv += sb;
2479  }
2480  }
2481  delete[] atemp;
2482  delete[] btemp;
2483  return;
2484 }
2486 
2525 template <typename T, class GG, int SIZE1>
2526 void convolution_timestep_les_retles(int n, std::complex<T> *cles, GG &C, GG &A, GG &Acc,
2527  std::complex<T> *ft, GG &B, GG &Bcc,
2528  integration::Integrator<T> &I, T beta, T h) {
2529  typedef std::complex<T> cplx;
2530  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1, size1 = C.size1();
2531  int sf, sa, sb, sc, ntau, j, m, l, n1;
2532  T weight, dtau;
2533  cplx *aret, *atemp, *btemp, *btemp1, *cles1, *bles, idtau;
2534 
2535  ntau = A.ntau();
2536  sa = A.element_size();
2537  sb = B.element_size();
2538  sc = C.element_size();
2539  sf = size1 * size1;
2540  dtau = beta / ntau;
2541  n1 = (n < k ? k : n);
2542  // check consistency:
2543  assert(sa * sb * sc != 0);
2544  assert(ntau >= k);
2545  assert(Acc.element_size() == sa);
2546  assert(Bcc.element_size() == sb);
2547  assert(Acc.ntau() == ntau);
2548  assert(B.ntau() == ntau);
2549  assert(Bcc.ntau() == ntau);
2550  assert(C.ntau() == ntau);
2551  assert(A.nt() >= n1);
2552  assert(Acc.nt() >= n1);
2553  assert(B.nt() >= n1);
2554  assert(Bcc.nt() >= n1);
2555 
2556  // contribution from Aret*Bles
2557  btemp = new cplx[(n1 + 1) * sb];
2558  btemp1 = new cplx[sb];
2559  atemp = new cplx[sa];
2560  for (m = 0; m <= n1; m++) { // btemp(m) --> B^<(m,n)
2561  if (m <= n) {
2562  for (l = 0; l < sb; l++)
2563  btemp1[l] = h * B.lesptr(m, n)[l];
2564  } else {
2565  element_conj<T, SIZE1>(size1, btemp1, Bcc.lesptr(n, m));
2566  for (l = 0; l < sb; l++)
2567  btemp1[l] *= -h;
2568  }
2569  element_mult<T, SIZE1>(size1, btemp + m * sb, ft + m * sf, btemp1);
2570  }
2571  delete[] btemp1;
2572  for (j = 0; j <= n; j++) { // compute -> cles1(j,n)
2573  cles1 = cles + j * sc;
2574  // CONTRINBUTION FROM A_RET
2575  if (j >= k2 - 1) {
2576  aret = A.retptr(j, 0);
2577  bles = btemp;
2578  for (m = 0; m <= k; m++) {
2579  weight = I.gregory_omega(m);
2580  element_incr<T, SIZE1>(size1, cles1, weight, aret, bles);
2581  bles += sb;
2582  aret += sa;
2583  }
2584  for (m = k1; m < j - k; m++) {
2585  element_incr<T, SIZE1>(size1, cles1, aret, bles);
2586  bles += sb;
2587  aret += sa;
2588  }
2589  for (m = j - k; m <= j; m++) {
2590  weight = I.gregory_omega(j - m);
2591  element_incr<T, SIZE1>(size1, cles1, weight, aret, bles);
2592  bles += sb;
2593  aret += sa;
2594  }
2595  } else {
2596  aret = A.retptr(j, 0);
2597  bles = btemp;
2598  for (m = 0; m <= j; m++) {
2599  weight = I.gregory_weights(j, m);
2600  element_incr<T, SIZE1>(size1, cles1, weight, aret, bles);
2601  bles += sb;
2602  aret += sa;
2603  }
2604  }
2605  // CONTRIBUTION FROM ACC_RET
2606  if (j < k) {
2607  for (m = j + 1; m <= k; m++) {
2608  element_conj<T, SIZE1>(size1, atemp, Acc.retptr(m, j));
2609  weight = -I.gregory_weights(j, m);
2610  bles = btemp + m * sb;
2611  element_incr<T, SIZE1>(size1, cles1, weight, atemp, bles);
2612  }
2613  }
2614  }
2615  delete[] btemp;
2616  delete[] atemp;
2617  return;
2618 }
2620 
2662 template <typename T, class GG, int SIZE1>
2663 void convolution_timestep_les(int n, GG &C, GG &A, GG &Acc, std::complex<T> *f0,
2664  std::complex<T> *ft, GG &B, GG &Bcc,
2665  integration::Integrator<T> &I, T beta, T h) {
2666  typedef std::complex<T> cplx;
2667  cplx *cles;
2668  int m, sc, n1, k = I.get_k(), size1 = C.size1();
2669  sc = C.element_size();
2670  assert(sc > 0);
2671  n1 = (k > n ? k : n);
2672  cles = new cplx[(n1 + 1) * sc];
2673  for (m = 0; m < sc * (n1 + 1); m++)
2674  cles[m] = 0;
2675  convolution_timestep_les_tvvt<T, GG, SIZE1>(n, cles, C, A, Acc, f0, B, Bcc, I, beta, h);
2676  convolution_timestep_les_lesadv<T, GG, SIZE1>(n, cles, C, A, Acc, ft, B, Bcc, I, beta,
2677  h);
2678  convolution_timestep_les_retles<T, GG, SIZE1>(n, cles, C, A, Acc, ft, B, Bcc, I, beta,
2679  h);
2680  for (m = 0; m <= n; m++)
2681  element_set<T, SIZE1>(size1, C.lesptr(m, n), cles + m * sc);
2682  delete[] cles;
2683  return;
2684 }
2686 
2723 template <typename T>
2724 void convolution_timestep(int n, herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &Acc,
2725  std::complex<T> *f0, std::complex<T> *ft, herm_matrix<T> &B,
2726  herm_matrix<T> &Bcc, integration::Integrator<T> &I, T beta, T h) {
2727  int size1 = C.size1(), ntau = C.ntau(), k = I.k(), n1 = (n < k ? k : n);
2728  if (n == -1) {
2729  convolution_matsubara(C, A, f0, B, I, beta);
2730  return;
2731  }
2732  assert(n >= 0);
2733  assert(A.size1() == size1);
2734  assert(Acc.size1() == size1);
2735  assert(B.size1() == size1);
2736  assert(Bcc.size1() == size1);
2737  assert(A.ntau() == ntau);
2738  assert(Acc.ntau() == ntau);
2739  assert(B.ntau() == ntau);
2740  assert(Bcc.ntau() == ntau);
2741  assert(A.nt() >= n1);
2742  assert(Acc.nt() >= n1);
2743  assert(B.nt() >= n1);
2744  assert(Bcc.nt() >= n1);
2745  assert(C.nt() >= n);
2746  if (size1 == 1) {
2747  convolution_timestep_ret<T, herm_matrix<T>, 1>(n, C, A, Acc, ft, B, Bcc, I, h);
2748  convolution_timestep_tv<T, herm_matrix<T>, 1>(n, C, A, Acc, f0, ft, B, Bcc, I, beta,
2749  h);
2750  convolution_timestep_les<T, herm_matrix<T>, 1>(n, C, A, Acc, f0, ft, B, Bcc, I, beta,
2751  h);
2752  } else {
2753  convolution_timestep_ret<T, herm_matrix<T>, LARGESIZE>(n, C, A, Acc, ft, B, Bcc, I,
2754  h);
2755  convolution_timestep_tv<T, herm_matrix<T>, LARGESIZE>(n, C, A, Acc, f0, ft, B, Bcc,
2756  I, beta, h);
2757  convolution_timestep_les<T, herm_matrix<T>, LARGESIZE>(n, C, A, Acc, f0, ft, B, Bcc,
2758  I, beta, h);
2759  }
2760 }
2762 
2796 template <typename T>
2797 void convolution_timestep(int n, herm_matrix<T> &C, herm_matrix<T> &A, std::complex<T> *f0,
2798  std::complex<T> *ft, herm_matrix<T> &B,
2799  integration::Integrator<T> &I, T beta, T h) {
2800  convolution_timestep<T>(n, C, A, A, f0, ft, B, B, I, beta, h);
2801 }
2802 
2804 
2841 template <typename T>
2842 void convolution(herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &Acc,
2843  std::complex<T> *f0, std::complex<T> *ft, herm_matrix<T> &B,
2844  herm_matrix<T> &Bcc, integration::Integrator<T> &I, T beta, T h) {
2845  int tstp;
2846  convolution_matsubara(C, A, f0, B, I, beta);
2847  for (tstp = 0; tstp <= C.nt(); tstp++)
2848  convolution_timestep<T>(tstp, C, A, Acc, f0, ft, B, Bcc, I, beta, h);
2849 }
2850 
2852 
2888 template <typename T>
2889 void convolution_timestep(int n, herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &Acc,
2890  function<T> &ft, herm_matrix<T> &B, herm_matrix<T> &Bcc,
2891  integration::Integrator<T> &I, T beta, T h) {
2892  int size1 = C.size1(), ntau = C.ntau(), k = I.k(), n1 = (n < k ? k : n);
2893  assert(ft.size1() == size1 && ft.nt() >= -1);
2894  if (n == -1) {
2895  convolution_matsubara(C, A, ft.ptr(-1), B, I, beta);
2896  return;
2897  }
2898  assert(n >= 0);
2899  assert(A.size1() == size1);
2900  assert(Acc.size1() == size1);
2901  assert(B.size1() == size1);
2902  assert(Bcc.size1() == size1);
2903  assert(A.ntau() == ntau);
2904  assert(Acc.ntau() == ntau);
2905  assert(B.ntau() == ntau);
2906  assert(Bcc.ntau() == ntau);
2907  assert(A.nt() >= n1);
2908  assert(Acc.nt() >= n1);
2909  assert(B.nt() >= n1);
2910  assert(Bcc.nt() >= n1);
2911  assert(ft.nt() >= n1);
2912  assert(C.nt() >= n);
2913  if (size1 == 1) {
2914  convolution_timestep_ret<T, herm_matrix<T>, 1>(n, C, A, Acc, ft.ptr(0), B, Bcc, I,
2915  h);
2916  convolution_timestep_tv<T, herm_matrix<T>, 1>(n, C, A, Acc, ft.ptr(-1), ft.ptr(0), B,
2917  Bcc, I, beta, h);
2918  convolution_timestep_les<T, herm_matrix<T>, 1>(n, C, A, Acc, ft.ptr(-1), ft.ptr(0),
2919  B, Bcc, I, beta, h);
2920  } else {
2921  convolution_timestep_ret<T, herm_matrix<T>, LARGESIZE>(n, C, A, Acc, ft.ptr(0), B,
2922  Bcc, I, h);
2923  convolution_timestep_tv<T, herm_matrix<T>, LARGESIZE>(n, C, A, Acc, ft.ptr(-1),
2924  ft.ptr(0), B, Bcc, I, beta, h);
2925  convolution_timestep_les<T, herm_matrix<T>, LARGESIZE>(
2926  n, C, A, Acc, ft.ptr(-1), ft.ptr(0), B, Bcc, I, beta, h);
2927  }
2928 }
2929 
2930 
2967 template <typename T>
2968 void convolution_timestep(int n, herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &Acc,
2969  function<T> &ft, herm_matrix<T> &B, herm_matrix<T> &Bcc,
2970  T beta, T h, int SolveOrder) {
2971  int size1 = C.size1(), ntau = C.ntau(), n1 = (n < SolveOrder ? SolveOrder : n);
2972  assert(ft.size1() == size1 && ft.nt() >= -1);
2973  if (n == -1) {
2974  convolution_matsubara(C, A, ft.ptr(-1), B, integration::I<T>(SolveOrder), beta);
2975  return;
2976  }
2977  assert(n >= 0);
2978  assert(A.size1() == size1);
2979  assert(Acc.size1() == size1);
2980  assert(B.size1() == size1);
2981  assert(Bcc.size1() == size1);
2982  assert(A.ntau() == ntau);
2983  assert(Acc.ntau() == ntau);
2984  assert(B.ntau() == ntau);
2985  assert(Bcc.ntau() == ntau);
2986  assert(A.nt() >= n1);
2987  assert(Acc.nt() >= n1);
2988  assert(B.nt() >= n1);
2989  assert(Bcc.nt() >= n1);
2990  assert(ft.nt() >= n1);
2991  assert(C.nt() >= n);
2992  if (size1 == 1) {
2993  convolution_timestep_ret<T, herm_matrix<T>, 1>(n, C, A, Acc, ft.ptr(0), B, Bcc, integration::I<T>(SolveOrder),
2994  h);
2995  convolution_timestep_tv<T, herm_matrix<T>, 1>(n, C, A, Acc, ft.ptr(-1), ft.ptr(0), B,
2996  Bcc, integration::I<T>(SolveOrder), beta, h);
2997  convolution_timestep_les<T, herm_matrix<T>, 1>(n, C, A, Acc, ft.ptr(-1), ft.ptr(0),
2998  B, Bcc, integration::I<T>(SolveOrder), beta, h);
2999  } else {
3000  convolution_timestep_ret<T, herm_matrix<T>, LARGESIZE>(n, C, A, Acc, ft.ptr(0), B,
3001  Bcc, integration::I<T>(SolveOrder), h);
3002  convolution_timestep_tv<T, herm_matrix<T>, LARGESIZE>(n, C, A, Acc, ft.ptr(-1),
3003  ft.ptr(0), B, Bcc, integration::I<T>(SolveOrder), beta, h);
3004  convolution_timestep_les<T, herm_matrix<T>, LARGESIZE>(
3005  n, C, A, Acc, ft.ptr(-1), ft.ptr(0), B, Bcc, integration::I<T>(SolveOrder), beta, h);
3006  }
3007 }
3008 
3009 
3011 
3043 template <typename T>
3044 void convolution_timestep(int n, herm_matrix<T> &C, herm_matrix<T> &A, function<T> &ft,
3045  herm_matrix<T> &B, integration::Integrator<T> &I, T beta, T h) {
3046  convolution_timestep<T>(n, C, A, A, ft, B, B, I, beta, h);
3047 }
3048 
3081 template <typename T>
3082 void convolution_timestep(int n, herm_matrix<T> &C, herm_matrix<T> &A, function<T> &ft,
3083  herm_matrix<T> &B, T beta, T h, int SolveOrder) {
3084  convolution_timestep<T>(n, C, A, A, ft, B, B, beta, h, SolveOrder);
3085 }
3087 
3122 template <typename T>
3123 void convolution(herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &Acc, function<T> &ft,
3124  herm_matrix<T> &B, herm_matrix<T> &Bcc, integration::Integrator<T> &I,
3125  T beta, T h) {
3126  int tstp;
3127  convolution_matsubara(C, A, ft.ptr(-1), B, I, beta);
3128  for (tstp = 0; tstp <= C.nt(); tstp++)
3129  convolution_timestep<T>(tstp, C, A, Acc, ft, B, Bcc, I, beta, h);
3130 }
3131 
3132 
3169 template <typename T>
3170 void convolution(herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &Acc, function<T> &ft,
3171  herm_matrix<T> &B, herm_matrix<T> &Bcc,
3172  T beta, T h, int SolveOrder) {
3173  int tstp;
3174  convolution_matsubara(C, A, ft.ptr(-1), B, integration::I<T>(SolveOrder), beta);
3175  for (tstp = 0; tstp <= C.nt(); tstp++)
3176  convolution_timestep<T>(tstp, C, A, Acc, ft, B, Bcc, beta, h, SolveOrder);
3177 }
3178 
3212 template <typename T>
3213 void convolution(herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &Acc,
3214  herm_matrix<T> &B, herm_matrix<T> &Bcc,
3215  T beta, T h, int SolveOrder) {
3216  int tstp;
3217  convolution_matsubara(C, A, B, integration::I<T>(SolveOrder), beta);
3218  for (tstp = 0; tstp <= C.nt(); tstp++)
3219  convolution_timestep<T>(tstp, C, A, Acc, B, Bcc, beta, h, SolveOrder);
3220 }
3221 
3224 // same funcions, but compute only at one timepoint
3225 
3227 
3257 template <typename T, class GG, int SIZE1>
3258 void convolution_matsubara_tau_dispatch(int m1, std::complex<T> *cc, int sizec, GG &A,
3259  std::complex<T> *f0, GG &B,
3260  integration::Integrator<T> &I, T beta) {
3261  int ntau, l, m, size1 = sizec, sb = B.element_size();
3262  std::complex<T> *bmat;
3263  T dtau;
3264  ntau = A.ntau();
3265  bmat = new std::complex<T>[sb * (ntau + 1)];
3266  if (f0 != NULL) {
3267  for (m = 0; m <= ntau; m++)
3268  element_mult<T, SIZE1>(size1, bmat + m * sb, f0, B.matptr(m));
3269  } else {
3270  for (m = 0; m <= ntau; m++)
3271  element_set<T, SIZE1>(size1, bmat + m * sb, B.matptr(m));
3272  }
3273  // compute cmat(m*dtau) = int_0^beta dx amat(tau-x) b(x)
3274  matsubara_integral_1<T, SIZE1>(size1, m1, ntau, cc, A.matptr(0), bmat, I, A.sig());
3275  delete[] bmat;
3276  // multiply by dtau:
3277  dtau = beta / ntau;
3278  for (l = 0; l < sizec * sizec; l++)
3279  cc[l] *= dtau;
3280  return;
3281 }
3283 
3320 template <typename T, class GG, int SIZE1>
3321 void convolution_timestep_les_jn_lesadv(int j, int n, std::complex<T> *cc, int sizec, GG &A,
3322  GG &Acc, std::complex<T> *ft, GG &B, GG &Bcc,
3323  integration::Integrator<T> &I, T beta, T h) {
3324  typedef std::complex<T> cplx;
3325  int k = I.get_k();
3326  int sf, sa, sb, sc, ntau, m, l, n1, size1 = sizec;
3327  T weight, dtau;
3328  cplx *atemp, *btemp, idtau, *ales, *badv, *btemp1;
3329 
3330  ntau = A.ntau();
3331  sa = A.element_size();
3332  sb = B.element_size();
3333  sc = sizec * sizec;
3334  sf = size1 * size1;
3335  dtau = beta / ntau;
3336  n1 = (n < k ? k : n);
3337 
3338  atemp = new cplx[(n1 + 1) * sa];
3339  btemp = new cplx[(n1 + 1) * sb];
3340  btemp1 = new cplx[sb];
3341  for (m = 0; m <= n1; m++) {
3342  weight = h * I.gregory_weights(n, m);
3343  if (m <= n) {
3344  element_conj<T, SIZE1>(size1, btemp1, Bcc.retptr(n, m));
3345  for (l = 0; l < sb; l++)
3346  btemp1[l] *= weight;
3347  } else {
3348  for (l = 0; l < sb; l++)
3349  btemp1[l] = -weight * B.retptr(m, n)[l];
3350  }
3351  if (ft != NULL) {
3352  element_mult<T, SIZE1>(size1, btemp + m * sb, ft + m * sf, btemp1);
3353  } else {
3354  element_set<T, SIZE1>(size1, btemp + m * sb, btemp1);
3355  }
3356  }
3357  delete[] btemp1;
3358  // for(j=0;j<=n1;j++){
3359  {
3360  for (m = 0; m <= n1; m++) {
3361  if (m < j) {
3362  element_conj<T, SIZE1>(size1, atemp + m * sa, Acc.lesptr(m, j));
3363  for (l = 0; l < sa; l++)
3364  atemp[m * sa + l] *= -1;
3365  } else {
3366  for (l = 0; l < sa; l++)
3367  atemp[m * sa + l] = A.lesptr(j, m)[l];
3368  }
3369  }
3370  ales = atemp;
3371  badv = btemp;
3372  for (m = 0; m <= n1; m++) {
3373  element_incr<T, SIZE1>(size1, cc, ales, badv);
3374  ales += sa;
3375  badv += sb;
3376  }
3377  }
3378  delete[] atemp;
3379  delete[] btemp;
3380  return;
3381 }
3383 
3420 template <typename T, class GG, int SIZE1>
3421 void convolution_timestep_les_jn_tvvt(int j, int n, std::complex<T> *cc, int sizec, GG &A,
3422  GG &Acc, std::complex<T> *f0, GG &B, GG &Bcc,
3423  integration::Integrator<T> &I, T beta, T h) {
3424  typedef std::complex<T> cplx;
3425  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1;
3426  int sa, sb, sc, ntau, m, l, n1, sig, size1 = sizec;
3427  T weight, dtau;
3428  cplx *atv, *btv, *btemp, *btemp1, *cles1, idtau;
3429 
3430  ntau = A.ntau();
3431  sa = A.element_size();
3432  sb = B.element_size();
3433  sc = sizec * sizec;
3434  dtau = beta / ntau;
3435  sig = A.sig();
3436  n1 = (n < k ? k : n);
3437  // contribution from Atv*Bvt = Atv(jh,tau) * Bcc^tv(nh,beta-tau)^* * (-Bose/Fermi)
3438  btemp = new cplx[(ntau + 1) * sb];
3439  btemp1 = new cplx[sb];
3440  idtau = cplx(0, -dtau);
3441  for (m = 0; m <= ntau; m++) {
3442  element_conj<T, SIZE1>(size1, btemp1, Bcc.tvptr(n, ntau - m));
3443  if (f0 != NULL) {
3444  element_mult<T, SIZE1>(size1, btemp + m * sb, f0, btemp1);
3445  } else {
3446  element_set<T, SIZE1>(size1, btemp + m * sb, btemp1);
3447  }
3448  }
3449  delete[] btemp1;
3450  for (l = 0; l < (ntau + 1) * sb; l++)
3451  btemp[l] *= idtau * (-(T)sig);
3452  // for(j=0;j<=n1;j++){
3453  {
3454  btv = btemp;
3455  atv = A.tvptr(j, 0);
3456  cles1 = cc;
3457  if (ntau < k2 - 1) {
3458  for (m = 0; m <= ntau; m++) {
3459  weight = I.gregory_weights(ntau, m);
3460  element_incr<T, SIZE1>(size1, cles1, weight, atv, btv);
3461  atv += sa;
3462  btv += sb;
3463  }
3464  } else {
3465  for (m = 0; m <= k; m++) {
3466  weight = I.gregory_omega(m);
3467  element_incr<T, SIZE1>(size1, cles1, weight, atv, btv);
3468  atv += sa;
3469  btv += sb;
3470  }
3471  for (m = k1; m < ntau - k; m++) {
3472  element_incr<T, SIZE1>(size1, cles1, atv, btv);
3473  atv += sa;
3474  btv += sb;
3475  }
3476  for (m = ntau - k; m <= ntau; m++) {
3477  weight = I.gregory_omega(ntau - m);
3478  element_incr<T, SIZE1>(size1, cles1, weight, atv, btv);
3479  atv += sa;
3480  btv += sb;
3481  }
3482  }
3483  }
3484  delete[] btemp;
3485  return;
3486 }
3488 
3525 template <typename T, class GG, int SIZE1>
3526 void convolution_timestep_les_jn_retles(int j, int n, std::complex<T> *cc, int sizec, GG &A,
3527  GG &Acc, std::complex<T> *ft, GG &B, GG &Bcc,
3528  integration::Integrator<T> &I, T beta, T h) {
3529  typedef std::complex<T> cplx;
3530  int k = I.get_k(), k1 = k + 1, k2 = 2 * k1, size1 = sizec;
3531  int sf, sa, sb, sc, ntau, m, l, n1;
3532  T weight, dtau;
3533  cplx *aret, *atemp, *btemp, *btemp1, *cles1, *bles, idtau;
3534 
3535  ntau = A.ntau();
3536  sa = A.element_size();
3537  sb = B.element_size();
3538  sc = sizec * sizec;
3539  sf = size1 * size1;
3540  dtau = beta / ntau;
3541  n1 = (n < k ? k : n);
3542  // contribution from Aret*Bles
3543  btemp = new cplx[(n1 + 1) * sb];
3544  btemp1 = new cplx[sb];
3545  atemp = new cplx[sa];
3546  for (m = 0; m <= n1; m++) { // btemp(m) --> B^<(m,n)
3547  if (m <= n) {
3548  for (l = 0; l < sb; l++)
3549  btemp1[l] = h * B.lesptr(m, n)[l];
3550  } else {
3551  element_conj<T, SIZE1>(size1, btemp1, Bcc.lesptr(n, m));
3552  for (l = 0; l < sb; l++)
3553  btemp1[l] *= -h;
3554  }
3555  if (ft != NULL) {
3556  element_mult<T, SIZE1>(size1, btemp + m * sb, ft + m * sf, btemp1);
3557  } else {
3558  element_set<T, SIZE1>(size1, btemp + m * sb, btemp1);
3559  }
3560  }
3561  delete[] btemp1;
3562  // for(j=0;j<=n;j++){ // compute -> cles1(j,n)
3563  {
3564  cles1 = cc;
3565  // CONTRINBUTION FROM A_RET
3566  if (j >= k2 - 1) {
3567  aret = A.retptr(j, 0);
3568  bles = btemp;
3569  for (m = 0; m <= k; m++) {
3570  weight = I.gregory_omega(m);
3571  element_incr<T, SIZE1>(size1, cles1, weight, aret, bles);
3572  bles += sb;
3573  aret += sa;
3574  }
3575  for (m = k1; m < j - k; m++) {
3576  element_incr<T, SIZE1>(size1, cles1, aret, bles);
3577  bles += sb;
3578  aret += sa;
3579  }
3580  for (m = j - k; m <= j; m++) {
3581  weight = I.gregory_omega(j - m);
3582  element_incr<T, SIZE1>(size1, cles1, weight, aret, bles);
3583  bles += sb;
3584  aret += sa;
3585  }
3586  } else {
3587  aret = A.retptr(j, 0);
3588  bles = btemp;
3589  for (m = 0; m <= j; m++) {
3590  weight = I.gregory_weights(j, m);
3591  element_incr<T, SIZE1>(size1, cles1, weight, aret, bles);
3592  bles += sb;
3593  aret += sa;
3594  }
3595  }
3596  // CONTRIBUTION FROM ACC_RET
3597  if (j < k) {
3598  for (m = j + 1; m <= k; m++) {
3599  element_conj<T, SIZE1>(size1, atemp, Acc.retptr(m, j));
3600  weight = -I.gregory_weights(j, m);
3601  bles = btemp + m * sb;
3602  element_incr<T, SIZE1>(size1, cles1, weight, atemp, bles);
3603  }
3604  }
3605  }
3606  delete[] btemp;
3607  delete[] atemp;
3608  return;
3609 }
3611 
3650 template <typename T, class GG, int SIZE1>
3651 void convolution_timestep_les_jn(int j, int n, std::complex<T> *cc, int sizec, GG &A,
3652  GG &Acc, std::complex<T> *f0, std::complex<T> *ft, GG &B,
3653  GG &Bcc, integration::Integrator<T> &I, T beta, T h) {
3654  element_set_zero<T, LARGESIZE>(sizec, cc);
3655  convolution_timestep_les_jn_tvvt<T, GG, SIZE1>(j, n, cc, sizec, A, Acc, f0, B, Bcc, I,
3656  beta, h);
3657  convolution_timestep_les_jn_lesadv<T, GG, SIZE1>(j, n, cc, sizec, A, Acc, ft, B, Bcc, I,
3658  beta, h);
3659  convolution_timestep_les_jn_retles<T, GG, SIZE1>(j, n, cc, sizec, A, Acc, ft, B, Bcc, I,
3660  beta, h);
3661  return;
3662 }
3664 // density matrix
3665 // note: !! works for both (GG=herm_matrix and GG=herm_pseudo
3666 // (because for density matrix, restricted convolution is the same as full convolution)
3667 
3669 
3706 template <typename T, class GG>
3707 void convolution_density_matrix(int n, std::complex<T> *rho, GG &A, GG &Acc, function<T> &ft,
3708  GG &B, GG &Bcc, integration::Integrator<T> &I, T beta, T h) {
3709  int size1 = A.size1(), ntau = A.ntau(), k = I.k(), n1 = (n < k && n >= 0 ? k : n);
3710  assert(ft.size1() == size1 && ft.nt() >= n1);
3711  assert(A.size1() == size1);
3712  assert(Acc.size1() == size1);
3713  assert(B.size1() == size1);
3714  assert(Bcc.size1() == size1);
3715  assert(A.ntau() == ntau);
3716  assert(Acc.ntau() == ntau);
3717  assert(B.ntau() == ntau);
3718  assert(Bcc.ntau() == ntau);
3719  assert(A.nt() >= n1);
3720  assert(Acc.nt() >= n1);
3721  assert(B.nt() >= n1);
3722  assert(Bcc.nt() >= n1);
3723  assert(ft.nt() >= n1);
3724  if (n == -1) {
3725  if (size1 == 1) {
3726  // note:
3727  convolution_matsubara_tau_dispatch<T, GG, 1>(ntau, rho, size1, A, ft.ptr(-1), B,
3728  I, beta);
3729  } else {
3730  convolution_matsubara_tau_dispatch<T, GG, LARGESIZE>(ntau, rho, size1, A,
3731  ft.ptr(-1), B, I, beta);
3732  }
3733  element_smul<T, LARGESIZE>(size1, rho, -1.0);
3734  } else {
3735  if (size1 == 1) {
3736  convolution_timestep_les_jn<T, GG, 1>(n, n, rho, size1, A, Acc, ft.ptr(-1),
3737  ft.ptr(0), B, Bcc, I, beta, h);
3738  } else {
3739  convolution_timestep_les_jn<T, GG, LARGESIZE>(
3740  n, n, rho, size1, A, Acc, ft.ptr(-1), ft.ptr(0), B, Bcc, I, beta, h);
3741  }
3742  element_smul<T, LARGESIZE>(size1, rho, std::complex<T>(0, -1.0));
3743  }
3744 }
3745 
3783 template <typename T, class GG>
3784 void convolution_density_matrix(int n, cdmatrix &rho, GG &A, GG &Acc, function<T> &ft,
3785  GG &B, GG &Bcc, T beta, T h, int SolveOrder) {
3786  int size1 = A.size1();
3787  int element_size = size1;
3788  std::complex<T> *rho_ptr;
3789  rho_ptr = new std::complex<T>[element_size];
3790 
3791  convolution_density_matrix<T, GG>(n, rho_ptr, A, Acc, ft, B, Bcc, integration::I<T>(SolveOrder), beta, h);
3792  map_ptr2matrix<T>(size1, size1, rho_ptr, rho);
3793 
3794  delete[] rho_ptr;
3795 }
3797 
3829 template <typename T, class GG>
3830 void convolution_density_matrix(int tstp, std::complex<T> *rho, GG &A, function<T> &ft,
3831  GG &B, integration::Integrator<T> &I, T beta, T h) {
3832  convolution_density_matrix<T, GG>(tstp, rho, A, A, ft, B, B, I, beta, h);
3833 }
3834 
3835 
3869 template <typename T, class GG>
3870 void convolution_density_matrix(int n, cdmatrix &rho, GG &A, function<T> &ft,
3871  GG &B, T beta, T h, int SolveOrder) {
3872  int size1 = A.size1();
3873  int element_size = size1;
3874  std::complex<T> *rho_ptr;
3875  rho_ptr = new std::complex<T>[element_size];
3876 
3877  convolution_density_matrix<T, GG>(n, rho_ptr, A, A, ft, B, B, integration::I<T>(SolveOrder), beta, h);
3878  map_ptr2matrix<T>(size1, size1, rho_ptr, rho);
3879 
3880  delete[] rho_ptr;
3881 }
3882 
3884 
3917 template <typename T, class GG>
3918 void convolution_density_matrix(int n, std::complex<T> *rho, GG &A, GG &Acc, GG &B, GG &Bcc,
3919  integration::Integrator<T> &I, T beta, T h) {
3920  int size1 = A.size1(), ntau = A.ntau(), k = I.k(), n1 = (n < k && n >= 0 ? k : n);
3921  assert(A.size1() == size1);
3922  assert(Acc.size1() == size1);
3923  assert(B.size1() == size1);
3924  assert(Bcc.size1() == size1);
3925  assert(A.ntau() == ntau);
3926  assert(Acc.ntau() == ntau);
3927  assert(B.ntau() == ntau);
3928  assert(Bcc.ntau() == ntau);
3929  assert(A.nt() >= n1);
3930  assert(Acc.nt() >= n1);
3931  assert(B.nt() >= n1);
3932  assert(Bcc.nt() >= n1);
3933  if (n == -1) {
3934  if (size1 == 1) {
3935  // note:
3936  convolution_matsubara_tau_dispatch<T, GG, 1>(ntau, rho, size1, A, NULL, B, I,
3937  beta);
3938  } else {
3939  convolution_matsubara_tau_dispatch<T, GG, LARGESIZE>(ntau, rho, size1, A, NULL,
3940  B, I, beta);
3941  }
3942  element_smul<T, LARGESIZE>(size1, rho, -1.0);
3943  } else {
3944  if (size1 == 1) {
3945  convolution_timestep_les_jn<T, GG, 1>(n, n, rho, size1, A, Acc, NULL, NULL, B,
3946  Bcc, I, beta, h);
3947  } else {
3948  convolution_timestep_les_jn<T, GG, LARGESIZE>(n, n, rho, size1, A, Acc, NULL,
3949  NULL, B, Bcc, I, beta, h);
3950  }
3951  element_smul<T, LARGESIZE>(size1, rho, std::complex<T>(0, -1.0));
3952  }
3953 }
3954 
3955 
3989 template <typename T, class GG>
3990 void convolution_density_matrix(int n, cdmatrix &rho, GG &A, GG &Acc, GG &B, GG &Bcc,
3991  T beta, T h, int SolveOrder) {
3992 
3993  int size1 = A.size1();
3994  int element_size = size1;
3995  std::complex<T> *rho_ptr;
3996  rho_ptr = new std::complex<T>[element_size];
3997 
3998  convolution_density_matrix<T, GG>(n, rho_ptr, A, Acc, B, Bcc, integration::I<T>(SolveOrder), beta, h);
3999  map_ptr2matrix<T>(size1, size1, rho_ptr, rho);
4000 
4001  delete[] rho_ptr;
4002 
4003 }
4004 
4006 
4035 template <typename T, class GG>
4036 void convolution_density_matrix(int tstp, std::complex<T> *rho, GG &A, GG &B,
4037  integration::Integrator<T> &I, T beta, T h) {
4038  convolution_density_matrix<T, GG>(tstp, rho, A, A, B, B, I, beta, h);
4039 }
4040 
4041 
4042 
4072 template <typename T, class GG>
4073 void convolution_density_matrix(int tstp, cdmatrix &rho, GG &A, GG &B,
4074  T beta, T h, int SolveOrder) {
4075  int size1 = A.size1();
4076  int element_size = size1;
4077  std::complex<T> *rho_ptr;
4078  rho_ptr = new std::complex<T>[element_size];
4079 
4080  convolution_density_matrix<T, GG>(tstp, rho_ptr, A, A, B, B, integration::I<T>(SolveOrder), beta, h);
4081  map_ptr2matrix<T>(size1, size1, rho_ptr, rho);
4082 
4083  delete[] rho_ptr;
4084 }
4085 
4087 template <typename T, class GG>
4088 void convolution_les_timediag(int tstp, cdmatrix &Cles, GG &A, GG &B,
4089  integration::Integrator<T> &I, T beta, T h){
4090  int size1=A.size1();
4091  std::complex<T> *Cles_ptr = new std::complex<T>[size1*size1];
4092  assert(Cles.rows() == size1);
4093  assert(Cles.cols() == size1);
4094 
4095  convolution_density_matrix<T, GG>(tstp, Cles_ptr, A, A, B, B, I, beta, h);
4096 
4097  for(int n1=0; n1 < size1; n1++){
4098  for(int n2=0; n2 < size1; n2++){
4099  Cles(n1,n2) = std::complex<T>(0, 1.0) * Cles_ptr[n1 * size1 + n2];
4100  }
4101  }
4102 
4103  delete Cles_ptr;
4104 
4105 }
4106 
4108 template <typename T, class GG>
4109 void convolution_density_matrix(int tstp, cdmatrix &Cles, GG &A, GG &B,
4110  integration::Integrator<T> &I, T beta, T h){
4111  int size1=A.size1();
4112  std::complex<T> *Cles_ptr = new std::complex<T>[size1*size1];
4113  assert(Cles.rows() == size1);
4114  assert(Cles.cols() == size1);
4115 
4116  convolution_density_matrix<T, GG>(tstp, Cles_ptr, A, A, B, B, I, beta, h);
4117  for(int n1=0; n1 < size1; n1++){
4118  for(int n2=0; n2 < size1; n2++){
4119  Cles(n1,n2) = std::complex<T>(0, 1.0) * Cles_ptr[n1 * size1 + n2];
4120  }
4121  }
4122  delete Cles_ptr;
4123 
4124 }
4125 
4126 /* /////////////////////////////////////////////////////////////////////////////////////////
4127 // INCREMENETAL
4130 
4162 #define CPLX std::complex<T>
4163 template <typename T, class GG, int SIZE1>
4164 void incr_convolution_mat(std::vector<bool> &mask, CPLX alpha, GG &C, GG &A, CPLX *f0, GG &B,
4165  integration::Integrator<T> &I, T beta) {
4166  int ntau = A.ntau();
4167  int size1 = A.size1();
4168  int sc = size1 * size1;
4169  bool func = (f0 == NULL ? false : true);
4170  CPLX adtau = alpha * beta * (1.0 / ntau);
4171  {
4172  // CONVOLUTION OF MATSUBARA GREENFUNCTIONS
4173  int m, sfb = sc;
4174  CPLX *ctemp = new CPLX[sc], *bmat1 = 0;
4175  CPLX *btemp;
4176  if (func) {
4177  bmat1 = new CPLX[(ntau + 1) * sfb];
4178  for (m = 0; m <= ntau; m++) {
4179  element_mult<T, SIZE1>(size1, bmat1 + m * sfb, f0, B.matptr(m));
4180  }
4181  btemp = bmat1;
4182  } else {
4183  btemp = B.matptr(0);
4184  }
4185  // t1=omp_get_wtime();
4186  for (m = 0; m <= ntau; m++) {
4187  if (mask[m]) {
4188  // compute cmat(m*dtau) = int_0^beta dx amat(tau-x) f0 b(x)
4189  matsubara_integral_1<T, SIZE1>(size1, m, ntau, ctemp, A.matptr(0), btemp, I,
4190  A.sig());
4191  element_incr<T, SIZE1>(size1, C.matptr(m), adtau, ctemp);
4192  }
4193  }
4194  // t2=omp_get_wtime();
4195  if (func)
4196  delete[] bmat1;
4197  delete[] ctemp;
4198  }
4199  // std::cout << "tid " << omp_get_thread_num() << " parallel time " << t2-t1 <<
4200  // std::endl;
4201  return;
4202 }
4204 
4242 template <typename T, class GG, int SIZE1>
4243 void incr_convolution_ret(int tstp, std::vector<bool> &mask, CPLX alpha, GG &C, GG &A,
4244  GG &Acc, CPLX *ft, GG &B, GG &Bcc, integration::Integrator<T> &I,
4245  T h) {
4246  int size1 = A.size1();
4247  int sc = size1 * size1;
4248  bool func = (ft == NULL ? false : true);
4249  CPLX adt = alpha * h;
4250  int SolveOrder = I.get_k();
4251  int n1 = (tstp >= SolveOrder ? tstp : SolveOrder);
4252  {
4253  // CONVOLUTION OF RET SECTION
4254  int j, n, l;
4255  int saf = size1 * size1;
4256  CPLX *ctemp0 = new CPLX[sc];
4257  CPLX *aret;
4258  CPLX *btmp = new CPLX[sc];
4259  CPLX *aret1 = 0;
4260  T wt;
4261  // aret[j]=Aret(tstp,j)*f(j) j=0 ... n1
4262  {
4263  if (func) {
4264  CPLX *atemp = new CPLX[sc];
4265  aret1 = new CPLX[(n1 + 1) * saf];
4266  aret = aret1;
4267  for (j = 0; j <= tstp; j++) {
4268  element_mult<T, SIZE1>(size1, aret1 + saf * j, A.retptr(tstp, j),
4269  ft + j * sc);
4270  }
4271  for (j = tstp + 1; j <= n1; j++) {
4272  element_minusconj<T, SIZE1>(size1, atemp, Acc.retptr(j, tstp));
4273  element_mult<T, SIZE1>(size1, aret1 + saf * j, atemp, ft + j * sc);
4274  }
4275  delete[] atemp;
4276  } else {
4277  if (n1 == tstp) {
4278  aret = A.retptr(tstp, 0);
4279  } else {
4280  aret1 = new CPLX[(n1 + 1) * saf];
4281  aret = aret1;
4282  for (j = 0; j <= tstp; j++)
4283  element_set<T, SIZE1>(size1, aret1 + sc * j, A.retptr(tstp, j));
4284  for (j = tstp + 1; j <= n1; j++)
4285  element_minusconj<T, SIZE1>(size1, aret1 + sc * j,
4286  Acc.retptr(j, tstp));
4287  }
4288  }
4289  }
4290  for (n = 0; n <= tstp; n++) {
4291  if (mask[n]) {
4292  // int_n^tstp dj Aret(tstp,j)Bret(j,n) ==> Cret(tstp,n)
4293  int n2 = tstp - n;
4294  element_set_zero<T, SIZE1>(size1, ctemp0);
4295  if (tstp < SolveOrder) {
4296  for (j = 0; j <= SolveOrder; j++) {
4297  wt = I.poly_integration(n, tstp, j);
4298  if (j >= n) {
4299  element_incr<T, SIZE1>(size1, ctemp0, wt, aret + j * saf,
4300  B.retptr(j, n));
4301  } else {
4302  element_minusconj<T, SIZE1>(size1, btmp, Bcc.retptr(n, j));
4303  element_incr<T, SIZE1>(size1, ctemp0, wt, aret + j * saf, btmp);
4304  }
4305  }
4306  } else if (n2 < SolveOrder) {
4307  for (l = 0; l <= SolveOrder; l++) {
4308  j = tstp - l;
4309  wt = I.gregory_weights(n2, l);
4310  if (j >= n) {
4311  element_incr<T, SIZE1>(size1, ctemp0, wt, aret + j * saf,
4312  B.retptr(j, n));
4313  } else {
4314  element_minusconj<T, SIZE1>(size1, btmp, Bcc.retptr(n, j));
4315  element_incr<T, SIZE1>(size1, ctemp0, wt, aret + j * saf, btmp);
4316  }
4317  }
4318  } else if (n2 <= 2 * SolveOrder + 2) {
4319  for (l = 0; l <= n2; l++) {
4320  j = n + l;
4321  wt = I.gregory_weights(n2, l);
4322  element_incr<T, SIZE1>(size1, ctemp0, wt, aret + j * saf,
4323  B.retptr(j, n));
4324  }
4325  } else {
4326  for (j = n; j <= n + SolveOrder; j++) {
4327  wt = I.gregory_omega(j - n);
4328  element_incr<T, SIZE1>(size1, ctemp0, wt, aret + j * saf,
4329  B.retptr(j, n));
4330  }
4331  for (j = n + SolveOrder + 1; j < tstp - SolveOrder; j++) {
4332  element_incr<T, SIZE1>(size1, ctemp0, aret + j * saf,
4333  B.retptr(j, n));
4334  }
4335  for (l = 0; l <= SolveOrder; l++) {
4336  j = tstp - l;
4337  wt = I.gregory_omega(l);
4338  element_incr<T, SIZE1>(size1, ctemp0, wt, aret + j * saf,
4339  B.retptr(j, n));
4340  }
4341  }
4342  element_incr<T, SIZE1>(size1, C.retptr(tstp, n), adt, ctemp0);
4343  }
4344  }
4345  delete[] ctemp0;
4346  delete[] btmp;
4347  if (aret1 != 0)
4348  delete[] aret1;
4349  }
4350  return;
4351 }
4353 
4395 template <typename T, class GG, int SIZE1>
4396 void incr_convolution_tv(int tstp, std::vector<bool> &mask, CPLX alpha, GG &C, GG &A,
4397  GG &Acc, CPLX *f0, CPLX *ft, GG &B, GG &Bcc,
4398  integration::Integrator<T> &I, T beta, T h) {
4399  int ntau = A.ntau();
4400  int size1 = A.size1();
4401  int sc = size1 * size1;
4402  bool func = (ft == NULL ? false : true);
4403  CPLX adt = alpha * h;
4404  CPLX adtau = alpha * beta * (1.0 / ntau);
4405  int SolveOrder = I.get_k();
4406  int n1 = (tstp >= SolveOrder ? tstp : SolveOrder);
4407  {
4408  // CONVOLUTION OF TV SECTION
4409  int j, m, n, saf = size1 * size1, sfb = size1 * size1;
4410  CPLX *ctemp1 = new CPLX[sc];
4411  CPLX *ctemp2 = new CPLX[sc];
4412  CPLX *bmat;
4413  CPLX *aret;
4414  CPLX *bmat1 = 0;
4415  CPLX *aret1 = 0;
4416  {
4417  // aret[j]=Aret(tstp,j)*f(j) j=0 ... n1
4418  if (func) {
4419  CPLX *atemp = new CPLX[sc];
4420  aret1 = new CPLX[(n1 + 1) * saf];
4421  aret = aret1;
4422  for (j = 0; j <= tstp; j++) {
4423  element_mult<T, SIZE1>(size1, aret1 + saf * j, A.retptr(tstp, j),
4424  ft + j * sc);
4425  }
4426  for (j = tstp + 1; j <= n1; j++) {
4427  element_minusconj<T, SIZE1>(size1, atemp, Acc.retptr(j, tstp));
4428  element_mult<T, SIZE1>(size1, aret1 + saf * j, atemp, ft + j * sc);
4429  }
4430  delete[] atemp;
4431  } else {
4432  if (n1 == tstp) {
4433  aret = A.retptr(tstp, 0);
4434  } else {
4435  aret1 = new CPLX[(n1 + 1) * saf];
4436  aret = aret1;
4437  for (j = 0; j <= tstp; j++)
4438  element_set<T, SIZE1>(size1, aret1 + sc * j, A.retptr(tstp, j));
4439  for (j = tstp + 1; j <= n1; j++)
4440  element_minusconj<T, SIZE1>(size1, aret1 + sc * j,
4441  Acc.retptr(j, tstp));
4442  }
4443  }
4444  // bmat[m]=f0*Bmat(m)
4445  if (func) {
4446  bmat1 = new CPLX[(ntau + 1) * sfb];
4447  bmat = bmat1;
4448  for (m = 0; m <= ntau; m++)
4449  element_mult<T, SIZE1>(size1, bmat1 + m * sfb, f0, B.matptr(m));
4450  } else {
4451  bmat = B.matptr(0);
4452  }
4453  }
4454  for (m = 0; m <= ntau; m++) {
4455  if (mask[m]) {
4456  // CONTRIBUTION FROM Atv * Bmat : very similar to computing the matsubara
4457  // convolution
4458  matsubara_integral_2<T, SIZE1>(size1, m, ntau, ctemp1, A.tvptr(tstp, 0),
4459  bmat, I, B.sig());
4460  element_smul<T, SIZE1>(size1, ctemp1, adtau);
4461  // CONTRIBUTION FROM Aret * Btv
4462  element_set_zero<T, SIZE1>(size1, ctemp2);
4463  if (tstp <= 2 * SolveOrder + 2) {
4464  for (n = 0; n <= n1; n++) {
4465  element_incr<T, SIZE1>(size1, ctemp2, I.gregory_weights(tstp, n),
4466  aret + n * saf, B.tvptr(n, m));
4467  }
4468  } else {
4469  for (n = 0; n <= SolveOrder; n++) {
4470  element_incr<T, SIZE1>(size1, ctemp2, I.gregory_omega(n),
4471  aret + n * saf, B.tvptr(n, m));
4472  }
4473  for (n = SolveOrder + 1; n < tstp - SolveOrder; n++) {
4474  element_incr<T, SIZE1>(size1, ctemp2, aret + n * saf, B.tvptr(n, m));
4475  }
4476  for (n = tstp - SolveOrder; n <= tstp; n++) {
4477  element_incr<T, SIZE1>(size1, ctemp2, I.gregory_omega(tstp - n),
4478  aret + n * saf, B.tvptr(n, m));
4479  }
4480  }
4481  element_smul<T, SIZE1>(size1, ctemp2, adt);
4482  element_incr<T, SIZE1>(size1, C.tvptr(tstp, m), ctemp1);
4483  element_incr<T, SIZE1>(size1, C.tvptr(tstp, m), ctemp2);
4484  }
4485  }
4486  if (bmat1 != 0)
4487  delete[] bmat1;
4488  if (aret1 != 0)
4489  delete[] aret1;
4490  delete[] ctemp1;
4491  delete[] ctemp2;
4492  }
4493  return;
4494 }
4496 
4538 template <typename T, class GG, int SIZE1>
4539 void incr_convolution_les(int tstp, std::vector<bool> &mask, CPLX alpha, GG &C, GG &A,
4540  GG &Acc, CPLX *f0, CPLX *ft, GG &B, GG &Bcc,
4541  integration::Integrator<T> &I, T beta, T h) {
4542  int ntau = A.ntau();
4543  int size1 = A.size1();
4544  int sc = size1 * size1;
4545  bool func = (ft == NULL ? false : true);
4546  CPLX adt = alpha * h;
4547  CPLX adtau = alpha * beta * (1.0 / ntau);
4548  int SolveOrder = I.get_k();
4549  int n1 = (tstp >= SolveOrder ? tstp : SolveOrder);
4550  {
4551  // CONVOLUTION OF LES SECTION
4552  int m, j, n, sfb = size1 * size1;
4553  CPLX *badv = new CPLX[(n1 + 1) * sfb]; // badv[j]=f(j)*Badv(j,tstp) j=0 ... n1
4554  CPLX *bles = new CPLX[(n1 + 1) * sfb]; // bles[j]=f(j)*Bles(j,tstp) j=0 ... n1
4555  CPLX *bvt = new CPLX[(ntau + 1) * sfb]; // bvt[m]=f0*Bvt(m,tstp) m=0 ... ntau
4556  CPLX *ctemp1 = new CPLX[sc];
4557  CPLX *ctemp2 = new CPLX[sc];
4558  CPLX *ctemp3 = new CPLX[sc];
4559  CPLX *atemp = new CPLX[sc];
4560  {
4561  CPLX *btemp = new CPLX[sc];
4562  if (func) {
4563  for (j = 0; j <= tstp; j++) {
4564  element_conj<T, SIZE1>(size1, btemp, Bcc.retptr(tstp, j));
4565  element_mult<T, SIZE1>(size1, badv + sfb * j, ft + sc * j, btemp);
4566  element_mult<T, SIZE1>(size1, bles + sfb * j, ft + sc * j,
4567  B.lesptr(j, tstp));
4568  }
4569  for (j = tstp + 1; j <= n1; j++) {
4570  element_mult<T, SIZE1>(size1, badv + sfb * j, ft + sc * j,
4571  B.retptr(j, tstp));
4572  element_smul<T, SIZE1>(size1, badv + sfb * j, CPLX(-1, 0));
4573  element_minusconj<T, SIZE1>(size1, btemp, Bcc.lesptr(tstp, j));
4574  element_mult<T, SIZE1>(size1, bles + sfb * j, ft + sc * j, btemp);
4575  }
4576  if (B.sig() == -1) {
4577  for (m = 0; m <= ntau; m++) {
4578  element_conj<T, SIZE1>(size1, btemp, Bcc.tvptr(tstp, ntau - m));
4579  element_mult<T, SIZE1>(size1, bvt + sfb * m, f0, btemp);
4580  }
4581  } else {
4582  // BOSE !!
4583  for (m = 0; m <= ntau; m++) {
4584  element_minusconj<T, SIZE1>(size1, btemp, Bcc.tvptr(tstp, ntau - m));
4585  element_mult<T, SIZE1>(size1, bvt + sfb * m, f0, btemp);
4586  }
4587  }
4588  } else {
4589  for (j = 0; j <= tstp; j++) {
4590  element_conj<T, SIZE1>(size1, badv + sc * j, Bcc.retptr(tstp, j));
4591  element_set<T, SIZE1>(size1, bles + sc * j, B.lesptr(j, tstp));
4592  }
4593  for (j = tstp + 1; j <= n1; j++) {
4594  element_set<T, SIZE1>(size1, badv + sc * j, B.retptr(j, tstp));
4595  element_smul<T, SIZE1>(size1, badv + sc * j, CPLX(-1, 0));
4596  element_minusconj<T, SIZE1>(size1, bles + sc * j, Bcc.lesptr(tstp, j));
4597  }
4598  if (B.sig() == -1) {
4599  for (m = 0; m <= ntau; m++) {
4600  element_conj<T, SIZE1>(size1, bvt + sc * m,
4601  Bcc.tvptr(tstp, ntau - m));
4602  }
4603  } else {
4604  // BOSE !!
4605  for (m = 0; m <= ntau; m++) {
4606  element_minusconj<T, SIZE1>(size1, bvt + sc * m,
4607  Bcc.tvptr(tstp, ntau - m));
4608  }
4609  }
4610  }
4611  delete[] btemp;
4612  }
4613  for (n = 0; n <= tstp; n++) {
4614  if (mask[n]) {
4615  T wt;
4616  /* int_0^n dj Aret(n,j)Bles(j,tstp) ==> Cles(n,tstp) */
4617  {
4618  int nup = (n > SolveOrder ? n : SolveOrder);
4619  element_set_zero<T, SIZE1>(size1, ctemp1);
4620  if (n <= 2 * SolveOrder + 2) {
4621  for (j = 0; j <= nup; j++) {
4622  wt = I.gregory_weights(n, j);
4623  if (j <= n) {
4624  element_incr<T, SIZE1>(size1, ctemp1, wt, A.retptr(n, j),
4625  bles + j * sfb);
4626  } else {
4627  element_minusconj<T, SIZE1>(size1, atemp, Acc.retptr(j, n));
4628  element_incr<T, SIZE1>(size1, ctemp1, wt, atemp,
4629  bles + j * sfb);
4630  }
4631  }
4632  } else {
4633  for (j = 0; j <= SolveOrder; j++) {
4634  wt = I.gregory_omega(j);
4635  element_incr<T, SIZE1>(size1, ctemp1, wt, A.retptr(n, j),
4636  bles + j * sfb);
4637  }
4638  for (j = SolveOrder + 1; j < n - SolveOrder; j++) {
4639  element_incr<T, SIZE1>(size1, ctemp1, A.retptr(n, j),
4640  bles + j * sfb);
4641  }
4642  for (j = n - SolveOrder; j <= n; j++) {
4643  wt = I.gregory_omega(n - j);
4644  element_incr<T, SIZE1>(size1, ctemp1, wt, A.retptr(n, j),
4645  bles + j * sfb);
4646  }
4647  }
4648  }
4649  /* int_0^beta dj Atv(n,j)Bvt(j,tstp) ==> Cles(n,tstp) */
4650  {
4651  element_set_zero<T, SIZE1>(size1, ctemp2);
4652  for (m = 0; m <= SolveOrder; m++) {
4653  wt = I.gregory_omega(m);
4654  element_incr<T, SIZE1>(size1, ctemp2, wt, A.tvptr(n, m),
4655  bvt + sfb * m);
4656  }
4657  for (m = SolveOrder + 1; m < ntau - SolveOrder; m++) {
4658  element_incr<T, SIZE1>(size1, ctemp2, A.tvptr(n, m), bvt + sfb * m);
4659  }
4660  for (m = ntau - SolveOrder; m <= ntau; m++) {
4661  wt = I.gregory_omega(ntau - m);
4662  element_incr<T, SIZE1>(size1, ctemp2, wt, A.tvptr(n, m),
4663  bvt + sfb * m);
4664  }
4665  }
4666  /* int_0^tstp dj Ales(n,j)Badv(j,tstp) ==> Cles(n,tstp) */
4667  {
4668  element_set_zero<T, SIZE1>(size1, ctemp3);
4669  if (tstp <= 2 * SolveOrder + 2) {
4670  for (j = 0; j <= n; j++) {
4671  wt = I.gregory_weights(tstp, j);
4672  element_minusconj<T, SIZE1>(size1, atemp, Acc.lesptr(j, n));
4673  element_incr<T, SIZE1>(size1, ctemp3, wt, atemp, badv + sfb * j);
4674  }
4675  for (j = n + 1; j <= n1; j++) {
4676  wt = I.gregory_weights(tstp, j);
4677  element_incr<T, SIZE1>(size1, ctemp3, wt, A.lesptr(n, j),
4678  badv + sfb * j);
4679  }
4680  } else {
4681  for (j = 0; j <= SolveOrder; j++) {
4682  wt = I.gregory_omega(j);
4683  if (j < n) {
4684  element_minusconj<T, SIZE1>(size1, atemp, Acc.lesptr(j, n));
4685  element_incr<T, SIZE1>(size1, ctemp3, wt, atemp,
4686  badv + sfb * j);
4687  } else {
4688  element_incr<T, SIZE1>(size1, ctemp3, wt, A.lesptr(n, j),
4689  badv + sfb * j);
4690  }
4691  }
4692  if (n <= SolveOrder) {
4693  for (j = SolveOrder + 1; j < tstp - SolveOrder; j++) {
4694  element_incr<T, SIZE1>(size1, ctemp3, A.lesptr(n, j),
4695  badv + sfb * j);
4696  }
4697  } else if (n < tstp - SolveOrder) {
4698  for (j = SolveOrder + 1; j < n; j++) {
4699  element_minusconj<T, SIZE1>(size1, atemp, Acc.lesptr(j, n));
4700  element_incr<T, SIZE1>(size1, ctemp3, atemp, badv + sfb * j);
4701  }
4702  for (j = n; j < tstp - SolveOrder; j++) {
4703  element_incr<T, SIZE1>(size1, ctemp3, A.lesptr(n, j),
4704  badv + sfb * j);
4705  }
4706  } else {
4707  for (j = SolveOrder + 1; j < tstp - SolveOrder; j++) {
4708  element_minusconj<T, SIZE1>(size1, atemp, Acc.lesptr(j, n));
4709  element_incr<T, SIZE1>(size1, ctemp3, atemp, badv + sfb * j);
4710  }
4711  }
4712  for (j = tstp - SolveOrder; j <= tstp; j++) {
4713  wt = I.gregory_omega(tstp - j);
4714  if (j < n) {
4715  element_minusconj<T, SIZE1>(size1, atemp, Acc.lesptr(j, n));
4716  element_incr<T, SIZE1>(size1, ctemp3, wt, atemp,
4717  badv + sfb * j);
4718  } else {
4719  element_incr<T, SIZE1>(size1, ctemp3, wt, A.lesptr(n, j),
4720  badv + sfb * j);
4721  }
4722  }
4723  }
4724  }
4725  element_incr<T, SIZE1>(size1, C.lesptr(n, tstp), adt, ctemp1);
4726  element_incr<T, SIZE1>(size1, C.lesptr(n, tstp), adt, ctemp3);
4727  element_incr<T, SIZE1>(size1, C.lesptr(n, tstp), adtau * CPLX(0.0, -1.0),
4728  ctemp2);
4729  }
4730  }
4731  delete[] ctemp1;
4732  delete[] ctemp2;
4733  delete[] ctemp3;
4734  delete[] atemp;
4735  delete[] bles;
4736  delete[] badv;
4737  delete[] bvt;
4738  }
4739  return;
4740 }
4741 
4743 // SINGLE-PROCESSOR ROUTINES ... function calls like for old version
4744 // with function object
4745 
4747 
4785 template <typename T, class GG, int SIZE1>
4786 void incr_convolution(int tstp, CPLX alpha, GG &C, GG &A, GG &Acc, CPLX *f0, CPLX *ft, GG &B,
4787  GG &Bcc, integration::Integrator<T> &I, T beta, T h) {
4788  int ntau = A.ntau();
4789  // this function is still not on top level, so no asserts!
4790  if (tstp == -1) {
4791  std::vector<bool> mask(ntau + 1, true);
4792  incr_convolution_mat<T, GG, SIZE1>(mask, alpha, C, A, f0, B, I, beta);
4793  } else if (tstp >= 0) {
4794  std::vector<bool> mask_ret(tstp + 1, true);
4795  std::vector<bool> mask_tv(ntau + 1, true);
4796  std::vector<bool> mask_les(tstp + 1, true);
4797  incr_convolution_ret<T, GG, SIZE1>(tstp, mask_ret, alpha, C, A, Acc, ft, B, Bcc, I,
4798  h);
4799  incr_convolution_tv<T, GG, SIZE1>(tstp, mask_tv, alpha, C, A, Acc, f0, ft, B, Bcc, I,
4800  beta, h);
4801  incr_convolution_les<T, GG, SIZE1>(tstp, mask_les, alpha, C, A, Acc, f0, ft, B, Bcc,
4802  I, beta, h);
4803  }
4804  return;
4805 }
4806 
4808 // NEW VERSIONS
4810 template <typename T>
4811 void convolution_timestep_new(int tstp, herm_matrix<T> &C, herm_matrix<T> &A,
4812  herm_matrix<T> &Acc, function<T> &ft, herm_matrix<T> &B,
4813  herm_matrix<T> &Bcc, integration::Integrator<T> &I, T beta,
4814  T h) {
4815  int SolveOrder = I.k();
4816  int ntmin = (tstp == -1 || tstp > SolveOrder ? tstp : SolveOrder);
4817  if (tstp < -1)
4818  return;
4819  int size1 = A.size1();
4820  std::complex<T> *fttemp;
4821 
4822  assert(C.size1() == size1);
4823  assert(ft.size1() == size1);
4824  assert(B.size1() == size1);
4825  assert(Bcc.size1() == size1);
4826  assert(Acc.size1() == size1);
4827  assert(SolveOrder > 0 && SolveOrder <= 5);
4828  assert(SolveOrder <= C.ntau());
4829  assert(ntmin <= C.nt());
4830  assert(ntmin <= A.nt());
4831  assert(ntmin <= Acc.nt());
4832  assert(ntmin <= B.nt());
4833  assert(ntmin <= Bcc.nt());
4834  assert(ntmin <= ft.nt());
4835  assert(C.ntau() == A.ntau());
4836  assert(C.ntau() == Acc.ntau());
4837  assert(C.ntau() == B.ntau());
4838  assert(C.ntau() == Bcc.ntau());
4839 
4840  C.set_timestep_zero(tstp);
4841  fttemp = (tstp == -1 ? ft.ptr(-1) : ft.ptr(0));
4842 
4843  switch (size1) {
4844  case 1:
4845  incr_convolution<T, herm_matrix<T>, 1>(tstp, CPLX(1, 0), C, A, Acc, ft.ptr(-1),
4846  fttemp, B, Bcc, integration::I<T>(SolveOrder), beta,
4847  h);
4848  break;
4849  case 2:
4850  incr_convolution<T, herm_matrix<T>, 2>(tstp, CPLX(1, 0), C, A, Acc, ft.ptr(-1),
4851  fttemp, B, Bcc, integration::I<T>(SolveOrder), beta,
4852  h);
4853  break;
4854  case 3:
4855  incr_convolution<T, herm_matrix<T>, 3>(tstp, CPLX(1, 0), C, A, Acc, ft.ptr(-1),
4856  fttemp, B, Bcc, integration::I<T>(SolveOrder), beta,
4857  h);
4858  break;
4859  case 4:
4860  incr_convolution<T, herm_matrix<T>, 4>(tstp, CPLX(1, 0), C, A, Acc, ft.ptr(-1),
4861  fttemp, B, Bcc, integration::I<T>(SolveOrder), beta,
4862  h);
4863  break;
4864  case 5:
4865  incr_convolution<T, herm_matrix<T>, 5>(tstp, CPLX(1, 0), C, A, Acc, ft.ptr(-1),
4866  fttemp, B, Bcc, integration::I<T>(SolveOrder), beta,
4867  h);
4868  break;
4869  case 6:
4870  incr_convolution<T, herm_matrix<T>, 6>(tstp, CPLX(1, 0), C, A, Acc, ft.ptr(-1),
4871  fttemp, B, Bcc, integration::I<T>(SolveOrder), beta,
4872  h);
4873  break;
4874  case 8:
4875  incr_convolution<T, herm_matrix<T>, 8>(tstp, CPLX(1, 0), C, A, Acc, ft.ptr(-1),
4876  fttemp, B, Bcc, integration::I<T>(SolveOrder), beta,
4877  h);
4878  break;
4879  default:
4880  incr_convolution<T, herm_matrix<T>, LARGESIZE>(tstp, CPLX(1, 0), C, A, Acc,
4881  ft.ptr(-1), fttemp, B, Bcc,
4882  integration::I<T>(SolveOrder), beta, h);
4883  break;
4884  }
4885 }
4887 template <typename T>
4888 void convolution_timestep_new(int n, herm_matrix<T> &C, herm_matrix<T> &A, function<T> &ft,
4889  herm_matrix<T> &B, integration::Integrator<T> &I, T beta,
4890  T h) {
4891  convolution_timestep_new<T>(n, C, A, A, ft, B, B, I, beta, h);
4892 }
4894 template <typename T>
4895 void convolution_matsubara_new(herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &B,
4896  integration::Integrator<T> &I, T beta) {
4897  convolution_timestep_new<T>(-1, C, A, A, B, B, I, beta, 0.0);
4898 }
4900 template <typename T>
4901 void convolution_matsubara_new(herm_matrix<T> &C, herm_matrix<T> &A, function<T> &ft,
4902  herm_matrix<T> &B, integration::Integrator<T> &I, T beta) {
4903  convolution_timestep_new<T>(-1, C, A, A, ft, B, B, I, beta, 0.0);
4904 }
4906 template <typename T>
4907 void convolution_new(herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &Acc,
4908  function<T> &ft, herm_matrix<T> &B, herm_matrix<T> &Bcc,
4909  integration::Integrator<T> &I, T beta, T h) {
4910  int tstp;
4911  for (tstp = -1; tstp <= C.nt(); tstp++)
4912  convolution_timestep_new<T>(tstp, C, A, Acc, ft, B, Bcc, I, beta, h);
4913 }
4915 template <typename T>
4916 void convolution_timestep_new(int tstp, herm_matrix<T> &C, herm_matrix<T> &A,
4917  herm_matrix<T> &Acc, herm_matrix<T> &B, herm_matrix<T> &Bcc,
4918  integration::Integrator<T> &I, T beta, T h) {
4919  int SolveOrder = I.k();
4920  int ntmin = (tstp == -1 || tstp > SolveOrder ? tstp : SolveOrder);
4921  if (tstp < -1)
4922  return;
4923  int size1 = A.size1();
4924 
4925  assert(C.size1() == size1);
4926  assert(B.size1() == size1);
4927  assert(Bcc.size1() == size1);
4928  assert(Acc.size1() == size1);
4929  assert(SolveOrder > 0 && SolveOrder <= 5);
4930  assert(SolveOrder <= C.ntau());
4931  assert(ntmin <= C.nt());
4932  assert(ntmin <= A.nt());
4933  assert(ntmin <= Acc.nt());
4934  assert(ntmin <= B.nt());
4935  assert(ntmin <= Bcc.nt());
4936  assert(C.ntau() == A.ntau());
4937  assert(C.ntau() == Acc.ntau());
4938  assert(C.ntau() == B.ntau());
4939  assert(C.ntau() == Bcc.ntau());
4940 
4941  C.set_timestep_zero(tstp);
4942  switch (size1) {
4943  case 1:
4944  incr_convolution<T, herm_matrix<T>, 1>(tstp, CPLX(1, 0), C, A, Acc, NULL, NULL, B,
4945  Bcc, integration::I<T>(SolveOrder), beta, h);
4946  break;
4947  case 2:
4948  incr_convolution<T, herm_matrix<T>, 2>(tstp, CPLX(1, 0), C, A, Acc, NULL, NULL, B,
4949  Bcc, integration::I<T>(SolveOrder), beta, h);
4950  break;
4951  case 3:
4952  incr_convolution<T, herm_matrix<T>, 3>(tstp, CPLX(1, 0), C, A, Acc, NULL, NULL, B,
4953  Bcc, integration::I<T>(SolveOrder), beta, h);
4954  break;
4955  case 4:
4956  incr_convolution<T, herm_matrix<T>, 4>(tstp, CPLX(1, 0), C, A, Acc, NULL, NULL, B,
4957  Bcc, integration::I<T>(SolveOrder), beta, h);
4958  break;
4959  case 5:
4960  incr_convolution<T, herm_matrix<T>, 5>(tstp, CPLX(1, 0), C, A, Acc, NULL, NULL, B,
4961  Bcc, integration::I<T>(SolveOrder), beta, h);
4962  break;
4963  case 6:
4964  incr_convolution<T, herm_matrix<T>, 6>(tstp, CPLX(1, 0), C, A, Acc, NULL, NULL, B,
4965  Bcc, integration::I<T>(SolveOrder), beta, h);
4966  break;
4967  case 8:
4968  incr_convolution<T, herm_matrix<T>, 8>(tstp, CPLX(1, 0), C, A, Acc, NULL, NULL, B,
4969  Bcc, integration::I<T>(SolveOrder), beta, h);
4970  break;
4971  default:
4972  incr_convolution<T, herm_matrix<T>, LARGESIZE>(
4973  tstp, CPLX(1, 0), C, A, Acc, NULL, NULL, B, Bcc, integration::I<T>(SolveOrder), beta, h);
4974  break;
4975  }
4976 }
4978 template <typename T>
4979 void convolution_timestep_new(int n, herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &B,
4980  integration::Integrator<T> &I, T beta, T h) {
4981  convolution_timestep_new<T>(n, C, A, A, B, B, I, beta, h);
4982 }
4984 template <typename T>
4985 void convolution_new(herm_matrix<T> &C, herm_matrix<T> &A, herm_matrix<T> &Acc,
4986  herm_matrix<T> &B, herm_matrix<T> &Bcc, integration::Integrator<T> &I,
4987  T beta, T h) {
4988  int tstp;
4989  for (tstp = -1; tstp <= C.nt(); tstp++)
4990  convolution_timestep_new<T>(tstp, C, A, Acc, B, Bcc, I, beta, h);
4991 }
4992 #undef CPLX
4993 // OPEN-MP paralellized routines
4995 #if CNTR_USE_OMP == 1
4996 
4997 #define CPLX std::complex<T>
4998 
5040 template <typename T, class GG, int SIZE1>
5041 void incr_convolution_omp(int omp_num_threads, int tstp, CPLX alpha, GG &C, GG &A, GG &Acc,
5042  CPLX *f0, CPLX *ft, GG &B, GG &Bcc, integration::Integrator<T> &I,
5043  T beta, T h) {
5044 #pragma omp parallel num_threads(omp_num_threads)
5045  {
5046  int nomp = omp_get_num_threads();
5047  int tid = omp_get_thread_num();
5048  int ntau = A.ntau(), i;
5049  // this function is still not on top level, so no asserts!
5050  if (tstp == -1) {
5051  std::vector<bool> mask(ntau + 1, false);
5052  for (i = 0; i <= ntau; i++)
5053  if (i % nomp == tid)
5054  mask[i] = true;
5055  incr_convolution_mat<T, GG, SIZE1>(mask, alpha, C, A, f0, B, I, beta);
5056  } else if (tstp >= 0) {
5057  std::vector<bool> mask_ret(tstp + 1, false);
5058  std::vector<bool> mask_tv(ntau + 1, false);
5059  std::vector<bool> mask_les(tstp + 1, false);
5060  for (i = 0; i <= tstp; i++)
5061  if (i % nomp == tid)
5062  mask_ret[i] = true;
5063  for (i = 0; i <= ntau; i++)
5064  if (i % nomp == tid)
5065  mask_tv[i] = true;
5066  for (i = 0; i <= tstp; i++)
5067  if (mask_ret[i])
5068  mask_les[tstp - i] = true;
5069  incr_convolution_ret<T, GG, SIZE1>(tstp, mask_ret, alpha, C, A, Acc, ft, B, Bcc,
5070  I, h);
5071  incr_convolution_tv<T, GG, SIZE1>(tstp, mask_tv, alpha, C, A, Acc, f0, ft, B,
5072  Bcc, I, beta, h);
5073  incr_convolution_les<T, GG, SIZE1>(tstp, mask_les, alpha, C, A, Acc, f0, ft, B,
5074  Bcc, I, beta, h);
5075  }
5076  }
5077  return;
5078 }
5080 
5117 template <typename T>
5118 void convolution_timestep_omp(int omp_num_threads, int tstp, herm_matrix<T> &C,
5119  herm_matrix<T> &A, herm_matrix<T> &Acc, function<T> &ft,
5120  herm_matrix<T> &B, herm_matrix<T> &Bcc,
5121  integration::Integrator<T> &I, T beta, T h) {
5122  int SolveOrder = I.k();
5123  int ntmin = (tstp == -1 || tstp > SolveOrder ? tstp : SolveOrder);
5124  if (tstp < -1)
5125  return;
5126  int size1 = A.size1();
5127  std::complex<T> *fttemp;
5128  assert(C.size1()==size1);
5129  assert(ft.size1()==size1);
5130  assert(B.size1()==size1);
5131  assert(Bcc.size1()==size1);
5132  assert(Acc.size1()==size1);
5133  assert(SolveOrder>=0 && SolveOrder <=5);
5134  assert(C.ntau()>=SolveOrder);
5135  assert(C.nt()>=ntmin);
5136  assert(A.nt()>=ntmin);
5137  assert(Acc.nt()>=ntmin);
5138  assert(B.nt()>=ntmin);
5139  assert(Bcc.nt()>=ntmin);
5140  assert(ft.nt()>=ntmin);
5141  assert(C.ntau()==A.ntau());
5142  assert(C.ntau()==Acc.ntau());
5143  assert(C.ntau()==B.ntau());
5144  assert(C.ntau()==Bcc.ntau());
5145  C.set_timestep_zero(tstp);
5146  fttemp = (tstp == -1 ? ft.ptr(-1) : ft.ptr(0));
5147 
5148  switch (size1) {
5149  case 1:
5150  incr_convolution_omp<T, herm_matrix<T>, 1>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5151  Acc, ft.ptr(-1), fttemp, B, Bcc,
5152  integration::I<T>(SolveOrder), beta, h);
5153  break;
5154  case 2:
5155  incr_convolution_omp<T, herm_matrix<T>, 2>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5156  Acc, ft.ptr(-1), fttemp, B, Bcc,
5157  integration::I<T>(SolveOrder), beta, h);
5158  break;
5159  case 3:
5160  incr_convolution_omp<T, herm_matrix<T>, 3>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5161  Acc, ft.ptr(-1), fttemp, B, Bcc,
5162  integration::I<T>(SolveOrder), beta, h);
5163  break;
5164  case 4:
5165  incr_convolution_omp<T, herm_matrix<T>, 4>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5166  Acc, ft.ptr(-1), fttemp, B, Bcc,
5167  integration::I<T>(SolveOrder), beta, h);
5168  break;
5169  case 5:
5170  incr_convolution_omp<T, herm_matrix<T>, 5>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5171  Acc, ft.ptr(-1), fttemp, B, Bcc,
5172  integration::I<T>(SolveOrder), beta, h);
5173  break;
5174  case 6:
5175  incr_convolution_omp<T, herm_matrix<T>, 6>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5176  Acc, ft.ptr(-1), fttemp, B, Bcc,
5177  integration::I<T>(SolveOrder), beta, h);
5178  break;
5179  case 8:
5180  incr_convolution_omp<T, herm_matrix<T>, 8>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5181  Acc, ft.ptr(-1), fttemp, B, Bcc,
5182  integration::I<T>(SolveOrder), beta, h);
5183  break;
5184  default:
5185  incr_convolution_omp<T, herm_matrix<T>, LARGESIZE>(
5186  omp_num_threads, tstp, CPLX(1, 0), C, A, Acc, ft.ptr(-1), fttemp, B, Bcc,
5187  integration::I<T>(SolveOrder), beta, h);
5188  break;
5189  }
5190 }
5191 
5192 
5230 template <typename T>
5231 void convolution_timestep_omp(int omp_num_threads, int tstp, herm_matrix<T> &C,
5232  herm_matrix<T> &A, herm_matrix<T> &Acc, function<T> &ft,
5233  herm_matrix<T> &B, herm_matrix<T> &Bcc,
5234  T beta, T h, int SolveOrder) {
5235  int ntmin = (tstp == -1 || tstp > SolveOrder ? tstp : SolveOrder);
5236  if (tstp < -1)
5237  return;
5238  int size1 = A.size1();
5239  std::complex<T> *fttemp;
5240  assert(C.size1()==size1);
5241  assert(ft.size1()==size1);
5242  assert(B.size1()==size1);
5243  assert(Bcc.size1()==size1);
5244  assert(Acc.size1()==size1);
5245  assert(SolveOrder>=0 && SolveOrder <=5);
5246  assert(C.ntau()>=SolveOrder);
5247  assert(C.nt()>=ntmin);
5248  assert(A.nt()>=ntmin);
5249  assert(Acc.nt()>=ntmin);
5250  assert(B.nt()>=ntmin);
5251  assert(Bcc.nt()>=ntmin);
5252  assert(ft.nt()>=ntmin);
5253  assert(C.ntau()==A.ntau());
5254  assert(C.ntau()==Acc.ntau());
5255  assert(C.ntau()==B.ntau());
5256  assert(C.ntau()==Bcc.ntau());
5257  C.set_timestep_zero(tstp);
5258  fttemp = (tstp == -1 ? ft.ptr(-1) : ft.ptr(0));
5259 
5260  switch (size1) {
5261  case 1:
5262  incr_convolution_omp<T, herm_matrix<T>, 1>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5263  Acc, ft.ptr(-1), fttemp, B, Bcc,
5264  integration::I<T>(SolveOrder), beta, h);
5265  break;
5266  case 2:
5267  incr_convolution_omp<T, herm_matrix<T>, 2>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5268  Acc, ft.ptr(-1), fttemp, B, Bcc,
5269  integration::I<T>(SolveOrder), beta, h);
5270  break;
5271  case 3:
5272  incr_convolution_omp<T, herm_matrix<T>, 3>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5273  Acc, ft.ptr(-1), fttemp, B, Bcc,
5274  integration::I<T>(SolveOrder), beta, h);
5275  break;
5276  case 4:
5277  incr_convolution_omp<T, herm_matrix<T>, 4>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5278  Acc, ft.ptr(-1), fttemp, B, Bcc,
5279  integration::I<T>(SolveOrder), beta, h);
5280  break;
5281  case 5:
5282  incr_convolution_omp<T, herm_matrix<T>, 5>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5283  Acc, ft.ptr(-1), fttemp, B, Bcc,
5284  integration::I<T>(SolveOrder), beta, h);
5285  break;
5286  case 6:
5287  incr_convolution_omp<T, herm_matrix<T>, 6>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5288  Acc, ft.ptr(-1), fttemp, B, Bcc,
5289  integration::I<T>(SolveOrder), beta, h);
5290  break;
5291  case 8:
5292  incr_convolution_omp<T, herm_matrix<T>, 8>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5293  Acc, ft.ptr(-1), fttemp, B, Bcc,
5294  integration::I<T>(SolveOrder), beta, h);
5295  break;
5296  default:
5297  incr_convolution_omp<T, herm_matrix<T>, LARGESIZE>(
5298  omp_num_threads, tstp, CPLX(1, 0), C, A, Acc, ft.ptr(-1), fttemp, B, Bcc,
5299  integration::I<T>(SolveOrder), beta, h);
5300  break;
5301  }
5302 }
5303 
5305 
5339 template <typename T>
5340 void convolution_timestep_omp(int omp_num_threads, int tstp, herm_matrix<T> &C,
5341  herm_matrix<T> &A, function<T> &ft, herm_matrix<T> &B,
5342  integration::Integrator<T> &I, T beta, T h) {
5343  convolution_timestep_omp<T>(omp_num_threads, tstp, C, A, A, ft, B, B, I, beta, h);
5344 }
5345 
5346 
5381 template <typename T>
5382 void convolution_timestep_omp(int omp_num_threads, int tstp, herm_matrix<T> &C,
5383  herm_matrix<T> &A, function<T> &ft, herm_matrix<T> &B,
5384  T beta, T h, int SolveOrder) {
5385  convolution_timestep_omp<T>(omp_num_threads, tstp, C, A, A, ft, B, B, beta, h, SolveOrder);
5386 }
5387 
5389 
5419 template <typename T>
5420 void convolution_matsubara_omp(int omp_num_threads, herm_matrix<T> &C, herm_matrix<T> &A,
5421  herm_matrix<T> &B, integration::Integrator<T> &I, T beta) {
5422  convolution_timestep_omp<T>(omp_num_threads, -1, C, A, A, B, B, I, beta, 0.0);
5423 }
5456 template <typename T>
5457 void convolution_matsubara_omp(int omp_num_threads, herm_matrix<T> &C, herm_matrix<T> &A,
5458  function<T> &ft, herm_matrix<T> &B,
5459  integration::Integrator<T> &I, T beta) {
5460  convolution_timestep_omp<T>(omp_num_threads, -1, C, A, A, ft, B, B, I, beta, 0.0);
5461 }
5463 
5499 template <typename T>
5500 void convolution_omp(int omp_num_threads, herm_matrix<T> &C, herm_matrix<T> &A,
5501  herm_matrix<T> &Acc, function<T> &ft, herm_matrix<T> &B,
5502  herm_matrix<T> &Bcc, integration::Integrator<T> &I, T beta, T h) {
5503  int tstp;
5504  for (tstp = -1; tstp <= C.nt(); tstp++)
5505  convolution_timestep_omp<T>(omp_num_threads, tstp, C, A, Acc, ft, B, Bcc, I, beta,
5506  h);
5507 }
5508 
5509 
5510 
5547 template <typename T>
5548 void convolution_omp(int omp_num_threads, herm_matrix<T> &C, herm_matrix<T> &A,
5549  herm_matrix<T> &Acc, function<T> &ft, herm_matrix<T> &B,
5550  herm_matrix<T> &Bcc, T beta, T h, int SolveOrder) {
5551  int tstp;
5552  for (tstp = -1; tstp <= C.nt(); tstp++)
5553  convolution_timestep_omp<T>(omp_num_threads, tstp, C, A, Acc, ft, B, Bcc, beta,
5554  h, SolveOrder);
5555 }
5556 
5557 
5559 
5595 template <typename T>
5596 void convolution_timestep_omp(int omp_num_threads, int tstp, herm_matrix<T> &C,
5597  herm_matrix<T> &A, herm_matrix<T> &Acc, herm_matrix<T> &B,
5598  herm_matrix<T> &Bcc, integration::Integrator<T> &I, T beta,
5599  T h) {
5600  int SolveOrder = I.k();
5601  int ntmin = (tstp == -1 || tstp > SolveOrder ? tstp : SolveOrder);
5602  if (tstp < -1)
5603  return;
5604  int size1 = A.size1();
5605  assert(C.size1()==size1);
5606  assert(B.size1()==size1);
5607  assert(Bcc.size1()==size1);
5608  assert(Acc.size1()==size1);
5609  assert(SolveOrder>=0 && SolveOrder <=5);
5610  assert(C.ntau()>=SolveOrder);
5611  assert(C.nt()>=ntmin);
5612  assert(A.nt()>=ntmin);
5613  assert(Acc.nt()>=ntmin);
5614  assert(B.nt()>=ntmin);
5615  assert(Bcc.nt()>=ntmin);
5616  assert(C.ntau()==A.ntau());
5617  assert(C.ntau()==Acc.ntau());
5618  assert(C.ntau()==B.ntau());
5619  assert(C.ntau()==Bcc.ntau());
5620  C.set_timestep_zero(tstp);
5621 
5622  switch (size1) {
5623  case 1:
5624  incr_convolution_omp<T, herm_matrix<T>, 1>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5625  Acc, NULL, NULL, B, Bcc,
5626  integration::I<T>(SolveOrder), beta, h);
5627  break;
5628  case 2:
5629  incr_convolution_omp<T, herm_matrix<T>, 2>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5630  Acc, NULL, NULL, B, Bcc,
5631  integration::I<T>(SolveOrder), beta, h);
5632  break;
5633  case 3:
5634  incr_convolution_omp<T, herm_matrix<T>, 3>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5635  Acc, NULL, NULL, B, Bcc,
5636  integration::I<T>(SolveOrder), beta, h);
5637  break;
5638  case 4:
5639  incr_convolution_omp<T, herm_matrix<T>, 4>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5640  Acc, NULL, NULL, B, Bcc,
5641  integration::I<T>(SolveOrder), beta, h);
5642  break;
5643  case 5:
5644  incr_convolution_omp<T, herm_matrix<T>, 5>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5645  Acc, NULL, NULL, B, Bcc,
5646  integration::I<T>(SolveOrder), beta, h);
5647  break;
5648  case 6:
5649  incr_convolution_omp<T, herm_matrix<T>, 6>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5650  Acc, NULL, NULL, B, Bcc,
5651  integration::I<T>(SolveOrder), beta, h);
5652  break;
5653  case 8:
5654  incr_convolution_omp<T, herm_matrix<T>, 8>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5655  Acc, NULL, NULL, B, Bcc,
5656  integration::I<T>(SolveOrder), beta, h);
5657  break;
5658  default:
5659  incr_convolution_omp<T, herm_matrix<T>, LARGESIZE>(omp_num_threads, tstp, CPLX(1, 0),
5660  C, A, Acc, NULL, NULL, B, Bcc,
5661  integration::I<T>(SolveOrder), beta, h);
5662  break;
5663  }
5664 }
5665 
5666 
5703 template <typename T>
5704 void convolution_timestep_omp(int omp_num_threads, int tstp, herm_matrix<T> &C,
5705  herm_matrix<T> &A, herm_matrix<T> &Acc, herm_matrix<T> &B,
5706  herm_matrix<T> &Bcc, T beta,
5707  T h, int SolveOrder) {
5708  int ntmin = (tstp == -1 || tstp > SolveOrder ? tstp : SolveOrder);
5709  if (tstp < -1)
5710  return;
5711  int size1 = A.size1();
5712  assert(C.size1()==size1);
5713  assert(B.size1()==size1);
5714  assert(Bcc.size1()==size1);
5715  assert(Acc.size1()==size1);
5716  assert(SolveOrder>=0 && SolveOrder <=5);
5717  assert(C.ntau()>=SolveOrder);
5718  assert(C.nt()>=ntmin);
5719  assert(A.nt()>=ntmin);
5720  assert(Acc.nt()>=ntmin);
5721  assert(B.nt()>=ntmin);
5722  assert(Bcc.nt()>=ntmin);
5723  assert(C.ntau()==A.ntau());
5724  assert(C.ntau()==Acc.ntau());
5725  assert(C.ntau()==B.ntau());
5726  assert(C.ntau()==Bcc.ntau());
5727  C.set_timestep_zero(tstp);
5728 
5729  switch (size1) {
5730  case 1:
5731  incr_convolution_omp<T, herm_matrix<T>, 1>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5732  Acc, NULL, NULL, B, Bcc,
5733  integration::I<T>(SolveOrder), beta, h);
5734  break;
5735  case 2:
5736  incr_convolution_omp<T, herm_matrix<T>, 2>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5737  Acc, NULL, NULL, B, Bcc,
5738  integration::I<T>(SolveOrder), beta, h);
5739  break;
5740  case 3:
5741  incr_convolution_omp<T, herm_matrix<T>, 3>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5742  Acc, NULL, NULL, B, Bcc,
5743  integration::I<T>(SolveOrder), beta, h);
5744  break;
5745  case 4:
5746  incr_convolution_omp<T, herm_matrix<T>, 4>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5747  Acc, NULL, NULL, B, Bcc,
5748  integration::I<T>(SolveOrder), beta, h);
5749  break;
5750  case 5:
5751  incr_convolution_omp<T, herm_matrix<T>, 5>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5752  Acc, NULL, NULL, B, Bcc,
5753  integration::I<T>(SolveOrder), beta, h);
5754  break;
5755  case 6:
5756  incr_convolution_omp<T, herm_matrix<T>, 6>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5757  Acc, NULL, NULL, B, Bcc,
5758  integration::I<T>(SolveOrder), beta, h);
5759  break;
5760  case 8:
5761  incr_convolution_omp<T, herm_matrix<T>, 8>(omp_num_threads, tstp, CPLX(1, 0), C, A,
5762  Acc, NULL, NULL, B, Bcc,
5763  integration::I<T>(SolveOrder), beta, h);
5764  break;
5765  default:
5766  incr_convolution_omp<T, herm_matrix<T>, LARGESIZE>(omp_num_threads, tstp, CPLX(1, 0),
5767  C, A, Acc, NULL, NULL, B, Bcc,
5768  integration::I<T>(SolveOrder), beta, h);
5769  break;
5770  }
5771 }
5772 
5773 
5774 
5775 
5808 template <typename T>
5809 void convolution_timestep_omp(int omp_num_threads, int tstp, herm_matrix<T> &C,
5810  herm_matrix<T> &A, herm_matrix<T> &B,
5811  integration::Integrator<T> &I, T beta, T h) {
5812  convolution_timestep_omp<T>(omp_num_threads, tstp, C, A, A, B, B, I, beta, h);
5813 }
5814 
5847 template <typename T>
5848 void convolution_timestep_omp(int omp_num_threads, int tstp, herm_matrix<T> &C,
5849  herm_matrix<T> &A, herm_matrix<T> &B,
5850  T beta, T h, int SolveOrder) {
5851  convolution_timestep_omp<T>(omp_num_threads, tstp, C, A, A, B, B, beta, h, SolveOrder);
5852 }
5853 
5888 template <typename T>
5889 void convolution_omp(int omp_num_threads, herm_matrix<T> &C, herm_matrix<T> &A,
5890  herm_matrix<T> &Acc, herm_matrix<T> &B, herm_matrix<T> &Bcc,
5891  integration::Integrator<T> &I, T beta, T h) {
5892  int tstp;
5893  for (tstp = -1; tstp <= C.nt(); tstp++)
5894  convolution_timestep_omp<T>(omp_num_threads, tstp, C, A, Acc, B, Bcc, I, beta, h);
5895 }
5896 
5931 template <typename T>
5932 void convolution_omp(int omp_num_threads, herm_matrix<T> &C, herm_matrix<T> &A,
5933  herm_matrix<T> &Acc, herm_matrix<T> &B, herm_matrix<T> &Bcc,
5934  T beta, T h, int SolveOrder) {
5935  int tstp;
5936  for (tstp = -1; tstp <= C.nt(); tstp++)
5937  convolution_timestep_omp<T>(omp_num_threads, tstp, C, A, Acc, B, Bcc, beta, h, SolveOrder);
5938 }
5939 
5940 #undef CPLX
5941 
5942 #endif
5943 
5944 } // namespace cntr
5945 
5946 #endif // CNTR_CONVOLUTION_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_elements.hpp
cntr_function_decl.hpp
fourier::cplx
std::complex< double > cplx
Definition: fourier.cpp:11
cntr::herm_matrix::set_timestep_zero
void set_timestep_zero(int tstp)
Sets all components at time step tstp to zero.
Definition: cntr_herm_matrix_impl.hpp:2065
cntr::herm_matrix::nt
int nt(void) const
Definition: cntr_herm_matrix_decl.hpp:71
integration::I
Integrator< T > & I(int k)
Definition: integration.hpp:506
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_convolution_decl.hpp
cntr::herm_matrix::ntau
int ntau(void) const
Definition: cntr_herm_matrix_decl.hpp:70
cntr::herm_matrix::size1
int size1(void) const
Definition: cntr_herm_matrix_decl.hpp:68
eigen_map.hpp