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 "commonParameters.h"
#include "communicator.h"

#include "bridge_complex.h"


class Field
{
 public:
  enum element_type { REAL = 1, COMPLEX = 2 };

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

  std::valarray<double> field;

  inline
  int myindex(const int jin, const int site, const int jex) const
  {
    return jin + m_Nin * (site + m_Nvol * jex);
  }

  Bridge::VerboseLevel m_vl;

 public:

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

  Field(const int Nin, const int Nvol, const int Nex, const element_type cmpl = 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 = 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(); }

  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 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());
//     for (int site = 0; site < m_Nvol; ++site) {
//       for (int jin = 0; jin < m_Nin; ++jin) {
//         field[myindex(jin, site, ex)] = w.field[myindex(jin, site, exw)];
//       }
//     }
    return copy(*this, ex, w, exw);
  }

  void addpart_ex(int ex, const Field& w, int exw)
  {
    assert(ex < nex());
    assert(exw < w.nex());
//     for (int site = 0; site < m_Nvol; ++site) {
//       for (int jin = 0; jin < m_Nin; ++jin) {
//         field[myindex(jin, site, ex)] += w.field[myindex(jin, site, exw)];
//       }
//     }
    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());
//     for (int site = 0; site < m_Nvol; ++site) {
//       for (int jin = 0; jin < m_Nin; ++jin) {
//         field[myindex(jin, site, ex)]
//           += prf * w.field[myindex(jin, site, exw)];
//       }
//     }
    return axpy(*this, ex, prf, w, exw);
  }

  // linear algebra operations
  //  Field& operator-() { field *= -1; return *this; }

  /*
  Field& operator+=(const Field& v) { field += v.field; return *this; }
  Field& operator-=(const Field& v) { field -= v.field; return *this; }
  Field& operator*=(const double a) { field *= a; return *this; }
  Field& operator/=(const double a) { field *= 1.0 / a; return *this; }
  Field& operator*=(const dcomplex a)
  {
    if (imag(a) == 0.0) {
      field *= real(a);
    } else {
      scal(*this, a);
    }
    return *this;
  }
  */

  // XXX
  //double operator*(const Field& rhs) { return dot(*this, rhs); }

  // BLAS-like routine

  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
  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 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;
};

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

/*
inline Field operator*(const Field& v, const double s)
{
  Field w(v);

  w *= s;
  return w;
}


inline Field operator*(const double s, const Field& v)
{
  return operator*(v, s);
}


inline Field operator+(const Field& lhs, const Field& rhs)
{
  Field w(lhs);

  w += rhs;
  return w;
}


inline Field operator-(const Field& lhs, const Field& rhs)
{
  Field w(lhs);

  w -= rhs;
  return w;
}
*/

//----------------------------------------------------------------
// BLAS-like routines

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