1 #ifndef CNTR_FUNCTION_IMPL_H 2 #define CNTR_FUNCTION_IMPL_H 46 int len = (nt + 2) * size1 * size1;
47 assert(size1 >= 0 && nt >= -1);
51 data_ =
new cplx[len];
55 element_size_ = size1 * size1_;
57 total_size_ = (nt_ + 2) * size1_ * size2_;
82 int len=(nt+2)*size1*size2;
83 assert(size1>=0 && nt>=-1);
86 data_ =
new cplx [len];
87 memset(data_, 0,
sizeof(
cplx)*len);
91 element_size_=size1*size2_;
112 template <
typename T>
118 element_size_ = size1_ * size2_;
119 len = (nt_ + 2) * size1_ * size2_;
122 data_ =
new cplx[len];
123 memcpy(data_, g.data_,
sizeof(
cplx) * len);
128 #if __cplusplus >= 201103L 147 template <
typename T>
153 element_size_(g.element_size_) {
178 template <
typename T>
187 element_size_ = g.element_size_;
215 template <
typename T>
225 element_size_ = size1_ * size2_;
226 len = (nt_ + 2) * size1_ * size2_;
229 data_ =
new cplx[len];
230 memcpy(data_, g.data_,
sizeof(
cplx) * len);
254 template <
typename T>
256 int len = (nt + 2) * size1 * size1;
257 assert(nt >= -1 && size1 >= 0);
263 data_ =
new cplx[len];
267 element_size_ = size1_ * size1_;
285 template <
typename T>
287 int len = element_size_ * (nt_ + 2);
289 memset(data_, 0,
sizeof(
cplx) * len);
307 template <
typename T>
310 for (t = -1; t <= nt_; t++)
311 element_set<T, LARGESIZE>(size1_, ptr(t), f0);
328 template <
typename T>
template<
class EigenMatrix>
330 assert(M.rows() == size1_);
332 for(t=-1;t<=nt_;t++){
355 template <
typename T>
template<
class EigenMatrix>
381 template <
typename T>
382 template <
class EigenMatrix>
386 assert(M.rows() == size1_);
387 assert(M.cols() == size2_);
390 for (i1 = 0; i1 < size1_; i1++)
391 for (i2 = 0; i2 < size2_; i2++)
392 ft[i1 * size2_ + i2] = M(i1, i2);
412 template <
typename T>
413 template <
class EigenMatrix>
417 M.resize(size1_, size2_);
419 for (i1 = 0; i1 < size1_; i1++)
420 for (i2 = 0; i2 < size2_; i2++)
421 M(i1, i2) = ft[i1 * size2_ + i2];
439 template <
typename T>
442 for(
int m = 0; m < total_size_; m++)
463 template <
typename T>
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_;
471 element_mult<T, LARGESIZE>(size1_, xtemp,
473 element_smul<T, LARGESIZE>(size1_, xtemp, weight);
474 element_set<T, LARGESIZE>(size1_, mytemp, xtemp);
497 template <
typename T>
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_;
505 element_mult<T, LARGESIZE>(size1_, xtemp,
507 element_smul<T, LARGESIZE>(size1_, xtemp, weight);
508 element_set<T, LARGESIZE>(size1_, mytemp, xtemp);
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);
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);
589 template <
typename T>
template<
class EigenMatrix>
591 assert(0<=i1 && i1<size1_ && 0<=i2 && i2<size2_);
592 assert(0<=j1 && j1<M.rows() && 0<=j2 && j2<M.cols());
595 ft[i1*size2_+i2]=M(j1,j2);
619 template <
typename T>
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++)
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);
651 template <
typename T>
654 assert(
this -> nt_ == g.
nt_);
655 for(
int m = 0; m < total_size_; m++)
657 this -> data_[m] += g.
data_[m] * weight;
684 template <
typename T>
686 assert(0<=i1 && i1<size1_ && 0<=i2 && i2<size2_);
687 assert(0<=j1 && j1<g.size1() && 0<=j2 && j2<g.size2());
689 for(
int tstp=-1;tstp<=nt_;tstp++){
692 set_matrixelement(tstp,i1,i2,M,j1,j2);
713 template <
typename T>
715 int i, l, sg = element_size_;
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() <<
" ";
744 template <
typename T>
746 int i, n, l, size1, sg;
750 out.open(file, std::ios::in);
751 if (!(out >> s >> n >> size1)) {
752 std::cerr <<
"read function from file " << file <<
" error in file" 756 if (n > nt_ || size1 != size1_)
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;
766 ptr(i)[l] = std::complex<T>(real, imag);
789 template <
typename T>
791 int i, n, l, size1, sg;
795 out.open(file, std::ios::in);
796 if (!(out >> s >> n >> size1)) {
797 std::cerr <<
"read function from file " << file <<
" error in file" 802 assert(nt1 >= -1 && nt1 <= nt_);
804 assert(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;
815 ptr(i)[l] = std::complex<T>(real, imag);
821 #if CNTR_USE_HDF5 == 1 839 template <
typename T>
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"),
846 hsize_t len_shape = 3, shape[3];
851 store_cplx_array_to_hid(group_id, std::string(
"data"), this->data_,
875 template <
typename T>
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);
900 template <
typename T>
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);
924 template <
typename T>
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);
931 hsize_t data_size = (nt + 2) * element_size_;
932 read_primitive_type_array(group_id,
"data", data_size, this->data_);
954 template <
typename T>
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);
979 template <
typename T>
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);
1004 template <
typename T>
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_);
1014 memcpy(data_, ftmp.
data_,
sizeof(
cplx) * (nt1 + 2) * element_size_);
1036 template <
typename T>
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);
1063 template <
typename T>
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);
1096 template <
typename T>
1098 int numtasks,taskid;
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);
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);
1114 #endif // CNTR_FUNCTION_IMPL_H int size2_
Number of the rows in the Matrix form.
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 '='. 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 ( ) ; 'data_+ element_size' 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
void get_matrixelement(int i1, int i2, function< T > &g)
Get a matrix element of the function object
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