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]
  static const std::string class_name;
 
 protected:
  int m_Nin;                   
  int m_Nvol;                  
  int m_Nex;                   
  element_type m_element_type; 
 
  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 Field& v)
  {
    assert(check_size(v.nin(), v.nvol(), v.nex()));
    // 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()); }
 
  void xI();
 
  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);
 
//----------------------------------------------------------------
#endif