staggered contraction code

.github/workflows/cuda_githubactions_build.yml

@@ -48,6 +48,7 @@ jobs:
           -DCMAKE_CUDA_COMPILER_LAUNCHER=ccache
           -DCMAKE_CXX_COMPILER_LAUNCHER=ccache
           -DQUDA_GPU_ARCH=sm_80 -DQUDA_GPU_ARCH_SUFFIX=virtual -DQUDA_JITIFY=ON
+          -DQUDA_COVDEV=ON
           -DQUDA_MULTIGRID=ON
           -DQUDA_MULTIGRID_NVEC_LIST=24
           -DQUDA_MDW_FUSED_LS_LIST=4

.github/workflows/rocm-build-ci.yml

@@ -30,6 +30,7 @@ jobs:
             -DQUDA_DIRAC_WILSON=ON \
             -DQUDA_DIRAC_LAPLACE=ON \
             -DQUDA_CLOVER_DYNAMIC=ON \
+            -DQUDA_COVDEV=ON \
             -DQUDA_QDPJIT=OFF \
             -DQUDA_INTERFACE_QDPJIT=OFF \
             -DQUDA_INTERFACE_MILC=ON \

.gitignore

@@ -4,7 +4,6 @@
 *.a
 *~
 tests/*_test
-make.inc
 milc_interface/*
 *#*
 *.pyc

CMakeLists.txt

@@ -154,7 +154,6 @@ option(QUDA_DIRAC_LAPLACE "build laplace operator" ${QUDA_DIRAC_DEFAULT})
 
 option(QUDA_DIRAC_DISTANCE_PRECONDITIONING "build code for distance preconditioned Wilson/clover Dirac operators" OFF)
 
-option(QUDA_CONTRACT "build code for bilinear contraction" OFF)
 option(QUDA_COVDEV "build code for covariant derivative" OFF)
 
 option(QUDA_QIO "build QIO code for binary I/O" OFF)

README.md

@@ -283,6 +283,7 @@ Advanced Scientific Computing (PASC21) [arXiv:2104.05615[hep-lat]].
 *  Mario Schröck (INFN)
 *  Aniket Sen (HISKP, University of Bonn)
 *  Guochun Shi (NCSA)
+*  James Simone (Fermi National Accelerator Laboratory)
 *  Alexei Strelchenko (Fermi National Accelerator Laboratory)
 *  Jiqun Tu (NVIDIA)
 *  Carsten Urbach (HISKP, University of Bonn)

ci/docker/Dockerfile.build

@@ -45,6 +45,7 @@ RUN  QUDA_TEST_GRID_SIZE="1 1 1 2" cmake -S /quda/src \
     -DQUDA_DIRAC_TWISTED_CLOVER=ON \
     -DQUDA_DIRAC_STAGGERED=ON \
     -DQUDA_DIRAC_LAPLACE=ON \
+    -DQUDA_COVDEV=ON \
     -GNinja \
     -B /quda/build 
 

include/contract_quda.h

@@ -4,6 +4,32 @@
 
 namespace quda
 {
+  /**
+   * Interface function that launch contraction compute kernels,
+   * used in interface_quda.cpp
+   * @param[in] x               input source field
+   * @param[in] y               input source field
+   * @param[out] result         container of complex contraction results for
+   *                            all decay slices and spins
+   * @param[in] cType           contraction types as defined in QudaContractType enum
+   * @param[in] source_position 4d array of source position
+   * @param[in] mom_mode        4d array of momentum
+   * @param[in] fft_type        Fourier phase factor type
+   *                            as defined in QudaFFTSymmType enum
+   * @param[in] s1              spin component index (0 for staggered)
+   * @param[in] b1              spin component index (0 for staggered)
+   */
+
+  void contractSummedQuda(const ColorSpinorField &x, const ColorSpinorField &y, std::vector<Complex> &result,
+                          QudaContractType cType, const int *const source_position, const int *const mom_mode,
+                          const QudaFFTSymmType *const fft_type, const size_t s1, const size_t b1);
+  /**
+   * @param[in] x       input color spinor
+   * @param[in] y       input color spinor
+   * @param[out] result pointer to the spinxspin projections per lattice site
+   * @param[in] cType   contraction type
+   */
+
   void contractQuda(const ColorSpinorField &x, const ColorSpinorField &y, void *result, QudaContractType cType);
 
   /**

include/dslash_quda.h

@@ -927,4 +927,26 @@ namespace quda
   */
   void gamma5(ColorSpinorField &out, const ColorSpinorField &in);
 
+  /**
+     @brief Applies a (1 \pm gamma5)/2 projection matrix to a spinor
+     @param[out] out Output field
+     @param[in] in Input field
+     @param[in] proj Sign of \pm projection
+  */
+  void ApplyChiralProj(ColorSpinorField &out, const ColorSpinorField &in, const int proj);
+
+  /**
+     @brief Constructs the mid-point 4D propagator from a 5D domain wall propagator
+     @param[out] out Output field
+     @param[in] in Input field
+  */
+  void make4DMidPointProp(ColorSpinorField &out, ColorSpinorField &in);
+
+  /**
+     @brief Constructs the chiral 4D propagator from a 5D domain wall propagator
+     @param[out] out Output field
+     @param[in] in Input field
+  */
+  void make4DChiralProp(ColorSpinorField &out, ColorSpinorField &in);
+
 } // namespace quda

include/enum_quda.h

@@ -29,6 +29,7 @@ typedef enum QudaLinkType_s {
   QUDA_MOMENTUM_LINKS,
   QUDA_COARSE_LINKS,                  // used for coarse-gauge field with multigrid
   QUDA_SMEARED_LINKS,                 // used for loading and saving gaugeSmeared in the interface
+  QUDA_TWOLINK_LINKS,                 // used for staggered fermion smearing.
   QUDA_WILSON_LINKS = QUDA_SU3_LINKS, // used by wilson, clover, twisted mass, and domain wall
   QUDA_ASQTAD_FAT_LINKS = QUDA_GENERAL_LINKS,
   QUDA_ASQTAD_LONG_LINKS = QUDA_THREE_LINKS,
@@ -544,12 +545,47 @@ typedef enum QudaStaggeredPhase_s {
   QUDA_STAGGERED_PHASE_INVALID = QUDA_INVALID_ENUM
 } QudaStaggeredPhase;
 
+typedef enum QudaSpinTasteGamma_s {
+  QUDA_SPIN_TASTE_G1 = 0,
+  QUDA_SPIN_TASTE_GX = 1,
+  QUDA_SPIN_TASTE_GY = 2,
+  QUDA_SPIN_TASTE_GZ = 4,
+  QUDA_SPIN_TASTE_GT = 8,
+  QUDA_SPIN_TASTE_G5 = 15,
+  QUDA_SPIN_TASTE_GYGZ = 6,
+  QUDA_SPIN_TASTE_GZGX = 5,
+  QUDA_SPIN_TASTE_GXGY = 3,
+  QUDA_SPIN_TASTE_GXGT = 9,
+  QUDA_SPIN_TASTE_GYGT = 10,
+  QUDA_SPIN_TASTE_GZGT = 12,
+  QUDA_SPIN_TASTE_G5GX = 14,
+  QUDA_SPIN_TASTE_G5GY = 13,
+  QUDA_SPIN_TASTE_G5GZ = 11,
+  QUDA_SPIN_TASTE_G5GT = 7,
+  QUDA_SPIN_TASTE_INVALID = QUDA_INVALID_ENUM
+} QudaSpinTasteGamma;
+
 typedef enum QudaContractType_s {
-  QUDA_CONTRACT_TYPE_OPEN, // Open spin elementals
-  QUDA_CONTRACT_TYPE_DR,   // DegrandRossi
+  QUDA_CONTRACT_TYPE_STAGGERED_FT_T, // Staggered, FT in tdim
+  QUDA_CONTRACT_TYPE_DR_FT_T,        // DegrandRossi insertion, FT in tdim
+  QUDA_CONTRACT_TYPE_DR_FT_Z,        // DegrandRossi insertion, FT in zdim
+  QUDA_CONTRACT_TYPE_STAGGERED,      // Staggered, no summation (TODO: remove line)
+  QUDA_CONTRACT_TYPE_DR,             // DegrandRossi insertion, no summation
+  QUDA_CONTRACT_TYPE_OPEN,           // Open spin elementals, no summation
+  QUDA_CONTRACT_TYPE_OPEN_SUM_T,     // Open spin elementals, spatially summed over tdim
+  QUDA_CONTRACT_TYPE_OPEN_SUM_Z,     // Open spin elementals, spatially summed over zdim
+  QUDA_CONTRACT_TYPE_OPEN_FT_T,      // Open spin elementals, FT in tdim
+  QUDA_CONTRACT_TYPE_OPEN_FT_Z,      // Open spin elementals, FT in zdim
   QUDA_CONTRACT_TYPE_INVALID = QUDA_INVALID_ENUM
 } QudaContractType;
 
+typedef enum QudaFFTSymmType_t {
+  QUDA_FFT_SYMM_ODD = 1,  // sin(phase)
+  QUDA_FFT_SYMM_EVEN = 2, // cos(phase)
+  QUDA_FFT_SYMM_EO = 3,   // exp(-i phase)
+  QUDA_FFT_SYMM_INVALID = QUDA_INVALID_ENUM
+} QudaFFTSymmType;
+
 typedef enum QudaContractGamma_s {
   QUDA_CONTRACT_GAMMA_I = 0,
   QUDA_CONTRACT_GAMMA_G1 = 1,
@@ -580,6 +616,12 @@ typedef enum QudaGaugeSmearType_s {
   QUDA_GAUGE_SMEAR_INVALID = QUDA_INVALID_ENUM
 } QudaGaugeSmearType;
 
+typedef enum QudaWFlowType_s {
+  QUDA_WFLOW_TYPE_WILSON,
+  QUDA_WFLOW_TYPE_SYMANZIK,
+  QUDA_WFLOW_TYPE_INVALID = QUDA_INVALID_ENUM
+} QudaWFlowType;
+
 typedef enum QudaFermionSmearType_s {
   QUDA_FERMION_SMEAR_TYPE_GAUSSIAN,
   QUDA_FERMION_SMEAR_TYPE_WUPPERTAL,

include/enum_quda_fortran.h

@@ -23,6 +23,7 @@
 #define QUDA_MOMENTUM_LINKS 3
 #define QUDA_COARSE_LINKS   4
 #define QUDA_SMEARED_LINKS  5
+#define QUDA_TWOLINK_LINKS 6
 
 #define QUDA_WILSON_LINKS         QUDA_SU3_LINKS
 #define QUDA_ASQTAD_FAT_LINKS     QUDA_GENERAL_LINKS
@@ -477,10 +478,37 @@
 #define QUDA_STAGGERED_PHASE_TIFR   3
 #define QUDA_STAGGERED_PHASE_INVALID QUDA_INVALID_ENUM
 
+#define QudaSpinTasteGamma integer(4)
+#define QUDA_SPIN_TASTE_G1 0
+#define QUDA_SPIN_TASTE_GX 1
+#define QUDA_SPIN_TASTE_GY 2
+#define QUDA_SPIN_TASTE_GZ 4
+#define QUDA_SPIN_TASTE_GT 8
+#define QUDA_SPIN_TASTE_G5 15
+#define QUDA_SPIN_TASTE_GYGZ 6
+#define QUDA_SPIN_TASTE_GZGX 5
+#define QUDA_SPIN_TASTE_GXGY 3
+#define QUDA_SPIN_TASTE_GXGT 9
+#define QUDA_SPIN_TASTE_GYGT 10
+#define QUDA_SPIN_TASTE_GZGT 12
+#define QUDA_SPIN_TASTE_G5GX 14
+#define QUDA_SPIN_TASTE_G5GY 13
+#define QUDA_SPIN_TASTE_G5GZ 11
+#define QUDA_SPIN_TASTE_G5GT 7
+#define QUDA_SPIN_TASTE_INVALID QUDA_INVALID_ENUM
+
 #define QudaContractType integer(4)
-#define QUDA_CONTRACT_TYPE_OPEN ,
-#define QUDA_CONTRACT_TYPE_DR ,
-#define QUDA_CONTRACT_TYPE_INVALID = QUDA_INVALID_ENUM
+#define QUDA_CONTRACT_TYPE_STAGGERED_FT_T 0
+#define QUDA_CONTRACT_TYPE_DR_FT_T 1
+#define QUDA_CONTRACT_TYPE_DR_FT_Z 2
+#define QUDA_CONTRACT_TYPE_STAGGERED 3
+#define QUDA_CONTRACT_TYPE_DR 4
+#define QUDA_CONTRACT_TYPE_OPEN 5
+#define QUDA_CONTRACT_TYPE_OPEN_SUM_T 6
+#define QUDA_CONTRACT_TYPE_OPEN_SUM_Z 7
+#define QUDA_CONTRACT_TYPE_OPEN_FT_T 8
+#define QUDA_CONTRACT_TYPE_OPEN_FT_Z 9
+#define QUDA_CONTRACT_TYPE_INVALID QUDA_INVALID_ENUM
 
 #define QudaContractGamma integer(4)
 #define QUDA_CONTRACT_GAMMA_I 0

include/gamma.cuh

@@ -291,4 +291,123 @@ namespace quda {
     inline constexpr int Dir() const { return dir;  }
   };
 
+  // list of specialized structures used in the contraction kernels:
+
+  constexpr array<array<int, 4>, 16> get_dr_gm_i()
+  {
+    return {{// VECTORS
+             // G_idx = 1: \gamma_1
+             {3, 2, 1, 0},
+
+             // G_idx = 2: \gamma_2
+             {3, 2, 1, 0},
+
+             // G_idx = 3: \gamma_3
+             {2, 3, 0, 1},
+
+             // G_idx = 4: \gamma_4
+             {2, 3, 0, 1},
+
+             // PSEUDO-VECTORS
+             // G_idx = 6: \gamma_5\gamma_1
+             {3, 2, 1, 0},
+
+             // G_idx = 7: \gamma_5\gamma_2
+             {3, 2, 1, 0},
+
+             // G_idx = 8: \gamma_5\gamma_3
+             {2, 3, 0, 1},
+
+             // G_idx = 9: \gamma_5\gamma_4
+             {2, 3, 0, 1},
+
+             // SCALAR
+             // G_idx = 0: I
+             {0, 1, 2, 3},
+
+             // PSEUDO-SCALAR
+             // G_idx = 5: \gamma_5
+             {0, 1, 2, 3},
+
+             // TENSORS
+             // G_idx = 10: (i/2) * [\gamma_1, \gamma_2]
+             {0, 1, 2, 3},
+
+             // G_idx = 11: (i/2) * [\gamma_1, \gamma_3]. this matrix was corrected
+             {1, 0, 3, 2},
+
+             // G_idx = 12: (i/2) * [\gamma_1, \gamma_4]
+             {1, 0, 3, 2},
+
+             // G_idx = 13: (i/2) * [\gamma_2, \gamma_3]
+             {1, 0, 3, 2},
+
+             // G_idx = 14: (i/2) * [\gamma_2, \gamma_4]
+             {1, 0, 3, 2},
+
+             // G_idx = 15: (i/2) * [\gamma_3, \gamma_4]. this matrix was corrected
+             {0, 1, 2, 3}}};
+  }
+
+  template <typename T> constexpr array<array<complex<T>, 4>, 16> get_dr_g5gm_z()
+  {
+
+    constexpr complex<T> p_i = complex<T>(0., +1.);
+    constexpr complex<T> m_i = complex<T>(0., -1.);
+    constexpr complex<T> p_1 = complex<T>(+1., 0.);
+    constexpr complex<T> m_1 = complex<T>(-1., 0.);
+
+    return {{// VECTORS
+             // G_idx = 1: \gamma_1
+             {p_i, p_i, p_i, p_i},
+
+             // G_idx = 2: \gamma_2
+             {m_1, p_1, m_1, p_1},
+
+             // G_idx = 3: \gamma_3
+             {p_i, m_i, p_i, m_i},
+
+             // G_idx = 4: \gamma_4
+             {p_1, p_1, m_1, m_1},
+
+             // PSEUDO-VECTORS
+             // G_idx = 6: \gamma_5\gamma_1
+             {p_i, p_i, m_i, m_i},
+
+             // G_idx = 7: \gamma_5\gamma_2
+             {m_1, p_1, p_1, m_1},
+
+             // G_idx = 8: \gamma_5\gamma_3
+             {p_i, m_i, m_i, p_i},
+
+             // G_idx = 9: \gamma_5\gamma_4
+             {p_1, p_1, p_1, p_1},
+
+             // SCALAR
+             // G_idx = 0: I
+             {p_1, p_1, m_1, m_1},
+
+             // PSEUDO-SCALAR
+             // G_idx = 5: \gamma_5
+             {p_1, p_1, p_1, p_1},
+
+             // TENSORS
+             // G_idx = 10: (i/2) * [\gamma_1, \gamma_2]
+             {p_1, m_1, m_1, p_1},
+
+             // G_idx = 11: (i/2) * [\gamma_1, \gamma_3]. this matrix was corrected
+             {m_i, p_i, p_i, m_i},
+
+             // G_idx = 12: (i/2) * [\gamma_1, \gamma_4]
+             {m_1, m_1, m_1, m_1},
+
+             // G_idx = 13: (i/2) * [\gamma_2, \gamma_3]
+             {p_1, p_1, m_1, m_1},
+
+             // G_idx = 14: (i/2) * [\gamma_2, \gamma_4]
+             {m_i, p_i, m_i, p_i},
+
+             // G_idx = 15: (i/2) * [\gamma_3, \gamma_4]. this matrix was corrected
+             {m_1, p_1, m_1, p_1}}};
+  }
 } // namespace quda