gauge_field.h

Go to the documentation of this file.

#ifndef _GAUGE_QUDA_H
#define _GAUGE_QUDA_H
 
#include <quda_internal.h>
#include <quda.h>
#include <lattice_field.h>
 
#include <comm_key.h>
 
namespace quda {
 
  namespace gauge
  {
 
    inline bool isNative(QudaGaugeFieldOrder order, QudaPrecision precision, QudaReconstructType reconstruct)
    {
      if (precision == QUDA_DOUBLE_PRECISION) {
        if (order == QUDA_FLOAT2_GAUGE_ORDER) return true;
      } else if (precision == QUDA_SINGLE_PRECISION) {
        if (reconstruct == QUDA_RECONSTRUCT_NO || reconstruct == QUDA_RECONSTRUCT_10) {
          if (order == QUDA_FLOAT2_GAUGE_ORDER) return true;
        } else if (reconstruct == QUDA_RECONSTRUCT_12 || reconstruct == QUDA_RECONSTRUCT_13
                   || reconstruct == QUDA_RECONSTRUCT_8 || reconstruct == QUDA_RECONSTRUCT_9) {
          if (order == QUDA_FLOAT4_GAUGE_ORDER) return true;
        }
      } else if (precision == QUDA_HALF_PRECISION || precision == QUDA_QUARTER_PRECISION) {
        if (reconstruct == QUDA_RECONSTRUCT_NO || reconstruct == QUDA_RECONSTRUCT_10) {
          if (order == QUDA_FLOAT2_GAUGE_ORDER) return true;
        } else if (reconstruct == QUDA_RECONSTRUCT_12 || reconstruct == QUDA_RECONSTRUCT_13) {
          if (order == QUDA_FLOAT4_GAUGE_ORDER) return true;
        } else if (reconstruct == QUDA_RECONSTRUCT_8 || reconstruct == QUDA_RECONSTRUCT_9) {
#ifdef FLOAT8
          if (order == QUDA_FLOAT8_GAUGE_ORDER) return true;
#else
          if (order == QUDA_FLOAT4_GAUGE_ORDER) return true;
#endif
        }
      }
      return false;
    }
 
  } // namespace gauge
 
  struct GaugeFieldParam : public LatticeFieldParam {
 
    QudaFieldLocation location; // where are we storing the field (CUDA or GPU)?
    int nColor;
    int nFace;
 
    QudaReconstructType reconstruct;
    QudaGaugeFieldOrder order;
    QudaGaugeFixed fixed;
    QudaLinkType link_type;
    QudaTboundary t_boundary;
 
    double anisotropy;
    double tadpole;
    void *gauge; // used when we use a reference to an external field
 
    QudaFieldCreate create; // used to determine the type of field created
 
    QudaFieldGeometry geometry; // whether the field is a scale, vector or tensor
 
    // whether we need to compute the fat link maxima
    // FIXME temporary flag until we have a kernel that can do this, then we just do this in copy()
    // always set to false, requires external override
    bool compute_fat_link_max;
 
    QudaStaggeredPhase staggeredPhaseType;
 
    bool staggeredPhaseApplied;
 
    double i_mu;
 
    size_t site_offset;
 
    size_t site_size;
 
    // Default constructor
    GaugeFieldParam(void *const h_gauge = NULL) :
      LatticeFieldParam(),
      location(QUDA_INVALID_FIELD_LOCATION),
      nColor(3),
      nFace(0),
      reconstruct(QUDA_RECONSTRUCT_NO),
      order(QUDA_INVALID_GAUGE_ORDER),
      fixed(QUDA_GAUGE_FIXED_NO),
      link_type(QUDA_WILSON_LINKS),
      t_boundary(QUDA_INVALID_T_BOUNDARY),
      anisotropy(1.0),
      tadpole(1.0),
      gauge(h_gauge),
      create(QUDA_REFERENCE_FIELD_CREATE),
      geometry(QUDA_VECTOR_GEOMETRY),
      compute_fat_link_max(false),
      staggeredPhaseType(QUDA_STAGGERED_PHASE_NO),
      staggeredPhaseApplied(false),
      i_mu(0.0),
      site_offset(0),
      site_size(0)
    {
    }
 
    GaugeFieldParam(const GaugeField &u);
 
    GaugeFieldParam(const int *x, const QudaPrecision precision, const QudaReconstructType reconstruct, const int pad,
                    const QudaFieldGeometry geometry, const QudaGhostExchange ghostExchange = QUDA_GHOST_EXCHANGE_PAD) :
      LatticeFieldParam(4, x, pad, precision, ghostExchange),
      location(QUDA_INVALID_FIELD_LOCATION),
      nColor(3),
      nFace(0),
      reconstruct(reconstruct),
      order(QUDA_INVALID_GAUGE_ORDER),
      fixed(QUDA_GAUGE_FIXED_NO),
      link_type(QUDA_WILSON_LINKS),
      t_boundary(QUDA_INVALID_T_BOUNDARY),
      anisotropy(1.0),
      tadpole(1.0),
      gauge(0),
      create(QUDA_NULL_FIELD_CREATE),
      geometry(geometry),
      compute_fat_link_max(false),
      staggeredPhaseType(QUDA_STAGGERED_PHASE_NO),
      staggeredPhaseApplied(false),
      i_mu(0.0),
      site_offset(0),
      site_size(0)
    {
    }
 
    GaugeFieldParam(void *h_gauge, const QudaGaugeParam &param, QudaLinkType link_type_ = QUDA_INVALID_LINKS) :
      LatticeFieldParam(param),
      location(QUDA_CPU_FIELD_LOCATION),
      nColor(3),
      nFace(0),
      reconstruct(QUDA_RECONSTRUCT_NO),
      order(param.gauge_order),
      fixed(param.gauge_fix),
      link_type(link_type_ != QUDA_INVALID_LINKS ? link_type_ : param.type),
      t_boundary(param.t_boundary),
      anisotropy(param.anisotropy),
      tadpole(param.tadpole_coeff),
      gauge(h_gauge),
      create(QUDA_REFERENCE_FIELD_CREATE),
      geometry(QUDA_VECTOR_GEOMETRY),
      compute_fat_link_max(false),
      staggeredPhaseType(param.staggered_phase_type),
      staggeredPhaseApplied(param.staggered_phase_applied),
      i_mu(param.i_mu),
      site_offset(param.gauge_offset),
      site_size(param.site_size)
    {
      switch (link_type) {
      case QUDA_SU3_LINKS:
      case QUDA_GENERAL_LINKS:
      case QUDA_SMEARED_LINKS:
      case QUDA_MOMENTUM_LINKS: nFace = 1; break;
      case QUDA_THREE_LINKS: nFace = 3; break;
      default: errorQuda("Error: invalid link type(%d)\n", link_type);
      }
    }
 
    void setPrecision(QudaPrecision precision, bool force_native = false)
    {
      // is the current status in native field order?
      bool native = force_native ? true : gauge::isNative(order, this->precision, reconstruct);
      this->precision = precision;
      this->ghost_precision = precision;
 
      if (native) {
        if (precision == QUDA_DOUBLE_PRECISION || reconstruct == QUDA_RECONSTRUCT_NO
            || reconstruct == QUDA_RECONSTRUCT_10) {
          order = QUDA_FLOAT2_GAUGE_ORDER;
        } else if ((precision == QUDA_HALF_PRECISION || precision == QUDA_QUARTER_PRECISION)
                   && (reconstruct == QUDA_RECONSTRUCT_8 || reconstruct == QUDA_RECONSTRUCT_9)) {
#ifdef FLOAT8
          order = QUDA_FLOAT8_GAUGE_ORDER;
#else
          order = QUDA_FLOAT4_GAUGE_ORDER;
#endif
        } else {
          order = QUDA_FLOAT4_GAUGE_ORDER;
        }
      }
    }
  };
 
  std::ostream& operator<<(std::ostream& output, const GaugeFieldParam& param);
 
  class GaugeField : public LatticeField {
 
  protected:
      size_t bytes;        // bytes allocated per full field
      size_t phase_offset; // offset in bytes to gauge phases - useful to keep track of texture alignment
      size_t phase_bytes;  // bytes needed to store the phases
      size_t length;
      size_t real_length;
      int nColor;
      int nFace;
      QudaFieldGeometry geometry; // whether the field is a scale, vector or tensor
 
      QudaReconstructType reconstruct;
      int nInternal; // number of degrees of freedom per link matrix
      QudaGaugeFieldOrder order;
      QudaGaugeFixed fixed;
      QudaLinkType link_type;
      QudaTboundary t_boundary;
 
      double anisotropy;
      double tadpole;
      double fat_link_max;
 
      QudaFieldCreate create; // used to determine the type of field created
 
      mutable void *ghost[2 * QUDA_MAX_DIM]; // stores the ghost zone of the gauge field (non-native fields only)
 
      mutable int ghostFace[QUDA_MAX_DIM]; // the size of each face
 
      QudaStaggeredPhase staggeredPhaseType;
 
      bool staggeredPhaseApplied;
 
      void exchange(void **recv, void **send, QudaDirection dir) const;
 
      double i_mu;
 
      size_t site_offset;
 
      size_t site_size;
 
      void createGhostZone(const int *R, bool no_comms_fill, bool bidir = true) const;
 
      void setTuningString();
 
  public:
    GaugeField(const GaugeFieldParam &param);
    virtual ~GaugeField();
 
    virtual void exchangeGhost(QudaLinkDirection = QUDA_LINK_BACKWARDS) = 0;
    virtual void injectGhost(QudaLinkDirection = QUDA_LINK_BACKWARDS) = 0;
 
    size_t Length() const { return length; }
    int Ncolor() const { return nColor; }
    QudaReconstructType Reconstruct() const { return reconstruct; }
    QudaGaugeFieldOrder Order() const { return order; }
    double Anisotropy() const { return anisotropy; }
    double Tadpole() const { return tadpole; }
    QudaTboundary TBoundary() const { return t_boundary; }
    QudaLinkType LinkType() const { return link_type; }
    QudaGaugeFixed GaugeFixed() const { return fixed; }
    QudaGaugeFieldOrder FieldOrder() const { return order; }
    QudaFieldGeometry Geometry() const { return geometry; }
    QudaStaggeredPhase StaggeredPhase() const { return staggeredPhaseType; }
    bool StaggeredPhaseApplied() const { return staggeredPhaseApplied; }
 
    using param_type = GaugeFieldParam;
 
    void applyStaggeredPhase(QudaStaggeredPhase phase=QUDA_STAGGERED_PHASE_INVALID);
 
    void removeStaggeredPhase();
 
    double iMu() const { return i_mu; }
 
    const double& LinkMax() const { return fat_link_max; }
    int Nface() const { return nFace; }
 
    virtual void exchangeExtendedGhost(const int *R, bool no_comms_fill = false) = 0;
 
    virtual void exchangeExtendedGhost(const int *R, TimeProfile &profile, bool no_comms_fill = false) = 0;
 
    void checkField(const LatticeField &) const;
 
    bool isNative() const { return gauge::isNative(order, precision, reconstruct); }
 
    size_t Bytes() const { return bytes; }
    size_t PhaseBytes() const { return phase_bytes; }
    size_t PhaseOffset() const { return phase_offset; }
 
    size_t TotalBytes() const { return bytes; }
 
    virtual void* Gauge_p() { errorQuda("Not implemented"); return (void*)0;}
    virtual void* Even_p() { errorQuda("Not implemented"); return (void*)0;}
    virtual void* Odd_p() { errorQuda("Not implemented"); return (void*)0;}
 
    virtual const void* Gauge_p() const { errorQuda("Not implemented"); return (void*)0;}
    virtual const void* Even_p() const { errorQuda("Not implemented"); return (void*)0;}
    virtual const void* Odd_p() const { errorQuda("Not implemented"); return (void*)0;}
 
    virtual int full_dim(int d) const { return x[d]; }
 
    const void** Ghost() const {
      if ( isNative() ) errorQuda("No ghost zone pointer for quda-native gauge fields");
      return (const void**)ghost;
    }
 
    void** Ghost() {
      if ( isNative() ) errorQuda("No ghost zone pointer for quda-native gauge fields");
      return ghost;
    }
 
    size_t SiteOffset() const { return site_offset; }
 
    size_t SiteSize() const { return site_size; }
 
    virtual void zero() = 0;
 
    virtual void copy(const GaugeField &src) = 0;
 
    double norm1(int dim = -1, bool fixed = false) const;
 
    double norm2(int dim = -1, bool fixed = false) const;
 
    double abs_max(int dim = -1, bool fixed = false) const;
 
    double abs_min(int dim = -1, bool fixed = false) const;
 
    uint64_t checksum(bool mini = false) const;
 
    static GaugeField* Create(const GaugeFieldParam &param);
 
  };
 
  class cudaGaugeField : public GaugeField {
 
  private:
    void *gauge;
    void *gauge_h; // mapped-memory pointer when allocating on the host
    void *even;
    void *odd;
 
    void zeroPad();
 
  public:
    cudaGaugeField(const GaugeFieldParam &);
    virtual ~cudaGaugeField();
 
    void exchangeGhost(QudaLinkDirection link_direction = QUDA_LINK_BACKWARDS);
 
    void injectGhost(QudaLinkDirection link_direction = QUDA_LINK_BACKWARDS);
 
    void createComms(const int *R, bool no_comms_fill, bool bidir=true);
 
    void allocateGhostBuffer(const int *R, bool no_comms_fill, bool bidir=true) const;
 
    void recvStart(int dim, int dir);
 
    void sendStart(int dim, int dir, qudaStream_t *stream_p = nullptr);
 
    void commsComplete(int dim, int dir);
 
    void exchangeExtendedGhost(const int *R, bool no_comms_fill=false);
 
    void exchangeExtendedGhost(const int *R, TimeProfile &profile, bool no_comms_fill=false);
 
    void copy(const GaugeField &src);
 
    void loadCPUField(const cpuGaugeField &cpu);
 
    void loadCPUField(const cpuGaugeField &cpu, TimeProfile &profile);
 
    void saveCPUField(cpuGaugeField &cpu) const;
 
    void saveCPUField(cpuGaugeField &cpu, TimeProfile &profile) const;
 
    // (ab)use with care
    void* Gauge_p() { return gauge; }
    void* Even_p() { return even; }
    void* Odd_p() { return odd; }
 
    const void* Gauge_p() const { return gauge; }
    const void* Even_p() const { return even; }
    const void *Odd_p() const { return odd; }
 
    virtual void copy_to_buffer(void *buffer) const;
 
    virtual void copy_from_buffer(void *buffer);
 
    void setGauge(void* _gauge); //only allowed when create== QUDA_REFERENCE_FIELD_CREATE
 
    void zero();
 
    void backup() const;
 
    void restore() const;
 
    void prefetch(QudaFieldLocation mem_space, qudaStream_t stream = 0) const;
  };
 
  class cpuGaugeField : public GaugeField {
 
    friend void cudaGaugeField::copy(const GaugeField &cpu);
    friend void cudaGaugeField::loadCPUField(const cpuGaugeField &cpu);
    friend void cudaGaugeField::saveCPUField(cpuGaugeField &cpu) const;
 
  private:
    void **gauge; // the actual gauge field
 
  public:
    cpuGaugeField(const GaugeFieldParam &param);
    virtual ~cpuGaugeField();
 
    void exchangeGhost(QudaLinkDirection link_direction = QUDA_LINK_BACKWARDS);
 
    void injectGhost(QudaLinkDirection link_direction = QUDA_LINK_BACKWARDS);
 
    void exchangeExtendedGhost(const int *R, bool no_comms_fill=false);
 
    void exchangeExtendedGhost(const int *R, TimeProfile &profile, bool no_comms_fill=false);
 
    void copy(const GaugeField &src);
 
    void* Gauge_p() { return gauge; }
    const void* Gauge_p() const { return gauge; }
 
    virtual void copy_to_buffer(void *buffer) const;
 
    virtual void copy_from_buffer(void *buffer);
 
    void setGauge(void** _gauge); //only allowed when create== QUDA_REFERENCE_FIELD_CREATE
 
    void zero();
 
    void backup() const;
 
    void restore() const;
  };
 
  double norm1(const GaugeField &u);
 
  double norm2(const GaugeField &u);
 
  void ax(const double &a, GaugeField &u);
 
  void copyGenericGauge(GaugeField &out, const GaugeField &in, QudaFieldLocation location, void *Out = 0, void *In = 0,
                        void **ghostOut = 0, void **ghostIn = 0, int type = 0);
 
  void copyFieldOffset(GaugeField &out, const GaugeField &in, CommKey offset, QudaPCType pc_type);
 
  void copyExtendedGauge(GaugeField &out, const GaugeField &in,
                         QudaFieldLocation location, void *Out=0, void *In=0);
 
  cudaGaugeField *createExtendedGauge(cudaGaugeField &in, const int *R, TimeProfile &profile,
                                      bool redundant_comms = false, QudaReconstructType recon = QUDA_RECONSTRUCT_INVALID);
 
  cpuGaugeField *createExtendedGauge(void **gauge, QudaGaugeParam &gauge_param, const int *R);
 
  void extractGaugeGhost(const GaugeField &u, void **ghost, bool extract=true, int offset=0);
 
  void extractExtendedGaugeGhost(const GaugeField &u, int dim, const int *R, void **ghost, bool extract);
 
  void applyGaugePhase(GaugeField &u);
 
  uint64_t Checksum(const GaugeField &u, bool mini=false);
 
  inline QudaReconstructType Reconstruct_(const char *func, const char *file, int line, const GaugeField &a,
                                          const GaugeField &b)
  {
    QudaReconstructType reconstruct = QUDA_RECONSTRUCT_INVALID;
    if (a.Reconstruct() == b.Reconstruct())
      reconstruct = a.Reconstruct();
    else
      errorQuda("Reconstruct %d %d do not match (%s:%d in %s())\n", a.Reconstruct(), b.Reconstruct(), file, line, func);
    return reconstruct;
  }
 
  template <typename... Args>
  inline QudaReconstructType Reconstruct_(const char *func, const char *file, int line, const GaugeField &a,
                                          const GaugeField &b, const Args &... args)
  {
    return static_cast<QudaReconstructType>(Reconstruct_(func, file, line, a, b)
                                            & Reconstruct_(func, file, line, a, args...));
  }
 
#define checkReconstruct(...) Reconstruct_(__func__, __FILE__, __LINE__, __VA_ARGS__)
 
} // namespace quda
 
#endif // _GAUGE_QUDA_H