Update to static data

Grid/DisableWarnings.h

@@ -44,22 +44,14 @@ directory
 #ifdef __NVCC__
  //disables nvcc specific warning in json.hpp
 #pragma clang diagnostic ignored "-Wdeprecated-register"
-
-#ifdef __NVCC_DIAG_PRAGMA_SUPPORT__
- //disables nvcc specific warning in json.hpp
-#pragma nv_diag_suppress unsigned_compare_with_zero
-#pragma nv_diag_suppress cast_to_qualified_type
- //disables nvcc specific warning in many files
-#pragma nv_diag_suppress esa_on_defaulted_function_ignored
-#pragma nv_diag_suppress extra_semicolon
-#else
- //disables nvcc specific warning in json.hpp
 #pragma diag_suppress unsigned_compare_with_zero
 #pragma diag_suppress cast_to_qualified_type
+
  //disables nvcc specific warning in many files
 #pragma diag_suppress esa_on_defaulted_function_ignored
 #pragma diag_suppress extra_semicolon
-#endif
+
+//Eigen only
 #endif
 
 // Disable vectorisation in Eigen on the Power8/9 and PowerPC

Grid/GridStd.h

@@ -16,7 +16,6 @@
 #include <functional>
 #include <stdio.h>
 #include <stdlib.h>
-#include <strings.h>
 #include <stdio.h>
 #include <signal.h>
 #include <ctime>

Grid/Grid_Eigen_Dense.h

@@ -14,11 +14,7 @@
 /* NVCC save and restore compile environment*/
 #ifdef __NVCC__
 #pragma push
-#ifdef __NVCC_DIAG_PRAGMA_SUPPORT__
-#pragma nv_diag_suppress code_is_unreachable
-#else
 #pragma diag_suppress code_is_unreachable
-#endif
 #pragma push_macro("__CUDA_ARCH__")
 #pragma push_macro("__NVCC__")
 #pragma push_macro("__CUDACC__")

Grid/algorithms/CoarsenedMatrix.h

@@ -262,7 +262,7 @@ public:
 	autoView( Tnp_v , (*Tnp), AcceleratorWrite);
 	autoView( Tnm_v , (*Tnm), AcceleratorWrite);
 	const int Nsimd = CComplex::Nsimd();
-	accelerator_for(ss, FineGrid->oSites(), Nsimd, {
+	accelerator_forNB(ss, FineGrid->oSites(), Nsimd, {
 	  coalescedWrite(y_v[ss],xscale*y_v(ss)+mscale*Tn_v(ss));
 	  coalescedWrite(Tnp_v[ss],2.0*y_v(ss)-Tnm_v(ss));
         });

Grid/algorithms/approx/Chebyshev.h

@@ -264,7 +264,7 @@ public:
       auto Tnp_v = Tnp->View();
       auto Tnm_v = Tnm->View();
       constexpr int Nsimd = vector_type::Nsimd();
-      accelerator_for(ss, in.Grid()->oSites(), Nsimd, {
+      accelerator_forNB(ss, in.Grid()->oSites(), Nsimd, {
 	  coalescedWrite(y_v[ss],xscale*y_v(ss)+mscale*Tn_v(ss));
 	  coalescedWrite(Tnp_v[ss],2.0*y_v(ss)-Tnm_v(ss));
       });

Grid/algorithms/iterative/Deflation.h

@@ -113,43 +113,7 @@ public:
     blockPromote(guess_coarse,guess,subspace);
     guess.Checkerboard() = src.Checkerboard();
   };
-
+};
-  void operator()(const std::vector<FineField> &src,std::vector<FineField> &guess) {
-    int Nevec = (int)evec_coarse.size();
-    int Nsrc = (int)src.size();
-    // make temp variables
-    std::vector<CoarseField> src_coarse(Nsrc,evec_coarse[0].Grid());
-    std::vector<CoarseField> guess_coarse(Nsrc,evec_coarse[0].Grid());    
-    //Preporcessing
-    std::cout << GridLogMessage << "Start BlockProject for loop" << std::endl;
-    for (int j=0;j<Nsrc;j++)
-    {
-    guess_coarse[j] = Zero();
-    std::cout << GridLogMessage << "BlockProject iter: " << j << std::endl;
-    blockProject(src_coarse[j],src[j],subspace);
-    }
-    //deflation set up for eigen vector batchsize 1 and source batch size equal number of sources
-    std::cout << GridLogMessage << "Start ProjectAccum for loop" << std::endl;
-    for (int i=0;i<Nevec;i++)
-    {
-      std::cout << GridLogMessage << "ProjectAccum Nvec: " << i << std::endl;
-      const CoarseField & tmp = evec_coarse[i];
-      for (int j=0;j<Nsrc;j++)
-      {
-        axpy(guess_coarse[j],TensorRemove(innerProduct(tmp,src_coarse[j])) / eval_coarse[i],tmp,guess_coarse[j]);
-      }
-    }
-    //postprocessing
-    std::cout << GridLogMessage << "Start BlockPromote for loop" << std::endl;
-    for (int j=0;j<Nsrc;j++)
-    {
-    std::cout << GridLogMessage << "BlockProject iter: " << j << std::endl;
-    blockPromote(guess_coarse[j],guess[j],subspace);
-    guess[j].Checkerboard() = src[j].Checkerboard();
-    }
-  };
-
-  };
 
 
 

Grid/allocator/MemoryManager.cc

@@ -40,7 +40,7 @@ void MemoryManager::PrintBytes(void)
 //////////////////////////////////////////////////////////////////////
 MemoryManager::AllocationCacheEntry MemoryManager::Entries[MemoryManager::NallocType][MemoryManager::NallocCacheMax];
 int MemoryManager::Victim[MemoryManager::NallocType];
-int MemoryManager::Ncache[MemoryManager::NallocType] = { 2, 8, 8, 16, 8, 16 };
+int MemoryManager::Ncache[MemoryManager::NallocType] = { 2, 8, 2, 8, 2, 8 };
 uint64_t MemoryManager::CacheBytes[MemoryManager::NallocType];
 //////////////////////////////////////////////////////////////////////
 // Actual allocation and deallocation utils

Grid/allocator/MemoryManager.h

@@ -36,11 +36,6 @@ NAMESPACE_BEGIN(Grid);
 
 #define GRID_ALLOC_SMALL_LIMIT (4096)
 
-#define STRINGIFY(x) #x
-#define TOSTRING(x) STRINGIFY(x)
-#define FILE_LINE __FILE__ ":" TOSTRING(__LINE__)
-#define AUDIT(a) MemoryManager::Audit(FILE_LINE)
-
 /*Pinning pages is costly*/
 ////////////////////////////////////////////////////////////////////////////
 // Advise the LatticeAccelerator class
@@ -97,9 +92,8 @@ private:
   static void *Insert(void *ptr,size_t bytes,AllocationCacheEntry *entries,int ncache,int &victim,uint64_t &cbytes) ;
   static void *Lookup(size_t bytes,AllocationCacheEntry *entries,int ncache,uint64_t &cbytes) ;
 
- public:
   static void PrintBytes(void);
-  static void Audit(std::string s);
+ public:
   static void Init(void);
   static void InitMessage(void);
   static void *AcceleratorAllocate(size_t bytes);
@@ -119,8 +113,11 @@ private:
   static uint64_t     DeviceToHostBytes;
   static uint64_t     HostToDeviceXfer;
   static uint64_t     DeviceToHostXfer;
-  static uint64_t     DeviceEvictions;
+  
-  static uint64_t     DeviceDestroy;
+  static uint64_t     DeviceAccesses;
+  static uint64_t     HostAccesses;
+  static uint64_t     DeviceAccessBytes;
+  static uint64_t     HostAccessBytes;
  
  private:
 #ifndef GRID_UVM
@@ -160,6 +157,7 @@ private:
 
   //  static void  LRUupdate(AcceleratorViewEntry &AccCache);
   static void  LRUinsert(AcceleratorViewEntry &AccCache);
+  static void  LRUinsertback(AcceleratorViewEntry &AccCache);
   static void  LRUremove(AcceleratorViewEntry &AccCache);
   
   // manage entries in the table
@@ -178,7 +176,6 @@ private:
 
  public:
   static void Print(void);
-  static void PrintAll(void);
   static void PrintState( void* CpuPtr);
   static int   isOpen   (void* CpuPtr);
   static void  ViewClose(void* CpuPtr,ViewMode mode);

Grid/allocator/MemoryManagerCache.cc

@@ -3,13 +3,8 @@
 
 #warning "Using explicit device memory copies"
 NAMESPACE_BEGIN(Grid);
-
+//#define dprintf(...) printf ( __VA_ARGS__ ); fflush(stdout);
-#define MAXLINE 512
+#define dprintf(...)
-static char print_buffer [ MAXLINE ];
-
-#define mprintf(...) snprintf (print_buffer,MAXLINE, __VA_ARGS__ ); std::cout << GridLogMemory << print_buffer;
-#define dprintf(...) snprintf (print_buffer,MAXLINE, __VA_ARGS__ ); std::cout << GridLogMemory << print_buffer;
-//#define dprintf(...) 
 
 
 ////////////////////////////////////////////////////////////
@@ -28,8 +23,11 @@ uint64_t  MemoryManager::HostToDeviceBytes;
 uint64_t  MemoryManager::DeviceToHostBytes;
 uint64_t  MemoryManager::HostToDeviceXfer;
 uint64_t  MemoryManager::DeviceToHostXfer;
-uint64_t  MemoryManager::DeviceEvictions;
+uint64_t  MemoryManager::DeviceAccesses;
-uint64_t  MemoryManager::DeviceDestroy;
+uint64_t  MemoryManager::HostAccesses;
+uint64_t  MemoryManager::DeviceAccessBytes;
+uint64_t  MemoryManager::HostAccessBytes;
+ 
 
 ////////////////////////////////////
 // Priority ordering for unlocked entries
@@ -93,6 +91,14 @@ void  MemoryManager::LRUinsert(AcceleratorViewEntry &AccCache)
   AccCache.LRU_valid = 1;
   DeviceLRUBytes+=AccCache.bytes;
 }
+void  MemoryManager::LRUinsertback(AcceleratorViewEntry &AccCache)
+{
+  assert(AccCache.LRU_valid==0);
+  LRU.push_back(AccCache.CpuPtr);
+  AccCache.LRU_entry = --LRU.end();
+  AccCache.LRU_valid = 1;
+  DeviceLRUBytes+=AccCache.bytes;
+}
 void  MemoryManager::LRUremove(AcceleratorViewEntry &AccCache)
 {
   assert(AccCache.LRU_valid==1);
@@ -111,17 +117,15 @@ void MemoryManager::AccDiscard(AcceleratorViewEntry &AccCache)
   ///////////////////////////////////////////////////////////
   assert(AccCache.state!=Empty);
   
-  mprintf("MemoryManager: Discard(%lx) %lx\n",(uint64_t)AccCache.CpuPtr,(uint64_t)AccCache.AccPtr); 
+   dprintf("MemoryManager: Discard(%llx) %llx\n",(uint64_t)AccCache.CpuPtr,(uint64_t)AccCache.AccPtr); 
   assert(AccCache.accLock==0);
   assert(AccCache.cpuLock==0);
   assert(AccCache.CpuPtr!=(uint64_t)NULL);
   if(AccCache.AccPtr) {
     AcceleratorFree((void *)AccCache.AccPtr,AccCache.bytes);
-    DeviceDestroy++;
     DeviceBytes   -=AccCache.bytes;
     LRUremove(AccCache);
-    AccCache.AccPtr=(uint64_t) NULL;
+    dprintf("MemoryManager: Free(%llx) LRU %lld Total %lld\n",(uint64_t)AccCache.AccPtr,DeviceLRUBytes,DeviceBytes);  
-    dprintf("MemoryManager: Free(%lx) LRU %ld Total %ld\n",(uint64_t)AccCache.AccPtr,DeviceLRUBytes,DeviceBytes);  
   }
   uint64_t CpuPtr = AccCache.CpuPtr;
   EntryErase(CpuPtr);
@@ -130,36 +134,27 @@ void MemoryManager::AccDiscard(AcceleratorViewEntry &AccCache)
 void MemoryManager::Evict(AcceleratorViewEntry &AccCache)
 {
   ///////////////////////////////////////////////////////////////////////////
-  // Make CPU consistent, remove from Accelerator, remove from LRU, LEAVE CPU only entry
+  // Make CPU consistent, remove from Accelerator, remove entry
-  // Cannot be acclocked. If allocated must be in LRU pool.
+  // Cannot be locked. If allocated must be in LRU pool.
-  //
-  // Nov 2022... Felix issue: Allocating two CpuPtrs, can have an entry in LRU-q with CPUlock.
-  //                          and require to evict the AccPtr copy. Eviction was a mistake in CpuViewOpen
-  //                          but there is a weakness where CpuLock entries are attempted for erase
-  //                          Take these OUT LRU queue when CPU locked?
-  //                          Cannot take out the table as cpuLock data is important.
   ///////////////////////////////////////////////////////////////////////////
   assert(AccCache.state!=Empty);
   
-  mprintf("MemoryManager: Evict cpu %lx acc %lx cpuLock %ld accLock %ld\n",
+  dprintf("MemoryManager: Evict(%llx) %llx\n",(uint64_t)AccCache.CpuPtr,(uint64_t)AccCache.AccPtr); 
-	  (uint64_t)AccCache.CpuPtr,(uint64_t)AccCache.AccPtr,
+  assert(AccCache.accLock==0);
-	  (uint64_t)AccCache.cpuLock,(uint64_t)AccCache.accLock); 
+  assert(AccCache.cpuLock==0);
-  assert(AccCache.accLock==0); // Cannot evict so logic bomb
+  
-  assert(AccCache.CpuPtr!=(uint64_t)NULL);
   if(AccCache.state==AccDirty) {
     Flush(AccCache);
   }
+  assert(AccCache.CpuPtr!=(uint64_t)NULL);
   if(AccCache.AccPtr) {
     AcceleratorFree((void *)AccCache.AccPtr,AccCache.bytes);
-    LRUremove(AccCache);
-    AccCache.AccPtr=(uint64_t)NULL;
-    AccCache.state=CpuDirty; // CPU primary now
     DeviceBytes   -=AccCache.bytes;
-    dprintf("MemoryManager: Free(%lx) footprint now %ld \n",(uint64_t)AccCache.AccPtr,DeviceBytes);  
+    LRUremove(AccCache);
+    dprintf("MemoryManager: Free(%llx) footprint now %lld \n",(uint64_t)AccCache.AccPtr,DeviceBytes);  
   }
-  //  uint64_t CpuPtr = AccCache.CpuPtr;
+  uint64_t CpuPtr = AccCache.CpuPtr;
-  DeviceEvictions++;
+  EntryErase(CpuPtr);
-  //  EntryErase(CpuPtr);
 }
 void MemoryManager::Flush(AcceleratorViewEntry &AccCache)
 {
@@ -169,7 +164,7 @@ void MemoryManager::Flush(AcceleratorViewEntry &AccCache)
   assert(AccCache.AccPtr!=(uint64_t)NULL);
   assert(AccCache.CpuPtr!=(uint64_t)NULL);
   acceleratorCopyFromDevice((void *)AccCache.AccPtr,(void *)AccCache.CpuPtr,AccCache.bytes);
-  mprintf("MemoryManager: Flush  %lx -> %lx\n",(uint64_t)AccCache.AccPtr,(uint64_t)AccCache.CpuPtr); fflush(stdout);
+  dprintf("MemoryManager: Flush  %llx -> %llx\n",(uint64_t)AccCache.AccPtr,(uint64_t)AccCache.CpuPtr); fflush(stdout);
   DeviceToHostBytes+=AccCache.bytes;
   DeviceToHostXfer++;
   AccCache.state=Consistent;
@@ -184,7 +179,7 @@ void MemoryManager::Clone(AcceleratorViewEntry &AccCache)
     AccCache.AccPtr=(uint64_t)AcceleratorAllocate(AccCache.bytes);
     DeviceBytes+=AccCache.bytes;
   }
-  mprintf("MemoryManager: Clone %lx <- %lx\n",(uint64_t)AccCache.AccPtr,(uint64_t)AccCache.CpuPtr); fflush(stdout);
+  dprintf("MemoryManager: Clone %llx <- %llx\n",(uint64_t)AccCache.AccPtr,(uint64_t)AccCache.CpuPtr); fflush(stdout);
   acceleratorCopyToDevice((void *)AccCache.CpuPtr,(void *)AccCache.AccPtr,AccCache.bytes);
   HostToDeviceBytes+=AccCache.bytes;
   HostToDeviceXfer++;
@@ -210,7 +205,6 @@ void MemoryManager::CpuDiscard(AcceleratorViewEntry &AccCache)
 void MemoryManager::ViewClose(void* Ptr,ViewMode mode)
 {
   if( (mode==AcceleratorRead)||(mode==AcceleratorWrite)||(mode==AcceleratorWriteDiscard) ){
-    dprintf("AcceleratorViewClose %lx\n",(uint64_t)Ptr);
     AcceleratorViewClose((uint64_t)Ptr);
   } else if( (mode==CpuRead)||(mode==CpuWrite)){
     CpuViewClose((uint64_t)Ptr);
@@ -222,7 +216,6 @@ void *MemoryManager::ViewOpen(void* _CpuPtr,size_t bytes,ViewMode mode,ViewAdvis
 {
   uint64_t CpuPtr = (uint64_t)_CpuPtr;
   if( (mode==AcceleratorRead)||(mode==AcceleratorWrite)||(mode==AcceleratorWriteDiscard) ){
-    dprintf("AcceleratorViewOpen %lx\n",(uint64_t)CpuPtr);
     return (void *) AcceleratorViewOpen(CpuPtr,bytes,mode,hint);
   } else if( (mode==CpuRead)||(mode==CpuWrite)){
     return (void *)CpuViewOpen(CpuPtr,bytes,mode,hint);
@@ -233,16 +226,13 @@ void *MemoryManager::ViewOpen(void* _CpuPtr,size_t bytes,ViewMode mode,ViewAdvis
 }
 void  MemoryManager::EvictVictims(uint64_t bytes)
 {
-  assert(bytes<DeviceMaxBytes);
   while(bytes+DeviceLRUBytes > DeviceMaxBytes){
     if ( DeviceLRUBytes > 0){
       assert(LRU.size()>0);
-      uint64_t victim = LRU.back(); // From the LRU
+      uint64_t victim = LRU.back();
       auto AccCacheIterator = EntryLookup(victim);
       auto & AccCache = AccCacheIterator->second;
       Evict(AccCache);
-    } else {
-      return;
     }
   }
 }
@@ -255,6 +245,9 @@ uint64_t MemoryManager::AcceleratorViewOpen(uint64_t CpuPtr,size_t bytes,ViewMod
     EntryCreate(CpuPtr,bytes,mode,hint);
   }
 
+  DeviceAccesses++;
+  DeviceAccessBytes+=bytes;
+
   auto AccCacheIterator = EntryLookup(CpuPtr);
   auto & AccCache = AccCacheIterator->second;
   if (!AccCache.AccPtr) {
@@ -265,12 +258,11 @@ uint64_t MemoryManager::AcceleratorViewOpen(uint64_t CpuPtr,size_t bytes,ViewMod
   assert(AccCache.cpuLock==0);  // Programming error
 
   if(AccCache.state!=Empty) {
-    dprintf("ViewOpen found entry %lx %lx : %ld %ld accLock %ld\n",
+    dprintf("ViewOpen found entry %llx %llx : %lld %lld\n",
 		    (uint64_t)AccCache.CpuPtr,
 		    (uint64_t)CpuPtr,
 		    (uint64_t)AccCache.bytes,
-	            (uint64_t)bytes,
+		    (uint64_t)bytes);
-		    (uint64_t)AccCache.accLock);
     assert(AccCache.CpuPtr == CpuPtr);
     assert(AccCache.bytes  ==bytes);
   }
@@ -305,7 +297,6 @@ uint64_t MemoryManager::AcceleratorViewOpen(uint64_t CpuPtr,size_t bytes,ViewMod
       AccCache.state  = Consistent; // Empty + AccRead => Consistent
     }
     AccCache.accLock= 1;
-    dprintf("Copied Empty entry into device accLock= %d\n",AccCache.accLock);
   } else if(AccCache.state==CpuDirty ){
     if(mode==AcceleratorWriteDiscard) {
       CpuDiscard(AccCache);
@@ -318,30 +309,28 @@ uint64_t MemoryManager::AcceleratorViewOpen(uint64_t CpuPtr,size_t bytes,ViewMod
       AccCache.state  = Consistent; // CpuDirty + AccRead => Consistent
     }
     AccCache.accLock++;
-    dprintf("CpuDirty entry into device ++accLock= %d\n",AccCache.accLock);
+    dprintf("Copied CpuDirty entry into device accLock %d\n",AccCache.accLock);
   } else if(AccCache.state==Consistent) {
     if((mode==AcceleratorWrite)||(mode==AcceleratorWriteDiscard))
       AccCache.state  = AccDirty;   // Consistent + AcceleratorWrite=> AccDirty
     else
       AccCache.state  = Consistent; // Consistent + AccRead => Consistent
     AccCache.accLock++;
-    dprintf("Consistent entry into device ++accLock= %d\n",AccCache.accLock);
+    dprintf("Consistent entry into device accLock %d\n",AccCache.accLock);
   } else if(AccCache.state==AccDirty) {
     if((mode==AcceleratorWrite)||(mode==AcceleratorWriteDiscard))
       AccCache.state  = AccDirty; // AccDirty + AcceleratorWrite=> AccDirty
     else
       AccCache.state  = AccDirty; // AccDirty + AccRead => AccDirty
     AccCache.accLock++;
-    dprintf("AccDirty entry ++accLock= %d\n",AccCache.accLock);
+    dprintf("AccDirty entry into device accLock %d\n",AccCache.accLock);
   } else {
     assert(0);
   }
 
-  assert(AccCache.accLock>0);
+  // If view is opened on device remove from LRU
-  // If view is opened on device must remove from LRU
   if(AccCache.LRU_valid==1){
     // must possibly remove from LRU as now locked on GPU
-    dprintf("AccCache entry removed from LRU \n");
     LRUremove(AccCache);
   }
 
@@ -362,12 +351,10 @@ void MemoryManager::AcceleratorViewClose(uint64_t CpuPtr)
   assert(AccCache.accLock>0);
 
   AccCache.accLock--;
+
   // Move to LRU queue if not locked and close on device
   if(AccCache.accLock==0) {
-    dprintf("AccleratorViewClose %lx AccLock decremented to %ld move to LRU queue\n",(uint64_t)CpuPtr,(uint64_t)AccCache.accLock);
     LRUinsert(AccCache);
-  } else {
-    dprintf("AccleratorViewClose %lx AccLock decremented to %ld\n",(uint64_t)CpuPtr,(uint64_t)AccCache.accLock);
   }
 }
 void MemoryManager::CpuViewClose(uint64_t CpuPtr)
@@ -379,6 +366,10 @@ void MemoryManager::CpuViewClose(uint64_t CpuPtr)
   assert(AccCache.accLock==0);
 
   AccCache.cpuLock--;
+
+  if(AccCache.cpuLock==0) {
+    LRUinsertback(AccCache);
+  }
 }
 /*
  *  Action  State   StateNext         Flush    Clone
@@ -401,13 +392,15 @@ uint64_t MemoryManager::CpuViewOpen(uint64_t CpuPtr,size_t bytes,ViewMode mode,V
     EntryCreate(CpuPtr,bytes,mode,transient);
   }
 
+  HostAccesses++;
+  HostAccessBytes+=bytes;
+  
   auto AccCacheIterator = EntryLookup(CpuPtr);
   auto & AccCache = AccCacheIterator->second;
 
-  // CPU doesn't need to free space
+  if (!AccCache.AccPtr) {
-  //  if (!AccCache.AccPtr) {
+     EvictVictims(bytes);
-  //    EvictVictims(bytes);
+  }
-  //  }
 
   assert((mode==CpuRead)||(mode==CpuWrite));
   assert(AccCache.accLock==0);  // Programming error
@@ -447,6 +440,12 @@ uint64_t MemoryManager::CpuViewOpen(uint64_t CpuPtr,size_t bytes,ViewMode mode,V
 
   AccCache.transient= transient? EvictNext : 0;
 
+  // If view is opened on host remove from LRU
+  // Host close says evict next from device
+  if(AccCache.LRU_valid==1){
+    LRUremove(AccCache);
+  }
+  
   return AccCache.CpuPtr;
 }
 void  MemoryManager::NotifyDeletion(void *_ptr)
@@ -461,28 +460,20 @@ void  MemoryManager::NotifyDeletion(void *_ptr)
 void  MemoryManager::Print(void)
 {
   PrintBytes();
-  std::cout << GridLogMessage << "--------------------------------------------" << std::endl;
+  std::cout << GridLogDebug << "--------------------------------------------" << std::endl;
-  std::cout << GridLogMessage << "Memory Manager                             " << std::endl;
+  std::cout << GridLogDebug << "Memory Manager                             " << std::endl;
-  std::cout << GridLogMessage << "--------------------------------------------" << std::endl;
+  std::cout << GridLogDebug << "--------------------------------------------" << std::endl;
-  std::cout << GridLogMessage << DeviceBytes   << " bytes allocated on device " << std::endl;
+  std::cout << GridLogDebug << DeviceBytes   << " bytes allocated on device " << std::endl;
-  std::cout << GridLogMessage << DeviceLRUBytes<< " bytes evictable on device " << std::endl;
+  std::cout << GridLogDebug << DeviceLRUBytes<< " bytes evictable on device " << std::endl;
-  std::cout << GridLogMessage << DeviceMaxBytes<< " bytes max on device       " << std::endl;
+  std::cout << GridLogDebug << DeviceMaxBytes<< " bytes max on device       " << std::endl;
-  std::cout << GridLogMessage << HostToDeviceXfer << " transfers        to   device " << std::endl;
+  std::cout << GridLogDebug << HostToDeviceXfer << " transfers        to   device " << std::endl;
-  std::cout << GridLogMessage << DeviceToHostXfer << " transfers        from device " << std::endl;
+  std::cout << GridLogDebug << DeviceToHostXfer << " transfers        from device " << std::endl;
-  std::cout << GridLogMessage << HostToDeviceBytes<< " bytes transfered to   device " << std::endl;
+  std::cout << GridLogDebug << HostToDeviceBytes<< " bytes transfered to   device " << std::endl;
-  std::cout << GridLogMessage << DeviceToHostBytes<< " bytes transfered from device " << std::endl;
+  std::cout << GridLogDebug << DeviceToHostBytes<< " bytes transfered from device " << std::endl;
-  std::cout << GridLogMessage << DeviceEvictions  << " Evictions from device " << std::endl;
+  std::cout << GridLogDebug << AccViewTable.size()<< " vectors " << LRU.size()<<" evictable"<< std::endl;
-  std::cout << GridLogMessage << DeviceDestroy    << " Destroyed vectors on device " << std::endl;
+  std::cout << GridLogDebug << "--------------------------------------------" << std::endl;
-  std::cout << GridLogMessage << AccViewTable.size()<< " vectors " << LRU.size()<<" evictable"<< std::endl;
+  std::cout << GridLogDebug << "CpuAddr\t\tAccAddr\t\tState\t\tcpuLock\taccLock\tLRU_valid "<<std::endl;
-  std::cout << GridLogMessage << "--------------------------------------------" << std::endl;
+  std::cout << GridLogDebug << "--------------------------------------------" << std::endl;
-}
-void  MemoryManager::PrintAll(void)
-{
-  Print();
-  std::cout << GridLogMessage << std::endl;
-  std::cout << GridLogMessage << "--------------------------------------------" << std::endl;
-  std::cout << GridLogMessage << "CpuAddr\t\tAccAddr\t\tState\t\tcpuLock\taccLock\tLRU_valid "<<std::endl;
-  std::cout << GridLogMessage << "--------------------------------------------" << std::endl;
   for(auto it=AccViewTable.begin();it!=AccViewTable.end();it++){
     auto &AccCache = it->second;
     
@@ -492,13 +483,13 @@ void  MemoryManager::PrintAll(void)
     if ( AccCache.state==AccDirty ) str = std::string("AccDirty");
     if ( AccCache.state==Consistent)str = std::string("Consistent");
 
-    std::cout << GridLogMessage << "0x"<<std::hex<<AccCache.CpuPtr<<std::dec
+    std::cout << GridLogDebug << "0x"<<std::hex<<AccCache.CpuPtr<<std::dec
 	      << "\t0x"<<std::hex<<AccCache.AccPtr<<std::dec<<"\t" <<str
 	      << "\t" << AccCache.cpuLock
 	      << "\t" << AccCache.accLock
 	      << "\t" << AccCache.LRU_valid<<std::endl;
   }
-  std::cout << GridLogMessage << "--------------------------------------------" << std::endl;
+  std::cout << GridLogDebug << "--------------------------------------------" << std::endl;
 
 };
 int   MemoryManager::isOpen   (void* _CpuPtr) 
@@ -512,61 +503,6 @@ int   MemoryManager::isOpen   (void* _CpuPtr)
     return 0;
   }
 }
-void MemoryManager::Audit(std::string s)
-{
-  uint64_t CpuBytes=0;
-  uint64_t AccBytes=0;
-  uint64_t LruBytes1=0;
-  uint64_t LruBytes2=0;
-  uint64_t LruCnt=0;
-  uint64_t LockedBytes=0;
-  
-  std::cout << " Memory Manager::Audit() from "<<s<<std::endl;
-  for(auto it=LRU.begin();it!=LRU.end();it++){
-    uint64_t cpuPtr = *it;
-    assert(EntryPresent(cpuPtr));
-    auto AccCacheIterator = EntryLookup(cpuPtr);
-    auto & AccCache = AccCacheIterator->second;
-    LruBytes2+=AccCache.bytes;
-    assert(AccCache.LRU_valid==1);
-    assert(AccCache.LRU_entry==it);
-  }
-  std::cout << " Memory Manager::Audit() LRU queue matches table entries "<<std::endl;
-  for(auto it=AccViewTable.begin();it!=AccViewTable.end();it++){
-    auto &AccCache = it->second;
-    
-    std::string str;
-    if ( AccCache.state==Empty    ) str = std::string("Empty");
-    if ( AccCache.state==CpuDirty ) str = std::string("CpuDirty");
-    if ( AccCache.state==AccDirty ) str = std::string("AccDirty");
-    if ( AccCache.state==Consistent)str = std::string("Consistent");
-
-    CpuBytes+=AccCache.bytes;
-    if( AccCache.AccPtr )    AccBytes+=AccCache.bytes;
-    if( AccCache.LRU_valid ) LruBytes1+=AccCache.bytes;
-    if( AccCache.LRU_valid ) LruCnt++;
-    
-    if ( AccCache.cpuLock || AccCache.accLock ) {
-      assert(AccCache.LRU_valid==0);
-      std::cout << GridLogError << s<< "\n\t 0x"<<std::hex<<AccCache.CpuPtr<<std::dec
-		<< "\t0x"<<std::hex<<AccCache.AccPtr<<std::dec<<"\t" <<str
-		<< "\t cpuLock  " << AccCache.cpuLock
-		<< "\t accLock  " << AccCache.accLock
-		<< "\t LRUvalid " << AccCache.LRU_valid<<std::endl;
-    }
-
-    assert( AccCache.cpuLock== 0 ) ;
-    assert( AccCache.accLock== 0 ) ;
-  }
-  std::cout << " Memory Manager::Audit() no locked table entries "<<std::endl;
-  assert(LruBytes1==LruBytes2);
-  assert(LruBytes1==DeviceLRUBytes);
-  std::cout << " Memory Manager::Audit() evictable bytes matches sum over table "<<std::endl;
-  assert(AccBytes==DeviceBytes);
-  std::cout << " Memory Manager::Audit() device bytes matches sum over table "<<std::endl;
-  assert(LruCnt == LRU.size());
-  std::cout << " Memory Manager::Audit() LRU entry count matches "<<std::endl;
-}
 
 void MemoryManager::PrintState(void* _CpuPtr)
 {
@@ -583,8 +519,8 @@ void MemoryManager::PrintState(void* _CpuPtr)
     if ( AccCache.state==EvictNext) str = std::string("EvictNext");
 
     std::cout << GridLogMessage << "CpuAddr\t\tAccAddr\t\tState\t\tcpuLock\taccLock\tLRU_valid "<<std::endl;
-    std::cout << GridLogMessage << "\tx"<<std::hex<<AccCache.CpuPtr<<std::dec
+    std::cout << GridLogMessage << "0x"<<std::hex<<AccCache.CpuPtr<<std::dec
-    << "\tx"<<std::hex<<AccCache.AccPtr<<std::dec<<"\t" <<str
+    << "\t0x"<<std::hex<<AccCache.AccPtr<<std::dec<<"\t" <<str
     << "\t" << AccCache.cpuLock
     << "\t" << AccCache.accLock
     << "\t" << AccCache.LRU_valid<<std::endl;

Grid/allocator/MemoryManagerShared.cc

@@ -12,10 +12,11 @@ uint64_t  MemoryManager::HostToDeviceBytes;
 uint64_t  MemoryManager::DeviceToHostBytes;
 uint64_t  MemoryManager::HostToDeviceXfer;
 uint64_t  MemoryManager::DeviceToHostXfer;
-uint64_t  MemoryManager::DeviceEvictions;
+uint64_t  MemoryManager::DeviceAccesses;
-uint64_t  MemoryManager::DeviceDestroy;
+uint64_t  MemoryManager::HostAccesses;
+uint64_t  MemoryManager::DeviceAccessBytes;
+uint64_t  MemoryManager::HostAccessBytes;
 
-void  MemoryManager::Audit(std::string s){};
 void  MemoryManager::ViewClose(void* AccPtr,ViewMode mode){};
 void *MemoryManager::ViewOpen(void* CpuPtr,size_t bytes,ViewMode mode,ViewAdvise hint){ return CpuPtr; };
 int   MemoryManager::isOpen   (void* CpuPtr) { return 0;}
@@ -24,7 +25,6 @@ void  MemoryManager::PrintState(void* CpuPtr)
 std::cout << GridLogMessage << "Host<->Device memory movement not currently managed by Grid." << std::endl;
 };
 void  MemoryManager::Print(void){};
-void  MemoryManager::PrintAll(void){};
 void  MemoryManager::NotifyDeletion(void *ptr){};
 
 NAMESPACE_END(Grid);

Grid/communicator/Communicator_mpi3.cc

@@ -392,9 +392,9 @@ double CartesianCommunicator::StencilSendToRecvFromBegin(std::vector<CommsReques
     acceleratorCopyDeviceToDeviceAsynch(xmit,shm,bytes);
   }
 
-  //  if ( CommunicatorPolicy == CommunicatorPolicySequential ) {
+  if ( CommunicatorPolicy == CommunicatorPolicySequential ) {
-  //    this->StencilSendToRecvFromComplete(list,dir);
+    this->StencilSendToRecvFromComplete(list,dir);
-  //  }
+  }
 
   return off_node_bytes;
 }

Grid/lattice/Lattice.h

@@ -46,4 +46,3 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/lattice/Lattice_unary.h>
 #include <Grid/lattice/Lattice_transfer.h>
 #include <Grid/lattice/Lattice_basis.h>
-#include <Grid/lattice/Lattice_crc.h>

Grid/lattice/Lattice_base.h

@@ -129,7 +129,7 @@ public:
     
     auto exprCopy = expr;
     ExpressionViewOpen(exprCopy);
-    auto me  = View(AcceleratorWrite);
+    auto me  = View(AcceleratorWriteDiscard);
     accelerator_for(ss,me.size(),vobj::Nsimd(),{
       auto tmp = eval(ss,exprCopy);
       coalescedWrite(me[ss],tmp);
@@ -152,7 +152,7 @@ public:
 
     auto exprCopy = expr;
     ExpressionViewOpen(exprCopy);
-    auto me  = View(AcceleratorWrite);
+    auto me  = View(AcceleratorWriteDiscard);
     accelerator_for(ss,me.size(),vobj::Nsimd(),{
       auto tmp = eval(ss,exprCopy);
       coalescedWrite(me[ss],tmp);
@@ -174,7 +174,7 @@ public:
     this->checkerboard=cb;
     auto exprCopy = expr;
     ExpressionViewOpen(exprCopy);
-    auto me  = View(AcceleratorWrite);
+    auto me  = View(AcceleratorWriteDiscard);
     accelerator_for(ss,me.size(),vobj::Nsimd(),{
       auto tmp = eval(ss,exprCopy);
       coalescedWrite(me[ss],tmp);
@@ -245,7 +245,7 @@ public:
   ///////////////////////////////////////////
   // user defined constructor
   ///////////////////////////////////////////
-  Lattice(GridBase *grid,ViewMode mode=AcceleratorWrite) { 
+  Lattice(GridBase *grid,ViewMode mode=AcceleratorWriteDiscard) { 
     this->_grid = grid;
     resize(this->_grid->oSites());
     assert((((uint64_t)&this->_odata[0])&0xF) ==0);
@@ -288,7 +288,7 @@ public:
     typename std::enable_if<!std::is_same<robj,vobj>::value,int>::type i=0;
     conformable(*this,r);
     this->checkerboard = r.Checkerboard();
-    auto me =   View(AcceleratorWrite);
+    auto me =   View(AcceleratorWriteDiscard);
     auto him= r.View(AcceleratorRead);
     accelerator_for(ss,me.size(),vobj::Nsimd(),{
       coalescedWrite(me[ss],him(ss));
@@ -303,7 +303,7 @@ public:
   inline Lattice<vobj> & operator = (const Lattice<vobj> & r){
     this->checkerboard = r.Checkerboard();
     conformable(*this,r);
-    auto me =   View(AcceleratorWrite);
+    auto me =   View(AcceleratorWriteDiscard);
     auto him= r.View(AcceleratorRead);
     accelerator_for(ss,me.size(),vobj::Nsimd(),{
       coalescedWrite(me[ss],him(ss));

Grid/lattice/Lattice_crc.h

@@ -1,55 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/lattice/Lattice_crc.h
-
-    Copyright (C) 2021
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#pragma once
-
-NAMESPACE_BEGIN(Grid);
-
-template<class vobj> void DumpSliceNorm(std::string s,Lattice<vobj> &f,int mu=-1)
-{
-  auto ff = localNorm2(f);
-  if ( mu==-1 ) mu = f.Grid()->Nd()-1;
-  typedef typename vobj::tensor_reduced normtype;
-  typedef typename normtype::scalar_object scalar;
-  std::vector<scalar> sff;
-  sliceSum(ff,sff,mu);
-  for(int t=0;t<sff.size();t++){
-    std::cout << s<<" "<<t<<" "<<sff[t]<<std::endl;
-  }
-}
-
-template<class vobj> uint32_t crc(Lattice<vobj> & buf)
-{
-  autoView( buf_v , buf, CpuRead);
-  return ::crc32(0L,(unsigned char *)&buf_v[0],(size_t)sizeof(vobj)*buf.oSites());
-}
-
-#define CRC(U) std::cout << "FingerPrint "<<__FILE__ <<" "<< __LINE__ <<" "<< #U <<" "<<crc(U)<<std::endl;
-
-NAMESPACE_END(Grid);
-
-

Grid/lattice/Lattice_reduction.h

@@ -28,9 +28,6 @@ Author: Christoph Lehner <christoph@lhnr.de>
 #if defined(GRID_CUDA)||defined(GRID_HIP)
 #include <Grid/lattice/Lattice_reduction_gpu.h>
 #endif
-#if defined(GRID_SYCL)
-#include <Grid/lattice/Lattice_reduction_sycl.h>
-#endif
 
 NAMESPACE_BEGIN(Grid);
 
@@ -130,7 +127,7 @@ inline Double max(const Double *arg, Integer osites)
 template<class vobj>
 inline typename vobj::scalar_object sum(const vobj *arg, Integer osites)
 {
-#if defined(GRID_CUDA)||defined(GRID_HIP)||defined(GRID_SYCL)
+#if defined(GRID_CUDA)||defined(GRID_HIP)
   return sum_gpu(arg,osites);
 #else
   return sum_cpu(arg,osites);
@@ -139,54 +136,29 @@ inline typename vobj::scalar_object sum(const vobj *arg, Integer osites)
 template<class vobj>
 inline typename vobj::scalar_objectD sumD(const vobj *arg, Integer osites)
 {
-#if defined(GRID_CUDA)||defined(GRID_HIP)||defined(GRID_SYCL)
+#if defined(GRID_CUDA)||defined(GRID_HIP)
   return sumD_gpu(arg,osites);
 #else
   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)||defined(GRID_SYCL)
-  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)
 {
-  Integer osites = arg.Grid()->oSites();
+#if defined(GRID_CUDA)||defined(GRID_HIP)
-#if defined(GRID_CUDA)||defined(GRID_HIP)||defined(GRID_SYCL)
-  typename vobj::scalar_object ssum;
   autoView( arg_v, arg, AcceleratorRead);
-  ssum= sum_gpu(&arg_v[0],osites);
+  Integer osites = arg.Grid()->oSites();
+  auto ssum= sum_gpu(&arg_v[0],osites);
 #else
   autoView(arg_v, arg, CpuRead);
+  Integer osites = arg.Grid()->oSites();
   auto ssum= sum_cpu(&arg_v[0],osites);
 #endif  
   arg.Grid()->GlobalSum(ssum);
   return ssum;
 }
 
-template<class vobj>
-inline typename vobj::scalar_object sum_large(const Lattice<vobj> &arg)
-{
-#if defined(GRID_CUDA)||defined(GRID_HIP)||defined(GRID_SYCL)
-  autoView( arg_v, arg, AcceleratorRead);
-  Integer osites = arg.Grid()->oSites();
-  auto ssum= sum_gpu_large(&arg_v[0],osites);
-#else
-  autoView(arg_v, arg, CpuRead);
-  Integer osites = arg.Grid()->oSites();
-  auto ssum= sum_cpu(&arg_v[0],osites);
-#endif
-  arg.Grid()->GlobalSum(ssum);
-  return ssum;
-}
-
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Deterministic Reduction operations
 ////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -238,10 +210,11 @@ inline ComplexD rankInnerProduct(const Lattice<vobj> &left,const Lattice<vobj> &
   typedef decltype(innerProductD(vobj(),vobj())) inner_t;
   Vector<inner_t> inner_tmp(sites);
   auto inner_tmp_v = &inner_tmp[0];
+    
   {
     autoView( left_v , left, AcceleratorRead);
     autoView( right_v,right, AcceleratorRead);
-    // This code could read coalesce
+
     // GPU - SIMT lane compliance...
     accelerator_for( ss, sites, 1,{
 	auto x_l = left_v[ss];

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);
+  getNumBlocksAndThreads(size, sizeof(sobj), numThreads, numBlocks);
-  assert(ok);
-
   Integer smemSize = numThreads * sizeof(sobj);
 
   Vector<sobj> buffer(numBlocks);
@@ -220,54 +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);
-
-  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;
-  
-  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
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -280,13 +230,6 @@ inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites)
   return result;
 }
 
-template <class vobj>
+
-inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osites)
-{
-  typedef typename vobj::scalar_object sobj;
-  sobj result;
-  result = sumD_gpu_large(lat,osites);
-  return result;
-}
 
 NAMESPACE_END(Grid);

Grid/lattice/Lattice_reduction_sycl.h

@@ -1,125 +0,0 @@
-NAMESPACE_BEGIN(Grid);
-
-/////////////////////////////////////////////////////////////////////////////////////////////////////////
-// Possibly promote to double and sum
-/////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_tensor(const vobj *lat, Integer osites) 
-{
-  typedef typename vobj::scalar_object sobj;
-  typedef typename vobj::scalar_objectD sobjD;
-  sobj *mysum =(sobj *) malloc_shared(sizeof(sobj),*theGridAccelerator);
-  sobj identity; zeroit(identity);
-  sobj ret ; 
-
-  Integer nsimd= vobj::Nsimd();
-  
-  theGridAccelerator->submit([&](cl::sycl::handler &cgh) {
-     auto Reduction = cl::sycl::reduction(mysum,identity,std::plus<>());
-     cgh.parallel_for(cl::sycl::range<1>{osites},
-		      Reduction,
-		      [=] (cl::sycl::id<1> item, auto &sum) {
-      auto osite   = item[0];
-      sum +=Reduce(lat[osite]);
-     });
-   });
-  theGridAccelerator->wait();
-  ret = mysum[0];
-  free(mysum,*theGridAccelerator);
-  sobjD dret; convertType(dret,ret);
-  return dret;
-}
-
-template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osites)
-{
-  return sumD_gpu_tensor(lat,osites);
-}
-template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites)
-{
-  return sumD_gpu_large(lat,osites);
-}
-
-template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
-{
-  return sumD_gpu_large(lat,osites);
-}
-
-/////////////////////////////////////////////////////////////////////////////////////////////////////////
-// Return as same precision as input performing reduction in double precision though
-/////////////////////////////////////////////////////////////////////////////////////////////////////////
-template <class vobj>
-inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites) 
-{
-  typedef typename vobj::scalar_object sobj;
-  sobj result;
-  result = sumD_gpu(lat,osites);
-  return result;
-}
-
-template <class vobj>
-inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osites)
-{
-  typedef typename vobj::scalar_object sobj;
-  sobj result;
-  result = sumD_gpu_large(lat,osites);
-  return result;
-}
-
-NAMESPACE_END(Grid);
-
-/*
-template<class Double> Double svm_reduce(Double *vec,uint64_t L)
-{
-  Double sumResult; zeroit(sumResult);
-  Double *d_sum =(Double *)cl::sycl::malloc_shared(sizeof(Double),*theGridAccelerator);
-  Double identity;  zeroit(identity);
-  theGridAccelerator->submit([&](cl::sycl::handler &cgh) {
-     auto Reduction = cl::sycl::reduction(d_sum,identity,std::plus<>());
-     cgh.parallel_for(cl::sycl::range<1>{L},
-		      Reduction,
-		      [=] (cl::sycl::id<1> index, auto &sum) {
-	 sum +=vec[index];
-     });
-   });
-  theGridAccelerator->wait();
-  Double ret = d_sum[0];
-  free(d_sum,*theGridAccelerator);
-  std::cout << " svm_reduce finished "<<L<<" sites sum = " << ret <<std::endl;
-  return ret;
-}
-
-template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_repack(const vobj *lat, Integer osites)
-{
-  typedef typename vobj::vector_type  vector;
-  typedef typename vobj::scalar_type  scalar;
-
-  typedef typename vobj::scalar_typeD scalarD;
-  typedef typename vobj::scalar_objectD sobjD;
-
-  sobjD ret;
-  scalarD *ret_p = (scalarD *)&ret;
-  
-  const int nsimd = vobj::Nsimd();
-  const int words = sizeof(vobj)/sizeof(vector);
-
-  Vector<scalar> buffer(osites*nsimd);
-  scalar *buf = &buffer[0];
-  vector *dat = (vector *)lat;
-
-  for(int w=0;w<words;w++) {
-
-    accelerator_for(ss,osites,nsimd,{
-	int lane = acceleratorSIMTlane(nsimd);
-	buf[ss*nsimd+lane] = dat[ss*words+w].getlane(lane);
-    });
-    //Precision change at this point is to late to gain precision
-    ret_p[w] = svm_reduce(buf,nsimd*osites);
-  }
-  return ret;
-}
-*/

Grid/log/Log.cc

@@ -65,40 +65,29 @@ GridLogger GridLogSolver (1, "Solver", GridLogColours, "NORMAL");
 GridLogger GridLogError  (1, "Error" , GridLogColours, "RED");
 GridLogger GridLogWarning(1, "Warning", GridLogColours, "YELLOW");
 GridLogger GridLogMessage(1, "Message", GridLogColours, "NORMAL");
-GridLogger GridLogMemory (1, "Memory", GridLogColours, "NORMAL");
-GridLogger GridLogTracing(1, "Tracing", GridLogColours, "NORMAL");
 GridLogger GridLogDebug  (1, "Debug", GridLogColours, "PURPLE");
 GridLogger GridLogPerformance(1, "Performance", GridLogColours, "GREEN");
-GridLogger GridLogDslash     (1, "Dslash", GridLogColours, "BLUE");
 GridLogger GridLogIterative  (1, "Iterative", GridLogColours, "BLUE");
 GridLogger GridLogIntegrator (1, "Integrator", GridLogColours, "BLUE");
-GridLogger GridLogHMC (1, "HMC", GridLogColours, "BLUE");
 
 void GridLogConfigure(std::vector<std::string> &logstreams) {
-  GridLogError.Active(1);
+  GridLogError.Active(0);
   GridLogWarning.Active(0);
   GridLogMessage.Active(1); // at least the messages should be always on
-  GridLogMemory.Active(0); 
-  GridLogTracing.Active(0); 
   GridLogIterative.Active(0);
   GridLogDebug.Active(0);
   GridLogPerformance.Active(0);
-  GridLogDslash.Active(0);
   GridLogIntegrator.Active(1);
   GridLogColours.Active(0);
-  GridLogHMC.Active(1);
 
   for (int i = 0; i < logstreams.size(); i++) {
-    if (logstreams[i] == std::string("Tracing"))     GridLogTracing.Active(1);
+    if (logstreams[i] == std::string("Error"))       GridLogError.Active(1);
-    if (logstreams[i] == std::string("Memory"))      GridLogMemory.Active(1);
     if (logstreams[i] == std::string("Warning"))     GridLogWarning.Active(1);
     if (logstreams[i] == std::string("NoMessage"))   GridLogMessage.Active(0);
     if (logstreams[i] == std::string("Iterative"))   GridLogIterative.Active(1);
     if (logstreams[i] == std::string("Debug"))       GridLogDebug.Active(1);
     if (logstreams[i] == std::string("Performance")) GridLogPerformance.Active(1);
-    if (logstreams[i] == std::string("Dslash"))      GridLogDslash.Active(1);
+    if (logstreams[i] == std::string("Integrator"))  GridLogIntegrator.Active(1);
-    if (logstreams[i] == std::string("NoIntegrator"))GridLogIntegrator.Active(0);
-    if (logstreams[i] == std::string("NoHMC"))       GridLogHMC.Active(0);
     if (logstreams[i] == std::string("Colours"))     GridLogColours.Active(1);
   }
 }

Grid/log/Log.h

@@ -138,8 +138,7 @@ public:
         stream << std::setw(log.topWidth);
       }
       stream << log.topName << log.background()<< " : ";
-      //      stream << log.colour() <<  std::left;
+      stream << log.colour() <<  std::left;
-      stream <<  std::left;
       if (log.chanWidth > 0)
       {
         stream << std::setw(log.chanWidth);
@@ -154,9 +153,9 @@ public:
 	stream << log.evidence()
 	       << now	       << log.background() << " : " ;
       }
-      //      stream << log.colour();
+      stream << log.colour();
-      stream <<  std::right;
       stream.flags(f);
+
       return stream;
     } else { 
       return devnull;
@@ -181,12 +180,8 @@ extern GridLogger GridLogWarning;
 extern GridLogger GridLogMessage;
 extern GridLogger GridLogDebug  ;
 extern GridLogger GridLogPerformance;
-extern GridLogger GridLogDslash;
 extern GridLogger GridLogIterative  ;
 extern GridLogger GridLogIntegrator  ;
-extern GridLogger GridLogHMC;
-extern GridLogger GridLogMemory;
-extern GridLogger GridLogTracing;
 extern Colours    GridLogColours;
 
 std::string demangle(const char* name) ;

Grid/parallelIO/IldgIO.h

@@ -31,7 +31,6 @@ directory
 #include <fstream>
 #include <iomanip>
 #include <iostream>
-#include <string>
 #include <map>
 
 #include <pwd.h>
@@ -655,8 +654,7 @@ class IldgWriter : public ScidacWriter {
     // Fill ILDG header data struct
     //////////////////////////////////////////////////////
     ildgFormat ildgfmt ;
-    const std::string stNC = std::to_string( Nc ) ;
+    ildgfmt.field     = std::string("su3gauge");
-    ildgfmt.field          = std::string("su"+stNC+"gauge");
 
     if ( format == std::string("IEEE32BIG") ) { 
       ildgfmt.precision = 32;
@@ -873,8 +871,7 @@ class IldgReader : public GridLimeReader {
     } else { 
 
       assert(found_ildgFormat);
-      const std::string stNC = std::to_string( Nc ) ;
+      assert ( ildgFormat_.field == std::string("su3gauge") );
-      assert ( ildgFormat_.field == std::string("su"+stNC+"gauge") );
 
       ///////////////////////////////////////////////////////////////////////////////////////
       // Populate our Grid metadata as best we can
@@ -882,7 +879,7 @@ class IldgReader : public GridLimeReader {
 
       std::ostringstream vers; vers << ildgFormat_.version;
       FieldMetaData_.hdr_version = vers.str();
-      FieldMetaData_.data_type = std::string("4D_SU"+stNC+"_GAUGE_"+stNC+"x"+stNC);
+      FieldMetaData_.data_type = std::string("4D_SU3_GAUGE_3X3");
 
       FieldMetaData_.nd=4;
       FieldMetaData_.dimension.resize(4);

Grid/parallelIO/MetaData.h

@@ -6,8 +6,8 @@
 
     Copyright (C) 2015
 
+
     Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-    Author: Jamie Hudspith <renwick.james.hudspth@gmail.com>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -182,8 +182,8 @@ class GaugeStatistics
 public:
   void operator()(Lattice<vLorentzColourMatrixD> & data,FieldMetaData &header)
   {
-    header.link_trace = WilsonLoops<Impl>::linkTrace(data);
+    header.link_trace=WilsonLoops<Impl>::linkTrace(data);
-    header.plaquette  = WilsonLoops<Impl>::avgPlaquette(data);
+    header.plaquette =WilsonLoops<Impl>::avgPlaquette(data);
   }
 };
 typedef GaugeStatistics<PeriodicGimplD> PeriodicGaugeStatistics;
@@ -203,24 +203,20 @@ template<> inline void PrepareMetaData<vLorentzColourMatrixD>(Lattice<vLorentzCo
 //////////////////////////////////////////////////////////////////////
 inline void reconstruct3(LorentzColourMatrix & cm)
 {
-  assert( Nc < 4 && Nc > 1 ) ;
+  const int x=0;
+  const int y=1;
+  const int z=2;
   for(int mu=0;mu<Nd;mu++){
-    #if Nc == 2
+    cm(mu)()(2,x) = adj(cm(mu)()(0,y)*cm(mu)()(1,z)-cm(mu)()(0,z)*cm(mu)()(1,y)); //x= yz-zy
-      cm(mu)()(1,0) = -adj(cm(mu)()(0,y)) ;
+    cm(mu)()(2,y) = adj(cm(mu)()(0,z)*cm(mu)()(1,x)-cm(mu)()(0,x)*cm(mu)()(1,z)); //y= zx-xz
-      cm(mu)()(1,1) =  adj(cm(mu)()(0,x)) ;
+    cm(mu)()(2,z) = adj(cm(mu)()(0,x)*cm(mu)()(1,y)-cm(mu)()(0,y)*cm(mu)()(1,x)); //z= xy-yx
-    #else
-      const int x=0 , y=1 , z=2 ; // a little disinenuous labelling
-      cm(mu)()(2,x) = adj(cm(mu)()(0,y)*cm(mu)()(1,z)-cm(mu)()(0,z)*cm(mu)()(1,y)); //x= yz-zy
-      cm(mu)()(2,y) = adj(cm(mu)()(0,z)*cm(mu)()(1,x)-cm(mu)()(0,x)*cm(mu)()(1,z)); //y= zx-xz
-      cm(mu)()(2,z) = adj(cm(mu)()(0,x)*cm(mu)()(1,y)-cm(mu)()(0,y)*cm(mu)()(1,x)); //z= xy-yx
-    #endif
   }
 }
 
 ////////////////////////////////////////////////////////////////////////////////
 // Some data types for intermediate storage
 ////////////////////////////////////////////////////////////////////////////////
-template<typename vtype> using iLorentzColour2x3 = iVector<iVector<iVector<vtype, Nc>, Nc-1>, Nd >;
+template<typename vtype> using iLorentzColour2x3 = iVector<iVector<iVector<vtype, Nc>, 2>, Nd >;
 
 typedef iLorentzColour2x3<Complex>  LorentzColour2x3;
 typedef iLorentzColour2x3<ComplexF> LorentzColour2x3F;
@@ -282,6 +278,7 @@ struct GaugeSimpleMunger{
 
 template <class fobj, class sobj>
 struct GaugeSimpleUnmunger {
+
   void operator()(sobj &in, fobj &out) {
     for (int mu = 0; mu < Nd; mu++) {
       for (int i = 0; i < Nc; i++) {
@@ -320,8 +317,8 @@ template<class fobj,class sobj>
 struct Gauge3x2munger{
   void operator() (fobj &in,sobj &out){
     for(int mu=0;mu<Nd;mu++){
-      for(int i=0;i<Nc-1;i++){
+      for(int i=0;i<2;i++){
-	for(int j=0;j<Nc;j++){
+	for(int j=0;j<3;j++){
 	  out(mu)()(i,j) = in(mu)(i)(j);
 	}}
     }
@@ -333,8 +330,8 @@ template<class fobj,class sobj>
 struct Gauge3x2unmunger{
   void operator() (sobj &in,fobj &out){
     for(int mu=0;mu<Nd;mu++){
-      for(int i=0;i<Nc-1;i++){
+      for(int i=0;i<2;i++){
-	for(int j=0;j<Nc;j++){
+	for(int j=0;j<3;j++){
 	  out(mu)(i)(j) = in(mu)()(i,j);
 	}}
     }

Grid/parallelIO/NerscIO.h

@@ -9,7 +9,6 @@
     Author: Matt Spraggs <matthew.spraggs@gmail.com>
     Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     Author: paboyle <paboyle@ph.ed.ac.uk>
-    Author: Jamie Hudspith <renwick.james.hudspth@gmail.com>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -31,8 +30,6 @@
 #ifndef GRID_NERSC_IO_H
 #define GRID_NERSC_IO_H
 
-#include <string>
-
 NAMESPACE_BEGIN(Grid);
 
 using namespace Grid;
@@ -148,17 +145,15 @@ public:
 
     std::string format(header.floating_point);
 
-    const int ieee32big = (format == std::string("IEEE32BIG"));
+    int ieee32big = (format == std::string("IEEE32BIG"));
-    const int ieee32    = (format == std::string("IEEE32"));
+    int ieee32    = (format == std::string("IEEE32"));
-    const int ieee64big = (format == std::string("IEEE64BIG"));
+    int ieee64big = (format == std::string("IEEE64BIG"));
-    const int ieee64    = (format == std::string("IEEE64") || \
+    int ieee64    = (format == std::string("IEEE64") || format == std::string("IEEE64LITTLE"));
-			   format == std::string("IEEE64LITTLE"));
 
     uint32_t nersc_csum,scidac_csuma,scidac_csumb;
     // depending on datatype, set up munger;
     // munger is a function of <floating point, Real, data_type>
-    const std::string stNC = std::to_string( Nc ) ;
+    if ( header.data_type == std::string("4D_SU3_GAUGE") ) {
-    if ( header.data_type == std::string("4D_SU"+stNC+"_GAUGE") ) {
       if ( ieee32 || ieee32big ) {
 	BinaryIO::readLatticeObject<vLorentzColourMatrixD, LorentzColour2x3F> 
 	  (Umu,file,Gauge3x2munger<LorentzColour2x3F,LorentzColourMatrix>(), offset,format,
@@ -169,7 +164,7 @@ public:
 	  (Umu,file,Gauge3x2munger<LorentzColour2x3D,LorentzColourMatrix>(),offset,format,
 	   nersc_csum,scidac_csuma,scidac_csumb);
       }
-    } else if ( header.data_type == std::string("4D_SU"+stNC+"_GAUGE_"+stNC+"x"+stNC) ) {
+    } else if ( header.data_type == std::string("4D_SU3_GAUGE_3x3") ) {
       if ( ieee32 || ieee32big ) {
 	BinaryIO::readLatticeObject<vLorentzColourMatrixD,LorentzColourMatrixF>
 	  (Umu,file,GaugeSimpleMunger<LorentzColourMatrixF,LorentzColourMatrix>(),offset,format,
@@ -214,29 +209,27 @@ public:
   template<class GaugeStats=PeriodicGaugeStatistics>
   static inline void writeConfiguration(Lattice<vLorentzColourMatrixD > &Umu,
 					std::string file, 
-					std::string ens_label = std::string("DWF"),
+					std::string ens_label = std::string("DWF"))
-					std::string ens_id = std::string("UKQCD"),
-					unsigned int sequence_number = 1)
   {
-    writeConfiguration(Umu,file,0,1,ens_label,ens_id,sequence_number);
+    writeConfiguration(Umu,file,0,1,ens_label);
   }
   template<class GaugeStats=PeriodicGaugeStatistics>
   static inline void writeConfiguration(Lattice<vLorentzColourMatrixD > &Umu,
 					std::string file, 
 					int two_row,
 					int bits32,
-					std::string ens_label = std::string("DWF"),
+					std::string ens_label = std::string("DWF"))
-					std::string ens_id = std::string("UKQCD"),
-					unsigned int sequence_number = 1)
   {
     typedef vLorentzColourMatrixD vobj;
     typedef typename vobj::scalar_object sobj;
 
     FieldMetaData header;
-    header.sequence_number = sequence_number;
+    ///////////////////////////////////////////
-    header.ensemble_id     = ens_id;
+    // Following should become arguments
+    ///////////////////////////////////////////
+    header.sequence_number = 1;
+    header.ensemble_id     = std::string("UKQCD");
     header.ensemble_label  = ens_label;
-    header.hdr_version     = "1.0" ;
 
     typedef LorentzColourMatrixD fobj3D;
     typedef LorentzColour2x3D    fobj2D;
@@ -250,14 +243,10 @@ public:
 
     uint64_t offset;
 
-    // Sod it -- always write NcxNc double
+    // Sod it -- always write 3x3 double
-    header.floating_point  = std::string("IEEE64BIG");
+    header.floating_point = std::string("IEEE64BIG");
-    const std::string stNC = std::to_string( Nc ) ;
+    header.data_type      = std::string("4D_SU3_GAUGE_3x3");
-    if( two_row ) {
+    GaugeSimpleUnmunger<fobj3D,sobj> munge;
-      header.data_type = std::string("4D_SU" + stNC + "_GAUGE" );
-    } else {
-      header.data_type = std::string("4D_SU" + stNC + "_GAUGE_" + stNC + "x" + stNC );
-    }
     if ( grid->IsBoss() ) { 
       truncate(file);
       offset = writeHeader(header,file);
@@ -265,15 +254,8 @@ public:
     grid->Broadcast(0,(void *)&offset,sizeof(offset));
 
     uint32_t nersc_csum,scidac_csuma,scidac_csumb;
-    if( two_row ) {
+    BinaryIO::writeLatticeObject<vobj,fobj3D>(Umu,file,munge,offset,header.floating_point,
-      Gauge3x2unmunger<fobj2D,sobj> munge;
+					      nersc_csum,scidac_csuma,scidac_csumb);
-      BinaryIO::writeLatticeObject<vobj,fobj2D>(Umu,file,munge,offset,header.floating_point,
-						nersc_csum,scidac_csuma,scidac_csumb);
-    } else {
-      GaugeSimpleUnmunger<fobj3D,sobj> munge;
-      BinaryIO::writeLatticeObject<vobj,fobj3D>(Umu,file,munge,offset,header.floating_point,
-						nersc_csum,scidac_csuma,scidac_csumb);
-    }
     header.checksum = nersc_csum;
     if ( grid->IsBoss() ) { 
       writeHeader(header,file);
@@ -305,7 +287,8 @@ public:
     header.plaquette=0.0;
     MachineCharacteristics(header);
 
-    uint64_t offset;
+	uint64_t offset;
+  
 #ifdef RNG_RANLUX
     header.floating_point = std::string("UINT64");
     header.data_type      = std::string("RANLUX48");
@@ -345,7 +328,7 @@ public:
 
     GridBase *grid = parallel.Grid();
 
-    uint64_t offset = readHeader(file,grid,header);
+	uint64_t offset = readHeader(file,grid,header);
 
     FieldMetaData clone(header);
 

Grid/perfmon/PerfCount.cc

@@ -27,12 +27,9 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 /*  END LEGAL */
 
 #include <Grid/GridCore.h>
-
-#include <Grid/perfmon/Timer.h>
 #include <Grid/perfmon/PerfCount.h>
-NAMESPACE_BEGIN(Grid);
 
-GridTimePoint theProgramStart = GridClock::now();
+NAMESPACE_BEGIN(Grid);
 
 #define CacheControl(L,O,R) ((PERF_COUNT_HW_CACHE_##L)|(PERF_COUNT_HW_CACHE_OP_##O<<8)| (PERF_COUNT_HW_CACHE_RESULT_##R<<16))
 #define RawConfig(A,B) (A<<8|B)

Grid/perfmon/PerfCount.h

@@ -30,12 +30,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #ifndef GRID_PERFCOUNT_H
 #define GRID_PERFCOUNT_H
 
-
-#ifndef __SSC_START
-#define __SSC_START
-#define __SSC_STOP
-#endif
-
 #include <sys/time.h>
 #include <ctime>
 #include <chrono>
@@ -78,9 +72,17 @@ static long perf_event_open(struct perf_event_attr *hw_event, pid_t pid,
 inline uint64_t cyclecount(void){ 
   return 0;
 }
+#define __SSC_MARK(mark) __asm__ __volatile__ ("movl %0, %%ebx; .byte 0x64, 0x67, 0x90 " ::"i"(mark):"%ebx")
+#define __SSC_STOP  __SSC_MARK(0x110)
+#define __SSC_START __SSC_MARK(0x111)
+
 
 #else
 
+#define __SSC_MARK(mark) 
+#define __SSC_STOP  
+#define __SSC_START 
+
 /*
  * cycle counters arch dependent
  */

Grid/perfmon/Timer.h

@@ -35,8 +35,17 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 NAMESPACE_BEGIN(Grid)
 
-//typedef  std::chrono::system_clock          GridClock;
+// Dress the output; use std::chrono
-typedef  std::chrono::high_resolution_clock   GridClock;
+// C++11 time facilities better?
+inline double usecond(void) {
+  struct timeval tv;
+#ifdef TIMERS_ON
+  gettimeofday(&tv,NULL);
+#endif
+  return 1.0*tv.tv_usec + 1.0e6*tv.tv_sec;
+}
+
+typedef  std::chrono::system_clock          GridClock;
 typedef  std::chrono::time_point<GridClock> GridTimePoint;
 
 typedef  std::chrono::seconds               GridSecs;
@@ -44,15 +53,6 @@ typedef  std::chrono::milliseconds          GridMillisecs;
 typedef  std::chrono::microseconds          GridUsecs;
 typedef  std::chrono::microseconds          GridTime;
 
-extern GridTimePoint theProgramStart;
-// Dress the output; use std::chrono
-// C++11 time facilities better?
-inline double usecond(void) {
-  auto usecs = std::chrono::duration_cast<GridUsecs>(GridClock::now()-theProgramStart); 
-  return 1.0*usecs.count();
-}
-
-
 inline std::ostream& operator<< (std::ostream & stream, const GridSecs & time)
 {
   stream << time.count()<<" s";

Grid/perfmon/Tracing.h

@@ -1,70 +0,0 @@
-#pragma once
-
-NAMESPACE_BEGIN(Grid);
-
-#ifdef GRID_TRACING_NVTX
-#include <nvToolsExt.h>
-class GridTracer {
-public:
-  GridTracer(const char* name) {
-    nvtxRangePushA(name);
-  }
-  ~GridTracer() {
-    nvtxRangePop();
-  }
-};
-inline void tracePush(const char *name) { nvtxRangePushA(name); }
-inline void tracePop(const char *name) { nvtxRangePop(); }
-inline int  traceStart(const char *name) {  }
-inline void traceStop(int ID) {  }
-#endif
-
-#ifdef GRID_TRACING_ROCTX
-#include <roctracer/roctx.h>
-class GridTracer {
- public:
-  GridTracer(const char* name) {
-    roctxRangePushA(name);
-    std::cout << "roctxRangePush "<<name<<std::endl;
-  }
-  ~GridTracer() {
-    roctxRangePop();
-    std::cout << "roctxRangePop "<<std::endl;
-  }
-};
-inline void tracePush(const char *name) { roctxRangePushA(name); }
-inline void tracePop(const char *name) { roctxRangePop(); }
-inline int  traceStart(const char *name) { roctxRangeStart(name); }
-inline void traceStop(int ID) { roctxRangeStop(ID); }
-#endif
-
-#ifdef GRID_TRACING_TIMER
-class GridTracer {
- public:
-  const char *name;
-  double elapsed;
-  GridTracer(const char* _name) {
-    name = _name;
-    elapsed=-usecond();
-  }
-  ~GridTracer() {
-    elapsed+=usecond();
-    std::cout << GridLogTracing << name << " took " <<elapsed<< " us" <<std::endl;
-  }
-};
-inline void tracePush(const char *name) {  }
-inline void tracePop(const char *name) {  }
-inline int  traceStart(const char *name) { return 0; }
-inline void traceStop(int ID) {  }
-#endif
-
-#ifdef GRID_TRACING_NONE
-#define GRID_TRACE(name) 
-inline void tracePush(const char *name) {  }
-inline void tracePop(const char *name) {  }
-inline int  traceStart(const char *name) { return 0;  }
-inline void traceStop(int ID) {  }
-#else
-#define GRID_TRACE(name) GridTracer uniq_name_using_macros##__COUNTER__(name);
-#endif
-NAMESPACE_END(Grid);

Grid/pugixml/pugixml.cc

@@ -16,12 +16,8 @@
 
 #ifdef __NVCC__
 #pragma push
-#ifdef __NVCC_DIAG_PRAGMA_SUPPORT__
-#pragma nv_diag_suppress declared_but_not_referenced // suppress "function was declared but never referenced warning"
-#else
 #pragma diag_suppress declared_but_not_referenced // suppress "function was declared but never referenced warning"
 #endif
-#endif
 
 #include "pugixml.h"
 

Grid/qcd/QCD.h

@@ -451,20 +451,9 @@ template<class vobj> void pokeLorentz(vobj &lhs,const decltype(peekIndex<Lorentz
 // Fermion <-> propagator assignements
 //////////////////////////////////////////////
 //template <class Prop, class Ferm>
-#define FAST_FERM_TO_PROP
 template <class Fimpl>
 void FermToProp(typename Fimpl::PropagatorField &p, const typename Fimpl::FermionField &f, const int s, const int c)
 {
-#ifdef FAST_FERM_TO_PROP
-  autoView(p_v,p,CpuWrite);
-  autoView(f_v,f,CpuRead);
-  thread_for(idx,p_v.oSites(),{
-      for(int ss = 0; ss < Ns; ++ss) {
-      for(int cc = 0; cc < Fimpl::Dimension; ++cc) {
-	p_v[idx]()(ss,s)(cc,c) = f_v[idx]()(ss)(cc); // Propagator sink index is LEFT, suitable for left mult by gauge link (e.g.)
-      }}
-    });
-#else
   for(int j = 0; j < Ns; ++j)
     {
       auto pjs = peekSpin(p, j, s);
@@ -476,23 +465,12 @@ void FermToProp(typename Fimpl::PropagatorField &p, const typename Fimpl::Fermio
 	}
       pokeSpin(p, pjs, j, s);
     }
-#endif
 }
     
 //template <class Prop, class Ferm>
 template <class Fimpl>
 void PropToFerm(typename Fimpl::FermionField &f, const typename Fimpl::PropagatorField &p, const int s, const int c)
 {
-#ifdef FAST_FERM_TO_PROP
-  autoView(p_v,p,CpuRead);
-  autoView(f_v,f,CpuWrite);
-  thread_for(idx,p_v.oSites(),{
-      for(int ss = 0; ss < Ns; ++ss) {
-      for(int cc = 0; cc < Fimpl::Dimension; ++cc) {
-	f_v[idx]()(ss)(cc) = p_v[idx]()(ss,s)(cc,c); // LEFT index is copied across for s,c right index
-      }}
-    });
-#else
   for(int j = 0; j < Ns; ++j)
     {
       auto pjs = peekSpin(p, j, s);
@@ -504,7 +482,6 @@ void PropToFerm(typename Fimpl::FermionField &f, const typename Fimpl::Propagato
 	}
       pokeSpin(f, fj, j);
     }
-#endif
 }
     
 //////////////////////////////////////////////

Grid/qcd/action/fermion/CayleyFermion5D.h

@@ -68,16 +68,9 @@ public:
   ///////////////////////////////////////////////////////////////
   // Support for MADWF tricks
   ///////////////////////////////////////////////////////////////
-  RealD Mass(void) { return (mass_plus + mass_minus) / 2.0; };
+  RealD Mass(void) { return mass; };
-  RealD MassPlus(void) { return mass_plus; };
-  RealD MassMinus(void) { return mass_minus; };
   void  SetMass(RealD _mass) { 
-    mass_plus=mass_minus=_mass; 
+    mass=_mass; 
-    SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c);  // Reset coeffs
-  } ;
-  void  SetMass(RealD _mass_plus, RealD _mass_minus) { 
-    mass_plus=_mass_plus;
-    mass_minus=_mass_minus;
     SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c);  // Reset coeffs
   } ;
   void  P(const FermionField &psi, FermionField &chi);
@@ -115,7 +108,7 @@ public:
   void   MeooeDag5D    (const FermionField &in, FermionField &out);
 
   //    protected:
-  RealD mass_plus, mass_minus;
+  RealD mass;
 
   // Save arguments to SetCoefficientsInternal
   Vector<Coeff_t> _gamma;

Grid/qcd/action/fermion/CloverHelpers.h

@@ -1,333 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./lib/qcd/action/fermion/WilsonCloverFermionImplementation.h
-
-    Copyright (C) 2017 - 2022
-
-    Author: paboyle <paboyle@ph.ed.ac.uk>
-    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
-    Author: Mattia Bruno <mattia.bruno@cern.ch>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-/*  END LEGAL */
-
-#pragma once
-
-#include <Grid/Grid.h>
-#include <Grid/qcd/spin/Dirac.h>
-#include <Grid/qcd/action/fermion/WilsonCloverHelpers.h>
-
-////////////////////////////////////////////
-// Standard Clover
-//   (4+m0) + csw * clover_term
-// Exp Clover
-//   (4+m0) * exp(csw/(4+m0) clover_term)
-//   = (4+m0) + csw * clover_term + ...
-////////////////////////////////////////////
-
-NAMESPACE_BEGIN(Grid);
-
-
-//////////////////////////////////
-// Generic Standard Clover
-//////////////////////////////////
-
-template<class Impl>
-class CloverHelpers: public WilsonCloverHelpers<Impl> {
-public:
-
-  INHERIT_IMPL_TYPES(Impl);
-  INHERIT_CLOVER_TYPES(Impl);
-
-  typedef WilsonCloverHelpers<Impl> Helpers;
-
-  static void Instantiate(CloverField& CloverTerm, CloverField& CloverTermInv, RealD csw_t, RealD diag_mass) {
-    GridBase *grid = CloverTerm.Grid();
-    CloverTerm += diag_mass;
-
-    int lvol = grid->lSites();
-    int DimRep = Impl::Dimension;
-    {
-      autoView(CTv,CloverTerm,CpuRead);
-      autoView(CTIv,CloverTermInv,CpuWrite);
-      thread_for(site, lvol, {
-        Coordinate lcoor;
-        grid->LocalIndexToLocalCoor(site, lcoor);
-        Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
-        Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
-        typename SiteClover::scalar_object Qx = Zero(), Qxinv = Zero();
-        peekLocalSite(Qx, CTv, lcoor);
-
-        for (int j = 0; j < Ns; j++)
-          for (int k = 0; k < Ns; k++)
-            for (int a = 0; a < DimRep; a++)
-              for (int b = 0; b < DimRep; b++){
-                auto zz =  Qx()(j, k)(a, b);
-                EigenCloverOp(a + j * DimRep, b + k * DimRep) = std::complex<double>(zz);
-              }
-
-        EigenInvCloverOp = EigenCloverOp.inverse();
-        for (int j = 0; j < Ns; j++)
-          for (int k = 0; k < Ns; k++)
-            for (int a = 0; a < DimRep; a++)
-              for (int b = 0; b < DimRep; b++)
-                Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep);
-               pokeLocalSite(Qxinv, CTIv, lcoor);
-      });
-    }
-  }
-
-  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu) {
-    return Helpers::Cmunu(U, lambda, mu, nu);
-  }
-
-};
-
-
-//////////////////////////////////
-// Generic Exp Clover
-//////////////////////////////////
-
-template<class Impl>
-class ExpCloverHelpers: public WilsonCloverHelpers<Impl> {
-public:
-
-  INHERIT_IMPL_TYPES(Impl);
-  INHERIT_CLOVER_TYPES(Impl);
-
-  template <typename vtype> using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
-  typedef WilsonCloverHelpers<Impl> Helpers;
-
-  // Can this be avoided?
-  static void IdentityTimesC(const CloverField& in, RealD c) {
-    int DimRep = Impl::Dimension;
-
-    autoView(in_v, in, AcceleratorWrite);
-
-    accelerator_for(ss, in.Grid()->oSites(), 1, {
-      for (int sa=0; sa<Ns; sa++)
-        for (int ca=0; ca<DimRep; ca++)
-          in_v[ss]()(sa,sa)(ca,ca) = c;
-    });
-  }
-
-  static int getNMAX(RealD prec, RealD R) {
-    /* compute stop condition for exponential */
-    int NMAX=1;
-    RealD cond=R*R/2.;
-
-    while (cond*std::exp(R)>prec) {
-      NMAX++;
-      cond*=R/(double)(NMAX+1);
-    }
-    return NMAX;
-  }
-
-  static int getNMAX(Lattice<iImplClover<vComplexD>> &t, RealD R) {return getNMAX(1e-12,R);}
-  static int getNMAX(Lattice<iImplClover<vComplexF>> &t, RealD R) {return getNMAX(1e-6,R);}
-
-  static void Instantiate(CloverField& Clover, CloverField& CloverInv, RealD csw_t, RealD diag_mass) {
-    GridBase* grid = Clover.Grid();
-    CloverField ExpClover(grid);
-
-    int NMAX = getNMAX(Clover, 3.*csw_t/diag_mass);
-
-    Clover *= (1.0/diag_mass);
-
-    // Taylor expansion, slow but generic
-    // Horner scheme: a0 + a1 x + a2 x^2 + .. = a0 + x (a1 + x(...))
-    // qN = cN
-    // qn = cn + qn+1 X
-    std::vector<RealD> cn(NMAX+1);
-    cn[0] = 1.0;
-    for (int i=1; i<=NMAX; i++)
-      cn[i] = cn[i-1] / RealD(i);
-
-    ExpClover = Zero();
-    IdentityTimesC(ExpClover, cn[NMAX]);
-    for (int i=NMAX-1; i>=0; i--)
-      ExpClover = ExpClover * Clover + cn[i];
-
-    // prepare inverse
-    CloverInv = (-1.0)*Clover;
-
-    Clover = ExpClover * diag_mass;
-
-    ExpClover = Zero();
-    IdentityTimesC(ExpClover, cn[NMAX]);
-    for (int i=NMAX-1; i>=0; i--)
-      ExpClover = ExpClover * CloverInv + cn[i];
-
-    CloverInv = ExpClover * (1.0/diag_mass);
-
-  }
-
-  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu) {
-    assert(0);
-    return lambda;
-  }
-
-};
-
-
-//////////////////////////////////
-// Compact Standard Clover
-//////////////////////////////////
-
-
-template<class Impl>
-class CompactCloverHelpers: public CompactWilsonCloverHelpers<Impl>,
-                            public WilsonCloverHelpers<Impl> {
-public:
-
-  INHERIT_IMPL_TYPES(Impl);
-  INHERIT_CLOVER_TYPES(Impl);
-  INHERIT_COMPACT_CLOVER_TYPES(Impl);
-
-  typedef WilsonCloverHelpers<Impl> Helpers;
-  typedef CompactWilsonCloverHelpers<Impl> CompactHelpers;
-
-  static void InstantiateClover(CloverField& Clover, CloverField& CloverInv, RealD csw_t, RealD diag_mass) {
-    Clover += diag_mass;
-  }
-
-  static void InvertClover(CloverField& InvClover,
-                            const CloverDiagonalField& diagonal,
-                            const CloverTriangleField& triangle,
-                            CloverDiagonalField&       diagonalInv,
-                            CloverTriangleField&       triangleInv,
-                            bool fixedBoundaries) {
-
-    CompactHelpers::Invert(diagonal, triangle, diagonalInv, triangleInv);
-  }
-
-  // TODO: implement Cmunu for better performances with compact layout, but don't do it
-  // here, but rather in WilsonCloverHelpers.h -> CompactWilsonCloverHelpers
-  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu) {
-    return Helpers::Cmunu(U, lambda, mu, nu);
-  }
-};
-
-//////////////////////////////////
-// Compact Exp Clover
-//////////////////////////////////
-
-template<class Impl>
-class CompactExpCloverHelpers: public CompactWilsonCloverHelpers<Impl> {
-public:
-
-  INHERIT_IMPL_TYPES(Impl);
-  INHERIT_CLOVER_TYPES(Impl);
-  INHERIT_COMPACT_CLOVER_TYPES(Impl);
-
-  template <typename vtype> using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
-  typedef CompactWilsonCloverHelpers<Impl> CompactHelpers;
-
-  // Can this be avoided?
-  static void IdentityTimesC(const CloverField& in, RealD c) {
-    int DimRep = Impl::Dimension;
-
-    autoView(in_v, in, AcceleratorWrite);
-
-    accelerator_for(ss, in.Grid()->oSites(), 1, {
-      for (int sa=0; sa<Ns; sa++)
-        for (int ca=0; ca<DimRep; ca++)
-          in_v[ss]()(sa,sa)(ca,ca) = c;
-    });
-  }
-
-  static int getNMAX(RealD prec, RealD R) {
-    /* compute stop condition for exponential */
-    int NMAX=1;
-    RealD cond=R*R/2.;
-
-    while (cond*std::exp(R)>prec) {
-      NMAX++;
-      cond*=R/(double)(NMAX+1);
-    }
-    return NMAX;
-  }
-
-  static int getNMAX(Lattice<iImplClover<vComplexD>> &t, RealD R) {return getNMAX(1e-12,R);}
-  static int getNMAX(Lattice<iImplClover<vComplexF>> &t, RealD R) {return getNMAX(1e-6,R);}
-
-  static void InstantiateClover(CloverField& Clover, CloverField& CloverInv, RealD csw_t, RealD diag_mass) {
-
-    GridBase* grid = Clover.Grid();
-    CloverField ExpClover(grid);
-
-    int NMAX = getNMAX(Clover, 3.*csw_t/diag_mass);
-
-    Clover *= (1.0/diag_mass);
-
-    // Taylor expansion, slow but generic
-    // Horner scheme: a0 + a1 x + a2 x^2 + .. = a0 + x (a1 + x(...))
-    // qN = cN
-    // qn = cn + qn+1 X
-    std::vector<RealD> cn(NMAX+1);
-    cn[0] = 1.0;
-    for (int i=1; i<=NMAX; i++)
-      cn[i] = cn[i-1] / RealD(i);
-
-    ExpClover = Zero();
-    IdentityTimesC(ExpClover, cn[NMAX]);
-    for (int i=NMAX-1; i>=0; i--)
-      ExpClover = ExpClover * Clover + cn[i];
-
-    // prepare inverse
-    CloverInv = (-1.0)*Clover;
-
-    Clover = ExpClover * diag_mass;
-
-    ExpClover = Zero();
-    IdentityTimesC(ExpClover, cn[NMAX]);
-    for (int i=NMAX-1; i>=0; i--)
-      ExpClover = ExpClover * CloverInv + cn[i];
-
-    CloverInv = ExpClover * (1.0/diag_mass);
-
-  }
-
-  static void InvertClover(CloverField& InvClover,
-                            const CloverDiagonalField& diagonal,
-                            const CloverTriangleField& triangle,
-                            CloverDiagonalField&       diagonalInv,
-                            CloverTriangleField&       triangleInv,
-                            bool fixedBoundaries) {
-
-    if (fixedBoundaries)
-    {
-      CompactHelpers::Invert(diagonal, triangle, diagonalInv, triangleInv);
-    }
-    else
-    {
-      CompactHelpers::ConvertLayout(InvClover, diagonalInv, triangleInv);
-    }
-  }
-
-  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu) {
-    assert(0);
-    return lambda;
-  }
-
-};
-
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/CompactWilsonCloverFermion.h

@@ -1,241 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./lib/qcd/action/fermion/CompactWilsonCloverFermion.h
-
-    Copyright (C) 2020 - 2022
-
-    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
-    Author: Nils Meyer <nils.meyer@ur.de>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-/*  END LEGAL */
-
-#pragma once
-
-#include <Grid/qcd/action/fermion/WilsonCloverTypes.h>
-#include <Grid/qcd/action/fermion/WilsonCloverHelpers.h>
-#include <Grid/qcd/action/fermion/CloverHelpers.h>
-
-NAMESPACE_BEGIN(Grid);
-
-// see Grid/qcd/action/fermion/WilsonCloverFermion.h for description
-//
-// Modifications done here:
-//
-// Original: clover term = 12x12 matrix per site
-//
-// But: Only two diagonal 6x6 hermitian blocks are non-zero (also true for original, verified by running)
-// Sufficient to store/transfer only the real parts of the diagonal and one triangular part
-// 2 * (6 + 15 * 2) = 72 real or 36 complex words to be stored/transfered
-//
-// Here: Above but diagonal as complex numbers, i.e., need to store/transfer
-// 2 * (6 * 2 + 15 * 2) = 84 real or 42 complex words
-//
-// Words per site and improvement compared to original (combined with the input and output spinors):
-//
-// - Original: 2*12 + 12*12 = 168 words -> 1.00 x less
-// - Minimal:  2*12 + 36    =  60 words -> 2.80 x less
-// - Here:     2*12 + 42    =  66 words -> 2.55 x less
-//
-// These improvements directly translate to wall-clock time
-//
-// Data layout:
-//
-// - diagonal and triangle part as separate lattice fields,
-//   this was faster than as 1 combined field on all tested machines
-// - diagonal: as expected
-// - triangle: store upper right triangle in row major order
-// - graphical:
-//        0  1  2  3  4
-//           5  6  7  8
-//              9 10 11 = upper right triangle indices
-//                12 13
-//                   14
-//     0
-//        1
-//           2
-//              3       = diagonal indices
-//                 4
-//                    5
-//     0
-//     1  5
-//     2  6  9          = lower left triangle indices
-//     3  7 10 12
-//     4  8 11 13 14
-//
-// Impact on total memory consumption:
-// - Original: (2 * 1 + 8 * 1/2) 12x12 matrices = 6 12x12 matrices = 864 complex words per site
-// - Here:     (2 * 1 + 4 * 1/2) diagonal parts = 4 diagonal parts =  24 complex words per site
-//           + (2 * 1 + 4 * 1/2) triangle parts = 4 triangle parts =  60 complex words per site
-//                                                                 =  84 complex words per site
-
-template<class Impl, class CloverHelpers>
-class CompactWilsonCloverFermion : public WilsonFermion<Impl>,
-                                   public WilsonCloverHelpers<Impl>,
-                                   public CompactWilsonCloverHelpers<Impl> {
-  /////////////////////////////////////////////
-  // Sizes
-  /////////////////////////////////////////////
-
-public:
-
-  INHERIT_COMPACT_CLOVER_SIZES(Impl);
-
-  /////////////////////////////////////////////
-  // Type definitions
-  /////////////////////////////////////////////
-
-public:
-
-  INHERIT_IMPL_TYPES(Impl);
-  INHERIT_CLOVER_TYPES(Impl);
-  INHERIT_COMPACT_CLOVER_TYPES(Impl);
-
-  typedef WilsonFermion<Impl>              WilsonBase;
-  typedef WilsonCloverHelpers<Impl>        Helpers;
-  typedef CompactWilsonCloverHelpers<Impl> CompactHelpers;
-
-  /////////////////////////////////////////////
-  // Constructors
-  /////////////////////////////////////////////
-
-public:
-
-  CompactWilsonCloverFermion(GaugeField& _Umu,
-			    GridCartesian& Fgrid,
-			    GridRedBlackCartesian& Hgrid,
-			    const RealD _mass,
-			    const RealD _csw_r = 0.0,
-			    const RealD _csw_t = 0.0,
-			    const RealD _cF = 1.0,
-			    const WilsonAnisotropyCoefficients& clover_anisotropy = WilsonAnisotropyCoefficients(),
-			    const ImplParams& impl_p = ImplParams());
-
-  /////////////////////////////////////////////
-  // Member functions (implementing interface)
-  /////////////////////////////////////////////
-
-public:
-
-  virtual void Instantiatable() {};
-  int          ConstEE()     override { return 0; };
-  int          isTrivialEE() override { return 0; };
-
-  void Dhop(const FermionField& in, FermionField& out, int dag) override;
-
-  void DhopOE(const FermionField& in, FermionField& out, int dag) override;
-
-  void DhopEO(const FermionField& in, FermionField& out, int dag) override;
-
-  void DhopDir(const FermionField& in, FermionField& out, int dir, int disp) override;
-
-  void DhopDirAll(const FermionField& in, std::vector<FermionField>& out) /* override */;
-
-  void M(const FermionField& in, FermionField& out) override;
-
-  void Mdag(const FermionField& in, FermionField& out) override;
-
-  void Meooe(const FermionField& in, FermionField& out) override;
-
-  void MeooeDag(const FermionField& in, FermionField& out) override;
-
-  void Mooee(const FermionField& in, FermionField& out) override;
-
-  void MooeeDag(const FermionField& in, FermionField& out) override;
-
-  void MooeeInv(const FermionField& in, FermionField& out) override;
-
-  void MooeeInvDag(const FermionField& in, FermionField& out) override;
-
-  void Mdir(const FermionField& in, FermionField& out, int dir, int disp) override;
-
-  void MdirAll(const FermionField& in, std::vector<FermionField>& out) override;
-
-  void MDeriv(GaugeField& force, const FermionField& X, const FermionField& Y, int dag) override;
-
-  void MooDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) override;
-
-  void MeeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) override;
-
-  /////////////////////////////////////////////
-  // Member functions (internals)
-  /////////////////////////////////////////////
-
-  void MooeeInternal(const FermionField&        in,
-                     FermionField&              out,
-                     const CloverDiagonalField& diagonal,
-                     const CloverTriangleField& triangle);
-
-  /////////////////////////////////////////////
-  // Helpers
-  /////////////////////////////////////////////
-
-  void ImportGauge(const GaugeField& _Umu) override;
-
-  /////////////////////////////////////////////
-  // Helpers
-  /////////////////////////////////////////////
-
-private:
-
-  template<class Field>
-  const MaskField* getCorrectMaskField(const Field &in) const {
-    if(in.Grid()->_isCheckerBoarded) {
-      if(in.Checkerboard() == Odd) {
-        return &this->BoundaryMaskOdd;
-      } else {
-        return &this->BoundaryMaskEven;
-      }
-    } else {
-      return &this->BoundaryMask;
-    }
-  }
-
-  template<class Field>
-  void ApplyBoundaryMask(Field& f) {
-    const MaskField* m = getCorrectMaskField(f); assert(m != nullptr);
-    assert(m != nullptr);
-    CompactHelpers::ApplyBoundaryMask(f, *m);
-  }
-
-  /////////////////////////////////////////////
-  // Member Data
-  /////////////////////////////////////////////
-
-public:
-
-  RealD csw_r;
-  RealD csw_t;
-  RealD cF;
-
-  bool fixedBoundaries;
-
-  CloverDiagonalField Diagonal,    DiagonalEven,    DiagonalOdd;
-  CloverDiagonalField DiagonalInv, DiagonalInvEven, DiagonalInvOdd;
-
-  CloverTriangleField Triangle,    TriangleEven,    TriangleOdd;
-  CloverTriangleField TriangleInv, TriangleInvEven, TriangleInvOdd;
-
-  FermionField Tmp;
-
-  MaskField BoundaryMask, BoundaryMaskEven, BoundaryMaskOdd;
-};
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/Fermion.h

@@ -53,7 +53,6 @@ NAMESPACE_CHECK(Wilson);
 #include <Grid/qcd/action/fermion/WilsonTMFermion.h>       // 4d wilson like
 NAMESPACE_CHECK(WilsonTM);
 #include <Grid/qcd/action/fermion/WilsonCloverFermion.h> // 4d wilson clover fermions
-#include <Grid/qcd/action/fermion/CompactWilsonCloverFermion.h> // 4d compact wilson clover fermions
 NAMESPACE_CHECK(WilsonClover);
 #include <Grid/qcd/action/fermion/WilsonFermion5D.h>     // 5d base used by all 5d overlap types
 NAMESPACE_CHECK(Wilson5D);
@@ -138,52 +137,21 @@ typedef WilsonTMFermion<WilsonImplF> WilsonTMFermionF;
 typedef WilsonTMFermion<WilsonImplD> WilsonTMFermionD;
 
 // Clover fermions
-template <typename WImpl> using WilsonClover = WilsonCloverFermion<WImpl, CloverHelpers<WImpl>>;
+typedef WilsonCloverFermion<WilsonImplR> WilsonCloverFermionR;
-template <typename WImpl> using WilsonExpClover = WilsonCloverFermion<WImpl, ExpCloverHelpers<WImpl>>;
+typedef WilsonCloverFermion<WilsonImplF> WilsonCloverFermionF;
+typedef WilsonCloverFermion<WilsonImplD> WilsonCloverFermionD;
 
-typedef WilsonClover<WilsonImplR> WilsonCloverFermionR;
+typedef WilsonCloverFermion<WilsonAdjImplR> WilsonCloverAdjFermionR;
-typedef WilsonClover<WilsonImplF> WilsonCloverFermionF;
+typedef WilsonCloverFermion<WilsonAdjImplF> WilsonCloverAdjFermionF;
-typedef WilsonClover<WilsonImplD> WilsonCloverFermionD;
+typedef WilsonCloverFermion<WilsonAdjImplD> WilsonCloverAdjFermionD;
 
-typedef WilsonExpClover<WilsonImplR> WilsonExpCloverFermionR;
+typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplR> WilsonCloverTwoIndexSymmetricFermionR;
-typedef WilsonExpClover<WilsonImplF> WilsonExpCloverFermionF;
+typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplF> WilsonCloverTwoIndexSymmetricFermionF;
-typedef WilsonExpClover<WilsonImplD> WilsonExpCloverFermionD;
+typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplD> WilsonCloverTwoIndexSymmetricFermionD;
 
-typedef WilsonClover<WilsonAdjImplR> WilsonCloverAdjFermionR;
+typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplR> WilsonCloverTwoIndexAntiSymmetricFermionR;
-typedef WilsonClover<WilsonAdjImplF> WilsonCloverAdjFermionF;
+typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplF> WilsonCloverTwoIndexAntiSymmetricFermionF;
-typedef WilsonClover<WilsonAdjImplD> WilsonCloverAdjFermionD;
+typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplD> WilsonCloverTwoIndexAntiSymmetricFermionD;
-
-typedef WilsonClover<WilsonTwoIndexSymmetricImplR> WilsonCloverTwoIndexSymmetricFermionR;
-typedef WilsonClover<WilsonTwoIndexSymmetricImplF> WilsonCloverTwoIndexSymmetricFermionF;
-typedef WilsonClover<WilsonTwoIndexSymmetricImplD> WilsonCloverTwoIndexSymmetricFermionD;
-
-typedef WilsonClover<WilsonTwoIndexAntiSymmetricImplR> WilsonCloverTwoIndexAntiSymmetricFermionR;
-typedef WilsonClover<WilsonTwoIndexAntiSymmetricImplF> WilsonCloverTwoIndexAntiSymmetricFermionF;
-typedef WilsonClover<WilsonTwoIndexAntiSymmetricImplD> WilsonCloverTwoIndexAntiSymmetricFermionD;
-
-// Compact Clover fermions
-template <typename WImpl> using CompactWilsonClover = CompactWilsonCloverFermion<WImpl, CompactCloverHelpers<WImpl>>;
-template <typename WImpl> using CompactWilsonExpClover = CompactWilsonCloverFermion<WImpl, CompactExpCloverHelpers<WImpl>>;
-
-typedef CompactWilsonClover<WilsonImplR> CompactWilsonCloverFermionR;
-typedef CompactWilsonClover<WilsonImplF> CompactWilsonCloverFermionF;
-typedef CompactWilsonClover<WilsonImplD> CompactWilsonCloverFermionD;
-
-typedef CompactWilsonExpClover<WilsonImplR> CompactWilsonExpCloverFermionR;
-typedef CompactWilsonExpClover<WilsonImplF> CompactWilsonExpCloverFermionF;
-typedef CompactWilsonExpClover<WilsonImplD> CompactWilsonExpCloverFermionD;
-
-typedef CompactWilsonClover<WilsonAdjImplR> CompactWilsonCloverAdjFermionR;
-typedef CompactWilsonClover<WilsonAdjImplF> CompactWilsonCloverAdjFermionF;
-typedef CompactWilsonClover<WilsonAdjImplD> CompactWilsonCloverAdjFermionD;
-
-typedef CompactWilsonClover<WilsonTwoIndexSymmetricImplR> CompactWilsonCloverTwoIndexSymmetricFermionR;
-typedef CompactWilsonClover<WilsonTwoIndexSymmetricImplF> CompactWilsonCloverTwoIndexSymmetricFermionF;
-typedef CompactWilsonClover<WilsonTwoIndexSymmetricImplD> CompactWilsonCloverTwoIndexSymmetricFermionD;
-
-typedef CompactWilsonClover<WilsonTwoIndexAntiSymmetricImplR> CompactWilsonCloverTwoIndexAntiSymmetricFermionR;
-typedef CompactWilsonClover<WilsonTwoIndexAntiSymmetricImplF> CompactWilsonCloverTwoIndexAntiSymmetricFermionF;
-typedef CompactWilsonClover<WilsonTwoIndexAntiSymmetricImplD> CompactWilsonCloverTwoIndexAntiSymmetricFermionD;
 
 // Domain Wall fermions
 typedef DomainWallFermion<WilsonImplR> DomainWallFermionR;

Grid/qcd/action/fermion/WilsonCloverFermion.h

@@ -4,11 +4,10 @@
 
     Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.h
 
-    Copyright (C) 2017 - 2022
+    Copyright (C) 2017
 
     Author: Guido Cossu <guido.cossu@ed.ac.uk>
     Author: David Preti <>
-    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -30,9 +29,7 @@
 
 #pragma once
 
-#include <Grid/qcd/action/fermion/WilsonCloverTypes.h>
+#include <Grid/Grid.h>
-#include <Grid/qcd/action/fermion/WilsonCloverHelpers.h>
-#include <Grid/qcd/action/fermion/CloverHelpers.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -52,16 +49,19 @@ NAMESPACE_BEGIN(Grid);
 // csw_r = csw_t to recover the isotropic version
 //////////////////////////////////////////////////////////////////
 
-template<class Impl, class CloverHelpers>
+template <class Impl>
-class WilsonCloverFermion : public WilsonFermion<Impl>,
+class WilsonCloverFermion : public WilsonFermion<Impl>
-                            public WilsonCloverHelpers<Impl>
 {
 public:
+  // Types definitions
   INHERIT_IMPL_TYPES(Impl);
-  INHERIT_CLOVER_TYPES(Impl);
+  template <typename vtype>
+  using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
+  typedef iImplClover<Simd> SiteCloverType;
+  typedef Lattice<SiteCloverType> CloverFieldType;
 
-  typedef WilsonFermion<Impl>       WilsonBase;
+public:
-  typedef WilsonCloverHelpers<Impl> Helpers;
+  typedef WilsonFermion<Impl> WilsonBase;
 
   virtual int    ConstEE(void)     { return 0; };
   virtual void Instantiatable(void){};
@@ -72,7 +72,42 @@ public:
                       const RealD _csw_r = 0.0,
                       const RealD _csw_t = 0.0,
                       const WilsonAnisotropyCoefficients &clover_anisotropy = WilsonAnisotropyCoefficients(),
-                      const ImplParams &impl_p = ImplParams());
+                      const ImplParams &impl_p = ImplParams()) : WilsonFermion<Impl>(_Umu,
+                                                                                     Fgrid,
+                                                                                     Hgrid,
+                                                                                     _mass, impl_p, clover_anisotropy),
+                                                                 CloverTerm(&Fgrid),
+                                                                 CloverTermInv(&Fgrid),
+                                                                 CloverTermEven(&Hgrid),
+                                                                 CloverTermOdd(&Hgrid),
+                                                                 CloverTermInvEven(&Hgrid),
+                                                                 CloverTermInvOdd(&Hgrid),
+                                                                 CloverTermDagEven(&Hgrid),
+                                                                 CloverTermDagOdd(&Hgrid),
+                                                                 CloverTermInvDagEven(&Hgrid),
+                                                                 CloverTermInvDagOdd(&Hgrid)
+  {
+    assert(Nd == 4); // require 4 dimensions
+
+    if (clover_anisotropy.isAnisotropic)
+    {
+      csw_r = _csw_r * 0.5 / clover_anisotropy.xi_0;
+      diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0);
+    }
+    else
+    {
+      csw_r = _csw_r * 0.5;
+      diag_mass = 4.0 + _mass;
+    }
+    csw_t = _csw_t * 0.5;
+
+    if (csw_r == 0)
+      std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_r = 0" << std::endl;
+    if (csw_t == 0)
+      std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_t = 0" << std::endl;
+
+    ImportGauge(_Umu);
+  }
 
   virtual void M(const FermionField &in, FermionField &out);
   virtual void Mdag(const FermionField &in, FermionField &out);
@@ -89,21 +124,250 @@ public:
   void ImportGauge(const GaugeField &_Umu);
 
   // Derivative parts unpreconditioned pseudofermions
-  void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag);
+  void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag)
+  {
+    conformable(X.Grid(), Y.Grid());
+    conformable(X.Grid(), force.Grid());
+    GaugeLinkField force_mu(force.Grid()), lambda(force.Grid());
+    GaugeField clover_force(force.Grid());
+    PropagatorField Lambda(force.Grid());
 
-public:
+    // Guido: Here we are hitting some performance issues:
+    // need to extract the components of the DoubledGaugeField
+    // for each call
+    // Possible solution
+    // Create a vector object to store them? (cons: wasting space)
+    std::vector<GaugeLinkField> U(Nd, this->Umu.Grid());
+
+    Impl::extractLinkField(U, this->Umu);
+
+    force = Zero();
+    // Derivative of the Wilson hopping term
+    this->DhopDeriv(force, X, Y, dag);
+
+    ///////////////////////////////////////////////////////////
+    // Clover term derivative
+    ///////////////////////////////////////////////////////////
+    Impl::outerProductImpl(Lambda, X, Y);
+    //std::cout << "Lambda:" << Lambda << std::endl;
+
+    Gamma::Algebra sigma[] = {
+        Gamma::Algebra::SigmaXY,
+        Gamma::Algebra::SigmaXZ,
+        Gamma::Algebra::SigmaXT,
+        Gamma::Algebra::MinusSigmaXY,
+        Gamma::Algebra::SigmaYZ,
+        Gamma::Algebra::SigmaYT,
+        Gamma::Algebra::MinusSigmaXZ,
+        Gamma::Algebra::MinusSigmaYZ,
+        Gamma::Algebra::SigmaZT,
+        Gamma::Algebra::MinusSigmaXT,
+        Gamma::Algebra::MinusSigmaYT,
+        Gamma::Algebra::MinusSigmaZT};
+
+    /*
+      sigma_{\mu \nu}=
+      | 0         sigma[0]  sigma[1]  sigma[2] |
+      | sigma[3]    0       sigma[4]  sigma[5] |
+      | sigma[6]  sigma[7]     0      sigma[8] |
+      | sigma[9]  sigma[10] sigma[11]   0      |
+    */
+
+    int count = 0;
+    clover_force = Zero();
+    for (int mu = 0; mu < 4; mu++)
+    {
+      force_mu = Zero();
+      for (int nu = 0; nu < 4; nu++)
+      {
+        if (mu == nu)
+        continue;
+        
+        RealD factor;
+        if (nu == 4 || mu == 4)
+        {
+          factor = 2.0 * csw_t;
+        }
+        else
+        {
+          factor = 2.0 * csw_r;
+        }
+        PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
+        Impl::TraceSpinImpl(lambda, Slambda);                   // traceSpin ok
+        force_mu -= factor*Cmunu(U, lambda, mu, nu);                   // checked
+        count++;
+      }
+
+      pokeLorentz(clover_force, U[mu] * force_mu, mu);
+    }
+    //clover_force *= csw;
+    force += clover_force;
+  }
+
+  // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis
+  GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu)
+  {
+    conformable(lambda.Grid(), U[0].Grid());
+    GaugeLinkField out(lambda.Grid()), tmp(lambda.Grid());
+    // insertion in upper staple
+    // please check redundancy of shift operations
+
+    // C1+
+    tmp = lambda * U[nu];
+    out = Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
+
+    // C2+
+    tmp = U[mu] * Impl::ShiftStaple(adj(lambda), mu);
+    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
+
+    // C3+
+    tmp = U[nu] * Impl::ShiftStaple(adj(lambda), nu);
+    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu);
+
+    // C4+
+    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * lambda;
+
+    // insertion in lower staple
+    // C1-
+    out -= Impl::ShiftStaple(lambda, mu) * Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
+
+    // C2-
+    tmp = adj(lambda) * U[nu];
+    out -= Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
+
+    // C3-
+    tmp = lambda * U[nu];
+    out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu);
+
+    // C4-
+    out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * lambda;
+
+    return out;
+  }
+
+protected:
   // here fixing the 4 dimensions, make it more general?
 
   RealD csw_r;                                               // Clover coefficient - spatial
   RealD csw_t;                                               // Clover coefficient - temporal
   RealD diag_mass;                                           // Mass term
-  CloverField CloverTerm, CloverTermInv;                     // Clover term
+  CloverFieldType CloverTerm, CloverTermInv;                 // Clover term
-  CloverField CloverTermEven, CloverTermOdd;                 // Clover term EO
+  CloverFieldType CloverTermEven, CloverTermOdd;             // Clover term EO
-  CloverField CloverTermInvEven, CloverTermInvOdd;           // Clover term Inv EO
+  CloverFieldType CloverTermInvEven, CloverTermInvOdd;       // Clover term Inv EO
-  CloverField CloverTermDagEven, CloverTermDagOdd;           // Clover term Dag EO
+  CloverFieldType CloverTermDagEven, CloverTermDagOdd;       // Clover term Dag EO
-  CloverField CloverTermInvDagEven, CloverTermInvDagOdd;     // Clover term Inv Dag EO
+  CloverFieldType CloverTermInvDagEven, CloverTermInvDagOdd; // Clover term Inv Dag EO
-};
 
+ public:
+  // eventually these can be compressed into 6x6 blocks instead of the 12x12
+  // using the DeGrand-Rossi basis for the gamma matrices
+  CloverFieldType fillCloverYZ(const GaugeLinkField &F)
+  {
+    CloverFieldType T(F.Grid());
+    T = Zero();
+    autoView(T_v,T,AcceleratorWrite);
+    autoView(F_v,F,AcceleratorRead);
+    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 1) = timesMinusI(F_v[i]()());
+      T_v[i]()(1, 0) = timesMinusI(F_v[i]()());
+      T_v[i]()(2, 3) = timesMinusI(F_v[i]()());
+      T_v[i]()(3, 2) = timesMinusI(F_v[i]()());
+    });
+
+    return T;
+  }
+
+  CloverFieldType fillCloverXZ(const GaugeLinkField &F)
+  {
+    CloverFieldType T(F.Grid());
+    T = Zero();
+    
+    autoView(T_v, T,AcceleratorWrite);
+    autoView(F_v, F,AcceleratorRead);
+    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 1) = -F_v[i]()();
+      T_v[i]()(1, 0) = F_v[i]()();
+      T_v[i]()(2, 3) = -F_v[i]()();
+      T_v[i]()(3, 2) = F_v[i]()();
+    });
+
+    return T;
+  }
+
+  CloverFieldType fillCloverXY(const GaugeLinkField &F)
+  {
+    CloverFieldType T(F.Grid());
+    T = Zero();
+
+    autoView(T_v,T,AcceleratorWrite);
+    autoView(F_v,F,AcceleratorRead);
+    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 0) = timesMinusI(F_v[i]()());
+      T_v[i]()(1, 1) = timesI(F_v[i]()());
+      T_v[i]()(2, 2) = timesMinusI(F_v[i]()());
+      T_v[i]()(3, 3) = timesI(F_v[i]()());
+    });
+
+    return T;
+  }
+
+  CloverFieldType fillCloverXT(const GaugeLinkField &F)
+  {
+    CloverFieldType T(F.Grid());
+    T = Zero();
+
+    autoView( T_v , T, AcceleratorWrite);
+    autoView( F_v , F, AcceleratorRead);
+    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 1) = timesI(F_v[i]()());
+      T_v[i]()(1, 0) = timesI(F_v[i]()());
+      T_v[i]()(2, 3) = timesMinusI(F_v[i]()());
+      T_v[i]()(3, 2) = timesMinusI(F_v[i]()());
+    });
+
+    return T;
+  }
+
+  CloverFieldType fillCloverYT(const GaugeLinkField &F)
+  {
+    CloverFieldType T(F.Grid());
+    T = Zero();
+    
+    autoView( T_v ,T,AcceleratorWrite);
+    autoView( F_v ,F,AcceleratorRead);
+    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 1) = -(F_v[i]()());
+      T_v[i]()(1, 0) = (F_v[i]()());
+      T_v[i]()(2, 3) = (F_v[i]()());
+      T_v[i]()(3, 2) = -(F_v[i]()());
+    });
+
+    return T;
+  }
+
+  CloverFieldType fillCloverZT(const GaugeLinkField &F)
+  {
+    CloverFieldType T(F.Grid());
+
+    T = Zero();
+
+    autoView( T_v , T,AcceleratorWrite);
+    autoView( F_v , F,AcceleratorRead);
+    accelerator_for(i, CloverTerm.Grid()->oSites(),1,
+    {
+      T_v[i]()(0, 0) = timesI(F_v[i]()());
+      T_v[i]()(1, 1) = timesMinusI(F_v[i]()());
+      T_v[i]()(2, 2) = timesMinusI(F_v[i]()());
+      T_v[i]()(3, 3) = timesI(F_v[i]()());
+    });
+
+    return T;
+  }
+};
 NAMESPACE_END(Grid);
 
 

Grid/qcd/action/fermion/WilsonCloverHelpers.h

@@ -1,763 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./lib/qcd/action/fermion/WilsonCloverHelpers.h
-
-    Copyright (C) 2021 - 2022
-
-    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-
-#pragma once
-
-// Helper routines that implement common clover functionality
-
-NAMESPACE_BEGIN(Grid);
-
-template<class Impl> class WilsonCloverHelpers {
-public:
-
-  INHERIT_IMPL_TYPES(Impl);
-  INHERIT_CLOVER_TYPES(Impl);
-
-  // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis
-  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu)
-  {
-    conformable(lambda.Grid(), U[0].Grid());
-    GaugeLinkField out(lambda.Grid()), tmp(lambda.Grid());
-    // insertion in upper staple
-    // please check redundancy of shift operations
-
-    // C1+
-    tmp = lambda * U[nu];
-    out = Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
-
-    // C2+
-    tmp = U[mu] * Impl::ShiftStaple(adj(lambda), mu);
-    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
-
-    // C3+
-    tmp = U[nu] * Impl::ShiftStaple(adj(lambda), nu);
-    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu);
-
-    // C4+
-    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * lambda;
-
-    // insertion in lower staple
-    // C1-
-    out -= Impl::ShiftStaple(lambda, mu) * Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
-
-    // C2-
-    tmp = adj(lambda) * U[nu];
-    out -= Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
-
-    // C3-
-    tmp = lambda * U[nu];
-    out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu);
-
-    // C4-
-    out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * lambda;
-
-    return out;
-  }
-
-  static CloverField fillCloverYZ(const GaugeLinkField &F)
-  {
-    CloverField T(F.Grid());
-    T = Zero();
-    autoView(T_v,T,AcceleratorWrite);
-    autoView(F_v,F,AcceleratorRead);
-    accelerator_for(i, T.Grid()->oSites(),CloverField::vector_type::Nsimd(),
-    {
-      coalescedWrite(T_v[i]()(0, 1), coalescedRead(timesMinusI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(1, 0), coalescedRead(timesMinusI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(2, 3), coalescedRead(timesMinusI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(3, 2), coalescedRead(timesMinusI(F_v[i]()())));
-    });
-
-    return T;
-  }
-
-  static CloverField fillCloverXZ(const GaugeLinkField &F)
-  {
-    CloverField T(F.Grid());
-    T = Zero();
-    
-    autoView(T_v, T,AcceleratorWrite);
-    autoView(F_v, F,AcceleratorRead);
-    accelerator_for(i, T.Grid()->oSites(),CloverField::vector_type::Nsimd(),
-    {
-      coalescedWrite(T_v[i]()(0, 1), coalescedRead(-F_v[i]()()));
-      coalescedWrite(T_v[i]()(1, 0), coalescedRead(F_v[i]()()));
-      coalescedWrite(T_v[i]()(2, 3), coalescedRead(-F_v[i]()()));
-      coalescedWrite(T_v[i]()(3, 2), coalescedRead(F_v[i]()()));
-    });
-
-    return T;
-  }
-
-  static CloverField fillCloverXY(const GaugeLinkField &F)
-  {
-    CloverField T(F.Grid());
-    T = Zero();
-
-    autoView(T_v,T,AcceleratorWrite);
-    autoView(F_v,F,AcceleratorRead);
-    accelerator_for(i, T.Grid()->oSites(),CloverField::vector_type::Nsimd(),
-    {
-      coalescedWrite(T_v[i]()(0, 0), coalescedRead(timesMinusI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(1, 1), coalescedRead(timesI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(2, 2), coalescedRead(timesMinusI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(3, 3), coalescedRead(timesI(F_v[i]()())));
-    });
-
-    return T;
-  }
-
-  static CloverField fillCloverXT(const GaugeLinkField &F)
-  {
-    CloverField T(F.Grid());
-    T = Zero();
-
-    autoView( T_v , T, AcceleratorWrite);
-    autoView( F_v , F, AcceleratorRead);
-    accelerator_for(i, T.Grid()->oSites(),CloverField::vector_type::Nsimd(),
-    {
-      coalescedWrite(T_v[i]()(0, 1), coalescedRead(timesI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(1, 0), coalescedRead(timesI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(2, 3), coalescedRead(timesMinusI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(3, 2), coalescedRead(timesMinusI(F_v[i]()())));
-    });
-
-    return T;
-  }
-
-  static CloverField fillCloverYT(const GaugeLinkField &F)
-  {
-    CloverField T(F.Grid());
-    T = Zero();
-    
-    autoView( T_v ,T,AcceleratorWrite);
-    autoView( F_v ,F,AcceleratorRead);
-    accelerator_for(i, T.Grid()->oSites(),CloverField::vector_type::Nsimd(),
-    {
-      coalescedWrite(T_v[i]()(0, 1), coalescedRead(-(F_v[i]()())));
-      coalescedWrite(T_v[i]()(1, 0), coalescedRead((F_v[i]()())));
-      coalescedWrite(T_v[i]()(2, 3), coalescedRead((F_v[i]()())));
-      coalescedWrite(T_v[i]()(3, 2), coalescedRead(-(F_v[i]()())));
-    });
-
-    return T;
-  }
-
-  static CloverField fillCloverZT(const GaugeLinkField &F)
-  {
-    CloverField T(F.Grid());
-
-    T = Zero();
-
-    autoView( T_v , T,AcceleratorWrite);
-    autoView( F_v , F,AcceleratorRead);
-    accelerator_for(i, T.Grid()->oSites(),CloverField::vector_type::Nsimd(),
-    {
-      coalescedWrite(T_v[i]()(0, 0), coalescedRead(timesI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(1, 1), coalescedRead(timesMinusI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(2, 2), coalescedRead(timesMinusI(F_v[i]()())));
-      coalescedWrite(T_v[i]()(3, 3), coalescedRead(timesI(F_v[i]()())));
-    });
-
-    return T;
-  }
-
-  template<class _Spinor>
-  static accelerator_inline void multClover(_Spinor& phi, const SiteClover& C, const _Spinor& chi) {
-    auto CC = coalescedRead(C);
-    mult(&phi, &CC, &chi);
-  }
-
-  template<class _SpinorField>
-  inline void multCloverField(_SpinorField& out, const CloverField& C, const _SpinorField& phi) {
-    const int Nsimd = SiteSpinor::Nsimd();
-    autoView(out_v, out, AcceleratorWrite);
-    autoView(phi_v, phi, AcceleratorRead);
-    autoView(C_v,   C,   AcceleratorRead);
-    typedef decltype(coalescedRead(out_v[0])) calcSpinor;
-    accelerator_for(sss,out.Grid()->oSites(),Nsimd,{
-      calcSpinor tmp;
-      multClover(tmp,C_v[sss],phi_v(sss));
-      coalescedWrite(out_v[sss],tmp);
-    });
-  }
-};
-
-
-////////////////////////////////////////////////////////
-
-template<class Impl> class CompactWilsonCloverHelpers {
-public:
-
-  INHERIT_COMPACT_CLOVER_SIZES(Impl);
-
-  INHERIT_IMPL_TYPES(Impl);
-  INHERIT_CLOVER_TYPES(Impl);
-  INHERIT_COMPACT_CLOVER_TYPES(Impl);
-
-  #if 0
-  static accelerator_inline typename SiteCloverTriangle::vector_type triangle_elem(const SiteCloverTriangle& triangle, int block, int i, int j) {
-    assert(i != j);
-    if(i < j) {
-      return triangle()(block)(triangle_index(i, j));
-    } else { // i > j
-      return conjugate(triangle()(block)(triangle_index(i, j)));
-    }
-  }
-  #else
-  template<typename vobj>
-  static accelerator_inline vobj triangle_elem(const iImplCloverTriangle<vobj>& triangle, int block, int i, int j) {
-    assert(i != j);
-    if(i < j) {
-      return triangle()(block)(triangle_index(i, j));
-    } else { // i > j
-      return conjugate(triangle()(block)(triangle_index(i, j)));
-    }
-  }
-  #endif
-
-  static accelerator_inline int triangle_index(int i, int j) {
-    if(i == j)
-      return 0;
-    else if(i < j)
-      return Nred * (Nred - 1) / 2 - (Nred - i) * (Nred - i - 1) / 2 + j - i - 1;
-    else // i > j
-      return Nred * (Nred - 1) / 2 - (Nred - j) * (Nred - j - 1) / 2 + i - j - 1;
-  }
-
-  static void MooeeKernel_gpu(int                        Nsite,
-                              int                        Ls,
-                              const FermionField&        in,
-                              FermionField&              out,
-                              const CloverDiagonalField& diagonal,
-                              const CloverTriangleField& triangle) {
-    autoView(diagonal_v, diagonal, AcceleratorRead);
-    autoView(triangle_v, triangle, AcceleratorRead);
-    autoView(in_v,       in,       AcceleratorRead);
-    autoView(out_v,      out,      AcceleratorWrite);
-
-    typedef decltype(coalescedRead(out_v[0])) CalcSpinor;
-
-    const uint64_t NN = Nsite * Ls;
-
-    accelerator_for(ss, NN, Simd::Nsimd(), {
-      int sF = ss;
-      int sU = ss/Ls;
-      CalcSpinor res;
-      CalcSpinor in_t = in_v(sF);
-      auto diagonal_t = diagonal_v(sU);
-      auto triangle_t = triangle_v(sU);
-      for(int block=0; block<Nhs; block++) {
-        int s_start = block*Nhs;
-        for(int i=0; i<Nred; i++) {
-          int si = s_start + i/Nc, ci = i%Nc;
-          res()(si)(ci) = diagonal_t()(block)(i) * in_t()(si)(ci);
-          for(int j=0; j<Nred; j++) {
-            if (j == i) continue;
-            int sj = s_start + j/Nc, cj = j%Nc;
-            res()(si)(ci) = res()(si)(ci) + triangle_elem(triangle_t, block, i, j) * in_t()(sj)(cj);
-          };
-        };
-      };
-      coalescedWrite(out_v[sF], res);
-    });
-  }
-
-  static void MooeeKernel_cpu(int                        Nsite,
-                              int                        Ls,
-                              const FermionField&        in,
-                              FermionField&              out,
-                              const CloverDiagonalField& diagonal,
-                              const CloverTriangleField& triangle) {
-    autoView(diagonal_v, diagonal, CpuRead);
-    autoView(triangle_v, triangle, CpuRead);
-    autoView(in_v,       in,       CpuRead);
-    autoView(out_v,      out,      CpuWrite);
-
-    typedef SiteSpinor CalcSpinor;
-
-#if defined(A64FX) || defined(A64FXFIXEDSIZE)
-#define PREFETCH_CLOVER(BASE) {                                     \
-    uint64_t base;                                                  \
-    int pf_dist_L1 = 1;                                             \
-    int pf_dist_L2 = -5; /* -> penalty -> disable */                \
-                                                                    \
-    if ((pf_dist_L1 >= 0) && (sU + pf_dist_L1 < Nsite)) {           \
-      base = (uint64_t)&diag_t()(pf_dist_L1+BASE)(0);               \
-      svprfd(svptrue_b64(), (int64_t*)(base +    0), SV_PLDL1STRM); \
-      svprfd(svptrue_b64(), (int64_t*)(base +  256), SV_PLDL1STRM); \
-      svprfd(svptrue_b64(), (int64_t*)(base +  512), SV_PLDL1STRM); \
-      svprfd(svptrue_b64(), (int64_t*)(base +  768), SV_PLDL1STRM); \
-      svprfd(svptrue_b64(), (int64_t*)(base + 1024), SV_PLDL1STRM); \
-      svprfd(svptrue_b64(), (int64_t*)(base + 1280), SV_PLDL1STRM); \
-    }                                                               \
-                                                                    \
-    if ((pf_dist_L2 >= 0) && (sU + pf_dist_L2 < Nsite)) {           \
-      base = (uint64_t)&diag_t()(pf_dist_L2+BASE)(0);               \
-      svprfd(svptrue_b64(), (int64_t*)(base +    0), SV_PLDL2STRM); \
-      svprfd(svptrue_b64(), (int64_t*)(base +  256), SV_PLDL2STRM); \
-      svprfd(svptrue_b64(), (int64_t*)(base +  512), SV_PLDL2STRM); \
-      svprfd(svptrue_b64(), (int64_t*)(base +  768), SV_PLDL2STRM); \
-      svprfd(svptrue_b64(), (int64_t*)(base + 1024), SV_PLDL2STRM); \
-      svprfd(svptrue_b64(), (int64_t*)(base + 1280), SV_PLDL2STRM); \
-    }                                                               \
-  }
-// TODO: Implement/generalize this for other architectures
-// I played around a bit on KNL (see below) but didn't bring anything
-// #elif defined(AVX512)
-// #define PREFETCH_CLOVER(BASE) {                              \
-//     uint64_t base;                                           \
-//     int pf_dist_L1 = 1;                                      \
-//     int pf_dist_L2 = +4;                                     \
-//                                                              \
-//     if ((pf_dist_L1 >= 0) && (sU + pf_dist_L1 < Nsite)) {    \
-//       base = (uint64_t)&diag_t()(pf_dist_L1+BASE)(0);        \
-//       _mm_prefetch((const char*)(base +    0), _MM_HINT_T0); \
-//       _mm_prefetch((const char*)(base +   64), _MM_HINT_T0); \
-//       _mm_prefetch((const char*)(base +  128), _MM_HINT_T0); \
-//       _mm_prefetch((const char*)(base +  192), _MM_HINT_T0); \
-//       _mm_prefetch((const char*)(base +  256), _MM_HINT_T0); \
-//       _mm_prefetch((const char*)(base +  320), _MM_HINT_T0); \
-//     }                                                        \
-//                                                              \
-//     if ((pf_dist_L2 >= 0) && (sU + pf_dist_L2 < Nsite)) {    \
-//       base = (uint64_t)&diag_t()(pf_dist_L2+BASE)(0);        \
-//       _mm_prefetch((const char*)(base +    0), _MM_HINT_T1); \
-//       _mm_prefetch((const char*)(base +   64), _MM_HINT_T1); \
-//       _mm_prefetch((const char*)(base +  128), _MM_HINT_T1); \
-//       _mm_prefetch((const char*)(base +  192), _MM_HINT_T1); \
-//       _mm_prefetch((const char*)(base +  256), _MM_HINT_T1); \
-//       _mm_prefetch((const char*)(base +  320), _MM_HINT_T1); \
-//     }                                                        \
-//   }
-#else
-#define PREFETCH_CLOVER(BASE)
-#endif
-
-    const uint64_t NN = Nsite * Ls;
-
-    thread_for(ss, NN, {
-      int sF = ss;
-      int sU = ss/Ls;
-      CalcSpinor res;
-      CalcSpinor in_t = in_v[sF];
-      auto diag_t     = diagonal_v[sU]; // "diag" instead of "diagonal" here to make code below easier to read
-      auto triangle_t = triangle_v[sU];
-
-      // upper half
-      PREFETCH_CLOVER(0);
-
-      auto in_cc_0_0 = conjugate(in_t()(0)(0)); // Nils: reduces number
-      auto in_cc_0_1 = conjugate(in_t()(0)(1)); // of conjugates from
-      auto in_cc_0_2 = conjugate(in_t()(0)(2)); // 30 to 20
-      auto in_cc_1_0 = conjugate(in_t()(1)(0));
-      auto in_cc_1_1 = conjugate(in_t()(1)(1));
-
-      res()(0)(0) =               diag_t()(0)( 0) * in_t()(0)(0)
-                  +           triangle_t()(0)( 0) * in_t()(0)(1)
-                  +           triangle_t()(0)( 1) * in_t()(0)(2)
-                  +           triangle_t()(0)( 2) * in_t()(1)(0)
-                  +           triangle_t()(0)( 3) * in_t()(1)(1)
-                  +           triangle_t()(0)( 4) * in_t()(1)(2);
-
-      res()(0)(1) =           triangle_t()(0)( 0) * in_cc_0_0;
-      res()(0)(1) =               diag_t()(0)( 1) * in_t()(0)(1)
-                  +           triangle_t()(0)( 5) * in_t()(0)(2)
-                  +           triangle_t()(0)( 6) * in_t()(1)(0)
-                  +           triangle_t()(0)( 7) * in_t()(1)(1)
-                  +           triangle_t()(0)( 8) * in_t()(1)(2)
-                  + conjugate(       res()(0)( 1));
-
-      res()(0)(2) =           triangle_t()(0)( 1) * in_cc_0_0
-                  +           triangle_t()(0)( 5) * in_cc_0_1;
-      res()(0)(2) =               diag_t()(0)( 2) * in_t()(0)(2)
-                  +           triangle_t()(0)( 9) * in_t()(1)(0)
-                  +           triangle_t()(0)(10) * in_t()(1)(1)
-                  +           triangle_t()(0)(11) * in_t()(1)(2)
-                  + conjugate(       res()(0)( 2));
-
-      res()(1)(0) =           triangle_t()(0)( 2) * in_cc_0_0
-                  +           triangle_t()(0)( 6) * in_cc_0_1
-                  +           triangle_t()(0)( 9) * in_cc_0_2;
-      res()(1)(0) =               diag_t()(0)( 3) * in_t()(1)(0)
-                  +           triangle_t()(0)(12) * in_t()(1)(1)
-                  +           triangle_t()(0)(13) * in_t()(1)(2)
-                  + conjugate(       res()(1)( 0));
-
-      res()(1)(1) =           triangle_t()(0)( 3) * in_cc_0_0
-                  +           triangle_t()(0)( 7) * in_cc_0_1
-                  +           triangle_t()(0)(10) * in_cc_0_2
-                  +           triangle_t()(0)(12) * in_cc_1_0;
-      res()(1)(1) =               diag_t()(0)( 4) * in_t()(1)(1)
-                  +           triangle_t()(0)(14) * in_t()(1)(2)
-                  + conjugate(       res()(1)( 1));
-
-      res()(1)(2) =           triangle_t()(0)( 4) * in_cc_0_0
-                  +           triangle_t()(0)( 8) * in_cc_0_1
-                  +           triangle_t()(0)(11) * in_cc_0_2
-                  +           triangle_t()(0)(13) * in_cc_1_0
-                  +           triangle_t()(0)(14) * in_cc_1_1;
-      res()(1)(2) =               diag_t()(0)( 5) * in_t()(1)(2)
-                  + conjugate(       res()(1)( 2));
-
-      vstream(out_v[sF]()(0)(0), res()(0)(0));
-      vstream(out_v[sF]()(0)(1), res()(0)(1));
-      vstream(out_v[sF]()(0)(2), res()(0)(2));
-      vstream(out_v[sF]()(1)(0), res()(1)(0));
-      vstream(out_v[sF]()(1)(1), res()(1)(1));
-      vstream(out_v[sF]()(1)(2), res()(1)(2));
-
-      // lower half
-      PREFETCH_CLOVER(1);
-
-      auto in_cc_2_0 = conjugate(in_t()(2)(0));
-      auto in_cc_2_1 = conjugate(in_t()(2)(1));
-      auto in_cc_2_2 = conjugate(in_t()(2)(2));
-      auto in_cc_3_0 = conjugate(in_t()(3)(0));
-      auto in_cc_3_1 = conjugate(in_t()(3)(1));
-
-      res()(2)(0) =               diag_t()(1)( 0) * in_t()(2)(0)
-                  +           triangle_t()(1)( 0) * in_t()(2)(1)
-                  +           triangle_t()(1)( 1) * in_t()(2)(2)
-                  +           triangle_t()(1)( 2) * in_t()(3)(0)
-                  +           triangle_t()(1)( 3) * in_t()(3)(1)
-                  +           triangle_t()(1)( 4) * in_t()(3)(2);
-
-      res()(2)(1) =           triangle_t()(1)( 0) * in_cc_2_0;
-      res()(2)(1) =               diag_t()(1)( 1) * in_t()(2)(1)
-                  +           triangle_t()(1)( 5) * in_t()(2)(2)
-                  +           triangle_t()(1)( 6) * in_t()(3)(0)
-                  +           triangle_t()(1)( 7) * in_t()(3)(1)
-                  +           triangle_t()(1)( 8) * in_t()(3)(2)
-                  + conjugate(       res()(2)( 1));
-
-      res()(2)(2) =           triangle_t()(1)( 1) * in_cc_2_0
-                  +           triangle_t()(1)( 5) * in_cc_2_1;
-      res()(2)(2) =               diag_t()(1)( 2) * in_t()(2)(2)
-                  +           triangle_t()(1)( 9) * in_t()(3)(0)
-                  +           triangle_t()(1)(10) * in_t()(3)(1)
-                  +           triangle_t()(1)(11) * in_t()(3)(2)
-                  + conjugate(       res()(2)( 2));
-
-      res()(3)(0) =           triangle_t()(1)( 2) * in_cc_2_0
-                  +           triangle_t()(1)( 6) * in_cc_2_1
-                  +           triangle_t()(1)( 9) * in_cc_2_2;
-      res()(3)(0) =               diag_t()(1)( 3) * in_t()(3)(0)
-                  +           triangle_t()(1)(12) * in_t()(3)(1)
-                  +           triangle_t()(1)(13) * in_t()(3)(2)
-                  + conjugate(       res()(3)( 0));
-
-      res()(3)(1) =           triangle_t()(1)( 3) * in_cc_2_0
-                  +           triangle_t()(1)( 7) * in_cc_2_1
-                  +           triangle_t()(1)(10) * in_cc_2_2
-                  +           triangle_t()(1)(12) * in_cc_3_0;
-      res()(3)(1) =               diag_t()(1)( 4) * in_t()(3)(1)
-                  +           triangle_t()(1)(14) * in_t()(3)(2)
-                  + conjugate(       res()(3)( 1));
-
-      res()(3)(2) =           triangle_t()(1)( 4) * in_cc_2_0
-                  +           triangle_t()(1)( 8) * in_cc_2_1
-                  +           triangle_t()(1)(11) * in_cc_2_2
-                  +           triangle_t()(1)(13) * in_cc_3_0
-                  +           triangle_t()(1)(14) * in_cc_3_1;
-      res()(3)(2) =               diag_t()(1)( 5) * in_t()(3)(2)
-                  + conjugate(       res()(3)( 2));
-
-      vstream(out_v[sF]()(2)(0), res()(2)(0));
-      vstream(out_v[sF]()(2)(1), res()(2)(1));
-      vstream(out_v[sF]()(2)(2), res()(2)(2));
-      vstream(out_v[sF]()(3)(0), res()(3)(0));
-      vstream(out_v[sF]()(3)(1), res()(3)(1));
-      vstream(out_v[sF]()(3)(2), res()(3)(2));
-    });
-  }
-
-  static void MooeeKernel(int                        Nsite,
-                          int                        Ls,
-                          const FermionField&        in,
-                          FermionField&              out,
-                          const CloverDiagonalField& diagonal,
-                          const CloverTriangleField& triangle) {
-#if defined(GRID_CUDA) || defined(GRID_HIP)
-    MooeeKernel_gpu(Nsite, Ls, in, out, diagonal, triangle);
-#else
-    MooeeKernel_cpu(Nsite, Ls, in, out, diagonal, triangle);
-#endif
-  }
-
-  static void Invert(const CloverDiagonalField& diagonal,
-                     const CloverTriangleField& triangle,
-                     CloverDiagonalField&       diagonalInv,
-                     CloverTriangleField&       triangleInv) {
-    conformable(diagonal, diagonalInv);
-    conformable(triangle, triangleInv);
-    conformable(diagonal, triangle);
-
-    diagonalInv.Checkerboard() = diagonal.Checkerboard();
-    triangleInv.Checkerboard() = triangle.Checkerboard();
-
-    GridBase* grid = diagonal.Grid();
-
-    long lsites = grid->lSites();
-
-    typedef typename SiteCloverDiagonal::scalar_object scalar_object_diagonal;
-    typedef typename SiteCloverTriangle::scalar_object scalar_object_triangle;
-
-    autoView(diagonal_v,  diagonal,  CpuRead);
-    autoView(triangle_v,  triangle,  CpuRead);
-    autoView(diagonalInv_v, diagonalInv, CpuWrite);
-    autoView(triangleInv_v, triangleInv, CpuWrite);
-
-    thread_for(site, lsites, { // NOTE: Not on GPU because of Eigen & (peek/poke)LocalSite
-      Eigen::MatrixXcd clover_inv_eigen = Eigen::MatrixXcd::Zero(Ns*Nc, Ns*Nc);
-      Eigen::MatrixXcd clover_eigen = Eigen::MatrixXcd::Zero(Ns*Nc, Ns*Nc);
-
-      scalar_object_diagonal diagonal_tmp     = Zero();
-      scalar_object_diagonal diagonal_inv_tmp = Zero();
-      scalar_object_triangle triangle_tmp     = Zero();
-      scalar_object_triangle triangle_inv_tmp = Zero();
-
-      Coordinate lcoor;
-      grid->LocalIndexToLocalCoor(site, lcoor);
-
-      peekLocalSite(diagonal_tmp, diagonal_v, lcoor);
-      peekLocalSite(triangle_tmp, triangle_v, lcoor);
-
-      // TODO: can we save time here by inverting the two 6x6 hermitian matrices separately?
-      for (long s_row=0;s_row<Ns;s_row++) {
-        for (long s_col=0;s_col<Ns;s_col++) {
-          if(abs(s_row - s_col) > 1 || s_row + s_col == 3) continue;
-          int block       = s_row / Nhs;
-          int s_row_block = s_row % Nhs;
-          int s_col_block = s_col % Nhs;
-          for (long c_row=0;c_row<Nc;c_row++) {
-            for (long c_col=0;c_col<Nc;c_col++) {
-              int i = s_row_block * Nc + c_row;
-              int j = s_col_block * Nc + c_col;
-              if(i == j)
-                clover_eigen(s_row*Nc+c_row, s_col*Nc+c_col) = static_cast<ComplexD>(TensorRemove(diagonal_tmp()(block)(i)));
-              else
-                clover_eigen(s_row*Nc+c_row, s_col*Nc+c_col) = static_cast<ComplexD>(TensorRemove(triangle_elem(triangle_tmp, block, i, j)));
-            }
-          }
-        }
-      }
-
-      clover_inv_eigen = clover_eigen.inverse();
-
-      for (long s_row=0;s_row<Ns;s_row++) {
-        for (long s_col=0;s_col<Ns;s_col++) {
-          if(abs(s_row - s_col) > 1 || s_row + s_col == 3) continue;
-          int block       = s_row / Nhs;
-          int s_row_block = s_row % Nhs;
-          int s_col_block = s_col % Nhs;
-          for (long c_row=0;c_row<Nc;c_row++) {
-            for (long c_col=0;c_col<Nc;c_col++) {
-              int i = s_row_block * Nc + c_row;
-              int j = s_col_block * Nc + c_col;
-              if(i == j)
-                diagonal_inv_tmp()(block)(i) = clover_inv_eigen(s_row*Nc+c_row, s_col*Nc+c_col);
-              else if(i < j)
-                triangle_inv_tmp()(block)(triangle_index(i, j)) = clover_inv_eigen(s_row*Nc+c_row, s_col*Nc+c_col);
-              else
-                continue;
-            }
-          }
-        }
-      }
-
-      pokeLocalSite(diagonal_inv_tmp, diagonalInv_v, lcoor);
-      pokeLocalSite(triangle_inv_tmp, triangleInv_v, lcoor);
-    });
-  }
-
-  static void ConvertLayout(const CloverField&   full,
-                            CloverDiagonalField& diagonal,
-                            CloverTriangleField& triangle) {
-    conformable(full, diagonal);
-    conformable(full, triangle);
-
-    diagonal.Checkerboard() = full.Checkerboard();
-    triangle.Checkerboard() = full.Checkerboard();
-
-    autoView(full_v,     full,     AcceleratorRead);
-    autoView(diagonal_v, diagonal, AcceleratorWrite);
-    autoView(triangle_v, triangle, AcceleratorWrite);
-
-    // NOTE: this function cannot be 'private' since nvcc forbids this for kernels
-    accelerator_for(ss, full.Grid()->oSites(), 1, {
-      for(int s_row = 0; s_row < Ns; s_row++) {
-        for(int s_col = 0; s_col < Ns; s_col++) {
-          if(abs(s_row - s_col) > 1 || s_row + s_col == 3) continue;
-          int block       = s_row / Nhs;
-          int s_row_block = s_row % Nhs;
-          int s_col_block = s_col % Nhs;
-          for(int c_row = 0; c_row < Nc; c_row++) {
-            for(int c_col = 0; c_col < Nc; c_col++) {
-              int i = s_row_block * Nc + c_row;
-              int j = s_col_block * Nc + c_col;
-              if(i == j)
-                diagonal_v[ss]()(block)(i) = full_v[ss]()(s_row, s_col)(c_row, c_col);
-              else if(i < j)
-                triangle_v[ss]()(block)(triangle_index(i, j)) = full_v[ss]()(s_row, s_col)(c_row, c_col);
-              else
-                continue;
-            }
-          }
-        }
-      }
-    });
-  }
-
-
-  static void ConvertLayout(const CloverDiagonalField& diagonal,
-                            const CloverTriangleField& triangle,
-                            CloverField&               full) {
-    conformable(full, diagonal);
-    conformable(full, triangle);
-
-    full.Checkerboard() = diagonal.Checkerboard();
-
-    full = Zero();
-
-    autoView(diagonal_v, diagonal, AcceleratorRead);
-    autoView(triangle_v, triangle, AcceleratorRead);
-    autoView(full_v,     full,     AcceleratorWrite);
-
-    // NOTE: this function cannot be 'private' since nvcc forbids this for kernels
-    accelerator_for(ss, full.Grid()->oSites(), 1, {
-      for(int s_row = 0; s_row < Ns; s_row++) {
-        for(int s_col = 0; s_col < Ns; s_col++) {
-          if(abs(s_row - s_col) > 1 || s_row + s_col == 3) continue;
-          int block       = s_row / Nhs;
-          int s_row_block = s_row % Nhs;
-          int s_col_block = s_col % Nhs;
-          for(int c_row = 0; c_row < Nc; c_row++) {
-            for(int c_col = 0; c_col < Nc; c_col++) {
-              int i = s_row_block * Nc + c_row;
-              int j = s_col_block * Nc + c_col;
-              if(i == j)
-                full_v[ss]()(s_row, s_col)(c_row, c_col) = diagonal_v[ss]()(block)(i);
-              else
-                full_v[ss]()(s_row, s_col)(c_row, c_col) = triangle_elem(triangle_v[ss], block, i, j);
-            }
-          }
-        }
-      }
-    });
-  }
-
-  static void ModifyBoundaries(CloverDiagonalField& diagonal, CloverTriangleField& triangle, RealD csw_t, RealD cF, RealD diag_mass) {
-    // Checks/grid
-    double t0 = usecond();
-    conformable(diagonal, triangle);
-    GridBase* grid = diagonal.Grid();
-
-    // Determine the boundary coordinates/sites
-    double t1 = usecond();
-    int t_dir = Nd - 1;
-    Lattice<iScalar<vInteger>> t_coor(grid);
-    LatticeCoordinate(t_coor, t_dir);
-    int T = grid->GlobalDimensions()[t_dir];
-
-    // Set off-diagonal parts at boundary to zero -- OK
-    double t2 = usecond();
-    CloverTriangleField zeroTriangle(grid);
-    zeroTriangle.Checkerboard() = triangle.Checkerboard();
-    zeroTriangle = Zero();
-    triangle = where(t_coor == 0,   zeroTriangle, triangle);
-    triangle = where(t_coor == T-1, zeroTriangle, triangle);
-
-    // Set diagonal to unity (scaled correctly) -- OK
-    double t3 = usecond();
-    CloverDiagonalField tmp(grid);
-    tmp.Checkerboard() = diagonal.Checkerboard();
-    tmp                = -1.0 * csw_t + diag_mass;
-    diagonal           = where(t_coor == 0,   tmp, diagonal);
-    diagonal           = where(t_coor == T-1, tmp, diagonal);
-
-    // Correct values next to boundary
-    double t4 = usecond();
-    if(cF != 1.0) {
-      tmp = cF - 1.0;
-      tmp += diagonal;
-      diagonal = where(t_coor == 1,   tmp, diagonal);
-      diagonal = where(t_coor == T-2, tmp, diagonal);
-    }
-
-    // Report timings
-    double t5 = usecond();
-#if 0
-    std::cout << GridLogMessage << "CompactWilsonCloverHelpers::ModifyBoundaries timings:"
-              << " checks = "          << (t1 - t0) / 1e6
-              << ", coordinate = "     << (t2 - t1) / 1e6
-              << ", off-diag zero = "  << (t3 - t2) / 1e6
-              << ", diagonal unity = " << (t4 - t3) / 1e6
-              << ", near-boundary = "  << (t5 - t4) / 1e6
-              << ", total = "          << (t5 - t0) / 1e6
-              << std::endl;
-#endif
-  }
-
-  template<class Field, class Mask>
-  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();
-    autoView(f_v, f, AcceleratorWrite);
-    autoView(m_v, m, AcceleratorRead);
-    // NOTE: this function cannot be 'private' since nvcc forbids this for kernels
-    accelerator_for(ss, Nsite, Nsimd, {
-      coalescedWrite(f_v[ss], m_v(ss) * f_v(ss));
-    });
-  }
-
-  template<class MaskField>
-  static void SetupMasks(MaskField& full, MaskField& even, MaskField& odd) {
-    assert(even.Grid()->_isCheckerBoarded && even.Checkerboard() == Even);
-    assert(odd.Grid()->_isCheckerBoarded  && odd.Checkerboard()  == Odd);
-    assert(!full.Grid()->_isCheckerBoarded);
-
-    GridBase* grid = full.Grid();
-    int t_dir = Nd-1;
-    Lattice<iScalar<vInteger>> t_coor(grid);
-    LatticeCoordinate(t_coor, t_dir);
-    int T = grid->GlobalDimensions()[t_dir];
-
-    MaskField zeroMask(grid); zeroMask = Zero();
-    full = 1.0;
-    full = where(t_coor == 0,   zeroMask, full);
-    full = where(t_coor == T-1, zeroMask, full);
-
-    pickCheckerboard(Even, even, full);
-    pickCheckerboard(Odd,  odd,  full);
-  }
-};
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/WilsonCloverTypes.h

@@ -1,90 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./lib/qcd/action/fermion/WilsonCloverTypes.h
-
-    Copyright (C) 2021 - 2022
-
-    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-
-#pragma once
-
-NAMESPACE_BEGIN(Grid);
-
-template<class Impl>
-class WilsonCloverTypes {
-public:
-  INHERIT_IMPL_TYPES(Impl);
-
-  template <typename vtype> using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
-
-  typedef iImplClover<Simd> SiteClover;
-
-  typedef Lattice<SiteClover> CloverField;
-};
-
-template<class Impl>
-class CompactWilsonCloverTypes {
-public:
-  INHERIT_IMPL_TYPES(Impl);
-
-  static constexpr int Nred      = Nc * Nhs;        // 6
-  static constexpr int Nblock    = Nhs;             // 2
-  static constexpr int Ndiagonal = Nred;            // 6
-  static constexpr int Ntriangle = (Nred - 1) * Nc; // 15
-
-  template<typename vtype> using iImplCloverDiagonal = iScalar<iVector<iVector<vtype, Ndiagonal>, Nblock>>;
-  template<typename vtype> using iImplCloverTriangle = iScalar<iVector<iVector<vtype, Ntriangle>, Nblock>>;
-
-  typedef iImplCloverDiagonal<Simd> SiteCloverDiagonal;
-  typedef iImplCloverTriangle<Simd> SiteCloverTriangle;
-  typedef iSinglet<Simd>            SiteMask;
-
-  typedef Lattice<SiteCloverDiagonal> CloverDiagonalField;
-  typedef Lattice<SiteCloverTriangle> CloverTriangleField;
-  typedef Lattice<SiteMask>           MaskField;
-};
-
-#define INHERIT_CLOVER_TYPES(Impl)                                 \
-  typedef typename WilsonCloverTypes<Impl>::SiteClover SiteClover; \
-  typedef typename WilsonCloverTypes<Impl>::CloverField CloverField;
-
-#define INHERIT_COMPACT_CLOVER_TYPES(Impl) \
-  typedef typename CompactWilsonCloverTypes<Impl>::SiteCloverDiagonal  SiteCloverDiagonal; \
-  typedef typename CompactWilsonCloverTypes<Impl>::SiteCloverTriangle  SiteCloverTriangle; \
-  typedef typename CompactWilsonCloverTypes<Impl>::SiteMask            SiteMask; \
-  typedef typename CompactWilsonCloverTypes<Impl>::CloverDiagonalField CloverDiagonalField; \
-  typedef typename CompactWilsonCloverTypes<Impl>::CloverTriangleField CloverTriangleField; \
-  typedef typename CompactWilsonCloverTypes<Impl>::MaskField           MaskField; \
-  /* ugly duplication but needed inside functionality classes */ \
-  template<typename vtype> using iImplCloverDiagonal = \
-    iScalar<iVector<iVector<vtype, CompactWilsonCloverTypes<Impl>::Ndiagonal>, CompactWilsonCloverTypes<Impl>::Nblock>>; \
-  template<typename vtype> using iImplCloverTriangle = \
-    iScalar<iVector<iVector<vtype, CompactWilsonCloverTypes<Impl>::Ntriangle>, CompactWilsonCloverTypes<Impl>::Nblock>>;
-
-#define INHERIT_COMPACT_CLOVER_SIZES(Impl)                                    \
-  static constexpr int Nred      = CompactWilsonCloverTypes<Impl>::Nred;      \
-  static constexpr int Nblock    = CompactWilsonCloverTypes<Impl>::Nblock;    \
-  static constexpr int Ndiagonal = CompactWilsonCloverTypes<Impl>::Ndiagonal; \
-  static constexpr int Ntriangle = CompactWilsonCloverTypes<Impl>::Ntriangle;
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/WilsonCompressor.h

@@ -117,19 +117,19 @@ public:
     typedef decltype(coalescedRead(*in))    sobj;
     typedef decltype(coalescedRead(*out0)) hsobj;
 
-    constexpr unsigned int Nsimd = vobj::Nsimd();
+    unsigned int Nsimd = vobj::Nsimd();
     unsigned int mask = Nsimd >> (type + 1);
     int lane = acceleratorSIMTlane(Nsimd);
     int j0 = lane &(~mask); // inner coor zero
     int j1 = lane |(mask) ; // inner coor one
-    const vobj *vp0 = &in[k];  // out0[j] = merge low bit of type from in[k] and in[m] 
+    const vobj *vp0 = &in[k];
-    const vobj *vp1 = &in[m];  // out1[j] = merge hi  bit of type from in[k] and in[m]
+    const vobj *vp1 = &in[m];
-    const vobj *vp = (lane&mask) ? vp1:vp0;// if my lane has high bit take vp1, low bit take vp0
+    const vobj *vp = (lane&mask) ? vp1:vp0;
-    auto sa = coalescedRead(*vp,j0); // lane to read for out 0, NB 50% read coalescing
+    auto sa = coalescedRead(*vp,j0);
-    auto sb = coalescedRead(*vp,j1); // lane to read for out 1
+    auto sb = coalescedRead(*vp,j1);
     hsobj psa, psb;
-    projector::Proj(psa,sa,mu,dag);  // spin project the result0
+    projector::Proj(psa,sa,mu,dag);
-    projector::Proj(psb,sb,mu,dag);  // spin project the result1
+    projector::Proj(psb,sb,mu,dag);
     coalescedWrite(out0[j],psa);
     coalescedWrite(out1[j],psb);
 #else

Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h

@@ -47,7 +47,7 @@ CayleyFermion5D<Impl>::CayleyFermion5D(GaugeField &_Umu,
 			FiveDimRedBlackGrid,
 			FourDimGrid,
 			FourDimRedBlackGrid,_M5,p),
-  mass_plus(_mass), mass_minus(_mass)
+  mass(_mass)
 { 
 }
 
@@ -209,8 +209,8 @@ void CayleyFermion5D<Impl>::M5D   (const FermionField &psi, FermionField &chi)
 {
   int Ls=this->Ls;
   Vector<Coeff_t> diag (Ls,1.0);
-  Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1]=mass_minus;
+  Vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1]=mass;
-  Vector<Coeff_t> lower(Ls,-1.0); lower[0]   =mass_plus;
+  Vector<Coeff_t> lower(Ls,-1.0); lower[0]   =mass;
   M5D(psi,chi,chi,lower,diag,upper);
 }
 template<class Impl>
@@ -220,8 +220,8 @@ void CayleyFermion5D<Impl>::Meooe5D    (const FermionField &psi, FermionField &D
   Vector<Coeff_t> diag = bs;
   Vector<Coeff_t> upper= cs;
   Vector<Coeff_t> lower= cs; 
-  upper[Ls-1]=-mass_minus*upper[Ls-1];
+  upper[Ls-1]=-mass*upper[Ls-1];
-  lower[0]   =-mass_plus*lower[0];
+  lower[0]   =-mass*lower[0];
   M5D(psi,psi,Din,lower,diag,upper);
 }
 // FIXME Redunant with the above routine; check this and eliminate
@@ -235,8 +235,8 @@ template<class Impl> void CayleyFermion5D<Impl>::Meo5D     (const FermionField &
     upper[i]=-ceo[i];
     lower[i]=-ceo[i];
   }
-  upper[Ls-1]=-mass_minus*upper[Ls-1];
+  upper[Ls-1]=-mass*upper[Ls-1];
-  lower[0]   =-mass_plus*lower[0];
+  lower[0]   =-mass*lower[0];
   M5D(psi,psi,chi,lower,diag,upper);
 }
 template<class Impl>
@@ -250,8 +250,8 @@ void CayleyFermion5D<Impl>::Mooee       (const FermionField &psi, FermionField &
     upper[i]=-cee[i];
     lower[i]=-cee[i];
   }
-  upper[Ls-1]=-mass_minus*upper[Ls-1];
+  upper[Ls-1]=-mass*upper[Ls-1];
-  lower[0]   =-mass_plus*lower[0];
+  lower[0]   =-mass*lower[0];
   M5D(psi,psi,chi,lower,diag,upper);
 }
 template<class Impl>
@@ -266,9 +266,9 @@ void CayleyFermion5D<Impl>::MooeeDag    (const FermionField &psi, FermionField &
     // Assemble the 5d matrix
     if ( s==0 ) {
       upper[s] = -cee[s+1] ;
-      lower[s] = mass_minus*cee[Ls-1];
+      lower[s] = mass*cee[Ls-1];
     } else if ( s==(Ls-1)) { 
-      upper[s] = mass_plus*cee[0];
+      upper[s] = mass*cee[0];
       lower[s] = -cee[s-1];
     } else {
       upper[s]=-cee[s+1];
@@ -291,8 +291,8 @@ void CayleyFermion5D<Impl>::M5Ddag (const FermionField &psi, FermionField &chi)
   Vector<Coeff_t> diag(Ls,1.0);
   Vector<Coeff_t> upper(Ls,-1.0);
   Vector<Coeff_t> lower(Ls,-1.0);
-  upper[Ls-1]=-mass_plus*upper[Ls-1];
+  upper[Ls-1]=-mass*upper[Ls-1];
-  lower[0]   =-mass_minus*lower[0];
+  lower[0]   =-mass*lower[0];
   M5Ddag(psi,chi,chi,lower,diag,upper);
 }
 
@@ -307,9 +307,9 @@ void CayleyFermion5D<Impl>::MeooeDag5D    (const FermionField &psi, FermionField
   for (int s=0;s<Ls;s++){
     if ( s== 0 ) {
       upper[s] = cs[s+1];
-      lower[s] =-mass_minus*cs[Ls-1];
+      lower[s] =-mass*cs[Ls-1];
     } else if ( s==(Ls-1) ) { 
-      upper[s] =-mass_plus*cs[0];
+      upper[s] =-mass*cs[0];
       lower[s] = cs[s-1];
     } else { 
       upper[s] = cs[s+1];
@@ -552,7 +552,7 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,Vector<Coeff_t
       
       lee[i] =-cee[i+1]/bee[i]; // sub-diag entry on the ith column
       
-      leem[i]=mass_minus*cee[Ls-1]/bee[0];
+      leem[i]=mass*cee[Ls-1]/bee[0];
       for(int j=0;j<i;j++) {
 	assert(bee[j+1]!=Coeff_t(0.0));
 	leem[i]*= aee[j]/bee[j+1];
@@ -560,7 +560,7 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,Vector<Coeff_t
       
       uee[i] =-aee[i]/bee[i];   // up-diag entry on the ith row
       
-      ueem[i]=mass_plus;
+      ueem[i]=mass;
       for(int j=1;j<=i;j++) ueem[i]*= cee[j]/bee[j];
       ueem[i]*= aee[0]/bee[0];
       
@@ -573,7 +573,7 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,Vector<Coeff_t
   }
 	
   { 
-    Coeff_t delta_d=mass_minus*cee[Ls-1];
+    Coeff_t delta_d=mass*cee[Ls-1];
     for(int j=0;j<Ls-1;j++) {
       assert(bee[j] != Coeff_t(0.0));
       delta_d *= cee[j]/bee[j];
@@ -642,10 +642,6 @@ void CayleyFermion5D<Impl>::ContractConservedCurrent( PropagatorField &q_in_1,
 						      Current curr_type,
 						      unsigned int mu)
 {
-
-  assert(mass_plus == mass_minus);
-  RealD mass = mass_plus;
-  
 #if (!defined(GRID_HIP))
   Gamma::Algebra Gmu [] = {
     Gamma::Algebra::GammaX,
@@ -781,8 +777,6 @@ void CayleyFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
   assert(mu>=0);
   assert(mu<Nd);
 
-  assert(mass_plus == mass_minus);
-  RealD mass = mass_plus;
 
 #if 0
   int tshift = (mu == Nd-1) ? 1 : 0;

Grid/qcd/action/fermion/implementation/CompactWilsonCloverFermionImplementation.h

@@ -1,377 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./lib/qcd/action/fermion/CompactWilsonCloverFermionImplementation.h
-
-    Copyright (C) 2017 - 2022
-
-    Author: paboyle <paboyle@ph.ed.ac.uk>
-    Author: Guido Cossu <guido.cossu@ed.ac.uk>
-    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-/*  END LEGAL */
-
-#include <Grid/Grid.h>
-#include <Grid/qcd/spin/Dirac.h>
-#include <Grid/qcd/action/fermion/CompactWilsonCloverFermion.h>
-
-
-NAMESPACE_BEGIN(Grid);
-template<class Impl, class CloverHelpers>
-CompactWilsonCloverFermion<Impl, CloverHelpers>::CompactWilsonCloverFermion(GaugeField& _Umu,
-                                                                            GridCartesian& Fgrid,
-                                                                            GridRedBlackCartesian& Hgrid,
-                                                                            const RealD _mass,
-                                                                            const RealD _csw_r,
-                                                                            const RealD _csw_t,
-                                                                            const RealD _cF,
-                                                                            const WilsonAnisotropyCoefficients& clover_anisotropy,
-                                                                            const ImplParams& impl_p)
-  : WilsonBase(_Umu, Fgrid, Hgrid, _mass, impl_p, clover_anisotropy)
-  , csw_r(_csw_r)
-  , csw_t(_csw_t)
-  , cF(_cF)
-  , fixedBoundaries(impl_p.boundary_phases[Nd-1] == 0.0)
-  , Diagonal(&Fgrid),        Triangle(&Fgrid)
-  , DiagonalEven(&Hgrid),    TriangleEven(&Hgrid)
-  , DiagonalOdd(&Hgrid),     TriangleOdd(&Hgrid)
-  , DiagonalInv(&Fgrid),     TriangleInv(&Fgrid)
-  , DiagonalInvEven(&Hgrid), TriangleInvEven(&Hgrid)
-  , DiagonalInvOdd(&Hgrid),  TriangleInvOdd(&Hgrid)
-  , Tmp(&Fgrid)
-  , BoundaryMask(&Fgrid)
-  , BoundaryMaskEven(&Hgrid), BoundaryMaskOdd(&Hgrid)
-{
-  assert(Nd == 4 && Nc == 3 && Ns == 4 && Impl::Dimension == 3);
-
-  csw_r *= 0.5;
-  csw_t *= 0.5;
-  if (clover_anisotropy.isAnisotropic)
-    csw_r /= clover_anisotropy.xi_0;
-
-  ImportGauge(_Umu);
-  if (fixedBoundaries) {
-    this->BoundaryMaskEven.Checkerboard() = Even;
-    this->BoundaryMaskOdd.Checkerboard() = Odd;
-    CompactHelpers::SetupMasks(this->BoundaryMask, this->BoundaryMaskEven, this->BoundaryMaskOdd);
-  }
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::Dhop(const FermionField& in, FermionField& out, int dag) {
-  WilsonBase::Dhop(in, out, dag);
-  if(fixedBoundaries) ApplyBoundaryMask(out);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopOE(const FermionField& in, FermionField& out, int dag) {
-  WilsonBase::DhopOE(in, out, dag);
-  if(fixedBoundaries) ApplyBoundaryMask(out);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopEO(const FermionField& in, FermionField& out, int dag) {
-  WilsonBase::DhopEO(in, out, dag);
-  if(fixedBoundaries) ApplyBoundaryMask(out);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopDir(const FermionField& in, FermionField& out, int dir, int disp) {
-  WilsonBase::DhopDir(in, out, dir, disp);
-  if(this->fixedBoundaries) ApplyBoundaryMask(out);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopDirAll(const FermionField& in, std::vector<FermionField>& out) {
-  WilsonBase::DhopDirAll(in, out);
-  if(this->fixedBoundaries) {
-    for(auto& o : out) ApplyBoundaryMask(o);
-  }
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::M(const FermionField& in, FermionField& out) {
-  out.Checkerboard() = in.Checkerboard();
-  WilsonBase::Dhop(in, out, DaggerNo); // call base to save applying bc
-  Mooee(in, Tmp);
-  axpy(out, 1.0, out, Tmp);
-  if(fixedBoundaries) ApplyBoundaryMask(out);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::Mdag(const FermionField& in, FermionField& out) {
-  out.Checkerboard() = in.Checkerboard();
-  WilsonBase::Dhop(in, out, DaggerYes);  // call base to save applying bc
-  MooeeDag(in, Tmp);
-  axpy(out, 1.0, out, Tmp);
-  if(fixedBoundaries) ApplyBoundaryMask(out);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::Meooe(const FermionField& in, FermionField& out) {
-  WilsonBase::Meooe(in, out);
-  if(fixedBoundaries) ApplyBoundaryMask(out);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::MeooeDag(const FermionField& in, FermionField& out) {
-  WilsonBase::MeooeDag(in, out);
-  if(fixedBoundaries) ApplyBoundaryMask(out);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::Mooee(const FermionField& in, FermionField& out) {
-  if(in.Grid()->_isCheckerBoarded) {
-    if(in.Checkerboard() == Odd) {
-      MooeeInternal(in, out, DiagonalOdd, TriangleOdd);
-    } else {
-      MooeeInternal(in, out, DiagonalEven, TriangleEven);
-    }
-  } else {
-    MooeeInternal(in, out, Diagonal, Triangle);
-  }
-  if(fixedBoundaries) ApplyBoundaryMask(out);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeDag(const FermionField& in, FermionField& out) {
-  Mooee(in, out); // blocks are hermitian
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeInv(const FermionField& in, FermionField& out) {
-  if(in.Grid()->_isCheckerBoarded) {
-    if(in.Checkerboard() == Odd) {
-      MooeeInternal(in, out, DiagonalInvOdd, TriangleInvOdd);
-    } else {
-      MooeeInternal(in, out, DiagonalInvEven, TriangleInvEven);
-    }
-  } else {
-    MooeeInternal(in, out, DiagonalInv, TriangleInv);
-  }
-  if(fixedBoundaries) ApplyBoundaryMask(out);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeInvDag(const FermionField& in, FermionField& out) {
-  MooeeInv(in, out); // blocks are hermitian
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::Mdir(const FermionField& in, FermionField& out, int dir, int disp) {
-  DhopDir(in, out, dir, disp);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::MdirAll(const FermionField& in, std::vector<FermionField>& out) {
-  DhopDirAll(in, out);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::MDeriv(GaugeField& force, const FermionField& X, const FermionField& Y, int dag) {
-  assert(!fixedBoundaries); // TODO check for changes required for open bc
-
-  // NOTE: code copied from original clover term
-  conformable(X.Grid(), Y.Grid());
-  conformable(X.Grid(), force.Grid());
-  GaugeLinkField force_mu(force.Grid()), lambda(force.Grid());
-  GaugeField clover_force(force.Grid());
-  PropagatorField Lambda(force.Grid());
-
-  // Guido: Here we are hitting some performance issues:
-  // need to extract the components of the DoubledGaugeField
-  // for each call
-  // Possible solution
-  // Create a vector object to store them? (cons: wasting space)
-  std::vector<GaugeLinkField> U(Nd, this->Umu.Grid());
-
-  Impl::extractLinkField(U, this->Umu);
-
-  force = Zero();
-  // Derivative of the Wilson hopping term
-  this->DhopDeriv(force, X, Y, dag);
-
-  ///////////////////////////////////////////////////////////
-  // Clover term derivative
-  ///////////////////////////////////////////////////////////
-  Impl::outerProductImpl(Lambda, X, Y);
-  //std::cout << "Lambda:" << Lambda << std::endl;
-
-  Gamma::Algebra sigma[] = {
-      Gamma::Algebra::SigmaXY,
-      Gamma::Algebra::SigmaXZ,
-      Gamma::Algebra::SigmaXT,
-      Gamma::Algebra::MinusSigmaXY,
-      Gamma::Algebra::SigmaYZ,
-      Gamma::Algebra::SigmaYT,
-      Gamma::Algebra::MinusSigmaXZ,
-      Gamma::Algebra::MinusSigmaYZ,
-      Gamma::Algebra::SigmaZT,
-      Gamma::Algebra::MinusSigmaXT,
-      Gamma::Algebra::MinusSigmaYT,
-      Gamma::Algebra::MinusSigmaZT};
-
-  /*
-    sigma_{\mu \nu}=
-    | 0         sigma[0]  sigma[1]  sigma[2] |
-    | sigma[3]    0       sigma[4]  sigma[5] |
-    | sigma[6]  sigma[7]     0      sigma[8] |
-    | sigma[9]  sigma[10] sigma[11]   0      |
-  */
-
-  int count = 0;
-  clover_force = Zero();
-  for (int mu = 0; mu < 4; mu++)
-  {
-    force_mu = Zero();
-    for (int nu = 0; nu < 4; nu++)
-    {
-      if (mu == nu)
-        continue;
-
-      RealD factor;
-      if (nu == 4 || mu == 4)
-      {
-        factor = 2.0 * csw_t;
-      }
-      else
-      {
-        factor = 2.0 * csw_r;
-      }
-      PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
-      Impl::TraceSpinImpl(lambda, Slambda);                   // traceSpin ok
-      force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu);   // checked
-      count++;
-    }
-
-    pokeLorentz(clover_force, U[mu] * force_mu, mu);
-  }
-  //clover_force *= csw;
-  force += clover_force;
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
-  assert(0);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::MeeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
-  assert(0);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeInternal(const FermionField&        in,
-                    FermionField&              out,
-                    const CloverDiagonalField& diagonal,
-                    const CloverTriangleField& triangle) {
-  assert(in.Checkerboard() == Odd || in.Checkerboard() == Even);
-  out.Checkerboard() = in.Checkerboard();
-  conformable(in, out);
-  conformable(in, diagonal);
-  conformable(in, triangle);
-
-  CompactHelpers::MooeeKernel(diagonal.oSites(), 1, in, out, diagonal, triangle);
-}
-
-template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl, CloverHelpers>::ImportGauge(const GaugeField& _Umu) {
-  // NOTE: parts copied from original implementation
-
-  // Import gauge into base class
-  double t0 = usecond();
-  WilsonBase::ImportGauge(_Umu); // NOTE: called here and in wilson constructor -> performed twice, but can't avoid that
-
-  // Initialize temporary variables
-  double t1 = usecond();
-  conformable(_Umu.Grid(), this->GaugeGrid());
-  GridBase* grid = _Umu.Grid();
-  typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid);
-  CloverField TmpOriginal(grid);
-  CloverField TmpInverse(grid);
-
-  // Compute the field strength terms mu>nu
-  double t2 = usecond();
-  WilsonLoops<Impl>::FieldStrength(Bx, _Umu, Zdir, Ydir);
-  WilsonLoops<Impl>::FieldStrength(By, _Umu, Zdir, Xdir);
-  WilsonLoops<Impl>::FieldStrength(Bz, _Umu, Ydir, Xdir);
-  WilsonLoops<Impl>::FieldStrength(Ex, _Umu, Tdir, Xdir);
-  WilsonLoops<Impl>::FieldStrength(Ey, _Umu, Tdir, Ydir);
-  WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
-
-  // Compute the Clover Operator acting on Colour and Spin
-  // multiply here by the clover coefficients for the anisotropy
-  double t3 = usecond();
-  TmpOriginal  = Helpers::fillCloverYZ(Bx) * csw_r;
-  TmpOriginal += Helpers::fillCloverXZ(By) * csw_r;
-  TmpOriginal += Helpers::fillCloverXY(Bz) * csw_r;
-  TmpOriginal += Helpers::fillCloverXT(Ex) * csw_t;
-  TmpOriginal += Helpers::fillCloverYT(Ey) * csw_t;
-  TmpOriginal += Helpers::fillCloverZT(Ez) * csw_t;
-
-  // Instantiate the clover term
-  // - In case of the standard clover the mass term is added
-  // - In case of the exponential clover the clover term is exponentiated
-  double t4 = usecond();
-  CloverHelpers::InstantiateClover(TmpOriginal, TmpInverse, csw_t, this->diag_mass);
-
-  // Convert the data layout of the clover term
-  double t5 = usecond();
-  CompactHelpers::ConvertLayout(TmpOriginal, Diagonal, Triangle);
-
-  // Modify the clover term at the temporal boundaries in case of open boundary conditions
-  double t6 = usecond();
-  if(fixedBoundaries) CompactHelpers::ModifyBoundaries(Diagonal, Triangle, csw_t, cF, this->diag_mass);
-
-  // Invert the Clover term
-  // In case of the exponential clover with (anti-)periodic boundary conditions exp(-Clover) saved
-  // in TmpInverse can be used. In all other cases the clover term has to be explictly inverted.
-  // TODO: For now this inversion is explictly done on the CPU
-  double t7 = usecond();
-  CloverHelpers::InvertClover(TmpInverse, Diagonal, Triangle, DiagonalInv, TriangleInv, fixedBoundaries);
-
-  // Fill the remaining clover fields
-  double t8 = usecond();
-  pickCheckerboard(Even, DiagonalEven,    Diagonal);
-  pickCheckerboard(Even, TriangleEven,    Triangle);
-  pickCheckerboard(Odd,  DiagonalOdd,     Diagonal);
-  pickCheckerboard(Odd,  TriangleOdd,     Triangle);
-  pickCheckerboard(Even, DiagonalInvEven, DiagonalInv);
-  pickCheckerboard(Even, TriangleInvEven, TriangleInv);
-  pickCheckerboard(Odd,  DiagonalInvOdd,  DiagonalInv);
-  pickCheckerboard(Odd,  TriangleInvOdd,  TriangleInv);
-
-  // Report timings
-  double t9 = usecond();
-
-  std::cout << GridLogDebug << "CompactWilsonCloverFermion::ImportGauge timings:" << std::endl;
-  std::cout << GridLogDebug << "WilsonFermion::Importgauge = " << (t1 - t0) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "allocations =                " << (t2 - t1) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "field strength =             " << (t3 - t2) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "fill clover =                " << (t4 - t3) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "instantiate clover =         " << (t5 - t4) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "convert layout =             " << (t6 - t5) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "modify boundaries =          " << (t7 - t6) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "invert clover =              " << (t8 - t7) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "pick cbs =                   " << (t9 - t8) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "total =                      " << (t9 - t0) / 1e6 << std::endl;
-}
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h

@@ -2,13 +2,12 @@
 
     Grid physics library, www.github.com/paboyle/Grid
 
-    Source file: ./lib/qcd/action/fermion/WilsonCloverFermionImplementation.h
+    Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.cc
 
-    Copyright (C) 2017 - 2022
+    Copyright (C) 2017
 
     Author: paboyle <paboyle@ph.ed.ac.uk>
     Author: Guido Cossu <guido.cossu@ed.ac.uk>
-    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -34,48 +33,9 @@
 
 NAMESPACE_BEGIN(Grid);
 
-template<class Impl, class CloverHelpers>
-WilsonCloverFermion<Impl, CloverHelpers>::WilsonCloverFermion(GaugeField&                         _Umu,
-                                               GridCartesian&                      Fgrid,
-                                               GridRedBlackCartesian&              Hgrid,
-                                               const RealD                         _mass,
-                                               const RealD                         _csw_r,
-                                               const RealD                         _csw_t,
-                                               const WilsonAnisotropyCoefficients& clover_anisotropy,
-                                               const ImplParams&                   impl_p)
-  : WilsonFermion<Impl>(_Umu, Fgrid, Hgrid, _mass, impl_p, clover_anisotropy)
-  , CloverTerm(&Fgrid)
-  , CloverTermInv(&Fgrid)
-  , CloverTermEven(&Hgrid)
-  , CloverTermOdd(&Hgrid)
-  , CloverTermInvEven(&Hgrid)
-  , CloverTermInvOdd(&Hgrid)
-  , CloverTermDagEven(&Hgrid)
-  , CloverTermDagOdd(&Hgrid)
-  , CloverTermInvDagEven(&Hgrid)
-  , CloverTermInvDagOdd(&Hgrid) {
-  assert(Nd == 4); // require 4 dimensions
-
-  if(clover_anisotropy.isAnisotropic) {
-    csw_r     = _csw_r * 0.5 / clover_anisotropy.xi_0;
-    diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0);
-  } else {
-    csw_r     = _csw_r * 0.5;
-    diag_mass = 4.0 + _mass;
-  }
-  csw_t = _csw_t * 0.5;
-
-  if(csw_r == 0)
-    std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_r = 0" << std::endl;
-  if(csw_t == 0)
-    std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_t = 0" << std::endl;
-
-  ImportGauge(_Umu);
-}
-
 // *NOT* EO
-template<class Impl, class CloverHelpers>
+template <class Impl>
-void WilsonCloverFermion<Impl, CloverHelpers>::M(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
 {
   FermionField temp(out.Grid());
 
@@ -89,8 +49,8 @@ void WilsonCloverFermion<Impl, CloverHelpers>::M(const FermionField &in, Fermion
   out += temp;
 }
 
-template<class Impl, class CloverHelpers>
+template <class Impl>
-void WilsonCloverFermion<Impl, CloverHelpers>::Mdag(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
 {
   FermionField temp(out.Grid());
 
@@ -104,16 +64,13 @@ void WilsonCloverFermion<Impl, CloverHelpers>::Mdag(const FermionField &in, Ferm
   out += temp;
 }
 
-template<class Impl, class CloverHelpers>
+template <class Impl>
-void WilsonCloverFermion<Impl, CloverHelpers>::ImportGauge(const GaugeField &_Umu)
+void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
 {
-  double t0 = usecond();
   WilsonFermion<Impl>::ImportGauge(_Umu);
-  double t1 = usecond();
   GridBase *grid = _Umu.Grid();
   typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid);
 
-  double t2 = usecond();
   // Compute the field strength terms mu>nu
   WilsonLoops<Impl>::FieldStrength(Bx, _Umu, Zdir, Ydir);
   WilsonLoops<Impl>::FieldStrength(By, _Umu, Zdir, Xdir);
@@ -122,20 +79,52 @@ void WilsonCloverFermion<Impl, CloverHelpers>::ImportGauge(const GaugeField &_Um
   WilsonLoops<Impl>::FieldStrength(Ey, _Umu, Tdir, Ydir);
   WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
 
-  double t3 = usecond();
   // Compute the Clover Operator acting on Colour and Spin
   // multiply here by the clover coefficients for the anisotropy
-  CloverTerm  = Helpers::fillCloverYZ(Bx) * csw_r;
+  CloverTerm  = fillCloverYZ(Bx) * csw_r;
-  CloverTerm += Helpers::fillCloverXZ(By) * csw_r;
+  CloverTerm += fillCloverXZ(By) * csw_r;
-  CloverTerm += Helpers::fillCloverXY(Bz) * csw_r;
+  CloverTerm += fillCloverXY(Bz) * csw_r;
-  CloverTerm += Helpers::fillCloverXT(Ex) * csw_t;
+  CloverTerm += fillCloverXT(Ex) * csw_t;
-  CloverTerm += Helpers::fillCloverYT(Ey) * csw_t;
+  CloverTerm += fillCloverYT(Ey) * csw_t;
-  CloverTerm += Helpers::fillCloverZT(Ez) * csw_t;
+  CloverTerm += fillCloverZT(Ez) * csw_t;
-   
+  CloverTerm += diag_mass;
-  double t4 = usecond();
+
-  CloverHelpers::Instantiate(CloverTerm, CloverTermInv, csw_t, this->diag_mass);
+  int lvol = _Umu.Grid()->lSites();
+  int DimRep = Impl::Dimension;
+
+  {
+    autoView(CTv,CloverTerm,CpuRead);
+    autoView(CTIv,CloverTermInv,CpuWrite);
+    thread_for(site, lvol, {
+      Coordinate lcoor;
+      grid->LocalIndexToLocalCoor(site, lcoor);
+      Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
+      Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
+      typename SiteCloverType::scalar_object Qx = Zero(), Qxinv = Zero();
+      peekLocalSite(Qx, CTv, lcoor);
+      //if (csw!=0){
+      for (int j = 0; j < Ns; j++)
+	for (int k = 0; k < Ns; k++)
+	  for (int a = 0; a < DimRep; a++)
+	    for (int b = 0; b < DimRep; b++){
+	      auto zz =  Qx()(j, k)(a, b);
+	      EigenCloverOp(a + j * DimRep, b + k * DimRep) = std::complex<double>(zz);
+	    }
+      //   if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl;
+      
+      EigenInvCloverOp = EigenCloverOp.inverse();
+      //std::cout << EigenInvCloverOp << std::endl;
+      for (int j = 0; j < Ns; j++)
+	for (int k = 0; k < Ns; k++)
+	  for (int a = 0; a < DimRep; a++)
+	    for (int b = 0; b < DimRep; b++)
+	      Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep);
+      //    if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl;
+      //  }
+      pokeLocalSite(Qxinv, CTIv, lcoor);
+    });
+  }
 
-  double t5 = usecond();
   // Separate the even and odd parts
   pickCheckerboard(Even, CloverTermEven, CloverTerm);
   pickCheckerboard(Odd, CloverTermOdd, CloverTerm);
@@ -148,47 +137,37 @@ void WilsonCloverFermion<Impl, CloverHelpers>::ImportGauge(const GaugeField &_Um
 
   pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv));
   pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv));
-  double t6 = usecond();
-
-  std::cout << GridLogDebug << "WilsonCloverFermion::ImportGauge timings:" << std::endl;
-  std::cout << GridLogDebug << "WilsonFermion::Importgauge = " << (t1 - t0) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "allocations =                " << (t2 - t1) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "field strength =             " << (t3 - t2) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "fill clover =                " << (t4 - t3) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "instantiation =              " << (t5 - t4) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "pick cbs =                   " << (t6 - t5) / 1e6 << std::endl;
-  std::cout << GridLogDebug << "total =                      " << (t6 - t0) / 1e6 << std::endl;
 }
 
-template<class Impl, class CloverHelpers>
+template <class Impl>
-void WilsonCloverFermion<Impl, CloverHelpers>::Mooee(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl>::Mooee(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerNo, InverseNo);
 }
 
-template<class Impl, class CloverHelpers>
+template <class Impl>
-void WilsonCloverFermion<Impl, CloverHelpers>::MooeeDag(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerYes, InverseNo);
 }
 
-template<class Impl, class CloverHelpers>
+template <class Impl>
-void WilsonCloverFermion<Impl, CloverHelpers>::MooeeInv(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerNo, InverseYes);
 }
 
-template<class Impl, class CloverHelpers>
+template <class Impl>
-void WilsonCloverFermion<Impl, CloverHelpers>::MooeeInvDag(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerYes, InverseYes);
 }
 
-template<class Impl, class CloverHelpers>
+template <class Impl>
-void WilsonCloverFermion<Impl, CloverHelpers>::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv)
+void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv)
 {
   out.Checkerboard() = in.Checkerboard();
-  CloverField *Clover;
+  CloverFieldType *Clover;
   assert(in.Checkerboard() == Odd || in.Checkerboard() == Even);
 
   if (dag)
@@ -203,12 +182,12 @@ void WilsonCloverFermion<Impl, CloverHelpers>::MooeeInternal(const FermionField
       {
         Clover = (inv) ? &CloverTermInvDagEven : &CloverTermDagEven;
       }
-      Helpers::multCloverField(out, *Clover, in);
+      out = *Clover * in;
     }
     else
     {
       Clover = (inv) ? &CloverTermInv : &CloverTerm;
-      Helpers::multCloverField(out, *Clover, in); // don't bother with adj, hermitian anyway
+      out = adj(*Clover) * in;
     }
   }
   else
@@ -226,109 +205,29 @@ void WilsonCloverFermion<Impl, CloverHelpers>::MooeeInternal(const FermionField
         //  std::cout << "Calling clover term Even" << std::endl;
         Clover = (inv) ? &CloverTermInvEven : &CloverTermEven;
       }
-      Helpers::multCloverField(out, *Clover, in);
+      out = *Clover * in;
       //  std::cout << GridLogMessage << "*Clover.Checkerboard() "  << (*Clover).Checkerboard() << std::endl;
     }
     else
     {
       Clover = (inv) ? &CloverTermInv : &CloverTerm;
-      Helpers::multCloverField(out, *Clover, in);
+      out = *Clover * in;
     }
   }
+
 } // MooeeInternal
 
-// Derivative parts unpreconditioned pseudofermions
-template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl, CloverHelpers>::MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag)
-{
-  conformable(X.Grid(), Y.Grid());
-  conformable(X.Grid(), force.Grid());
-  GaugeLinkField force_mu(force.Grid()), lambda(force.Grid());
-  GaugeField clover_force(force.Grid());
-  PropagatorField Lambda(force.Grid());
-
-  // Guido: Here we are hitting some performance issues:
-  // need to extract the components of the DoubledGaugeField
-  // for each call
-  // Possible solution
-  // Create a vector object to store them? (cons: wasting space)
-  std::vector<GaugeLinkField> U(Nd, this->Umu.Grid());
-
-  Impl::extractLinkField(U, this->Umu);
-
-  force = Zero();
-  // Derivative of the Wilson hopping term
-  this->DhopDeriv(force, X, Y, dag);
-
-  ///////////////////////////////////////////////////////////
-  // Clover term derivative
-  ///////////////////////////////////////////////////////////
-  Impl::outerProductImpl(Lambda, X, Y);
-  //std::cout << "Lambda:" << Lambda << std::endl;
-
-  Gamma::Algebra sigma[] = {
-      Gamma::Algebra::SigmaXY,
-      Gamma::Algebra::SigmaXZ,
-      Gamma::Algebra::SigmaXT,
-      Gamma::Algebra::MinusSigmaXY,
-      Gamma::Algebra::SigmaYZ,
-      Gamma::Algebra::SigmaYT,
-      Gamma::Algebra::MinusSigmaXZ,
-      Gamma::Algebra::MinusSigmaYZ,
-      Gamma::Algebra::SigmaZT,
-      Gamma::Algebra::MinusSigmaXT,
-      Gamma::Algebra::MinusSigmaYT,
-      Gamma::Algebra::MinusSigmaZT};
-
-  /*
-    sigma_{\mu \nu}=
-    | 0         sigma[0]  sigma[1]  sigma[2] |
-    | sigma[3]    0       sigma[4]  sigma[5] |
-    | sigma[6]  sigma[7]     0      sigma[8] |
-    | sigma[9]  sigma[10] sigma[11]   0      |
-  */
-
-  int count = 0;
-  clover_force = Zero();
-  for (int mu = 0; mu < 4; mu++)
-  {
-    force_mu = Zero();
-    for (int nu = 0; nu < 4; nu++)
-    {
-      if (mu == nu)
-      continue;
-
-      RealD factor;
-      if (nu == 4 || mu == 4)
-      {
-        factor = 2.0 * csw_t;
-      }
-      else
-      {
-        factor = 2.0 * csw_r;
-      }
-      PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
-      Impl::TraceSpinImpl(lambda, Slambda);                   // traceSpin ok
-      force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu);                   // checked
-      count++;
-    }
-
-    pokeLorentz(clover_force, U[mu] * force_mu, mu);
-  }
-  //clover_force *= csw;
-  force += clover_force;
-}
 
 // Derivative parts
-template<class Impl, class CloverHelpers>
+template <class Impl>
-void WilsonCloverFermion<Impl, CloverHelpers>::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag)
+void WilsonCloverFermion<Impl>::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag)
 {
   assert(0);
 }
 
 // Derivative parts
-template<class Impl, class CloverHelpers>
+template <class Impl>
-void WilsonCloverFermion<Impl, CloverHelpers>::MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
+void WilsonCloverFermion<Impl>::MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
 {
   assert(0); // not implemented yet
 }

Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h

@@ -4,13 +4,12 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: ./lib/qcd/action/fermion/WilsonFermion.cc
 
-Copyright (C) 2022
+Copyright (C) 2015
 
 Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
 Author: paboyle <paboyle@ph.ed.ac.uk>
-Author: Fabian Joswig <fabian.joswig@ed.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -600,47 +599,11 @@ void WilsonFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
                                                    Current curr_type,
                                                    unsigned int mu)
 {
-  if(curr_type != Current::Vector)
-  {
-    std::cout << GridLogError << "Only the conserved vector current is implemented so far." << std::endl;
-    exit(1);
-  }
-
   Gamma g5(Gamma::Algebra::Gamma5);
   conformable(_grid, q_in_1.Grid());
   conformable(_grid, q_in_2.Grid());
   conformable(_grid, q_out.Grid());
-  auto UGrid= this->GaugeGrid();
+  assert(0);
-
-  PropagatorField tmp_shifted(UGrid);
-  PropagatorField g5Lg5(UGrid);
-  PropagatorField R(UGrid);
-  PropagatorField gmuR(UGrid);
-
-    Gamma::Algebra Gmu [] = {
-    Gamma::Algebra::GammaX,
-    Gamma::Algebra::GammaY,
-    Gamma::Algebra::GammaZ,
-    Gamma::Algebra::GammaT,
-  };
-  Gamma gmu=Gamma(Gmu[mu]);
-
-  g5Lg5=g5*q_in_1*g5;
-  tmp_shifted=Cshift(q_in_2,mu,1);
-  Impl::multLinkField(R,this->Umu,tmp_shifted,mu);
-  gmuR=gmu*R;
-
-  q_out=adj(g5Lg5)*R;
-  q_out-=adj(g5Lg5)*gmuR;
-
-  tmp_shifted=Cshift(q_in_1,mu,1);
-  Impl::multLinkField(g5Lg5,this->Umu,tmp_shifted,mu);
-  g5Lg5=g5*g5Lg5*g5;
-  R=q_in_2;
-  gmuR=gmu*R;
-
-  q_out-=adj(g5Lg5)*R;
-  q_out-=adj(g5Lg5)*gmuR;
 }
 
 
@@ -654,51 +617,9 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
                                               unsigned int tmax,
 					      ComplexField &lattice_cmplx)
 {
-  if(curr_type != Current::Vector)
-  {
-    std::cout << GridLogError << "Only the conserved vector current is implemented so far." << std::endl;
-    exit(1);
-  }
-
-  int tshift = (mu == Nd-1) ? 1 : 0;
-  unsigned int LLt    = GridDefaultLatt()[Tp];
   conformable(_grid, q_in.Grid());
   conformable(_grid, q_out.Grid());
-  auto UGrid= this->GaugeGrid();
+  assert(0);
-
-  PropagatorField tmp(UGrid);
-  PropagatorField Utmp(UGrid);
-  PropagatorField L(UGrid);
-  PropagatorField zz (UGrid);
-  zz=Zero();
-  LatticeInteger lcoor(UGrid); LatticeCoordinate(lcoor,Nd-1);
-
-    Gamma::Algebra Gmu [] = {
-    Gamma::Algebra::GammaX,
-    Gamma::Algebra::GammaY,
-    Gamma::Algebra::GammaZ,
-    Gamma::Algebra::GammaT,
-  };
-  Gamma gmu=Gamma(Gmu[mu]);
-
-  tmp = Cshift(q_in,mu,1);
-  Impl::multLinkField(Utmp,this->Umu,tmp,mu);
-  tmp = ( Utmp*lattice_cmplx - gmu*Utmp*lattice_cmplx ); // Forward hop
-  tmp = where((lcoor>=tmin),tmp,zz); // Mask the time
-  q_out = where((lcoor<=tmax),tmp,zz); // Position of current complicated
-
-  tmp = q_in *lattice_cmplx;
-  tmp = Cshift(tmp,mu,-1);
-  Impl::multLinkField(Utmp,this->Umu,tmp,mu+Nd); // Adjoint link
-  tmp = -( Utmp + gmu*Utmp );
-  // Mask the time
-  if (tmax == LLt - 1 && tshift == 1){ // quick fix to include timeslice 0 if tmax + tshift is over the last timeslice
-    unsigned int t0 = 0;
-    tmp = where(((lcoor==t0) || (lcoor>=tmin+tshift)),tmp,zz);
-  } else {
-    tmp = where((lcoor>=tmin+tshift),tmp,zz);
-  }
-  q_out+= where((lcoor<=tmax+tshift),tmp,zz); // Position of current complicated
 }
 
 NAMESPACE_END(Grid);

Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h

@@ -498,7 +498,6 @@ void WilsonKernels<Impl>::DhopKernel(int Opt,StencilImpl &st,  DoubledGaugeField
 #ifndef GRID_CUDA
      if (Opt == WilsonKernelsStatic::OptInlineAsm  ) {  ASM_CALL(AsmDhopSiteDag);     return;}
 #endif
-     acceleratorFenceComputeStream();
    } else if( interior ) {
      if (Opt == WilsonKernelsStatic::OptGeneric    ) { KERNEL_CALL(GenericDhopSiteDagInt); return;}
      if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDagInt);    return;}
@@ -506,13 +505,11 @@ void WilsonKernels<Impl>::DhopKernel(int Opt,StencilImpl &st,  DoubledGaugeField
      if (Opt == WilsonKernelsStatic::OptInlineAsm  ) {  ASM_CALL(AsmDhopSiteDagInt);     return;}
 #endif
    } else if( exterior ) {
-     acceleratorFenceComputeStream();
      if (Opt == WilsonKernelsStatic::OptGeneric    ) { KERNEL_CALL(GenericDhopSiteDagExt); return;}
      if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDagExt);    return;}
 #ifndef GRID_CUDA
      if (Opt == WilsonKernelsStatic::OptInlineAsm  ) {  ASM_CALL(AsmDhopSiteDagExt);     return;}
 #endif
-     acceleratorFenceComputeStream();
    }
    assert(0 && " Kernel optimisation case not covered ");
   }

Grid/qcd/action/fermion/instantiation/CompactWilsonCloverFermionInstantiation.cc.master

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./lib/ qcd/action/fermion/instantiation/CompactWilsonCloverFermionInstantiation.cc.master
-
-    Copyright (C) 2017 - 2022
-
-    Author: paboyle <paboyle@ph.ed.ac.uk>
-    Author: Guido Cossu <guido.cossu@ed.ac.uk>
-    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
-    Author: Mattia Bruno <mattia.bruno@cern.ch>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-/*  END LEGAL */
-
-#include <Grid/Grid.h>
-#include <Grid/qcd/spin/Dirac.h>
-#include <Grid/qcd/action/fermion/CompactWilsonCloverFermion.h>
-#include <Grid/qcd/action/fermion/implementation/CompactWilsonCloverFermionImplementation.h>
-#include <Grid/qcd/action/fermion/CloverHelpers.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class CompactWilsonCloverFermion<IMPLEMENTATION, CompactCloverHelpers<IMPLEMENTATION>>; 
-template class CompactWilsonCloverFermion<IMPLEMENTATION, CompactExpCloverHelpers<IMPLEMENTATION>>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/WilsonAdjImplD/WilsonKernelsInstantiationWilsonAdjImplD.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/WilsonAdjImplD/WilsonKernelsInstantiationWilsonAdjImplD.cc

@@ -0,0 +1,51 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
+
+Copyright (C) 2015, 2020
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Nils Meyer <nils.meyer@ur.de> Regensburg University
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
+
+#ifndef AVX512
+#ifndef QPX
+#ifndef A64FX
+#ifndef A64FXFIXEDSIZE
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
+#endif
+#endif
+#endif
+#endif
+
+NAMESPACE_BEGIN(Grid);
+
+#include "impl.h"
+template class WilsonKernels<IMPLEMENTATION>;
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/WilsonAdjImplF/WilsonKernelsInstantiationWilsonAdjImplF.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/WilsonAdjImplF/WilsonKernelsInstantiationWilsonAdjImplF.cc

@@ -0,0 +1,51 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
+
+Copyright (C) 2015, 2020
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Nils Meyer <nils.meyer@ur.de> Regensburg University
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
+
+#ifndef AVX512
+#ifndef QPX
+#ifndef A64FX
+#ifndef A64FXFIXEDSIZE
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
+#endif
+#endif
+#endif
+#endif
+
+NAMESPACE_BEGIN(Grid);
+
+#include "impl.h"
+template class WilsonKernels<IMPLEMENTATION>;
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/WilsonCloverFermionInstantiation.cc.master

@@ -8,8 +8,7 @@
 
     Author: paboyle <paboyle@ph.ed.ac.uk>
     Author: Guido Cossu <guido.cossu@ed.ac.uk>
-    Author: Mattia Bruno <mattia.bruno@cern.ch>
+
-    
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
     the Free Software Foundation; either version 2 of the License, or
@@ -32,12 +31,10 @@
 #include <Grid/qcd/spin/Dirac.h>
 #include <Grid/qcd/action/fermion/WilsonCloverFermion.h>
 #include <Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h>
-#include <Grid/qcd/action/fermion/CloverHelpers.h>
 
 NAMESPACE_BEGIN(Grid);
 
 #include "impl.h"
-template class WilsonCloverFermion<IMPLEMENTATION, CloverHelpers<IMPLEMENTATION>>; 
+template class WilsonCloverFermion<IMPLEMENTATION>; 
-template class WilsonCloverFermion<IMPLEMENTATION, ExpCloverHelpers<IMPLEMENTATION>>; 
 
 NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/WilsonImplD/CompactWilsonCloverFermionInstantiationWilsonImplD.cc

@@ -1 +0,0 @@
-../CompactWilsonCloverFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/WilsonImplD/WilsonKernelsInstantiationWilsonImplD.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/WilsonImplD/WilsonKernelsInstantiationWilsonImplD.cc

@@ -0,0 +1,51 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
+
+Copyright (C) 2015, 2020
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Nils Meyer <nils.meyer@ur.de> Regensburg University
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
+
+#ifndef AVX512
+#ifndef QPX
+#ifndef A64FX
+#ifndef A64FXFIXEDSIZE
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
+#endif
+#endif
+#endif
+#endif
+
+NAMESPACE_BEGIN(Grid);
+
+#include "impl.h"
+template class WilsonKernels<IMPLEMENTATION>;
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/WilsonImplF/CompactWilsonCloverFermionInstantiationWilsonImplF.cc

@@ -1 +0,0 @@
-../CompactWilsonCloverFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/WilsonImplF/WilsonKernelsInstantiationWilsonImplF.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/WilsonImplF/WilsonKernelsInstantiationWilsonImplF.cc

@@ -0,0 +1,51 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
+
+Copyright (C) 2015, 2020
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Nils Meyer <nils.meyer@ur.de> Regensburg University
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
+
+#ifndef AVX512
+#ifndef QPX
+#ifndef A64FX
+#ifndef A64FXFIXEDSIZE
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
+#endif
+#endif
+#endif
+#endif
+
+NAMESPACE_BEGIN(Grid);
+
+#include "impl.h"
+template class WilsonKernels<IMPLEMENTATION>;
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplD.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplD.cc

@@ -0,0 +1,51 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
+
+Copyright (C) 2015, 2020
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Nils Meyer <nils.meyer@ur.de> Regensburg University
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
+
+#ifndef AVX512
+#ifndef QPX
+#ifndef A64FX
+#ifndef A64FXFIXEDSIZE
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
+#endif
+#endif
+#endif
+#endif
+
+NAMESPACE_BEGIN(Grid);
+
+#include "impl.h"
+template class WilsonKernels<IMPLEMENTATION>;
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplF.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplF.cc

@@ -0,0 +1,51 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
+
+Copyright (C) 2015, 2020
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Nils Meyer <nils.meyer@ur.de> Regensburg University
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
+
+#ifndef AVX512
+#ifndef QPX
+#ifndef A64FX
+#ifndef A64FXFIXEDSIZE
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
+#endif
+#endif
+#endif
+#endif
+
+NAMESPACE_BEGIN(Grid);
+
+#include "impl.h"
+template class WilsonKernels<IMPLEMENTATION>;
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplD.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplD.cc

@@ -0,0 +1,51 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
+
+Copyright (C) 2015, 2020
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Nils Meyer <nils.meyer@ur.de> Regensburg University
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
+
+#ifndef AVX512
+#ifndef QPX
+#ifndef A64FX
+#ifndef A64FXFIXEDSIZE
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
+#endif
+#endif
+#endif
+#endif
+
+NAMESPACE_BEGIN(Grid);
+
+#include "impl.h"
+template class WilsonKernels<IMPLEMENTATION>;
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplF.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplF.cc

@@ -0,0 +1,51 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
+
+Copyright (C) 2015, 2020
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Nils Meyer <nils.meyer@ur.de> Regensburg University
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
+
+#ifndef AVX512
+#ifndef QPX
+#ifndef A64FX
+#ifndef A64FXFIXEDSIZE
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
+#endif
+#endif
+#endif
+#endif
+
+NAMESPACE_BEGIN(Grid);
+
+#include "impl.h"
+template class WilsonKernels<IMPLEMENTATION>;
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/ZWilsonImplD/WilsonKernelsInstantiationZWilsonImplD.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/ZWilsonImplD/WilsonKernelsInstantiationZWilsonImplD.cc

@@ -0,0 +1,51 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
+
+Copyright (C) 2015, 2020
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Nils Meyer <nils.meyer@ur.de> Regensburg University
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
+
+#ifndef AVX512
+#ifndef QPX
+#ifndef A64FX
+#ifndef A64FXFIXEDSIZE
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
+#endif
+#endif
+#endif
+#endif
+
+NAMESPACE_BEGIN(Grid);
+
+#include "impl.h"
+template class WilsonKernels<IMPLEMENTATION>;
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/ZWilsonImplF/WilsonKernelsInstantiationZWilsonImplF.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/ZWilsonImplF/WilsonKernelsInstantiationZWilsonImplF.cc

@@ -0,0 +1,51 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
+
+Copyright (C) 2015, 2020
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Nils Meyer <nils.meyer@ur.de> Regensburg University
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
+
+#ifndef AVX512
+#ifndef QPX
+#ifndef A64FX
+#ifndef A64FXFIXEDSIZE
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
+#endif
+#endif
+#endif
+#endif
+
+NAMESPACE_BEGIN(Grid);
+
+#include "impl.h"
+template class WilsonKernels<IMPLEMENTATION>;
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/generate_instantiations.sh

@@ -18,10 +18,6 @@ WILSON_IMPL_LIST=" \
 	   GparityWilsonImplF \
 	   GparityWilsonImplD "
 
-COMPACT_WILSON_IMPL_LIST=" \
-	   WilsonImplF \
-	   WilsonImplD "
-
 DWF_IMPL_LIST=" \
 	   WilsonImplF \
 	   WilsonImplD \
@@ -50,17 +46,7 @@ for impl in $WILSON_IMPL_LIST
 do
 for f in $CC_LIST
 do
-  ln -f -s ../$f.cc.master $impl/$f$impl.cc
+  ln -f -s ../$f.cc.master $impl/$f$impl.cc 
-done
-done
-
-CC_LIST="CompactWilsonCloverFermionInstantiation"
-
-for impl in $COMPACT_WILSON_IMPL_LIST
-do
-for f in $CC_LIST
-do
-  ln -f -s ../$f.cc.master $impl/$f$impl.cc
 done
 done
 
@@ -77,14 +63,14 @@ for impl in $DWF_IMPL_LIST $GDWF_IMPL_LIST
 do
 for f in $CC_LIST
 do
-  ln -f -s ../$f.cc.master $impl/$f$impl.cc
+  ln -f -s ../$f.cc.master $impl/$f$impl.cc 
 done
 done
 
 # overwrite the .cc file in Gparity directories
 for impl in $GDWF_IMPL_LIST
 do
-  ln -f -s ../WilsonKernelsInstantiationGparity.cc.master $impl/WilsonKernelsInstantiation$impl.cc
+  ln -f -s ../WilsonKernelsInstantiationGparity.cc.master $impl/WilsonKernelsInstantiation$impl.cc 
 done
 
 
@@ -98,7 +84,7 @@ for impl in $STAG_IMPL_LIST
 do
 for f in $CC_LIST
 do
-  ln -f -s ../$f.cc.master $impl/$f$impl.cc
+  ln -f -s ../$f.cc.master $impl/$f$impl.cc 
 done
 done
 

Grid/qcd/action/gauge/Photon.h

@@ -49,7 +49,7 @@ NAMESPACE_BEGIN(Grid);
     
     typedef Lattice<SiteLink>  LinkField;
     typedef Lattice<SiteField> Field;
-    typedef LinkField          ComplexField;
+    typedef Field              ComplexField;
   };
   
   typedef QedGImpl<vComplex> QedGImplR;

Grid/qcd/utils/GaugeFix.h

@@ -55,12 +55,12 @@ public:
     }
   }  
 
-  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) {
+  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
     GridBase *grid = Umu.Grid();
     GaugeMat xform(grid);
-    SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog,err_on_no_converge);
+    SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog);
   }
-  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) {
+  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) {
 
     GridBase *grid = Umu.Grid();
 
@@ -122,8 +122,6 @@ public:
 
       }
     }
-    std::cout << GridLogError << "Gauge fixing did not converge in " << maxiter << " iterations." << std::endl;
-    if (err_on_no_converge) assert(0);
   };
   static Real SteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu,int orthog) {
     GridBase *grid = U[0].Grid();

Grid/qcd/utils/WilsonLoops.h

@@ -125,6 +125,7 @@ public:
     return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME
   }
 
+
   //////////////////////////////////////////////////
   // average over all x,y,z the temporal loop
   //////////////////////////////////////////////////
@@ -164,7 +165,7 @@ public:
 
     double vol = Umu.Grid()->gSites();
 
-    return p.real() / vol / (4.0 * Nc ) ;
+    return p.real() / vol / 4.0 / 3.0;
   };
 
   //////////////////////////////////////////////////

Grid/serialisation/JSON_IO.cc

@@ -26,7 +26,7 @@
     *************************************************************************************/
     /*  END LEGAL */
 #include <Grid/Grid.h>
-#ifndef GRID_HIP
+#if (!defined(GRID_CUDA)) && (!defined(GRID_HIP))
 
 NAMESPACE_BEGIN(Grid);
 
@@ -82,7 +82,7 @@ void JSONWriter::writeDefault(const std::string &s,	const std::string &x)
   if (s.size())
     ss_ << "\""<< s << "\" : \"" << os.str() << "\" ," ;
   else
-    ss_ << "\""<< os.str() << "\" ," ;
+    ss_ << os.str() << " ," ;
 }
 
 // Reader implementation ///////////////////////////////////////////////////////

Grid/serialisation/JSON_IO.h

@@ -54,7 +54,7 @@ namespace Grid
     void pop(void);
     template <typename U>
     void writeDefault(const std::string &s, const U &x);
-#if defined(GRID_CUDA) || defined(GRID_HIP)
+#ifdef __NVCC__
     void writeDefault(const std::string &s, const Grid::ComplexD &x) 
     { 
       std::complex<double> z(real(x),imag(x));
@@ -101,7 +101,7 @@ namespace Grid
     void readDefault(const std::string &s, std::vector<U> &output);
     template <typename U, typename P>
     void readDefault(const std::string &s, std::pair<U,P> &output);
-#if defined(GRID_CUDA) || defined(GRID_HIP)
+#ifdef __NVCC__
     void readDefault(const std::string &s, ComplexD &output)
     { 
       std::complex<double> z;

Grid/serialisation/Serialisation.h

@@ -36,7 +36,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include "BinaryIO.h"
 #include "TextIO.h"
 #include "XmlIO.h"
-#ifndef GRID_HIP
+#if (!defined(GRID_CUDA)) && (!defined(GRID_HIP))
 #include "JSON_IO.h"
 #endif
 

Grid/stencil/Stencil.h

@@ -80,14 +80,11 @@ void Gather_plane_simple_table (commVector<std::pair<int,int> >& table,const Lat
 ///////////////////////////////////////////////////////////////////
 template<class cobj,class vobj,class compressor>
 void Gather_plane_exchange_table(const Lattice<vobj> &rhs,
-				 commVector<cobj *> pointers,
+				 commVector<cobj *> pointers,int dimension,int plane,int cbmask,compressor &compress,int type) __attribute__((noinline));
-				 int dimension,int plane,
-				 int cbmask,compressor &compress,int type) __attribute__((noinline));
 
 template<class cobj,class vobj,class compressor>
-void Gather_plane_exchange_table(commVector<std::pair<int,int> >& table,
+void Gather_plane_exchange_table(commVector<std::pair<int,int> >& table,const Lattice<vobj> &rhs,
-				 const Lattice<vobj> &rhs,
+				 Vector<cobj *> pointers,int dimension,int plane,int cbmask,
-				 std::vector<cobj *> &pointers,int dimension,int plane,int cbmask,
 				 compressor &compress,int type)
 {
   assert( (table.size()&0x1)==0);
@@ -95,15 +92,14 @@ void Gather_plane_exchange_table(commVector<std::pair<int,int> >& table,
   int so  = plane*rhs.Grid()->_ostride[dimension]; // base offset for start of plane
 
   auto rhs_v = rhs.View(AcceleratorRead);
-  auto rhs_p = &rhs_v[0];
   auto p0=&pointers[0][0];
   auto p1=&pointers[1][0];
   auto tp=&table[0];
   accelerator_forNB(j, num, vobj::Nsimd(), {
-      compress.CompressExchange(p0,p1, rhs_p, j,
+      compress.CompressExchange(p0,p1, &rhs_v[0], j,
-				so+tp[2*j  ].second,
+			      so+tp[2*j  ].second,
-				so+tp[2*j+1].second,
+			      so+tp[2*j+1].second,
-				type);
+			      type);
   });
   rhs_v.ViewClose();
 }
@@ -234,8 +230,8 @@ public:
   };
   struct Merge {
     cobj * mpointer;
-    //    std::vector<scalar_object *> rpointers;
+    Vector<scalar_object *> rpointers;
-    std::vector<cobj *> vpointers;
+    Vector<cobj *> vpointers;
     Integer buffer_size;
     Integer type;
   };
@@ -410,7 +406,6 @@ public:
       comms_bytes+=bytes;
       shm_bytes  +=2*Packets[i].bytes-bytes;
     }
-    _grid->StencilBarrier();// Synch shared memory on a single nodes
   }
 
   void CommunicateComplete(std::vector<std::vector<CommsRequest_t> > &reqs)
@@ -425,7 +420,7 @@ public:
   ////////////////////////////////////////////////////////////////////////
   void Communicate(void)
   {
-    if ( 0 ){
+    if ( CartesianCommunicator::CommunicatorPolicy == CartesianCommunicator::CommunicatorPolicySequential ){
       thread_region {
 	// must be called in parallel region
 	int mythread  = thread_num();
@@ -574,7 +569,7 @@ public:
     d.buffer_size = buffer_size;
     dv.push_back(d);
   }
-  void AddMerge(cobj *merge_p,std::vector<cobj *> &rpointers,Integer buffer_size,Integer type,std::vector<Merge> &mv) {
+  void AddMerge(cobj *merge_p,Vector<cobj *> &rpointers,Integer buffer_size,Integer type,std::vector<Merge> &mv) {
     Merge m;
     m.type     = type;
     m.mpointer = merge_p;
@@ -587,7 +582,6 @@ public:
   }
   template<class decompressor>  void CommsMergeSHM(decompressor decompress) {
     mpi3synctime-=usecond();
-    accelerator_barrier();
     _grid->StencilBarrier();// Synch shared memory on a single nodes
     mpi3synctime+=usecond();
     shmmergetime-=usecond();
@@ -1120,8 +1114,8 @@ public:
     int bytes = (reduced_buffer_size*datum_bytes)/simd_layout;
     assert(bytes*simd_layout == reduced_buffer_size*datum_bytes);
 
-    std::vector<cobj *> rpointers(maxl);
+    Vector<cobj *> rpointers(maxl);
-    std::vector<cobj *> spointers(maxl);
+    Vector<cobj *> spointers(maxl);
 
     ///////////////////////////////////////////
     // Work out what to send where

Grid/threads/Accelerator.cc

@@ -195,15 +195,12 @@ void acceleratorInit(void)
 #ifdef GRID_SYCL
 
 cl::sycl::queue *theGridAccelerator;
-cl::sycl::queue *theCopyAccelerator;
 void acceleratorInit(void)
 {
   int nDevices = 1;
   cl::sycl::gpu_selector selector;
   cl::sycl::device selectedDevice { selector };
   theGridAccelerator = new sycl::queue (selectedDevice);
-  //  theCopyAccelerator = new sycl::queue (selectedDevice);
-  theCopyAccelerator = theGridAccelerator; // Should proceed concurrenlty anyway.
 
 #ifdef GRID_SYCL_LEVEL_ZERO_IPC
   zeInit(0);

Grid/threads/Accelerator.h

@@ -247,6 +247,7 @@ inline int  acceleratorIsCommunicable(void *ptr)
 //////////////////////////////////////////////
 // SyCL acceleration
 //////////////////////////////////////////////
+
 #ifdef GRID_SYCL
 NAMESPACE_END(Grid);
 #include <CL/sycl.hpp>
@@ -261,7 +262,6 @@ NAMESPACE_END(Grid);
 NAMESPACE_BEGIN(Grid);
 
 extern cl::sycl::queue *theGridAccelerator;
-extern cl::sycl::queue *theCopyAccelerator;
 
 #ifdef __SYCL_DEVICE_ONLY__
 #define GRID_SIMT
@@ -289,7 +289,7 @@ accelerator_inline int acceleratorSIMTlane(int Nsimd) {
       cgh.parallel_for(					\
       cl::sycl::nd_range<3>(global,local), \
       [=] (cl::sycl::nd_item<3> item) /*mutable*/     \
-      [[intel::reqd_sub_group_size(16)]]	      \
+      [[intel::reqd_sub_group_size(8)]]	      \
       {						      \
       auto iter1    = item.get_global_id(0);	      \
       auto iter2    = item.get_global_id(1);	      \
@@ -298,19 +298,19 @@ accelerator_inline int acceleratorSIMTlane(int Nsimd) {
      });	   			              \
     });
 
-#define accelerator_barrier(dummy) { theGridAccelerator->wait(); }
+#define accelerator_barrier(dummy) theGridAccelerator->wait();
 
 inline void *acceleratorAllocShared(size_t bytes){ return malloc_shared(bytes,*theGridAccelerator);};
 inline void *acceleratorAllocDevice(size_t bytes){ return malloc_device(bytes,*theGridAccelerator);};
 inline void acceleratorFreeShared(void *ptr){free(ptr,*theGridAccelerator);};
 inline void acceleratorFreeDevice(void *ptr){free(ptr,*theGridAccelerator);};
-
+inline void acceleratorCopyDeviceToDeviceAsynch(void *from,void *to,size_t bytes)  {
-inline void acceleratorCopySynchronise(void) {  theCopyAccelerator->wait(); }
+  theGridAccelerator->memcpy(to,from,bytes);
-inline void acceleratorCopyDeviceToDeviceAsynch(void *from,void *to,size_t bytes)  {  theCopyAccelerator->memcpy(to,from,bytes);}
+}
-inline void acceleratorCopyToDevice(void *from,void *to,size_t bytes)  { theCopyAccelerator->memcpy(to,from,bytes); theCopyAccelerator->wait();}
+inline void acceleratorCopySynchronise(void) {  theGridAccelerator->wait(); std::cout<<"acceleratorCopySynchronise() wait "<<std::endl; }
-inline void acceleratorCopyFromDevice(void *from,void *to,size_t bytes){ theCopyAccelerator->memcpy(to,from,bytes); theCopyAccelerator->wait();}
+inline void acceleratorCopyToDevice(void *from,void *to,size_t bytes)  { theGridAccelerator->memcpy(to,from,bytes); theGridAccelerator->wait();}
-inline void acceleratorMemSet(void *base,int value,size_t bytes) { theCopyAccelerator->memset(base,value,bytes); theCopyAccelerator->wait();}
+inline void acceleratorCopyFromDevice(void *from,void *to,size_t bytes){ theGridAccelerator->memcpy(to,from,bytes); theGridAccelerator->wait();}
-
+inline void acceleratorMemSet(void *base,int value,size_t bytes) { theGridAccelerator->memset(base,value,bytes); theGridAccelerator->wait();}
 inline int  acceleratorIsCommunicable(void *ptr)
 {
 #if 0
@@ -342,7 +342,7 @@ extern hipStream_t copyStream;
 /*These routines define mapping from thread grid to loop & vector lane indexing */
 accelerator_inline int acceleratorSIMTlane(int Nsimd) {
 #ifdef GRID_SIMT
-  return hipThreadIdx_x; 
+  return hipThreadIdx_z; 
 #else
   return 0;
 #endif
@@ -356,41 +356,19 @@ accelerator_inline int acceleratorSIMTlane(int Nsimd) {
       { __VA_ARGS__;}							\
     };									\
     int nt=acceleratorThreads();					\
-    dim3 hip_threads(nsimd, nt, 1);					 \
+    dim3 hip_threads(nt,1,nsimd);					\
-    dim3 hip_blocks ((num1+nt-1)/nt,num2,1); \
+    dim3 hip_blocks ((num1+nt-1)/nt,num2,1);				\
-    if(hip_threads.x * hip_threads.y * hip_threads.z <= 64){ \
+    hipLaunchKernelGGL(LambdaApply,hip_blocks,hip_threads,		\
-      hipLaunchKernelGGL(LambdaApply64,hip_blocks,hip_threads,		\
+		       0,0,						\
-            0,0,						\
+		       num1,num2,nsimd,lambda);				\
-            num1,num2,nsimd, lambda);				\
-    } else { \
-      hipLaunchKernelGGL(LambdaApply,hip_blocks,hip_threads,		\
-            0,0,						\
-            num1,num2,nsimd, lambda);				\
-    } \
   }
 
-
 template<typename lambda>  __global__
-__launch_bounds__(64,1)
-void LambdaApply64(uint64_t numx, uint64_t numy, uint64_t numz, lambda Lambda)
-{
-  // Following the same scheme as CUDA for now
-  uint64_t x = threadIdx.y + blockDim.y*blockIdx.x;
-  uint64_t y = threadIdx.z + blockDim.z*blockIdx.y;
-  uint64_t z = threadIdx.x;
-  if ( (x < numx) && (y<numy) && (z<numz) ) {
-    Lambda(x,y,z);
-  }
-}
-
-template<typename lambda>  __global__
-__launch_bounds__(1024,1)
 void LambdaApply(uint64_t numx, uint64_t numy, uint64_t numz, lambda Lambda)
 {
-  // Following the same scheme as CUDA for now
+  uint64_t x = hipThreadIdx_x + hipBlockDim_x*hipBlockIdx_x;
-  uint64_t x = threadIdx.y + blockDim.y*blockIdx.x;
+  uint64_t y = hipThreadIdx_y + hipBlockDim_y*hipBlockIdx_y;
-  uint64_t y = threadIdx.z + blockDim.z*blockIdx.y;
+  uint64_t z = hipThreadIdx_z ;//+ hipBlockDim_z*hipBlockIdx_z;
-  uint64_t z = threadIdx.x;
   if ( (x < numx) && (y<numy) && (z<numz) ) {
     Lambda(x,y,z);
   }
@@ -481,10 +459,9 @@ inline void acceleratorCopySynchronise(void) { hipStreamSynchronize(copyStream);
 #define accelerator_for2d(iter1, num1, iter2, num2, nsimd, ... ) thread_for2d(iter1,num1,iter2,num2,{ __VA_ARGS__ });
 
 accelerator_inline int acceleratorSIMTlane(int Nsimd) { return 0; } // CUDA specific
-
+inline void acceleratorCopyToDevice(void *from,void *to,size_t bytes)  { memcpy(to,from,bytes);}
-inline void acceleratorCopyToDevice(void *from,void *to,size_t bytes)  { thread_bcopy(from,to,bytes); }
+inline void acceleratorCopyFromDevice(void *from,void *to,size_t bytes){ memcpy(to,from,bytes);}
-inline void acceleratorCopyFromDevice(void *from,void *to,size_t bytes){ thread_bcopy(from,to,bytes);}
+inline void acceleratorCopyDeviceToDeviceAsynch(void *from,void *to,size_t bytes)  { memcpy(to,from,bytes);}
-inline void acceleratorCopyDeviceToDeviceAsynch(void *from,void *to,size_t bytes)  { thread_bcopy(from,to,bytes);}
 inline void acceleratorCopySynchronise(void) {};
 
 inline int  acceleratorIsCommunicable(void *ptr){ return 1; }
@@ -511,16 +488,7 @@ inline void *acceleratorAllocCpu(size_t bytes){return memalign(GRID_ALLOC_ALIGN,
 inline void acceleratorFreeCpu  (void *ptr){free(ptr);};
 #endif
 
-//////////////////////////////////////////////
-// Fencing needed ONLY for SYCL
-//////////////////////////////////////////////
 
-#ifdef GRID_SYCL
-inline void acceleratorFenceComputeStream(void){ accelerator_barrier();};
-#else
-// Ordering within a stream guaranteed on Nvidia & AMD
-inline void acceleratorFenceComputeStream(void){ };
-#endif
 
 ///////////////////////////////////////////////////
 // Synchronise across local threads for divergence resynch

Grid/threads/Threads.h

@@ -72,20 +72,3 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #define thread_region                                       DO_PRAGMA(omp parallel)
 #define thread_critical                                     DO_PRAGMA(omp critical)
 
-#ifdef GRID_OMP
-inline void thread_bcopy(void *from, void *to,size_t bytes)
-{
-  uint64_t *ufrom = (uint64_t *)from;
-  uint64_t *uto   = (uint64_t *)to;
-  assert(bytes%8==0);
-  uint64_t words=bytes/8;
-  thread_for(w,words,{
-      uto[w] = ufrom[w];
-  });
-}
-#else
-inline void thread_bcopy(void *from, void *to,size_t bytes)
-{
-  bcopy(from,to,bytes);
-}
-#endif

Grid/util/Init.cc

@@ -167,13 +167,6 @@ void GridCmdOptionInt(std::string &str,int & val)
   return;
 }
 
-void GridCmdOptionFloat(std::string &str,float & val)
-{
-  std::stringstream ss(str);
-  ss>>val;
-  return;
-}
-
 
 void GridParseLayout(char **argv,int argc,
 		     Coordinate &latt_c,
@@ -534,7 +527,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

Grid/util/Init.h

@@ -57,7 +57,6 @@ void GridCmdOptionCSL(std::string str,std::vector<std::string> & vec);
 template<class VectorInt>
 void GridCmdOptionIntVector(const std::string &str,VectorInt & vec);
 void GridCmdOptionInt(std::string &str,int & val);
-void GridCmdOptionFloat(std::string &str,float & val);
 
 
 void GridParseLayout(char **argv,int argc,

Grid/util/Sha.h

@@ -27,7 +27,6 @@
     /*  END LEGAL */
 extern "C" {
 #include <openssl/sha.h>
-#include <openssl/evp.h>
 }
 #ifdef USE_IPP
 #include "ipp.h"
@@ -71,8 +70,10 @@ public:
   static inline std::vector<unsigned char> sha256(const void *data,size_t bytes)
   {
     std::vector<unsigned char> hash(SHA256_DIGEST_LENGTH);
-    auto digest = EVP_get_digestbyname("SHA256");
+    SHA256_CTX sha256;
-    EVP_Digest(data, bytes, &hash[0], NULL, digest, NULL);
+    SHA256_Init  (&sha256);
+    SHA256_Update(&sha256, data,bytes);
+    SHA256_Final (&hash[0], &sha256);
     return hash;
   }
   static inline std::vector<int> sha256_seeds(const std::string &s)

README.md

@@ -148,7 +148,7 @@ If you want to build all the tests at once just use `make tests`.
 - `--enable-mkl[=<path>]`: use Intel MKL for FFT (and LAPACK if enabled) routines. A UNIX prefix containing the library can be specified (optional).
 - `--enable-numa`: enable NUMA first touch optimisation
 - `--enable-simd=<code>`: setup Grid for the SIMD target `<code>` (default: `GEN`). A list of possible SIMD targets is detailed in a section below.
-- `--enable-gen-simd-width=<size>`: select the size (in bytes) of the generic SIMD vector type (default: 64 bytes).
+- `--enable-gen-simd-width=<size>`: select the size (in bytes) of the generic SIMD vector type (default: 32 bytes).
 - `--enable-comms=<comm>`: Use `<comm>` for message passing (default: `none`). A list of possible SIMD targets is detailed in a section below.
 - `--enable-rng={sitmo|ranlux48|mt19937}`: choose the RNG (default: `sitmo `).
 - `--disable-timers`: disable system dependent high-resolution timers.

benchmarks/Benchmark_halo.cc

@@ -1,131 +0,0 @@
- /*************************************************************************************
-    Grid physics library, www.github.com/paboyle/Grid
-    Source file: ./benchmarks/Benchmark_dwf.cc
-    Copyright (C) 2015
-
-    Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-    Author: paboyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
-#include <Grid/Grid.h>
-#ifdef GRID_CUDA
-#define CUDA_PROFILE
-#endif
-
-#ifdef CUDA_PROFILE
-#include <cuda_profiler_api.h>
-#endif
-
-using namespace std;
-using namespace Grid;
-
-template<class d>
-struct scal {
-  d internal;
-};
-
-  Gamma::Algebra Gmu [] = {
-    Gamma::Algebra::GammaX,
-    Gamma::Algebra::GammaY,
-    Gamma::Algebra::GammaZ,
-    Gamma::Algebra::GammaT
-  };
-
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  Coordinate latt4= GridDefaultLatt();
-  Coordinate mpi  = GridDefaultMpi();
-  Coordinate simd = GridDefaultSimd(Nd,vComplexF::Nsimd());
-
-  GridLogLayout();
-
-  int Ls=16;
-  for(int i=0;i<argc;i++)
-    if(std::string(argv[i]) == "-Ls"){
-      std::stringstream ss(argv[i+1]); ss >> Ls;
-    }
-
-  
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(latt4,simd ,mpi);
-  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-
-  std::cout << GridLogMessage << "Making s innermost grids"<<std::endl;
-  GridCartesian         * sUGrid   = SpaceTimeGrid::makeFourDimDWFGrid(GridDefaultLatt(),GridDefaultMpi());
-  GridRedBlackCartesian * sUrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(sUGrid);
-  GridCartesian         * sFGrid   = SpaceTimeGrid::makeFiveDimDWFGrid(Ls,UGrid);
-  GridRedBlackCartesian * sFrbGrid = SpaceTimeGrid::makeFiveDimDWFRedBlackGrid(Ls,UGrid);
-
-  std::vector<int> seeds4({1,2,3,4});
-  std::vector<int> seeds5({5,6,7,8});
-
-  std::cout << GridLogMessage << "Initialising 4d RNG" << std::endl;
-  GridParallelRNG          RNG4(UGrid);  RNG4.SeedUniqueString(std::string("The 4D RNG"));
-  std::cout << GridLogMessage << "Initialising 5d RNG" << std::endl;
-  GridParallelRNG          RNG5(FGrid);  RNG5.SeedUniqueString(std::string("The 5D RNG"));
-  std::cout << GridLogMessage << "Initialised RNGs" << std::endl;
-
-  LatticeFermionF src   (FGrid); random(RNG5,src);
-  RealD N2 = 1.0/::sqrt(norm2(src));
-  src = src*N2;
-
-  std::cout << GridLogMessage << "Drawing gauge field" << std::endl;
-  LatticeGaugeFieldF Umu(UGrid);
-  SU<Nc>::HotConfiguration(RNG4,Umu);
-  std::cout << GridLogMessage << "Random gauge initialised " << std::endl;
-
-  RealD mass=0.1;
-  RealD M5  =1.8;
-
-  RealD NP = UGrid->_Nprocessors;
-  RealD NN = UGrid->NodeCount();
-
-  DomainWallFermionF Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-
-  const int ncall = 500;
-  std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
-  std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionF::HaloGatherOpt         "<<std::endl;
-  std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
-  {
-    typename DomainWallFermionF::Compressor compressor(0);
-    FGrid->Barrier();
-    Dw.Stencil.HaloExchangeOptGather(src,compressor);
-    double t0=usecond();
-    for(int i=0;i<ncall;i++){
-      Dw.Stencil.HaloExchangeOptGather(src,compressor);
-    }
-    double t1=usecond();
-    FGrid->Barrier();
-
-    double bytes=0.0;
-    if(mpi[0]) bytes+=latt4[1]*latt4[2]*latt4[3];
-    if(mpi[1]) bytes+=latt4[0]*latt4[2]*latt4[3];
-    if(mpi[2]) bytes+=latt4[0]*latt4[1]*latt4[3];
-    if(mpi[3]) bytes+=latt4[0]*latt4[1]*latt4[2];
-    bytes = bytes * Ls * 8.* (24.+12.)* 2.0;
-
-    std::cout<<GridLogMessage << "Gather us /call =   "<< (t1-t0)/ncall<<std::endl;
-    std::cout<<GridLogMessage << "Gather MBs /call =   "<< bytes*ncall/(t1-t0)<<std::endl;
-
-  }
-
-  Grid_finalize();
-  exit(0);
-}

benchmarks/Benchmark_mooee.cc

@@ -81,8 +81,8 @@ int main (int argc, char ** argv)
     Vector<Coeff_t> diag = Dw.bs;
     Vector<Coeff_t> upper= Dw.cs;
     Vector<Coeff_t> lower= Dw.cs;
-    upper[Ls-1]=-Dw.mass_minus*upper[Ls-1];
+    upper[Ls-1]=-Dw.mass*upper[Ls-1];
-    lower[0]   =-Dw.mass_plus*lower[0];
+    lower[0]   =-Dw.mass*lower[0];
     
     LatticeFermion r_eo(FGrid);
     LatticeFermion src_e (FrbGrid);

benchmarks/Benchmark_wilson_sweep.cc

@@ -44,13 +44,6 @@ void bench_wilson (
 		   double const     volume,
 		   int const           dag );
 
-void bench_wilson_eo (
-       LatticeFermion &    src,
-       LatticeFermion & result,
-       WilsonFermionR &     Dw,
-       double const     volume,
-       int const           dag );
-
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
@@ -117,8 +110,8 @@ int main (int argc, char ** argv)
 	  bench_wilson(src,result,Dw,volume,DaggerYes);
 	  std::cout << "\t";
     // EO
-    bench_wilson_eo(src_o,result_e,Dw,volume,DaggerNo);
+	  bench_wilson(src,result,Dw,volume,DaggerNo);
-    bench_wilson_eo(src_o,result_e,Dw,volume,DaggerYes);
+	  bench_wilson(src,result,Dw,volume,DaggerYes);
 	  std::cout << std::endl;
 	}
     }

configure.ac

@@ -159,7 +159,7 @@ case ${ac_ZMOBIUS} in
 esac
 ############### Nc
 AC_ARG_ENABLE([Nc],
-    [AC_HELP_STRING([--enable-Nc=2|3|4|5], [enable number of colours])],
+    [AC_HELP_STRING([--enable-Nc=2|3|4], [enable number of colours])],
     [ac_Nc=${enable_Nc}], [ac_Nc=3])
 
 case ${ac_Nc} in

systems/Crusher/config-command

@@ -1,12 +0,0 @@
-../../configure --enable-comms=mpi-auto \
---enable-unified=no \
---enable-shm=nvlink \
---enable-accelerator=hip \
---enable-gen-simd-width=64 \
---enable-simd=GPU \
---disable-fermion-reps \
---disable-gparity \
-CXX=hipcc MPICXX=mpicxx \
-CXXFLAGS="-fPIC -I/opt/rocm-4.5.0/include/ -std=c++14 -I${MPICH_DIR}/include " \
- LDFLAGS=" -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa "
-HIPFLAGS = --amdgpu-target=gfx90a

systems/Crusher/dwf.slurm

@@ -1,30 +0,0 @@
-#!/bin/bash
-# Begin LSF Directives
-#SBATCH -A LGT104
-#SBATCH -t 01:00:00
-##SBATCH -U openmpThu
-##SBATCH -p ecp
-#SBATCH -J DWF
-#SBATCH -o DWF.%J
-#SBATCH -e DWF.%J
-#SBATCH -N 1
-#SBATCH -n 1
-#SBATCH --exclusive  
-
-DIR=.
-module list
-#export MPIR_CVAR_GPU_EAGER_DEVICE_MEM=0
-export MPICH_GPU_SUPPORT_ENABLED=1
-export MPICH_SMP_SINGLE_COPY_MODE=XPMEM
-#export MPICH_SMP_SINGLE_COPY_MODE=NONE
-#export MPICH_SMP_SINGLE_COPY_MODE=CMA
-export OMP_NUM_THREADS=1
-
-AT=8
-echo MPICH_SMP_SINGLE_COPY_MODE $MPICH_SMP_SINGLE_COPY_MODE
-
-PARAMS=" --accelerator-threads ${AT} --grid 24.24.24.24 --shm-mpi 0 --mpi 1.1.1.1"
-
-srun --gpus-per-task 1 -n1 ./benchmarks/Benchmark_dwf_fp32 $PARAMS
-
-

systems/Crusher/dwf4.slurm

@@ -1,27 +0,0 @@
-#!/bin/bash
-# Begin LSF Directives
-#SBATCH -A LGT104
-#SBATCH -t 01:00:00
-##SBATCH -U openmpThu
-#SBATCH -J DWF
-#SBATCH -o DWF.%J
-#SBATCH -e DWF.%J
-#SBATCH -N 1
-#SBATCH -n 4
-#SBATCH --exclusive
-
-DIR=.
-module list
-export MPIR_CVAR_GPU_EAGER_DEVICE_MEM=0
-export MPICH_GPU_SUPPORT_ENABLED=1
-#export MPICH_SMP_SINGLE_COPY_MODE=XPMEM
-export MPICH_SMP_SINGLE_COPY_MODE=NONE
-#export MPICH_SMP_SINGLE_COPY_MODE=CMA
-export OMP_NUM_THREADS=4
-
-echo MPICH_SMP_SINGLE_COPY_MODE $MPICH_SMP_SINGLE_COPY_MODE
-PARAMS=" --accelerator-threads 8 --grid 32.32.64.64 --mpi 1.1.2.2 --comms-overlap --shm 2048 --shm-mpi 0"
-
-srun --gpus-per-task 1 -n4 ./mpiwrapper.sh ./benchmarks/Benchmark_dwf_fp32 $PARAMS
-
-

systems/Crusher/dwf8.slurm

@@ -1,27 +0,0 @@
-#!/bin/bash
-# Begin LSF Directives
-#SBATCH -A LGT104
-#SBATCH -t 01:00:00
-##SBATCH -U openmpThu
-#SBATCH -J DWF
-#SBATCH -o DWF.%J
-#SBATCH -e DWF.%J
-#SBATCH -N 1
-#SBATCH -n 8
-#SBATCH --exclusive
-
-DIR=.
-module list
-export MPIR_CVAR_GPU_EAGER_DEVICE_MEM=0
-export MPICH_GPU_SUPPORT_ENABLED=1
-export MPICH_SMP_SINGLE_COPY_MODE=XPMEM
-#export MPICH_SMP_SINGLE_COPY_MODE=NONE
-#export MPICH_SMP_SINGLE_COPY_MODE=CMA
-export OMP_NUM_THREADS=1
-
-echo MPICH_SMP_SINGLE_COPY_MODE $MPICH_SMP_SINGLE_COPY_MODE
-PARAMS=" --accelerator-threads 8 --grid 32.64.64.64 --mpi 1.2.2.2 --comms-overlap --shm 2048 --shm-mpi 0"
-
-srun --gpus-per-task 1 -n8 ./mpiwrapper.sh ./benchmarks/Benchmark_dwf_fp32 $PARAMS
-
-

systems/Crusher/mpiwrapper.sh

@@ -1,12 +0,0 @@
-#!/bin/bash
-
-lrank=$SLURM_LOCALID
-
-export ROCR_VISIBLE_DEVICES=$SLURM_LOCALID
-
-echo "`hostname` - $lrank device=$ROCR_VISIBLE_DEVICES binding=$BINDING"
-
-$*
-
-
-

systems/Crusher/sourceme.sh

@@ -1,5 +0,0 @@
-module load PrgEnv-gnu
-module load rocm/4.5.0
-module load gmp
-module load cray-fftw
-module load craype-accel-amd-gfx90a

systems/PVC/benchmarks/run-1tile.sh

@@ -1,62 +0,0 @@
-#!/bin/sh
-##SBATCH -p PVC-SPR-QZEH 
-##SBATCH -p PVC-ICX-QZNW
-#SBATCH -p QZ1J-ICX-PVC
-##SBATCH -p QZ1J-SPR-PVC-2C
-
-source /nfs/site/home/paboylex/ATS/GridNew/Grid/systems/PVC-nightly/setup.sh
-
-export NT=8
-
-export I_MPI_OFFLOAD=1
-export I_MPI_OFFLOAD_TOPOLIB=level_zero
-export I_MPI_OFFLOAD_DOMAIN_SIZE=-1
-
-# export IGC_EnableLSCFenceUGMBeforeEOT=0
-# export SYCL_PROGRAM_COMPILE_OPTIONS="-ze-opt-large-register-file=False"
-export SYCL_DEVICE_FILTER=gpu,level_zero
-#export IGC_ShaderDumpEnable=1 
-#export IGC_DumpToCurrentDir=1
-export I_MPI_OFFLOAD_CELL=tile
-export EnableImplicitScaling=0
-export EnableWalkerPartition=0
-export ZE_AFFINITY_MASK=0.0
-mpiexec -launcher ssh -n 1 -host localhost  ./Benchmark_dwf_fp32 --mpi 1.1.1.1 --grid 32.32.32.32 --accelerator-threads $NT --comms-sequential --shm-mpi 1 --device-mem 32768
-
-export ZE_AFFINITY_MASK=0
-export I_MPI_OFFLOAD_CELL=device
-export EnableImplicitScaling=1
-export EnableWalkerPartition=1
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-#mpiexec -launcher ssh -n 2 -host localhost  vtune -collect gpu-hotspots -knob gpu-sampling-interval=1 -data-limit=0 -r ./vtune_run4 -- ./wrap.sh ./Benchmark_dwf_fp32 --mpi 2.1.1.1 --grid 64.32.32.32 --accelerator-threads $NT --comms-overlap --shm-mpi 1
-
-#mpiexec  -launcher ssh -n 1 -host localhost ./wrap.sh ./Benchmark_dwf_fp32 --mpi 1.1.1.1 --grid 64.32.32.32 --accelerator-threads $NT --comms-overlap --shm-mpi 1
-
-#mpiexec  -launcher ssh -n 2 -host localhost ./wrap.sh ./Benchmark_dwf_fp32 --mpi 2.1.1.1 --grid 64.32.32.32 --accelerator-threads $NT --comms-sequential --shm-mpi 1
-
-#mpiexec  -launcher ssh -n 2 -host localhost ./wrap.sh ./Benchmark_dwf_fp32 --mpi 2.1.1.1 --grid 64.32.32.32 --accelerator-threads $NT --comms-overlap --shm-mpi 1
-
-#mpiexec  -launcher ssh -n 2 -host localhost ./wrap.sh ./Benchmark_dwf_fp32 --mpi 2.1.1.1 --grid 64.32.32.32 --accelerator-threads $NT --comms-sequential --shm-mpi 0
-
-#mpirun -np 2 ./wrap.sh ./Benchmark_dwf_fp32 --mpi 1.1.1.2 --grid 16.32.32.64 --accelerator-threads $NT --comms-sequential --shm-mpi 0
-#mpirun -np 2 ./wrap.sh ./Benchmark_dwf_fp32 --mpi 1.1.1.2 --grid 32.32.32.64 --accelerator-threads $NT --comms-sequential --shm-mpi 1
-

systems/PVC/benchmarks/run-2tile-mpi.sh

@@ -1,34 +0,0 @@
-#!/bin/bash
-##SBATCH -p PVC-SPR-QZEH 
-##SBATCH -p PVC-ICX-QZNW
-
-#SBATCH -p QZ1J-ICX-PVC
-
-source /nfs/site/home/paboylex/ATS/GridNew/Grid/systems/PVC-nightly/setup.sh
-
-export NT=16
-
-# export IGC_EnableLSCFenceUGMBeforeEOT=0
-# export SYCL_PROGRAM_COMPILE_OPTIONS="-ze-opt-large-register-file=False"
-#export IGC_ShaderDumpEnable=1 
-#export IGC_DumpToCurrentDir=1
-export I_MPI_OFFLOAD=1
-export I_MPI_OFFLOAD_TOPOLIB=level_zero
-export I_MPI_OFFLOAD_DOMAIN_SIZE=-1
-export SYCL_DEVICE_FILTER=gpu,level_zero
-export I_MPI_OFFLOAD_CELL=tile
-export EnableImplicitScaling=0
-export EnableWalkerPartition=0
-export SYCL_PI_LEVEL_ZERO_DEVICE_SCOPE_EVENTS=1
-export SYCL_PI_LEVEL_ZERO_USE_IMMEDIATE_COMMANDLISTS=1
-export SYCL_PI_LEVEL_ZERO_USE_COPY_ENGINE=0
-
-for i in 0 
-do
-mpiexec -launcher ssh -n 2 -host localhost  ./wrap4gpu.sh ./Benchmark_dwf_fp32 --mpi 1.1.1.2 --grid 32.32.32.64 --accelerator-threads $NT  --shm-mpi 1  --device-mem 32768
-mpiexec -launcher ssh -n 2 -host localhost  ./wrap4gpu.sh ./Benchmark_dwf_fp32 --mpi 2.1.1.1 --grid 64.32.32.32 --accelerator-threads $NT  --shm-mpi 1  --device-mem 32768
-done
-#mpiexec -launcher ssh -n 2 -host localhost  ./wrap4gpu.sh ./Benchmark_halo --mpi 1.1.1.2 --grid 32.32.32.64 --accelerator-threads $NT  --shm-mpi 1 > halo.2tile.1x2.log
-#mpiexec -launcher ssh -n 2 -host localhost  ./wrap4gpu.sh ./Benchmark_halo --mpi 2.1.1.1 --grid 64.32.32.32 --accelerator-threads $NT  --shm-mpi 1 > halo.2tile.2x1.log
-
-

systems/PVC/benchmarks/wrap.sh

@@ -1,14 +0,0 @@
-#!/bin/sh
-
-export ZE_AFFINITY_MASK=0.$MPI_LOCALRANKID
-
-echo Ranke $MPI_LOCALRANKID ZE_AFFINITY_MASK is $ZE_AFFINITY_MASK
-
-
-if [ $MPI_LOCALRANKID = "0" ] 
-then
-#  ~psteinbr/build_pti/ze_tracer -h $@
-  onetrace --chrome-device-timeline $@
-else
-  $@
-fi

systems/PVC/config-command

@@ -1,16 +0,0 @@
-INSTALL=/nfs/site/home/azusayax/install
-../../configure \
-	--enable-simd=GPU \
-	--enable-gen-simd-width=64 \
-	--enable-comms=mpi-auto \
-	--disable-accelerator-cshift \
-	--disable-gparity \
-	--disable-fermion-reps \
-	--enable-shm=nvlink \
-	--enable-accelerator=sycl \
-	--enable-unified=no \
-	MPICXX=mpicxx \
-	CXX=dpcpp \
-	LDFLAGS="-fsycl-device-code-split=per_kernel -fsycl-device-lib=all -lze_loader -L$INSTALL/lib" \
-	CXXFLAGS="-fsycl-unnamed-lambda -fsycl -no-fma -I$INSTALL/include -Wno-tautological-compare"
-

systems/PVC/setup.sh

@@ -1,11 +0,0 @@
-export https_proxy=http://proxy-chain.intel.com:911
-export LD_LIBRARY_PATH=/nfs/site/home/azusayax/install/lib:$LD_LIBRARY_PATH
-
-module load intel-release
-source /opt/intel/oneapi/PVC_setup.sh
-#source /opt/intel/oneapi/ATS_setup.sh
-module load intel/mpich/pvc45.3
-export PATH=~/ATS/pti-gpu/tools/onetrace/:$PATH
-
-#clsh embargo-ci-neo-022845
-#source /opt/intel/vtune_amplifier/amplxe-vars.sh

systems/mac-arm/config-command-mpi

@@ -1 +0,0 @@
-CXX=mpicxx-openmpi-mp CXXFLAGS=-I/opt/local/include/ LDFLAGS=-L/opt/local/lib/ ../../configure --enable-simd=GEN --enable-debug --enable-comms=mpi --enable-unified=no

tests/IO/Test_nersc_io.cc

@@ -147,7 +147,7 @@ int main (int argc, char ** argv)
   Complex p  = TensorRemove(Tp);
   std::cout<<GridLogMessage << "calculated plaquettes " <<p*PlaqScale<<std::endl;
 
-  Complex LinkTraceScale(1.0/vol/4.0/(Real)Nc);
+  Complex LinkTraceScale(1.0/vol/4.0/3.0);
   TComplex Tl = sum(LinkTrace);
   Complex l  = TensorRemove(Tl);
   std::cout<<GridLogMessage << "calculated link trace " <<l*LinkTraceScale<<std::endl;
@@ -157,10 +157,8 @@ int main (int argc, char ** argv)
   Complex ll= TensorRemove(TcP);
   std::cout<<GridLogMessage << "coarsened plaquettes sum to " <<ll*PlaqScale<<std::endl;
 
-  const string stNc   = to_string( Nc   ) ;
+  std::string clone2x3("./ckpoint_clone2x3.4000");
-  const string stNcM1 = to_string( Nc-1 ) ;
+  std::string clone3x3("./ckpoint_clone3x3.4000");
-  std::string clone2x3("./ckpoint_clone"+stNcM1+"x"+stNc+".4000");
-  std::string clone3x3("./ckpoint_clone"+stNc+"x"+stNc+".4000");
 
   NerscIO::writeConfiguration(Umu,clone3x3,0,precision32);
   NerscIO::writeConfiguration(Umu,clone2x3,1,precision32);

tests/Test_simd.cc

@@ -793,7 +793,6 @@ int main (int argc, char ** argv)
     }
     std::cout <<" OK ! "<<std::endl;
 
-#ifdef USE_FP16
     // Double to Half
     std::cout << GridLogMessage<< "Double to half" ;
     precisionChange(&H[0],&D[0],Ndp);
@@ -823,7 +822,6 @@ int main (int argc, char ** argv)
       assert( tmp < 1.0e-3 );
     }
     std::cout <<" OK ! "<<std::endl;
-#endif
 
   }
   Grid_finalize();

tests/core/Test_compact_wilson_clover_speedup.cc

@@ -1,226 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./tests/core/Test_compact_wilson_clover_speedup.cc
-
-    Copyright (C) 2020 - 2022
-
-    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
-    Author: Nils Meyer       <nils.meyer@ur.de>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-
-#include <Grid/Grid.h>
-
-using namespace Grid;
-
-NAMESPACE_BEGIN(CommandlineHelpers);
-
-static bool checkPresent(int* argc, char*** argv, const std::string& option) {
-  return GridCmdOptionExists(*argv, *argv + *argc, option);
-}
-
-static std::string getContent(int* argc, char*** argv, const std::string& option) {
-  return GridCmdOptionPayload(*argv, *argv + *argc, option);
-}
-
-static int readInt(int* argc, char*** argv, std::string&& option, int defaultValue) {
-  std::string arg;
-  int         ret = defaultValue;
-  if(checkPresent(argc, argv, option)) {
-    arg = getContent(argc, argv, option);
-    GridCmdOptionInt(arg, ret);
-  }
-  return ret;
-}
-
-static float readFloat(int* argc, char*** argv, std::string&& option, float defaultValue) {
-  std::string arg;
-  float       ret = defaultValue;
-  if(checkPresent(argc, argv, option)) {
-    arg = getContent(argc, argv, option);
-    GridCmdOptionFloat(arg, ret);
-  }
-  return ret;
-}
-
-NAMESPACE_END(CommandlineHelpers);
-
-
-#define _grid_printf(LOGGER, ...) \
-  { \
-    if((LOGGER).isActive()) { /* this makes it safe to put, e.g., norm2 in the calling code w.r.t. performance */ \
-      char _printf_buf[1024]; \
-      std::sprintf(_printf_buf, __VA_ARGS__); \
-      std::cout << (LOGGER) << _printf_buf; \
-      fflush(stdout); \
-    } \
-  }
-#define grid_printf_msg(...) _grid_printf(GridLogMessage, __VA_ARGS__)
-
-
-template<typename Field>
-bool resultsAgree(const Field& ref, const Field& res, const std::string& name) {
-  RealD checkTolerance = (getPrecision<Field>::value == 2) ? 1e-15 : 1e-7;
-  Field diff(ref.Grid());
-  diff = ref - res;
-  auto absDev = norm2(diff);
-  auto relDev = absDev / norm2(ref);
-  std::cout << GridLogMessage
-            << "norm2(reference), norm2(" << name << "), abs. deviation, rel. deviation: " << norm2(ref) << " "
-            << norm2(res) << " " << absDev << " " << relDev << " -> check "
-            << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
-
-  return relDev <= checkTolerance;
-}
-
-
-template<typename vCoeff_t>
-void runBenchmark(int* argc, char*** argv) {
-  // read from command line
-  const int   nIter        = CommandlineHelpers::readInt(     argc, argv, "--niter", 1000);
-  const RealD mass         = CommandlineHelpers::readFloat(   argc, argv, "--mass",  0.5);
-  const RealD csw          = CommandlineHelpers::readFloat(   argc, argv, "--csw",   1.0);
-  const RealD cF           = CommandlineHelpers::readFloat(   argc, argv, "--cF",    1.0);
-  const bool  antiPeriodic = CommandlineHelpers::checkPresent(argc, argv, "--antiperiodic");
-
-  // precision
-  static_assert(getPrecision<vCoeff_t>::value == 2 || getPrecision<vCoeff_t>::value == 1, "Incorrect precision"); // double or single
-  std::string precision = (getPrecision<vCoeff_t>::value == 2 ? "double" : "single");
-
-  // setup grids
-  GridCartesian*         UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vCoeff_t::Nsimd()), GridDefaultMpi());
-  GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-  // clang-format on
-
-  // setup rng
-  std::vector<int> seeds({1, 2, 3, 4});
-  GridParallelRNG  pRNG(UGrid);
-  pRNG.SeedFixedIntegers(seeds);
-
-  // type definitions
-  typedef WilsonImpl<vCoeff_t, FundamentalRepresentation, CoeffReal> WImpl;
-  typedef WilsonCloverFermion<WImpl, CloverHelpers<WImpl>> WilsonCloverOperator;
-  typedef CompactWilsonCloverFermion<WImpl, CompactCloverHelpers<WImpl>> CompactWilsonCloverOperator;
-  typedef typename WilsonCloverOperator::FermionField Fermion;
-  typedef typename WilsonCloverOperator::GaugeField Gauge;
-
-  // setup fields
-  Fermion src(UGrid); random(pRNG, src);
-  Fermion ref(UGrid); ref = Zero();
-  Fermion res(UGrid); res = Zero();
-  Fermion hop(UGrid); hop = Zero();
-  Fermion diff(UGrid); diff = Zero();
-  Gauge   Umu(UGrid); SU3::HotConfiguration(pRNG, Umu);
-
-  // setup boundary phases
-  typename WilsonCloverOperator::ImplParams implParams;
-  std::vector<Complex> boundary_phases(Nd, 1.);
-  if(antiPeriodic) boundary_phases[Nd-1] = -1.;
-  implParams.boundary_phases = boundary_phases;
-  WilsonAnisotropyCoefficients anisParams;
-
-  // misc stuff needed for benchmarks
-  double volume=1.0; for(int mu=0; mu<Nd; mu++) volume*=UGrid->_fdimensions[mu];
-
-  // setup fermion operators
-  WilsonCloverOperator        Dwc(        Umu, *UGrid, *UrbGrid, mass, csw, csw,     anisParams, implParams);
-  CompactWilsonCloverOperator Dwc_compact(Umu, *UGrid, *UrbGrid, mass, csw, csw, cF, anisParams, implParams);
-
-  // now test the conversions
-  typename CompactWilsonCloverOperator::CloverField         tmp_ref(UGrid);  tmp_ref  = Dwc.CloverTerm;
-  typename CompactWilsonCloverOperator::CloverField         tmp_res(UGrid);  tmp_res  = Zero();
-  typename CompactWilsonCloverOperator::CloverField         tmp_diff(UGrid); tmp_diff = Zero();
-  typename CompactWilsonCloverOperator::CloverDiagonalField diagonal(UGrid); diagonal = Zero();
-  typename CompactWilsonCloverOperator::CloverTriangleField triangle(UGrid); diagonal = Zero();
-  CompactWilsonCloverOperator::CompactHelpers::ConvertLayout(tmp_ref, diagonal, triangle);
-  CompactWilsonCloverOperator::CompactHelpers::ConvertLayout(diagonal, triangle, tmp_res);
-  tmp_diff = tmp_ref - tmp_res;
-  std::cout << GridLogMessage << "conversion: ref, res, diff, eps"
-            << " " << norm2(tmp_ref)
-            << " " << norm2(tmp_res)
-            << " " << norm2(tmp_diff)
-            << " " << norm2(tmp_diff) / norm2(tmp_ref)
-            << std::endl;
-
-  // performance per site (use minimal values necessary)
-  double hop_flop_per_site            = 1320; // Rich's Talk + what Peter uses
-  double hop_byte_per_site            = (8 * 9 + 9 * 12) * 2 * getPrecision<vCoeff_t>::value * 4;
-  double clov_flop_per_site           = 504; // Rich's Talk and 1412.2629
-  double clov_byte_per_site           = (2 * 18 + 12 + 12) * 2 * getPrecision<vCoeff_t>::value * 4;
-  double clov_flop_per_site_performed = 1128;
-  double clov_byte_per_site_performed = (12 * 12 + 12 + 12) * 2 * getPrecision<vCoeff_t>::value * 4;
-
-  // total performance numbers
-  double hop_gflop_total            = volume * nIter * hop_flop_per_site / 1e9;
-  double hop_gbyte_total            = volume * nIter * hop_byte_per_site / 1e9;
-  double clov_gflop_total           = volume * nIter * clov_flop_per_site / 1e9;
-  double clov_gbyte_total           = volume * nIter * clov_byte_per_site / 1e9;
-  double clov_gflop_performed_total = volume * nIter * clov_flop_per_site_performed / 1e9;
-  double clov_gbyte_performed_total = volume * nIter * clov_byte_per_site_performed / 1e9;
-
-  // warmup + measure dhop
-  for(auto n : {1, 2, 3, 4, 5}) Dwc.Dhop(src, hop, 0);
-  double t0 = usecond();
-  for(int n = 0; n < nIter; n++) Dwc.Dhop(src, hop, 0);
-  double t1 = usecond();
-  double secs_hop = (t1-t0)/1e6;
-  grid_printf_msg("Performance(%35s, %s): %2.4f s, %6.0f GFlop/s, %6.0f GByte/s, speedup vs ref = %.2f, fraction of hop = %.2f\n",
-              "hop", precision.c_str(), secs_hop, hop_gflop_total/secs_hop, hop_gbyte_total/secs_hop, 0.0, secs_hop/secs_hop);
-
-#define BENCH_CLOVER_KERNEL(KERNEL) { \
-  /* warmup + measure reference clover */ \
-  for(auto n : {1, 2, 3, 4, 5}) Dwc.KERNEL(src, ref); \
-  double t2 = usecond(); \
-  for(int n = 0; n < nIter; n++) Dwc.KERNEL(src, ref); \
-  double t3 = usecond(); \
-  double secs_ref = (t3-t2)/1e6; \
-  grid_printf_msg("Performance(%35s, %s): %2.4f s, %6.0f GFlop/s, %6.0f GByte/s, speedup vs ref = %.2f, fraction of hop = %.2f\n", \
-                  "reference_"#KERNEL, precision.c_str(), secs_ref, clov_gflop_total/secs_ref, clov_gbyte_total/secs_ref, secs_ref/secs_ref, secs_ref/secs_hop); \
-  grid_printf_msg("Performance(%35s, %s): %2.4f s, %6.0f GFlop/s, %6.0f GByte/s, speedup vs ref = %.2f, fraction of hop = %.2f\n", /* to see how well the ET performs */  \
-                  "reference_"#KERNEL"_performed", precision.c_str(), secs_ref, clov_gflop_performed_total/secs_ref, clov_gbyte_performed_total/secs_ref, secs_ref/secs_ref, secs_ref/secs_hop); \
-\
-  /* warmup + measure compact clover */ \
-  for(auto n : {1, 2, 3, 4, 5}) Dwc_compact.KERNEL(src, res); \
-  double t4 = usecond(); \
-  for(int n = 0; n < nIter; n++) Dwc_compact.KERNEL(src, res); \
-  double t5 = usecond(); \
-  double secs_res = (t5-t4)/1e6; \
-  grid_printf_msg("Performance(%35s, %s): %2.4f s, %6.0f GFlop/s, %6.0f GByte/s, speedup vs ref = %.2f, fraction of hop = %.2f\n", \
-                  "compact_"#KERNEL, precision.c_str(), secs_res, clov_gflop_total/secs_res, clov_gbyte_total/secs_res, secs_ref/secs_res, secs_res/secs_hop); \
-  assert(resultsAgree(ref, res, #KERNEL)); \
-}
-
-  BENCH_CLOVER_KERNEL(Mooee);
-  BENCH_CLOVER_KERNEL(MooeeDag);
-  BENCH_CLOVER_KERNEL(MooeeInv);
-  BENCH_CLOVER_KERNEL(MooeeInvDag);
-
-  grid_printf_msg("finalize %s\n", precision.c_str());
-}
-
-int main(int argc, char** argv) {
-  Grid_init(&argc, &argv);
-
-  runBenchmark<vComplexD>(&argc, &argv);
-  runBenchmark<vComplexF>(&argc, &argv);
-
-  Grid_finalize();
-}