NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_herm_matrix_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_HERM_MATRIX_IMPL_H
2 #define CNTR_HERM_MATRIX_IMPL_H
3 
5 //#include "cntr_exception.hpp"
6 #include "cntr_elements.hpp"
7 #include "cntr_function_decl.hpp"
10 
11 namespace cntr {
12 
13 /* #######################################################################################
14 #
15 # CONSTRUCTION/DESTRUCTION
16 #
17 ########################################################################################*/
18 template <typename T>
20  les_ = 0;
21  tv_ = 0;
22  ret_ = 0;
23  mat_ = 0;
24  ntau_ = 0;
25  nt_ = 0;
26  size1_ = 0;
27  size2_ = 0;
28  element_size_ = 0;
29  sig_ = -1;
30 }
31 template <typename T>
33  delete[] les_;
34  delete[] ret_;
35  delete[] tv_;
36  delete[] mat_;
37 }
38 
61 template <typename T>
62 herm_matrix<T>::herm_matrix(int nt, int ntau, int size1, int sig) {
63  assert(size1 >= 0 && nt >= -1 && sig * sig == 1 && ntau >= 0);
64  nt_ = nt;
65  ntau_ = ntau;
66  sig_ = sig;
67  size1_ = size1;
68  size2_ = size1;
69  element_size_ = size1 * size1;
70  if (size1 > 0) {
71  mat_ = new cplx[(ntau_ + 1) * element_size_];
72  } else {
73  mat_ = 0;
74  }
75  if (nt >= 0 && size1 > 0) {
76  les_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
77  ret_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
78  tv_ = new cplx[(nt_ + 1) * (ntau_ + 1) * element_size_];
79  } else {
80  les_ = 0;
81  tv_ = 0;
82  ret_ = 0;
83  }
84 }
85 
111 template <typename T> herm_matrix<T>::herm_matrix(int nt,int ntau,int size1,int size2,int sig){
112  assert(size1>=0 && size2>=0 && nt>=-1 && sig*sig==1 && ntau>=0);
113  nt_=nt;
114  ntau_=ntau;
115  sig_=sig;
116  size1_=size1;
117  size2_=size2;
118  element_size_=size1*size2;
119  if(size1>0){
120  mat_ = new cplx [(ntau_+1)*element_size_];
121  memset(mat_, 0, sizeof(cplx)*(ntau_+1)*element_size_);
122  }else{
123  mat_=0;
124  }
125  if(nt>=0 && size1>0){
126  les_ = new cplx [((nt_+1)*(nt_+2))/2*element_size_];
127  ret_ = new cplx [((nt_+1)*(nt_+2))/2*element_size_];
128  tv_ = new cplx [(nt_+1)*(ntau_+1)*element_size_];
129  memset(les_, 0, sizeof(cplx)*((nt_+1)*(nt_+2))/2*element_size_);
130  memset(ret_, 0, sizeof(cplx)*((nt_+1)*(nt_+2))/2*element_size_);
131  memset(tv_, 0, sizeof(cplx)*(nt_+1)*(ntau_+1)*element_size_);
132  }else{
133  les_=0;
134  tv_=0;
135  ret_=0;
136  }
137 }
138 
156 template <typename T>
158  nt_ = g.nt_;
159  ntau_ = g.ntau_;
160  sig_ = g.sig_;
161  size1_ = g.size1_;
162  size2_ = g.size1_;
163  element_size_ = size1_ * size1_;
164  if (size1_ > 0) {
165  mat_ = new cplx[(ntau_ + 1) * element_size_];
166  memcpy(mat_, g.mat_, sizeof(cplx) * (ntau_ + 1) * element_size_);
167  } else {
168  mat_ = 0;
169  }
170  if (nt_ >= 0 && size1_ > 0) {
171  les_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
172  ret_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
173  tv_ = new cplx[(nt_ + 1) * (ntau_ + 1) * element_size_];
174  memcpy(les_, g.les_,
175  sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
176  memcpy(ret_, g.ret_,
177  sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
178  memcpy(tv_, g.tv_,
179  sizeof(cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
180  } else {
181  les_ = 0;
182  ret_ = 0;
183  tv_ = 0;
184  }
185 }
186 template <typename T>
188  if (this == &g)
189  return *this;
190  sig_ = g.sig_;
191  if (nt_ != g.nt_ || ntau_ != g.ntau_ || size1_ != g.size1_) {
192  delete[] les_;
193  delete[] ret_;
194  delete[] tv_;
195  delete[] mat_;
196  nt_ = g.nt_;
197  ntau_ = g.ntau_;
198  size1_ = g.size1_;
199  size2_ = g.size1_;
200  element_size_ = size1_ * size1_;
201  if (size1_ > 0) {
202  mat_ = new cplx[(ntau_ + 1) * element_size_];
203  } else {
204  mat_ = 0;
205  }
206  if (size1_ > 0 && nt_ >= 0) {
207  les_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
208  ret_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
209  tv_ = new cplx[(nt_ + 1) * (ntau_ + 1) * element_size_];
210  } else {
211  les_ = 0;
212  ret_ = 0;
213  tv_ = 0;
214  }
215  }
216  if (size1_ > 0) {
217  memcpy(mat_, g.mat_, sizeof(cplx) * (ntau_ + 1) * element_size_);
218  if (nt_ >= 0) {
219  memcpy(les_, g.les_, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 *
220  element_size_);
221  memcpy(ret_, g.ret_, sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 *
222  element_size_);
223  memcpy(tv_, g.tv_,
224  sizeof(cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
225  }
226  }
227  return *this;
228 }
229 #if __cplusplus >= 201103L
230 template <typename T>
232  : les_(g.les_),
233  ret_(g.ret_),
234  tv_(g.tv_),
235  mat_(g.mat_),
236  nt_(g.nt_),
237  ntau_(g.ntau_),
238  size1_(g.size1_),
239  size2_(g.size2_),
240  element_size_(g.element_size_),
241  sig_(g.sig_) {
242  g.les_ = nullptr;
243  g.tv_ = nullptr;
244  g.ret_ = nullptr;
245  g.mat_ = nullptr;
246  g.ntau_ = 0;
247  g.nt_ = 0;
248  g.size1_ = 0;
249  g.size2_ = 0;
250  g.element_size_ = 0;
251 }
252 template <typename T>
254  if (&g == this)
255  return *this;
256 
257  les_ = g.les_;
258  ret_ = g.ret_;
259  tv_ = g.tv_;
260  mat_ = g.mat_;
261  nt_ = g.nt_;
262  ntau_ = g.ntau_;
263  size1_ = g.size1_;
264  size2_ = g.size2_;
265  element_size_ = g.element_size_;
266  sig_ = g.sig_;
267 
268  g.les_ = nullptr;
269  g.tv_ = nullptr;
270  g.ret_ = nullptr;
271  g.mat_ = nullptr;
272  g.ntau_ = 0;
273  g.nt_ = 0;
274  g.size1_ = 0;
275  g.size2_ = 0;
276  g.element_size_ = 0;
277 
278  return *this;
279 }
280 #endif
281 /* #######################################################################################
282 #
283 # RESIZE
284 #
285 ########################################################################################*/
286 
312 template <typename T>
313 void herm_matrix<T>::resize_discard(int nt, int ntau, int size1) {
314  assert(ntau >= 0 && nt >= -1 && size1 >= 0);
315  delete[] les_;
316  delete[] ret_;
317  delete[] tv_;
318  delete[] mat_;
319  nt_ = nt;
320  ntau_ = ntau;
321  size1_ = size1;
322  size2_ = size1;
323  element_size_ = size1 * size1;
324  if (size1 > 0) {
325  mat_ = new cplx[(ntau_ + 1) * element_size_];
326  } else {
327  mat_ = 0;
328  }
329  if (nt_ >= 0 && size1_ > 0) {
330  les_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
331  ret_ = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
332  tv_ = new cplx[(nt_ + 1) * (ntau_ + 1) * element_size_];
333  } else {
334  les_ = 0;
335  tv_ = 0;
336  ret_ = 0;
337  }
338 }
339 
360 template <typename T>
362  int nt1 = (nt_ > nt ? nt : nt_);
363  cplx *ret, *les, *tv;
364  assert(nt >= -1);
365  nt_ = nt;
366  if (size1_ == 0)
367  return;
368  if (nt_ >= 0) {
369  les = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
370  ret = new cplx[((nt_ + 1) * (nt_ + 2)) / 2 * element_size_];
371  tv = new cplx[(nt_ + 1) * (ntau_ + 1) * element_size_];
372  if (nt1 >= 0) {
373  memcpy(les, les_, sizeof(cplx) * ((nt1 + 1) * (nt1 + 2)) / 2 *
374  element_size_);
375  memcpy(ret, ret_, sizeof(cplx) * ((nt1 + 1) * (nt1 + 2)) / 2 *
376  element_size_);
377  memcpy(tv, tv_,
378  sizeof(cplx) * (nt1 + 1) * (ntau_ + 1) * element_size_);
379  }
380  } else {
381  les = 0;
382  ret = 0;
383  tv = 0;
384  }
385  delete[] les_;
386  delete[] ret_;
387  delete[] tv_;
388  les_ = les;
389  ret_ = ret;
390  tv_ = tv;
391 }
416 template <typename T>
417 void herm_matrix<T>::resize(int nt, int ntau, int size1) {
418  // std::cout << "herm_matrix<T>::" << __FUNCTION__ << " " << nt << " " <<
419  // ntau << " " << size1 << std::endl;
420  if (ntau == ntau_ && size1_ == size1)
421  resize_nt(nt);
422  else
423  resize_discard(nt, ntau, size1);
424 }
436 template <typename T>
438  if (size1_ == 0)
439  return;
440  memset(mat_, 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
441  if (nt_ >= 0) {
442  memset(les_, 0,
443  sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
444  memset(ret_, 0,
445  sizeof(cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
446  memset(tv_, 0,
447  sizeof(cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
448  }
449 }
450 /* #######################################################################################
451 #
452 # RAW POINTERS TO ELEMENTS
453 #
454 ########################################################################################*/
456 template <typename T>
457 int herm_matrix<T>::les_offset(int t, int t1) const {
458  assert(t >= 0 && t1 >= 0 && t <= t1 && t1 <= nt_);
459  return ((t1 * (t1 + 1)) / 2 + t) * element_size_;
460 }
462 template <typename T>
463 int herm_matrix<T>::ret_offset(int t, int t1) const {
464  assert(t >= 0 && t1 >= 0 && t <= nt_ && t1 <= t);
465  return ((t * (t + 1)) / 2 + t1) * element_size_;
466 }
468 template <typename T>
469 int herm_matrix<T>::tv_offset(int t, int tau) const {
470  assert(t >= 0 && tau >= 0 && t <= nt_ && tau <= ntau_);
471  return (t * (ntau_ + 1) + tau) * element_size_;
472 }
474 template <typename T>
475 int herm_matrix<T>::mat_offset(int tau) const {
476  assert(tau >= 0 && tau <= ntau_);
477  return tau * element_size_;
478 }
480 template <typename T>
481 inline std::complex<T> *herm_matrix<T>::lesptr(int t, int t1) {
482  return les_ + les_offset(t, t1);
483 }
485 template <typename T>
486 inline std::complex<T> *herm_matrix<T>::retptr(int t, int t1) {
487  return ret_ + ret_offset(t, t1);
488 }
490 template <typename T>
491 inline std::complex<T> *herm_matrix<T>::tvptr(int t, int tau) {
492  return tv_ + tv_offset(t, tau);
493 }
495 template <typename T>
496 inline std::complex<T> *herm_matrix<T>::matptr(int tau) {
497  return mat_ + mat_offset(tau);
498 }
500 template <typename T>
501 inline const std::complex<T> *herm_matrix<T>::lesptr(int t, int t1) const {
502  return les_ + les_offset(t, t1);
503 }
505 template <typename T>
506 inline const std::complex<T> *herm_matrix<T>::retptr(int t, int t1) const {
507  return ret_ + ret_offset(t, t1);
508 }
510 template <typename T>
511 inline const std::complex<T> *herm_matrix<T>::tvptr(int t, int tau) const {
512  return tv_ + tv_offset(t, tau);
513 }
515 template <typename T>
516 inline const std::complex<T> *herm_matrix<T>::matptr(int tau) const {
517  return mat_ + mat_offset(tau);
518 }
519 /* #######################################################################################
520 #
521 # READING ELEMENTS TO ANY MATRIX TYPE
522 # OR TO COMPLEX NUMBERS (then only the (0,0) element is addressed for dim>0)
523 # note: these are not efficient, in particular the conjugation
524 #
525 ########################################################################################*/
526 #define herm_matrix_READ_ELEMENT \
527  { \
528  int r, s; \
529  M.resize(size1_, size2_); \
530  for (r = 0; r < size1_; r++) \
531  for (s = 0; s < size2_; s++) \
532  M(r, s) = x[r * size2_ + s]; \
533  }
534 // TODO: Please think about this function, what happens for rectangular
535 #define herm_matrix_READ_ELEMENT_MINUS_CONJ \
536  { \
537  cplx w; \
538  int r, s, dim = size1_; \
539  M.resize(dim, dim); \
540  for (r = 0; r < dim; r++) \
541  for (s = 0; s < dim; s++) { \
542  w = x[s * dim + r]; \
543  M(r, s) = std::complex<T>(-w.real(), w.imag()); \
544  } \
545  }
546 
569  template <typename T> template <class Matrix>
570  void herm_matrix<T>::get_les(int i, int j, Matrix &M) const {
571  assert(i <= nt_ && j <= nt_);
572  const cplx *x;
573  if (i <= j) {
574  x = lesptr(i, j);
575  herm_matrix_READ_ELEMENT
576  } else {
577  x = lesptr(j, i);
578  herm_matrix_READ_ELEMENT_MINUS_CONJ
579  }
580  }
582 
605 template <typename T>
606 template <class Matrix>
607 void herm_matrix<T>::get_ret(int i, int j, Matrix &M) const {
608  assert(i <= nt_ && j <= nt_);
609  const cplx *x;
610  if (i >= j) {
611  x = retptr(i, j);
612  herm_matrix_READ_ELEMENT
613  } else {
614  x = retptr(j, i);
615  herm_matrix_READ_ELEMENT_MINUS_CONJ
616  }
617 }
619 
639 template <typename T>
640 template <class Matrix>
641 void herm_matrix<T>::get_tv(int i, int j, Matrix &M) const {
642  assert(i <= nt_ && j <= ntau_);
643  const cplx *x = tvptr(i, j);
644  herm_matrix_READ_ELEMENT
645 }
647 
668 template <typename T>
669 template <class Matrix>
670 void herm_matrix<T>::get_vt(int i, int j, Matrix &M) const {
671  assert(i <= ntau_ && j <= nt_);
672  const cplx *x = tvptr(j, ntau_ - i);
673  herm_matrix_READ_ELEMENT_MINUS_CONJ if (sig_ == -1) M = -M;
674 }
676 
694 template <typename T>
695 template <class Matrix>
696 void herm_matrix<T>::get_mat(int i, Matrix &M) const {
697  assert(i <= ntau_);
698  const cplx *x = matptr(i);
699  herm_matrix_READ_ELEMENT
700 }
702 
721 template <typename T>
722 template <class Matrix>
723 void herm_matrix<T>::get_matminus(int i, Matrix &M) const {
724  assert(i <= ntau_);
725  const cplx *x = matptr(ntau_ - i);
726  herm_matrix_READ_ELEMENT if (sig_ == -1) M = -M;
727 }
728 #undef herm_matrix_READ_ELEMENT
729 #undef herm_matrix_READ_ELEMENT_MINUS_CONJ
730 
732 
752 template <typename T>
753 template <class Matrix>
754 void herm_matrix<T>::get_gtr(int i, int j, Matrix &M) const {
755  assert(i <= nt_ && j <= nt_);
756  Matrix M1;
757  get_ret(i, j, M);
758  get_les(i, j, M1);
759  M += M1;
760 }
762 
785 template <typename T>
786 inline void herm_matrix<T>::get_ret(int i, int j, cplx &x) const {
787  assert(i <= nt_ && j <= nt_);
788  if (i >= j)
789  x = *retptr(i, j);
790  else {
791  x = *retptr(j, i);
792  x = -std::conj(x);
793  }
794 }
796 
817 template <typename T>
818 inline void herm_matrix<T>::get_les(int i, int j, cplx &x) const {
819  assert(i <= nt_ && j <= nt_);
820  if (i <= j)
821  x = *lesptr(i, j);
822  else {
823  x = *lesptr(j, i);
824  x = -std::conj(x);
825  }
826 }
828 
848 template <typename T>
849 inline void herm_matrix<T>::get_tv(int i, int j, cplx &x) const {
850  assert(i <= nt_ && j <= ntau_);
851  x = *tvptr(i, j);
852 }
854 
875 template <typename T>
876 inline void herm_matrix<T>::get_vt(int i, int j, cplx &x) const {
877  assert(i <= ntau_ && j <= nt_);
878  x = *tvptr(j, ntau_ - i);
879  if (sig_ == -1)
880  x = std::conj(x);
881  else
882  x = -std::conj(x);
883 }
885 
903 template <typename T>
904 inline void herm_matrix<T>::get_mat(int i, cplx &x) const {
905  assert(i <= ntau_);
906  x = *matptr(i);
907 }
909 
928 template <typename T>
929 inline void herm_matrix<T>::get_matminus(int i, cplx &x) const {
930  assert(i <= ntau_);
931  x = *matptr(ntau_ - i);
932  if (sig_ == -1)
933  x = -x;
934 }
936 
956 template <typename T>
957 inline void herm_matrix<T>::get_gtr(int i, int j, cplx &x) const {
958  assert(i <= nt_ && j <= nt_);
959  cplx x1;
960  get_ret(i, j, x);
961  get_les(i, j, x1);
962  x += x1;
963 }
965 
983 template <typename T>
984 std::complex<T> herm_matrix<T>::density_matrix(int tstp) const {
985  assert(tstp >= -1 && tstp <= nt_);
986  cplx x1;
987  if (tstp == -1) {
988  get_mat(ntau_, x1);
989  return -x1;
990  } else {
991  get_les(tstp, tstp, x1);
992  return std::complex<T>(0.0, sig_) * x1;
993  }
994 }
995 
1016 template <typename T>
1017 template <class Matrix>
1018 void herm_matrix<T>::density_matrix(int tstp, Matrix &M) const {
1019  assert(M.rows() == size1_ && M.cols() == size2_);
1020  assert(tstp >= -1 && tstp <= nt_);
1021  if (tstp == -1) {
1022  get_mat(ntau_, M);
1023  M *= (-1.0);
1024  } else {
1025  get_les(tstp, tstp, M);
1026  M *= std::complex<T>(0.0, 1.0 * sig_);
1027  }
1028 }
1029 
1030 /* #######################################################################################
1031 #
1032 # WRITING ELEMENTS FROM ANY MATRIX TYPE
1033 # OR FROM COMPLEX NUMBERS (then only the (0,0) element is addressed for
1034 dim>0)
1035 #
1036 ########################################################################################*/
1037 #define herm_matrix_SET_ELEMENT_MATRIX \
1038  { \
1039  int r, s; \
1040  for (r = 0; r < size1_; r++) \
1041  for (s = 0; s < size2_; s++) \
1042  x[r * size2_ + s] = M(r, s); \
1043  }
1044 
1046 
1063 template <typename T>
1064 template <class Matrix>
1065 void herm_matrix<T>::set_ret(int i, int j, Matrix &M) {
1066  assert(i <= nt_ && j <= nt_);
1067  cplx *x = retptr(i, j);
1068  herm_matrix_SET_ELEMENT_MATRIX
1069 }
1071 
1089 template <typename T>
1090 template <class Matrix>
1091 void herm_matrix<T>::set_les(int i, int j, Matrix &M) {
1092  assert(i <= nt_ && j <= nt_);
1093  cplx *x = lesptr(i, j);
1094  herm_matrix_SET_ELEMENT_MATRIX
1095 }
1097 
1115 template <typename T>
1116 template <class Matrix>
1117 void herm_matrix<T>::set_tv(int i, int j, Matrix &M) {
1118  assert(i <= nt_ && j <= ntau_);
1119  cplx *x = tvptr(i, j);
1120  herm_matrix_SET_ELEMENT_MATRIX
1121 }
1123 
1138 template <typename T>
1139 template <class Matrix>
1140 void herm_matrix<T>::set_mat(int i, Matrix &M) {
1141  assert(i <= ntau_);
1142  cplx *x = matptr(i);
1143  herm_matrix_SET_ELEMENT_MATRIX
1144 }
1145 
1146 
1161 template<typename T>
1162 void herm_matrix<T>::set_mat_herm(int i){
1163  cdmatrix tmp(size1_,size2_);
1164  this->get_mat(i,tmp);
1165  tmp=(tmp+tmp.adjoint())/2.0;
1166  this->set_mat(i,tmp);
1167 }
1169 
1179 template<typename T>
1180 void herm_matrix<T>::set_mat_herm(void){
1181  cdmatrix tmp(size1_,size2_);
1182  for(int i=0;i<ntau_;i++){
1183  this->get_mat(i,tmp);
1184  tmp=(tmp+tmp.adjoint())/2.0;
1185  this->set_mat(i,tmp);
1186  }
1187 }
1189 
1206 template <typename T>
1207 inline void herm_matrix<T>::set_les(int i, int j, cplx x) {
1208  assert(i <= nt_ && j <= nt_);
1209  *lesptr(i, j) = x;
1210 }
1211 
1213 
1230 template <typename T>
1231 inline void herm_matrix<T>::set_ret(int i, int j, cplx x) {
1232  assert(i <= nt_ && j <= nt_);
1233  *retptr(i, j) = x;
1234 }
1236 
1253 template <typename T>
1254 inline void herm_matrix<T>::set_tv(int i, int j, cplx x) {
1255  assert(i <= nt_ && j <= ntau_);
1256  *tvptr(i, j) = x;
1257 }
1259 
1274 template <typename T>
1275 inline void herm_matrix<T>::set_mat(int i, cplx x) {
1276  assert(i <= ntau_ );
1277  *matptr(i) = x;
1278 }
1279 
1280 
1281 
1282 /* #######################################################################################
1283 #
1284 # INPUT/OUTPUT FROM/TO FILES
1285 #
1286 ########################################################################################*/
1287 
1307 template <typename T>
1308 void herm_matrix<T>::print_to_file(const char *file, int precision) const {
1309  int i, j, l, sg = element_size_;
1310  std::ofstream out;
1311  out.open(file, std::ios::out);
1312  out.precision(precision);
1313  out << "# " << nt_ << " " << ntau_ << " " << size1_ << " "
1314  << " " << sig_ << std::endl;
1315  for (j = 0; j <= ntau_; j++) {
1316  out << "mat: " << j;
1317  for (l = 0; l < sg; l++)
1318  out << " " << matptr(j)[l].real() << " " << matptr(j)[l].imag();
1319  out << std::endl;
1320  }
1321  out << std::endl;
1322  if (nt_ >= 0) {
1323  for (i = 0; i <= nt_; i++) {
1324  for (j = 0; j <= i; j++) {
1325  out << "ret: " << i << " " << j;
1326  for (l = 0; l < sg; l++)
1327  out << " " << retptr(i, j)[l].real() << " "
1328  << retptr(i, j)[l].imag();
1329  out << std::endl;
1330  }
1331  out << std::endl;
1332  }
1333  out << std::endl;
1334  for (i = 0; i <= nt_; i++) {
1335  for (j = 0; j <= ntau_; j++) {
1336  out << "tv: " << i << " " << j;
1337  for (l = 0; l < sg; l++)
1338  out << " " << tvptr(i, j)[l].real() << " "
1339  << tvptr(i, j)[l].imag();
1340  out << std::endl;
1341  }
1342  out << std::endl;
1343  }
1344  out << std::endl;
1345  for (j = 0; j <= nt_; j++) {
1346  for (i = 0; i <= j; i++) {
1347  out << "les: " << i << " " << j;
1348  for (l = 0; l < sg; l++)
1349  out << " " << lesptr(i, j)[l].real() << " "
1350  << lesptr(i, j)[l].imag();
1351  out << std::endl;
1352  }
1353  out << std::endl;
1354  }
1355  out << std::endl;
1356  }
1357  out.close();
1358 }
1375 template <typename T>
1376 void herm_matrix<T>::read_from_file(const char *file) {
1377  int i, n, m, j, l, size1, sg, sig;
1378  double real, imag;
1379  std::string s;
1380  std::ifstream out;
1381  out.open(file, std::ios::in);
1382  if (!(out >> s >> n >> m >> size1 >> sig)) {
1383  std::cerr << "read G from file " << file << " error in file"
1384  << std::endl;
1385  abort();
1386  }
1387  if (n > nt_ || m != ntau_ || size1 != size1_)
1388  resize(n, m, size1);
1389  sig_ = sig;
1390  sg = element_size_;
1391  for (j = 0; j <= ntau_; j++) {
1392  out >> s >> s;
1393  for (l = 0; l < sg; l++) {
1394  if (!(out >> real >> imag)) {
1395  std::cerr << "read G from file " << file << " error at mat ("
1396  << j << ")" << std::endl;
1397  abort();
1398  }
1399  matptr(j)[l] = std::complex<T>(real, imag);
1400  }
1401  }
1402  if (n >= 0) {
1403  for (i = 0; i <= n; i++) {
1404  for (j = 0; j <= i; j++) {
1405  out >> s >> s >> s;
1406  for (l = 0; l < sg; l++) {
1407  if (!(out >> real >> imag)) {
1408  std::cerr << "read G from file " << file
1409  << " error at ret (" << i << "," << j << ")"
1410  << std::endl;
1411  abort();
1412  }
1413  retptr(i, j)[l] = std::complex<T>(real, imag);
1414  }
1415  }
1416  }
1417  for (i = 0; i <= n; i++) {
1418  for (j = 0; j <= ntau_; j++) {
1419  out >> s >> s >> s;
1420  for (l = 0; l < sg; l++) {
1421  if (!(out >> real >> imag)) {
1422  std::cerr << "read G from file " << file
1423  << " error at tv (" << i << "," << j << ")"
1424  << std::endl;
1425  abort();
1426  }
1427  tvptr(i, j)[l] = std::complex<T>(real, imag);
1428  }
1429  }
1430  }
1431  for (j = 0; j <= n; j++) {
1432  for (i = 0; i <= j; i++) {
1433  out >> s >> s >> s;
1434  for (l = 0; l < sg; l++) {
1435  if (!(out >> real >> imag)) {
1436  std::cerr << "read G from file " << file
1437  << " error at les (" << i << "," << j << ")"
1438  << std::endl;
1439  abort();
1440  }
1441  lesptr(i, j)[l] = std::complex<T>(real, imag);
1442  }
1443  }
1444  }
1445  }
1446  out.close();
1447 }
1468 template <typename T>
1469 void herm_matrix<T>::read_from_file(int nt1, const char *file) {
1470  int i, n, m, j, l, size1, sg, sig;
1471  double real, imag;
1472  std::string s;
1473  std::ifstream out;
1474  out.open(file, std::ios::in);
1475 
1476  if (!(out >> s >> n >> m >> size1 >> sig)) {
1477  std::cerr << "read G from file " << file << " error in file"
1478  << std::endl;
1479  abort();
1480  }
1481 
1482  assert(nt1 <= nt_ && "nt1 > nt_");
1483  assert(nt1 <= n && "nt1 > n");
1484  assert(m == ntau_ && "m /= ntau_");
1485  assert(size1 == size1_ && "size1 /= size1_");
1486 
1487  sg = size1 * size1;
1488  if (nt1 >= -1) {
1489  for (j = 0; j <= ntau_; j++) {
1490  out >> s >> s;
1491  for (l = 0; l < sg; l++) {
1492  if (!(out >> real >> imag)) {
1493  std::cerr << "read G from file " << file
1494  << " error at mat (" << j << ")" << std::endl;
1495  abort();
1496  }
1497  matptr(j)[l] = std::complex<T>(real, imag);
1498  }
1499  }
1500  }
1501  if (n >= 0) {
1502  for (i = 0; i <= n; i++) {
1503  for (j = 0; j <= i; j++) {
1504  out >> s >> s >> s;
1505  for (l = 0; l < sg; l++) {
1506  if (!(out >> real >> imag)) {
1507  std::cerr << "read G from file " << file
1508  << " error at ret (" << i << "," << j << ")"
1509  << std::endl;
1510  abort();
1511  }
1512  if (i <= nt1)
1513  retptr(i, j)[l] = std::complex<T>(real, imag);
1514  }
1515  }
1516  }
1517  for (i = 0; i <= n; i++) {
1518  for (j = 0; j <= ntau_; j++) {
1519  out >> s >> s >> s;
1520  for (l = 0; l < sg; l++) {
1521  if (!(out >> real >> imag)) {
1522  std::cerr << "read G from file " << file
1523  << " error at tv (" << i << "," << j << ")"
1524  << std::endl;
1525  abort();
1526  }
1527  if (i <= nt1)
1528  tvptr(i, j)[l] = std::complex<T>(real, imag);
1529  }
1530  }
1531  }
1532  for (j = 0; j <= n; j++) {
1533  for (i = 0; i <= j; i++) {
1534  out >> s >> s >> s;
1535  for (l = 0; l < sg; l++) {
1536  if (!(out >> real >> imag)) {
1537  std::cerr << "read G from file " << file
1538  << " error at les (" << i << "," << j << ")"
1539  << std::endl;
1540  abort();
1541  }
1542  if (j <= nt1)
1543  lesptr(i, j)[l] = std::complex<T>(real, imag);
1544  }
1545  }
1546  }
1547  }
1548  out.close();
1549 }
1550 
1551 #if CNTR_USE_HDF5 == 1
1552 
1570 template <typename T>
1571 void herm_matrix<T>::write_to_hdf5(hid_t group_id) {
1572  store_int_attribute_to_hid(group_id, std::string("ntau"), ntau_);
1573  store_int_attribute_to_hid(group_id, std::string("nt"), nt_);
1574  store_int_attribute_to_hid(group_id, std::string("sig"), sig_);
1575  store_int_attribute_to_hid(group_id, std::string("size1"), size1_);
1576  store_int_attribute_to_hid(group_id, std::string("size2"), size2_);
1577  store_int_attribute_to_hid(group_id, std::string("element_size"),
1578  element_size_);
1579  hsize_t len_shape = 3, shape[3];
1580  shape[1] = size1_;
1581  shape[2] = size2_;
1582  if (nt_ > -2) {
1583  shape[0] = ntau_ + 1;
1584  store_cplx_array_to_hid(group_id, std::string("mat"), matptr(0),
1585  shape, len_shape);
1586  }
1587  if (nt_ > -1) {
1588  shape[0] = ((nt_ + 1) * (nt_ + 2)) / 2;
1589  // CHECK: implement store_cplx_array_to_hid with template typename T
1590  store_cplx_array_to_hid(group_id, std::string("ret"), retptr(0, 0),
1591  shape, len_shape);
1592  store_cplx_array_to_hid(group_id, std::string("les"), lesptr(0, 0),
1593  shape, len_shape);
1594  shape[0] = (nt_ + 1) * (ntau_ + 1);
1595  store_cplx_array_to_hid(group_id, std::string("tv"), tvptr(0, 0),
1596  shape, len_shape);
1597  }
1598 }
1618 template <typename T>
1619 void herm_matrix<T>::write_to_hdf5(hid_t group_id, const char *groupname) {
1620  hid_t sub_group_id = create_group(group_id, groupname);
1621  this->write_to_hdf5(sub_group_id);
1622  close_group(sub_group_id);
1623 }
1643 template <typename T>
1644 void herm_matrix<T>::write_to_hdf5(const char *filename,
1645  const char *groupname) {
1646  hid_t file_id = open_hdf5_file(filename);
1647  this->write_to_hdf5(file_id, groupname);
1648  close_hdf5_file(file_id);
1649 }
1667 template <typename T>
1668 void herm_matrix<T>::read_from_hdf5(hid_t group_id) {
1669  // -- Read dimensions
1670  int nt = read_primitive_type<int>(group_id, "nt");
1671  int ntau = read_primitive_type<int>(group_id, "ntau");
1672  int sig = read_primitive_type<int>(group_id, "sig");
1673  int size1 = read_primitive_type<int>(group_id, "size1");
1674  // RESIZE G
1675  this->resize(nt, ntau, size1);
1676  sig_ = sig;
1677  if (nt > -2) {
1678  hsize_t mat_size = (ntau + 1) * element_size_;
1679  read_primitive_type_array(group_id, "mat", mat_size, matptr(0));
1680  }
1681  if (nt > -1) {
1682  hsize_t ret_size = ((nt + 1) * (nt + 2)) / 2 * element_size_;
1683  hsize_t les_size = ((nt + 1) * (nt + 2)) / 2 * element_size_;
1684  hsize_t tv_size = ((nt + 1) * (ntau + 1)) * element_size_;
1685  read_primitive_type_array(group_id, "ret", ret_size, retptr(0, 0));
1686  read_primitive_type_array(group_id, "les", les_size, lesptr(0, 0));
1687  read_primitive_type_array(group_id, "tv", tv_size, tvptr(0, 0));
1688  }
1689 }
1709 template <typename T>
1710 void herm_matrix<T>::read_from_hdf5(hid_t group_id, const char *groupname) {
1711  hid_t sub_group_id = open_group(group_id, groupname);
1712  this->read_from_hdf5(sub_group_id);
1713  close_group(sub_group_id);
1714 }
1734 template <typename T>
1735 void herm_matrix<T>::read_from_hdf5(const char *filename,
1736  const char *groupname) {
1737  hid_t file_id = read_hdf5_file(filename);
1738  this->read_from_hdf5(file_id, groupname);
1739  close_hdf5_file(file_id);
1740 }
1761 template <typename T>
1762 void herm_matrix<T>::read_from_hdf5(int nt1, hid_t group_id) {
1763  herm_matrix<T> gtmp;
1764  gtmp.read_from_hdf5(group_id);
1765 
1766  assert(nt1 >= -1 && nt1 <= gtmp.nt() && "nt1 >= -1 && nt1 <= gtmp.nt()");
1767  assert(nt1 >= -1 && nt1 <= nt_ && "nt1 >= -1 && nt1 <= nt_");
1768  assert(gtmp.size1() == size1_ && "gtmp.size1() == size1_");
1769  assert(gtmp.element_size() == element_size_ && "gtmp.element_size() == element_size_");
1770  assert(gtmp.ntau() == ntau_ && "gtmp.ntau() == ntau_");
1771 
1772  for (int n = -1; n <= nt1; n++)
1773  this->set_timestep(n, gtmp);
1774 }
1797 template <typename T>
1798 void herm_matrix<T>::read_from_hdf5(int nt1, hid_t group_id,
1799  const char *groupname) {
1800  hid_t sub_group_id = open_group(group_id, groupname);
1801  this->read_from_hdf5(nt1, sub_group_id);
1802  close_group(sub_group_id);
1803 }
1826 template <typename T>
1827 void herm_matrix<T>::read_from_hdf5(int nt1, const char *filename,
1828  const char *groupname) {
1829  hid_t file_id = read_hdf5_file(filename);
1830  this->read_from_hdf5(nt1, file_id, groupname);
1831  close_hdf5_file(file_id);
1832 }
1851 template <typename T>
1852 void herm_matrix<T>::write_to_hdf5_slices(hid_t group_id, int dt) {
1853  assert(dt >= 1);
1854 
1855  char groupname[20];
1856  hid_t subgroup_id;
1857  for (int tstp = -1; tstp <= nt_; tstp++) {
1858  if (tstp == -1 || tstp % dt == 0) {
1859  herm_matrix_timestep_view<T> tmp(tstp, *this);
1860  std::sprintf(groupname, "t%d", tstp);
1861  subgroup_id = create_group(group_id, std::string(groupname));
1862  tmp.write_to_hdf5(subgroup_id);
1863  }
1864  }
1865 }
1886 template <typename T>
1887 void herm_matrix<T>::write_to_hdf5_slices(hid_t group_id,
1888  const char *groupname, int dt) {
1889  hid_t sub_group_id = create_group(group_id, groupname);
1890  this->write_to_hdf5_slices(sub_group_id, dt);
1891  close_group(sub_group_id);
1892 }
1913 template <typename T>
1914 void herm_matrix<T>::write_to_hdf5_slices(const char *filename,
1915  const char *groupname, int dt) {
1916  hid_t file_id = open_hdf5_file(filename);
1917  this->write_to_hdf5_slices(file_id, groupname, dt);
1918  close_hdf5_file(file_id);
1919 }
1938 template <typename T>
1939 void herm_matrix<T>::write_to_hdf5_tavtrel(hid_t group_id, int dt) {
1940  assert(dt >= 1);
1941 
1942  typedef std::complex<double> complex;
1943  char name[100];
1944  if (nt_ < 0)
1945  return;
1946  complex *ggtr = new complex[(nt_ + 1) * element_size_]; // temp storage
1947  complex *gles = new complex[(nt_ + 1) * element_size_]; // temp storage
1948  store_int_attribute_to_hid(group_id, std::string("nt"), nt_);
1949  store_int_attribute_to_hid(group_id, std::string("size1"), size1_);
1950  store_int_attribute_to_hid(group_id, std::string("size2"), size2_);
1951  store_int_attribute_to_hid(group_id, std::string("element_size"),
1952  element_size_);
1953  hid_t group_id_les = create_group(group_id, std::string("les"));
1954  hid_t group_id_gtr = create_group(group_id, std::string("gtr"));
1955  hsize_t len_shape = 3, shape[3];
1956  shape[1] = size1_;
1957  shape[2] = size2_;
1958  for (int tav = 0; tav <= nt_; tav += dt) {
1959  int len = (tav <= nt_ - tav ? tav : nt_ - tav);
1960  std::sprintf(name, "%d", tav);
1961  shape[0] = len + 1;
1962  // read data: trel2 is trel/2
1963  for (int trel2 = 0; trel2 <= len; trel2++) {
1964  cdmatrix ggtrtt, glestt;
1965  this->get_les(tav + trel2, tav - trel2, glestt);
1966  this->get_gtr(tav + trel2, tav - trel2, ggtrtt);
1967  for (int i = 0; i < element_size_; i++) {
1968  gles[trel2 * element_size_ + i] =
1969  glestt(i / size2_, i % size2_);
1970  ggtr[trel2 * element_size_ + i] =
1971  ggtrtt(i / size2_, i % size2_);
1972  }
1973  }
1974  store_cplx_array_to_hid(group_id_les, std::string(name), gles, shape,
1975  len_shape);
1976  store_cplx_array_to_hid(group_id_gtr, std::string(name), ggtr, shape,
1977  len_shape);
1978  }
1979  delete[] ggtr;
1980  delete[] gles;
1981 }
2003 template <typename T>
2004 void herm_matrix<T>::write_to_hdf5_tavtrel(hid_t group_id,
2005  const char *groupname, int dt) {
2006  hid_t sub_group_id = create_group(group_id, groupname);
2007  this->write_to_hdf5_tavtrel(sub_group_id, dt);
2008  close_group(sub_group_id);
2009 }
2031 template <typename T>
2032 void herm_matrix<T>::write_to_hdf5_tavtrel(const char *filename,
2033  const char *groupname, int dt) {
2034  hid_t file_id = open_hdf5_file(filename);
2035  this->write_to_hdf5_tavtrel(file_id, groupname, dt);
2036  close_hdf5_file(file_id);
2037 }
2038 
2039 #endif
2040 
2041 /* #######################################################################################
2042 #
2043 # SIMPLE OPERATIONS ON TIMESTEPS
2044 # NOTE: tstp IS A PHYSICAL TIME, tstp=-1 is the matsubara branch
2045 #
2046 ########################################################################################*/
2064 template <typename T>
2065 void herm_matrix<T>::set_timestep_zero(int tstp) {
2066  assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_");
2067  if (tstp == -1) {
2068  memset(matptr(0), 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
2069  } else {
2070  memset(retptr(tstp, 0), 0, sizeof(cplx) * (tstp + 1) * element_size_);
2071  memset(tvptr(tstp, 0), 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
2072  memset(lesptr(0, tstp), 0, sizeof(cplx) * (tstp + 1) * element_size_);
2073  }
2074 }
2097 template <typename T>
2098 void herm_matrix<T>::set_timestep(int tstp, herm_matrix<T> &g1) {
2099  assert(tstp >= -1 && tstp <= nt_ && tstp <= g1.nt() && "tstp >= -1 && tstp <= nt_ && tstp <= g1.nt()");
2100  assert(g1.size1() == size1_ && "g1.size1() == size1_");
2101  assert(g1.ntau() == ntau_ && "g1.ntau() == ntau_");
2102  if (tstp == -1) {
2103  memcpy(mat_, g1.mat_, sizeof(cplx) * (ntau_ + 1) * element_size_);
2104  } else {
2105  memcpy(retptr(tstp, 0), g1.retptr(tstp, 0),
2106  sizeof(cplx) * (tstp + 1) * element_size_);
2107  memcpy(tvptr(tstp, 0), g1.tvptr(tstp, 0),
2108  sizeof(cplx) * (ntau_ + 1) * element_size_);
2109  memcpy(lesptr(0, tstp), g1.lesptr(0, tstp),
2110  sizeof(cplx) * (tstp + 1) * element_size_);
2111  }
2112 }
2113 
2136 template <typename T>
2137 void herm_matrix<T>::set_timestep(int tstp,
2138  herm_matrix_timestep<T> &timestep) {
2139  cplx *x = timestep.data_;
2140  assert(tstp >= -1 && tstp <= nt_ && "(tstp >= -1 && tstp <= nt_");
2141  assert(timestep.tstp_ == tstp && timestep.ntau_ == ntau_ &&
2142  timestep.size1_ == size1_ &&
2143  "timestep.tstp_ == tstp && timestep.ntau_ == ntau_ && timestep.size1_ == size1_");
2144  if (tstp == -1) {
2145  memcpy(mat_, x, sizeof(cplx) * (ntau_ + 1) * element_size_);
2146  } else {
2147  memcpy(retptr(tstp, 0), x, sizeof(cplx) * (tstp + 1) * element_size_);
2148  memcpy(tvptr(tstp, 0), x + (tstp + 1) * element_size_,
2149  sizeof(cplx) * (ntau_ + 1) * element_size_);
2150  memcpy(lesptr(0, tstp), x + (tstp + 1 + ntau_ + 1) * element_size_,
2151  sizeof(cplx) * (tstp + 1) * element_size_);
2152  }
2153 }
2154 
2155 template <typename T>
2156 void herm_matrix<T>::get_timestep(int tstp,
2157  herm_matrix_timestep<T> &timestep) const {
2158  int len = (2 * (tstp + 1) + ntau_ + 1) * element_size_;
2159  cplx *x;
2160  assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_");
2161  if (timestep.total_size_ < len)
2162  timestep.resize(tstp, ntau_, size1_);
2163  x = timestep.data_;
2164  timestep.tstp_ = tstp;
2165  timestep.ntau_ = ntau_;
2166  timestep.size1_ = size1_;
2167  if (tstp == -1) {
2168  memcpy(x, mat_, sizeof(cplx) * (ntau_ + 1) * element_size_);
2169  } else {
2170  memcpy(x, retptr(tstp, 0), sizeof(cplx) * (tstp + 1) * element_size_);
2171  memcpy(x + (tstp + 1) * element_size_, tvptr(tstp, 0),
2172  sizeof(cplx) * (ntau_ + 1) * element_size_);
2173  memcpy(x + (tstp + 1 + ntau_ + 1) * element_size_, lesptr(0, tstp),
2174  sizeof(cplx) * (tstp + 1) * element_size_);
2175  }
2176 }
2177 
2178 template <typename T>
2179 void herm_matrix<T>::get_timestep(int tstp, herm_matrix<T> &g1) const {
2180 
2181  assert(tstp >= -1 && tstp <= nt_ && tstp <= g1.nt() && "tstp >= -1 && tstp <= nt_ && tstp <= g1.nt()");
2182  assert(g1.size1() == size1_ && "g1.size1() == size1_");
2183  assert(g1.ntau() == ntau_ && "g1.ntau() == ntau_");
2184  if (tstp == -1) {
2185  memcpy(g1.mat_, mat_, sizeof(cplx) * (ntau_ + 1) * element_size_);
2186  } else {
2187  memcpy(g1.retptr(tstp, 0), retptr(tstp, 0),
2188  sizeof(cplx) * (tstp + 1) * element_size_);
2189  memcpy(g1.tvptr(tstp, 0), tvptr(tstp, 0),
2190  sizeof(cplx) * (ntau_ + 1) * element_size_);
2191  memcpy(g1.lesptr(0, tstp), lesptr(0, tstp),
2192  sizeof(cplx) * (tstp + 1) * element_size_);
2193  }
2194 }
2195 
2196 
2221 template <typename T>
2222 void herm_matrix<T>::set_matrixelement(int tstp, int i1, int i2,
2223  herm_matrix_timestep<T> &g) {
2224  int i, sij = i1 * size2_ + i2;
2225  assert(tstp == g.tstp_ && tstp <= nt_);
2226  assert(0 <= i1 && i1 < size1_ && 0 <= i2 && i2 < size2_ && g.size1() == 1);
2227  if (tstp == -1) {
2228  for (i = 0; i <= ntau_; i++)
2229  matptr(i)[sij] = *(g.matptr(i));
2230  } else {
2231  for (i = 0; i <= tstp; i++)
2232  retptr(tstp, i)[sij] = *(g.retptr(i));
2233  for (i = 0; i <= ntau_; i++)
2234  tvptr(tstp, i)[sij] = *(g.tvptr(i));
2235  for (i = 0; i <= tstp; i++)
2236  lesptr(i, tstp)[sij] = *(g.lesptr(i));
2237  }
2238 }
2263 template <typename T>
2264 void herm_matrix<T>::set_matrixelement(int tstp, int i1, int i2,
2265  herm_matrix<T> &g) {
2266  int i, sij = i1 * size2_ + i2;
2267  assert(tstp <= g.nt() && tstp <= nt_ && "tstp <= g.nt() && tstp <= nt_");
2268  assert(0 <= i1 && i1 < size1_ && 0 <= i2 && i2 < size2_ && g.size1() == 1
2269  && "0 <= i1 && i1 < size1_ && 0 <= i2 && i2 < size2_ && g.size1() == 1");
2270  if (tstp == -1) {
2271  for (i = 0; i <= ntau_; i++)
2272  matptr(i)[sij] = *(g.matptr(i));
2273  } else {
2274  for (i = 0; i <= tstp; i++)
2275  retptr(tstp, i)[sij] = *(g.retptr(tstp, i));
2276  for (i = 0; i <= ntau_; i++)
2277  tvptr(tstp, i)[sij] = *(g.tvptr(tstp, i));
2278  for (i = 0; i <= tstp; i++)
2279  lesptr(i, tstp)[sij] = *(g.lesptr(i, tstp));
2280  }
2281 }
2282 
2312 template <typename T>
2313 void herm_matrix<T>::set_matrixelement(int tstp, int i1, int i2,
2314  herm_matrix_timestep<T> &g, int j1,
2315  int j2) {
2316  int i, sij = i1 * size2_ + i2, tij = j1 * g.size2() + j2;
2317  assert(tstp == g.tstp_ && tstp <= nt_ && "tstp == g.tstp_ && tstp <= nt_");
2318  assert(0 <= i1 && i1 < size1_ && 0 <= i2 && i2 < size2_ && "0 <= i1 && i1 < size1_ && 0 <= i2 && i2 < size2_");
2319  assert(0 <= j1 && j1 < g.size1() && 0 <= j2 && j2 < g.size1()
2320  && "0 <= j1 && j1 < g.size1() && 0 <= j2 && j2 < g.size1()");
2321  if (tstp == -1) {
2322  for (i = 0; i <= ntau_; i++)
2323  matptr(i)[sij] = g.matptr(i)[tij];
2324  } else {
2325  for (i = 0; i <= tstp; i++)
2326  retptr(tstp, i)[sij] = g.retptr(i)[tij];
2327  for (i = 0; i <= ntau_; i++)
2328  tvptr(tstp, i)[sij] = g.tvptr(i)[tij];
2329  for (i = 0; i <= tstp; i++)
2330  lesptr(i, tstp)[sij] = g.lesptr(i)[tij];
2331  }
2332 }
2333 
2364 template <typename T>
2365 void herm_matrix<T>::set_matrixelement(int tstp, int i1, int i2,
2366  herm_matrix<T> &g, int j1, int j2) {
2367  int i, sij = i1 * size2_ + i2, tij = j1 * g.size2() + j2;
2368  assert(tstp <= g.nt_ && tstp <= nt_ && "tstp <= g.nt_ && tstp <= nt_");
2369  assert(0 <= i1 && i1 < size1_ && 0 <= i2 && i2 < size2_
2370  && "0 <= i1 && i1 < size1_ && 0 <= i2 && i2 < size2_");
2371  assert(0 <= j1 && j1 < g.size1() && 0 <= j2 && j2 < g.size2()
2372  && "0 <= j1 && j1 < g.size1() && 0 <= j2 && j2 < g.size2()");
2373  if (tstp == -1) {
2374  for (i = 0; i <= ntau_; i++)
2375  matptr(i)[sij] = g.matptr(i)[tij];
2376  } else {
2377  for (i = 0; i <= tstp; i++)
2378  retptr(tstp, i)[sij] = g.retptr(tstp, i)[tij];
2379  for (i = 0; i <= ntau_; i++)
2380  tvptr(tstp, i)[sij] = g.tvptr(tstp, i)[tij];
2381  for (i = 0; i <= tstp; i++)
2382  lesptr(i, tstp)[sij] = g.lesptr(i, tstp)[tij];
2383  }
2384 }
2385 
2386 
2410 template<typename T> void herm_matrix<T>::set_matrixelement(int i1,int i2,herm_matrix<T> &g){
2411  assert(nt_ == g.nt_);
2412  for(int tstp=-1;tstp<=nt_;tstp++){
2413  this->set_matrixelement(tstp,i1,i2,g);
2414  }
2415 }
2416 
2445 template<typename T> void herm_matrix<T>::set_matrixelement(int i1,int i2,herm_matrix<T> &g,int j1,int j2){
2446  assert(nt_ == g.nt_);
2447  for(int tstp=-1;tstp<=nt_;tstp++){
2448  this->set_matrixelement(tstp,i1,i2,g,j1,j2);
2449  }
2450 }
2451 
2452 // Set the submatrix of the Green'd function
2453 // G_{j_1(k) j_2(k)}=G_{i_1(k) i_2(k) }, where k is an iterator over subspace
2454 
2455 
2483 template<typename T> void herm_matrix<T>::set_submatrix(std::vector<int> &i1,
2484  std::vector<int> &i2,herm_matrix<T> &g,std::vector<int> &j1,std::vector<int> &j2){
2485  assert(nt_ == g.nt_);
2486  assert(i1.size()==i2.size() && i1.size()==j1.size() && j1.size()==j2.size());
2487  assert(size1_*size2_==i1.size());
2488  for(int k1=0;k1<i1.size();k1++){
2489  this->set_matrixelement(i1[k1],i2[k1],g,j1[k1],j2[k1]);
2490  }
2491 }
2492 
2522 template<typename T> void herm_matrix<T>::set_submatrix(int tstp, std::vector<int> &i1,
2523  std::vector<int> &i2,herm_matrix<T> &g,std::vector<int> &j1,std::vector<int> &j2){
2524  assert(tstp <= nt_);
2525  assert(nt_ == g.nt_);
2526  assert(i1.size()==i2.size() && i1.size()==j1.size() && j1.size()==j2.size());
2527  assert(size1_*size2_==i1.size());
2528  for(int k1=0;k1<i1.size();k1++){
2529  this->set_matrixelement(tstp,i1[k1],i2[k1],g,j1[k1],j2[k1]);
2530  }
2531 }
2532 
2533 
2534 #define HERM_MATRIX_INCR_TSTP \
2535  if (alpha == cplx(1.0, 0.0)) { \
2536  for (i = 0; i < len; i++) \
2537  x0[i] += x[i]; \
2538  } else { \
2539  for (i = 0; i < len; i++) \
2540  x0[i] += alpha * x[i]; \
2541  }
2542 
2543 
2544 
2566 template <typename T>
2567 void herm_matrix<T>::incr_timestep(int tstp,
2568  herm_matrix_timestep<T> &timestep,
2569  std::complex<T> alpha) {
2570  int i, len;
2571  cplx *x, *x0;
2572  assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_");
2573  assert(timestep.tstp_ == tstp && timestep.ntau_ == ntau_ && timestep.size1_ == size1_
2574  && "timestep.tstp_ == tstp && timestep.ntau_ == ntau_ && timestep.size1_ == size1_");
2575  if (tstp == -1) {
2576  len = (ntau_ + 1) * element_size_;
2577  x0 = matptr(0);
2578  x = timestep.data_;
2579  HERM_MATRIX_INCR_TSTP
2580  } else {
2581  len = (tstp + 1) * element_size_;
2582  x0 = retptr(tstp, 0);
2583  x = timestep.data_;
2584  HERM_MATRIX_INCR_TSTP
2585  len = (ntau_ + 1) * element_size_;
2586  x0 = tvptr(tstp, 0);
2587  x = timestep.data_ + (tstp + 1) * element_size_;
2588  HERM_MATRIX_INCR_TSTP
2589  len = (tstp + 1) * element_size_;
2590  x0 = lesptr(0, tstp);
2591  x = timestep.data_ + (tstp + 1 + ntau_ + 1) * element_size_;
2592  HERM_MATRIX_INCR_TSTP
2593  }
2594 }
2595 
2615 template <typename T>
2616 void herm_matrix<T>::incr_timestep(int tstp,
2617  herm_matrix_timestep<T> &timestep) {
2618  int i, len;
2619  cplx *x, *x0;
2620  cplx alpha = cplx(1.0, 0.0);
2621  assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_");
2622  assert(timestep.tstp_ == tstp && timestep.ntau_ == ntau_ && timestep.size1_ == size1_
2623  && "timestep.tstp_ == tstp && timestep.ntau_ == ntau_ && timestep.size1_ == size1_");
2624  if (tstp == -1) {
2625  len = (ntau_ + 1) * element_size_;
2626  x0 = matptr(0);
2627  x = timestep.data_;
2628  HERM_MATRIX_INCR_TSTP
2629  } else {
2630  len = (tstp + 1) * element_size_;
2631  x0 = retptr(tstp, 0);
2632  x = timestep.data_;
2633  HERM_MATRIX_INCR_TSTP
2634  len = (ntau_ + 1) * element_size_;
2635  x0 = tvptr(tstp, 0);
2636  x = timestep.data_ + (tstp + 1) * element_size_;
2637  HERM_MATRIX_INCR_TSTP
2638  len = (tstp + 1) * element_size_;
2639  x0 = lesptr(0, tstp);
2640  x = timestep.data_ + (tstp + 1 + ntau_ + 1) * element_size_;
2641  HERM_MATRIX_INCR_TSTP
2642  }
2643 }
2644 
2645 #undef HERM_MATRIX_INCR_TSTP
2646 
2668 template <typename T>
2669 void herm_matrix<T>::incr_timestep(int tstp, herm_matrix<T> &g,
2670  std::complex<T> alpha) {
2671  int m;
2672  cplx *x, *x0;
2673  assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_");
2674  assert(g.nt_ >= tstp && g.ntau_ == ntau_ && g.size1_ == size1_
2675  && "g.nt_ >= tstp && g.ntau_ == ntau_ && g.size1_ == size1_");
2676  if (tstp == -1) {
2677  x0 = matptr(0);
2678  x = g.matptr(0);
2679  for (m = 0; m <= ntau_; m++) {
2680  element_incr<T, LARGESIZE>(size1_, x0 + m * element_size_, alpha,
2681  x + m * element_size_);
2682  }
2683  } else {
2684  x0 = retptr(tstp, 0);
2685  x = g.retptr(tstp, 0);
2686  for (m = 0; m <= tstp; m++) {
2687  element_incr<T, LARGESIZE>(size1_, x0 + m * element_size_, alpha,
2688  x + m * element_size_);
2689  }
2690  x0 = tvptr(tstp, 0);
2691  x = g.tvptr(tstp, 0);
2692  for (m = 0; m <= ntau_; m++) {
2693  element_incr<T, LARGESIZE>(size1_, x0 + m * element_size_, alpha,
2694  x + m * element_size_);
2695  }
2696  x0 = lesptr(0, tstp);
2697  x = g.lesptr(0, tstp);
2698  for (m = 0; m <= tstp; m++) {
2699  element_incr<T, LARGESIZE>(size1_, x0 + m * element_size_, alpha,
2700  x + m * element_size_);
2701  }
2702  }
2703 }
2704 
2724 template <typename T>
2725 void herm_matrix<T>::incr_timestep(int tstp, herm_matrix<T> &g) {
2726  int m;
2727  cplx *x, *x0;
2728  cplx alpha = cplx(1.0, 0.0);
2729  assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_");
2730  assert(g.nt_ >= tstp && g.ntau_ == ntau_ && g.size1_ == size1_
2731  && "g.nt_ >= tstp && g.ntau_ == ntau_ && g.size1_ == size1_");
2732  if (tstp == -1) {
2733  x0 = matptr(0);
2734  x = g.matptr(0);
2735  for (m = 0; m <= ntau_; m++) {
2736  element_incr<T, LARGESIZE>(size1_, x0 + m * element_size_, alpha,
2737  x + m * element_size_);
2738  }
2739  } else {
2740  x0 = retptr(tstp, 0);
2741  x = g.retptr(tstp, 0);
2742  for (m = 0; m <= tstp; m++) {
2743  element_incr<T, LARGESIZE>(size1_, x0 + m * element_size_, alpha,
2744  x + m * element_size_);
2745  }
2746  x0 = tvptr(tstp, 0);
2747  x = g.tvptr(tstp, 0);
2748  for (m = 0; m <= ntau_; m++) {
2749  element_incr<T, LARGESIZE>(size1_, x0 + m * element_size_, alpha,
2750  x + m * element_size_);
2751  }
2752  x0 = lesptr(0, tstp);
2753  x = g.lesptr(0, tstp);
2754  for (m = 0; m <= tstp; m++) {
2755  element_incr<T, LARGESIZE>(size1_, x0 + m * element_size_, alpha,
2756  x + m * element_size_);
2757  }
2758  }
2759 }
2760 
2761 
2781 template <typename T>
2782 void herm_matrix<T>::incr(herm_matrix<T> &g, std::complex<T> alpha) {
2783  assert(g.nt_ >= nt_ && g.ntau_ == ntau_ && g.size1_ == size1_
2784  && "g.nt_ >= nt_ && g.ntau_ == ntau_ && g.size1_ == size1_");
2785  for (int m = -1; m <= nt_; m++)
2786  this->incr_timestep(m, g, alpha);
2787 }
2788 
2805 template <typename T>
2806 void herm_matrix<T>::incr(herm_matrix<T> &g) {
2807  assert(g.nt_ >= nt_ && g.ntau_ == ntau_ && g.size1_ == size1_
2808  && "g.nt_ >= nt_ && g.ntau_ == ntau_ && g.size1_ == size1_");
2809  cplx alpha = (1.0,0.0);
2810  for (int m = -1; m <= nt_; m++)
2811  this->incr_timestep(m, g, alpha);
2812 }
2813 
2815 /** \brief <b> Left-multiplies the `herm_matrix` with contour function. </b>
2816 *
2817 * <!-- ====== DOCUMENTATION ====== -->
2818 *
2819 * \par Purpose
2820 * <!-- ========= -->
2821 *
2822 * Performs the operation \f$C(t,t^\prime) \rightarrow w F(t)C(t,t^\prime)\f$, where \f$C(t,t^\prime)\f$ is the
2823 * `herm_matrix`, \f$F(t)\f$ is a `function` given in pointer format and \f$w\f$ is a real weight. The operation is performed
2824 * at given time step `tstp`. This is a lower-level routine. The interface where \f$F(t)\f$ is supplied as
2825 * contour function is preferred.
2826 *
2827 * <!-- ARGUMENTS
2828 * ========= -->
2829 *
2830 * @param tstp
2831 * > [int] The time step at which \f$F(t)\f$ and \f$C(t,t^\prime)\f$ are multiplied.
2832 * @param *f0
2833 * > [complex<T>] Pointer to \f$F(-\mathrm{i}\beta)\f$ on the Matsubara axis (this is just a constant matrix).
2834 * @param *ft
2835 * > [complex<T>] Pointer to \f$F(t)\f$ on the real axis. It is assumed that ft+t*element_size_ points to \f$ F(t) \f$.
2836 * @param weight
2837 * > [T] The weight as above.
2838 */
2839 template <typename T>
2840 void herm_matrix<T>::left_multiply(int tstp, std::complex<T> *f0,
2841  std::complex<T> *ft, T weight) {
2842  int m;
2843  cplx *xtemp, *ftemp, *x0;
2844  xtemp = new cplx[element_size_];
2845  assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_");
2846  if (tstp == -1) {
2847  x0 = matptr(0);
2848  for (m = 0; m <= ntau_; m++) {
2849  element_mult<T, LARGESIZE>(size1_, xtemp, f0,
2850  x0 + m * element_size_);
2851  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
2852  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
2853  }
2854  } else {
2855  ftemp = ft + tstp * element_size_;
2856  x0 = retptr(tstp, 0);
2857  for (m = 0; m <= tstp; m++) {
2858  element_mult<T, LARGESIZE>(size1_, xtemp, ftemp,
2859  x0 + m * element_size_);
2860  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
2861  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
2862  }
2863  x0 = tvptr(tstp, 0);
2864  for (m = 0; m <= ntau_; m++) {
2865  element_mult<T, LARGESIZE>(size1_, xtemp, ftemp,
2866  x0 + m * element_size_);
2867  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
2868  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
2869  }
2870  x0 = lesptr(0, tstp);
2871  for (m = 0; m <= tstp; m++) {
2872  element_mult<T, LARGESIZE>(size1_, xtemp, ft + m * element_size_,
2873  x0 + m * element_size_);
2874  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
2875  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
2876  }
2877  }
2878  delete[] xtemp;
2879 }
2880 
2902 template <typename T>
2903 void herm_matrix<T>::left_multiply(int tstp, function<T> &ft, T weight) {
2904  switch (size1_) {
2905  case 1:
2906  left_multiply_(tstp, ft, weight, std::integral_constant<int, 1>{});
2907  break;
2908  case 2:
2909  left_multiply_(tstp, ft, weight, std::integral_constant<int, 2>{});
2910  break;
2911  case 3:
2912  left_multiply_(tstp, ft, weight, std::integral_constant<int, 3>{});
2913  break;
2914  case 4:
2915  left_multiply_(tstp, ft, weight, std::integral_constant<int, 4>{});
2916  break;
2917  case 5:
2918  left_multiply_(tstp, ft, weight, std::integral_constant<int, 5>{});
2919  break;
2920  case 8:
2921  left_multiply_(tstp, ft, weight, std::integral_constant<int, 8>{});
2922  break;
2923  default:
2924  left_multiply_(tstp, ft, weight, std::integral_constant<int, LARGESIZE>{});
2925  break;
2926  }
2927 }
2929 /** \brief <b> Left-multiplies the `herm_matrix` with contour function. Internal version, not intended
2930 * to be called directly. </b>
2931 *
2932 * <!-- ====== DOCUMENTATION ====== -->
2933 *
2934 * \par Purpose
2935 * <!-- ========= -->
2936 *
2937 * Performs the operation \f$C(t,t^\prime) \rightarrow w F(t)C(t,t^\prime)\f$, where \f$C(t,t^\prime)\f$ is the
2938 * `herm_matrix`, \f$F(t)\f$ is a contour `function` and \f$w\f$ is a real weight. The operation is performed
2939 * at given time step `tstp`. This is an internal version with specializations with respect to the matrix size.
2940 *
2941 * \note parameter SIZE determines the matrix size specified by the top-level interface.
2942 *
2943 * <!-- ARGUMENTS
2944 * ========= -->
2945 *
2946 * @param tstp
2947 * > [int] The time step at which \f$F(t)\f$ and \f$C(t,t^\prime)\f$ are multiplied.
2948 * @param ft
2949 * > [function] The contour function.
2950 * @param weight
2951 * > [T] The weight as above.
2952 */
2953 template <typename T>
2954 template <int SIZE>
2955 void herm_matrix<T>::left_multiply_(int tstp, function<T> &ft, T weight,
2956  std::integral_constant<int, SIZE>) {
2957  int m;
2958  cplx *xtemp, *ftemp, *x0, *f0;
2959  xtemp = new cplx[element_size_];
2960  assert(tstp >= -1 && tstp <= nt_ && ft.nt() >= tstp &&
2961  ft.size1() == size1_ && ft.size2() == size2_
2962  && "tstp >= -1 && tstp <= nt_ && ft.nt() >= tstp && ft.size1() == size1_ && ft.size2() == size2_");
2963  f0 = ft.ptr(-1);
2964  if (tstp == -1) {
2965  x0 = matptr(0);
2966  for (m = 0; m <= ntau_; m++) {
2967  element_mult<T, SIZE>(size1_, xtemp, f0, x0 + m * element_size_);
2968  element_smul<T, SIZE>(size1_, xtemp, weight);
2969  element_set<T, SIZE>(size1_, x0 + m * element_size_, xtemp);
2970  }
2971  } else {
2972  ftemp = ft.ptr(tstp);
2973  x0 = retptr(tstp, 0);
2974  for (m = 0; m <= tstp; m++) {
2975  element_mult<T, SIZE>(size1_, xtemp, ftemp, x0 + m * element_size_);
2976  element_smul<T, SIZE>(size1_, xtemp, weight);
2977  element_set<T, SIZE>(size1_, x0 + m * element_size_, xtemp);
2978  }
2979  x0 = tvptr(tstp, 0);
2980  for (m = 0; m <= ntau_; m++) {
2981  element_mult<T, SIZE>(size1_, xtemp, ftemp, x0 + m * element_size_);
2982  element_smul<T, SIZE>(size1_, xtemp, weight);
2983  element_set<T, SIZE>(size1_, x0 + m * element_size_, xtemp);
2984  }
2985  x0 = lesptr(0, tstp);
2986  for (m = 0; m <= tstp; m++) {
2987  element_mult<T, SIZE>(size1_, xtemp, ft.ptr(m), x0 + m * element_size_);
2988  element_smul<T, SIZE>(size1_, xtemp, weight);
2989  element_set<T, SIZE>(size1_, x0 + m * element_size_, xtemp);
2990  }
2991  }
2992  delete[] xtemp;
2993 }
2994 
3016 template <typename T>
3017 void herm_matrix<T>::left_multiply_hermconj(int tstp, function<T> &ft,
3018  T weight) {
3019  int m;
3020  cplx *xtemp, *ftemp, *x0, *f0, *fcc;
3021  xtemp = new cplx[element_size_];
3022  fcc = new cplx[element_size_];
3023  assert(tstp >= -1 && tstp <= nt_ && ft.nt() >= tstp &&
3024  ft.size1() == size1_ && ft.size2() == size2_ &&
3025  "tstp >= -1 && tstp <= nt_ && ft.nt() >= tstp && ft.size1() == size1_ && ft.size2() == size2_");
3026 
3027  f0 = ft.ptr(-1);
3028  element_conj<T, LARGESIZE>(size1_, fcc, f0);
3029  if (tstp == -1) {
3030  x0 = matptr(0);
3031  for (m = 0; m <= ntau_; m++) {
3032  element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
3033  x0 + m * element_size_);
3034  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3035  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3036  }
3037  } else {
3038  ftemp = ft.ptr(tstp);
3039  element_conj<T, LARGESIZE>(size1_, fcc, ftemp);
3040  x0 = retptr(tstp, 0);
3041  for (m = 0; m <= tstp; m++) {
3042  element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
3043  x0 + m * element_size_);
3044  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3045  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3046  }
3047  x0 = tvptr(tstp, 0);
3048  for (m = 0; m <= ntau_; m++) {
3049  element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
3050  x0 + m * element_size_);
3051  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3052  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3053  }
3054  x0 = lesptr(0, tstp);
3055  for (m = 0; m <= tstp; m++) {
3056  element_conj<T, LARGESIZE>(size1_, fcc, ft.ptr(m));
3057  element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
3058  x0 + m * element_size_);
3059  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3060  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3061  }
3062  }
3063  delete[] xtemp;
3064  delete[] fcc;
3065 }
3067 /** \brief <b> Right-multiplies the `herm_matrix` with contour function. </b>
3068 *
3069 * <!-- ====== DOCUMENTATION ====== -->
3070 *
3071 * \par Purpose
3072 * <!-- ========= -->
3073 *
3074 * Performs the operation \f$C(t,t^\prime) \rightarrow w C(t,t^\prime) F(t^\prime)\f$, where \f$C(t,t^\prime)\f$ is the
3075 * `herm_matrix`, \f$F(t)\f$ is a `function` given in pointer format and \f$w\f$ is a real weight. The operation is performed
3076 * at given time step `tstp`. This is a lower-level routine. The interface where \f$F(t)\f$ is supplied as
3077 * contour function is preferred.
3078 *
3079 * <!-- ARGUMENTS
3080 * ========= -->
3081 *
3082 * @param tstp
3083 * > [int] The time step at which \f$F(t)\f$ and \f$C(t,t^\prime)\f$ are multiplied.
3084 * @param *f0
3085 * > [complex<T>] Pointer to \f$F(-\mathrm{i}\beta)\f$ on the Matsubara axis (this is just a constant matrix).
3086 * @param *ft
3087 * > [complex<T>] Pointer to \f$F(t)\f$ on the real axis. It is assumed that ft+t*element_size_ points to \f$ F(t) \f$.
3088 * @param weight
3089 * > [T] The weight as above.
3090 */
3091 template <typename T>
3092 void herm_matrix<T>::right_multiply(int tstp, std::complex<T> *f0,
3093  std::complex<T> *ft, T weight) {
3094  int m;
3095  cplx *xtemp, *ftemp, *x0;
3096  xtemp = new cplx[element_size_];
3097  assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_");
3098  if (tstp == -1) {
3099  x0 = matptr(0);
3100  for (m = 0; m <= ntau_; m++) {
3101  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
3102  f0);
3103  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3104  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3105  }
3106  } else {
3107  x0 = retptr(tstp, 0);
3108  for (m = 0; m <= tstp; m++) {
3109  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
3110  ft + m * element_size_);
3111  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3112  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3113  }
3114  x0 = tvptr(tstp, 0);
3115  for (m = 0; m <= ntau_; m++) {
3116  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
3117  f0);
3118  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3119  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3120  }
3121  ftemp = ft + tstp * element_size_;
3122  x0 = lesptr(0, tstp);
3123  for (m = 0; m <= tstp; m++) {
3124  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
3125  ftemp);
3126  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3127  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3128  }
3129  }
3130  delete[] xtemp;
3131 }
3153 template <typename T>
3154 void herm_matrix<T>::right_multiply(int tstp, function<T> &ft, T weight) {
3155  switch (size1_) {
3156  case 1:
3157  right_multiply_(tstp, ft, weight, std::integral_constant<int, 1>{});
3158  break;
3159  case 2:
3160  right_multiply_(tstp, ft, weight, std::integral_constant<int, 2>{});
3161  break;
3162  case 3:
3163  right_multiply_(tstp, ft, weight, std::integral_constant<int, 3>{});
3164  break;
3165  case 4:
3166  right_multiply_(tstp, ft, weight, std::integral_constant<int, 4>{});
3167  break;
3168  case 5:
3169  right_multiply_(tstp, ft, weight, std::integral_constant<int, 5>{});
3170  break;
3171  case 8:
3172  right_multiply_(tstp, ft, weight, std::integral_constant<int, 8>{});
3173  break;
3174  default:
3175  right_multiply_(tstp, ft, weight, std::integral_constant<int, LARGESIZE>{});
3176  break;
3177  }
3178 }
3180 /** \brief <b> Right-multiplies the `herm_matrix` with contour function.
3181 * Internal version, not intended to be called directly. </b>
3182 *
3183 * <!-- ====== DOCUMENTATION ====== -->
3184 *
3185 * \par Purpose
3186 * <!-- ========= -->
3187 *
3188 * Performs the operation \f$C(t,t^\prime) \rightarrow w C(t,t^\prime)F(t^\prime)\f$, where \f$C(t,t^\prime)\f$ is the
3189 * `herm_matrix`, \f$F(t)\f$ is a contour `function` and \f$w\f$ is a real weight. The operation is performed
3190 * at given time step `tstp`. This is an internal version with specializations with respect to the matrix size.
3191 *
3192 * \note Parameter SIZE determines the matrix size specified by the top-level interface.
3193 *
3194 * <!-- ARGUMENTS
3195 * ========= -->
3196 *
3197 * @param tstp
3198 * > The time step at which \f$F(t)\f$ and \f$C(t,t^\prime)\f$ are multiplied.
3199 * @param ft
3200 * > The contour function.
3201 * @param weight
3202 * > The weight as above.
3203 */
3204 template <typename T>
3205 template <int SIZE>
3206 void herm_matrix<T>::right_multiply_(int tstp, function<T> &ft, T weight,
3207  std::integral_constant<int, SIZE>) {
3208  int m;
3209  cplx *xtemp, *ftemp, *x0, *f0;
3210  xtemp = new cplx[element_size_];
3211  assert(ft.size1() == size1_ && "ft.size1() == size1_");
3212  assert(ft.nt() >= tstp && "ft.nt() >= tstp");
3213  assert(tstp <= nt_ && "tstp <= nt_");
3214  assert(tstp >= -1 && "tstp >= -1");
3215  f0 = ft.ptr(-1);
3216  if (tstp == -1) {
3217  x0 = matptr(0);
3218  for (m = 0; m <= ntau_; m++) {
3219  element_mult<T, SIZE>(size1_, xtemp, x0 + m * element_size_, f0);
3220  element_smul<T, SIZE>(size1_, xtemp, weight);
3221  element_set<T, SIZE>(size1_, x0 + m * element_size_, xtemp);
3222  }
3223  } else {
3224  x0 = retptr(tstp, 0);
3225  for (m = 0; m <= tstp; m++) {
3226  element_mult<T, SIZE>(size1_, xtemp, x0 + m * element_size_, ft.ptr(m));
3227  element_smul<T, SIZE>(size1_, xtemp, weight);
3228  element_set<T, SIZE>(size1_, x0 + m * element_size_, xtemp);
3229  }
3230  x0 = tvptr(tstp, 0);
3231  for (m = 0; m <= ntau_; m++) {
3232  element_mult<T, SIZE>(size1_, xtemp, x0 + m * element_size_, f0);
3233  element_smul<T, SIZE>(size1_, xtemp, weight);
3234  element_set<T, SIZE>(size1_, x0 + m * element_size_, xtemp);
3235  }
3236  ftemp = ft.ptr(tstp);
3237  x0 = lesptr(0, tstp);
3238  for (m = 0; m <= tstp; m++) {
3239  element_mult<T, SIZE>(size1_, xtemp, x0 + m * element_size_, ftemp);
3240  element_smul<T, SIZE>(size1_, xtemp, weight);
3241  element_set<T, SIZE>(size1_, x0 + m * element_size_, xtemp);
3242  }
3243  }
3244  delete[] xtemp;
3245 }
3246 
3269 template <typename T>
3270 void herm_matrix<T>::right_multiply_hermconj(int tstp, function<T> &ft,
3271  T weight) {
3272  int m;
3273  cplx *xtemp, *ftemp, *x0, *f0, *fcc;
3274  xtemp = new cplx[element_size_];
3275  fcc = new cplx[element_size_];
3276  assert(ft.size1() == size1_ && "ft.size1() == size1_");
3277  assert(ft.nt() >= tstp && "ft.nt() >= tstp");
3278  assert(tstp <= nt_&& "ft.nt() >= tstp");
3279  assert(tstp >= -1&& "ft.nt() >= tstp");
3280  f0 = ft.ptr(-1);
3281  element_conj<T, LARGESIZE>(size1_, fcc, f0);
3282  if (tstp == -1) {
3283  x0 = matptr(0);
3284  for (m = 0; m <= ntau_; m++) {
3285  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
3286  fcc);
3287  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3288  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3289  }
3290  } else {
3291  x0 = retptr(tstp, 0);
3292  for (m = 0; m <= tstp; m++) {
3293  element_conj<T, LARGESIZE>(size1_, fcc, ft.ptr(m));
3294  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
3295  fcc);
3296  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3297  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3298  }
3299  x0 = tvptr(tstp, 0);
3300  element_conj<T, LARGESIZE>(size1_, fcc, f0);
3301  for (m = 0; m <= ntau_; m++) {
3302  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
3303  fcc);
3304  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3305  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3306  }
3307  ftemp = ft.ptr(tstp);
3308  element_conj<T, LARGESIZE>(size1_, fcc, ftemp);
3309  x0 = lesptr(0, tstp);
3310  for (m = 0; m <= tstp; m++) {
3311  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
3312  fcc);
3313  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
3314  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
3315  }
3316  }
3317  delete[] fcc;
3318  delete[] xtemp;
3319 }
3320 
3338 template <typename T>
3339 void herm_matrix<T>::smul(int tstp, T weight) {
3340  int m;
3341  cplx *x0;
3342  assert(tstp >= -1 && tstp <= nt_);
3343  if (tstp == -1) {
3344  x0 = matptr(0);
3345  for (m = 0; m <= ntau_; m++) {
3346  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
3347  weight);
3348  }
3349  } else {
3350  x0 = retptr(tstp, 0);
3351  for (m = 0; m <= tstp; m++) {
3352  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
3353  weight);
3354  }
3355  x0 = tvptr(tstp, 0);
3356  for (m = 0; m <= ntau_; m++) {
3357  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
3358  weight);
3359  }
3360  x0 = lesptr(0, tstp);
3361  for (m = 0; m <= tstp; m++) {
3362  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
3363  weight);
3364  }
3365  }
3366 }
3367 
3386 template <typename T>
3387 void herm_matrix<T>::smul(int tstp, std::complex<T> weight) {
3388  int m;
3389  cplx *x0;
3390  assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_");
3391  if (tstp == -1) {
3392  x0 = matptr(0);
3393  for (m = 0; m <= ntau_; m++) {
3394  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
3395  weight);
3396  }
3397  } else {
3398  x0 = retptr(tstp, 0);
3399  for (m = 0; m <= tstp; m++) {
3400  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
3401  weight);
3402  }
3403  x0 = tvptr(tstp, 0);
3404  for (m = 0; m <= ntau_; m++) {
3405  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
3406  weight);
3407  }
3408  x0 = lesptr(0, tstp);
3409  for (m = 0; m <= tstp; m++) {
3410  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
3411  weight);
3412  }
3413  }
3414 }
3415 
3416 /* #######################################################################################
3417 #
3418 # MPI UTILS
3419 #
3420 ########################################################################################*/
3421 #if CNTR_USE_MPI == 1
3422 
3441 template <typename T>
3442 void herm_matrix<T>::Reduce_timestep(int tstp, int root) {
3443  assert(tstp <= nt_);
3444 
3445  herm_matrix_timestep<T> Gtemp;
3446  Gtemp.resize(tstp, ntau_, size1_);
3447  this->get_timestep(tstp, Gtemp);
3448 
3449  Gtemp.Reduce_timestep(tstp, root);
3450 
3451  this->set_timestep(tstp, Gtemp);
3452 }
3453 
3454 
3472 template <typename T>
3473 void herm_matrix<T>::Bcast_timestep(int tstp, int root) {
3474  int numtasks, taskid;
3475  herm_matrix_timestep<T> Gtemp;
3476  MPI_Comm_size(MPI_COMM_WORLD, &numtasks);
3477  MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
3478  Gtemp.resize(tstp, ntau_, size1_);
3479  if (taskid == root)
3480  this->get_timestep(tstp, Gtemp);
3481  if (sizeof(T) == sizeof(double))
3482  MPI_Bcast(Gtemp.data_, Gtemp.total_size_,
3483  MPI_DOUBLE_COMPLEX, root, MPI_COMM_WORLD);
3484  else
3485  MPI_Bcast(Gtemp.data_, Gtemp.total_size_, MPI_COMPLEX,
3486  root, MPI_COMM_WORLD);
3487  if (taskid != root)
3488  this->set_timestep(tstp, Gtemp);
3489 }
3490 
3510 template <typename T>
3511 void herm_matrix<T>::Send_timestep(int tstp, int dest, int tag) {
3512  int taskid;
3513  MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
3514  if (!(taskid == dest)) {
3515  herm_matrix_timestep<T> Gtemp;
3516  Gtemp.resize(tstp, ntau_, size1_);
3517  this->get_timestep(tstp, Gtemp);
3518  if (sizeof(T) == sizeof(double))
3519  MPI_Send(Gtemp.data_, Gtemp.total_size_,
3520  MPI_DOUBLE_COMPLEX, dest, tag, MPI_COMM_WORLD);
3521  else
3522  MPI_Send(Gtemp.data_, Gtemp.total_size_, MPI_COMPLEX,
3523  dest, tag, MPI_COMM_WORLD);
3524  }
3525 }
3526 
3546 template <typename T>
3547 void herm_matrix<T>::Recv_timestep(int tstp, int root, int tag) {
3548  int taskid;
3549  MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
3550  if (!(taskid == root)) {
3551  herm_matrix_timestep<T> Gtemp;
3552  Gtemp.resize(tstp, ntau_, size1_);
3553  if (sizeof(T) == sizeof(double))
3554  MPI_Recv(Gtemp.data_, Gtemp.total_size_,
3555  MPI_DOUBLE_COMPLEX, root, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
3556  else
3557  MPI_Recv(Gtemp.data_, Gtemp.total_size_, MPI_COMPLEX,
3558  root, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
3559  this->set_timestep(tstp, Gtemp);
3560  }
3561 }
3562 #endif
3563 
3564 } // namespace cntr
3565 
3566 #endif // CNTR_HERM_MATRIX_IMPL_H
void clear(void)
Sets all values to zero.
void get_mat(const int m, std::complex< T > &G_mat, herm_matrix< T > &G, herm_matrix< T > &Gcc)
std::complex< T > cplx
void resize_nt(int nt)
Resizes herm_matrix object with respect to the number of time points nt. Internal routine; see top-l...
void resize_discard(int nt, int ntau, int size1)
Discards the herm_matrix object and resize with respect to the number of time points nt...
std::complex< double > cplx
Definition: fourier.cpp:11
void get_ret(const int i, const int j, std::complex< T > &G_ret, herm_matrix< T > &G, herm_matrix< T > &Gcc)
Class herm_matrix for two-time contour objects with hermitian symmetry.
void print_to_file(const char *file, int precision=16) const
Saves the herm_matrix to file in text format.
herm_matrix & operator=(const herm_matrix &g)
void get_les(const int i, const int j, std::complex< T > &G_les, herm_matrix< T > &G, herm_matrix< T > &Gcc)
void resize(int nt, int ntau, int size1)
Resizes herm_matrix object with respect to the number of time points nt, points on the Matsubara bra...
void read_from_file(const char *file)
Reads the herm_matrix from file in text format.