1 #ifndef CNTR_HERM_MATRIX_IMPL_H 2 #define CNTR_HERM_MATRIX_IMPL_H 63 assert(size1 >= 0 && nt >= -1 && sig * sig == 1 && ntau >= 0);
69 element_size_ = size1 * size1;
71 mat_ =
new cplx[(ntau_ + 1) * element_size_];
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_];
112 assert(size1>=0 && size2>=0 && nt>=-1 && sig*sig==1 && ntau>=0);
118 element_size_=size1*size2;
120 mat_ =
new cplx [(ntau_+1)*element_size_];
121 memset(mat_, 0,
sizeof(
cplx)*(ntau_+1)*element_size_);
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_);
156 template <
typename T>
163 element_size_ = size1_ * size1_;
165 mat_ =
new cplx[(ntau_ + 1) * element_size_];
166 memcpy(mat_, g.mat_,
sizeof(
cplx) * (ntau_ + 1) * element_size_);
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_];
175 sizeof(
cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
177 sizeof(
cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
179 sizeof(
cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
186 template <
typename T>
191 if (nt_ != g.nt_ || ntau_ != g.ntau_ || size1_ != g.size1_) {
200 element_size_ = size1_ * size1_;
202 mat_ =
new cplx[(ntau_ + 1) * element_size_];
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_];
217 memcpy(mat_, g.mat_,
sizeof(
cplx) * (ntau_ + 1) * element_size_);
219 memcpy(les_, g.les_,
sizeof(
cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 *
221 memcpy(ret_, g.ret_,
sizeof(
cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 *
224 sizeof(
cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
229 #if __cplusplus >= 201103L 230 template <
typename T>
240 element_size_(g.element_size_),
252 template <
typename T>
265 element_size_ = g.element_size_;
312 template <
typename T>
314 assert(ntau >= 0 && nt >= -1 && size1 >= 0);
323 element_size_ = size1 * size1;
325 mat_ =
new cplx[(ntau_ + 1) * element_size_];
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_];
360 template <
typename T>
362 int nt1 = (nt_ > nt ? nt : nt_);
363 cplx *ret, *les, *tv;
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_];
373 memcpy(les, les_,
sizeof(
cplx) * ((nt1 + 1) * (nt1 + 2)) / 2 *
375 memcpy(ret, ret_,
sizeof(
cplx) * ((nt1 + 1) * (nt1 + 2)) / 2 *
378 sizeof(
cplx) * (nt1 + 1) * (ntau_ + 1) * element_size_);
416 template <
typename T>
420 if (ntau == ntau_ && size1_ == size1)
423 resize_discard(nt, ntau, size1);
436 template <
typename T>
440 memset(mat_, 0,
sizeof(
cplx) * (ntau_ + 1) * element_size_);
443 sizeof(
cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
445 sizeof(
cplx) * ((nt_ + 1) * (nt_ + 2)) / 2 * element_size_);
447 sizeof(
cplx) * (nt_ + 1) * (ntau_ + 1) * element_size_);
456 template <
typename T>
458 assert(t >= 0 && t1 >= 0 && t <= t1 && t1 <= nt_);
459 return ((t1 * (t1 + 1)) / 2 + t) * element_size_;
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_;
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_;
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_;
480 template <
typename T>
481 inline std::complex<T> *herm_matrix<T>::lesptr(
int t,
int t1) {
482 return les_ + les_offset(t, t1);
485 template <
typename T>
486 inline std::complex<T> *herm_matrix<T>::retptr(
int t,
int t1) {
487 return ret_ + ret_offset(t, t1);
490 template <
typename T>
491 inline std::complex<T> *herm_matrix<T>::tvptr(
int t,
int tau) {
492 return tv_ + tv_offset(t, tau);
495 template <
typename T>
496 inline std::complex<T> *herm_matrix<T>::matptr(
int tau) {
497 return mat_ + mat_offset(tau);
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);
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);
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);
515 template <
typename T>
516 inline const std::complex<T> *herm_matrix<T>::matptr(
int tau)
const {
517 return mat_ + mat_offset(tau);
526 #define herm_matrix_READ_ELEMENT \ 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]; \ 535 #define herm_matrix_READ_ELEMENT_MINUS_CONJ \ 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()); \ 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_);
575 herm_matrix_READ_ELEMENT
578 herm_matrix_READ_ELEMENT_MINUS_CONJ
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_);
612 herm_matrix_READ_ELEMENT
615 herm_matrix_READ_ELEMENT_MINUS_CONJ
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
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;
694 template <
typename T>
695 template <
class Matrix>
696 void herm_matrix<T>::get_mat(
int i, Matrix &M)
const {
698 const cplx *x = matptr(i);
699 herm_matrix_READ_ELEMENT
721 template <
typename T>
722 template <
class Matrix>
723 void herm_matrix<T>::get_matminus(
int i, Matrix &M)
const {
725 const cplx *x = matptr(ntau_ - i);
726 herm_matrix_READ_ELEMENT
if (sig_ == -1) M = -M;
728 #undef herm_matrix_READ_ELEMENT 729 #undef herm_matrix_READ_ELEMENT_MINUS_CONJ 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_);
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_);
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_);
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_);
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);
903 template <
typename T>
904 inline void herm_matrix<T>::get_mat(
int i,
cplx &x)
const {
928 template <
typename T>
929 inline void herm_matrix<T>::get_matminus(
int i,
cplx &x)
const {
931 x = *matptr(ntau_ - i);
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_);
983 template <
typename T>
984 std::complex<T> herm_matrix<T>::density_matrix(
int tstp)
const {
985 assert(tstp >= -1 && tstp <= nt_);
992 return std::complex<T>(0.0, sig_) * x1;
1016 template <
typename T>
1017 template <
class Matrix>
1019 assert(M.rows() == size1_ && M.cols() == size2_);
1020 assert(tstp >= -1 && tstp <= nt_);
1026 M *= std::complex<T>(0.0, 1.0 * sig_);
1037 #define herm_matrix_SET_ELEMENT_MATRIX \ 1040 for (r = 0; r < size1_; r++) \ 1041 for (s = 0; s < size2_; s++) \ 1042 x[r * size2_ + s] = M(r, s); \ 1063 template <
typename T>
1064 template <
class Matrix>
1066 assert(i <= nt_ && j <= nt_);
1067 cplx *x = retptr(i, j);
1068 herm_matrix_SET_ELEMENT_MATRIX
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
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
1138 template <
typename T>
1139 template <
class Matrix>
1140 void herm_matrix<T>::set_mat(
int i, Matrix &M) {
1142 cplx *x = matptr(i);
1143 herm_matrix_SET_ELEMENT_MATRIX
1161 template<
typename T>
1162 void herm_matrix<T>::set_mat_herm(
int i){
1163 cdmatrix tmp(size1_,size2_);
1165 tmp=(tmp+tmp.adjoint())/2.0;
1166 this->set_mat(i,tmp);
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++){
1184 tmp=(tmp+tmp.adjoint())/2.0;
1185 this->set_mat(i,tmp);
1206 template <
typename T>
1207 inline void herm_matrix<T>::set_les(
int i,
int j,
cplx x) {
1208 assert(i <= nt_ && j <= nt_);
1230 template <
typename T>
1231 inline void herm_matrix<T>::set_ret(
int i,
int j,
cplx x) {
1232 assert(i <= nt_ && j <= nt_);
1253 template <
typename T>
1254 inline void herm_matrix<T>::set_tv(
int i,
int j,
cplx x) {
1255 assert(i <= nt_ && j <= ntau_);
1274 template <
typename T>
1275 inline void herm_matrix<T>::set_mat(
int i,
cplx x) {
1276 assert(i <= ntau_ );
1307 template <
typename T>
1309 int i, j, l, sg = element_size_;
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();
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();
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();
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();
1375 template <
typename T>
1377 int i, n, m, j, l, size1, sg, sig;
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" 1387 if (n > nt_ || m != ntau_ || size1 != size1_)
1388 resize(n, m, size1);
1391 for (j = 0; j <= ntau_; j++) {
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;
1399 matptr(j)[l] = std::complex<T>(real, imag);
1403 for (i = 0; i <= n; i++) {
1404 for (j = 0; j <= i; j++) {
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 <<
")" 1413 retptr(i, j)[l] = std::complex<T>(real, imag);
1417 for (i = 0; i <= n; i++) {
1418 for (j = 0; j <= ntau_; j++) {
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 <<
")" 1427 tvptr(i, j)[l] = std::complex<T>(real, imag);
1431 for (j = 0; j <= n; j++) {
1432 for (i = 0; i <= j; i++) {
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 <<
")" 1441 lesptr(i, j)[l] = std::complex<T>(real, imag);
1468 template <
typename T>
1470 int i, n, m, j, l, size1, sg, sig;
1474 out.open(file, std::ios::in);
1476 if (!(out >> s >> n >> m >> size1 >> sig)) {
1477 std::cerr <<
"read G from file " << file <<
" error in file" 1482 assert(nt1 <= nt_ && "nt1 > nt_
"); 1483 assert(nt1 <= n && "nt1 > n
"); 1484 assert(m == ntau_ && "m /= ntau_
"); 1485 assert(size1 == size1_ && "size1 /= size1_
"); 1489 for (j = 0; j <= ntau_; j++) { 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; 1497 matptr(j)[l] = std::complex<T>(real, imag); 1502 for (i = 0; i <= n; i++) { 1503 for (j = 0; j <= i; j++) { 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 << ")
" 1513 retptr(i, j)[l] = std::complex<T>(real, imag); 1517 for (i = 0; i <= n; i++) { 1518 for (j = 0; j <= ntau_; j++) { 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 << ")
" 1528 tvptr(i, j)[l] = std::complex<T>(real, imag); 1532 for (j = 0; j <= n; j++) { 1533 for (i = 0; i <= j; i++) { 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 << ")
" 1543 lesptr(i, j)[l] = std::complex<T>(real, imag); 1551 #if CNTR_USE_HDF5 == 1 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
"), 1579 hsize_t len_shape = 3, shape[3]; 1583 shape[0] = ntau_ + 1; 1584 store_cplx_array_to_hid(group_id, std::string("mat
"), matptr(0), 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), 1592 store_cplx_array_to_hid(group_id, std::string("les
"), lesptr(0, 0), 1594 shape[0] = (nt_ + 1) * (ntau_ + 1); 1595 store_cplx_array_to_hid(group_id, std::string("tv
"), tvptr(0, 0), 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); 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); 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
"); 1675 this->resize(nt, ntau, size1); 1678 hsize_t mat_size = (ntau + 1) * element_size_; 1679 read_primitive_type_array(group_id, "mat
", mat_size, matptr(0)); 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)); 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); 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); 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); 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_
"); 1772 for (int n = -1; n <= nt1; n++) 1773 this->set_timestep(n, gtmp); 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); 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); 1851 template <typename T> 1852 void herm_matrix<T>::write_to_hdf5_slices(hid_t group_id, int dt) { 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); 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); 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); 1938 template <typename T> 1939 void herm_matrix<T>::write_to_hdf5_tavtrel(hid_t group_id, int dt) { 1942 typedef std::complex<double> complex; 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
"), 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]; 1958 for (int tav = 0; tav <= nt_; tav += dt) { 1959 int len = (tav <= nt_ - tav ? tav : nt_ - tav); 1960 std::sprintf(name, "%d
", tav); 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_); 1974 store_cplx_array_to_hid(group_id_les, std::string(name), gles, shape, 1976 store_cplx_array_to_hid(group_id_gtr, std::string(name), ggtr, shape, 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); 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); 2041 /* ####################################################################################### 2043 # SIMPLE OPERATIONS ON TIMESTEPS 2044 # NOTE: tstp IS A PHYSICAL TIME, tstp=-1 is the matsubara branch 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_
"); 2068 memset(matptr(0), 0, sizeof(cplx) * (ntau_ + 1) * element_size_); 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_); 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_
"); 2103 memcpy(mat_, g1.mat_, sizeof(cplx) * (ntau_ + 1) * element_size_); 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_); 2136 template <typename T> 2137 void herm_matrix<T>::set_timestep(int tstp, 2138 herm_matrix_timestep<T> ×tep) { 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_
"); 2145 memcpy(mat_, x, sizeof(cplx) * (ntau_ + 1) * element_size_); 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_); 2155 template <typename T> 2156 void herm_matrix<T>::get_timestep(int tstp, 2157 herm_matrix_timestep<T> ×tep) const { 2158 int len = (2 * (tstp + 1) + ntau_ + 1) * element_size_; 2160 assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_
"); 2161 if (timestep.total_size_ < len) 2162 timestep.resize(tstp, ntau_, size1_); 2164 timestep.tstp_ = tstp; 2165 timestep.ntau_ = ntau_; 2166 timestep.size1_ = size1_; 2168 memcpy(x, mat_, sizeof(cplx) * (ntau_ + 1) * element_size_); 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_); 2178 template <typename T> 2179 void herm_matrix<T>::get_timestep(int tstp, herm_matrix<T> &g1) const { 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_
"); 2185 memcpy(g1.mat_, mat_, sizeof(cplx) * (ntau_ + 1) * element_size_); 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_); 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); 2228 for (i = 0; i <= ntau_; i++) 2229 matptr(i)[sij] = *(g.matptr(i)); 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)); 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
"); 2271 for (i = 0; i <= ntau_; i++) 2272 matptr(i)[sij] = *(g.matptr(i)); 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)); 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, 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()
"); 2322 for (i = 0; i <= ntau_; i++) 2323 matptr(i)[sij] = g.matptr(i)[tij]; 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]; 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()
"); 2374 for (i = 0; i <= ntau_; i++) 2375 matptr(i)[sij] = g.matptr(i)[tij]; 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]; 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); 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); 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 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]); 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]); 2534 #define HERM_MATRIX_INCR_TSTP \ 2535 if (alpha == cplx(1.0, 0.0)) { \ 2536 for (i = 0; i < len; i++) \ 2539 for (i = 0; i < len; i++) \ 2540 x0[i] += alpha * x[i]; \ 2566 template <typename T> 2567 void herm_matrix<T>::incr_timestep(int tstp, 2568 herm_matrix_timestep<T> ×tep, 2569 std::complex<T> alpha) { 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_
"); 2576 len = (ntau_ + 1) * element_size_; 2579 HERM_MATRIX_INCR_TSTP 2581 len = (tstp + 1) * element_size_; 2582 x0 = retptr(tstp, 0); 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 2615 template <typename T> 2616 void herm_matrix<T>::incr_timestep(int tstp, 2617 herm_matrix_timestep<T> ×tep) { 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_
"); 2625 len = (ntau_ + 1) * element_size_; 2628 HERM_MATRIX_INCR_TSTP 2630 len = (tstp + 1) * element_size_; 2631 x0 = retptr(tstp, 0); 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 2645 #undef HERM_MATRIX_INCR_TSTP 2668 template <typename T> 2669 void herm_matrix<T>::incr_timestep(int tstp, herm_matrix<T> &g, 2670 std::complex<T> alpha) { 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_
"); 2679 for (m = 0; m <= ntau_; m++) { 2680 element_incr<T, LARGESIZE>(size1_, x0 + m * element_size_, alpha, 2681 x + m * element_size_); 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_); 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_); 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_); 2724 template <typename T> 2725 void herm_matrix<T>::incr_timestep(int tstp, herm_matrix<T> &g) { 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_
"); 2735 for (m = 0; m <= ntau_; m++) { 2736 element_incr<T, LARGESIZE>(size1_, x0 + m * element_size_, alpha, 2737 x + m * element_size_); 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_); 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_); 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_); 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); 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); 2815 /** \brief <b> Left-multiplies the `herm_matrix` with contour function. </b> 2817 * <!-- ====== DOCUMENTATION ====== --> 2820 * <!-- ========= --> 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. 2831 * > [int] The time step at which \f$F(t)\f$ and \f$C(t,t^\prime)\f$ are multiplied. 2833 * > [complex<T>] Pointer to \f$F(-\mathrm{i}\beta)\f$ on the Matsubara axis (this is just a constant matrix). 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$. 2837 * > [T] The weight as above. 2839 template <typename T> 2840 void herm_matrix<T>::left_multiply(int tstp, std::complex<T> *f0, 2841 std::complex<T> *ft, T weight) { 2843 cplx *xtemp, *ftemp, *x0; 2844 xtemp = new cplx[element_size_]; 2845 assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_
"); 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); 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); 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); 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); 2902 template <typename T> 2903 void herm_matrix<T>::left_multiply(int tstp, function<T> &ft, T weight) { 2906 left_multiply_(tstp, ft, weight, std::integral_constant<int, 1>{}); 2909 left_multiply_(tstp, ft, weight, std::integral_constant<int, 2>{}); 2912 left_multiply_(tstp, ft, weight, std::integral_constant<int, 3>{}); 2915 left_multiply_(tstp, ft, weight, std::integral_constant<int, 4>{}); 2918 left_multiply_(tstp, ft, weight, std::integral_constant<int, 5>{}); 2921 left_multiply_(tstp, ft, weight, std::integral_constant<int, 8>{}); 2924 left_multiply_(tstp, ft, weight, std::integral_constant<int, LARGESIZE>{}); 2929 /** \brief <b> Left-multiplies the `herm_matrix` with contour function. Internal version, not intended 2930 * to be called directly. </b> 2932 * <!-- ====== DOCUMENTATION ====== --> 2935 * <!-- ========= --> 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. 2941 * \note parameter SIZE determines the matrix size specified by the top-level interface. 2947 * > [int] The time step at which \f$F(t)\f$ and \f$C(t,t^\prime)\f$ are multiplied. 2949 * > [function] The contour function. 2951 * > [T] The weight as above. 2953 template <typename T> 2955 void herm_matrix<T>::left_multiply_(int tstp, function<T> &ft, T weight, 2956 std::integral_constant<int, SIZE>) { 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_
"); 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); 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); 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); 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); 3016 template <typename T> 3017 void herm_matrix<T>::left_multiply_hermconj(int tstp, function<T> &ft, 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_
"); 3028 element_conj<T, LARGESIZE>(size1_, fcc, f0); 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); 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); 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); 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); 3067 /** \brief <b> Right-multiplies the `herm_matrix` with contour function. </b> 3069 * <!-- ====== DOCUMENTATION ====== --> 3072 * <!-- ========= --> 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. 3083 * > [int] The time step at which \f$F(t)\f$ and \f$C(t,t^\prime)\f$ are multiplied. 3085 * > [complex<T>] Pointer to \f$F(-\mathrm{i}\beta)\f$ on the Matsubara axis (this is just a constant matrix). 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$. 3089 * > [T] The weight as above. 3091 template <typename T> 3092 void herm_matrix<T>::right_multiply(int tstp, std::complex<T> *f0, 3093 std::complex<T> *ft, T weight) { 3095 cplx *xtemp, *ftemp, *x0; 3096 xtemp = new cplx[element_size_]; 3097 assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_
"); 3100 for (m = 0; m <= ntau_; m++) { 3101 element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_, 3103 element_smul<T, LARGESIZE>(size1_, xtemp, weight); 3104 element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp); 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); 3114 x0 = tvptr(tstp, 0); 3115 for (m = 0; m <= ntau_; m++) { 3116 element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_, 3118 element_smul<T, LARGESIZE>(size1_, xtemp, weight); 3119 element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp); 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_, 3126 element_smul<T, LARGESIZE>(size1_, xtemp, weight); 3127 element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp); 3153 template <typename T> 3154 void herm_matrix<T>::right_multiply(int tstp, function<T> &ft, T weight) { 3157 right_multiply_(tstp, ft, weight, std::integral_constant<int, 1>{}); 3160 right_multiply_(tstp, ft, weight, std::integral_constant<int, 2>{}); 3163 right_multiply_(tstp, ft, weight, std::integral_constant<int, 3>{}); 3166 right_multiply_(tstp, ft, weight, std::integral_constant<int, 4>{}); 3169 right_multiply_(tstp, ft, weight, std::integral_constant<int, 5>{}); 3172 right_multiply_(tstp, ft, weight, std::integral_constant<int, 8>{}); 3175 right_multiply_(tstp, ft, weight, std::integral_constant<int, LARGESIZE>{}); 3180 /** \brief <b> Right-multiplies the `herm_matrix` with contour function. 3181 * Internal version, not intended to be called directly. </b> 3183 * <!-- ====== DOCUMENTATION ====== --> 3186 * <!-- ========= --> 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. 3192 * \note Parameter SIZE determines the matrix size specified by the top-level interface. 3198 * > The time step at which \f$F(t)\f$ and \f$C(t,t^\prime)\f$ are multiplied. 3200 * > The contour function. 3202 * > The weight as above. 3204 template <typename T> 3206 void herm_matrix<T>::right_multiply_(int tstp, function<T> &ft, T weight, 3207 std::integral_constant<int, SIZE>) { 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
"); 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); 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); 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); 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); 3269 template <typename T> 3270 void herm_matrix<T>::right_multiply_hermconj(int tstp, function<T> &ft, 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
"); 3281 element_conj<T, LARGESIZE>(size1_, fcc, f0); 3284 for (m = 0; m <= ntau_; m++) { 3285 element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_, 3287 element_smul<T, LARGESIZE>(size1_, xtemp, weight); 3288 element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp); 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_, 3296 element_smul<T, LARGESIZE>(size1_, xtemp, weight); 3297 element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp); 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_, 3304 element_smul<T, LARGESIZE>(size1_, xtemp, weight); 3305 element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp); 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_, 3313 element_smul<T, LARGESIZE>(size1_, xtemp, weight); 3314 element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp); 3338 template <typename T> 3339 void herm_matrix<T>::smul(int tstp, T weight) { 3342 assert(tstp >= -1 && tstp <= nt_); 3345 for (m = 0; m <= ntau_; m++) { 3346 element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, 3350 x0 = retptr(tstp, 0); 3351 for (m = 0; m <= tstp; m++) { 3352 element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, 3355 x0 = tvptr(tstp, 0); 3356 for (m = 0; m <= ntau_; m++) { 3357 element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, 3360 x0 = lesptr(0, tstp); 3361 for (m = 0; m <= tstp; m++) { 3362 element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, 3386 template <typename T> 3387 void herm_matrix<T>::smul(int tstp, std::complex<T> weight) { 3390 assert(tstp >= -1 && tstp <= nt_ && "tstp >= -1 && tstp <= nt_
"); 3393 for (m = 0; m <= ntau_; m++) { 3394 element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, 3398 x0 = retptr(tstp, 0); 3399 for (m = 0; m <= tstp; m++) { 3400 element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, 3403 x0 = tvptr(tstp, 0); 3404 for (m = 0; m <= ntau_; m++) { 3405 element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, 3408 x0 = lesptr(0, tstp); 3409 for (m = 0; m <= tstp; m++) { 3410 element_smul<T, LARGESIZE>(size1_, x0 + m * element_size_, 3416 /* ####################################################################################### 3420 ########################################################################################*/ 3421 #if CNTR_USE_MPI == 1 3441 template <typename T> 3442 void herm_matrix<T>::Reduce_timestep(int tstp, int root) { 3443 assert(tstp <= nt_); 3445 herm_matrix_timestep<T> Gtemp; 3446 Gtemp.resize(tstp, ntau_, size1_); 3447 this->get_timestep(tstp, Gtemp); 3449 Gtemp.Reduce_timestep(tstp, root); 3451 this->set_timestep(tstp, Gtemp); 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_); 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); 3485 MPI_Bcast(Gtemp.data_, Gtemp.total_size_, MPI_COMPLEX, 3486 root, MPI_COMM_WORLD); 3488 this->set_timestep(tstp, Gtemp); 3510 template <typename T> 3511 void herm_matrix<T>::Send_timestep(int tstp, int dest, int tag) { 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); 3522 MPI_Send(Gtemp.data_, Gtemp.total_size_, MPI_COMPLEX, 3523 dest, tag, MPI_COMM_WORLD); 3546 template <typename T> 3547 void herm_matrix<T>::Recv_timestep(int tstp, int root, int tag) { 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); 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); 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)
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
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.