Staggered fix finished

Grid/lattice/Lattice_reduction.h

@@ -142,15 +142,6 @@ inline typename vobj::scalar_objectD sumD(const vobj *arg, Integer osites)
   return sumD_cpu(arg,osites);
 #endif  
 }
-template<class vobj>
-inline typename vobj::scalar_objectD sumD_large(const vobj *arg, Integer osites)
-{
-#if defined(GRID_CUDA)||defined(GRID_HIP)
-  return sumD_gpu_large(arg,osites);
-#else
-  return sumD_cpu(arg,osites);
-#endif  
-}
 
 template<class vobj>
 inline typename vobj::scalar_object sum(const Lattice<vobj> &arg)

Grid/lattice/Lattice_reduction_gpu.h

@@ -23,7 +23,7 @@ unsigned int nextPow2(Iterator x) {
 }
 
 template <class Iterator>
-int getNumBlocksAndThreads(const Iterator n, const size_t sizeofsobj, Iterator &threads, Iterator &blocks) {
+void getNumBlocksAndThreads(const Iterator n, const size_t sizeofsobj, Iterator &threads, Iterator &blocks) {
   
   int device;
 #ifdef GRID_CUDA
@@ -37,13 +37,13 @@ int getNumBlocksAndThreads(const Iterator n, const size_t sizeofsobj, Iterator &
   Iterator sharedMemPerBlock   = gpu_props[device].sharedMemPerBlock;
   Iterator maxThreadsPerBlock  = gpu_props[device].maxThreadsPerBlock;
   Iterator multiProcessorCount = gpu_props[device].multiProcessorCount;
-  /*  
+  
   std::cout << GridLogDebug << "GPU has:" << std::endl;
   std::cout << GridLogDebug << "\twarpSize            = " << warpSize << std::endl;
   std::cout << GridLogDebug << "\tsharedMemPerBlock   = " << sharedMemPerBlock << std::endl;
   std::cout << GridLogDebug << "\tmaxThreadsPerBlock  = " << maxThreadsPerBlock << std::endl;
   std::cout << GridLogDebug << "\tmultiProcessorCount = " << multiProcessorCount << std::endl;
-  */  
+  
   if (warpSize != WARP_SIZE) {
     std::cout << GridLogError << "The warp size of the GPU in use does not match the warp size set when compiling Grid." << std::endl;
     exit(EXIT_FAILURE);
@@ -53,12 +53,12 @@ int getNumBlocksAndThreads(const Iterator n, const size_t sizeofsobj, Iterator &
   threads = warpSize;
   if ( threads*sizeofsobj > sharedMemPerBlock ) {
     std::cout << GridLogError << "The object is too large for the shared memory." << std::endl;
-    return 0;
+    exit(EXIT_FAILURE);
   }
   while( 2*threads*sizeofsobj < sharedMemPerBlock && 2*threads <= maxThreadsPerBlock ) threads *= 2;
   // keep all the streaming multiprocessors busy
   blocks = nextPow2(multiProcessorCount);
-  return 1;
+  
 }
 
 template <class sobj, class Iterator>
@@ -198,7 +198,7 @@ __global__ void reduceKernel(const vobj *lat, sobj *buffer, Iterator n) {
 // Possibly promote to double and sum
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites) 
+inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites) 
 {
   typedef typename vobj::scalar_objectD sobj;
   typedef decltype(lat) Iterator;
@@ -207,9 +207,7 @@ inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osi
   Integer size = osites*nsimd;
 
   Integer numThreads, numBlocks;
-  int ok = getNumBlocksAndThreads(size, sizeof(sobj), numThreads, numBlocks);
-  assert(ok);
-
+  getNumBlocksAndThreads(size, sizeof(sobj), numThreads, numBlocks);
   Integer smemSize = numThreads * sizeof(sobj);
 
   Vector<sobj> buffer(numBlocks);
@@ -220,62 +218,6 @@ inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osi
   auto result = buffer_v[0];
   return result;
 }
-
-template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osites)
-{
-  typedef typename vobj::vector_type  vector;
-  typedef typename vobj::scalar_typeD scalarD;
-  typedef typename vobj::scalar_objectD sobj;
-  sobj ret;
-  scalarD *ret_p = (scalarD *)&ret;
-  
-  const int words = sizeof(vobj)/sizeof(vector);
-
-  Integer nsimd= vobj::Nsimd();
-  Integer size = osites*nsimd;
-  Integer numThreads, numBlocks;
-
-  Vector<vector> buffer(osites);
-  vector *dat = (vector *)lat;
-  vector *buf = &buffer[0];
-  iScalar<vector> *tbuf =(iScalar<vector> *)  &buffer[0];
-  for(int w=0;w<words;w++) {
-
-    accelerator_for(ss,osites,1,{
-	buf[ss] = dat[ss*words+w];
-      });
-      
-    ret_p[w] = sumD_gpu_small(tbuf,osites);
-  }
-  return ret;
-}
-
-template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
-{
-  typedef typename vobj::vector_type  vector;
-  typedef typename vobj::scalar_typeD scalarD;
-  typedef typename vobj::scalar_objectD sobj;
-  sobj ret;
-  scalarD *ret_p = (scalarD *)&ret;
-  
-  const int words = sizeof(vobj)/sizeof(vector);
-
-  Integer nsimd= vobj::Nsimd();
-  Integer size = osites*nsimd;
-  Integer numThreads, numBlocks;
-  int ok = getNumBlocksAndThreads(size, sizeof(sobj), numThreads, numBlocks);
-  
-  if ( ok ) {
-    ret = sumD_gpu_small(lat,osites);
-  } else {
-    ret = sumD_gpu_large(lat,osites);
-  }
-  return ret;
-}
-
-
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Return as same precision as input performing reduction in double precision though
 /////////////////////////////////////////////////////////////////////////////////////////////////////////

Grid/qcd/action/fermion/WilsonCloverHelpers.h

@@ -726,8 +726,8 @@ public:
   static strong_inline void ApplyBoundaryMask(Field& f, const Mask& m) {
     conformable(f, m);
     auto grid  = f.Grid();
-    const uint32_t Nsite = grid->oSites();
-    const uint32_t Nsimd = grid->Nsimd();
+    const int Nsite = grid->oSites();
+    const int Nsimd = grid->Nsimd();
     autoView(f_v, f, AcceleratorWrite);
     autoView(m_v, m, AcceleratorRead);
     // NOTE: this function cannot be 'private' since nvcc forbids this for kernels

Grid/stencil/Stencil.h

@@ -240,6 +240,20 @@ public:
     cobj * mpi_p;
     Integer buffer_size;
   };
+  struct CopyReceiveBuffer {
+    void * from_p;
+    void * to_p;
+    Integer bytes;
+  };
+  struct CachedTransfer {
+    Integer direction;
+    Integer OrthogPlane;
+    Integer DestProc;
+    Integer bytes;
+    Integer lane;
+    Integer cb;
+    void *recv_buf;
+  };
 
 
 protected:
@@ -271,7 +285,8 @@ public:
   std::vector<Merge> MergersSHM;
   std::vector<Decompress> Decompressions;
   std::vector<Decompress> DecompressionsSHM;
-
+  std::vector<CopyReceiveBuffer> CopyReceiveBuffers ;
+  std::vector<CachedTransfer> CachedTransfers;
   ///////////////////////////////////////////////////////////
   // Unified Comms buffers for all directions
   ///////////////////////////////////////////////////////////
@@ -551,8 +566,62 @@ public:
     Mergers.resize(0);
     MergersSHM.resize(0);
     Packets.resize(0);
+    CopyReceiveBuffers.resize(0);
+    CachedTransfers.resize(0);
     calls++;
   }
+  void AddCopy(void *from,void * to, Integer bytes)
+  {
+    //    std::cout << "Adding CopyReceiveBuffer "<<std::hex<<from<<" "<<to<<std::dec<<" "<<bytes<<std::endl;
+    CopyReceiveBuffer obj;
+    obj.from_p = from;
+    obj.to_p = to;
+    obj.bytes= bytes;
+    CopyReceiveBuffers.push_back(obj);
+  }
+  void CommsCopy()
+  {
+    //    These are device resident MPI buffers.
+    for(int i=0;i<CopyReceiveBuffers.size();i++){
+      cobj *from=(cobj *)CopyReceiveBuffers[i].from_p;
+      cobj *to  =(cobj *)CopyReceiveBuffers[i].to_p;
+      Integer words = CopyReceiveBuffers[i].bytes/sizeof(cobj);
+      //    std::cout << "CopyReceiveBuffer "<<std::hex<<from<<" "<<to<<std::dec<<" "<<words*sizeof(cobj)<<std::endl;
+      accelerator_forNB(j, words, cobj::Nsimd(), {
+	  coalescedWrite(to[j] ,coalescedRead(from [j]));
+      });
+    }
+  }
+  
+  Integer CheckForDuplicate(Integer direction, Integer OrthogPlane, Integer DestProc, void *recv_buf,Integer lane,Integer bytes,Integer cb)
+  {
+    CachedTransfer obj;
+    obj.direction   = direction;
+    obj.OrthogPlane = OrthogPlane;
+    obj.DestProc    = DestProc;
+    obj.recv_buf    = recv_buf;
+    obj.lane        = lane;
+    obj.bytes       = bytes;
+    obj.cb          = cb;
+
+    for(int i=0;i<CachedTransfers.size();i++){
+      if (   (CachedTransfers[i].direction  ==direction)
+	   &&(CachedTransfers[i].OrthogPlane==OrthogPlane)
+	   &&(CachedTransfers[i].DestProc   ==DestProc)
+	   &&(CachedTransfers[i].bytes      ==bytes)
+	   &&(CachedTransfers[i].lane       ==lane)
+	   &&(CachedTransfers[i].cb         ==cb)
+	     ){
+	//	std::cout << "Found duplicate plane dir "<<direction<<" plane "<< OrthogPlane<< " simd "<<lane << " relproc "<<DestProc<< " bytes "<<bytes <<std::endl;
+	AddCopy(CachedTransfers[i].recv_buf,recv_buf,bytes);
+	return 1;
+      }
+    }
+
+    //    std::cout << "No duplicate plane dir "<<direction<<" plane "<< OrthogPlane<< " simd "<<lane << " relproc "<<DestProc<<"  bytes "<<bytes<<std::endl;
+    CachedTransfers.push_back(obj);
+    return 0;
+  }
   void AddPacket(void *xmit,void * rcv, Integer to,Integer from,Integer bytes){
     Packet p;
     p.send_buf = xmit;
@@ -578,6 +647,7 @@ public:
     mv.push_back(m);
   }
   template<class decompressor>  void CommsMerge(decompressor decompress)    {
+    CommsCopy();
     CommsMerge(decompress,Mergers,Decompressions);
   }
   template<class decompressor>  void CommsMergeSHM(decompressor decompress) {
@@ -590,8 +660,8 @@ public:
   }
 
   template<class decompressor>
-  void CommsMerge(decompressor decompress,std::vector<Merge> &mm,std::vector<Decompress> &dd) {
-
+  void CommsMerge(decompressor decompress,std::vector<Merge> &mm,std::vector<Decompress> &dd)
+  {
     
     mergetime-=usecond();
     for(int i=0;i<mm.size();i++){
@@ -1011,9 +1081,11 @@ public:
 
       int sx        = (x+sshift)%rd;
       int comm_proc = ((x+sshift)/rd)%pd;
-
+      
       if (comm_proc) {
 
+	
+	
 	int words = buffer_size;
 	if (cbmask != 0x3) words=words>>1;
 
@@ -1045,9 +1117,10 @@ public:
 	  recv_buf=this->u_recv_buf_p;
 	}
 
+
 	cobj *send_buf;
 	send_buf = this->u_send_buf_p; // Gather locally, must send
-
+	
 	////////////////////////////////////////////////////////
 	// Gather locally
 	////////////////////////////////////////////////////////
@@ -1056,23 +1129,27 @@ public:
 	Gather_plane_simple_table(face_table[face_idx],rhs,send_buf,compress,u_comm_offset,so); face_idx++;
 	gathertime+=usecond();
 
-	///////////////////////////////////////////////////////////
-	// Build a list of things to do after we synchronise GPUs
-	// Start comms now???
-	///////////////////////////////////////////////////////////
-	AddPacket((void *)&send_buf[u_comm_offset],
-		  (void *)&recv_buf[u_comm_offset],
-		  xmit_to_rank,
-		  recv_from_rank,
-		  bytes);
+	int duplicate = CheckForDuplicate(dimension,sx,comm_proc,(void *)&recv_buf[u_comm_offset],0,bytes,cbmask);
+	if ( (!duplicate) ) { // Force comms for now
 
+	  ///////////////////////////////////////////////////////////
+	  // Build a list of things to do after we synchronise GPUs
+	  // Start comms now???
+	  ///////////////////////////////////////////////////////////
+	  AddPacket((void *)&send_buf[u_comm_offset],
+		    (void *)&recv_buf[u_comm_offset],
+		    xmit_to_rank,
+		    recv_from_rank,
+		    bytes);
+	}
+	
 	if ( compress.DecompressionStep() ) {
 	  AddDecompress(&this->u_recv_buf_p[u_comm_offset],
 			&recv_buf[u_comm_offset],
 			words,Decompressions);
 	}
 	u_comm_offset+=words;
-      }
+	}
     }
     return 0;
   }
@@ -1181,8 +1258,10 @@ public:
 
 	    rpointers[i] = rp;
 
-	    AddPacket((void *)sp,(void *)rp,xmit_to_rank,recv_from_rank,bytes);
-
+	    int duplicate = CheckForDuplicate(dimension,sx,nbr_proc,(void *)rp,i,bytes,cbmask);
+	    if ( (!duplicate) ) { // Force comms for now
+	      AddPacket((void *)sp,(void *)rp,xmit_to_rank,recv_from_rank,bytes);
+	    }
 
 	  } else {
 

Grid/util/Init.cc

@@ -534,7 +534,6 @@ void Grid_init(int *argc,char ***argv)
 void Grid_finalize(void)
 {
 #if defined (GRID_COMMS_MPI) || defined (GRID_COMMS_MPI3) || defined (GRID_COMMS_MPIT)
-  MPI_Barrier(MPI_COMM_WORLD);
   MPI_Finalize();
   Grid_unquiesce_nodes();
 #endif