staples.cpp

Go to the documentation of this file.

#include "staples.h"

#include "shiftField_lex.h"

//====================================================================
Staples::Staples()
{
  m_vl = CommonParameters::Vlevel();
  m_shift = new ShiftField_lex;
}

//====================================================================
Staples::~Staples()
{
  delete m_shift;
}

//====================================================================
double Staples::plaquette(const Field_G& U)
{
  return (plaq_s(U) + plaq_t(U)) / 2;
}

//====================================================================
double Staples::plaq_s(const Field_G& U)
{
  int Nc       = CommonParameters::Nc();
  int Lvol     = CommonParameters::Lvol();
  int Ndim     = CommonParameters::Ndim();
  int Ndim_spc = Ndim - 1;

  double plaq = 0.0;

  Field_G Umu;

  for (int mu = 0; mu < Ndim_spc; ++mu) {
    int nu = (mu + 1) % Ndim_spc;

    Umu.setpart_ex(0, U, mu);

    upper(m_staple, U, mu, nu);
    //    lower(m_staple,U,mu,nu);

    //    #pragma omp parallel
    {
      //  double plaq1 = m_staple.ddotc_ex(0, U, mu);
      double plaq1 = real(dotc(m_staple, Umu));
#pragma omp master
      plaq += plaq1;
    }
  }

  return plaq / (Lvol * Nc * Ndim_spc);
}

//====================================================================
double Staples::plaq_t(const Field_G& U)
{
  int Nc       = CommonParameters::Nc();
  int Lvol     = CommonParameters::Lvol();
  int Ndim     = CommonParameters::Ndim();
  int Ndim_spc = Ndim - 1;

  double plaq = 0.0;

  Field_G Umu;

  int mu = Ndim_spc;

  Umu.setpart_ex(0, U, mu);

  for (int nu = 0; nu < Ndim_spc; ++nu) {
    upper(m_staple, U, mu, nu);
    //    lower(m_staple, U, mu, nu);

    //    #pragma omp parallel
    {
      //   double plaq1 = m_staple.ddotc_ex(0, U, mu);
      double plaq1 = real(dotc(m_staple, Umu));
#pragma omp master
      plaq += plaq1;
    }
  }

  return plaq / (Lvol * Nc * Ndim_spc);
}

//====================================================================
void Staples::staple(Field_G& W, const Field_G& U, const int mu)
{
  int Ndim = CommonParameters::Ndim();

  W = 0.0;

  for (int nu = 0; nu < Ndim; ++nu) {
    if (nu != mu) {
      W += upper(U, mu, nu);
      W += lower(U, mu, nu);
    }
  }

}

//====================================================================
Field_G Staples::upper(const Field_G& U, const int mu, const int nu)
{
  Field_G c;

  upper(c, U, mu, nu);

  return c;
}

//====================================================================
Field_G Staples::lower(const Field_G& U, const int mu, const int nu)
{
  Field_G c;

  lower(c, U, mu, nu);

  return c;
}

//====================================================================
void Staples::upper(Field_G& c, const Field_G& U,
                    const int mu, const int nu)
{
  // (1)  mu (2)
  //    +-->--+
  // nu |     |
  //   i+     +

  Field_G Umu, Unu;

  //#pragma omp parallel
  {
    //  m_shift->backward_ex(m_v, 0, U, nu, mu);
    //  m_shift->backward_ex(c,   0, U, mu, nu);
    Umu.setpart_ex(0, U, mu);
    Unu.setpart_ex(0, U, nu);
    m_shift->backward(m_v, Unu, mu);
    m_shift->backward(c, Umu, nu);

    mult_Field_Gnd(m_w, 0, c, 0, m_v, 0);
    mult_Field_Gnn(c, 0, U, nu, m_w, 0);
  }

}

//====================================================================
void Staples::lower(Field_G& c, const Field_G& U,
                    const int mu, const int nu)
{
  //    +     +
  // nu |     |
  //   i+-->--+
  //  (1)  mu (2)

  Field_G Umu, Unu;

  //#pragma omp parallel
  {
    //    m_shift->backward_ex(m_w, 0, U, nu, mu);
    Unu.setpart_ex(0, U, nu);
    m_shift->backward(m_w, Unu, mu);

    mult_Field_Gnn(m_v, 0, U, mu, m_w, 0);
    mult_Field_Gdn(m_w, 0, U, nu, m_v, 0);

    //    m_shift->forward_ex(c, 0, m_w, 0, nu);
    m_shift->forward(c, m_w, nu);
  }

}

//====================================================================
//============================================================END=====