field.h

Go to the documentation of this file.

#ifndef FIELD_INCLUDED
#define FIELD_INCLUDED

#include <iostream>
#include <valarray>
#include <string>
#include <assert.h>

#include "Parameters/commonParameters.h"
#include "Communicator/communicator.h"

#include "bridge_complex.h"
#include "bridge_defs.h"


class Field
{
 public:
  //  enum element_type { REAL = 1, COMPLEX = 2 };
  typedef Element_type::type   element_type;
  typedef double               real_t; // [23 Apr 2018 H.M. added]

 protected:
  int m_Nin;           // internal d.o.f.
  int m_Nvol;          // lattice volume
  int m_Nex;           // external d.o.f.
  // Element_type::type  m_element_type;
  element_type m_element_type;
  // total degree of freedom is Nin * Nsite * Nex.

  std::valarray<double> field;

  // as valarray takes size_t, we do not use long_t here
  inline
  size_t myindex(const int jin, const int site, const int jex) const
  {
    return jin + m_Nin * ((size_t)site + m_Nvol * jex);
  }

  Bridge::VerboseLevel m_vl;

 public:

  Field() :
    m_Nin(0), m_Nvol(0), m_Nex(0), m_element_type(Element_type::COMPLEX),
    field(0),
    m_vl(CommonParameters::Vlevel())
  {
    check();
  }

  Field(const int Nin, const int Nvol, const int Nex,
        const element_type cmpl = Element_type::COMPLEX) :
    m_Nin(Nin), m_Nvol(Nvol), m_Nex(Nex), m_element_type(cmpl),
    field(ntot()),
    m_vl(CommonParameters::Vlevel())
  {
    check();
  }

  Field clone() const
  {
    return Field(m_Nin, m_Nvol, m_Nex, m_element_type);
  }

  void reset(const int Nin, const int Nvol, const int Nex,
             const element_type cmpl = Element_type::COMPLEX)
  {
    if ((m_Nin == Nin) &&
        (m_Nvol == Nvol) &&
        (m_Nex == Nex) &&
        (m_element_type == cmpl)) return;

    m_Nin          = Nin;
    m_Nvol         = Nvol;
    m_Nex          = Nex;
    m_element_type = cmpl;

    field.resize(ntot());
  }

  // assignment
  //Field& operator=(const double a) { field = a; return *this; }
  Field& operator=(const Field& v)
  {
    assert(m_Nin == v.nin());
    assert(m_Nvol == v.nvol());
    assert(m_Nex == v.nex());
    copy(*this, v);
    return *this;
  }

 private:
  void check();

 public:
  int nin() const { return m_Nin; }
  int nvol() const { return m_Nvol; }
  int nex() const { return m_Nex; }
  element_type field_element_type() const { return m_element_type; }

  int ntot() const { return m_Nin * m_Nvol * m_Nex; }
  int size() const { return ntot(); }

  bool check_size(const int nin, const int nvol, const int nex) const
  {
    bool chk = true;

    if ((m_Nin != nin) || (m_Nvol != nvol) || (m_Nex != nex)) chk = false;
    return chk;
  }

  double cmp(const int jin, const int site, const int jex) const
  {
    return field[myindex(jin, site, jex)];
  }

  double cmp(const int i) const
  {
    return field[i];
  }

  const double *ptr(const int jin, const int site, const int jex) const
  {
    // when using c++03, const_cast is necessary.
    return &(const_cast<std::valarray<double>&>(field)[myindex(jin, site, jex)]);
  }

  double *ptr(const int jin, const int site, const int jex)
  {
    return &field[myindex(jin, site, jex)];
  }

  const double *ptr(const int i) const
  {
    // when using c++03, const_cast is necessary.
    return &(const_cast<std::valarray<double>&>(field)[i]);
  }

  double *ptr(const int i)
  {
    return &field[i];
  }

  void set(const int jin, const int site, const int jex, double v)
  {
    field[myindex(jin, site, jex)] = v;
  }

  void set(const int i, double v)
  {
    field[i] = v;
  }

  void set(double a);

  void setc(double a);

  void setc(dcomplex a);

  void add(const int jin, const int site, const int jex, double v)
  {
    field[myindex(jin, site, jex)] += v;
  }

  void add(const int i, double v)
  {
    field[i] += v;
  }

  void setpart_ex(int ex, const Field& w, int exw)
  {
    assert(ex < nex());
    assert(exw < w.nex());
    return copy(*this, ex, w, exw);
  }

  void addpart_ex(int ex, const Field& w, int exw)
  {
    assert(ex < nex());
    assert(exw < w.nex());
    return axpy(*this, ex, 1.0, w, exw);
  }

  void addpart_ex(int ex, const Field& w, int exw, double prf)
  {
    assert(ex < nex());
    assert(exw < w.nex());
    return axpy(*this, ex, prf, w, exw);
  }

  // BLAS-like routine (friend declaration)

  double norm2() const;

  double norm() const { return sqrt(norm2()); }

  friend
  double dot(const Field& y, const Field& x);

  friend
  double dot(const Field& y, const int exy, const Field& x, const int exx);

  friend
  void dot_and_norm2(double& yx, double& y2, double& x2, const Field& y, const Field& x);

  friend
  void dot_and_norm2(double& yx, double& y2, double& x2, const Field& y, const int exy, const Field& x, const int exx);

  friend
  dcomplex dotc(const Field& y, const Field& x);

  friend
  dcomplex dotc(const Field& y, const int exy, const Field& x, const int exx);

  friend
  void dotc_and_norm2(double& yx, double& y2, double& x2, const Field& y, const Field& x);

  friend
  void dotc_and_norm2(double& yx, double& y2, double& x2, const Field& y, const int exy, const Field& x, const int exx);

  friend
  void axpy(Field& y, const double a, const Field& x);

  friend
  void axpy(Field& y, const int exy, const double a, const Field& x, const int exx);

  friend
  void axpy(Field& y, const dcomplex a, const Field& x);

  friend
  void axpy(Field& y, const int exy, const dcomplex a, const Field& x, const int exx);

  friend
  void scal(Field& x, const double a);

  friend
  void scal(Field& x, const int exx, const double a);

  friend
  void scal(Field& x, const dcomplex a);

  friend
  void scal(Field& x, const int exx, const dcomplex a);

  friend
  void copy(Field& y, const Field& x);

  friend
  void copy(Field& y, const int exy, const Field& x, const int exx);

  friend
  void aypx(const double a, Field& y, const Field& x);

  friend
  void aypx(const dcomplex a, Field& y, const Field& x);


  void stat(double& Fave, double& Fmax, double& Fdev) const;
};

//----------------------------------------------------------------
// BLAS-like routines defined outside the Field class.

double dot(const Field& y, const Field& x);

double dot(const Field& y, const int exy, const Field& x, const int exx);

dcomplex dotc(const Field& y, const Field& x);

dcomplex dotc(const Field& y, const int exy, const Field& x, const int nexx);

void axpy(Field& y, const double a, const Field& x);

void axpy(Field& y, const int exy, const double a, const Field& x, const int exx);

void axpy(Field& y, const dcomplex a, const Field& x);

void axpy(Field& y, const int exy, const dcomplex a, const Field& x, const int exx);

void scal(Field& x, const double a);

void scal(Field& x, const int exx, const double a);

void scal(Field& x, const dcomplex a);

void scal(Field& x, const int exx, const dcomplex a);

void copy(Field& y, const Field& x);

void copy(Field& y, const int exy, const Field& x, const int exx);

void aypx(const double a, Field& y, const Field& x);

void aypx(const dcomplex a, Field& y, const Field& x);

//----------------------------------------------------------------

inline int exchange(int count, Field *recv_buf, Field *send_buf, int idir, int ipm, int tag)
{
  return Communicator::exchange(count, recv_buf->ptr(0), send_buf->ptr(0), idir, ipm, tag);
}


inline int send_1to1(int count, Field *recv_buf, Field *send_buf, int p_to, int p_from, int tag)
{
  return Communicator::send_1to1(count, recv_buf->ptr(0), send_buf->ptr(0), p_to, p_from, tag);
}


//----------------------------------------------------------------
// utility: report statistics.
void report_field_stat(const Bridge::VerboseLevel vl, const std::string& msg, const Field& f);

//----------------------------------------------------------------
#endif