cntr_herm_matrix_decl.hpp

Go to the documentation of this file.

#ifndef CNTR_HERM_MATRIX_DECL_H
#define CNTR_HERM_MATRIX_DECL_H

#include "cntr_global_settings.hpp"

namespace cntr {

template <typename T> class function;
template <typename T> class herm_matrix_timestep;

template <typename T>
class herm_matrix {
  public:
    typedef std::complex<T> cplx;
    typedef T scalar_type;

    /* construction, destruction */
    herm_matrix();
    ~herm_matrix();
    herm_matrix(int nt, int ntau, int size1 = 1, int sig = -1);
    herm_matrix(int nt, int ntau, int size1, int size2, int sig);
    herm_matrix(const herm_matrix &g);
    herm_matrix &operator=(const herm_matrix &g);
#if __cplusplus >= 201103L
    herm_matrix(herm_matrix &&g) noexcept;
    herm_matrix &operator=(herm_matrix &&g) noexcept;
#endif
    /* resize */
    void resize_discard(int nt, int ntau, int size1);
    void resize_nt(int nt);
    void resize(int nt, int ntau, int size1);
    void clear(void);
    /* access size etc ... */
    int element_size(void) const { return element_size_; }
    int size1(void) const { return size1_; }
    int size2(void) const { return size2_; }
    int ntau(void) const { return ntau_; }
    int nt(void) const { return nt_; }
    int sig(void) const { return sig_; }
    void set_sig(int sig) { sig_ = sig; }
    /* conversion from other types */
    // herm_matrix<T> & herm_matrix(const matrix<T> &g1);
    // herm_matrix<T> & herm_matrix(const scalar<T> &g1);
    // raw pointer to elements ... to be used with care
    inline cplx *lesptr(int i, int j);
    inline cplx *retptr(int i, int j);
    inline cplx *tvptr(int i, int j);
    inline cplx *matptr(int i);
    inline const cplx *lesptr(int i, int j) const;
    inline const cplx *retptr(int i, int j) const;
    inline const cplx *tvptr(int i, int j) const;
    inline const cplx *matptr(int i) const;
    // reading basic and derived elements to any Matrix type
    template <class Matrix>
    void get_les(int i, int j, Matrix &M) const;
    template <class Matrix>
    void get_gtr(int i, int j, Matrix &M) const;
    template <class Matrix>
    void get_ret(int i, int j, Matrix &M) const;
    template <class Matrix>
    void get_vt(int i, int j, Matrix &M) const;
    template <class Matrix>
    void get_tv(int i, int j, Matrix &M) const;
    template <class Matrix>
    void get_mat(int i, Matrix &M) const;
    template <class Matrix>
    void get_matminus(int i, Matrix &M) const;
    // reading complex numbers:
    // these will adress only (0,0) element for dim>1:
    inline void get_les(int i, int j, cplx &x) const;
    inline void get_gtr(int i, int j, cplx &x) const;
    inline void get_ret(int i, int j, cplx &x) const;
    inline void get_vt(int i, int j, cplx &x) const;
    inline void get_tv(int i, int j, cplx &x) const;
    inline void get_mat(int i, cplx &x) const;
    inline void get_matminus(int i, cplx &x) const;
    cplx density_matrix(int i) const;
    // should work with eigen
    template <class Matrix>
    void density_matrix(int tstp, Matrix &M) const;

    // writing basic elements:
    template <class Matrix>
    void set_les(int i, int j, Matrix &M);
    template <class Matrix>
    void set_ret(int i, int j, Matrix &M);
    template <class Matrix>
    void set_tv(int i, int j, Matrix &M);
    template <class Matrix>
    void set_mat(int i, Matrix &M);
    template<class Matrix> void set_mat_real(int i,Matrix &M);
    void set_mat_herm(int i);
    void set_mat_herm(void);
    inline void set_les(int i, int j, cplx x);
    inline void set_ret(int i, int j, cplx x);
    inline void set_tv(int i, int j, cplx x);
    inline void set_mat(int i, cplx x);
    // INPUT/OUTPUT
    void print_to_file(const char *file, int precision = 16) const;
    void read_from_file(const char *file);
    // read up to timestep nt1: NO RESIZE
    void read_from_file(int nt1, const char *file);
// HDF5 I/O
#if CNTR_USE_HDF5 == 1
    void write_to_hdf5(hid_t group_id);
    void write_to_hdf5(hid_t group_id, const char *groupname);
    void write_to_hdf5(const char *filename, const char *groupname);
    void read_from_hdf5(hid_t group_id);
    void read_from_hdf5(hid_t group_id, const char *groupname);
    void read_from_hdf5(const char *filename, const char *groupname);
    // read up to timestep nt1: NO RESIZE
    void read_from_hdf5(int nt1, hid_t group_id);
    void read_from_hdf5(int nt1, hid_t group_id, const char *groupname);
    void read_from_hdf5(int nt1, const char *filename, const char *groupname);
    // writing some other compressed format: just a collection of timeslices
    // -1,0,dt,2*dt,...
    void write_to_hdf5_slices(hid_t group_id, int dt);
    void write_to_hdf5_slices(hid_t group_id, const char *groupname, int dt);
    void write_to_hdf5_slices(const char *filename, const char *groupname,
                              int dt);
    void write_to_hdf5_tavtrel(hid_t group_id, int dt);
    void write_to_hdf5_tavtrel(hid_t group_id, const char *groupname, int dt);
    void write_to_hdf5_tavtrel(const char *filename, const char *groupname,
                               int dt);
#endif
    // simple operations on the timesteps (n<0: Matsubara)
    void set_timestep_zero(int tstp);
    void set_timestep(int tstp, herm_matrix &g1);
    void set_timestep(int tstp, herm_matrix_timestep<T> &timestep);
    void get_timestep(int tstp, herm_matrix_timestep<T> &timestep) const;
    void get_timestep(int tstp, herm_matrix<T> &timestep) const;
    void incr_timestep(int tstp, herm_matrix_timestep<T> &timestep,
                       cplx alpha);
    void incr_timestep(int tstp, herm_matrix_timestep<T> &timestep);
    void incr_timestep(int tstp, herm_matrix<T> &g, cplx alpha);
    void incr_timestep(int tstp, herm_matrix<T> &g);
    // check what this is actually doing
    void incr(herm_matrix<T> &g, cplx alpha);
    void incr(herm_matrix<T> &g);

    void set_matrixelement(int tstp, int i1, int i2,
                           herm_matrix_timestep<T> &timestep);
    void set_matrixelement(int tstp, int i1, int i2, herm_matrix<T> &g);
    // this(i1,i2) <- g(j1,j2)
    void set_matrixelement(int tstp, int i1, int i2,
                           herm_matrix_timestep<T> &g, int j1, int j2);
    void set_matrixelement(int tstp, int i1, int i2, herm_matrix<T> &g,
                           int j1, int j2);
    void set_matrixelement(int i1,int i2,herm_matrix<T> &g);
    void set_matrixelement(int i1,int i2,herm_matrix<T> &g,int j1,int j2);

    void set_submatrix(std::vector<int> &i1,std::vector<int> &i2,
      herm_matrix<T> &g,std::vector<int> &j1,std::vector<int> &j2);
    void set_submatrix(int tstp, std::vector<int> &i1,std::vector<int> &i2,
      herm_matrix<T> &g,std::vector<int> &j1,std::vector<int> &j2);

    // multiplication with function

    // check if multiplication with pointer can be removed (there might be
    //  internal dependencies)
    void
    left_multiply(int tstp, cplx *f0, cplx *ft,
                  T weight = 1.0); // version with function object is safer
    void
    right_multiply(int tstp, cplx *f0, cplx *ft,
                   T weight = 1.0); // version with function object is safer
    void left_multiply(int tstp, function<T> &ft, T weight = 1.0);
    // - check why underscore (private routine?)
    // - check what std::integral_constant is doing
    // - add interface to multiply with a matrix from left/right
    template <int SIZE>
    void left_multiply_(int tstp, function<T> &ft, T weight,
                        std::integral_constant<int, SIZE>);
    void right_multiply(int tstp, function<T> &ft, T weight = 1.0);
    template <int SIZE>
    void right_multiply_(int tstp, function<T> &ft, T weight,
                         std::integral_constant<int, SIZE>);
    void left_multiply_hermconj(
        int tstp, function<T> &ft,
        T weight = 1.0); // G(t,t') -> weight*f(t)^dag * G(t,t')
    void right_multiply_hermconj(
        int tstp, function<T> &ft,
        T weight = 1.0); // G(t,t') -> weight*G(t,t') * f(t')^dag
    void smul(int tstp, T weight);
    void smul(int tstp, cplx weight);
// MPI UTILS
#if CNTR_USE_MPI == 1
    void Reduce_timestep(int tstp, int root);
    void Bcast_timestep(int tstp, int root);
    void Send_timestep(int tstp, int dest, int tag);
    void Recv_timestep(int tstp, int root, int tag);
#endif
  private:
    int les_offset(int t, int t1) const;
    int ret_offset(int t, int t1) const;
    int tv_offset(int t, int tau) const;
    int mat_offset(int tau) const;

  private:

    cplx *les_;

    cplx *ret_;

    cplx *tv_; 

    cplx *mat_;

    int nt_;

    int ntau_;

    int size1_;

    int size2_;

    int element_size_;

    int sig_; // Bose = +1, Fermi =-1
};

}  // namespace cntr

#endif  // CNTR_HERM_MATRIX_DECL_H