Merging with current QUDA develop

CMakeLists.txt

@@ -5,7 +5,7 @@ set(CMAKE_INCLUDE_DIRECTORIES_PROJECT_BEFORE ON)
 set(CMAKE_COLOR_MAKEFILE ON)
 set(CMAKE_CUDA_SEPARABLE_COMPILATION ON)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
-set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD 17)
 include(CheckIncludeFileCXX)
 include(CheckSymbolExists)
 
@@ -33,6 +33,16 @@ endif()
 FIND_LIBRARY(QUDA_LIB quda ${QUDA_HOME}/lib)
 include_directories(SYSTEM ${QUDA_HOME}/include)
 
+set(QUDA_PRECISION
+  "14"
+  CACHE STRING "which precisions to instantiate in QUDA (4-bit number - double, single, half, quarter)")
+set(QUDA_RECONSTRUCT
+  "7"
+  CACHE STRING "which reconstructs to instantiate in QUDA (3-bit number - 18, 13/12, 9/8)")
+add_definitions(-DQUDA_PRECISION=${QUDA_PRECISION})
+add_definitions(-DQUDA_RECONSTRUCT=${QUDA_RECONSTRUCT})
+
+
 set(QUDA_SRC "" CACHE PATH "path to QUDA src directory")
 if("${QUDA_SRC}" STREQUAL "")
   message( FATAL_ERROR "QUDA_SRC must be defined" )
@@ -149,7 +159,6 @@ set(PLEGMA_LIGHT_BARYONS FALSE CACHE BOOL "compile contractions for all light ba
 if(PLEGMA_LIGHT_BARYONS)
   add_definitions(-DPLEGMA_LIGHT_BARYONS)
 endif()
-
 set(PLEGMA_SCATTERING_CONTRACTIONS FALSE CACHE BOOL "compile contractions for pi-N scattering")
 if (PLEGMA_SCATTERING_CONTRACTIONS)
   add_definitions(-DPLEGMA_SCATTERING_CONTRACTIONS)
@@ -264,7 +273,7 @@ find_package(PythonInterp)
 if(${CMAKE_VERSION} VERSION_GREATER 3.7.99)
   find_package(CUDAWrapper)
   set(USING_CUDA_LANG_SUPPORT True)
-  set(CMAKE_CUDA_STANDARD 14)
+  set(CMAKE_CUDA_STANDARD 17)
   set(CMAKE_CUDA_STANDARD_REQUIRED True)
 else()
   set(CUDA_HOST_COMPILER "${CMAKE_CXX_COMPILER}" CACHE FILEPATH "Host side compiler used by NVCC")
@@ -415,7 +424,7 @@ add_definitions(-D__COMPUTE_CAPABILITY__=${COMP_CAP})
 
 # NVCC FLAGS independet off build type
 if(NOT USING_CUDA_LANG_SUPPORT)
-  set(CUDA_NVCC_FLAGS "-std c++11 -arch=${GPU_ARCH} -ftz=true -prec-div=false -prec-sqrt=false")
+  set(CUDA_NVCC_FLAGS "-std c++17 -arch=${GPU_ARCH} -ftz=true -prec-div=false -prec-sqrt=false")
 else()
   set(CUDA_NVCC_FLAGS "-ftz=true -prec-div=false -prec-sqrt=false")
   set(CMAKE_CUDA_FLAGS "-arch=${GPU_ARCH}" CACHE STRING "Flags used by the CUDA compiler" FORCE)
@@ -471,17 +480,17 @@ set(CMAKE_CUDA_FLAGS_DEBUG "${CUDA_NVCC_FLAGS} -g -DHOST_DEBUG -G" CACHE STRING
     FORCE )
 
 #define CXX FLAGS
-set(CMAKE_CXX_FLAGS_DEVEL  "${OpenMP_CXX_FLAGS} -O3 -Wall ${CLANG_FORCE_COLOR}" CACHE STRING
+set(CMAKE_CXX_FLAGS_DEVEL "${OpenMP_CXX_FLAGS} -std=c++17 -O3 -Wall ${CLANG_FORCE_COLOR}" CACHE STRING
 "Flags used by the C++ compiler during regular development builds.")
-set(CMAKE_CXX_FLAGS_STRICT  "${OpenMP_CXX_FLAGS} -O3 -Wall -Werror ${CLANG_NOERROR}" CACHE STRING
+set(CMAKE_CXX_FLAGS_STRICT  "${OpenMP_CXX_FLAGS} -std=c++17 -O3 -Wall -Werror ${CLANG_NOERROR}" CACHE STRING
 "Flags used by the C++ compiler during strict jenkins builds.")
-set(CMAKE_CXX_FLAGS_RELEASE "${OpenMP_CXX_FLAGS} -O3 -w" CACHE STRING
+set(CMAKE_CXX_FLAGS_RELEASE "${OpenMP_CXX_FLAGS} -std=c++17 -O3 -w" CACHE STRING
     "Flags used by the C++ compiler during release builds.")
-set(CMAKE_CXX_FLAGS_HOSTDEBUG "${OpenMP_CXX_FLAGS} -Wall -Wno-unknown-pragmas -g -fno-inline -DHOST_DEBUG ${CLANG_FORCE_COLOR}" CACHE STRING
+set(CMAKE_CXX_FLAGS_HOSTDEBUG "${OpenMP_CXX_FLAGS} -std=c++17 -Wall -Wno-unknown-pragmas -g -fno-inline -DHOST_DEBUG ${CLANG_FORCE_COLOR}" CACHE STRING
     "Flags used by the C++ compiler during host-debug builds.")
-set(CMAKE_CXX_FLAGS_DEVICEDEBUG "${OpenMP_CXX_FLAGS} -Wall -Wno-unknown-pragmas ${CLANG_FORCE_COLOR}" CACHE STRING
+set(CMAKE_CXX_FLAGS_DEVICEDEBUG "${OpenMP_CXX_FLAGS} -std=c++17 -Wall -Wno-unknown-pragmas ${CLANG_FORCE_COLOR}" CACHE STRING
     "Flags used by the C++ compiler during device-debug builds.")
-set(CMAKE_CXX_FLAGS_DEBUG "${OpenMP_CXX_FLAGS} -Wall -Wno-unknown-pragmas -g -fno-inline -DHOST_DEBUG ${CLANG_FORCE_COLOR}" CACHE STRING
+set(CMAKE_CXX_FLAGS_DEBUG "${OpenMP_CXX_FLAGS} -std=c++17 -Wall -Wno-unknown-pragmas -g -fno-inline -DHOST_DEBUG ${CLANG_FORCE_COLOR}" CACHE STRING
     "Flags used by the C++ compiler during full (host+device) debug builds.")
 
 OPTION(USE_OpenMP "Use OpenMP" ON)
@@ -555,4 +564,3 @@ add_subdirectory(tests)
 ### add PLEGMA directory
 add_subdirectory(utils)
 add_subdirectory(plegma)
-

include/PLEGMA_BLAS.h

@@ -13,6 +13,7 @@
 
 #include <cublas_v2.h>
 #include <mpi.h>
+#include <PLEGMA_utils.h>
 #pragma once
 enum OPER_MATR_BLAS {NOTRANS, TRANS, DAGGER};
 namespace cBLAS{
@@ -233,7 +234,7 @@ namespace cuBLAS{
     if(comm == MPI_COMM_NULL) PLEGMA_error("Communicator is NULL and cannot be used for MPI reduction");
     cuBLAS::gemv_(trans, m, n, alpha, A, x, beta, y);
     cudaMemcpy(yHost,y,n*2*sizeof(Float),cudaMemcpyDeviceToHost);
-    checkCudaError();
+    checkQudaError();
     int mpiErr = MPI_Allreduce(MPI_IN_PLACE,yHost,n*2,MPI_Type(yHost),MPI_SUM,comm);
     if(mpiErr != MPI_SUCCESS) PLEGMA_error("MPI_Allreduce failed with error %d\n", mpiErr);
   }
@@ -245,7 +246,7 @@ namespace cuBLAS{
     Float yHost[n*2];
     cuBLAS::gemv(trans,m, n, alpha, A, x, beta, y, yHost,comm);
     cudaMemcpy(y,yHost,sizeof(yHost),cudaMemcpyHostToDevice);
-    checkCudaError();
+    checkQudaError();
   }
 }
 //=================================================================//

include/PLEGMA_Correlator.h

@@ -204,9 +204,13 @@ namespace plegma {
     void contractMesonsAll(PLEGMA_Propagator<Float> &prop1,
 			   PLEGMA_Propagator<Float> &prop2);
 
-
     void contractBaryons(PLEGMA_Propagator<Float> &prop1,
 			 PLEGMA_Propagator<Float> &prop2);
+
+    void contractBaryonsWall(PLEGMA_Propagator<Float> &prop1,
+			     PLEGMA_Propagator<Float> &prop2,
+			     PLEGMA_Propagator<Float> &prop3,
+			     PLEGMA_Propagator<Float> &prop4);
 
     void contractBaryonsUDSC(PLEGMA_Propagator<Float> &propUP,
 			     PLEGMA_Propagator<Float> &propDN, 
@@ -319,12 +323,31 @@ namespace plegma {
 					int signProps, std::vector<GAMMAS> gammas,
 					int z, std::string quark);
 
+    void contractNucleonThrp_staple(PLEGMA_Propagator<Float> &bwdProp,
+				    PLEGMA_Propagator<Float> &fwdProp,
+				    PLEGMA_Su3field<Float> &su3,
+				    int signProps, std::vector<GAMMAS> gammas, bool isZfac, int b, int l, int z);
+
     void contractNucleonThrp_qgq(PLEGMA_Propagator<Float> &bwdProp,
 				 PLEGMA_Propagator<Float> &fwdProp,
 				 PLEGMA_Su3field<Float> &su3_1,
 				 PLEGMA_Fmunu<Float> &Fmunu, std::pair<int,int> munu,
 				 PLEGMA_Su3field<Float> &su3_2,
 				 int signProps, std::vector<GAMMAS> gammas, int z2, int z1);
+    void contractMesonsFourp_ultralocal_oneendtrick(PLEGMA_Propagator<Float> &prop1,
+						    PLEGMA_Propagator<Float> &prop2,int b);
+
+    void contractMesonsFourp_ultralocal(PLEGMA_Propagator<Float> &prop1,
+					PLEGMA_Propagator<Float> &prop2,
+					PLEGMA_Propagator<Float> &prop3,
+					PLEGMA_Propagator<Float> &prop4,int b);
+
+    void contractTMDWFMesonsTrick_Zfac(PLEGMA_Propagator<Float> &prop1,
+				       PLEGMA_Su3field<float> &staple, int b, int l);
+
+    void contractTMDWFMesons_Zfac(PLEGMA_Propagator<Float> &prop1,
+				  PLEGMA_Propagator<Float> &prop2,
+                                  PLEGMA_Su3field<float> &staple, int b, int l);
 
     virtual void writeASCII(std::string filename) const;
     virtual void writeHDF5(std::string filename) const;

include/PLEGMA_Field.h

@@ -341,6 +341,9 @@ namespace plegma {
      * @param munu: the mu,nu pair of lorentz indices of the Plq
      */
     void trPmunu(PLEGMA_Gauge<Float> &gauge, std::pair<int,int> munu);
+
+    void SU3Trace(PLEGMA_Su3field<Float> &su3field);
+
     /**
        @brief Absorbs all  elements from a 3D field and puts it at a specific global time of the 4D field
        @param PLEGMA_Field3D<Float> prop, The 3D field

include/PLEGMA_Hprobing.h

@@ -1,4 +1,6 @@
 #pragma once
+//#include <PLEGMA_utils.h>
+#include <tune_quda.h>
 
 namespace plegma {
   inline int getVecToInd(std::vector<int> x, std::vector<int> L){
@@ -100,15 +102,15 @@ namespace plegma {
       createColLattice();
       //  if(check)checkColoring();
       cudaMalloc((void**)&d_arrVc, HGC_localVolume*sizeof(int));
-      checkCudaError();
+//      checkQudaError();
       cudaMemcpy(d_arrVc, h_arrVc, HGC_localVolume*sizeof(int), cudaMemcpyHostToDevice);
-      checkCudaError();    
+//      checkQudaError();    
     }
     ~PLEGMA_Hprobing(){
       delete[] h_arrVc;
       delete[] arrlc;
       cudaFree(d_arrVc);
-      checkCudaError();
+      //checkQudaError();
     }
     int* H_arrVc() const{return h_arrVc;}
     int* D_arrVc() const{return d_arrVc;}

include/PLEGMA_Propagator.h

@@ -57,6 +57,7 @@ namespace plegma {
     void applyBoundaries_device(int t0);
     void rotateToPhysicalBase_host(int sign);
     void rotateToPhysicalBase_device(int sign);
+    void PropmulVVdag(PLEGMA_Vector<Float> &vec1,PLEGMA_Vector<Float> &vec2);
   };
 
   ///////////////////////////////

include/PLEGMA_ScattCorrelator.h

@@ -190,7 +190,7 @@ namespace plegma {
     void T1( std::vector<GAMMAS_SCATT> &Gammas_i, std::vector<GAMMAS_SCATT> &Gammas_f, PLEGMA_Propagator<Float> &S1, PLEGMA_Propagator<Float> &S2, PLEGMA_Propagator<Float> &S3 );
     /**
      *  @brief performs T2 type reduction to compute baryon 2pt functions
-     *  T2_{alpha,beta}=\epsilon_{a,b,c}\epsilon_{l,m,n}S2^{c,l}_{alpha,beta}\Gamma_{i}_{alpha0,alpha1}S2^{b,m}_{beta0,alpha1}\Gamma_{f}_{beta0,beta1}S3^{a,n}_{beta1,alpha0} 
+     *  T2_{alpha,beta}=\epsilon_{a,b,c}\epsilon_{l,m,n}S1^{c,l}_{alpha,beta}\Gamma_{i}_{alpha0,alpha1}S2^{b,m}_{beta0,alpha1}\Gamma_{f}_{beta0,beta1}S3^{a,n}_{beta1,alpha0} 
      *  @param std::vector<GAMMAS_SCATT> &Gammas_i: list of gammas at the source for the contractions
      *  @param std::vector<GAMMAS_SCATT> &Gammas_f: list of gammas at the sink for the contractions
      *  @param PLEGMA_Propagator<Float> &S1: propagator to the sink spin from the source spin index

include/PLEGMA_global.h

@@ -43,6 +43,9 @@ using namespace std::chrono_literals;
 
 //======== Preprocessor macros =========//
 #include <global/PLEGMA_macros.hpp>
+#include <tune_quda.h>
+#include <comm_quda.h>
+using namespace quda;
 
 namespace plegma {
   //======== PLEGMA_printf, PLEGMA_error, PLEGMA_warning =========//

include/PLEGMA_utils.h

@@ -33,6 +33,7 @@ void setInvertParam(QudaInvertParam &inv_param);
 void setEigMultigridParam(QudaMultigridParam &mg_param, QudaEigParam *mg_eig_param);
 #endif
 //============== PLEGMA_utils.cpp =======================================//
+void checkQudaError();
 void createMom(int *Nmom, int momElem[][3], int Q_qs);
 extern const std::vector<std::string> listAvailOptPLEGMA;
 void initializeOptions(int argc, char **argv, bool withQuda=true, std::vector<std::string> listOptPLEGMA = listAvailOptPLEGMA);

include/global/PLEGMA_Options.h

@@ -1,5 +1,6 @@
 #pragma once
 #include "PLEGMA_templates.h"
+#include <comm_quda.h>
 struct argument{std::string name, value;};
 
 class Arguments{

include/global/PLEGMA_prints.hpp

@@ -8,7 +8,6 @@ extern bool HGC_hold_exit; // used to hold exit until all the errors have been p
 extern bool HGC_init_PLEGMA_flag;
 extern struct global_vars HGC_global_vars;
 extern class Options * HGC_options;
-
 #define PLEGMA_exit(value) do {			\
     if(! HGC_hold_exit) {			\
       if (HGC_init_PLEGMA_flag) {		\
@@ -53,7 +52,7 @@ extern class Options * HGC_options;
       fprintf(getOutputFile(), "%sERROR: ", getOutputPrefix());		\
       fprintf(getOutputFile(), __VA_ARGS__);				\
       fprintf(getOutputFile(), " (rank %d, host %s, " __FILE__ ":%d in %s())\n", \
-	      comm_rank(), comm_hostname(), __LINE__, __func__);	\
+             comm_rank(), comm_hostname(), __LINE__, __func__);         \
       fprintf(getOutputFile(), "%s       last kernel called was (name=%s,volume=%s,aux=%s)\n", \
 	      getOutputPrefix(), getLastTuneKey().name,			\
 	      getLastTuneKey().volume, getLastTuneKey().aux);		\

include/global/PLEGMA_structs.h

@@ -63,14 +63,18 @@ struct pointer_holder {
 
   void copyToDeviceConstant() {
     if(devPointer != nullptr) {
-      cudaMemcpyToSymbol( *devPointer, hostPointer, bytes*size);
-      checkCudaError();
+      cudaError_t err = cudaMemcpyToSymbol( *devPointer, hostPointer, bytes*size);
+      if (err != cudaSuccess) {
+        errorQuda("Failed to copy constant host memory of size to device %zu \n", size);
+      }
     }
   }
   void copyFromDeviceConstant() {
     if(false and devPointer != nullptr) {
-      cudaMemcpyFromSymbol(hostPointer, *devPointer, bytes*size, 0, cudaMemcpyDeviceToHost);
-      checkCudaError();
+      cudaError_t err = cudaMemcpyFromSymbol(hostPointer, *devPointer, bytes*size, 0, cudaMemcpyDeviceToHost);
+      if (err != cudaSuccess) {
+        errorQuda("Failed to copy constant host memory of size to device %zu \n", size);
+      }
     }
   }
   bool checkDeviceConstant() {

include/global/PLEGMA_templates.h

@@ -46,14 +46,18 @@ template<typename T> inline void hostFree(T &ptr) {
 
 // Pinned memory allocation
 template<typename T> inline void hostMallocPinned(T &ptr, size_t size){
-  cudaMallocHost((void**)&ptr, size);
-  checkCudaError();
+  cudaError_t err = cudaMallocHost((void**)&ptr, size);
+  if (err != cudaSuccess) {
+    errorQuda("Failed to allocate host memory of size %zu \n", size);
+  }
   HGC_used_memory += size;
 }
 
 template<typename T> inline void hostFreePinned(T &ptr, size_t size) {
-  cudaFreeHost(ptr);
-  checkCudaError();
+  cudaError_t err = cudaFreeHost(ptr);
+  if (err != cudaSuccess) {
+    errorQuda("Failed to free host memory of size %zu \n", size);
+  }
   ptr=NULL;
   HGC_used_memory -= size;
 }

include/io/PLEGMA_lime.h

@@ -2,7 +2,7 @@
 extern "C" {
 #include <lime.h>
 }
-
+#include <tune_quda.h>
 #include <PLEGMA_global.h>
 #include <unistd.h>
 #include <stdio.h>

include/utils/QUDA_interface.h

@@ -17,7 +17,7 @@ namespace quda {
     QudaMultigridParam mg_param;
     Dirac *D, *DSloppy, *DPre;
     DiracMatrix *M, *MSloppy, *MPre;
-    cudaColorSpinorField *b, *x;
+    ColorSpinorField *b, *x;
 #ifdef QUDA_INCLUDES_COMMIT_775a033
     QudaEigParam *mg_eig_param;
 #endif
@@ -28,9 +28,9 @@ namespace quda {
     void UpdateSolver();
     QUDA_solver(double mu);
     virtual ~QUDA_solver();
-    cudaColorSpinorField* solve(cudaColorSpinorField * rhs);
+    ColorSpinorField* solve(ColorSpinorField * rhs);
     template<typename Float>
-    cudaColorSpinorField* solve(PLEGMA_Vector<Float> &vectorIn);
+    ColorSpinorField* solve(PLEGMA_Vector<Float> &vectorIn);
     template<typename Float>
     void solve(PLEGMA_Vector<Float> &out, PLEGMA_Vector<Float> &in);
     template<typename Float>
@@ -43,7 +43,7 @@ namespace quda {
     DiracParam dParam;
     QudaInvertParam inv_param;
     Dirac *D;
-    cudaColorSpinorField *in, *out;
+    ColorSpinorField *in, *out;
     template<APP_TYPE type> void apply();
   public:
      //only QUDA_WILSON_DSLASH, QUDA_CLOVER_WILSON_DSLASH, QUDA_TWISTED_MASS_DSLASH, QUDA_TWISTED_CLOVER_DSLASH

include/utils/QUDA_types.h

@@ -284,8 +284,8 @@ inline QudaSchwarzType get_schwarz_type(std::string s) {
 inline QudaTwistFlavorType get_flavor_type(std::string s) {
   if(s=="singlet")
     return QUDA_TWIST_SINGLET;
-  else if(s=="deg-doublet")
-    return QUDA_TWIST_DEG_DOUBLET;
+//  else if(s=="deg-doublet")
+//    return QUDA_TWIST_DEG_DOUBLET;
   else if(s=="nondeg-doublet")
     return QUDA_TWIST_NONDEG_DOUBLET;
   else if(s=="no")
@@ -300,8 +300,8 @@ inline std::string get_flavor_str(QudaTwistFlavorType type) {
   switch(type) {
   case QUDA_TWIST_SINGLET:
     return "singlet";
-  case QUDA_TWIST_DEG_DOUBLET:
-    return "deg-doublet";
+//  case QUDA_TWIST_DEG_DOUBLET:
+//    return "deg-doublet";
   case QUDA_TWIST_NONDEG_DOUBLET:
     return "nondeg-doublet";
   case QUDA_TWIST_NO:
@@ -321,10 +321,10 @@ inline QudaInverterType get_solver_type(std::string s) {
     return QUDA_GCR_INVERTER;
   else if(s=="pcg")
     return QUDA_PCG_INVERTER;
-  else if(s=="mpcg")
-    return QUDA_MPCG_INVERTER; 
-  else if(s=="mpbicgstab")
-    return QUDA_MPBICGSTAB_INVERTER;
+//  else if(s=="mpcg")
+//    return QUDA_MPCG_INVERTER; 
+//  else if(s=="mpbicgstab")
+//    return QUDA_MPBICGSTAB_INVERTER;
   else if(s=="mr")
     return QUDA_MR_INVERTER;
   else if(s=="sd")
@@ -377,10 +377,10 @@ inline std::string get_solver_str(QudaInverterType type) {
     return "gcr";
   case QUDA_PCG_INVERTER:
     return "pcg";
-  case QUDA_MPCG_INVERTER:
-    return "mpcg";
-  case QUDA_MPBICGSTAB_INVERTER:
-    return "mpbicgstab";
+//  case QUDA_MPCG_INVERTER:
+//    return "mpcg";
+//  case QUDA_MPBICGSTAB_INVERTER:
+//    return "mpbicgstab";
   case QUDA_MR_INVERTER:
     return "mr";
   case QUDA_SD_INVERTER:

lib/CMakeLists.txt

@@ -18,6 +18,7 @@ set (PLEGMA_OBJS
   kernels/PLEGMA_threep_threeD_part3.cu
   kernels/PLEGMA_threep_threeD_part4.cu
   kernels/PLEGMA_threep_wilsonLine.cu
+  kernels/PLEGMA_threep_staple.cu
   kernels/PLEGMA_threep_qgq.cu	
   PLEGMA_QLoops.cu
   PLEGMA_FT.cu