Hip compile

CMakeLists.txt

@@ -122,8 +122,8 @@ set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS DEVEL RELEASE STRICT DEBUG
 string(REPLACE ":" ";" LIBRARY_DIRS $ENV{LD_LIBRARY_PATH})
 
 set(DEFAULT_GPU_ARCH sm_60)
-set(GPU_ARCH ${DEFAULT_GPU_ARCH} CACHE STRING "set the GPU architecture (sm_20, sm_21, sm_30, sm_35, sm_37, sm_50, sm_52, sm_60, sm_70)")
-set_property(CACHE GPU_ARCH PROPERTY STRINGS sm_20 sm_21 sm_30 sm_35 sm_37 sm_50 sm_52 sm_60 sm_70)
+set(GPU_ARCH ${DEFAULT_GPU_ARCH} CACHE STRING "set the GPU architecture (sm_20, sm_21, sm_30, sm_35, sm_37, sm_50, sm_52, sm_60, sm_70, sm_80, sm_90)")
+set_property(CACHE GPU_ARCH PROPERTY STRINGS sm_20 sm_21 sm_30 sm_35 sm_37 sm_50 sm_52 sm_60 sm_70 sm_80 sm_90)
 
 # max threads
 set(PLEGMA_MAX_THREADS 256 CACHE INT "maximum number of threads per block in tuner")

include/PLEGMA_BLAS.h

@@ -14,6 +14,7 @@
 #include <cublas_v2.h>
 #include <mpi.h>
 #include <PLEGMA_utils.h>
+#include <quda_api.h>
 #pragma once
 enum OPER_MATR_BLAS {NOTRANS, TRANS, DAGGER};
 namespace cBLAS{
@@ -233,7 +234,7 @@ namespace cuBLAS{
     if(trans == NOTRANS) PLEGMA_error("Use gemv without MPI comm");
     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);
+    qudaMemcpy(yHost,y,n*2*sizeof(Float),qudaMemcpyDeviceToHost);
     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{
   inline void gemv(OPER_MATR_BLAS trans, int m, int n, Float alpha[2], Float* A, Float* x, Float beta[2], Float* y, MPI_Comm comm){
     Float yHost[n*2];
     cuBLAS::gemv(trans,m, n, alpha, A, x, beta, y, yHost,comm);
-    cudaMemcpy(y,yHost,sizeof(yHost),cudaMemcpyHostToDevice);
+    qudaMemcpy(y,yHost,sizeof(yHost),qudaMemcpyHostToDevice);
     checkQudaError();
   }
 }

include/PLEGMA_Correlator.h

@@ -8,6 +8,7 @@ namespace plegma {
 
   // forward declaration
   template<typename Float>  class PLEGMA_Vector;
+  template<typename Float>  class PLEGMA_Vector3D;
   template<typename Float>  class PLEGMA_Propagator;
   template<typename Float>  class PLEGMA_Propagator3D;
   template<typename Float>  class PLEGMA_Fmunu;
@@ -191,6 +192,21 @@ namespace plegma {
       std::vector<std::string> d = {s};
       setGroups(d);
     }
+    void contractEigVecs(Float *ptr, int nvecs, size_t vec_size, bool dev_ptr);
+
+    void contractPropEigVecs(PLEGMA_Propagator<Float> &prop1, Float *spinVals,
+			     Float *ptr, int nvecs, size_t vec_size, bool dev_ptr);
+
+    void contractPropEigVecsClosed(PLEGMA_Propagator<Float> &prop1, Float *spinVals,
+		    Float *ptr, int nvecs, size_t vec_size, bool dev_ptr);
+
+    void contractLoop(PLEGMA_Vector3D<Float> &vec,
+		    PLEGMA_Propagator3D<Float> &prop);
+
+    void contractLoopSIB(PLEGMA_Vector3D<Float> &vec,
+		    PLEGMA_Propagator3D<Float> &prop);
+
+
     void contractMesons(PLEGMA_Propagator<Float> &prop1,
 			PLEGMA_Propagator<Float> &prop2);
 
@@ -201,6 +217,47 @@ namespace plegma {
 
     void contractMesonsNew(PLEGMA_Propagator<Float> &prop1,
 			   PLEGMA_Propagator<Float> &prop2);
+
+    void contractMesonsOpen(PLEGMA_Propagator<Float> &prop1,
+			    PLEGMA_Propagator<Float> &prop2,
+			    bool all_cols=true);
+
+    void contractMesons1ps(PLEGMA_Propagator<Float> &prop1,
+                           PLEGMA_Propagator<Float> &prop2,
+                           PLEGMA_Gauge<Float> &gauge,
+                           bool all_cols=true);
+
+    void contractMesonsSIB(PLEGMA_Propagator<Float> &prop1,
+                           PLEGMA_Propagator<Float> &prop2);
+
+    void contractMesonsOpenSIB(PLEGMA_Propagator<Float> &prop1,
+                               PLEGMA_Propagator<Float> &prop2);
+
+    void contractMesonsSIR(PLEGMA_Propagator<Float> &prop1,
+                           PLEGMA_Propagator<Float> &prop2);
+
+    void contractMesonsOpenSIR(PLEGMA_Propagator<Float> &prop1,
+                               PLEGMA_Propagator<Float> &prop2);
+
+    void contractMesonsOpenDefl(PLEGMA_Propagator<Float> &prop1,
+                                PLEGMA_Propagator<Float> &prop2);
+
+    void contractMesons1psSIB(PLEGMA_Propagator<Float> &prop1,
+                              PLEGMA_Propagator<Float> &prop2,
+                              PLEGMA_Gauge<Float> &gauge);
+
+    void contractMesonsLIBE(PLEGMA_Propagator<Float> &prop1,
+                           PLEGMA_Propagator<Float> &prop2);
+
+    void contractMesonsOpenLIBE(PLEGMA_Propagator<Float> &prop1,
+                               PLEGMA_Propagator<Float> &prop2);
+
+    void contractMesons1psLIBE(PLEGMA_Propagator<Float> &prop1,
+                              PLEGMA_Propagator<Float> &prop2,
+                              PLEGMA_Gauge<Float> &gauge);
+
+
+
     void contractMesonsAll(PLEGMA_Propagator<Float> &prop1,
 			   PLEGMA_Propagator<Float> &prop2);
 
@@ -211,13 +268,17 @@ namespace plegma {
 			     PLEGMA_Propagator<Float> &prop2,
 			     PLEGMA_Propagator<Float> &prop3,
 			     PLEGMA_Propagator<Float> &prop4);
-
+
+    void contractBaryonsEEE(Float *ptr, Float* evs, int nvecs, size_t vec_size, bool dev_ptr);
+
+
     void contractBaryonsUDSC(PLEGMA_Propagator<Float> &propUP,
-			     PLEGMA_Propagator<Float> &propDN, 
-			     PLEGMA_Propagator<Float> &propST, 
-			     PLEGMA_Propagator<Float> &propCH,
-			     bool only_st=false, bool only_ch=false);
-
+                             PLEGMA_Propagator<Float> &propDN,
+                             PLEGMA_Propagator<Float> &propST,
+                             PLEGMA_Propagator<Float> &propCH,
+           bool only_up = false, bool only_dn = false,
+           bool only_st = false, bool only_ch = false, bool excludeHeavyOnly = false);
+
     template<typename Float2>
     void contractTetraquarks(PLEGMA_Propagator<Float2> &propLT,
 			     PLEGMA_Propagator<Float2> &propST, 

include/PLEGMA_Field.h

@@ -27,9 +27,12 @@ namespace plegma {
     size_t ghost_length; /*!< Member variable to hold the size of ghosts that involve in the communication (only side ghosts) */
     size_t ghost_corner_length; /*!< Member variable to hold the size of ghosts that involve in the communication (ghosts which are on the corners) */
     size_t ghost_vertex_length;
+    size_t single_ghost_length;
+    size_t single_corner_length;
 
     Float *h_elem; /*!< Member variable pointer to the elements of the field on CPU */
     Float *d_elem; /*!< Member variable pointer to the elements of the field on GPU */
+    Float *d_ext_ghost;
     Float *h_ext_ghost_r; /*!< Member variable pointer to the side ghost elements of the field on CPU (receive version)*/
     Float *h_ext_ghost_s; /*!< Member variable pointer to the side ghost elements of the field on CPU (send version)*/
     Float *h_ext_ghost_corner_r; /*!< Member variable pointer to the corner ghost elements of the field on CPU (receive version)*/
@@ -72,6 +75,7 @@ namespace plegma {
        @param field_l: see "field_length"
        @param vol_l: see "total_length"
      */
+
     void initialize(ALLOCATION_FLAG alloc_flag, int field_l, size_t vol_l);
   public:
     /**
@@ -126,6 +130,8 @@ namespace plegma {
      * @return a pointer to access device elements of the field
      */
     Float* D_elem() const { return d_elem; }
+    void D_elem(Float* ptr) { d_elem = ptr; }
+
     /**
      * @return a boolean if the host memory is allocated
      */
@@ -134,6 +140,7 @@ namespace plegma {
      * @return a boolean if the device memory is allocated
      */
     bool IsAllocDevice() const { return isAllocDevice;}
+
     /**
      * @return the kind of allocation we have for the field
      */
@@ -151,8 +158,11 @@ namespace plegma {
      * @return the length of the ghost part of the field in lattice points (local)
      */
     size_t Ghost_length() const { return ghost_length;} // the length of the ghost
+    size_t SingleGhost_length() const { return single_ghost_length;} // the length of a single ghost
     size_t GhostCorner_length() const { return ghost_corner_length;} // the length of the ghost for corners
+    size_t SingleCorner_length() const { return single_corner_length;} // the length of a single ghost for corners
     size_t GhostVertex_length() const { return ghost_vertex_length;} // the length of the ghost for vertex
+    size_t TotalGhost_length() const { return Ghost_length()+GhostCorner_length()+GhostVertex_length();}
     /**
      * @return the length of the field + ghost in lattice points (local)
      */
@@ -162,9 +172,12 @@ namespace plegma {
      * @return the bytes the length of the field including d.o.f
      */
     size_t Bytes_total() const { return this->Total_length()*this->Field_length()*2*sizeof(Float); }
+    size_t Bytes_total_ghost() const { return this->TotalGhost_length()*this->Field_length()*2*sizeof(Float); }
     /**
      * @return the bytes the length of the ghost including d.o.f
      */
+    size_t Bytes_singleghost() const { return this->SingleGhost_length()*this->Field_length()*2*sizeof(Float); }
+    size_t Bytes_singleCorner() const { return this->SingleCorner_length()*this->Field_length()*2*sizeof(Float); }
     size_t Bytes_ghost() const { return this->Ghost_length()*this->Field_length()*2*sizeof(Float); }
     size_t Bytes_ghostCorner() const { return this->GhostCorner_length()*this->Field_length()*2*sizeof(Float); }
     size_t Bytes_ghostVertex() const { return this->GhostVertex_length()*this->Field_length()*2*sizeof(Float); }
@@ -207,11 +220,16 @@ namespace plegma {
      * @param dirOr: choose the dir,orien. If negative does all dir, orien. If >=0 then (0,1,2,3,4,5,6,7,8) -> (+x,+y,+z,+t,-x,-y,-z,-t)
      */
     void communicateSideGhost(short dir=-1, ORIENTATION sign=DIR_BOTH, ACTION action=DO_ALL);
+    void communicateSecondSideGhost(short dir=-1, ORIENTATION sign=DIR_BOTH, ACTION action=DO_ALL);
+    void communicateThirdSideGhost(short dir=-1, ORIENTATION sign=DIR_BOTH, ACTION action=DO_ALL);
+
     /**
      * @brief Communicates side+corner ghosts in chosen direction,orientation
      * @param dirOr: choose the dir,orien. If negative does all dir, orien. If >=0 then (0,1,2,3,4,5,6,7,8) -> (+x,+y,+z,+t,-x,-y,-z,-t)
      */
     void communicateCornerGhost(short dir=-1, ORIENTATION sign=DIR_BOTH, ACTION action=DO_ALL);
+    void communicateSecondCornerGhost(short dir=-1, ORIENTATION sign=DIR_BOTH, ACTION action=DO_ALL);
+
     void communicateVertexGhost(short dir=-1, ORIENTATION sign=DIR_BOTH, ACTION action=DO_ALL);
     /**
      * @brief Communicates side or side+corner ghosts in chosen direction,orientation give the ghost type

include/PLEGMA_Gauge.h

@@ -1,4 +1,5 @@
 #include <PLEGMA_Field.h>
+#include <PLEGMA_GaugeU1.h>
 #include "complex"
 
 #ifndef _PLEGMA_GAUGE_H
@@ -120,6 +121,8 @@ namespace plegma {
      **/
     void gluonField(PLEGMA_Gauge<Float> &uIn);
     void U3xU1(PLEGMA_Gauge<Float> &u3, PLEGMA_U1Gauge<Float> &u1);
+    void mul_dag(PLEGMA_Gauge<Float> &uIn);
+    void qedPhase(PLEGMA_GaugeU1<Float> &uIn, Float phase);
   };
 
   /////////////////////////////////////

include/PLEGMA_GaugeU1.h

@@ -0,0 +1,40 @@
+#include <PLEGMA_Field.h>
+#include "complex"
+
+#ifndef _PLEGMA_GAUGEU1_H
+#define _PLEGMA_GAUGEU1_H
+
+namespace plegma {
+  template<typename Float> class PLEGMA_Su3field;
+  ////////////////////////
+  // CLASS: PLEGMA_GaugeU1 //
+  ////////////////////////
+
+  template<typename Float>
+  class PLEGMA_GaugeU1 : virtual public PLEGMA_Field<Float> {
+  public:
+    PLEGMA_GaugeU1(ALLOCATION_FLAG alloc_flag=BOTH, GHOST_FLAG ghost_flag=FIRST_CORNER);
+    ~PLEGMA_GaugeU1(){;}
+    Float calculatePlaq(Float phase);
+   /* 
+    void absorbDir_device(PLEGMA_Su3field<Float> &su,int dir);
+    void absorbDir_host(PLEGMA_Su3field<Float> &su,int dir);
+
+
+    //PLEGMA_topocharge.cuh
+    Float calculateTopo(TOPO_CHARGE_DEF charge_def);
+
+    //PLEGMA_WFlow.cuh
+    void GFlow_step( PLEGMA_GaugeU1<Float> &Z, double eps );
+    void applyGradientFlow( PLEGMA_GaugeU1<Float> &Z, int N, double eps);
+    void unitarize();
+
+    void scaleDirWise(std::complex<Float> scale[N_DIMS]);
+    void momPhase(Float phase[N_DIMS],int mom[N_DIMS]);
+    void gFixingLandau(PLEGMA_GaugeU1<Float> &uIn,Float overelaxPar=0.2,Float tolerance=1.0e-12,int maxIter=10000, int seedOverRelax=123456);
+    void gluonField(PLEGMA_GaugeU1<Float> &uIn);
+  */
+  };
+}
+
+#endif

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();
@@ -53,6 +58,36 @@ namespace plegma {
        @return void
      **/    
     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
+       @param PLEGMA_Propagator<Float> in initial propagator to be broadcasted
+       @param int sinktimeslice  the sinktime slice to be broadcasted
+     **/ 
+    void pack_propagator_as_sink(PLEGMA_Propagator<Float> &in, int sinktimeslice, int max_source_sink_separations, bool initial);
+
+    /**
+       @brief packing the sinktime slice in all other time slices
+       @param the initial propagator to be broadcasted
+       @param PLEGMA_Propagator<Float> in initial propagator to be broadcasted
+       @param int sinktimeslice  the sinktime slice to be broadcasted
+     **/
+    void pack_propagator_from_source_to_sink(PLEGMA_Propagator<Float> &in, int sinktimeslice, int source_sink_separations, bool initial);
+
+
+
 
     void applyBoundaries_device(int t0);
     void rotateToPhysicalBase_host(int sign);