NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_bubble_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_BUBBLE_IMPL_H
2 #define CNTR_BUBBLE_IMPL_H
3 
4 #include "cntr_bubble_decl.hpp"
5 #include "cntr_herm_matrix_timestep_view_decl.hpp"
6 
7 namespace cntr {
8 
10 
11 // BUBBLE 1 : C(t,t') = ii * A(t,t') * B(t',t)
12 // BUBBLE 2 : C(t,t') = ii * A(t,t') * B(t,t')
13 
15 
16 // ------------- Auxiliary routines: -------------
17 
19 // Matsubara:
20 // C_{c1,c2}(tau) = - A_{a1,a2}(tau) * B_{b2,b1}(-tau)
21 // = - A_{a1,a2}(tau) * B_{b2,b1}(beta-tau) (no cc needed !!)
22 template <typename T>
23 void get_bubble_1_mat(std::complex<T> *cmat, int sc, int c1, int c2, std::complex<T> *amat,
24  int sa, int a1, int a2, std::complex<T> *bmat, int sb, int b1, int b2,
25  int sigb, int ntau) {
26  int m;
27  int a12 = a1 * sa + a2, sa2 = sa * sa;
28  int b21 = b2 * sb + b1, sb2 = sb * sb;
29  int c12 = c1 * sc + c2, sc2 = sc * sc;
30  double sig = -1.0 * sigb;
31  for (m = 0; m <= ntau; m++)
32  cmat[c12 + m * sc2] = sig * amat[m * sa2 + a12] * bmat[(ntau - m) * sb2 + b21];
33 }
35 template <typename T>
36 void
37 get_bubble_1_timestep(int tstp, std::complex<T> *cret, std::complex<T> *ctv,
38  std::complex<T> *cles, int sc, int c1, int c2, std::complex<T> *aret,
39  std::complex<T> *atv, std::complex<T> *ales, std::complex<T> *accret,
40  std::complex<T> *acctv, std::complex<T> *accles, int sa, int a1,
41  int a2, std::complex<T> *bret, std::complex<T> *btv,
42  std::complex<T> *bles, std::complex<T> *bccret, std::complex<T> *bcctv,
43  std::complex<T> *bccles, int sb, int b1, int b2, int sigb, int ntau) {
44  int m;
45  int a12 = a1 * sa + a2, a21 = a2 * sa + a1, sa2 = sa * sa;
46  int b12 = b1 * sb + b2, b21 = b2 * sb + b1, sb2 = sb * sb;
47  int c12 = c1 * sc + c2, sc2 = sc * sc;
48  std::complex<T> msigb = -1.0 * sigb, ii = std::complex<T>(0, 1.0);
49  std::complex<T> bgtr21_tt1, cgtr_tt1, cles_tt1, ales12_tt1, bgtr21_t1t, agtr12_tt1,
50  bles21_t1t, bvt21, bles21_tt1;
51  for (m = 0; m <= tstp; m++) {
52  // Bles_{21}(tstp,m) = - Bccles_{12}(m,tstp)^*
53  bles21_tt1 = -conj(bccles[m * sb2 + b12]);
54  // Bgtr_{21}(tstp,m) = Bret_{21}(tstp,m) - bles21_tt1;
55  bgtr21_tt1 = bret[m * sb2 + b21] + bles21_tt1;
56  // bgtr21_t1t = - Bccgtr_{12}(t,t1)^*
57  // = - [ Bccret_{12}(t,t1) + Bccles_{12}(t,t1) ]^*
58  // = - Bccret_{12}(t,t1)^* + Bles_{21}(t1,t)
59  bles21_t1t = bles[m * sb2 + b21];
60  bgtr21_t1t = -conj(bccret[m * sb2 + b12]) + bles21_t1t;
61  // Ales_{12}(tstp,m) = -Accles_{21}(m,tstp)^*
62  ales12_tt1 = -conj(accles[sa2 * m + a21]);
63  // Agtr_{a1,a2}(tstp,m) = Aret_{a1,a2}(tstp,m) - Ales_{a1,a2}(tstp,m)
64  agtr12_tt1 = aret[m * sa2 + a12] + ales12_tt1;
65  // Cgtr_{12}(tstp,m) = ii * Agtr_{12}(tstp,m)*Bles_{21}(m,tstp)
66  cgtr_tt1 = ii * agtr12_tt1 * bles21_t1t;
67  // Cles_{12}(tstp,m) = ii * Ales_{12}(tstp,m)*Bgtr_{21}(m,tstp)
68  cles_tt1 = ii * ales12_tt1 * bgtr21_t1t;
69  // Cret_{12}(tstp,m) = Cgtr_{12}(tstp,m) - Cles_{12}(tstp,m)
70  cret[m * sc2 + c12] = cgtr_tt1 - cles_tt1;
71  // Cles_{12}(m,tstp) = ii * Ales_{12}(m,tstp)*Bgtr_{21}(tstp,m)
72  cles[m * sc2 + c12] = ii * ales[m * sa2 + a12] * bgtr21_tt1;
73  }
74  for (m = 0; m <= ntau; m++) {
75  bvt21 = msigb * conj(bcctv[(ntau - m) * sb2 + b12]);
76  ctv[m * sc2 + c12] = ii * atv[m * sa2 + a12] * bvt21;
77  }
78 }
79 
80 // BUBBLE 2 :
81 // C(t,t') = ii * A(t,t') * B(t,t')
82 
84 // Matsubara:
85 // C_{c1,c2}(tau) = - A_{a1,a2}(tau) * B_{b1,b2}(tau)
86 template <typename T>
87 void get_bubble_2_mat(std::complex<T> *cmat, int sc, int c1, int c2, std::complex<T> *amat,
88  int sa, int a1, int a2, std::complex<T> *bmat, int sb, int b1, int b2,
89  int ntau) {
90  int m;
91  int a12 = a1 * sa + a2, sa2 = sa * sa;
92  int b12 = b1 * sb + b2, sb2 = sb * sb;
93  int c12 = c1 * sc + c2, sc2 = sc * sc;
94  for (m = 0; m <= ntau; m++)
95  cmat[c12 + m * sc2] = -amat[m * sa2 + a12] * bmat[m * sb2 + b12];
96 }
98 template <typename T>
99 void
100 get_bubble_2_timestep(int tstp, std::complex<T> *cret, std::complex<T> *ctv,
101  std::complex<T> *cles, int sc, int c1, int c2, std::complex<T> *aret,
102  std::complex<T> *atv, std::complex<T> *ales, std::complex<T> *accret,
103  std::complex<T> *acctv, std::complex<T> *accles, int sa, int a1,
104  int a2, std::complex<T> *bret, std::complex<T> *btv,
105  std::complex<T> *bles, std::complex<T> *bccret, std::complex<T> *bcctv,
106  std::complex<T> *bccles, int sb, int b1, int b2, int ntau) {
107  int m;
108  int a12 = a1 * sa + a2, a21 = a2 * sa + a1, sa2 = sa * sa;
109  int b12 = b1 * sb + b2, b21 = b2 * sb + b1, sb2 = sb * sb;
110  int c12 = c1 * sc + c2, sc2 = sc * sc;
111  std::complex<T> ii = std::complex<T>(0, 1.0);
112  std::complex<T> bgtr12_tt1, agtr12_tt1, cgtr_tt1, cles_tt1;
113  for (m = 0; m <= tstp; m++) {
114  bgtr12_tt1 = bret[m * sb2 + b12] - conj(bccles[m * sb2 + b21]);
115  agtr12_tt1 = aret[m * sa2 + a12] - conj(accles[m * sa2 + a21]);
116  // Cgtr_{12}(tstp,m) = ii * Agtr_{12}(tstp,m)*Bles_{21}(m,tstp)
117  cgtr_tt1 = ii * agtr12_tt1 * bgtr12_tt1;
118  // Cles_{12}(tstp,m) = ii * Ales_{12}(tstp,m)*Bles_{12}(tstp,m)
119  cles_tt1 = ii * conj(accles[m * sa2 + a21]) * conj(bccles[m * sb2 + b21]);
120  // Cret_{12}(tstp,m) = Cgtr_{12}(tstp,m) - Cles_{12}(tstp,m)
121  cret[m * sc2 + c12] = cgtr_tt1 - cles_tt1;
122  // Cles_{12}(m,tstp) = ii * Ales_{12}(m,tstp)*Bgtr_{12}(m,tstp)
123  cles[m * sc2 + c12] = ii * ales[m * sa2 + a12] * bles[m * sb2 + b12];
124  }
125  for (m = 0; m <= ntau; m++) {
126  ctv[m * sc2 + c12] = ii * atv[m * sa2 + a12] * btv[m * sb2 + b12];
127  }
128 }
129 
132 
134 // A,B hermitian, orbital dimension > 1:
135 // C_{c1,c2}(t1,t2) = ii * A_{a1,a2}(t1,t2) * B_{b2,b1}(t2,t1)
136 
137 template <typename T>
138 void Bubble1(int tstp, herm_matrix_timestep_view<T> &C, int c1, int c2,
139  herm_matrix_timestep_view<T> &A, herm_matrix_timestep_view<T> &Acc, int a1,
140  int a2, herm_matrix_timestep_view<T> &B, herm_matrix_timestep_view<T> &Bcc,
141  int b1, int b2) {
142  int ntau = C.ntau();
143  assert(ntau == A.ntau());
144  assert(ntau == B.ntau());
145  assert(ntau == Acc.ntau());
146  assert(ntau == Bcc.ntau());
147  assert(tstp == C.tstp_);
148  assert(tstp == B.tstp_);
149  assert(tstp == A.tstp_);
150  assert(tstp == Acc.tstp_);
151  assert(tstp == Bcc.tstp_);
152  assert(a1 <= A.size1());
153  assert(a2 <= A.size1());
154  assert(b1 <= B.size1());
155  assert(b2 <= B.size1());
156  assert(c1 <= C.size1());
157  assert(c2 <= C.size1());
158  assert(Acc.size1() == A.size1());
159  assert(Bcc.size1() == B.size1());
160 
161  if (tstp == -1) {
162  get_bubble_1_mat(C.matptr(0), C.size1(), c1, c2, A.matptr(0), A.size1(), a1, a2,
163  B.matptr(0), B.size1(), b1, b2, B.sig(), ntau);
164  } else {
165  get_bubble_1_timestep(tstp, C.retptr(0), C.tvptr(0), C.lesptr(0), C.size1(), c1, c2,
166  A.retptr(0), A.tvptr(0), A.lesptr(0), Acc.retptr(0),
167  Acc.tvptr(0), Acc.lesptr(0), A.size1(), a1, a2, B.retptr(0),
168  B.tvptr(0), B.lesptr(0), Bcc.retptr(0), Bcc.tvptr(0),
169  Bcc.lesptr(0), B.size1(), b1, b2, B.sig(), ntau);
170  }
171 }
173 template <typename T>
174 void Bubble1(int tstp, herm_matrix_timestep_view<T> &C, int c1, int c2,
175  herm_matrix_timestep_view<T> &A, int a1, int a2,
176  herm_matrix_timestep_view<T> &B, int b1, int b2) {
177  return Bubble1(tstp, C, c1, c2, A, A, a1, a2, B, B, b1, b2);
178 }
180 template <typename T>
181 void Bubble1(int tstp, herm_matrix_timestep_view<T> &C, herm_matrix_timestep_view<T> &A,
182  herm_matrix_timestep_view<T> &Acc, herm_matrix_timestep_view<T> &B,
183  herm_matrix_timestep_view<T> &Bcc) {
184  return Bubble1(tstp, C, 0, 0, A, Acc, 0, 0, B, Bcc, 0, 0);
185 }
187 template <typename T>
188 void Bubble1(int tstp, herm_matrix_timestep_view<T> &C, herm_matrix_timestep_view<T> &A,
189  herm_matrix_timestep_view<T> &B) {
190  return Bubble1(tstp, C, A, A, B, B);
191 }
193 template <class GGC, class GGA, class GGB>
194 void Bubble1(int tstp, GGC &C, int c1, int c2, GGA &A, GGA &Acc, int a1, int a2, GGB &B,
195  GGB &Bcc, int b1, int b2) {
196  herm_matrix_timestep_view<typename GGC::scalar_type> ctmp(tstp, C);
197  herm_matrix_timestep_view<typename GGA::scalar_type> atmp(tstp, A);
198  herm_matrix_timestep_view<typename GGA::scalar_type> acctmp(tstp, Acc);
199  herm_matrix_timestep_view<typename GGB::scalar_type> btmp(tstp, B);
200  herm_matrix_timestep_view<typename GGB::scalar_type> bcctmp(tstp, Bcc);
201  Bubble1(tstp, ctmp, c1, c2, atmp, acctmp, a1, a2, btmp, bcctmp, b1, b2);
202 }
204 template <class GGC, class GGA, class GGB>
205 void Bubble1(int tstp, GGC &C, int c1, int c2, GGA &A, int a1, int a2, GGB &B, int b1,
206  int b2) {
207  herm_matrix_timestep_view<typename GGC::scalar_type> ctmp(tstp, C);
208  herm_matrix_timestep_view<typename GGA::scalar_type> atmp(tstp, A);
209  herm_matrix_timestep_view<typename GGB::scalar_type> btmp(tstp, B);
210  Bubble1(tstp, ctmp, c1, c2, atmp, a1, a2, btmp, b1, b2);
211 }
213 template <class GGC, class GGA, class GGB>
214 void Bubble1(int tstp, GGC &C, GGA &A, GGA &Acc, GGB &B, GGB &Bcc) {
215  herm_matrix_timestep_view<typename GGC::scalar_type> ctmp(tstp, C);
216  herm_matrix_timestep_view<typename GGA::scalar_type> atmp(tstp, A);
217  herm_matrix_timestep_view<typename GGA::scalar_type> acctmp(tstp, Acc);
218  herm_matrix_timestep_view<typename GGB::scalar_type> btmp(tstp, B);
219  herm_matrix_timestep_view<typename GGB::scalar_type> bcctmp(tstp, Bcc);
220  Bubble1(tstp, ctmp, atmp, acctmp, btmp, bcctmp);
221 }
223 template <class GGC, class GGA, class GGB> void Bubble1(int tstp, GGC &C, GGA &A, GGB &B) {
224  herm_matrix_timestep_view<typename GGC::scalar_type> ctmp(tstp, C);
225  herm_matrix_timestep_view<typename GGA::scalar_type> atmp(tstp, A);
226  herm_matrix_timestep_view<typename GGB::scalar_type> btmp(tstp, B);
227  Bubble1(tstp, ctmp, atmp, btmp);
228 }
229 
231 // A,B hermitian, orbital dimension > 1:
232 // C_{c1,c2}(t1,t2) = ii * A_{a1,a2}(t1,t2) * B_{b1,b2}(t1,t2)
233 template <typename T>
234 void Bubble2(int tstp, herm_matrix_timestep_view<T> &C, int c1, int c2,
235  herm_matrix_timestep_view<T> &A, herm_matrix_timestep_view<T> &Acc, int a1,
236  int a2, herm_matrix_timestep_view<T> &B, herm_matrix_timestep_view<T> &Bcc,
237  int b1, int b2) {
238  int ntau = C.ntau();
239 
240  assert(ntau == A.ntau());
241  assert(ntau == B.ntau());
242  assert(ntau == Acc.ntau());
243  assert(ntau == Bcc.ntau());
244  assert(tstp == C.tstp_);
245  assert(tstp == B.tstp_);
246  assert(tstp == A.tstp_);
247  assert(tstp == Acc.tstp_);
248  assert(tstp == Bcc.tstp_);
249  assert(a1 <= A.size1());
250  assert(a2 <= A.size1());
251  assert(b1 <= B.size1());
252  assert(b2 <= B.size1());
253  assert(c1 <= C.size1());
254  assert(c2 <= C.size1());
255  assert(Acc.size1() == A.size1());
256  assert(Bcc.size1() == B.size1());
257 
258  if (tstp == -1) {
259  get_bubble_2_mat(C.matptr(0), C.size1(), c1, c2, A.matptr(0), A.size1(), a1, a2,
260  B.matptr(0), B.size1(), b1, b2, ntau);
261  } else {
262  get_bubble_2_timestep(tstp, C.retptr(0), C.tvptr(0), C.lesptr(0), C.size1(), c1, c2,
263  A.retptr(0), A.tvptr(0), A.lesptr(0), Acc.retptr(0),
264  Acc.tvptr(0), Acc.lesptr(0), A.size1(), a1, a2, B.retptr(0),
265  B.tvptr(0), B.lesptr(0), Bcc.retptr(0), Bcc.tvptr(0),
266  Bcc.lesptr(0), B.size1(), b1, b2, ntau);
267  }
268 }
270 template <typename T>
271 void Bubble2(int tstp, herm_matrix_timestep_view<T> &C, int c1, int c2,
272  herm_matrix_timestep_view<T> &A, int a1, int a2,
273  herm_matrix_timestep_view<T> &B, int b1, int b2) {
274  return Bubble2(tstp, C, c1, c2, A, A, a1, a2, B, B, b1, b2);
275 }
277 template <typename T>
278 void Bubble2(int tstp, herm_matrix_timestep_view<T> &C, herm_matrix_timestep_view<T> &A,
279  herm_matrix_timestep_view<T> &Acc, herm_matrix_timestep_view<T> &B,
280  herm_matrix_timestep_view<T> &Bcc) {
281  return Bubble2(tstp, C, 0, 0, A, Acc, 0, 0, B, Bcc, 0, 0);
282 }
284 template <typename T>
285 void Bubble2(int tstp, herm_matrix_timestep_view<T> &C, herm_matrix_timestep_view<T> &A,
286  herm_matrix_timestep_view<T> &B) {
287  return Bubble2(tstp, C, A, A, B, B);
288 }
290 template <class GGC, class GGA, class GGB>
291 void Bubble2(int tstp, GGC &C, int c1, int c2, GGA &A, GGA &Acc, int a1, int a2, GGB &B,
292  GGB &Bcc, int b1, int b2) {
293  herm_matrix_timestep_view<typename GGC::scalar_type> ctmp(tstp, C);
294  herm_matrix_timestep_view<typename GGA::scalar_type> atmp(tstp, A);
295  herm_matrix_timestep_view<typename GGA::scalar_type> acctmp(tstp, Acc);
296  herm_matrix_timestep_view<typename GGB::scalar_type> btmp(tstp, B);
297  herm_matrix_timestep_view<typename GGB::scalar_type> bcctmp(tstp, Bcc);
298  Bubble2(tstp, ctmp, c1, c2, atmp, acctmp, a1, a2, btmp, bcctmp, b1, b2);
299 }
301 template <class GGC, class GGA, class GGB>
302 void Bubble2(int tstp, GGC &C, int c1, int c2, GGA &A, int a1, int a2, GGB &B, int b1,
303  int b2) {
304  herm_matrix_timestep_view<typename GGC::scalar_type> ctmp(tstp, C);
305  herm_matrix_timestep_view<typename GGA::scalar_type> atmp(tstp, A);
306  herm_matrix_timestep_view<typename GGB::scalar_type> btmp(tstp, B);
307  Bubble2(tstp, ctmp, c1, c2, atmp, a1, a2, btmp, b1, b2);
308 }
310 template <class GGC, class GGA, class GGB>
311 void Bubble2(int tstp, GGC &C, GGA &A, GGA &Acc, GGB &B, GGB &Bcc) {
312  herm_matrix_timestep_view<typename GGC::scalar_type> ctmp(tstp, C);
313  herm_matrix_timestep_view<typename GGA::scalar_type> atmp(tstp, A);
314  herm_matrix_timestep_view<typename GGA::scalar_type> acctmp(tstp, Acc);
315  herm_matrix_timestep_view<typename GGB::scalar_type> btmp(tstp, B);
316  herm_matrix_timestep_view<typename GGB::scalar_type> bcctmp(tstp, Bcc);
317  Bubble2(tstp, ctmp, atmp, acctmp, btmp, bcctmp);
318 }
320 template <class GGC, class GGA, class GGB> void Bubble2(int tstp, GGC &C, GGA &A, GGB &B) {
321  herm_matrix_timestep_view<typename GGC::scalar_type> ctmp(tstp, C);
322  herm_matrix_timestep_view<typename GGA::scalar_type> atmp(tstp, A);
323  herm_matrix_timestep_view<typename GGB::scalar_type> btmp(tstp, B);
324  Bubble2(tstp, ctmp, atmp, btmp);
325 }
326 
327 } // namespace cntr
328 
329 #endif // CNTR_BUBBLE_IMPL_H
cntr::herm_matrix_timestep_view
Class herm_matrix_timestep_view serves for interfacing with class herm_matrix_timestep without copyi...
Definition: cntr_getset_decl.hpp:11
cntr_bubble_decl.hpp
cntr::Bubble2
void Bubble2(int tstp, GGC &C, int c1, int c2, GGA &A, GGA &Acc, int a1, int a2, GGB &B, GGB &Bcc, int b1, int b2)
Evaluate a bubble diagram ( ) from two-time contour functions at the time step; ...
Definition: cntr_bubble_impl.hpp:291
cntr_herm_matrix_timestep_view_decl.hpp
cntr::Bubble1
void Bubble1(int tstp, GGC &C, int c1, int c2, GGA &A, GGA &Acc, int a1, int a2, GGB &B, GGB &Bcc, int b1, int b2)
Evaluate a bubble diagram ( ) from two-time contour functions at the time step; ...
Definition: cntr_bubble_impl.hpp:194