cntr_herm_matrix_timestep_impl.hpp

Go to the documentation of this file.

#ifndef CNTR_HERM_MATRIX_TIMESTEP_IMPL_H
#define CNTR_HERM_MATRIX_TIMESTEP_IMPL_H

#include "cntr_herm_matrix_timestep_decl.hpp"
#include "cntr_elements.hpp"
#include "cntr_function_decl.hpp"
#include "cntr_herm_matrix_decl.hpp"
#include "cntr_herm_matrix_timestep_view_impl.hpp"

namespace cntr {

/* #######################################################################################
#
#   CONSTRUCTION/DESTRUCTION
#
########################################################################################*/
template <typename T>
   herm_matrix_timestep<T>::herm_matrix_timestep() {
      data_ = 0;
      ntau_ = 0;
      tstp_ = 0;
      size1_ = 0;
      size2_ = 0;
      element_size_ = 0;
      total_size_ = 0;
      sig_ = -1;
   }
template <typename T>
   herm_matrix_timestep<T>::~herm_matrix_timestep() {
      delete[] data_;
   }

template <typename T> herm_matrix_timestep<T>::herm_matrix_timestep(int tstp,int ntau,int size1){
   int len=((tstp+1)*2+(ntau+1))*size1*size1;
   assert(size1>=0 && tstp>=-1 && ntau>=0);
   if(len==0) data_=0;
   else{
      data_ = new cplx [len];
      memset(data_, 0, sizeof(cplx)*len);
   }
   size1_=size1;
   size2_=size1;
   element_size_=size1*size1_;
   tstp_=tstp;
   ntau_=ntau;
   total_size_=len;
   sig_=-1;
}

template <typename T> herm_matrix_timestep<T>::herm_matrix_timestep(int tstp,int ntau,int size1,int size2,int sig){
int len=((tstp+1)*2+(ntau+1))*size1*size2;
assert(size1>=0 && tstp>=-1 && ntau>=0);
if(len==0) data_=0;
else{
   data_ = new cplx [len];
   memset(data_, 0, sizeof(cplx)*len);
}
size1_=size1;
size2_=size2;
element_size_=size1*size2;
tstp_=tstp;
ntau_=ntau;
total_size_=len;
sig_=sig;
}

template <typename T> herm_matrix_timestep<T>::herm_matrix_timestep(int tstp,int ntau,int size1,int sig){
int len=((tstp+1)*2+(ntau+1))*size1*size1;
assert(size1>=0 && tstp>=-1 && ntau>=0 && sig*sig==1);
if(len==0) data_=0;
else{
   data_ = new cplx [len];
   memset(data_, 0, sizeof(cplx)*len);
}
size1_=size1;
size2_=size1;
element_size_=size1*size1_;
tstp_=tstp;
ntau_=ntau;
total_size_=len;
sig_=sig;
}

template <typename T>
herm_matrix_timestep<T>::herm_matrix_timestep(const herm_matrix_timestep &g) {
   tstp_ = g.tstp_;
   ntau_ = g.ntau_;
   size1_ = g.size1_;
   size2_ = g.size1_;
   element_size_ = size1_ * size1_;
   total_size_ = g.total_size_;
   if (total_size_ > 0) {
      data_ = new cplx[total_size_];
      memcpy(data_, g.data_, sizeof(cplx) * total_size_);
   } else {
      data_ = 0;
   }
   sig_ = g.sig_;
}
template <typename T>
herm_matrix_timestep<T> &herm_matrix_timestep<T>::
operator=(const herm_matrix_timestep &g) {
   if (this == &g)
      return *this;
   if (total_size_ != g.total_size_) {
      delete[] data_;
      total_size_ = g.total_size_;
      if (total_size_ > 0)
         data_ = new cplx[total_size_];
      else
         data_ = 0;
   }
   tstp_ = g.tstp_;
   ntau_ = g.ntau_;
   size1_ = g.size1_;
   size2_ = g.size1_;
   sig_ = g.sig_;
   element_size_ = size1_ * size1_;
   if (total_size_ > 0)
      memcpy(data_, g.data_, sizeof(cplx) * total_size_);
   return *this;
}
#if __cplusplus >= 201103L
template <typename T>
herm_matrix_timestep<T>::herm_matrix_timestep(
   herm_matrix_timestep &&g) noexcept : data_(g.data_),
tstp_(g.tstp_),
ntau_(g.ntau_),
size1_(g.size1_),
size2_(g.size2_),
element_size_(g.element_size_),
total_size_(g.total_size_),
sig_(g.sig_) {
   g.data_ = nullptr;
   g.tstp_ = 0;
   g.ntau_ = 0;
   g.size1_ = 0;
   g.size2_ = 0;
   g.element_size_ = 0;
   g.total_size_ = 0;
}
template <typename T>
herm_matrix_timestep<T> &herm_matrix_timestep<T>::
operator=(herm_matrix_timestep &&g) noexcept {
   if (&g == this)
      return *this;

   data_ = g.data_;
   tstp_ = g.tstp_;
   ntau_ = g.ntau_;
   size1_ = g.size1_;
   size2_ = g.size2_;
   element_size_ = g.element_size_;
   total_size_ = g.total_size_;
   sig_ = g.sig_;

   g.data_ = nullptr;
   g.tstp_ = 0;
   g.ntau_ = 0;
   g.size1_ = 0;
   g.size2_ = 0;
   g.element_size_ = 0;
   g.total_size_ = 0;

   return *this;
}
#endif

template <typename T>
void herm_matrix_timestep<T>::resize(int tstp, int ntau, int size1) {
   int len = ((tstp + 1) * 2 + (ntau + 1)) * size1 * size1;
   assert(ntau >= 0 && tstp >= -1 && size1 >= 0);
   delete[] data_;
   if (len == 0)
      data_ = 0;
   else {
      data_ = new cplx[len];
   }
   size1_ = size1;
   size2_ = size1;
   element_size_ = size1_ * size1_;
   tstp_ = tstp;
   ntau_ = ntau;
   total_size_ = len;
}

template <typename T>
void herm_matrix_timestep<T>::clear(void) {
   if (total_size_ > 0)
      memset(data_, 0, sizeof(cplx) * total_size_);
}


template <typename T>
void herm_matrix_timestep<T>::set_timestep_zero(int tstp) {
    assert(tstp == tstp_);
    assert(tstp >= -1);
    if (tstp == -1) {
        memset(matptr(0), 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
    } else {
        memset(retptr(0), 0, sizeof(cplx) * (tstp + 1) * element_size_);
        memset(tvptr(0), 0, sizeof(cplx) * (ntau_ + 1) * element_size_);
        memset(lesptr(0), 0, sizeof(cplx) * (tstp + 1) * element_size_);
    }
}


template <typename T>
void herm_matrix_timestep<T>::set_timestep(int tstp, herm_matrix<T> &g1) {
    assert(tstp == tstp_);
    assert(tstp >= -1 && tstp <= g1.nt() && "tstp >= -1 && tstp <= g1.nt()");
    assert(g1.size1() == size1_ && "g1.size1() == size1_");
    assert(g1.ntau() == ntau_ && "g1.ntau() == ntau_");
    if (tstp == -1) {
        memcpy(matptr(0), g1.matptr(0), sizeof(cplx) * (ntau_ + 1) * element_size_);
    } else {
        memcpy(retptr(0), g1.retptr(tstp, 0),
               sizeof(cplx) * (tstp + 1) * element_size_);
        memcpy(tvptr(0), g1.tvptr(tstp, 0),
               sizeof(cplx) * (ntau_ + 1) * element_size_);
        memcpy(lesptr(0), g1.lesptr(0, tstp),
               sizeof(cplx) * (tstp + 1) * element_size_);
    }
}

template <typename T>
void herm_matrix_timestep<T>::set_timestep(int tstp, herm_matrix_timestep<T> &g1) {
    assert(tstp == tstp_);
    assert(tstp == g1.tstp());
    assert(tstp >= -1 && tstp <= g1.tstp() && "tstp >= -1 && tstp <= g1.tstp()");
    assert(g1.size1() == size1_ && "g1.size1() == size1_");
    assert(g1.ntau() == ntau_ && "g1.ntau() == ntau_");
    if (tstp == -1) {
        memcpy(matptr(0), g1.matptr(0), sizeof(cplx) * (ntau_ + 1) * element_size_);
    } else {
        memcpy(retptr(0), g1.retptr(0),
               sizeof(cplx) * (tstp + 1) * element_size_);
        memcpy(tvptr(0), g1.tvptr(0),
               sizeof(cplx) * (ntau_ + 1) * element_size_);
        memcpy(lesptr(0), g1.lesptr(0),
               sizeof(cplx) * (tstp + 1) * element_size_);
    }
}



// the following routines are not very "efficent" but sometimes simple to
// implement
#define herm_matrix_READ_ELEMENT                                             \
{                                                                        \
   int r, s, dim = size1_;                                              \
   M.resize(dim, dim);                                                  \
   for (r = 0; r < dim; r++)                                            \
      for (s = 0; s < dim; s++)                                        \
         M(r, s) = x[r * dim + s];                                    \
   }
#define herm_matrix_READ_ELEMENT_MINUS_CONJ                                  \
   {                                                                        \
      cplx w;                                                              \
      int r, s, dim = size1_;                                              \
      M.resize(dim, dim);                                                  \
      for (r = 0; r < dim; r++)                                            \
         for (s = 0; s < dim; s++) {                                      \
            w = x[s * dim + r];                                          \
            M(r, s) = std::complex<T>(-w.real(), w.imag());              \
         }                                                                \
      }


template <typename T> template <class Matrix>
      void herm_matrix_timestep<T>::get_les(int i, int j, Matrix &M){
         assert(j == tstp_ && i <= tstp_);
         cplx *x;
         x = lesptr(i);
         herm_matrix_READ_ELEMENT
      }


template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_les_t_tstp(int i, Matrix &M) {
         cplx *x;
         x = lesptr(i);
         herm_matrix_READ_ELEMENT
      }

template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_les_tstp_t(int i, Matrix &M) {
         cplx *x;
         x = lesptr(i);
         herm_matrix_READ_ELEMENT_MINUS_CONJ
      }



template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_ret(int i, int j, Matrix &M) {
         assert(i == tstp_ && j <= tstp_);
         cplx *x;
         x = retptr(j);
         herm_matrix_READ_ELEMENT
      }



template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_ret_tstp_t(int j, Matrix &M) {
         cplx *x;
         x = retptr(j);
         herm_matrix_READ_ELEMENT
      }

template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_ret_t_tstp(int i, Matrix &M) {
         cplx *x;
         x = retptr(i);
         herm_matrix_READ_ELEMENT_MINUS_CONJ
      }


template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_tv(int i, int j, Matrix &M) {
         assert(i == tstp_);
         cplx *x = tvptr(j);
         herm_matrix_READ_ELEMENT
      }


template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_tv(int j, Matrix &M) {
         cplx *x = tvptr(j);
         herm_matrix_READ_ELEMENT
      }

template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_vt(int i, Matrix &M, int sig) {
         cplx *x = tvptr(ntau_ - i);
         herm_matrix_READ_ELEMENT_MINUS_CONJ if (sig == -1) M = -M;
      }

template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_mat(int i, Matrix &M) {
         cplx *x = matptr(i);
         herm_matrix_READ_ELEMENT
      }

template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_matminus(int i, Matrix &M, int sig) {
         cplx *x = matptr(ntau_ - i);
         herm_matrix_READ_ELEMENT if (sig == -1) M = -M;
      }


template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_matminus(int i, Matrix &M) {
         cplx *x = matptr(ntau_ - i);
         herm_matrix_READ_ELEMENT if (sig_ == -1) M = -M;
      }

template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_gtr_tstp_t(int i, Matrix &M) {
         Matrix M1;
         get_ret_tstp_t(i, M);
         get_les_tstp_t(i, M1);
         M += M1;
      }

template <typename T>
template <class Matrix>
      void herm_matrix_timestep<T>::get_gtr_t_tstp(int i, Matrix &M) {
         Matrix M1;
         get_ret_t_tstp(i, M);
         get_les_t_tstp(i, M1);
         M += M1;
      }

// Note: same syntax like for herm_matrix (with dummy argument tstp).
// The following routines are not very "efficent" but sometimes simple to implement

template <typename T>
      inline void herm_matrix_timestep<T>::get_ret(int i, int j, cplx &x) {
         assert(i == tstp_ || j == tstp_);

         if (i == tstp_) {
            x = *retptr(j);
         } else {
            x = *retptr(i);
            x = -std::conj(x);
         }
      }

template <typename T>
      inline void herm_matrix_timestep<T>::get_les(int i, int j, cplx &x) {
         assert(i == tstp_ || j == tstp_);

         if (j == tstp_) {
            x = *lesptr(i);
         } else {
            x = *lesptr(j);
            x = -std::conj(x);
         }
      }

template <typename T>
      inline void herm_matrix_timestep<T>::get_tv(int i, int j, cplx &x) {
         assert(i == tstp_);
         x = *tvptr(j);
      }

template <typename T>
      inline void herm_matrix_timestep<T>::get_vt(int i, int j, cplx &x) {
         assert(j == tstp_);
         x = *tvptr(ntau_ - i);
         if (sig_ == -1)
            x = std::conj(x);
         else
            x = -std::conj(x);
      }

template <typename T>
      inline void herm_matrix_timestep<T>::get_mat(int i, cplx &x) {
         assert(tstp_ == -1);
         x = *matptr(i);
      }

template <typename T>
      inline void herm_matrix_timestep<T>::get_matminus(int i, cplx &x) {
         assert(tstp_ == -1);
         x = *matptr(ntau_ - i);
         if (sig_ == -1)
            x = -x;
      }


template <typename T>
      std::complex<T> herm_matrix_timestep<T>::density_matrix(int tstp) {
         assert(tstp_ == tstp);
         cplx x1;
         if (tstp_ == -1) {
            get_mat(ntau_, x1);
            return -x1;
         } else {
            get_les(tstp_, tstp_, x1);
            return std::complex<T>(0.0, sig_) * x1;
         }
      }

template <typename T>
      inline void herm_matrix_timestep<T>::density_matrix(int tstp, cplx &rho) {
         assert(tstp_ == tstp);
         cplx x1;
         if (tstp_ == -1) {
            get_mat(ntau_, x1);
            rho = -x1;
         } else {
            get_les(tstp_, tstp_, x1);
            rho = std::complex<T>(0.0, sig_) * x1;
         }
      }


template<typename T> template<class Matrix>
      void herm_matrix_timestep<T>::density_matrix(Matrix &M){

         if(tstp_==-1){
            get_mat(ntau_,M);
            M *= (-1.0);
         }else{
            get_les_tstp_t(tstp_,M);
            M *= std::complex<T>(0.0,1.0*sig_);
         }
      }


template<typename T> template<class Matrix>
      void herm_matrix_timestep<T>::density_matrix(int tstp, Matrix &M){
         assert(tstp == tstp_);
         if(tstp_==-1){
            get_mat(ntau_,M);
            M *= (-1.0);
         }else{
            get_les_tstp_t(tstp_,M);
            M *= std::complex<T>(0.0,1.0*sig_);
         }
      }


/* #######################################################################################
#
#   WRITING ELEMENTS FROM ANY MATRIX TYPE
#   OR FROM COMPLEX NUMBERS (then only the (0,0) element is addressed for dim>0)
#
########################################################################################*/
#define herm_matrix_SET_ELEMENT_MATRIX                                       \
      {                                                                        \
         int r, s;                                                             \
         for (r = 0; r < size1_; r++)                                            \
            for (s = 0; s < size2_; s++)                                        \
               x[r * size2_ + s] = M(r, s);                                    \
         }

template<typename T> template <class Matrix> void herm_matrix_timestep<T>::set_ret(int j,Matrix &M){
         cplx *x=retptr(j);
         herm_matrix_SET_ELEMENT_MATRIX
      }


template<typename T> template <class Matrix> void herm_matrix_timestep<T>::set_ret(int i, int j,Matrix &M){
         assert(i == tstp_);
         assert(j <= i);
         cplx *x=retptr(j);
         herm_matrix_SET_ELEMENT_MATRIX
      }


template<typename T> inline void herm_matrix_timestep<T>::set_ret(int i, int j, cplx &x){
      assert(i == tstp_ && j <= i);
      *retptr(j) = x;
   }


template<typename T> template <class Matrix> void herm_matrix_timestep<T>::set_les(int j,Matrix &M){
      cplx *x=lesptr(j);
      herm_matrix_SET_ELEMENT_MATRIX
   }

template<typename T> template <class Matrix> void herm_matrix_timestep<T>::set_les(int i, int j, Matrix &M){
      assert(j == tstp_);
      assert(i <= j);
      cplx *x=lesptr(i);
      herm_matrix_SET_ELEMENT_MATRIX
   }

template<typename T> inline void herm_matrix_timestep<T>::set_les(int i, int j, cplx &x){
      assert(j == tstp_ && i <= j);
      *lesptr(i) = x;
   }


template<typename T> template <class Matrix> void herm_matrix_timestep<T>::set_tv(int j,Matrix &M){
  cplx *x=tvptr(j);
  herm_matrix_SET_ELEMENT_MATRIX
}


template<typename T> template <class Matrix> void herm_matrix_timestep<T>::set_tv(int i, int j, Matrix &M){
   assert(i == tstp_);
   assert(j <= ntau_);
   cplx *x=tvptr(j);
   herm_matrix_SET_ELEMENT_MATRIX
}


template<typename T> inline void herm_matrix_timestep<T>::set_tv(int i, int j, cplx &x){
   assert(i == tstp_);
   assert(j <= ntau_);
   *tvptr(j) = x;
}


template<typename T> template <class Matrix> void herm_matrix_timestep<T>::set_mat(int i,Matrix &M){
   assert(i <= ntau_ );
   cplx *x=matptr(i);
   herm_matrix_SET_ELEMENT_MATRIX
}

template<typename T> inline void herm_matrix_timestep<T>::set_mat(int i,cplx &x){
   assert(i <= ntau_ );
   *matptr(i) = x;
}















// G(t,t') ==> F(t)G(t,t')   ... ft+t*element_size_ points to F(t)

template <typename T>
void herm_matrix_timestep<T>::left_multiply(std::complex<T> *f0,
   std::complex<T> *ft, T weight) {
   int m;
   cplx *x0, *xtemp, *ftemp;
   xtemp = new cplx[element_size_];
   if (tstp_ == -1) {
      x0 = data_;
      for (m = 0; m <= ntau_; m++) {
         element_mult<T, LARGESIZE>(size1_, xtemp, f0,
            x0 + m * element_size_);
         element_smul<T, LARGESIZE>(size1_, xtemp, weight);
         element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
      }
   } else {
      ftemp = ft + tstp_ * element_size_;
      x0 = data_;
      for (m = 0; m <= tstp_; m++) {
         element_mult<T, LARGESIZE>(size1_, xtemp, ftemp,
            x0 + m * element_size_);
         element_smul<T, LARGESIZE>(size1_, xtemp, weight);
         element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
      }
      x0 = data_ + (tstp_ + 1) * element_size_;
      for (m = 0; m <= ntau_; m++) {
         element_mult<T, LARGESIZE>(size1_, xtemp, ftemp,
            x0 + m * element_size_);
         element_smul<T, LARGESIZE>(size1_, xtemp, weight);
         element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
      }
      x0 = data_ + (tstp_ + 1 + ntau_ + 1) * element_size_;
      for (m = 0; m <= tstp_; m++) {
         element_mult<T, LARGESIZE>(size1_, xtemp, ft + m * element_size_,
            x0 + m * element_size_);
         element_smul<T, LARGESIZE>(size1_, xtemp, weight);
         element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
      }
   }
   delete[] xtemp;
}

template <typename T>
void herm_matrix_timestep<T>::left_multiply(function<T> &ft, T weight) {
   assert( ft.nt() >= tstp_);

   this->left_multiply(ft.ptr(-1), ft.ptr(0), weight);
}

template <typename T>
void herm_matrix_timestep<T>::left_multiply(int tstp, function<T> &ft, T weight) {
   assert(tstp == tstp_);
   assert( ft.nt() >= tstp_);

   this->left_multiply(ft.ptr(-1), ft.ptr(0), weight);
}


template <typename T>
void herm_matrix_timestep<T>::left_multiply_hermconj(int tstp, function<T> &ft,
                                            T weight) {
    assert(tstp == tstp_);
    int m;
    cplx *xtemp, *ftemp, *x0, *f0, *fcc;
    xtemp = new cplx[element_size_];
    fcc = new cplx[element_size_];
    assert(tstp >= -1 && ft.nt() >= tstp &&
           ft.size1() == size1_ && ft.size2() == size2_ &&
      "tstp >= -1 && ft.nt() >= tstp && ft.size1() == size1_ && ft.size2() == size2_");

    f0 = ft.ptr(-1);
    element_conj<T, LARGESIZE>(size1_, fcc, f0);
    if (tstp == -1) {
        x0 = data_;
        for (m = 0; m <= ntau_; m++) {
            element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
                                       x0 + m * element_size_);
            element_smul<T, LARGESIZE>(size1_, xtemp, weight);
            element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
        }
    } else {
        ftemp = ft.ptr(tstp);
        element_conj<T, LARGESIZE>(size1_, fcc, ftemp);
        x0 = data_;
        for (m = 0; m <= tstp; m++) {
            element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
                                       x0 + m * element_size_);
            element_smul<T, LARGESIZE>(size1_, xtemp, weight);
            element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
        }
        x0 = data_ + (tstp_ + 1) * element_size_;
        for (m = 0; m <= ntau_; m++) {
            element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
                                       x0 + m * element_size_);
            element_smul<T, LARGESIZE>(size1_, xtemp, weight);
            element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
        }
        x0 = data_ + (tstp_ + 1 + ntau_ + 1) * element_size_;
        for (m = 0; m <= tstp; m++) {
            element_conj<T, LARGESIZE>(size1_, fcc, ft.ptr(m));
            element_mult<T, LARGESIZE>(size1_, xtemp, fcc,
                                       x0 + m * element_size_);
            element_smul<T, LARGESIZE>(size1_, xtemp, weight);
            element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
        }
    }
    delete[] xtemp;
    delete[] fcc;
}





// G(t,t') ==> G(t,t')F(t')

template <typename T>
void herm_matrix_timestep<T>::right_multiply(std::complex<T> *f0,
   std::complex<T> *ft, T weight) {
   int m;
   cplx *xtemp, *ftemp, *x0;
   xtemp = new cplx[element_size_];
   if (tstp_ == -1) {
      x0 = data_;
      for (m = 0; m <= ntau_; m++) {
         element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
            f0);
         element_smul<T, LARGESIZE>(size1_, xtemp, weight);
         element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
      }
   } else {
      x0 = data_;
      for (m = 0; m <= tstp_; m++) {
         element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
            ft + m * element_size_);
         element_smul<T, LARGESIZE>(size1_, xtemp, weight);
         element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
      }
      x0 = data_ + (tstp_ + 1) * element_size_;
      for (m = 0; m <= ntau_; m++) {
         element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
            f0);
         element_smul<T, LARGESIZE>(size1_, xtemp, weight);
         element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
      }
      ftemp = ft + tstp_ * element_size_;
      x0 = data_ + (tstp_ + 1 + ntau_ + 1) * element_size_;
      for (m = 0; m <= tstp_; m++) {
         element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
            ftemp);
         element_smul<T, LARGESIZE>(size1_, xtemp, weight);
         element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
      }
   }
   delete[] xtemp;
}

template <typename T>
void herm_matrix_timestep<T>::right_multiply(function<T> &ft, T weight) {
   assert(ft.nt() >= tstp_);

   this->right_multiply(ft.ptr(-1), ft.ptr(0), weight);
}


template <typename T>
void herm_matrix_timestep<T>::right_multiply(int tstp, function<T> &ft, T weight) {
   assert(tstp == tstp_);
   assert(ft.nt() >= tstp_);

   this->right_multiply(ft.ptr(-1), ft.ptr(0), weight);
}



template <typename T>
void herm_matrix_timestep<T>::right_multiply_hermconj(int tstp, function<T> &ft,
                                             T weight) {
    assert(tstp == tstp_);
    int m;
    cplx *xtemp, *ftemp, *x0, *f0, *fcc;
    xtemp = new cplx[element_size_];
    fcc = new cplx[element_size_];
    assert(ft.size1() == size1_ && "ft.size1() == size1_");
    assert(ft.nt() >= tstp && "ft.nt() >= tstp");
    assert(tstp >= -1&& "ft.nt() >= tstp");
    f0 = ft.ptr(-1);
    element_conj<T, LARGESIZE>(size1_, fcc, f0);
    if (tstp == -1) {
        x0 = data_;
        for (m = 0; m <= ntau_; m++) {
            element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
                                       fcc);
            element_smul<T, LARGESIZE>(size1_, xtemp, weight);
            element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
        }
    } else {
        x0 = data_;
        for (m = 0; m <= tstp; m++) {
            element_conj<T, LARGESIZE>(size1_, fcc, ft.ptr(m));
            element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
                                       fcc);
            element_smul<T, LARGESIZE>(size1_, xtemp, weight);
            element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
        }
        x0 = data_ + (tstp_ + 1) * element_size_;
        element_conj<T, LARGESIZE>(size1_, fcc, f0);
        for (m = 0; m <= ntau_; m++) {
            element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
                                       fcc);
            element_smul<T, LARGESIZE>(size1_, xtemp, weight);
            element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
        }
        ftemp = ft.ptr(tstp);
        element_conj<T, LARGESIZE>(size1_, fcc, ftemp);
        x0 = data_ + (tstp_ + 1 + ntau_ + 1) * element_size_;
        for (m = 0; m <= tstp; m++) {
            element_mult<T, LARGESIZE>(size1_, xtemp, x0 + m * element_size_,
                                       fcc);
            element_smul<T, LARGESIZE>(size1_, xtemp, weight);
            element_set<T, LARGESIZE>(size1_, x0 + m * element_size_, xtemp);
        }
    }
    delete[] fcc;
    delete[] xtemp;
}


template <typename T>
void herm_matrix_timestep<T>::incr(herm_matrix_timestep<T> &g1, T weight) {
   int m;
   assert(g1.size1_ == size1_ && g1.ntau_ == ntau_ && g1.tstp_ == tstp_);
   for (m = 0; m < total_size_; m++)
      data_[m] += weight * g1.data_[m];
}


template <typename T>
void herm_matrix_timestep<T>::incr_timestep(int tstp, herm_matrix_timestep<T> &g1, T weight) {
   assert(tstp == tstp_);
   assert(g1.size1_ == size1_ && g1.ntau_ == ntau_ && g1.tstp_ == tstp_);
   for (int m = 0; m < total_size_; m++)
      data_[m] += weight * g1.data_[m];
}



#define HERM_MATRIX_INCR_TSTP                                                \
if (alpha == cplx(1.0, 0.0)) {                                           \
   for (i = 0; i < len; i++)                                            \
      x0[i] += x[i];                                                   \
} else {                                                                 \
   for (i = 0; i < len; i++)                                            \
      x0[i] += alpha * x[i];                                           \
}

template <typename T>
void herm_matrix_timestep<T>::incr(herm_matrix<T> &g, T alpha) {
   int i, len;
   cplx *x, *x0;
   assert(tstp_ <= g.nt() && ntau_ == g.ntau() && size1_ == g.size1());
   if (tstp_ == -1) {
      len = (ntau_ + 1) * element_size_;
      x = g.matptr(0);
      x0 = data_;
      HERM_MATRIX_INCR_TSTP
   } else {
      len = (tstp_ + 1) * element_size_;
      x = g.retptr(tstp_, 0);
      x0 = data_;
      HERM_MATRIX_INCR_TSTP
      len = (ntau_ + 1) * element_size_;
      x = g.tvptr(tstp_, 0);
      x0 = data_ + (tstp_ + 1) * element_size_;
      HERM_MATRIX_INCR_TSTP
      len = (tstp_ + 1) * element_size_;
      x = g.lesptr(0, tstp_);
      x0 = data_ + (tstp_ + 1 + ntau_ + 1) * element_size_;
      HERM_MATRIX_INCR_TSTP
   }
}

template <typename T>
void herm_matrix_timestep<T>::incr_timestep(int tstp, herm_matrix<T> &g, T alpha) {
   assert(tstp == tstp_);
   assert(tstp_ <= g.nt() && ntau_ == g.ntau() && size1_ == g.size1());
   int i, len;
   cplx *x, *x0;
   if (tstp_ == -1) {
      len = (ntau_ + 1) * element_size_;
      x = g.matptr(0);
      x0 = data_;
      HERM_MATRIX_INCR_TSTP
   } else {
      len = (tstp_ + 1) * element_size_;
      x = g.retptr(tstp_, 0);
      x0 = data_;
      HERM_MATRIX_INCR_TSTP
      len = (ntau_ + 1) * element_size_;
      x = g.tvptr(tstp_, 0);
      x0 = data_ + (tstp_ + 1) * element_size_;
      HERM_MATRIX_INCR_TSTP
      len = (tstp_ + 1) * element_size_;
      x = g.lesptr(0, tstp_);
      x0 = data_ + (tstp_ + 1 + ntau_ + 1) * element_size_;
      HERM_MATRIX_INCR_TSTP
   }
}


#undef HERM_MATRIX_INCR_TSTP

template <typename T>
void herm_matrix_timestep<T>::smul(T weight) {
   herm_matrix_timestep_view<T> g_view(*this);

   g_view.smul(weight);

}


template <typename T>
void herm_matrix_timestep<T>::smul(int tstp, T weight) {
   assert(tstp == tstp_);
   herm_matrix_timestep_view<T> g_view(*this);

   g_view.smul(weight);

}



template <typename T>
void herm_matrix_timestep<T>::set_matrixelement(int i1, int i2,
   herm_matrix_timestep<T> &g, int j1, int j2) {
   herm_matrix_timestep_view<T> tmp(g);
   herm_matrix_timestep_view<T> tmp1(*this);
   tmp1.set_matrixelement(i1, i2, tmp, j1, j2);
}


template <typename T>
void herm_matrix_timestep<T>::set_matrixelement(int tstp, int i1, int i2,
   herm_matrix_timestep<T> &g, int j1, int j2) {
   assert(tstp == tstp_);
   herm_matrix_timestep_view<T> tmp(g);
   herm_matrix_timestep_view<T> tmp1(*this);
   tmp1.set_matrixelement(i1, i2, tmp, j1, j2);
}


template <typename T>
void herm_matrix_timestep<T>::set_matrixelement(
   int i1, int i2, herm_matrix_timestep_view<T> &g, int j1, int j2) {
   herm_matrix_timestep_view<T> tmp1(*this);
   tmp1.set_matrixelement(i1, i2, g, j1, j2);
}


template <typename T>
void herm_matrix_timestep<T>::set_matrixelement(int tstp,
   int i1, int i2, herm_matrix_timestep_view<T> &g, int j1, int j2) {
   herm_matrix_timestep_view<T> tmp1(*this);
   tmp1.set_matrixelement(i1, i2, g, j1, j2);
}


template <typename T>
void herm_matrix_timestep<T>::set_matrixelement(int i1, int i2,
   herm_matrix<T> &g, int j1,
   int j2) {
   herm_matrix_timestep_view<T> tmp(tstp_, g);
   herm_matrix_timestep_view<T> tmp1(*this);
   tmp1.set_matrixelement(i1, i2, tmp, j1, j2);
}


template <typename T>
void herm_matrix_timestep<T>::set_matrixelement(int tstp, int i1, int i2,
   herm_matrix<T> &g, int j1,
   int j2) {
   assert(tstp == tstp_);
   herm_matrix_timestep_view<T> tmp(tstp_, g);
   herm_matrix_timestep_view<T> tmp1(*this);
   tmp1.set_matrixelement(i1, i2, tmp, j1, j2);
}

template<typename T> void herm_matrix_timestep<T>::set_submatrix(int tstp, std::vector<int> &i1,
  std::vector<int> &i2,herm_matrix<T> &g,std::vector<int> &j1,std::vector<int> &j2){
    assert(tstp == tstp_);
    assert(tstp <= g.nt());
    assert(i1.size()==i2.size() && i1.size()==j1.size() && j1.size()==j2.size());
    assert(size1_*size2_==i1.size());
    for(int k1=0;k1<i1.size();k1++){
        this->set_matrixelement(tstp,i1[k1],i2[k1],g,j1[k1],j2[k1]);
    }
}


template<typename T> void herm_matrix_timestep<T>::set_submatrix(int tstp, std::vector<int> &i1,
  std::vector<int> &i2,herm_matrix_timestep<T> &g,std::vector<int> &j1,std::vector<int> &j2){
    assert(tstp == tstp_);
    assert(tstp == g.tstp());
    assert(i1.size()==i2.size() && i1.size()==j1.size() && j1.size()==j2.size());
    assert(size1_*size2_==i1.size());
    for(int k1=0;k1<i1.size();k1++){
        this->set_matrixelement(tstp,i1[k1],i2[k1],g,j1[k1],j2[k1]);
    }
}

/* #########################################################################
#
#   MPI UTILS
#
############################################################################ */


#if CNTR_USE_MPI == 1


template <typename T>
void herm_matrix_timestep<T>::Reduce_timestep(int root) {
   int taskid;
   int len = 2 * (2 * (tstp_ + 1) + ntau_ + 1) * element_size_;
   MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
   if (sizeof(T) == sizeof(double)) {
      if (taskid == root) {
         MPI_Reduce(MPI_IN_PLACE, (double *)this->data_, len, MPI_DOUBLE_PRECISION, MPI_SUM, root, MPI_COMM_WORLD);
      } else {
         MPI_Reduce((double *)this->data_,
            (double *)this->data_, len, MPI_DOUBLE_PRECISION, MPI_SUM, root, MPI_COMM_WORLD);
      }
   } else {
      std::cerr << "herm_matrix_timestep<T>::MPI_Reduce only for double "
      << std::endl;
      exit(0);
   }
}


template <typename T>
void herm_matrix_timestep<T>::Reduce_timestep(int tstp, int root) {
   assert(tstp == tstp_);
   int taskid;
   int len = 2 * (2 * (tstp_ + 1) + ntau_ + 1) * element_size_;
   MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
   if (sizeof(T) == sizeof(double)) {
      if (taskid == root) {
         MPI_Reduce(MPI_IN_PLACE, (double *)this->data_, len,
            MPI_DOUBLE_PRECISION, MPI_SUM, root, MPI_COMM_WORLD);
      } else {
         MPI_Reduce((double *)this->data_,
            (double *)this->data_, len, MPI_DOUBLE_PRECISION,
            MPI_SUM, root, MPI_COMM_WORLD);
      }
   } else {
      std::cerr << "herm_matrix_timestep<T>::MPI_Reduce only for double "
      << std::endl;
      exit(0);
   }
}



template <typename T>
void herm_matrix_timestep<T>::Bcast_timestep(int tstp, int ntau, int size1,
   int root) {
   int numtasks,taskid;
   MPI_Comm_size(MPI_COMM_WORLD, &numtasks);
   MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
   if (taskid != root)
      resize(tstp, ntau, size1);
   int len = (2 * (tstp_ + 1) + ntau_ + 1) * element_size_;
// test effective on root:
   assert(tstp == tstp_);
   assert(ntau == ntau_);
   assert(size1 == size1_);
   if (sizeof(T) == sizeof(double))
      MPI_Bcast(data_, len, MPI_DOUBLE_COMPLEX, root, MPI_COMM_WORLD);
   else
      MPI_Bcast(data_, len, MPI_COMPLEX, root, MPI_COMM_WORLD);
}

template <typename T>
void herm_matrix_timestep<T>::Bcast_timestep(int tstp, int root) {
   int numtasks,taskid;
   MPI_Comm_size(MPI_COMM_WORLD, &numtasks);
   MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
   if (taskid != root)
      resize(tstp, ntau_, size1_);
   int len = (2 * (tstp_ + 1) + ntau_ + 1) * element_size_;
// test effective on root:
   assert(tstp == tstp_);
   if (sizeof(T) == sizeof(double))
      MPI_Bcast(data_, len, MPI_DOUBLE_COMPLEX, root, MPI_COMM_WORLD);
   else
      MPI_Bcast(data_, len, MPI_COMPLEX, root, MPI_COMM_WORLD);
}

template <typename T>
void herm_matrix_timestep<T>::Send_timestep(int tstp, int ntau, int size1,
   int dest, int tag) {
   int taskid;
   MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
   int len = (2 * (tstp_ + 1) + ntau_ + 1) * element_size_;
   if (!(taskid == dest)) {
      assert(tstp == tstp_);
      assert(ntau == ntau_);
      assert(size1 == size1_);
      if (sizeof(T) == sizeof(double))
         MPI_Send(data_, len, MPI_DOUBLE_COMPLEX, dest, tag, MPI_COMM_WORLD);
      else
         MPI_Send(data_, len, MPI_COMPLEX, dest, tag, MPI_COMM_WORLD);
   } else {
   }
}

template <typename T>
void herm_matrix_timestep<T>::Send_timestep(int tstp, int dest, int tag) {
   int taskid;
   MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
   int len = (2 * (tstp_ + 1) + ntau_ + 1) * element_size_;
   if (!(taskid == dest)) {
      assert(tstp == tstp_);
      if (sizeof(T) == sizeof(double))
         MPI_Send(data_, len, MPI_DOUBLE_COMPLEX, dest, tag, MPI_COMM_WORLD);
      else
         MPI_Send(data_, len, MPI_COMPLEX, dest, tag, MPI_COMM_WORLD);
   } else {
   }
}

template <typename T>
void herm_matrix_timestep<T>::Recv_timestep(int tstp, int ntau, int size1,
   int root, int tag) {
   int taskid;
   MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
   if (!(taskid == root)) {
      resize(tstp, ntau, size1);
      int len = (2 * (tstp_ + 1) + ntau_ + 1) * element_size_;
      if (sizeof(T) == sizeof(double))
         MPI_Recv(data_, len, MPI_DOUBLE_COMPLEX, root, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
      else
         MPI_Recv(data_, len, MPI_COMPLEX, root, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
   }
}


template <typename T>
void herm_matrix_timestep<T>::Recv_timestep(int tstp, int root, int tag) {
   int taskid;
   MPI_Comm_rank(MPI_COMM_WORLD, &taskid);
   if (!(taskid == root)) {
      resize(tstp, ntau_, size1_);
      int len = (2 * (tstp_ + 1) + ntau_ + 1) * element_size_;
      if (sizeof(T) == sizeof(double))
         MPI_Recv(data_, len, MPI_DOUBLE_COMPLEX, root, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
      else
         MPI_Recv(data_, len, MPI_COMPLEX, root, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
   }
}


#endif

#undef herm_matrix_READ_ELEMENT
#undef herm_matrix_READ_ELEMENT_MINUS_CONJ

/* ##########################################################################
#
#   HDF5 I/O (using implementation in herm_matrix_timestep_view)
#
#############################################################################*/

#if CNTR_USE_HDF5 == 1

template <typename T>
void herm_matrix_timestep<T>::write_to_hdf5(hid_t group_id) {
   herm_matrix_timestep_view<T> tmp(*this);
   tmp.write_to_hdf5(group_id);
}

template <typename T>
void herm_matrix_timestep<T>::write_to_hdf5(hid_t group_id,
   const char *groupname) {
   herm_matrix_timestep_view<T> tmp(*this);
   tmp.write_to_hdf5(group_id, groupname);
}

template <typename T>
void herm_matrix_timestep<T>::write_to_hdf5(const char *filename,
   const char *groupname) {
   herm_matrix_timestep_view<T> tmp(*this);
   tmp.write_to_hdf5(filename, groupname);
}

template <typename T>
void herm_matrix_timestep<T>::read_from_hdf5(hid_t group_id) {
// -- implementation is different from herm_matrix_timestep_view, because G can be resized:
   int tstp = read_primitive_type<int>(group_id, "tstp");
   int ntau = read_primitive_type<int>(group_id, "ntau");
   int sig = read_primitive_type<int>(group_id, "sig");
   int size1 = read_primitive_type<int>(group_id, "size1");
// RESIZE G
   this->resize(tstp, ntau, size1);
   sig_ = sig;
   herm_matrix_timestep_view<T> tmp(*this);
   tmp.read_from_hdf5(group_id);
}

template <typename T>
void herm_matrix_timestep<T>::read_from_hdf5(hid_t group_id,
   const char *groupname) {
   hid_t sub_group_id = open_group(group_id, groupname);
   this->read_from_hdf5(sub_group_id);
   close_group(sub_group_id);
}

template <typename T>
void herm_matrix_timestep<T>::read_from_hdf5(const char *filename,
   const char *groupname) {
   hid_t file_id = read_hdf5_file(filename);
   this->read_from_hdf5(file_id, groupname);
   close_hdf5_file(file_id);
}
#endif

} // namespace cntr

#endif  // CNTR_HERM_MATRIX_TIMESTEP_IMPL_H