NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_herm_matrix_timestep_view_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_HERM_TIMESTEP_VIEW_IMPL_H
2 #define CNTR_HERM_TIMESTEP_VIEW_IMPL_H
3 
5 #include "cntr_elements.hpp"
6 //#include "cntr_exception.hpp"
7 
8 namespace cntr {
9 
10 #define CPLX std::complex<T>
11 template <typename T>
12 /* #######################################################################################
13 #
14 # CONSTRUCTION/DESTRUCTION
15 #
16 ########################################################################################*/
18  ret_ = 0;
19  les_ = 0;
20  tv_ = 0;
21  mat_ = 0;
22  ntau_ = 0;
23  tstp_ = -2;
24  size1_ = 0;
25  size2_ = 0;
26  element_size_ = 0;
27  sig_ = -1;
28 }
29 template <typename T>
31  // donothing
32 }
33 
53 template <typename T>
55  const herm_matrix_timestep_view &g) {
56  tstp_ = g.tstp_;
57  ntau_ = g.ntau_;
58  size1_ = g.size1_;
59  size2_ = g.size1_;
60  element_size_ = size1_ * size1_;
61  sig_ = g.sig_;
62  // copy pointers
63  ret_ = g.ret_;
64  les_ = g.les_;
65  tv_ = g.tv_;
66  mat_ = g.mat_;
67 }
68 
69 
70 
95 template <typename T> herm_matrix_timestep_view<T>::herm_matrix_timestep_view(int tstp,int ntau,int size1,int size2,int sig){
96 
97  assert(size1>=0 && tstp>=-1 && ntau>=0);
98  ret_ = 0;
99  les_ = 0;
100  tv_ = 0;
101  mat_ = 0;
102  size1_=size1;
103  size2_=size2;
104  element_size_=size1*size2;
105  tstp_=tstp;
106  ntau_=ntau;
107  sig_=sig;
108 
109 }
110 
111 
112 
132 template <typename T>
135  if (this == &g)
136  return *this;
137  tstp_ = g.tstp_;
138  ntau_ = g.ntau_;
139  size1_ = g.size1_;
140  size2_ = g.size1_;
141  sig_ = g.sig_;
142  element_size_ = size1_ * size1_;
143  // copy pointers
144  ret_ = g.ret_;
145  les_ = g.les_;
146  tv_ = g.tv_;
147  mat_ = g.mat_;
148  return *this;
149 }
150 
171 template <typename T>
174  tstp_ = g.tstp_;
175  ntau_ = g.ntau_;
176  size1_ = g.size1_;
177  size2_ = g.size2_;
178  element_size_ = size1_ * size2_;
179  sig_ = g.sig_;
180  if (tstp_ == -1) {
181  mat_ = g.matptr(0);
182  ret_ = 0;
183  les_ = 0;
184  tv_ = 0;
185  } else if (tstp_ >= 0) {
186  mat_ = 0;
187  ret_ = g.retptr(0);
188  les_ = g.lesptr(0);
189  tv_ = g.tvptr(0);
190  }
191 }
192 
217 template <typename T>
219  int tstp, herm_matrix_timestep<T> &g, bool check_assert) {
220  if(check_assert) assert(tstp==g.tstp_);
221 
222  tstp_ = g.tstp_;
223  ntau_ = g.ntau_;
224  size1_ = g.size1_;
225  size2_ = g.size2_;
226  element_size_ = size1_ * size2_;
227  sig_ = g.sig_;
228  if (tstp_ == -1) {
229  mat_ = g.matptr(0);
230  ret_ = 0;
231  les_ = 0;
232  tv_ = 0;
233  } else if (tstp_ >= 0) {
234  mat_ = 0;
235  ret_ = g.retptr(0);
236  les_ = g.lesptr(0);
237  tv_ = g.tvptr(0);
238  }
239 }
240 
262 template <typename T>
264  int tstp, herm_matrix_timestep_view<T> &g) {
265  assert(tstp == g.tstp_);
266 
267  tstp_ = g.tstp_;
268  ntau_ = g.ntau_;
269  size1_ = g.size1_;
270  size2_ = g.size1_;
271  element_size_ = size1_ * size1_;
272  sig_ = g.sig_;
273  // copy pointers
274  ret_ = g.ret_;
275  les_ = g.les_;
276  tv_ = g.tv_;
277  mat_ = g.mat_;
278 }
279 
304 template <typename T>
306  herm_matrix<T> &g, bool check_assert) {
307 
308  if(check_assert) assert(tstp>=-1 && tstp <=g.nt());
309  tstp_ = tstp;
310  ntau_ = g.ntau();
311  size1_ = g.size1();
312  size2_ = g.size2();
313  element_size_ = size1_ * size2_;
314  sig_ = g.sig();
315  if (tstp_ == -1) {
316  mat_ = g.matptr(0);
317  ret_ = 0;
318  les_ = 0;
319  tv_ = 0;
320  } else if (tstp_ >= 0) {
321  mat_ = 0;
322  ret_ = g.retptr(tstp_, 0);
323  les_ = g.lesptr(0, tstp_);
324  tv_ = g.tvptr(tstp_, 0);
325  }
326 }
327 
328 /* #######################################################################################
329 #
330 # WRITING ELEMENTS FROM ANY MATRIX TYPE
331 # OR FROM COMPLEX NUMBERS (then only the (0,0) element is addressed for dim>0)
332 #
333 ########################################################################################*/
334 #define herm_matrix_SET_ELEMENT_MATRIX \
335  { \
336  int r, s; \
337  for (r = 0; r < size1_; r++) \
338  for (s = 0; s < size2_; s++) \
339  x[r * size2_ + s] = M(r, s); \
340  }
341 
343 
363 template<typename T> template <class Matrix> void herm_matrix_timestep_view<T>::set_ret(int i, int j,Matrix &M){
364  assert(i == tstp_);
365  assert(j <= i);
366  cplx *x=retptr(j);
367  herm_matrix_SET_ELEMENT_MATRIX
368  }
369 
371 
391 template<typename T> template <class Matrix> void herm_matrix_timestep_view<T>::set_les(int i, int j, Matrix &M){
392  assert(j == tstp_);
393  assert(i <= j);
394  cplx *x=lesptr(i);
395  herm_matrix_SET_ELEMENT_MATRIX
396  }
397 
399 
419 template<typename T> template <class Matrix> void herm_matrix_timestep_view<T>::set_tv(int i, int j, Matrix &M){
420  assert(i == tstp_);
421  assert(j <= ntau_);
422  cplx *x=tvptr(j);
423  herm_matrix_SET_ELEMENT_MATRIX
424 }
425 
427 
445 template<typename T> template <class Matrix> void herm_matrix_timestep_view<T>::set_mat(int i,Matrix &M){
446  assert(i <= ntau_ );
447  cplx *x=matptr(i);
448  herm_matrix_SET_ELEMENT_MATRIX
449 }
450 
451 
452 
454 // the following routines are not very "efficent" but sometimes simple to
455 // implement
456 #define herm_matrix_READ_ELEMENT \
457 { \
458  int r, s, dim = size1_; \
459  M.resize(dim, dim); \
460  for (r = 0; r < dim; r++) \
461  for (s = 0; s < dim; s++) \
462  M(r, s) = x[r * dim + s]; \
463  }
464 #define herm_matrix_READ_ELEMENT_MINUS_CONJ \
465  { \
466  cplx w; \
467  int r, s, dim = size1_; \
468  M.resize(dim, dim); \
469  for (r = 0; r < dim; r++) \
470  for (s = 0; s < dim; s++) { \
471  w = x[s * dim + r]; \
472  M(r, s) = std::complex<T>(-w.real(), w.imag()); \
473  } \
474  }
475 
477 
497 template <typename T> template <class Matrix>
498  void herm_matrix_timestep_view<T>::get_les(int i, int j, Matrix &M){
499  assert(j == tstp_ && i <= tstp_);
500  cplx *x;
501  x = lesptr(i);
502  herm_matrix_READ_ELEMENT
503  }
504 
506 
526 template <typename T>
527 template <class Matrix>
528  void herm_matrix_timestep_view<T>::get_ret(int i, int j, Matrix &M) {
529  assert(i == tstp_ && j <= tstp_);
530  cplx *x;
531  x = retptr(j);
532  herm_matrix_READ_ELEMENT
533  }
534 
536 
556 template <typename T>
557 template <class Matrix>
558  void herm_matrix_timestep_view<T>::get_tv(int i, int j, Matrix &M) {
559  assert(i == tstp_);
560  cplx *x = tvptr(j);
561  herm_matrix_READ_ELEMENT
562  }
563 
565 
583 template <typename T>
584 template <class Matrix>
585  void herm_matrix_timestep_view<T>::get_mat(int i, Matrix &M) {
586  cplx *x = matptr(i);
587  herm_matrix_READ_ELEMENT
588  }
589 
591 
612 template <typename T>
613 template <class Matrix>
614  void herm_matrix_timestep_view<T>::get_matminus(int i, Matrix &M) {
615  cplx *x = matptr(ntau_ - i);
616  herm_matrix_READ_ELEMENT if (sig_ == -1) M = -M;
617  }
618 
619 
638 template <typename T>
640  assert(tstp == tstp_);
641  assert(tstp >= -1);
642  if (tstp == -1) {
643  memset(matptr(0), 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
644  } else {
645  memset(retptr(0), 0, sizeof(cplx) * (tstp + 1) * element_size_);
646  memset(tvptr(0), 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
647  memset(lesptr(0), 0, sizeof(cplx) * (tstp + 1) * element_size_);
648  }
649 }
650 
651 
674 template <typename T>
676  assert(tstp == tstp_);
677  assert(tstp >= -1 && tstp <= g1.nt() && "tstp >= -1 && tstp <= g1.nt()");
678  assert(g1.size1() == size1_ && "g1.size1() == size1_");
679  assert(g1.ntau() == ntau_ && "g1.ntau() == ntau_");
680  if (tstp == -1) {
681  memcpy(matptr(0), g1.matptr(0), sizeof(cplx) * (ntau_ + 1) * element_size_);
682  } else {
683  memcpy(retptr(0), g1.retptr(tstp, 0),
684  sizeof(cplx) * (tstp + 1) * element_size_);
685  memcpy(tvptr(0), g1.tvptr(tstp, 0),
686  sizeof(cplx) * (ntau_ + 1) * element_size_);
687  memcpy(lesptr(0), g1.lesptr(0, tstp),
688  sizeof(cplx) * (tstp + 1) * element_size_);
689  }
690 }
691 
714 template <typename T>
716  assert(tstp == tstp_);
717  assert(tstp == g1.tstp());
718  assert(tstp >= -1 && tstp <= g1.nt() && "tstp >= -1 && tstp <= g1.nt()");
719  assert(g1.size1() == size1_ && "g1.size1() == size1_");
720  assert(g1.ntau() == ntau_ && "g1.ntau() == ntau_");
721  if (tstp == -1) {
722  memcpy(matptr(0), g1.matptr(0), sizeof(cplx) * (ntau_ + 1) * element_size_);
723  } else {
724  memcpy(retptr(0), g1.retptr(0),
725  sizeof(cplx) * (tstp + 1) * element_size_);
726  memcpy(tvptr(0), g1.tvptr(0),
727  sizeof(cplx) * (ntau_ + 1) * element_size_);
728  memcpy(lesptr(0), g1.lesptr(0),
729  sizeof(cplx) * (tstp + 1) * element_size_);
730  }
731 }
732 
733 
734 
755 template <typename T>
757  assert(tstp == tstp_);
758  assert( ft.nt() >= tstp_);
759 
760  this->left_multiply(ft.ptr(-1), ft.ptr(0), weight);
761 }
762 
764 
788 template <typename T>
789 void herm_matrix_timestep_view<T>::left_multiply(int tstp, std::complex<T> *f0,
790  std::complex<T> *ft, T weight) {
791  int m;
792  cplx *xtemp, *ftemp, *x0;
793  xtemp = new cplx[element_size_];
794  assert(tstp >= -1 && "tstp >= -1");
795  if (tstp == -1) {
796  x0 = matptr(0);
797  for (m = 0; m <= ntau_; m++) {
798  element_mult<T, LARGESIZE>(size1_, xtemp, f0,
799  x0 + m * element_size_);
800  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
801  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
802  }
803  } else {
804  ftemp = ft + tstp * element_size_;
805  x0 = retptr(tstp, 0);
806  for (m = 0; m <= tstp; m++) {
807  element_mult<T, LARGESIZE>(size1_, xtemp, ftemp,
808  x0 + m * element_size_);
809  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
810  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
811  }
812  x0 = tvptr(tstp, 0);
813  for (m = 0; m <= ntau_; m++) {
814  element_mult<T, LARGESIZE>(size1_, xtemp, ftemp,
815  x0 + m * element_size_);
816  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
817  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
818  }
819  x0 = lesptr(0, tstp);
820  for (m = 0; m <= tstp; m++) {
821  element_mult<T, LARGESIZE>(size1_, xtemp, ft + m * element_size_,
822  x0 + m * element_size_);
823  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
824  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
825  }
826  }
827  delete[] xtemp;
828 }
829 
830 
851 template <typename T>
853  assert(tstp == tstp_);
854  assert( ft.nt() >= tstp_);
855 
856  this->right_multiply(ft.ptr(-1), ft.ptr(0), weight);
857 }
858 
860 
884 template <typename T>
885 void herm_matrix_timestep_view<T>::right_multiply(int tstp, std::complex<T> *f0,
886  std::complex<T> *ft, T weight) {
887  int m;
888  cplx *xtemp, *ftemp, *x0;
889  xtemp = new cplx[element_size_];
890  assert(tstp >= -1 && "tstp >= -1");
891  if (tstp == -1) {
892  x0 = matptr(0);
893  for (m = 0; m <= ntau_; m++) {
894  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
895  f0);
896  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
897  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
898  }
899  } else {
900  x0 = retptr(tstp, 0);
901  for (m = 0; m <= tstp; m++) {
902  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
903  ft + m * element_size_);
904  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
905  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
906  }
907  x0 = tvptr(tstp, 0);
908  for (m = 0; m <= ntau_; m++) {
909  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
910  f0);
911  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
912  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
913  }
914  ftemp = ft + tstp * element_size_;
915  x0 = lesptr(0, tstp);
916  for (m = 0; m <= tstp; m++) {
917  element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
918  ftemp);
919  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
920  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
921  }
922  }
923  delete[] xtemp;
924 }
925 
926 
927 
928 
929 
930 
952 template <typename T>
954  T weight) {
955  int m;
956  cplx *xtemp, *ftemp, *x0, *f0, *fcc;
957  xtemp = new cplx[element_size_];
958  fcc = new cplx[element_size_];
959  assert(tstp >= -1 && ft.nt() >= tstp &&
960  ft.size1() == size1_ && ft.size2() == size2_ &&
961  "tstp >= -1 && ft.nt() >= tstp && ft.size1() == size1_ && ft.size2() == size2_");
962 
963  f0 = ft.ptr(-1);
964  element_conj<T, LARGESIZE>(size1_, fcc, f0);
965  if (tstp == -1) {
966  x0 = matptr(0);
967  for (m = 0; m <= ntau_; m++) {
968  element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
969  x0 + m * element_size_);
970  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
971  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
972  }
973  } else {
974  ftemp = ft.ptr(tstp);
975  element_conj<T, LARGESIZE>(size1_, fcc, ftemp);
976  x0 = retptr(tstp, 0);
977  for (m = 0; m <= tstp; m++) {
978  element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
979  x0 + m * element_size_);
980  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
981  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
982  }
983  x0 = tvptr(tstp, 0);
984  for (m = 0; m <= ntau_; m++) {
985  element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
986  x0 + m * element_size_);
987  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
988  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
989  }
990  x0 = lesptr(0, tstp);
991  for (m = 0; m <= tstp; m++) {
992  element_conj<T, LARGESIZE>(size1_, fcc, ft.ptr(m));
993  element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
994  x0 + m * element_size_);
995  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
996  element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
997  }
998  }
999  delete[] xtemp;
1000  delete[] fcc;
1001 }
1002 
1003 
1004 
1005 
1006 
1007 
1037 template <typename T>
1038 void herm_matrix_timestep_view<T>::set_to_data(CPLX *data, int tstp, int ntau,
1039  int size, int sig) {
1040  tstp_ = tstp;
1041  ntau_ = ntau;
1042  size1_ = size;
1043  size2_ = size;
1044  element_size_ = size1_ * size2_;
1045  sig_ = sig;
1046  if (tstp == -1) {
1047  mat_ = data;
1048  ret_ = 0;
1049  les_ = 0;
1050  tv_ = 0;
1051  } else if (tstp >= 0) {
1052  mat_ = 0;
1053  ret_ = data;
1054  tv_ = data + (tstp_ + 1) * element_size_;
1055  les_ = data + (tstp_ + 1 + ntau_ + 1) * element_size_;
1056  }
1057 }
1058 
1082 template <typename T>
1084  assert(tstp>=-1 && tstp<=g.nt());
1085  tstp_ = tstp;
1086  ntau_ = g.ntau();
1087  size1_ = g.size1();
1088  size2_ = g.size2();
1089  element_size_ = size1_ * size2_;
1090  sig_ = g.sig();
1091  if (tstp == -1) {
1092  mat_ = g.matptr(0);
1093  ret_ = 0;
1094  les_ = 0;
1095  tv_ = 0;
1096  } else if (tstp >= 0) {
1097  mat_ = 0;
1098  ret_ = g.retptr(tstp, 0);
1099  tv_ = g.tvptr(tstp, 0);
1100  les_ = g.lesptr(0, tstp);
1101  }
1102 }
1103 
1125 template <typename T>
1127  tstp_ = g.tstp_;
1128  ntau_ = g.ntau();
1129  size1_ = g.size1();
1130  size2_ = g.size2();
1131  element_size_ = size1_ * size2_;
1132  sig_ = g.sig();
1133  if (tstp_ == -1) {
1134  mat_ = g.matptr(0);
1135  ret_ = 0;
1136  les_ = 0;
1137  tv_ = 0;
1138  } else if (tstp_ >= 0) {
1139  mat_ = 0;
1140  ret_ = g.retptr(0);
1141  tv_ = g.tvptr(0);
1142  les_ = g.lesptr(0);
1143  }
1144 }
1145 
1169 template <typename T>
1170 void herm_matrix_timestep_view<T>::get_data(CPLX *ret, CPLX *les, CPLX *tv,
1171  CPLX *mat) {
1172  if (tstp_ == -1) {
1173  memcpy(mat_, mat, sizeof(CPLX) * (ntau_ + 1) * element_size_);
1174  ret_ = 0;
1175  les_ = 0;
1176  tv_ = 0;
1177  } else if (tstp_ >= 0) {
1178  mat_ = 0;
1179  memcpy(ret_, ret, sizeof(CPLX) * (tstp_ + 1) * element_size_);
1180  memcpy(les_, les, sizeof(CPLX) * (tstp_ + 1) * element_size_);
1181  memcpy(tv_, tv, sizeof(CPLX) * (ntau_ + 1) * element_size_);
1182  }
1183 }
1184 
1203 template <typename T>
1205  assert(tstp_==g.tstp_ && size1_ == g.size1() && size2_ == g.size2() && ntau_ == g.ntau());
1206  if (tstp_ == -1)
1207  get_data(NULL, NULL, NULL, g.mat_);
1208  else
1209  get_data(g.ret_, g.les_, g.tv_, NULL);
1210 }
1211 
1212 
1231 template <typename T>
1232 template <class GG>
1234  herm_matrix_timestep_view<T> tmp(tstp_, g);
1235  get_data(tmp);
1236 }
1237 
1257 template <typename T>
1259  herm_matrix_timestep<T> &timestep) {
1260  int len = (2 * (tstp + 1) + ntau_ + 1) * element_size_;
1261  cplx *x;
1262  assert(tstp_ == tstp);
1263  if (timestep.total_size_ < len)
1264  timestep.resize(tstp, ntau_, size1_);
1265  x = timestep.data_;
1266  timestep.tstp_ = tstp;
1267  timestep.ntau_ = ntau_;
1268  timestep.size1_ = size1_;
1269  if (tstp == -1) {
1270  memcpy(x, mat_, sizeof(cplx) * (ntau_ + 1) * element_size_);
1271  } else {
1272  memcpy(x, ret_, sizeof(cplx) * (tstp + 1) * element_size_);
1273  memcpy(x + (tstp + 1) * element_size_, tv_,
1274  sizeof(cplx) * (ntau_ + 1) * element_size_);
1275  memcpy(x + (tstp + 1 + ntau_ + 1) * element_size_, les_,
1276  sizeof(cplx) * (tstp + 1) * element_size_);
1277  }
1278 }
1279 
1280 
1308 template <typename T>
1310  int i1, int i2, herm_matrix_timestep_view<T> &g, int j1, int j2) {
1311  int i, sij = i1 * size2_ + i2, tij = j1 * g.size2() + j2;
1312  assert(tstp_==g.tstp_ && ntau_ == g.ntau_);
1313 
1314  assert(i1>=0 && i1 <=size1_ -1);
1315  assert(i2>=0 && i2 <=size2_ -1);
1316  assert(j1>=0 && j1 <=g.size1_-1);
1317  assert(j2>=0 && j2 <=g.size2_-1);
1318 
1319  if (tstp_ == -1) {
1320  for (i = 0; i <= ntau_; i++)
1321  matptr(i)[sij] = g.matptr(i)[tij];
1322  } else {
1323  for (i = 0; i <= tstp_; i++)
1324  retptr(i)[sij] = g.retptr(i)[tij];
1325  for (i = 0; i <= ntau_; i++)
1326  tvptr(i)[sij] = g.tvptr(i)[tij];
1327  for (i = 0; i <= tstp_; i++)
1328  lesptr(i)[sij] = g.lesptr(i)[tij];
1329  }
1330 }
1331 
1332 
1360 template <typename T>
1361 template <class GG>
1363  int j1, int j2) {
1364  herm_matrix_timestep_view<T> tmp(tstp_, g);
1365  set_matrixelement(i1, i2, tmp, j1, j2);
1366 }
1367 
1368 
1369 #define HERM_MATRIX_INCR_TSTP \
1370  if (alpha == CPLX(1.0, 0.0)) { \
1371  for (i = 0; i < len; i++) \
1372  x0[i] += x[i]; \
1373  } else { \
1374  for (i = 0; i < len; i++) \
1375  x0[i] += alpha * x[i]; \
1376  }
1377 
1378 
1401 template <typename T>
1402 void
1404  CPLX alpha) {
1405  assert(tstp_==g.tstp_ && ntau_ ==g.ntau_ && size1_ == g.size1_ && size2_ == g.size2_);
1406 
1407  int i, len;
1408  CPLX *x, *x0;
1409  if (tstp_ == -1) {
1410  len = (ntau_ + 1) * element_size_;
1411  x0 = mat_;
1412  x = g.mat_;
1413  HERM_MATRIX_INCR_TSTP
1414  } else {
1415  len = (tstp_ + 1) * element_size_;
1416  x0 = ret_;
1417  x = g.ret_;
1418  HERM_MATRIX_INCR_TSTP
1419  x0 = les_;
1420  x = g.les_;
1421  HERM_MATRIX_INCR_TSTP
1422  len = (ntau_ + 1) * element_size_;
1423  x0 = tv_;
1424  x = g.tv_;
1425  HERM_MATRIX_INCR_TSTP
1426  }
1427 }
1428 #undef HERM_MATRIX_INCR_TSTP
1429 
1453 template <typename T>
1454 template <class GG>
1456  herm_matrix_timestep_view<T> tmp(tstp_, g);
1457  incr_timestep(tmp, alpha);
1458 }
1459 #define HERM_MATRIX_INCR_TSTP \
1460  if (alpha == 1.0) { \
1461  for (i = 0; i < len; i++) \
1462  x0[i] += x[i]; \
1463  } else { \
1464  for (i = 0; i < len; i++) \
1465  x0[i] += alpha * x[i]; \
1466  }
1467 
1468 
1491 template <typename T>
1492 void
1494  T alpha) {
1495  assert(tstp_==g.tstp_ && ntau_ == g.ntau_ && size1_ == g.size1_ && size2_ == g.size2_);
1496 
1497  int i, len;
1498  CPLX *x, *x0;
1499  if (tstp_ == -1) {
1500  len = (ntau_ + 1) * element_size_;
1501  x0 = mat_;
1502  x = g.mat_;
1503  HERM_MATRIX_INCR_TSTP
1504  } else {
1505  len = (tstp_ + 1) * element_size_;
1506  x0 = ret_;
1507  x = g.ret_;
1508  HERM_MATRIX_INCR_TSTP
1509  x0 = les_;
1510  x = g.les_;
1511  HERM_MATRIX_INCR_TSTP
1512  len = (ntau_ + 1) * element_size_;
1513  x0 = tv_;
1514  x = g.tv_;
1515  HERM_MATRIX_INCR_TSTP
1516  }
1517 }
1518 #undef HERM_MATRIX_INCR_TSTP
1519 
1542 template <typename T>
1543 template <class GG>
1545  herm_matrix_timestep_view<T> tmp(tstp_, g);
1546  incr_timestep(tmp, alpha);
1547 }
1548 
1568 template <typename T>
1570  int m;
1571  CPLX *x0;
1572  if (tstp_ == -1) {
1573  x0 = mat_;
1574  for (m = 0; m <= ntau_; m++) {
1575  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
1576  weight);
1577  }
1578  } else {
1579  x0 = ret_;
1580  for (m = 0; m <= tstp_; m++) {
1581  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
1582  weight);
1583  }
1584  x0 = tv_;
1585  for (m = 0; m <= ntau_; m++) {
1586  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
1587  weight);
1588  }
1589  x0 = les_;
1590  for (m = 0; m <= tstp_; m++) {
1591  element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_,
1592  weight);
1593  }
1594  }
1595 }
1597 
1598 #if CNTR_USE_MPI == 1
1599 
1601 template <typename T>
1602 void my_mpi_reduce(std::complex<T> *data, int len, int root) {
1603  std::cerr << __PRETTY_FUNCTION__ << ", LEN=" << len
1604  << " ... NOT DEFINED FOR THIS TYPE " << std::endl;
1605  exit(0);
1606  }
1607 
1608 
1610 template<>
1611 inline void my_mpi_reduce<double>(std::complex<double> *data, int len, int root) {
1612  int tid, ntasks;
1613  MPI_Comm_rank(MPI_COMM_WORLD, &tid);
1614  MPI_Comm_size(MPI_COMM_WORLD, &ntasks);
1615  assert(root>=0 && root <= ntasks -1);
1616  assert(len>=0);
1617  if (tid == root) {
1618  MPI_Reduce(MPI_IN_PLACE, (double *)data, 2 * len,
1619  MPI_DOUBLE, MPI_SUM, root, MPI_COMM_WORLD);
1620  } else {
1621  MPI_Reduce((double *)data, (double *)data, 2 * len,
1622  MPI_DOUBLE, MPI_SUM, root, MPI_COMM_WORLD);
1623  }
1624 }
1625 
1627 
1646 template <typename T>
1647 void herm_matrix_timestep_view<T>::MPI_Reduce(int root) {
1648  if (tstp_ == -1) {
1649  my_mpi_reduce<T>(mat_, (ntau_ + 1) * element_size_, root);
1650  } else {
1651  my_mpi_reduce<T>(les_, (tstp_ + 1) * element_size_, root);
1652  my_mpi_reduce<T>(ret_, (tstp_ + 1) * element_size_, root);
1653  my_mpi_reduce<T>(tv_, (ntau_ + 1) * element_size_, root);
1654  }
1655 }
1656 
1658 
1679 template <typename T>
1681  assert(tstp == tstp_);
1682  if (tstp_ == -1) {
1683  my_mpi_reduce<T>(mat_, (ntau_ + 1) * element_size_, root);
1684  } else {
1685  my_mpi_reduce<T>(les_, (tstp_ + 1) * element_size_, root);
1686  my_mpi_reduce<T>(ret_, (tstp_ + 1) * element_size_, root);
1687  my_mpi_reduce<T>(tv_, (ntau_ + 1) * element_size_, root);
1688  }
1689 }
1690 
1710 template <typename T>
1712  int numtasks,taskid;
1713  MPI_Comm_size(MPI_COMM_WORLD, &numtasks);
1714  MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
1715  // test effective on root:
1716  assert(tstp == tstp_);
1717  if (sizeof(T) == sizeof(double)){
1718  if (tstp_ == -1){
1719  MPI_Bcast(mat_, (ntau_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, root, MPI_COMM_WORLD);
1720  } else {
1721  MPI_Bcast(les_, (tstp_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, root, MPI_COMM_WORLD);
1722  MPI_Bcast(ret_, (tstp_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, root, MPI_COMM_WORLD);
1723  MPI_Bcast(tv_, (ntau_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, root, MPI_COMM_WORLD);
1724  }
1725  } else { // assuming single precision
1726  if (tstp_ == -1){
1727  MPI_Bcast(mat_, (ntau_ + 1) * element_size_, MPI_COMPLEX, root, MPI_COMM_WORLD);
1728  } else {
1729  MPI_Bcast(les_, (tstp_ + 1) * element_size_, MPI_COMPLEX, root, MPI_COMM_WORLD);
1730  MPI_Bcast(ret_, (tstp_ + 1) * element_size_, MPI_COMPLEX, root, MPI_COMM_WORLD);
1731  MPI_Bcast(tv_, (ntau_ + 1) * element_size_, MPI_COMPLEX, root, MPI_COMM_WORLD);
1732  }
1733  }
1734 
1735 }
1736 
1756 template <typename T>
1757 void herm_matrix_timestep_view<T>::Send_timestep(int tstp, int dest, int tag) {
1758  int taskid;
1759  MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
1760  assert(tstp == tstp_);
1761  if (!(taskid == dest)) {
1762  if (sizeof(T) == sizeof(double)){
1763  if (tstp_ == -1){
1764  MPI_Send(mat_, (ntau_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, dest, tag, MPI_COMM_WORLD);
1765  } else {
1766  MPI_Send(les_, (tstp_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, dest, tag, MPI_COMM_WORLD);
1767  MPI_Send(ret_, (tstp_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, dest, tag, MPI_COMM_WORLD);
1768  MPI_Send(tv_, (ntau_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, dest, tag, MPI_COMM_WORLD);
1769  }
1770  }
1771  else {
1772  if (tstp_ == -1){
1773  MPI_Send(mat_, (ntau_ + 1) * element_size_, MPI_COMPLEX, dest, tag, MPI_COMM_WORLD);
1774  } else {
1775  MPI_Send(les_, (tstp_ + 1) * element_size_, MPI_COMPLEX, dest, tag, MPI_COMM_WORLD);
1776  MPI_Send(ret_, (tstp_ + 1) * element_size_, MPI_COMPLEX, dest, tag, MPI_COMM_WORLD);
1777  MPI_Send(tv_, (ntau_ + 1) * element_size_, MPI_COMPLEX, dest, tag, MPI_COMM_WORLD);
1778  }
1779  }
1780  } else {
1781  // do nothing
1782  }
1783 }
1784 
1804 template <typename T>
1805 void herm_matrix_timestep_view<T>::Recv_timestep(int tstp, int root, int tag) {
1806  int taskid;
1807  MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
1808  assert(tstp == tstp_);
1809  if (!(taskid == root)) {
1810  if (sizeof(T) == sizeof(double))
1811  if (tstp_ == -1){
1812  MPI_Recv(mat_, (ntau_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, root, tag, MPI_COMM_WORLD);
1813  } else {
1814  MPI_Recv(les_, (tstp_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, root, tag, MPI_COMM_WORLD);
1815  MPI_Recv(ret_, (tstp_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, root, tag, MPI_COMM_WORLD);
1816  MPI_Recv(tv_, (ntau_ + 1) * element_size_, MPI_DOUBLE_COMPLEX, root, tag, MPI_COMM_WORLD);
1817  }
1818  else{
1819  if (tstp_ == -1){
1820  MPI_Recv(mat_, (ntau_ + 1) * element_size_, MPI_COMPLEX, root, tag, MPI_COMM_WORLD);
1821  } else {
1822  MPI_Recv(les_, (tstp_ + 1) * element_size_, MPI_COMPLEX, root, tag, MPI_COMM_WORLD);
1823  MPI_Recv(ret_, (tstp_ + 1) * element_size_, MPI_COMPLEX, root, tag, MPI_COMM_WORLD);
1824  MPI_Recv(tv_, (ntau_ + 1) * element_size_, MPI_COMPLEX, root, tag, MPI_COMM_WORLD);
1825  }
1826  }
1827  }
1828 }
1829 
1830 
1831 #endif
1832 
1852 #if CNTR_USE_HDF5 == 1
1853 template <typename T>
1855  store_int_attribute_to_hid(group_id, std::string("ntau"), ntau_);
1856  store_int_attribute_to_hid(group_id, std::string("tstp"), tstp_);
1857  store_int_attribute_to_hid(group_id, std::string("sig"), sig_);
1858  store_int_attribute_to_hid(group_id, std::string("size1"), size1_);
1859  store_int_attribute_to_hid(group_id, std::string("size2"), size2_);
1860  store_int_attribute_to_hid(group_id, std::string("element_size"),
1861  element_size_);
1862  hsize_t len_shape = 3, shape[3];
1863  shape[1] = size1_;
1864  shape[2] = size2_;
1865  if (tstp_ == -1) {
1866  shape[0] = ntau_ + 1;
1867  store_cplx_array_to_hid(group_id, std::string("mat"), matptr(0),
1868  shape, len_shape);
1869  } else if (tstp_ > -1) {
1870  shape[0] = (tstp_ + 1);
1871  // CHECK: implement store_cplx_array_to_hid with template typename T
1872  store_cplx_array_to_hid(group_id, std::string("ret"), retptr(0),
1873  shape, len_shape);
1874  store_cplx_array_to_hid(group_id, std::string("les"), lesptr(0),
1875  shape, len_shape);
1876  shape[0] = (ntau_ + 1);
1877  store_cplx_array_to_hid(group_id, std::string("tv"), tvptr(0), shape,
1878  len_shape);
1879  }
1880 }
1881 
1903 template <typename T>
1905  const char *groupname) {
1906  hid_t sub_group_id = create_group(group_id, groupname);
1907  this->write_to_hdf5(sub_group_id);
1908  close_group(sub_group_id);
1909 }
1910 
1932 template <typename T>
1934  const char *groupname) {
1935  hid_t file_id = open_hdf5_file(filename);
1936  this->write_to_hdf5(file_id, groupname);
1937  close_hdf5_file(file_id);
1938 }
1939 
1959 template <typename T>
1961  // -- Read data: NO RESIZE POSSIBLE
1962  int tstp = read_primitive_type<int>(group_id, "tstp");
1963  int ntau = read_primitive_type<int>(group_id, "ntau");
1964  int sig = read_primitive_type<int>(group_id, "sig");
1965  int size1 = read_primitive_type<int>(group_id, "size1");
1966  int size2 = read_primitive_type<int>(group_id, "size2");
1967  int element_size = read_primitive_type<int>(group_id, "element_size");
1968 
1969  assert(element_size == element_size_ && size1 == size1_ && size2 == size2_ && tstp == tstp_ && ntau ==ntau_);
1970 
1971  sig_ = sig;
1972  if (tstp == -1) {
1973  hsize_t mat_size = (ntau + 1) * element_size;
1974  read_primitive_type_array(group_id, "mat", mat_size, matptr(0));
1975  } else if (tstp_ > -1) {
1976  hsize_t ret_size = (tstp_ + 1) * element_size;
1977  read_primitive_type_array(group_id, "ret", ret_size, retptr(0));
1978  hsize_t les_size = (tstp_ + 1) * element_size;
1979  read_primitive_type_array(group_id, "les", les_size, lesptr(0));
1980  hsize_t tv_size = (ntau_ + 1) * element_size;
1981  read_primitive_type_array(group_id, "tv", tv_size, tvptr(0));
1982  }
1983 }
1984 
2006 template <typename T>
2008  const char *groupname) {
2009  hid_t sub_group_id = open_group(group_id, groupname);
2010  this->read_from_hdf5(sub_group_id);
2011  close_group(sub_group_id);
2012 }
2013 
2034 template <typename T>
2036  const char *groupname) {
2037  hid_t file_id = read_hdf5_file(filename);
2038  this->read_from_hdf5(file_id, groupname);
2039  close_hdf5_file(file_id);
2040 }
2041 #endif
2042 
2044 
2061 template <typename T>
2063  assert(tstp==tstp_);
2064 
2065  CPLX x1;
2066  if (tstp_ == -1) {
2067  x1 = *matptr(ntau_);
2068  return -x1;
2069  } else {
2070  x1 = *lesptr(tstp_);
2071  return CPLX(0.0, sig_) * x1;
2072  }
2073 }
2074 
2075 
2095 template <typename T>
2096 void herm_matrix_timestep_view<T>::density_matrix(int tstp, std::complex<T> &M) {
2097  assert(tstp==tstp_);
2098  if (tstp_ == -1) {
2099  M = -(*matptr(ntau_));
2100  } else {
2101  M = CPLX(0.0, sig_)*(*lesptr(tstp_));
2102  }
2103 }
2104 
2124 template <typename T>
2125 template <class Matrix>
2127  assert(tstp == tstp_);
2128 
2129  CPLX *x;
2130  if (tstp == -1) {
2131  x = matptr(ntau_);
2132  herm_matrix_READ_ELEMENT M *= (-1.0);
2133  } else {
2134  x = lesptr(tstp);
2135  herm_matrix_READ_ELEMENT M *= CPLX(0.0, 1.0 * sig_);
2136  }
2137 }
2138 
2139 #undef herm_matrix_READ_ELEMENT
2140 #undef herm_matrix_READ_ELEMENT_MINUS_CONJ
2141 #undef CPLX
2142 
2143 } // namespace cntr
2144 
2145 #endif // CNTR_HERM_TIMESTEP_VIEW_IMPL_H
void incr_timestep(herm_matrix_timestep_view< T > &g, std::complex< T > alpha)
Class herm_matrix_timestep_view serves for interfacing with class herm_matrix_timestep without copyi...
int size1(void) const
Class herm_matrix_timestep deals with contour objects at a particular timestep .
void smul(T alpha)
Multiply herm_matrix_timestep_view with scalar weight.
int size2(void) const
void Recv_timestep(int tstp, int root, int tag)
Recevies the herm_matrix_timestep_view at a given time step from a specific task.
void resize(int nt, int ntau, int size1)
Resizes herm_matrix_timestep object with respect to the number of points on the Matsubara branch and...
std::complex< double > cplx
Definition: fourier.cpp:11
herm_matrix_timestep_view & operator=(const herm_matrix_timestep_view &g)
Copy assignment operator for herm_matrix_timestep_view g class.
void set_matrixelement(int i1, int i2, herm_matrix_timestep_view< T > &g, int j1, int j2)
Set the matrix element of for each component from
Class function for objects with time on real axis.
void set_timestep(int tstp, herm_matrix< T > &g1)
Sets all components of herm_matrix_timestep_view to the components of a given herm_matrix at time st...
void get_timestep(int tstp, herm_matrix_timestep< T > &timestep)
Get timestep into herm_matrix_timestep
void get_data(herm_matrix_timestep_view< T > &g)
Return the data into herm_matrix_timestep_view
void Reduce_timestep(int tstp, int root)
MPI reduce for the herm_matrix_timestep_view.
void read_from_hdf5(hid_t group_id, const char *groupname)
Read herm_matrix_timestep_view from hdf5 group given by group ID group_id and group name groupname...
int nt(void) const
Class herm_matrix for two-time contour objects with hermitian symmetry.
void Send_timestep(int tstp, int dest, int tag)
Sends the herm_matrix_timestep_view at a given time step to a specific task.
void set_timestep_zero(int tstp)
Sets all components at time step tstp to zero.
void Bcast_timestep(int tstp, int root)
Broadcasts the herm_matrix_timestep_view at a given time step to all ranks.
void write_to_hdf5(hid_t group_id)
Write herm_matrix_timestep_view to hdf5 group given by group_id.
void set_to_data(std::complex< T > *data, int tstp, int ntau, int size, int sig)
void left_multiply_hermconj(int tstp, function< T > &ft, T weight=1.0)
Left-multiplies the herm_matrix_timestep_view with the hermitian conjugate of a contour function...