cylqcd · pittlerf · Apr 13, 2023 · Apr 13, 2023 · Apr 13, 2023 · May 2, 2023
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
 set(CMAKE_INCLUDE_CURRENT_DIR ON)
 set(CMAKE_INCLUDE_DIRECTORIES_PROJECT_BEFORE ON)
 set(CMAKE_COLOR_MAKEFILE ON)
-set(CMAKE_CUDA_SEPARABLE_COMPILATION ON)
+set(CMAKE_HIP_SEPARABLE_COMPILATION ON)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 set(CMAKE_CXX_STANDARD 17)
 include(CheckIncludeFileCXX)
@@ -424,13 +424,21 @@ if (CMAKE_CXX_COMPILER_ID MATCHES "Clang")
     endif()
 endif()
 
+set(CPU_ARCH ${CMAKE_HOST_SYSTEM_PROCESSOR})
+if(${CPU_ARCH} STREQUAL "x86_64")
+  set(CXX_OPT "-mtune=native")
+elseif(${CPU_ARCH} STREQUAL "ppc64le")
+  set(CXX_OPT "-Ofast -mcpu=native -mtune=native")
+endif()
+
+
 
 #define CXX FLAGS
 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} -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} -std=c++17 -O3 -w" CACHE STRING
+set(CMAKE_CXX_FLAGS_RELEASE "${OpenMP_CXX_FLAGS} -std=c++17 -O3 -w ${CXX_OPT}" CACHE STRING
     "Flags used by the C++ compiler during release builds.")
 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.")
@@ -443,10 +451,11 @@ OPTION(USE_OpenMP "Use OpenMP" ON)
 IF(USE_OpenMP)
     FIND_PACKAGE(OpenMP)
     IF( OPENMP_FOUND )
-        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+	 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -fPIC")
    	set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
     	set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
         set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -Xcompiler ${OpenMP_CXX_FLAGS}")
+        set(CMAKE_HIP_FLAGS "${CMAKE_HIP_FLAGS} -Xcompiler ${OpenMP_CXX_FLAGS}")
     ENDIF()
 ENDIF()
 

diff --git a/include/PLEGMA_Hprobing.h b/include/PLEGMA_Hprobing.h
@@ -97,7 +97,7 @@ namespace plegma {
 	arrlc = new int[Nc];
       }
       catch (const std::bad_alloc& err) {
-	PLEGMA_error(err.what());
+	PLEGMA_error("%s\n",err.what());
       }
       createElemColBlock();
       createColLattice();

diff --git a/include/PLEGMA_Propagator.h b/include/PLEGMA_Propagator.h
@@ -1,4 +1,5 @@
 #include <PLEGMA_Field.h>
+#include <color_spinor_field.h>
 
 #ifndef _PLEGMA_PROPAGATOR_H
 #define _PLEGMA_PROPAGATOR_H
@@ -17,6 +18,10 @@ namespace plegma {
   public:
     PLEGMA_Propagator(ALLOCATION_FLAG alloc_flag=BOTH, GHOST_FLAG ghost_flag=FIRST_SIDE);
     ~PLEGMA_Propagator(){;}
+
+    void copyToQUDA( std::vector<ColorSpinorField>& cudaVector,  bool isEv = false);
+    void copyFromQUDA( std::vector<ColorSpinorField>& cudaVector, bool isEv = false);
+
 
     void apply_gamma(GAMMAS gMat, LEFTRIGHT LR = LEFT);
     void apply_gamma5();
@@ -54,6 +59,17 @@ namespace plegma {
      **/    
     void absorb(PLEGMA_Vector3D<Float> &vec, int global_it, int nu, int c2);
 
+/**
+       @brief Applies N times Gaussian(Wuppertal) smearing operator on all time-slices of a vector. NOTE: works also for Vector3D
+       @param PLEGMA_Vector<Float> &vecIn, The 4D input vector (Exchange of boundaries happens inside the function)
+       @param PLEGMA_Gauge<Float> &gauge, The gauge field that will be used in the Gaussian smearing operator (Exchange of boundaries happens inside the function)
+       @param int nsmearGauss, The number of times to apply the operator (if zero copies inVec to outVec)
+       @param Float alphaGauss, alpha parameter of the Gaussian smearing
+    **/
+
+    void gaussianSmearing(PLEGMA_Propagator<Float> &propIn, PLEGMA_Gauge<Float> &gauge, int nsmearGauss, Float alphaGauss);
+
+
     /** 
        @brief packing the sinktime slice in all other time slices
        @param the initial propagator to be broadcasted

diff --git a/include/PLEGMA_Vector.h b/include/PLEGMA_Vector.h
@@ -25,8 +25,8 @@ namespace plegma {
     PLEGMA_Vector(ALLOCATION_FLAG alloc_flag=BOTH, GHOST_FLAG ghost_flag=FIRST_SIDE);
     ~PLEGMA_Vector(){;}
 
-    void copyToQUDA( quda::ColorSpinorField *cudaVector, bool isEv = false);
-    void copyFromQUDA( quda::ColorSpinorField *cudaVector, bool isEv = false);
+    void copyToQUDA( std::vector<quda::ColorSpinorField>& cudaVector, bool isEv = false);
+    void copyFromQUDA( std::vector<quda::ColorSpinorField>& cudaVector, bool isEv = false);
 
     /**
        @brief Applies N times Gaussian(Wuppertal) smearing operator on all time-slices of a vector. NOTE: works also for Vector3D
@@ -67,6 +67,9 @@ namespace plegma {
        @return void
      **/    
     void absorb(PLEGMA_Propagator<Float> &prop, int nu , int c2);
+
+    void absorb(PLEGMA_Propagator<Float> *prop, int nu , int c2);
+
 
     void dilutespin(PLEGMA_Vector<Float> &vecIn, int spin);
 
@@ -108,8 +111,8 @@ namespace plegma {
 
   };
 
-  template<typename Float> void copyToQUDA(quda::ColorSpinorField *cudaVector, Float* delem, bool isEv = false); // delem is a device pointer
-  template<typename Float> void copyFromQUDA(Float* delem, quda::ColorSpinorField *cudaVector, bool isEv = false);
+  template<typename Float> void copyToQUDA(std::vector<quda::ColorSpinorField>& cudaVector, Float* delem, bool isEv = false); // delem is a device pointer
+  template<typename Float> void copyFromQUDA(Float* delem, std::vector<quda::ColorSpinorField>& cudaVector,  bool isEv = false);
 
   /////////////////////////////////////
   // CLASS: PLEGMA_Vector3D ///////////

diff --git a/include/PLEGMA_global.h b/include/PLEGMA_global.h
@@ -15,7 +15,7 @@
 #include <quda.h>
 #include <quda_internal.h>
 #ifndef __HIP__
-#include <cuda_fp16.h>
+//##include <cuda_fp16.h>
 #endif
 //#include <cublas_v2.h>
 #ifdef __GNUG__ // gnu C++ compiler

diff --git a/include/PLEGMA_utils.h b/include/PLEGMA_utils.h
@@ -14,6 +14,32 @@
 
 using namespace plegma;
 
+   /**
+     @brief Helper function for signalling that DD is not supported
+     @param[in] a Input field
+     @return true if all fields are supported
+   */
+  template <typename T1>
+  inline bool Nrhs_(const char *func, const char *file, int line, const T1 &a_)
+  {
+    const unwrap_t<T1> &a(a_);
+    if (a->Nrhs() != 12) errorQuda("Propagator solve is not supported (%s:%d in %s())", file, line, func);
+    return true;
+  }
+  /**
+     @brief Helper function for signalling that DD is not supported
+     @param[in] a Input field
+     @param[in] args List of additional fields to check domain decomposition on
+     @return true if all fields match
+   */
+  template <typename T1, typename... Args>
+  inline bool Nrhs_(const char *func, const char *file, int line, const T1 &a, const Args &...args)
+  {
+    // checking all possible pairs
+    return (Nrhs_(func, file, line, a) && Nrhs_(func, file, line, args...));
+  }
+  #define assertNrhs(...) Nrhs_(__func__, __FILE__, __LINE__, __VA_ARGS__)
+
 //============ QUDA_interface.cpp ===============================//
 void initComms(int argc, char **argv, const int *commDims);
 void finalizeComms();
@@ -23,6 +49,8 @@ void finalizeGaugeQuda();
 void plaqQuda();
 void gFixingLandauOVR_QUDA(PLEGMA_Gauge<double> &gaugeOut, PLEGMA_Gauge<double> &gaugeIn, int type, double overelaxPar=1.5, double tolerance=1e-12,
 			   int maxiter=10000,int verbosePerSteps=1, int reunit_interval=1, int stop_theta=0);
+void gSmear_QUDA(PLEGMA_Gauge<double> &gaugeOut,PLEGMA_Gauge<double> &gaugeIn, bool antiperiodic);
+
 
 //============= QUDA_params.cpp ===================================//
 void infoQuda();

diff --git a/include/global/PLEGMA_global_constants.h b/include/global/PLEGMA_global_constants.h
@@ -22,18 +22,16 @@
 #endif
 #ifdef ADD_TO_GLOBAL
 
-void *symbolAddress;
 #define global_host(dtype, name, ...)					\
   HGC_global_vars.add<dtype>(#name,					\
 				   PRODUCT(__VA_ARGS__),		\
 				   &HGC_##name PARENTHESES(0,__VA_ARGS__))
 
 #define global_both(dtype, name, ...)					\
-  HIP_CHECK(hipGetSymbolAddress(&symbolAddress, HIP_SYMBOL(DGC_##name))); \
   HGC_global_vars.add<dtype>(#name,					\
 				   PRODUCT(__VA_ARGS__),		\
 				   &HGC_##name PARENTHESES(0,__VA_ARGS__), \
-			           (void**) &symbolAddress)
+			           (void**) &DGC_##name)
 
 #else
 #ifdef ALLOCATE

diff --git a/include/global/PLEGMA_prints.hpp b/include/global/PLEGMA_prints.hpp
@@ -50,7 +50,7 @@ extern class Options * HGC_options;
     HGC_options->checkErrors();						\
     if (HGC_init_PLEGMA_flag) {						\
       fprintf(getOutputFile(), "%sERROR: ", getOutputPrefix());		\
-      fprintf(getOutputFile(), __VA_ARGS__);				\
+      fprintf(getOutputFile(),  __VA_ARGS__);				\
       fprintf(getOutputFile(), " (rank %d, host %s, " __FILE__ ":%d in %s())\n", \
              comm_rank(), comm_hostname(), __LINE__, __func__);         \
       fprintf(getOutputFile(), "%s       last kernel called was (name=%s,volume=%s,aux=%s)\n", \

diff --git a/include/global/PLEGMA_structs.h b/include/global/PLEGMA_structs.h
@@ -2,6 +2,7 @@
 #define HIP_CHECK(error)                                                                           \
   {                                                                                                \
     hipError_t localError = error;                                                                 \
+    PLEGMA_printf("Error in MemcpytoSymbol %d\n", error);					   \
     if ((localError != hipSuccess) && (localError != hipErrorPeerAccessAlreadyEnabled)) {          \
       PLEGMA_printf("Error: %s, Code %d\n", hipGetErrorString(localError), localError);            \
       PLEGMA_printf("FILE: %s, LINE %d\n", __FILE__, __LINE__ );			           \
@@ -60,6 +61,8 @@ struct tex_mom_list {
   inline __device__ int4 get(const size_t &i) const {
 #ifdef __NVCC__
     return tex1Dfetch<int4>(tex,i);
+#elif defined (__HIP__)
+    return tex1Dfetch<int4>(tex,i);
 #else
     return make_int4(0,0,0,0);
 #endif
@@ -87,8 +90,15 @@ struct pointer_holder {
     if (err != cudaSuccess) {
         errorQuda("Failed to copy constant host memory of size to device %zu \n", size);
     }
-#else
-    HIP_CHECK(hipMemcpy(*devPointer, hostPointer, bytes*size, hipMemcpyHostToDevice));
+#elif defined(__HIP__)
+    hipError_t err = hipMemcpyToSymbol( *devPointer, hostPointer, bytes*size);
+    if (err != hipSuccess) {
+        errorQuda("Failed to copy constant host memory of size to device %zu \n", size);
+    }
+    else{
+        printf("No error in hipMemcpyToSymbol %d\n",err);
+    }
+  //HIP_CHECK(hipMemcpyToSymbol(*devPointer, hostPointer, bytes*size));
 #endif
     }
   }
@@ -100,13 +110,16 @@ struct pointer_holder {
         errorQuda("Failed to copy constant host memory of size to device %zu \n", size);
       }
 #else
-     qudaMemcpy(hostPointer, *devPointer,  bytes*size,qudaMemcpyDeviceToHost);
+      hipError_t err = hipMemcpyFromSymbol(hostPointer, *devPointer, bytes*size, 0, hipMemcpyDeviceToHost);
+      if (err != hipSuccess) {
+        errorQuda("Failed to copy constant host memory of size to device %zu \n", size);
+      }
 #endif
     }
   }
   bool checkDeviceConstant() {
-    if(false and devPointer != nullptr) {
-      char tmp[bytes*size];
+    if( devPointer != nullptr) {
+     char tmp[bytes*size];
       memcpy(tmp,hostPointer,bytes*size);
       copyFromDeviceConstant();
       bool check=true;
@@ -145,8 +158,10 @@ struct global_vars {
       globals[i].copyToDeviceConstant();
     }
   }
-  bool check() {
+  bool check() { 
     for(size_t i = 0; i < globals.size(); i++) {
+
+      std::string line = "HGC_" + globals[i].get_value();
       if(globals[i].checkDeviceConstant() == false) return false;
     }
     return true;
@@ -159,7 +174,7 @@ struct global_vars {
       PLEGMA_printf("%s",line.c_str());
     }
     PLEGMA_printf("\nGlobal constants available on both, host and device:\n");
-    for(int i = 0; i < globals.size(); i++) {
+    for(int i = 0; i < (int)globals.size(); i++) {
       if(globals[i].devPointer == nullptr) continue;
       if(globals[i].checkDeviceConstant()) {
 	std::string line = "H/DGC_" + globals[i].get_value();

diff --git a/include/io/PLEGMA_lime.h b/include/io/PLEGMA_lime.h
@@ -11,7 +11,7 @@ inline bool exists_file (const char* name) {
   return ( access( name, F_OK ) != -1 );
 }
 
-static void swap_8(double *Rd, int N)
+static inline void swap_8(double *Rd, int N)
 {
    char *i,*j,*k;
    char swap;
@@ -32,7 +32,7 @@ static void swap_8(double *Rd, int N)
    }
 }
 
-static void swap_4(float *Rd, int N)
+static inline void swap_4(float *Rd, int N)
 {
   char *i,*j,*k;
   char swap;
@@ -49,7 +49,7 @@ static void swap_4(float *Rd, int N)
   }
 }
 
-static int isBigEndian()
+static inline int isBigEndian()
 {
    union{
      char C[4];

diff --git a/include/utils/QUDA_interface.h b/include/utils/QUDA_interface.h
@@ -17,7 +17,7 @@ namespace quda {
     QudaMultigridParam mg_param;
     Dirac *D, *DSloppy, *DPre;
     DiracMatrix *M, *MSloppy, *MPre;
-    ColorSpinorField *b, *x;
+    std::vector<ColorSpinorField> b, x;
 #ifdef QUDA_INCLUDES_COMMIT_775a033
     QudaEigParam *mg_eig_param;
 #endif
@@ -26,15 +26,32 @@ namespace quda {
     QudaInvertParam getInvParams() const{return inv_param;}
     SolverParam* getSolverParam() const{return solverParam;}
     void UpdateSolver();
-    QUDA_solver(double mu);
+    void UpdateGaugeSolver();
+
+    double Mu() const { return inv_param.mu;}
+    int Nrhs()  const { return inv_param.num_src;}
+
+    QUDA_solver(double mu, int nsrc);
     virtual ~QUDA_solver();
-    ColorSpinorField* solve(ColorSpinorField * rhs);
-    template<typename Float>
-    ColorSpinorField* solve(PLEGMA_Vector<Float> &vectorIn);
-    template<typename Float>
-    void solve(PLEGMA_Vector<Float> &out, PLEGMA_Vector<Float> &in);
-    template<typename Float>
-    void runOneIter(PLEGMA_Vector<Float> &out, PLEGMA_Vector<Float> &in);
+    std::vector<ColorSpinorField> solve(std::vector<ColorSpinorField>& rhs);
+
+    template<typename Float> std::vector<ColorSpinorField> solve(PLEGMA_Vector<Float> &vectorIn);
+    template<typename Float> std::vector<ColorSpinorField> solve(PLEGMA_Propagator<Float> &vectorIn);
+
+
+    template<typename Float> void solve(PLEGMA_Propagator<Float> &out, PLEGMA_Propagator<Float> &in);
+    template<typename Float> void solve(PLEGMA_Vector<Float> &out, PLEGMA_Vector<Float> &in);
+
+    template<typename Float> void runOneIter( PLEGMA_Vector<Float> &out, PLEGMA_Vector<Float> &in);
+    template<typename Float> void runOneIter( PLEGMA_Propagator<Float> &out, PLEGMA_Propagator<Float> &in);
+
+
+    template<bool bl, typename Float>
+    std::vector<ColorSpinorField> solve(typename std::conditional<bl==true, PLEGMA_Vector<Float>,  PLEGMA_Propagator<Float>>::type &vectorIn);
+    template<bool bl, typename Float>
+    void solve(typename std::conditional<bl==true, PLEGMA_Vector<Float>,  PLEGMA_Propagator<Float>>::type &out, typename std::conditional<bl==true, PLEGMA_Vector<Float>,  PLEGMA_Propagator<Float>>::type &in);
+    template<bool bl, typename Float>
+    void runOneIter(typename std::conditional<bl==true, PLEGMA_Vector<Float>,  PLEGMA_Propagator<Float>>::type &out, typename std::conditional<bl==true, PLEGMA_Vector<Float>,  PLEGMA_Propagator<Float>>::type &in);
   };
 
   enum APP_TYPE {M,Mdag,MdagM,MMdag};
@@ -43,7 +60,7 @@ namespace quda {
     DiracParam dParam;
     QudaInvertParam inv_param;
     Dirac *D;
-    ColorSpinorField *in, *out;
+    std::vector<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
@@ -52,8 +69,10 @@ namespace quda {
     void print(){dParam.print();}
     void switchMu(double mu);
     void switchKappa(double kappa);
+    template<APP_TYPE type, bool bl, typename Float> void apply(typename std::conditional<bl==true, PLEGMA_Vector<Float>,  PLEGMA_Propagator<Float>>::type &out, typename std::conditional<bl==true, PLEGMA_Vector<Float>,  PLEGMA_Propagator<Float>>::type &in, QudaMassNormalization normType = QUDA_KAPPA_NORMALIZATION); // Default it is without any normalization
     template<APP_TYPE type, typename Float> void apply(PLEGMA_Vector<Float> &out, PLEGMA_Vector<Float> &in, QudaMassNormalization normType = QUDA_KAPPA_NORMALIZATION); // Default it is without any normalization
+    template<APP_TYPE type, typename Float> void apply(PLEGMA_Propagator<Float> &out, PLEGMA_Propagator<Float> &in, QudaMassNormalization normType = QUDA_KAPPA_NORMALIZATION); // Default it is without any normalization
+
     template<APP_TYPE type, typename Float> void apply(Float *dout, Float *din, QudaMassNormalization normType = QUDA_KAPPA_NORMALIZATION); // Default is without any normalization // Note that dout and din are device pointers
   };
 }
-