field_F.h

Go to the documentation of this file.

#ifndef FIELD_F_INCLUDED
#define FIELD_F_INCLUDED

#include "commonParameters.h"
#include "communicator.h"
#include "field.h"
#include "field_G.h"
#include "vec_SU_N.h"
#include "gammaMatrix.h"


class Field_F : public Field {
 private:
  int m_Nc;   // num of the color elements
  int m_Nc2;  // num of the double color elements
  int m_Nvol; // lattice volume
  int m_Nd;   // num of the spinor elements
  int m_Nex;  // num of the extra indices

  int myindex(const int c2, const int s, const int site, const int ex)
  const
  {
    return Field::myindex(c2 + m_Nc2 * s, site, ex);
  }

 public:
  Field_F(const int Nvol = CommonParameters::Nvol(), const int Nex = 1) :
    Field(), m_Nc(CommonParameters::Nc()),
    m_Nvol(Nvol),
    m_Nd(CommonParameters::Nd()),
    m_Nex(Nex)
  {
    m_Nc2 = 2 * m_Nc;
    reset(m_Nc2 * m_Nd, m_Nvol, m_Nex);
    check();
  }

  Field_F(const Field& x) :
    Field(x),
    m_Nc(CommonParameters::Nc()),
    m_Nc2(2 * m_Nc),
    m_Nvol(x.nvol()),
    m_Nd(CommonParameters::Nd()),
    m_Nex(x.nex())
  {
  }

  int nc() const { return m_Nc; }
  int nc2() const { return m_Nc2; }
  int nd() const { return m_Nd; }

  double cmp_r(const int cc, const int s, const int site, const int e = 0)
  const
  {
    return field[myindex(2 * cc, s, site, e)];
  }

  double cmp_i(const int cc, const int s, const int site, const int e = 0)
  const
  {
    return field[myindex(2 * cc + 1, s, site, e)];
  }

  void set_r(const int cc, const int s, const int site, const int e,
             const double re)
  {
    field[myindex(2 * cc, s, site, e)] = re;
  }

  void set_i(const int cc, const int s, const int site, const int e,
             const double im)
  {
    field[myindex(2 * cc + 1, s, site, e)] = im;
  }

  void set_ri(const int cc, const int s, const int site, const int e,
              const double re, const double im)
  {
    field[myindex(2 * cc, s, site, e)]     = re;
    field[myindex(2 * cc + 1, s, site, e)] = im;
  }

  Vec_SU_N vec(const int s, const int site, const int e = 0) const
  {
    Vec_SU_N Tmp;

    for (int cc = 0; cc < m_Nc; ++cc) {
      Tmp.set(cc, field[myindex(2 * cc, s, site, e)],
              field[myindex(2 * cc + 1, s, site, e)]);
    }
    return Tmp;
  }

  void set_vec(const int s, const int site, const int e, const Vec_SU_N& F)
  {
    for (int cc = 0; cc < m_Nc; ++cc) {
      field[myindex(2 * cc, s, site, e)]     = F.r(cc);
      field[myindex(2 * cc + 1, s, site, e)] = F.i(cc);
    }
  }

  void add_vec(const int s, const int site, const int e, const Vec_SU_N& F)
  {
    for (int cc = 0; cc < m_Nc; ++cc) {
      field[myindex(2 * cc, s, site, e)]     += F.r(cc);
      field[myindex(2 * cc + 1, s, site, e)] += F.i(cc);
    }
  }

  void clear_vec(const int s, const int site, const int e)
  {
    for (int cc = 0; cc < m_Nc2; ++cc) {
      field[myindex(cc, s, site, e)] = 0.0;
    }
  }

  void xI()
  {
    double tmp_r, tmp_i;

    for (int ex = 0; ex < m_Nex; ++ex) {
      for (int site = 0; site < m_Nvol; ++site) {
        for (int s = 0; s < m_Nd; ++s) {
          for (int cc = 0; cc < m_Nc; ++cc) {
            tmp_r = field[myindex(2 * cc, s, site, ex)];
            tmp_i = field[myindex(2 * cc + 1, s, site, ex)];
            field[myindex(2 * cc, s, site, ex)]     = -tmp_i;
            field[myindex(2 * cc + 1, s, site, ex)] = tmp_r;
          }
        }
      }
    }
  }

  void Ix(const Field_F& w)
  {
    for (int ex = 0; ex < m_Nex; ++ex) {
      for (int site = 0; site < m_Nvol; ++site) {
        for (int s = 0; s < m_Nd; ++s) {
          for (int cc = 0; cc < m_Nc; ++cc) {
            field[myindex(2 * cc, s, site, ex)]     = -w.cmp_i(cc, s, site, ex);
            field[myindex(2 * cc + 1, s, site, ex)] = w.cmp_r(cc, s, site, ex);
          }
        }
      }
    }
  }

  void mult_Field_Gn(int ex, const Field_G&, int ex1,
                     const Field_F&, int ex2);

  void mult_Field_Gd(int ex, const Field_G&, int ex1,
                     const Field_F&, int ex2);

  void multadd_Field_Gn(int ex, const Field_G&, int ex1,
                        const Field_F&, int ex2, double);

  void multadd_Field_Gd(int ex, const Field_G&, int ex1,
                        const Field_F&, int ex2, double);

  void mult_GM(const GammaMatrix&, const Field_F&);

  void mult_iGM(const GammaMatrix&, const Field_F&);

  void mult_GMproj(int, const GammaMatrix&, const Field_F&);

  void mult_GMproj2(int, const GammaMatrix&, const Field_F&);

  double operator*(const Field_F&);

  // check several assumptions for performance implementation.
  void check();

  template<typename T>
  Field_F& operator=(const T& rhs)
  {
    *this = rhs.eval();
    return *this;
  }

  template<typename T>
  Field_F& operator+=(const T& rhs)
  {
    *this += rhs.eval();
    return *this;
  }

  template<typename T>
  Field_F& operator-=(const T& rhs)
  {
    *this -= rhs.eval();
    return *this;
  }

  /*
  */

  Field_F& operator-();
  Field_F& operator=(const double&);
  Field_F& operator+=(const Field_F&);
  Field_F& operator-=(const Field_F&);
  Field_F& operator*=(const double&);
  Field_F& operator*=(const dcomplex&);
  Field_F& operator/=(const double&);
  Field_F& operator/=(const dcomplex&);
};

inline Field_F& Field_F::operator=(const double& r)
{
  field = r;
  return *this;
}


inline Field_F& Field_F::operator-()
{
  field = -field;
  return *this;
}


inline Field_F& Field_F::operator+=(const Field_F& rhs)
{
  field += rhs.field;
  return *this;
}


inline Field_F& Field_F::operator-=(const Field_F& rhs)
{
  field -= rhs.field;
  return *this;
}


inline Field_F& Field_F::operator*=(const double& rhs)
{
  field *= rhs;
  return *this;
}


inline Field_F& Field_F::operator/=(const double& rhs)
{
  field /= rhs;
  return *this;
}


inline double Field_F::operator*(const Field_F& rhs)
{
  double a = (field * rhs.field).sum();
  double b = Communicator::reduce_sum(a);

  return b;
  //  return (field*rhs.field).sum();
}
#endif