NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
cntr_function_impl.hpp
Go to the documentation of this file.
1 #ifndef CNTR_FUNCTION_IMPL_H
2 #define CNTR_FUNCTION_IMPL_H
3 
4 #include "cntr_function_decl.hpp"
5 #include "linalg.hpp"
6 //#include "cntr_exception.hpp"
7 #include "cntr_elements.hpp"
8 
9 namespace cntr {
10 
11 template <typename T>
13  data_ = 0;
14  nt_ = -2;
15  size1_ = 0;
16  size2_ = 0;
17  element_size_ = 0;
18 }
19 template <typename T>
21  if (data_ != 0)
22  delete[] data_;
23 }
44 template <typename T>
45 function<T>::function(int nt, int size1) {
46  int len = (nt + 2) * size1 * size1;
47  assert(size1 >= 0 && nt >= -1);
48  if (len == 0)
49  data_ = 0;
50  else {
51  data_ = new cplx[len];
52  }
53  size1_ = size1;
54  size2_ = size1;
55  element_size_ = size1 * size1_;
56  nt_ = nt;
57  total_size_ = (nt_ + 2) * size1_ * size2_;
58 }
59 
81 template <typename T> function<T>::function(int nt,int size1,int size2){
82  int len=(nt+2)*size1*size2;
83  assert(size1>=0 && nt>=-1);
84  if(len==0) data_=0;
85  else{
86  data_ = new cplx [len];
87  memset(data_, 0, sizeof(cplx)*len);
88  }
89  size1_=size1;
90  size2_=size2;
91  element_size_=size1*size2_;
92  nt_=nt;
93  total_size_ = len;
94 }
112 template <typename T>
113 function<T>::function(const function &g) {
114  int len;
115  nt_ = g.nt_;
116  size1_ = g.size1_;
117  size2_ = g.size2_;
118  element_size_ = size1_ * size2_;
119  len = (nt_ + 2) * size1_ * size2_;
120  if (len > 0) {
121  total_size_ = len;
122  data_ = new cplx[len];
123  memcpy(data_, g.data_, sizeof(cplx) * len);
124  } else {
125  data_ = 0;
126  }
127 }
128 #if __cplusplus >= 201103L
129 
147 template <typename T>
148 function<T>::function(function &&g) noexcept
149  : data_(g.data_),
150  nt_(g.nt_),
151  size1_(g.size1_),
152  size2_(g.size2_),
153  element_size_(g.element_size_) {
154  g.data_ = nullptr;
155  g.nt_ = -2;
156  g.size1_ = 0;
157  g.size2_ = 0;
158  g.element_size_ = 0;
159 }
178 template <typename T>
179 function<T> &function<T>::operator=(function &&g) noexcept {
180  if (&g == this)
181  return *this;
182 
183  data_ = g.data_;
184  nt_ = g.nt_;
185  size1_ = g.size1_;
186  size2_ = g.size2_;
187  element_size_ = g.element_size_;
188 
189  g.data_ = nullptr;
190  g.nt_ = -2;
191  g.size1_ = 0;
192  g.size2_ = 0;
193  g.element_size_ = 0;
194 
195  return *this;
196 }
197 #endif
198 
215 template <typename T>
216 function<T> &function<T>::operator=(const function &g) {
217  int len;
218  if (this == &g)
219  return *this;
220  if (data_ != 0)
221  delete[] data_;
222  nt_ = g.nt_;
223  size1_ = g.size1_;
224  size2_ = g.size2_;
225  element_size_ = size1_ * size2_;
226  len = (nt_ + 2) * size1_ * size2_;
227  if (len > 0) {
228  total_size_ = len;
229  data_ = new cplx[len];
230  memcpy(data_, g.data_, sizeof(cplx) * len);
231  } else {
232  data_ = 0;
233  }
234  return *this;
235 }
254 template <typename T>
255 void function<T>::resize(int nt, int size1) {
256  int len = (nt + 2) * size1 * size1;
257  assert(nt >= -1 && size1 >= 0);
258  if (data_ != 0)
259  delete[] data_;
260  if (len == 0)
261  data_ = 0;
262  else {
263  data_ = new cplx[len];
264  }
265  size1_ = size1;
266  size2_ = size1;
267  element_size_ = size1_ * size1_;
268  nt_ = nt;
269  total_size_ = len;
270 }
285 template <typename T>
287  int len = element_size_ * (nt_ + 2);
288  if (len > 0)
289  memset(data_, 0, sizeof(cplx) * len);
290 }
307 template <typename T>
308 void function<T>::set_constant(std::complex<T> *f0) {
309  int t;
310  for (t = -1; t <= nt_; t++)
311  element_set<T, LARGESIZE>(size1_, ptr(t), f0);
312 }
328 template <typename T> template<class EigenMatrix>
329 void function<T>::set_constant(EigenMatrix &M){
330  assert(M.rows() == size1_);
331  int t;
332  for(t=-1;t<=nt_;t++){
333  set_value(t,M);
334  }
335 }
336 
337 //template <class EigenMatrix>
355 template <typename T> template<class EigenMatrix>
356 void function<T>::set_value(int tstp, cplx x) {
357  int i1, i2;
358  cplx *ft;
359  assert(1 == size1_);
360  assert(1 == size2_);
361  *ptr(tstp) =x;
362 }
381 template <typename T>
382 template <class EigenMatrix>
383 void function<T>::set_value(int tstp, EigenMatrix &M) {
384  int i1, i2;
385  cplx *ft;
386  assert(M.rows() == size1_);
387  assert(M.cols() == size2_);
388 
389  ft = ptr(tstp);
390  for (i1 = 0; i1 < size1_; i1++)
391  for (i2 = 0; i2 < size2_; i2++)
392  ft[i1 * size2_ + i2] = M(i1, i2);
393 }
412 template <typename T>
413 template <class EigenMatrix>
414 void function<T>::get_value(int tstp, EigenMatrix &M) const {
415  int i1, i2;
416  const cplx *ft;
417  M.resize(size1_, size2_);
418  ft = ptr(tstp);
419  for (i1 = 0; i1 < size1_; i1++)
420  for (i2 = 0; i2 < size2_; i2++)
421  M(i1, i2) = ft[i1 * size2_ + i2];
422 }
439 template <typename T>
440 void function<T>::smul(T weight)
441 {
442  for(int m = 0; m < total_size_; m++)
443  data_[m] *= weight;
444 }
463 template <typename T>
465 {
466  cplx *mytemp, *xtemp, *ftemp;
467  xtemp = new cplx[element_size_];
468  for (int tstp = -1; tstp <= nt_; tstp++) {
469  ftemp = f + (tstp + 1) * element_size_;
470  mytemp = ptr(tstp);
471  element_mult<T, LARGESIZE>(size1_, xtemp,
472  ftemp, mytemp);
473  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
474  element_set<T, LARGESIZE>(size1_, mytemp, xtemp);
475  }
476  delete[] xtemp;
477 }
497 template <typename T>
499 {
500  cplx *mytemp, *xtemp, *ftemp;
501  xtemp = new cplx[element_size_];
502  for (int tstp = -1; tstp <= nt_; tstp++) {
503  ftemp = f + (tstp + 1) * element_size_;
504  mytemp = ptr(tstp);
505  element_mult<T, LARGESIZE>(size1_, xtemp,
506  mytemp, ftemp);
507  element_smul<T, LARGESIZE>(size1_, xtemp, weight);
508  element_set<T, LARGESIZE>(size1_, mytemp, xtemp);
509  }
510  delete[] xtemp;
511 }
530 template <typename T>
532  assert(this -> nt_ == f.nt_);
533  assert(this -> size1_ == f.size1_);
534  assert(this -> size2_ == f.size2_);
535  this -> left_multiply(f.ptr(-1), weight);
536 }
555 template <typename T>
557  assert(this -> nt_ == f.nt_);
558  assert(this -> size1_ == f.size1_);
559  assert(this -> size2_ == f.size2_);
560  this->right_multiply(f.ptr(-1), weight);
561 }
589 template <typename T> template<class EigenMatrix>
590 void function<T>::set_matrixelement(int tstp,int i1,int i2,EigenMatrix &M,int j1,int j2){
591  assert(0<=i1 && i1<size1_ && 0<=i2 && i2<size2_);
592  assert(0<=j1 && j1<M.rows() && 0<=j2 && j2<M.cols());
593  cplx *ft;
594  ft=ptr(tstp);
595  ft[i1*size2_+i2]=M(j1,j2);
596 }
619 template <typename T>
621 {
622  assert(this -> nt_ == g.nt_);
623  assert(i1 <= this -> size1_);
624  assert(i2 <= this -> size2_);
625  cplx *target_position, *my_data;
626  for(int tstp = -1; tstp <= nt_; tstp++)
627  {
628  my_data = this -> ptr(tstp) + i1 * size2_ + i2;
629  target_position = g.ptr(tstp);
630  element_set <T, 1> (1, target_position, my_data);
631  //*target_position = *my_data;
632  }
633 }
651 template <typename T>
652 void function<T>::incr(function<T> &g, T weight)
653 {
654  assert(this -> nt_ == g.nt_);
655  for(int m = 0; m < total_size_; m++)
656  {
657  this -> data_[m] += g.data_[m] * weight;
658  }
659 }
684 template <typename T>
685 void function<T>::set_matrixelement(int i1,int i2,function &g,int j1,int j2){
686  assert(0<=i1 && i1<size1_ && 0<=i2 && i2<size2_);
687  assert(0<=j1 && j1<g.size1() && 0<=j2 && j2<g.size2());
688 
689  for(int tstp=-1;tstp<=nt_;tstp++){
690  cdmatrix M;
691  g.get_value(tstp,M);
692  set_matrixelement(tstp,i1,i2,M,j1,j2);
693  }
694 
695 }
713 template <typename T>
714 void function<T>::print_to_file(const char *file, int precision) const {
715  int i, l, sg = element_size_;
716  std::ofstream out;
717  out.open(file, std::ios::out);
718  out.precision(precision);
719  out << "# " << nt_ << " " << size1_ << std::endl;
720  for (i = -1; i <= nt_; i++) {
721  for (l = 0; l < sg; l++)
722  out << ptr(i)[l].real() << " " << ptr(i)[l].imag() << " ";
723  out << std::endl;
724  }
725  out.close();
726 }
744 template <typename T>
745 void function<T>::read_from_file(const char *file) {
746  int i, n, l, size1, sg;
747  double real, imag;
748  std::string s;
749  std::ifstream out;
750  out.open(file, std::ios::in);
751  if (!(out >> s >> n >> size1)) {
752  std::cerr << "read function from file " << file << " error in file"
753  << std::endl;
754  abort();
755  }
756  if (n > nt_ || size1 != size1_)
757  resize(n, size1);
758  sg = element_size_;
759  for (i = -1; i <= n; i++) {
760  for (l = 0; l < sg; l++) {
761  if (!(out >> real >> imag)) {
762  std::cerr << "read function from file " << file
763  << " error at t= " << i << std::endl;
764  abort();
765  }
766  ptr(i)[l] = std::complex<T>(real, imag);
767  }
768  }
769  out.close();
770 }
789 template <typename T>
790 void function<T>::read_from_file(int nt1, const char *file) {
791  int i, n, l, size1, sg;
792  double real, imag;
793  std::string s;
794  std::ifstream out;
795  out.open(file, std::ios::in);
796  if (!(out >> s >> n >> size1)) {
797  std::cerr << "read function from file " << file << " error in file"
798  << std::endl;
799  abort();
800  }
801 
802  assert(nt1 >= -1 && nt1 <= nt_);
803  assert(nt1 <= n);
804  assert(size1 == size1_);
805 
806  sg = size1 * size1;
807  for (i = -1; i <= n; i++) {
808  for (l = 0; l < sg; l++) {
809  if (!(out >> real >> imag)) {
810  std::cerr << "read function from file " << file
811  << " error at t= " << i << std::endl;
812  abort();
813  }
814  if (i <= nt1)
815  ptr(i)[l] = std::complex<T>(real, imag);
816  }
817  }
818  out.close();
819 }
820 
821 #if CNTR_USE_HDF5 == 1
822 // identical from herm_matrix<T>
839 template <typename T>
840 void function<T>::write_to_hdf5(hid_t group_id) const {
841  store_int_attribute_to_hid(group_id, std::string("nt"), nt_);
842  store_int_attribute_to_hid(group_id, std::string("size1"), size1_);
843  store_int_attribute_to_hid(group_id, std::string("size2"), size2_);
844  store_int_attribute_to_hid(group_id, std::string("element_size"),
845  element_size_);
846  hsize_t len_shape = 3, shape[3];
847  shape[1] = size1_;
848  shape[2] = size2_;
849  if (nt_ > -2) {
850  shape[0] = nt_ + 2;
851  store_cplx_array_to_hid(group_id, std::string("data"), this->data_,
852  shape, len_shape);
853  }
854 }
855 // identical from herm_matrix<T>
875 template <typename T>
876 void function<T>::write_to_hdf5(hid_t group_id, const char *groupname) const {
877  hid_t sub_group_id = create_group(group_id, groupname);
878  this->write_to_hdf5(sub_group_id);
879  close_group(sub_group_id);
880 }
881 // identical from herm_matrix<T>
900 template <typename T>
901 void function<T>::write_to_hdf5(const char *filename,
902  const char *groupname) const {
903  hid_t file_id = open_hdf5_file(filename);
904  this->write_to_hdf5(file_id, groupname);
905  close_hdf5_file(file_id);
906 }
924 template <typename T>
925 void function<T>::read_from_hdf5(hid_t group_id) {
926  // -- Read dimensions
927  int nt = read_primitive_type<int>(group_id, "nt");
928  int size1 = read_primitive_type<int>(group_id, "size1");
929  this->resize(nt, size1);
930  if (nt > -2) {
931  hsize_t data_size = (nt + 2) * element_size_;
932  read_primitive_type_array(group_id, "data", data_size, this->data_);
933  }
934 }
954 template <typename T>
955 void function<T>::read_from_hdf5(hid_t group_id, const char *groupname) {
956  hid_t sub_group_id = open_group(group_id, groupname);
957  this->read_from_hdf5(sub_group_id);
958  close_group(sub_group_id);
959 }
979 template <typename T>
980 void function<T>::read_from_hdf5(const char *filename,
981  const char *groupname) {
982  hid_t file_id = read_hdf5_file(filename);
983  this->read_from_hdf5(file_id, groupname);
984  close_hdf5_file(file_id);
985 }
1004 template <typename T>
1005 void function<T>::read_from_hdf5(int nt1, hid_t group_id) {
1006  function<T> ftmp;
1007  ftmp.read_from_hdf5(group_id);
1008 
1009  assert(nt1 >= -1 && nt1 <= ftmp.nt());
1010  assert(nt1 >= -1 && nt1 <= nt_);
1011  assert(ftmp.size1() == size1_);
1012  assert(ftmp.element_size() == element_size_);
1013 
1014  memcpy(data_, ftmp.data_, sizeof(cplx) * (nt1 + 2) * element_size_);
1015 }
1036 template <typename T>
1037 void function<T>::read_from_hdf5(int nt1, hid_t group_id,
1038  const char *groupname) {
1039  hid_t sub_group_id = open_group(group_id, groupname);
1040  this->read_from_hdf5(nt1, sub_group_id);
1041  close_group(sub_group_id);
1042 }
1063 template <typename T>
1064 void function<T>::read_from_hdf5(int nt1, const char *filename,
1065  const char *groupname) {
1066  hid_t file_id = read_hdf5_file(filename);
1067  this->read_from_hdf5(nt1, file_id, groupname);
1068  close_hdf5_file(file_id);
1069 }
1070 #endif
1071 
1072 /* #######################################################################################
1073 #
1074 # MPI UTILS
1075 #
1076 ########################################################################################*/
1077 #if CNTR_USE_MPI==1
1078 
1096 template <typename T>
1097 void function<T>::Bcast_timestep(int tstp,int root){
1098  int numtasks,taskid;
1099  cdmatrix ftemp;
1100  ftemp.resize(size1_,size2_);
1101  MPI_Comm_size(MPI_COMM_WORLD, &numtasks);
1102  MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
1103  if(taskid==root) this->get_value(tstp,ftemp);
1104  // if(sizeof(T)==sizeof(double)) MPI::COMM_WORLD.Bcast(ftemp.data(),ftemp.size(),MPI::DOUBLE_COMPLEX,root);
1105  // else MPI::COMM_WORLD.Bcast(ftemp.data(),ftemp.size(),MPI::COMPLEX,root);
1106  if(sizeof(T)==sizeof(double)) MPI_Bcast(ftemp.data(),ftemp.size(),MPI_DOUBLE_COMPLEX,root,MPI_COMM_WORLD);
1107  if(taskid!=root) this->set_value(tstp,ftemp);
1108 }
1109 #endif
1110 
1111 
1112 } // namespace cntr
1113 
1114 #endif // CNTR_FUNCTION_IMPL_H
int size2_
Number of the rows in the Matrix form.
int size1(void) const
void Bcast_timestep(int tstp, int root)
Broadcast a time step of the function using MPI
void resize(int nt, int size1)
Resize the time length and/or matrix dimension of a square-matrix function
function & operator=(const function &f)
Overloaded operator &#39;=&#39;. It copies data from an existing function object
void set_constant(cplx *f0)
Set all data to a constant(scalar or matrix) for the function class
void set_zero(void)
Set all data to zero for the function class
cplx * data_
Pointer to the function in the Matrix form on the real-time axis ( ) ; &#39;data_+ element_size&#39; corresp...
int size1_
Number of the colums in the Matrix form.
Class function for objects with time on real axis.
void get_value(int tstp, EigenMatrix &M) const
Get matrix value of this function object at a specific time point
void left_multiply(cplx *f, T weight=1.0)
Multiply function values with a matrix-valued function from the left
void read_from_hdf5(hid_t group_id)
Read the function data from a hdf5 file
void incr(function< T > &g, T weight=1.0)
Increase the function object by another function.
void right_multiply(cplx *f, T weight=1.0)
Multiply the function with another matrix-valued function from the right
void set_matrixelement(int tstp, int i1, int i2, EigenMatrix &M, int j1, int j2)
Set a matrix element from an Eigen matrix
int nt(void) const
void get_matrixelement(int i1, int i2, function< T > &g)
Get a matrix element of the function object
std::complex< T > cplx
int nt_
Maximum number of the time steps.
void read_from_file(const char *file)
Read the function data from a txt file
void write_to_hdf5(hid_t group_id) const
Write the function data to a hdf5 file
void print_to_file(const char *file, int precision=16) const
Store the function data to a txt file
void set_value(int tstp, EigenMatrix &M)
Set matrix value at a specific time point
void smul(T weight)
Multiply the function with a scalar