NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_herm_pseudo_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_HERM_PSEUDO_IMPL_H
2 #define CNTR_HERM_PSEUDO_IMPL_H
3 
4 #include "cntr_herm_pseudo_decl.hpp"
5 //#include "cntr_exception.hpp"
6 #include "cntr_elements.hpp"
7 #include "cntr_herm_matrix_timestep_decl.hpp"
8 #include "cntr_function_decl.hpp"
9 
10 namespace cntr {
11 
12 /* #######################################################################################
13 #
14 # CONSTRUCTION/DESTRUCTION
15 #
16 ########################################################################################*/
18 template <typename T> herm_pseudo<T>::herm_pseudo() {
19  les_ = 0;
20  tv_ = 0;
21  ret_ = 0;
22  mat_ = 0;
23  ntau_ = 0;
24  nt_ = 0;
25  size1_ = 0;
26  size2_ = 0;
27  element_size_ = 0;
28  sig_ = -1;
29 }
31 template <typename T> herm_pseudo<T>::~herm_pseudo() {
32  delete[] les_;
33  delete[] ret_;
34  delete[] tv_;
35  delete[] mat_;
36 }
38 template <typename T> herm_pseudo<T>::herm_pseudo(int nt, int ntau, int size1, int sig) {
39  assert(size1 >= 0 && nt >= -1 && sig * sig == 1 && ntau >= 0);
40  nt_ = nt;
41  ntau_ = ntau;
42  sig_ = sig;
43  size1_ = size1;
44  size2_ = size1;
45  element_size_ = size1 * size1;
46  if (size1 > 0) {
47  mat_ = new cplx[(ntau_ + 1) * element_size_];
48  memset(mat_, 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
49  } else {
50  mat_ = 0;
51  }
52  if (nt >= 0 && size1 > 0) {
53  les_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
54  ret_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
55  tv_ = new cplx[(nt_ + 1) * (ntau_ + 1) * element_size_];
56  memset(les_, 0, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
57  memset(ret_, 0, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
58  memset(tv_, 0, sizeof(cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
59  } else {
60  les_ = 0;
61  tv_ = 0;
62  ret_ = 0;
63  }
64 }
66 template <typename T> unsigned long RAM_herm_pseudo(int nt, int ntau, int dim) {
67  int element_size = dim * dim;
68  size_t storage = 0;
69  storage += sizeof(typename herm_pseudo<T>::cplx) * (ntau + 1) * element_size;
70  storage +=
71  sizeof(typename herm_pseudo<T>::cplx) * ((nt + 1) * (nt + 2)) / 2 * element_size;
72  storage +=
73  sizeof(typename herm_pseudo<T>::cplx) * ((nt + 1) * (nt + 2)) / 2 * element_size;
74  storage += sizeof(typename herm_pseudo<T>::cplx) * (nt + 1) * (ntau + 1) * element_size;
75  return storage;
76 }
78 template <typename T> herm_pseudo<T>::herm_pseudo(const herm_pseudo &g) {
79  nt_ = g.nt_;
80  ntau_ = g.ntau_;
81  sig_ = g.sig_;
82  size1_ = g.size1_;
83  size2_ = g.size1_;
84  element_size_ = size1_ * size1_;
85  if (size1_ > 0) {
86  mat_ = new cplx[(ntau_ + 1) * element_size_];
87  memcpy(mat_, g.mat_, sizeof(cplx) * (ntau_ + 1) * element_size_);
88  } else {
89  mat_ = 0;
90  }
91  if (nt_ >= 0 && size1_ > 0) {
92  les_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
93  ret_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
94  tv_ = new cplx[(nt_ + 1) * (ntau_ + 1) * element_size_];
95  memcpy(les_, g.les_, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
96  memcpy(ret_, g.ret_, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
97  memcpy(tv_, g.tv_, sizeof(cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
98  } else {
99  les_ = 0;
100  ret_ = 0;
101  tv_ = 0;
102  }
103 }
104 template <typename T> herm_pseudo<T> &herm_pseudo<T>::operator=(const herm_pseudo &g) {
105  if (this == &g)
106  return *this;
107  sig_ = g.sig_;
108  if (nt_ != g.nt_ || ntau_ != g.ntau_ || size1_ != g.size1_) {
109  delete[] les_;
110  delete[] ret_;
111  delete[] tv_;
112  delete[] mat_;
113  nt_ = g.nt_;
114  ntau_ = g.ntau_;
115  size1_ = g.size1_;
116  size2_ = g.size1_;
117  element_size_ = size1_ * size1_;
118  if (size1_ > 0) {
119  mat_ = new cplx[(ntau_ + 1) * element_size_];
120  } else {
121  mat_ = 0;
122  }
123  if (size1_ > 0 && nt_ >= 0) {
124  les_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
125  ret_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
126  tv_ = new cplx[(nt_ + 1) * (ntau_ + 1) * element_size_];
127  } else {
128  les_ = 0;
129  ret_ = 0;
130  tv_ = 0;
131  }
132  }
133  if (size1_ > 0) {
134  memcpy(mat_, g.mat_, sizeof(cplx) * (ntau_ + 1) * element_size_);
135  if (nt_ >= 0) {
136  memcpy(les_, g.les_, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
137  memcpy(ret_, g.ret_, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
138  memcpy(tv_, g.tv_, sizeof(cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
139  }
140  }
141  return *this;
142 }
143 /* #######################################################################################
144 #
145 # RESIZE
146 #
147 ########################################################################################*/
149 template <typename T> void herm_pseudo<T>::resize_discard(int nt, int ntau, int size1) {
150  assert(ntau >= 0 && nt >= -1 && size1 >= 0);
151  delete[] les_;
152  delete[] ret_;
153  delete[] tv_;
154  delete[] mat_;
155  nt_ = nt;
156  ntau_ = ntau;
157  size1_ = size1;
158  size2_ = size1;
159  element_size_ = size1 * size1;
160  if (size1 > 0) {
161  mat_ = new cplx[(ntau_ + 1) * element_size_];
162  memset(mat_, 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
163  } else {
164  mat_ = 0;
165  }
166  if (nt_ >= 0 && size1_ > 0) {
167  les_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
168  ret_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
169  tv_ = new cplx[(nt_ + 1) * (ntau_ + 1) * element_size_];
170  memset(les_, 0, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
171  memset(ret_, 0, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
172  memset(tv_, 0, sizeof(cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
173  } else {
174  les_ = 0;
175  tv_ = 0;
176  ret_ = 0;
177  }
178 }
180 template <typename T> void herm_pseudo<T>::resize_nt(int nt) {
181  int nt1 = (nt_ > nt ? nt : nt_);
182  cplx *ret, *les, *tv;
183  assert(nt >= -1);
184  nt_ = nt;
185  if (size1_ == 0)
186  return;
187  if (nt_ >= 0) {
188  les = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
189  ret = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
190  tv = new cplx[(nt_ + 1) * (ntau_ + 1) * element_size_];
191  memset(les, 0, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
192  memset(ret, 0, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
193  memset(tv, 0, sizeof(cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
194  if (nt1 >= 0) {
195  memcpy(les, les_, sizeof(cplx) * ((nt1 + 1) * (nt1 + 2)) / 2 * element_size_);
196  memcpy(ret, ret_, sizeof(cplx) * ((nt1 + 1) * (nt1 + 2)) / 2 * element_size_);
197  memcpy(tv, tv_, sizeof(cplx) * (nt1 + 1) * (ntau_ + 1) * element_size_);
198  }
199  } else {
200  les = 0;
201  ret = 0;
202  tv = 0;
203  }
204  delete[] les_;
205  delete[] ret_;
206  delete[] tv_;
207  les_ = les;
208  ret_ = ret;
209  tv_ = tv;
210 }
212 template <typename T> void herm_pseudo<T>::resize(int nt, int ntau, int size1) {
213  // std::cout << "herm_pseudo<T>::" << __FUNCTION__ << " " << nt << " " << ntau << " " <<
214  // size1 << std::endl;
215  if (ntau == ntau_ && size1_ == size1)
216  resize_nt(nt);
217  else
218  resize_discard(nt, ntau, size1);
219 }
221 template <typename T> void herm_pseudo<T>::clear(void) {
222  if (size1_ == 0)
223  return;
224  memset(mat_, 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
225  if (nt_ >= 0) {
226  memset(les_, 0, sizeof(cplx) * ((nt_ + 2) * (nt_ + 1)) / 2 * element_size_);
227  memset(ret_, 0, sizeof(cplx) * ((nt_ + 2) * (nt_ + 1)) / 2 * element_size_);
228  memset(tv_, 0, sizeof(cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
229  }
230 }
231 /* #######################################################################################
232 #
233 # RAW POINTERS TO ELEMENTS
234 #
235 ########################################################################################*/
237 template <typename T> inline std::complex<T> *herm_pseudo<T>::lesptr(int t, int t1) {
238  assert(t >= 0 && t1 >= 0 && t <= t1 && t1 <= nt_);
239  return les_ + ((t1 * (t1 + 1)) / 2 + t) * element_size_;
240 }
242 template <typename T> inline std::complex<T> *herm_pseudo<T>::retptr(int t, int t1) {
243  assert(t >= 0 && t1 >= 0 && t <= nt_ && t1 <= t);
244  return ret_ + ((t * (t + 1)) / 2 + t1) * element_size_;
245 }
247 template <typename T> inline std::complex<T> *herm_pseudo<T>::tvptr(int t, int tau) {
248  assert(t >= 0 && tau >= 0 && t <= nt_ && tau <= ntau_);
249  return tv_ + (t * (ntau_ + 1) + tau) * element_size_;
250 }
252 template <typename T> inline std::complex<T> *herm_pseudo<T>::matptr(int tau) {
253  assert(tau >= 0 && tau <= ntau_);
254  return mat_ + tau * element_size_;
255 }
256 /* #######################################################################################
257 #
258 # READING ELEMENTS TO ANY MATRIX TYPE
259 # OR TO COMPLEX NUMBERS (then only the (0,0) element is addressed for dim>0)
260 # note: these are not efficient, in particular the conjugation
261 #
262 ########################################################################################*/
264 #define herm_pseudo_READ_ELEMENT \
265  { \
266  int r, s, dim = size1_; \
267  M.resize(dim, dim); \
268  for (r = 0; r < dim; r++) \
269  for (s = 0; s < dim; s++) \
270  M(r, s) = x[r * dim + s]; \
271  }
272 #define herm_pseudo_READ_ELEMENT_MINUS_CONJ \
274  { \
275  cplx w; \
276  int r, s, dim = size1_; \
277  M.resize(dim, dim); \
278  for (r = 0; r < dim; r++) \
279  for (s = 0; s < dim; s++) { \
280  w = x[s * dim + r]; \
281  M(r, s) = std::complex<T>(-w.real(), w.imag()); \
282  } \
283  }
284 template <typename T>
286 template <class Matrix>
287 void herm_pseudo<T>::get_les(int i, int j, Matrix &M) {
288  cplx *x;
289  if (i <= j) {
290  x = lesptr(i, j);
291  herm_pseudo_READ_ELEMENT
292  } else {
293  x = lesptr(j, i);
294  herm_pseudo_READ_ELEMENT_MINUS_CONJ
295  }
296 }
298 template <typename T>
299 template <class Matrix>
300 void herm_pseudo<T>::get_ret(int i, int j, Matrix &M) {
301  cplx *x;
302  if (i >= j) {
303  x = retptr(i, j);
304  herm_pseudo_READ_ELEMENT
305  } else {
306  x = retptr(j, i);
307  herm_pseudo_READ_ELEMENT_MINUS_CONJ
308  }
309 }
311 template <typename T>
312 template <class Matrix>
313 void herm_pseudo<T>::get_tv(int i, int j, Matrix &M) {
314  cplx *x = tvptr(i, j);
315  herm_pseudo_READ_ELEMENT
316 }
318 template <typename T>
319 template <class Matrix>
320 void herm_pseudo<T>::get_vt(int i, int j, Matrix &M) {
321  cplx *x = tvptr(j, ntau_ - i);
322  herm_pseudo_READ_ELEMENT_MINUS_CONJ if (sig_ == -1) M = -M;
323 }
325 template <typename T>
326 template <class Matrix>
327 void herm_pseudo<T>::get_mat(int i, Matrix &M) {
328  cplx *x = matptr(i);
329  herm_pseudo_READ_ELEMENT
330 }
332 template <typename T>
333 template <class Matrix>
334 void herm_pseudo<T>::get_matminus(int i, Matrix &M) {
335  cplx *x = matptr(ntau_ - i);
336  herm_pseudo_READ_ELEMENT if (sig_ == -1) M = -M;
337 }
339 template <typename T>
340 template <class Matrix>
341 void herm_pseudo<T>::get_gtr(int i, int j, Matrix &M) {
342  get_ret(i, j, M);
343 }
345 template <typename T> inline void herm_pseudo<T>::get_ret(int i, int j, cplx &x) {
346  if (i >= j)
347  x = *retptr(i, j);
348  else {
349  x = *retptr(j, i);
350  x = -std::conj(x);
351  }
352 }
354 template <typename T> inline void herm_pseudo<T>::get_les(int i, int j, cplx &x) {
355  if (i <= j)
356  x = *lesptr(i, j);
357  else {
358  x = *lesptr(j, i);
359  x = -std::conj(x);
360  }
361 }
363 template <typename T> inline void herm_pseudo<T>::get_tv(int i, int j, cplx &x) {
364  x = *tvptr(i, j);
365 }
367 template <typename T> inline void herm_pseudo<T>::get_vt(int i, int j, cplx &x) {
368  x = *tvptr(j, ntau_ - i);
369  if (sig_ == -1)
370  x = std::conj(x);
371  else
372  x = -std::conj(x);
373 }
375 template <typename T> inline void herm_pseudo<T>::get_mat(int i, cplx &x) { x = *matptr(i); }
377 template <typename T> inline void herm_pseudo<T>::get_matminus(int i, cplx &x) {
378  x = *matptr(ntau_ - i);
379  if (sig_ == -1)
380  x = -x;
381 }
383 template <typename T> inline void herm_pseudo<T>::get_gtr(int i, int j, cplx &x) {
384  get_ret(i, j, x);
385 }
386 /* #######################################################################################
387 #
388 # WRITING ELEMENTS FROM ANY MATRIX TYPE
389 # OR FROM COMPLEX NUMBERS (then only the (0,0) element is addressed for dim>0)
390 #
391 ########################################################################################*/
393 #define herm_pseudo_SET_ELEMENT_MATRIX \
394  { \
395  int r, s, dim = size1_; \
396  assert(M.rows() == dim && M.cols() == dim); \
397  for (r = 0; r < dim; r++) \
398  for (s = 0; s < dim; s++) \
399  x[r * dim + s] = M(r, s); \
400  }
401 template <typename T>
403 template <class Matrix>
404 void herm_pseudo<T>::set_ret(int i, int j, Matrix &M) {
405  cplx *x = retptr(i, j);
406  herm_pseudo_SET_ELEMENT_MATRIX
407 }
409 template <typename T>
410 template <class Matrix>
411 void herm_pseudo<T>::set_les(int i, int j, Matrix &M) {
412  cplx *x = lesptr(i, j);
413  herm_pseudo_SET_ELEMENT_MATRIX
414 }
416 template <typename T>
417 template <class Matrix>
418 void herm_pseudo<T>::set_tv(int i, int j, Matrix &M) {
419  cplx *x = tvptr(i, j);
420  herm_pseudo_SET_ELEMENT_MATRIX
421 }
423 template <typename T>
424 template <class Matrix>
425 void herm_pseudo<T>::set_mat(int i, Matrix &M) {
426  cplx *x = matptr(i);
427  herm_pseudo_SET_ELEMENT_MATRIX
428 }
430 template <typename T> inline void herm_pseudo<T>::set_les(int i, int j, cplx x) {
431  *lesptr(i, j) = x;
432 }
434 template <typename T> inline void herm_pseudo<T>::set_ret(int i, int j, cplx x) {
435  *retptr(i, j) = x;
436 }
438 template <typename T> inline void herm_pseudo<T>::set_tv(int i, int j, cplx x) {
439  *tvptr(i, j) = x;
440 }
442 template <typename T> inline void herm_pseudo<T>::set_mat(int i, cplx x) { *matptr(i) = x; }
444 template <typename T> std::complex<T> herm_pseudo<T>::density_matrix(int tstp) {
445  cplx x1;
446  if (tstp == -1) {
447  get_mat(ntau_, x1);
448  return -x1;
449  } else {
450  get_les(tstp, tstp, x1);
451  return std::complex<T>(0.0, sig_) * x1;
452  }
453 }
454 
455 // should work with Eigen
457 template <typename T>
458 template <class Matrix>
459 void herm_pseudo<T>::density_matrix(int tstp, Matrix &M) {
460  if (tstp == -1) {
461  get_mat(ntau_, M);
462  M *= (-1.0);
463  } else {
464  get_les(tstp, tstp, M);
465  M *= std::complex<T>(0.0, 1.0 * sig_);
466  }
467 }
468 
469 /* #######################################################################################
470 #
471 # INPUT/OUTPUT FROM/TO FILES
472 #
473 ########################################################################################*/
475 template <typename T> void herm_pseudo<T>::print_to_file(const char *file, int precision) {
476  int i, j, l, sg = element_size_;
477  std::ofstream out;
478  out.open(file, std::ios::out);
479  out.precision(precision);
480  out << "# " << nt_ << " " << ntau_ << " " << size1_ << " "
481  << " " << sig_ << std::endl;
482  for (j = 0; j <= ntau_; j++) {
483  out << "mat: " << j;
484  for (l = 0; l < sg; l++)
485  out << " " << matptr(j)[l].real() << " " << matptr(j)[l].imag();
486  out << std::endl;
487  }
488  out << std::endl;
489  if (nt_ >= 0) {
490  for (i = 0; i <= nt_; i++) {
491  for (j = 0; j <= i; j++) {
492  out << "ret: " << i << " " << j;
493  for (l = 0; l < sg; l++)
494  out << " " << retptr(i, j)[l].real() << " " << retptr(i, j)[l].imag();
495  out << std::endl;
496  }
497  out << std::endl;
498  }
499  out << std::endl;
500  for (i = 0; i <= nt_; i++) {
501  for (j = 0; j <= ntau_; j++) {
502  out << "tv: " << i << " " << j;
503  for (l = 0; l < sg; l++)
504  out << " " << tvptr(i, j)[l].real() << " " << tvptr(i, j)[l].imag();
505  out << std::endl;
506  }
507  out << std::endl;
508  }
509  out << std::endl;
510  for (j = 0; j <= nt_; j++) {
511  for (i = 0; i <= j; i++) {
512  out << "les: " << i << " " << j;
513  for (l = 0; l < sg; l++)
514  out << " " << lesptr(i, j)[l].real() << " " << lesptr(i, j)[l].imag();
515  out << std::endl;
516  }
517  out << std::endl;
518  }
519  out << std::endl;
520  }
521  out.close();
522 }
524 template <typename T> void herm_pseudo<T>::read_from_file(const char *file) {
525  int i, n, m, j, l, size1, sg, sig;
526  double real, imag;
527  std::string s;
528  std::ifstream out;
529  out.open(file, std::ios::in);
530  if (!(out >> s >> n >> m >> size1 >> sig)) {
531  std::cerr << "read G from file " << file << " error in file" << std::endl;
532  abort();
533  }
534  if (n > nt_ || m != ntau_ || size1 != size1_)
535  resize(n, m, size1);
536  sig_ = sig;
537  sg = element_size_;
538  for (j = 0; j <= ntau_; j++) {
539  out >> s >> s;
540  for (l = 0; l < sg; l++) {
541  if (!(out >> real >> imag)) {
542  std::cerr << "read G from file " << file << " error at mat (" << j << ")"
543  << std::endl;
544  abort();
545  }
546  matptr(j)[l] = std::complex<T>(real, imag);
547  }
548  }
549  if (n >= 0) {
550  for (i = 0; i <= n; i++) {
551  for (j = 0; j <= i; j++) {
552  out >> s >> s >> s;
553  for (l = 0; l < sg; l++) {
554  if (!(out >> real >> imag)) {
555  std::cerr << "read G from file " << file << " error at ret (" << i
556  << "," << j << ")" << std::endl;
557  abort();
558  }
559  retptr(i, j)[l] = std::complex<T>(real, imag);
560  }
561  }
562  }
563  for (i = 0; i <= n; i++) {
564  for (j = 0; j <= ntau_; j++) {
565  out >> s >> s >> s;
566  for (l = 0; l < sg; l++) {
567  if (!(out >> real >> imag)) {
568  std::cerr << "read G from file " << file << " error at tv (" << i
569  << "," << j << ")" << std::endl;
570  abort();
571  }
572  tvptr(i, j)[l] = std::complex<T>(real, imag);
573  }
574  }
575  }
576  for (j = 0; j <= n; j++) {
577  for (i = 0; i <= j; i++) {
578  out >> s >> s >> s;
579  for (l = 0; l < sg; l++) {
580  if (!(out >> real >> imag)) {
581  std::cerr << "read G from file " << file << " error at les (" << i
582  << "," << j << ")" << std::endl;
583  abort();
584  }
585  lesptr(i, j)[l] = std::complex<T>(real, imag);
586  }
587  }
588  }
589  }
590  out.close();
591 }
592 #if CNTR_USE_HDF5 == 1
593 // identical from herm_matrix<T>
595 template <typename T> void herm_pseudo<T>::write_to_hdf5(hid_t group_id) {
596  store_int_attribute_to_hid(group_id, std::string("ntau"), ntau_);
597  store_int_attribute_to_hid(group_id, std::string("nt"), nt_);
598  store_int_attribute_to_hid(group_id, std::string("sig"), sig_);
599  store_int_attribute_to_hid(group_id, std::string("size1"), size1_);
600  store_int_attribute_to_hid(group_id, std::string("size2"), size2_);
601  store_int_attribute_to_hid(group_id, std::string("element_size"), element_size_);
602  hsize_t len_shape = 3, shape[3];
603  shape[1] = size1_;
604  shape[2] = size2_;
605  if (nt_ > -2) {
606  shape[0] = ntau_ + 1;
607  store_cplx_array_to_hid(group_id, std::string("mat"), matptr(0), shape, len_shape);
608  }
609  if (nt_ > -1) {
610  shape[0] = ((nt_ + 1) * (nt_ + 2)) / 2;
611  // CHECK: implement store_cplx_array_to_hid with template typename T
612  store_cplx_array_to_hid(group_id, std::string("ret"), retptr(0, 0), shape,
613  len_shape);
614  store_cplx_array_to_hid(group_id, std::string("les"), lesptr(0, 0), shape,
615  len_shape);
616  shape[0] = (nt_ + 1) * (ntau_ + 1);
617  store_cplx_array_to_hid(group_id, std::string("tv"), tvptr(0, 0), shape, len_shape);
618  }
619 }
621 template <typename T>
622 void herm_pseudo<T>::write_to_hdf5(hid_t group_id, const char *groupname) {
623  hid_t sub_group_id = create_group(group_id, groupname);
624  this->write_to_hdf5(sub_group_id);
625  close_group(sub_group_id);
626 }
628 template <typename T>
629 void herm_pseudo<T>::write_to_hdf5(const char *filename, const char *groupname) {
630  hid_t file_id = open_hdf5_file(filename);
631  this->write_to_hdf5(file_id, groupname);
632  close_hdf5_file(file_id);
633 }
635 template <typename T> void herm_pseudo<T>::read_from_hdf5(hid_t group_id) {
636  // -- Read dimensions
637  int nt = read_primitive_type<int>(group_id, "nt");
638  int ntau = read_primitive_type<int>(group_id, "ntau");
639  int sig = read_primitive_type<int>(group_id, "sig");
640  int size1 = read_primitive_type<int>(group_id, "size1");
641  // RESIZE G
642  this->resize(nt, ntau, size1);
643  sig_ = sig;
644  if (nt > -2) {
645  hsize_t mat_size = (ntau + 1) * element_size_;
646  read_primitive_type_array(group_id, "mat", mat_size, matptr(0));
647  }
648  if (nt > -1) {
649  hsize_t ret_size = ((nt + 1) * (nt + 2)) / 2 * element_size_;
650  hsize_t les_size = ((nt + 1) * (nt + 2)) / 2 * element_size_;
651  hsize_t tv_size = ((nt + 1) * (ntau + 1)) * element_size_;
652  read_primitive_type_array(group_id, "ret", ret_size, retptr(0, 0));
653  read_primitive_type_array(group_id, "les", les_size, lesptr(0, 0));
654  read_primitive_type_array(group_id, "tv", tv_size, tvptr(0, 0));
655  }
656 }
658 template <typename T>
659 void herm_pseudo<T>::read_from_hdf5(hid_t group_id, const char *groupname) {
660  hid_t sub_group_id = open_group(group_id, groupname);
661  this->read_from_hdf5(sub_group_id);
662  close_group(sub_group_id);
663 }
665 template <typename T>
666 void herm_pseudo<T>::read_from_hdf5(const char *filename, const char *groupname) {
667  hid_t file_id = read_hdf5_file(filename);
668  this->read_from_hdf5(file_id, groupname);
669  close_hdf5_file(file_id);
670 }
672 template <typename T> void herm_pseudo<T>::read_from_hdf5(int nt1, hid_t group_id) {
673  herm_pseudo<T> gtmp;
674  gtmp.read_from_hdf5(group_id);
675  assert(nt1 >= -1 && nt1 <= gtmp.nt());
676  assert(nt1 >= -1 && nt1 <= nt_);
677  assert(gtmp.size1() == size1_);
678  assert(gtmp.element_size() == element_size_);
679  assert(gtmp.ntau() == ntau_);
680 
681  for (int n = -1; n <= nt1; n++)
682  this->set_timestep(n, gtmp);
683 }
685 template <typename T>
686 void herm_pseudo<T>::read_from_hdf5(int nt1, hid_t group_id, const char *groupname) {
687  hid_t sub_group_id = open_group(group_id, groupname);
688  this->read_from_hdf5(nt1, sub_group_id);
689  close_group(sub_group_id);
690 }
692 template <typename T>
693 void herm_pseudo<T>::read_from_hdf5(int nt1, const char *filename, const char *groupname) {
694  hid_t file_id = read_hdf5_file(filename);
695  this->read_from_hdf5(nt1, file_id, groupname);
696  close_hdf5_file(file_id);
697 }
698 #endif
699 
700 /* #######################################################################################
701 #
702 # SIMPLE OPERATIONS ON TIMESTEPS
703 # NOTE: tstp IS A PHYSICAL TIME, tstp=-1 is the matsubara branch
704 #
705 ########################################################################################*/
707 template <typename T> void herm_pseudo<T>::set_timestep_zero(int tstp) {
708  assert(tstp >= -1 && tstp <= nt_);
709  if (tstp == -1) {
710  memset(matptr(0), 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
711  } else {
712  memset(retptr(tstp, 0), 0, sizeof(cplx) * (tstp + 1) * element_size_);
713  memset(tvptr(tstp, 0), 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
714  memset(lesptr(0, tstp), 0, sizeof(cplx) * (tstp + 1) * element_size_);
715  }
716 }
718 template <typename T> void herm_pseudo<T>::set_timestep(int tstp, herm_pseudo<T> &g1) {
719  assert(tstp >= -1 && tstp <= nt_ && tstp <= g1.nt());
720  assert(g1.size1() == size1_);
721  assert(g1.ntau() == ntau_);
722  if (tstp == -1) {
723  memcpy(mat_, g1.mat_, sizeof(cplx) * (ntau_ + 1) * element_size_);
724  } else {
725  memcpy(retptr(tstp, 0), g1.retptr(tstp, 0),
726  sizeof(cplx) * (tstp + 1) * element_size_);
727  memcpy(tvptr(tstp, 0), g1.tvptr(tstp, 0),
728  sizeof(cplx) * (ntau_ + 1) * element_size_);
729  memcpy(lesptr(0, tstp), g1.lesptr(0, tstp),
730  sizeof(cplx) * (tstp + 1) * element_size_);
731  }
732 }
734 template <typename T>
735 void herm_pseudo<T>::set_timestep(int tstp, herm_matrix_timestep<T> &timestep) {
736  cplx *x = timestep.data_;
737  assert(tstp >= -1 && tstp <= nt_);
738  assert(timestep.tstp_ == tstp && timestep.ntau_ == ntau_ && timestep.size1_ == size1_);
739  if (tstp == -1) {
740  memcpy(mat_, x, sizeof(cplx) * (ntau_ + 1) * element_size_);
741  } else {
742  memcpy(retptr(tstp, 0), x, sizeof(cplx) * (tstp + 1) * element_size_);
743  memcpy(tvptr(tstp, 0), x + (tstp + 1) * element_size_,
744  sizeof(cplx) * (ntau_ + 1) * element_size_);
745  memcpy(lesptr(0, tstp), x + (tstp + 1 + ntau_ + 1) * element_size_,
746  sizeof(cplx) * (tstp + 1) * element_size_);
747  }
748 }
750 template <typename T>
751 void herm_pseudo<T>::get_timestep(int tstp, herm_matrix_timestep<T> &timestep) {
752  int len = (2 * (tstp + 1) + ntau_ + 1) * element_size_;
753  cplx *x;
754  assert(tstp >= -1 && tstp <= nt_);
755  if (timestep.total_size_ < len)
756  timestep.resize(tstp, ntau_, size1_);
757  x = timestep.data_;
758  timestep.tstp_ = tstp;
759  timestep.ntau_ = ntau_;
760  timestep.size1_ = size1_;
761  if (tstp == -1) {
762  memcpy(x, mat_, sizeof(cplx) * (ntau_ + 1) * element_size_);
763  } else {
764  memcpy(x, retptr(tstp, 0), sizeof(cplx) * (tstp + 1) * element_size_);
765  memcpy(x + (tstp + 1) * element_size_, tvptr(tstp, 0),
766  sizeof(cplx) * (ntau_ + 1) * element_size_);
767  memcpy(x + (tstp + 1 + ntau_ + 1) * element_size_, lesptr(0, tstp),
768  sizeof(cplx) * (tstp + 1) * element_size_);
769  }
770 }
772 #define herm_pseudo_INCR_TSTP \
773  if (alpha == cplx(1.0, 0.0)) { \
774  for (i = 0; i < len; i++) \
775  x0[i] += x[i]; \
776  } else { \
777  for (i = 0; i < len; i++) \
778  x0[i] += alpha * x[i]; \
779  }
780 template <typename T>
782 void herm_pseudo<T>::incr_timestep(int tstp, herm_matrix_timestep<T> &timestep,
783  std::complex<T> alpha) {
784  int i, len;
785  cplx *x, *x0;
786  assert(tstp >= -1 && tstp <= nt_);
787  assert(timestep.tstp_ == tstp && timestep.ntau_ == ntau_ && timestep.size1_ == size1_);
788  if (tstp == -1) {
789  len = (ntau_ + 1) * element_size_;
790  x0 = matptr(0);
791  x = timestep.data_;
792  herm_pseudo_INCR_TSTP
793  } else {
794  len = (tstp + 1) * element_size_;
795  x0 = retptr(tstp, 0);
796  x = timestep.data_;
797  herm_pseudo_INCR_TSTP len = (ntau_ + 1) * element_size_;
798  x0 = tvptr(tstp, 0);
799  x = timestep.data_ + (tstp + 1) * element_size_;
800  herm_pseudo_INCR_TSTP len = (tstp + 1) * element_size_;
801  x0 = lesptr(0, tstp);
802  x = timestep.data_ + (tstp + 1 + ntau_ + 1) * element_size_;
803  herm_pseudo_INCR_TSTP
804  }
805 }
806 #undef herm_pseudo_INCR_TSTP
807 template <typename T> void herm_pseudo<T>::smul(int tstp, T weight) {
809  int m;
810  cplx *x0;
811  assert(tstp >= -1 && tstp <= nt_);
812  if (tstp == -1) {
813  x0 = matptr(0);
814  for (m = 0; m <= ntau_; m++) {
815  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, weight);
816  }
817  } else {
818  x0 = retptr(tstp, 0);
819  for (m = 0; m <= tstp; m++) {
820  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, weight);
821  }
822  x0 = tvptr(tstp, 0);
823  for (m = 0; m <= ntau_; m++) {
824  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, weight);
825  }
826  x0 = lesptr(0, tstp);
827  for (m = 0; m <= tstp; m++) {
828  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, weight);
829  }
830  }
831 }
833 template <typename T> void herm_pseudo<T>::smul(int tstp, std::complex<T> weight) {
834  int m;
835  cplx *x0;
836  assert(tstp >= -1 && tstp <= nt_);
837  if (tstp == -1) {
838  x0 = matptr(0);
839  for (m = 0; m <= ntau_; m++) {
840  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, weight);
841  }
842  } else {
843  x0 = retptr(tstp, 0);
844  for (m = 0; m <= tstp; m++) {
845  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, weight);
846  }
847  x0 = tvptr(tstp, 0);
848  for (m = 0; m <= ntau_; m++) {
849  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, weight);
850  }
851  x0 = lesptr(0, tstp);
852  for (m = 0; m <= tstp; m++) {
853  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, weight);
854  }
855  }
856 }
857 
859 // multiply timestep with function ... same as for herm_matrix
860 // G(t,t') ==> F(t)G(t,t') ... ft+t*element_size_ points to F(t)
862 template <typename T>
863 void herm_pseudo<T>::left_multiply(int tstp, std::complex<T> *f0, std::complex<T> *ft,
864  T weight) {
865  int m;
866  cplx *xtemp, *ftemp, *x0;
867  xtemp = new cplx[element_size_];
868  assert(tstp >= -1 && tstp <= nt_);
869  if (tstp == -1) {
870  x0 = matptr(0);
871  for (m = 0; m <= ntau_; m++) {
872  element_mult<T, LARGESIZE>(size1_, xtemp, f0, x0 + m * element_size_);
873  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
874  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
875  }
876  } else {
877  ftemp = ft + tstp * element_size_;
878  x0 = retptr(tstp, 0);
879  for (m = 0; m <= tstp; m++) {
880  element_mult<T, LARGESIZE>(size1_, xtemp, ftemp, x0 + m * element_size_);
881  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
882  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
883  }
884  x0 = tvptr(tstp, 0);
885  for (m = 0; m <= ntau_; m++) {
886  element_mult<T, LARGESIZE>(size1_, xtemp, ftemp, x0 + m * element_size_);
887  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
888  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
889  }
890  x0 = lesptr(0, tstp);
891  for (m = 0; m <= tstp; m++) {
892  element_mult<T, LARGESIZE>(size1_, xtemp, ft + m * element_size_,
893  x0 + m * element_size_);
894  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
895  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
896  }
897  }
898  delete[] xtemp;
899 }
901 template <typename T>
902 void herm_pseudo<T>::left_multiply(int tstp, function<T> &ft, T weight) {
903  assert(tstp >= -1 && tstp <= nt_ && ft.nt() >= tstp && ft.size1() == size1_ &&
904  ft.size2() == size2_);
905  if (tstp >= 0)
906  left_multiply(tstp, ft.ptr(-1), ft.ptr(0), weight);
907  else
908  left_multiply(tstp, ft.ptr(-1), 0, weight);
909 }
910 // G(t,t') ==> F(t)G(t,t') ... ft+t*element_size_ points to F(t)
912 template <typename T>
913 void herm_pseudo<T>::right_multiply(int tstp, std::complex<T> *f0, std::complex<T> *ft,
914  T weight) {
915  int m;
916  cplx *xtemp, *ftemp, *x0;
917  xtemp = new cplx[element_size_];
918  assert(tstp >= -1 && tstp <= nt_);
919  if (tstp == -1) {
920  x0 = matptr(0);
921  for (m = 0; m <= ntau_; m++) {
922  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_, f0);
923  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
924  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
925  }
926  } else {
927  x0 = retptr(tstp, 0);
928  for (m = 0; m <= tstp; m++) {
929  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
930  ft + m * element_size_);
931  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
932  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
933  }
934  x0 = tvptr(tstp, 0);
935  for (m = 0; m <= ntau_; m++) {
936  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_, f0);
937  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
938  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
939  }
940  ftemp = ft + tstp * element_size_;
941  x0 = lesptr(0, tstp);
942  for (m = 0; m <= tstp; m++) {
943  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_, ftemp);
944  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
945  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
946  }
947  }
948  delete[] xtemp;
949 }
951 template <typename T>
952 void herm_pseudo<T>::right_multiply(int tstp, function<T> &ft, T weight) {
953  assert(tstp >= -1 && tstp <= nt_ && ft.nt() >= tstp && ft.size1() == size1_ &&
954  ft.size2() == size2_);
955  if (tstp >= 0)
956  right_multiply(tstp, ft.ptr(-1), ft.ptr(0), weight);
957  else
958  right_multiply(tstp, ft.ptr(-1), 0, weight);
959 }
960 
961 #if 0
962 
963 template<typename T> void herm_pseudo<T>::get_timestep(int tstp,std::vector<cplx> &data){
964  int len=2*(tstp+1)+ntau_,len1,j;
965  cplx *x;
966  len1=len*element_size_;
967  assert(tstp>=-1 && tstp<=nt_);
968  if((int)data.size()<len1) data.resize(len1);
969  j=(tstp>=0 ? IDXLOWER(tstp) : IDXIMAG(0));
970  x=elementptr(j,j);
971  for(i=0;i<len1;i++) data[i]=x[i];
972 }
973 template<typename T,class Matrix> void herm_pseudo<T>::get_timestep(int tstp,std::vector<Matrix> &data){
974  int len=2*(tstp+1)+ntau_,j;
975  cplx *x;
976  assert(tstp>=-1 && tstp<=nt_);
977  if((int)data.size()<len) data.resize(len);
978  j=(tstp>=0 ? IDXLOWER(tstp) : IDXIMAG(0));
979  x=elementptr(j,j);
980  for(i=0;i<len;i++){
981  if(data[i].size1()!=size1_ || data[i].size2()!=size1_) data[i].resize(size1_,size1_);
982 
983  data[i]=x[i*element_size_];
984  }
985 }
986 template<typename T> void herm_pseudo<T>::set_timestep(int tstp,std::vector<cplx> &data){
987  int len=2*(tstp+1)+ntau_;
988  assert((int)data.size()>=len*element_size_);
989  if(tstp>=0){
990  j=IDXLOWER(tstp);
991  for(i=0;i<len;i++){
992  i1=nt_-tstp+i;
993  }
994  }else{
995  }
996 }
997 template<typename T> void herm_pseudo<T>::set_timestep_to_zero(int tstp);
998 template<typename T> void herm_pseudo<T>::smul_timestep(int tstp,cplx lambda);
999 template<typename T> void herm_pseudo<T>::increment_timestep(int tstp,std::vector<cplx> &data,cplx alpha=1.0);
1000 template<typename T> void herm_pseudo<T>::set_timestep(int tstp,const herm_pseudo<T> &g1);
1001 template<typename T> void herm_pseudo<T>::set_timestep_pseudo(int tstp,const matrix<T> &g1);
1002 template<class Matrix> void get_timestep(int tstp,std::vector<Matrix> &data);
1003 template<class Matrix> void set_timestep(int tstp,std::vector<Matrix> &data);
1004 
1005 /* #######################################################################################
1006 #
1007 # MANIPULATING TIMESTEPS
1008 #
1009 ########################################################################################*/
1010 template <class GG>
1011 void set_timestep_pseudo(int tstp,GG &G,std::complex<double> *timestep){
1012  int sg=G.element_size(),ntau=G.ntau(),m,l=0;
1013  std::complex<double> *x1;
1014  x1=timestep;
1015  if(tstp>=0){
1016  for(m=0;m<=tstp;m++){
1017  // ret(tstp,m) = gtr(tstp,m)
1018  // = -G(m-,tstp-)^*
1019  G.element_set(G.retptr(tstp,tstp-m),x1,G);
1020  G.element_smul(G.retptr(tstp,tstp-m),-1);
1021  G.element_conj(G.retptr(tstp,tstp-m));
1022  x1+=sg;l++;
1023  }
1024  for(m=0;m<=ntau;m++){
1025  // x1 points to Gvt(m,tstp) = -bosefermi *Gtv(tstp,ntau-m)^*
1026  G.element_set(G.tvptr(tstp,ntau-m),x1,G);
1027  G.element_conj(G.tvptr(tstp,ntau-m));
1028  G.element_smul(G.tvptr(tstp,ntau-m),-G.sig());
1029  x1+=sg;l++;
1030  }
1031  for(m=0;m<=tstp;m++){
1032  // x1 points to G(m+,tstp+)=les(m,tstp)
1033  G.element_set(G.lesptr(m,tstp),x1,G);
1034  x1+=sg;l++;
1035  }
1036  }else{
1037  for(m=0;m<=ntau;m++){
1038  // x1 points to Gmat(m,0) = i*Gmat(m)
1039  G.element_set(G.matptr(m),x1,G);
1040  G.element_smul(G.matptr(m),std::complex<double>(0,-1));
1041  x1+=sg;l++;
1042  }
1043  }
1044 }
1045 template <class GG>
1046 void get_timestep_pseudo(int tstp,GG &G,std::complex<double> *timestep){
1047 int sg=G.element_size(),ntau=G.ntau(),m,l=0;
1048  std::complex<double> *x1;
1049  x1=timestep;
1050  //std::cout << "get tsp " << tstp << std::endl;
1051  if(tstp>=0){
1052  for(m=0;m<=tstp;m++){
1053  // ret(tstp,tstp-m) = gtr(tstp,tstp-m)
1054  // = -G^gtr(tstp-m,tstp)^*
1055  G.element_set(x1,G.retptr(tstp,tstp-m),G);
1056  G.element_smul(x1,-1);
1057  G.element_conj(x1);
1058  x1+=sg;l++;
1059  }
1060  for(m=0;m<=ntau;m++){
1061  // x1 points to Gvt(m,tstp) = -bosefermi *Gtv(tstp,ntau-m)^*
1062  G.element_set(x1,G.tvptr(tstp,ntau-m),G);
1063  G.element_conj(x1);
1064  G.element_smul(x1,-G.sig());
1065  x1+=sg;l++;
1066  }
1067  for(m=0;m<=tstp;m++){
1068  // x1 points to G(m+,tstp+)=les(m,tstp)
1069  G.element_set(x1,G.lesptr(m,tstp),G);
1070  x1+=sg;l++;
1071  }
1072  }else{
1073  for(m=0;m<=ntau;m++){
1074  // x1 points to Gmat(m,0) = i*Gmat(m)
1075  G.element_set(x1,G.matptr(m),G);
1076  G.element_smul(x1,std::complex<double>(0,1));
1077  x1+=sg;l++;
1078  }
1079  }
1080 }
1081 template <typename T,class Matrix>
1082 void set_timestep_matrix_pseudo(int tstp,matrix<T> &G,std::vector<Matrix> &data){
1083  int i,l,m,sg,dim,ntau=G.ntau(),nmax=(tstp<0 ? ntau+1 : ntau+2*tstp+3);
1084  std::complex<T> *data1;
1085  assert((int)data.size()>=nmax);
1086  dim=G.size1();
1087  sg=dim*dim;
1088  data1=new std::complex<T> [nmax*sg];
1089  for(i=0;i<nmax;i++) for(l=0;l<dim;l++) for(m=0;m<dim;m++) data1[i*sg+l*dim+m]=data[i](l,m);
1090  set_timestep_pseudo(tstp,G,data1);
1091  delete [] data1;
1092 }
1093 template <typename T,class Matrix>
1094 void increment_timestep_matrix_pseudo(int tstp,matrix<T> &G,std::complex<T> alpha,std::vector<Matrix> &data){
1095  int i,l,m,sg,dim,ntau=G.ntau(),nmax=(tstp<0 ? ntau+1 : ntau+2*tstp+3);
1096  std::complex<T> *data1;
1097  assert((int)data.size()>=nmax);
1098  dim=G.size1();
1099  sg=dim*dim;
1100  data1=new std::complex<T> [nmax*sg];
1101  get_timestep_pseudo(tstp,G,data1);
1102  if(alpha.real()==1 && alpha.imag()==0) for(i=0;i<nmax;i++) for(l=0;l<dim;l++) for(m=0;m<dim;m++) data1[i*sg+l*dim+m] += data[i](l,m);
1103  else for(i=0;i<nmax;i++) for(l=0;l<dim;l++) for(m=0;m<dim;m++) data1[i*sg+l*dim+m] += alpha*data[i](l,m);
1104  set_timestep_pseudo(tstp,G,data1);
1105  delete [] data1;
1106 }
1107 template <typename T,class Matrix>
1108 void get_timestep_matrix_pseudo(int tstp,matrix<T> &G,std::vector<Matrix> &data){
1109  int i,l,m,sg,dim,ntau=G.ntau(),nmax=(tstp<0 ? ntau+1 : ntau+2*tstp+3);
1110  std::complex<T> *data1;
1111  assert((int)data.size()>=nmax);
1112  dim=G.size1();
1113  sg=dim*dim;
1114  data1=new std::complex<T> [nmax*sg];
1115  get_timestep_pseudo(tstp,G,data1);
1116  for(i=0;i<nmax;i++) for(l=0;l<dim;l++) for(m=0;m<dim;m++) data[i](l,m)=data1[i*sg+l*dim+m];
1117  delete [] data1;
1118 }
1119 
1120 #endif // 0
1121 
1122 } // namespace cntr
1123 
1124 #endif // CNTR_HERM_PSEUDO_IMPL_H
cntr::get_mat
void get_mat(const int m, std::complex< T > &G_mat, herm_matrix< T > &G, herm_matrix< T > &Gcc)
cntr::get_timestep
void get_timestep(int tstp, std::vector< Matrix > &data)
cntr::set_timestep_pseudo
void set_timestep_pseudo(int tstp, GG &G, std::complex< double > *timestep)
Definition: cntr_herm_pseudo_impl.hpp:1011
cntr_herm_matrix_timestep_decl.hpp
cntr_herm_pseudo_decl.hpp
cntr_elements.hpp
cntr_function_decl.hpp
cntr::get_tv
void get_tv(const int i, const int m, std::complex< T > &G_tv, herm_matrix< T > &G, herm_matrix< T > &Gcc)
fourier::cplx
std::complex< double > cplx
Definition: fourier.cpp:11
cntr::get_vt
void get_vt(const int m, const int i, std::complex< T > &G_vt, herm_matrix< T > &G, herm_matrix< T > &Gcc)
cntr::increment_timestep_matrix_pseudo
void increment_timestep_matrix_pseudo(int tstp, matrix< T > &G, std::complex< T > alpha, std::vector< Matrix > &data)
Definition: cntr_herm_pseudo_impl.hpp:1094
cntr::get_ret
void get_ret(const int i, const int j, std::complex< T > &G_ret, herm_matrix< T > &G, herm_matrix< T > &Gcc)
cntr::set_timestep_matrix_pseudo
void set_timestep_matrix_pseudo(int tstp, matrix< T > &G, std::vector< Matrix > &data)
Definition: cntr_herm_pseudo_impl.hpp:1082
cntr::get_timestep_matrix_pseudo
void get_timestep_matrix_pseudo(int tstp, matrix< T > &G, std::vector< Matrix > &data)
Definition: cntr_herm_pseudo_impl.hpp:1108
cntr::get_les
void get_les(const int i, const int j, std::complex< T > &G_les, herm_matrix< T > &G, herm_matrix< T > &Gcc)
cntr::get_gtr
void get_gtr(const int i, const int j, std::complex< T > &G_gtr, herm_matrix< T > &G, herm_matrix< T > &Gcc)
cntr::set_timestep
void set_timestep(int tstp, std::vector< Matrix > &data)
cntr::get_timestep_pseudo
void get_timestep_pseudo(int tstp, GG &G, std::complex< double > *timestep)
Definition: cntr_herm_pseudo_impl.hpp:1046