Add HYP smearing. Also add 3D/4D version of STOUT and APE smearing.

include/enum_quda.h

@@ -571,6 +571,7 @@ typedef enum QudaGaugeSmearType_s {
   QUDA_GAUGE_SMEAR_APE,
   QUDA_GAUGE_SMEAR_STOUT,
   QUDA_GAUGE_SMEAR_OVRIMP_STOUT,
+  QUDA_GAUGE_SMEAR_HYP,
   QUDA_GAUGE_SMEAR_WILSON_FLOW,
   QUDA_GAUGE_SMEAR_SYMANZIK_FLOW,
   QUDA_GAUGE_SMEAR_INVALID = QUDA_INVALID_ENUM

include/enum_quda_fortran.h

@@ -504,8 +504,9 @@
 #define QUDA_GAUGE_SMEAR_APE 0
 #define QUDA_GAUGE_SMEAR_STOUT 1
 #define QUDA_GAUGE_SMEAR_OVRIMP_STOUT 2
-#define QUDA_GAUGE_SMEAR_WILSON_FLOW 3
-#define QUDA_GAUGE_SMEAR_SYMANZIK_FLOW 4
+#define QUDA_GAUGE_SMEAR_HYP 3
+#define QUDA_GAUGE_SMEAR_WILSON_FLOW 4
+#define QUDA_GAUGE_SMEAR_SYMANZIK_FLOW 5
 #define QUDA_GAUGE_SMEAR_INVALID QUDA_INVALID_ENUM
 
 #define QudaFermionSmearType integer(4)

include/gauge_field.h

@@ -638,6 +638,19 @@ namespace quda {
      */
     static bool are_compatible_weak(const GaugeField &a, const GaugeField &b);
 
+    /**
+     * @brief Helper function that returns the default tolerance used by SU(3) projection
+     * @return The default tolerance, which is ~10x epsilon
+     */
+    double toleranceSU3() const
+    {
+      switch (precision) {
+      case QUDA_DOUBLE_PRECISION: return 2e-15;
+      case QUDA_SINGLE_PRECISION: return 1e-6;
+      default: return 1e-6;
+      }
+    }
+
     friend struct GaugeFieldParam;
   };
 

include/gauge_tools.h

@@ -95,17 +95,19 @@ namespace quda
      @param[out] dataDs Output smeared field
      @param[in] dataOr Input gauge field
      @param[in] alpha smearing parameter
+     @param[in] dir_ignore ignored direction
   */
-  void APEStep(GaugeField &dataDs, GaugeField &dataOr, double alpha);
+  void APEStep(GaugeField &dataDs, GaugeField &dataOr, double alpha, int dir_ignore);
 
   /**
      @brief Apply STOUT smearing to the gauge field
 
      @param[out] dataDs Output smeared field
      @param[in] dataOr Input gauge field
      @param[in] rho smearing parameter
+     @param[in] dir_ignore ignored direction
   */
-  void STOUTStep(GaugeField &dataDs, GaugeField &dataOr, double rho);
+  void STOUTStep(GaugeField &dataDs, GaugeField &dataOr, double rho, int dir_ignore);
 
   /**
      @brief Apply Over Improved STOUT smearing to the gauge field
@@ -114,8 +116,20 @@ namespace quda
      @param[in] dataOr Input gauge field
      @param[in] rho smearing parameter
      @param[in] epsilon smearing parameter
+     @param[in] dir_ignore ignored direction
   */
-  void OvrImpSTOUTStep(GaugeField &dataDs, GaugeField &dataOr, double rho, double epsilon);
+  void OvrImpSTOUTStep(GaugeField &dataDs, GaugeField &dataOr, double rho, double epsilon, int dir_ignore);
+
+  /**
+     @brief Apply HYP smearing to the gauge field
+     @param[out] dataDs Output smeared field
+     @param[in] dataOr Input gauge field
+     @param[in] alpha1 smearing parameter
+     @param[in] alpha2 smearing parameter
+     @param[in] alpha3 smearing parameter
+     @param[in] dir_ignore ignored direction
+  */
+  void HYPStep(GaugeField &dataDs, GaugeField &dataOr, double alpha1, double alpha2, double alpha3, int dir_ignore);
 
   /**
      @brief Apply Wilson Flow steps W1, W2, Vt to the gauge field.

include/kernels/gauge_ape.cuh

@@ -8,9 +8,6 @@
 namespace quda
 {
 
-#define  DOUBLE_TOL	1e-15
-#define  SINGLE_TOL	2e-6
-
   template <typename Float_, int nColor_, QudaReconstructType recon_, int apeDim_>
   struct GaugeAPEArg : kernel_param<> {
     using Float = Float_;
@@ -26,14 +23,16 @@ namespace quda
     int X[4];    // grid dimensions
     int border[4];
     const Float alpha;
+    const int dir_ignore;
     const Float tolerance;
 
-    GaugeAPEArg(GaugeField &out, const GaugeField &in, double alpha) :
+    GaugeAPEArg(GaugeField &out, const GaugeField &in, double alpha, int dir_ignore) :
       kernel_param(dim3(in.LocalVolumeCB(), 2, apeDim)),
       out(out),
       in(in),
       alpha(alpha),
-      tolerance(in.Precision() == QUDA_DOUBLE_PRECISION ? DOUBLE_TOL : SINGLE_TOL)
+      dir_ignore(dir_ignore),
+      tolerance(in.toleranceSU3())
     {
       for (int dir = 0; dir < 4; ++dir) {
         border[dir] = in.R()[dir];
@@ -61,16 +60,17 @@ namespace quda
         x[dr] += arg.border[dr];
         X[dr] += 2 * arg.border[dr];
       }
+      dir = dir + (dir >= arg.dir_ignore);
 
       int dx[4] = {0, 0, 0, 0};
       Link U, Stap, TestU, I;
       // This function gets stap = S_{mu,nu} i.e., the staple of length 3,
-      computeStaple(arg, x, X, parity, dir, Stap, Arg::apeDim);
+      computeStaple(arg, x, X, parity, dir, Stap, arg.dir_ignore);
 
       // Get link U
       U = arg.in(dir, linkIndexShift(x, dx, X), parity);
-    
-      Stap = Stap * (arg.alpha / ((real)(2. * (3. - 1.))));
+
+      Stap = Stap * (arg.alpha / ((real)(2 * (Arg::apeDim - 1))));
       setIdentity(&I);
 
       TestU = I * (static_cast<real>(1.0) - arg.alpha) + Stap * conj(U);