Remove Zmobius for faster compile on this bracnh only

.travis.yml

@@ -0,0 +1,56 @@
+language: cpp
+
+cache:
+  directories:
+    - clang
+
+matrix:
+  include:
+    - os:        osx
+      osx_image: xcode8.3
+      compiler: clang
+      
+before_install:
+    - export GRIDDIR=`pwd`
+    - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]] && [ ! -e clang/bin ]; then wget $CLANG_LINK; tar -xf `basename $CLANG_LINK`; mkdir clang; mv clang+*/* clang/; fi
+    - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]]; then export PATH="${GRIDDIR}/clang/bin:${PATH}"; fi
+    - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]]; then export LD_LIBRARY_PATH="${GRIDDIR}/clang/lib:${LD_LIBRARY_PATH}"; fi
+    - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew update; fi
+    - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install libmpc openssl; fi
+    
+install:
+    - export CWD=`pwd`
+    - echo $CWD
+    - export CC=$CC$VERSION
+    - export CXX=$CXX$VERSION
+    - echo $PATH
+    - which autoconf
+    - autoconf  --version
+    - which automake
+    - automake  --version
+    - which $CC
+    - $CC  --version
+    - which $CXX
+    - $CXX --version
+    - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then export LDFLAGS='-L/usr/local/lib'; fi
+    - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then export EXTRACONF='--with-openssl=/usr/local/opt/openssl'; fi
+    
+script:
+    - ./bootstrap.sh
+    - mkdir build
+    - cd build
+    - mkdir lime
+    - cd lime
+    - mkdir build
+    - cd build
+    - wget http://usqcd-software.github.io/downloads/c-lime/lime-1.3.2.tar.gz
+    - tar xf lime-1.3.2.tar.gz
+    - cd lime-1.3.2
+    - ./configure --prefix=$CWD/build/lime/install
+    - make -j4
+    - make install
+    - cd $CWD/build
+    - ../configure --enable-simd=SSE4 --enable-comms=none --with-lime=$CWD/build/lime/install ${EXTRACONF}
+    - make -j4 
+    - ./benchmarks/Benchmark_dwf --threads 1 --debug-signals
+    - make check

Grid/Makefile.am

@@ -54,11 +54,9 @@ Version.h: version-cache
 include Make.inc
 include Eigen.inc
 
+#extra_sources+=$(ZWILS_FERMION_FILES)
 extra_sources+=$(WILS_FERMION_FILES)
 extra_sources+=$(STAG_FERMION_FILES)
-if BUILD_ZMOBIUS
-  extra_sources+=$(ZWILS_FERMION_FILES)
-endif
 if BUILD_GPARITY
   extra_sources+=$(GP_FERMION_FILES)
 endif

Grid/communicator/SharedMemoryMPI.cc

@@ -7,7 +7,6 @@
     Copyright (C) 2015
 
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Christoph Lehner <christoph@lhnr.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
@@ -170,23 +169,6 @@ static inline int divides(int a,int b)
 }
 void GlobalSharedMemory::GetShmDims(const Coordinate &WorldDims,Coordinate &ShmDims)
 {
-  ////////////////////////////////////////////////////////////////
-  // Allow user to configure through environment variable
-  ////////////////////////////////////////////////////////////////
-  char* str = getenv(("GRID_SHM_DIMS_" + std::to_string(ShmDims.size())).c_str());
-  if ( str ) {
-    std::vector<int> IntShmDims;
-    GridCmdOptionIntVector(std::string(str),IntShmDims);
-    assert(IntShmDims.size() == WorldDims.size());
-    long ShmSize = 1;
-    for (int dim=0;dim<WorldDims.size();dim++) {
-      ShmSize *= (ShmDims[dim] = IntShmDims[dim]);
-      assert(divides(ShmDims[dim],WorldDims[dim]));
-    }
-    assert(ShmSize == WorldShmSize);
-    return;
-  }
-  
   ////////////////////////////////////////////////////////////////
   // Powers of 2,3,5 only in prime decomposition for now
   ////////////////////////////////////////////////////////////////

Grid/cshift/Cshift_common.h

@@ -110,11 +110,9 @@ Gather_plane_extract(const Lattice<vobj> &rhs,
   int n1=rhs.Grid()->_slice_stride[dimension];
 
   if ( cbmask ==0x3){
-#ifdef ACCELERATOR_CSHIFT
+#ifdef ACCELERATOR_CSHIFT    
     autoView(rhs_v , rhs, AcceleratorRead);
-    accelerator_for(nn,e1*e2,1,{
+    accelerator_for2d(n,e1,b,e2,1,{
-	int n = nn%e1;
-	int b = nn/e1;
 	int o      =   n*n1;
 	int offset = b+n*e2;
 	
@@ -137,9 +135,7 @@ Gather_plane_extract(const Lattice<vobj> &rhs,
     std::cout << " Dense packed buffer WARNING " <<std::endl; // Does this get called twice once for each cb?
 #ifdef ACCELERATOR_CSHIFT    
     autoView(rhs_v , rhs, AcceleratorRead);
-    accelerator_for(nn,e1*e2,1,{
+    accelerator_for2d(n,e1,b,e2,1,{
-	int n = nn%e1;
-	int b = nn/e1;
 
 	Coordinate coor;
 
@@ -261,9 +257,7 @@ template<class vobj> void Scatter_plane_merge(Lattice<vobj> &rhs,ExtractPointerA
     int _slice_block = rhs.Grid()->_slice_block[dimension];
 #ifdef ACCELERATOR_CSHIFT    
     autoView( rhs_v , rhs, AcceleratorWrite);
-    accelerator_for(nn,e1*e2,1,{
+    accelerator_for2d(n,e1,b,e2,1,{
-	int n = nn%e1;
-	int b = nn/e1;
 	int o      = n*_slice_stride;
 	int offset = b+n*_slice_block;
 	merge(rhs_v[so+o+b],pointers,offset);
@@ -280,7 +274,7 @@ template<class vobj> void Scatter_plane_merge(Lattice<vobj> &rhs,ExtractPointerA
 
     // Case of SIMD split AND checker dim cannot currently be hit, except in 
     // Test_cshift_red_black code.
-    std::cout << "Scatter_plane merge assert(0); think this is buggy FIXME "<< std::endl;// think this is buggy FIXME
+    //    std::cout << "Scatter_plane merge assert(0); think this is buggy FIXME "<< std::endl;// think this is buggy FIXME
     std::cout<<" Unthreaded warning -- buffer is not densely packed ??"<<std::endl;
     assert(0); // This will fail if hit on GPU
     autoView( rhs_v, rhs, CpuWrite);

Grid/cshift/Cshift_mpi.h

@@ -122,8 +122,8 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
   assert(shift<fd);
   
   int buffer_size = rhs.Grid()->_slice_nblock[dimension]*rhs.Grid()->_slice_block[dimension];
-  static cshiftVector<vobj> send_buf; send_buf.resize(buffer_size);
+  cshiftVector<vobj> send_buf(buffer_size);
-  static cshiftVector<vobj> recv_buf; recv_buf.resize(buffer_size);
+  cshiftVector<vobj> recv_buf(buffer_size);
     
   int cb= (cbmask==0x2)? Odd : Even;
   int sshift= rhs.Grid()->CheckerBoardShiftForCB(rhs.Checkerboard(),dimension,shift,cb);
@@ -198,8 +198,8 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
   int buffer_size = grid->_slice_nblock[dimension]*grid->_slice_block[dimension];
   //  int words = sizeof(vobj)/sizeof(vector_type);
 
-  static std::vector<cshiftVector<scalar_object> >  send_buf_extract; send_buf_extract.resize(Nsimd);
+  std::vector<cshiftVector<scalar_object> >  send_buf_extract(Nsimd);
-  static std::vector<cshiftVector<scalar_object> >  recv_buf_extract; recv_buf_extract.resize(Nsimd);
+  std::vector<cshiftVector<scalar_object> >  recv_buf_extract(Nsimd);
   scalar_object *  recv_buf_extract_mpi;
   scalar_object *  send_buf_extract_mpi;
  
@@ -294,8 +294,8 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
   assert(shift<fd);
   
   int buffer_size = rhs.Grid()->_slice_nblock[dimension]*rhs.Grid()->_slice_block[dimension];
-  static cshiftVector<vobj> send_buf_v; send_buf_v.resize(buffer_size);
+  cshiftVector<vobj> send_buf_v(buffer_size);
-  static cshiftVector<vobj> recv_buf_v; recv_buf_v.resize(buffer_size);
+  cshiftVector<vobj> recv_buf_v(buffer_size);
   vobj *send_buf;
   vobj *recv_buf;
   {
@@ -381,8 +381,8 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
   int buffer_size = grid->_slice_nblock[dimension]*grid->_slice_block[dimension];
   //  int words = sizeof(vobj)/sizeof(vector_type);
 
-  static std::vector<cshiftVector<scalar_object> >  send_buf_extract; send_buf_extract.resize(Nsimd);
+  std::vector<cshiftVector<scalar_object> >  send_buf_extract(Nsimd);
-  static std::vector<cshiftVector<scalar_object> >  recv_buf_extract; recv_buf_extract.resize(Nsimd);
+  std::vector<cshiftVector<scalar_object> >  recv_buf_extract(Nsimd);
   scalar_object *  recv_buf_extract_mpi;
   scalar_object *  send_buf_extract_mpi;
   {

Grid/lattice/Lattice_transfer.h

@@ -97,20 +97,6 @@ accelerator_inline void convertType(ComplexF & out, const std::complex<float> &
   out = in;
 }
 
-template<typename T>
-accelerator_inline EnableIf<isGridFundamental<T>> convertType(T & out, const T & in) {
-  out = in;
-}
-
-// This would allow for conversions between GridFundamental types, but is not strictly needed as yet
-/*template<typename T1, typename T2>
-accelerator_inline typename std::enable_if<isGridFundamental<T1>::value && isGridFundamental<T2>::value>::type
-// Or to make this very broad, conversions between anything that's not a GridTensor could be allowed
-//accelerator_inline typename std::enable_if<!isGridTensor<T1>::value && !isGridTensor<T2>::value>::type
-convertType(T1 & out, const T2 & in) {
-  out = in;
-}*/
-
 #ifdef GRID_SIMT
 accelerator_inline void convertType(vComplexF & out, const ComplexF & in) {
   ((ComplexF*)&out)[acceleratorSIMTlane(vComplexF::Nsimd())] = in;
@@ -131,18 +117,23 @@ accelerator_inline void convertType(vComplexD2 & out, const vComplexF & in) {
   Optimization::PrecisionChange::StoD(in.v,out._internal[0].v,out._internal[1].v);
 }
 
-template<typename T1,typename T2>
+template<typename T1,typename T2,int N>
-accelerator_inline void convertType(iScalar<T1> & out, const iScalar<T2> & in) {
+  accelerator_inline void convertType(iMatrix<T1,N> & out, const iMatrix<T2,N> & in);
-  convertType(out._internal,in._internal);
+template<typename T1,typename T2,int N>
+  accelerator_inline void convertType(iVector<T1,N> & out, const iVector<T2,N> & in);
+
+template<typename T1,typename T2, typename std::enable_if<!isGridScalar<T1>::value, T1>::type* = nullptr>
+accelerator_inline void convertType(T1 & out, const iScalar<T2> & in) {
+  convertType(out,in._internal);
 }
 
-template<typename T1,typename T2>
+template<typename T1, typename std::enable_if<!isGridScalar<T1>::value, T1>::type* = nullptr>
-accelerator_inline NotEnableIf<isGridScalar<T1>> convertType(T1 & out, const iScalar<T2> & in) {
+accelerator_inline void convertType(T1 & out, const iScalar<T1> & in) {
   convertType(out,in._internal);
 }
 
 template<typename T1,typename T2>
-accelerator_inline NotEnableIf<isGridScalar<T2>> convertType(iScalar<T1> & out, const T2 & in) {
+accelerator_inline void convertType(iScalar<T1> & out, const T2 & in) {
   convertType(out._internal,in);
 }
 
@@ -159,6 +150,11 @@ accelerator_inline void convertType(iVector<T1,N> & out, const iVector<T2,N> & i
     convertType(out._internal[i],in._internal[i]);
 }
 
+template<typename T, typename std::enable_if<isGridFundamental<T>::value, T>::type* = nullptr>
+accelerator_inline void convertType(T & out, const T & in) {
+  out = in;
+}
+
 template<typename T1,typename T2>
 accelerator_inline void convertType(Lattice<T1> & out, const Lattice<T2> & in) {
   autoView( out_v , out,AcceleratorWrite);

Grid/lattice/Lattice_view.h

@@ -71,10 +71,10 @@ public:
   //  accelerator_inline const vobj & operator[](size_t i) const { return this->_odata[i]; };
   accelerator_inline vobj       & operator[](size_t i) const { return this->_odata[i]; };
 #else
-  accelerator_inline const vobj & operator[](size_t i) const { return this->_odata[i]; };
+  //  accelerator_inline const vobj & operator[](size_t i) const { return this->_odata[i]; };
-  accelerator_inline vobj       & operator[](size_t i)       { return this->_odata[i]; };
+  //  accelerator_inline vobj       & operator[](size_t i)       { return this->_odata[i]; };
 #endif
-  
+
   accelerator_inline uint64_t begin(void) const { return 0;};
   accelerator_inline uint64_t end(void)   const { return this->_odata_size; };
   accelerator_inline uint64_t size(void)  const { return this->_odata_size; };

Grid/lattice/Lattice_where.h

@@ -43,7 +43,7 @@ inline void whereWolf(Lattice<vobj> &ret,const Lattice<iobj> &predicate,Lattice<
   conformable(iftrue,predicate);
   conformable(iftrue,ret);
 
-  GridBase *grid=iftrue.Grid();
+  GridBase *grid=iftrue._grid;
 
   typedef typename vobj::scalar_object scalar_object;
   typedef typename vobj::scalar_type scalar_type;
@@ -52,23 +52,22 @@ inline void whereWolf(Lattice<vobj> &ret,const Lattice<iobj> &predicate,Lattice<
 
   const int Nsimd = grid->Nsimd();
 
-  autoView(iftrue_v,iftrue,CpuRead);
+  std::vector<Integer> mask(Nsimd);
-  autoView(iffalse_v,iffalse,CpuRead);
+  std::vector<scalar_object> truevals (Nsimd);
-  autoView(predicate_v,predicate,CpuRead);
+  std::vector<scalar_object> falsevals(Nsimd);
-  autoView(ret_v,ret,CpuWrite);
+
-  Integer NN= grid->oSites();
+  parallel_for(int ss=0;ss<iftrue._grid->oSites(); ss++){
-  thread_for(ss,NN,{
+
-    Integer mask;
+    extract(iftrue._odata[ss]   ,truevals);
-    scalar_object trueval;
+    extract(iffalse._odata[ss]  ,falsevals);
-    scalar_object falseval;
+    extract<vInteger,Integer>(TensorRemove(predicate._odata[ss]),mask);
-    for(int l=0;l<Nsimd;l++){
+
-      trueval =extractLane(l,iftrue_v[ss]);
+    for(int s=0;s<Nsimd;s++){
-      falseval=extractLane(l,iffalse_v[ss]);
+      if (mask[s]) falsevals[s]=truevals[s];
-      mask    =extractLane(l,predicate_v[ss]);
-      if (mask) falseval=trueval;
-      insertLane(l,ret_v[ss],falseval);
     }
-  });
+
+    merge(ret._odata[ss],falsevals);
+  }
 }
 
 template<class vobj,class iobj>
@@ -77,9 +76,9 @@ inline Lattice<vobj> whereWolf(const Lattice<iobj> &predicate,Lattice<vobj> &ift
   conformable(iftrue,iffalse);
   conformable(iftrue,predicate);
 
-  Lattice<vobj> ret(iftrue.Grid());
+  Lattice<vobj> ret(iftrue._grid);
 
-  whereWolf(ret,predicate,iftrue,iffalse);
+  where(ret,predicate,iftrue,iffalse);
 
   return ret;
 }

Grid/parallelIO/BinaryIO.h

@@ -271,7 +271,7 @@ class BinaryIO {
 			      uint32_t &scidac_csumb)
   {
     grid->Barrier();
-    GridStopWatch timer, insideTimer; 
+    GridStopWatch timer; 
     GridStopWatch bstimer;
     
     nersc_csum=0;
@@ -363,10 +363,7 @@ class BinaryIO {
 	std::cout<< GridLogMessage<<"IOobject: MPI read I/O "<< file<< std::endl;
 	ierr=MPI_File_open(grid->communicator,(char *) file.c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);    assert(ierr==0);
 	ierr=MPI_File_set_view(fh, disp, mpiObject, fileArray, "native", MPI_INFO_NULL);    assert(ierr==0);
-  insideTimer.Start();
+	ierr=MPI_File_read_all(fh, &iodata[0], 1, localArray, &status);    assert(ierr==0);
-	ierr=MPI_File_read_all(fh, &iodata[0], 1, localArray, &status);
-  insideTimer.Stop();
-  assert(ierr==0);
 	MPI_File_close(&fh);
 	MPI_Type_free(&fileArray);
 	MPI_Type_free(&localArray);
@@ -441,9 +438,7 @@ class BinaryIO {
         assert(ierr == 0);
 
         std::cout << GridLogDebug << "MPI write I/O write all " << file << std::endl;
-        insideTimer.Start();
         ierr = MPI_File_write_all(fh, &iodata[0], 1, localArray, &status);
-        insideTimer.Stop();
         assert(ierr == 0);
 
         MPI_Offset os;
@@ -521,13 +516,8 @@ class BinaryIO {
     if ( control & BINARYIO_READ) std::cout << " read  ";
     else                          std::cout << " write ";
     uint64_t bytes = sizeof(fobj)*iodata.size()*nrank;
-    std::cout<< lastPerf.size <<"bytes in "<< timer.Elapsed() <<" "
+    std::cout<< lastPerf.size <<" bytes in "<< timer.Elapsed() <<" "
 	     << lastPerf.mbytesPerSecond <<" MB/s "<<std::endl;
-    std::cout << GridLogMessage << "IOobject: pure MPI IO call " 
-              << lastPerf.size <<" bytes in " 
-              << insideTimer.Elapsed() << " "
-              << lastPerf.size/1024./1024./(insideTimer.useconds()/1.0e6)
-              <<" MB/s "<<std::endl;
 
     std::cout<<GridLogMessage<<"IOobject: endian and checksum overhead "<<bstimer.Elapsed()  <<std::endl;
 

Grid/parallelIO/MetaData.h

@@ -128,7 +128,7 @@ inline void MachineCharacteristics(FieldMetaData &header)
   std::time_t t = std::time(nullptr);
   std::tm tm_ = *std::localtime(&t);
   std::ostringstream oss; 
-  oss << std::put_time(&tm_, "%c %Z");
+  //      oss << std::put_time(&tm_, "%c %Z");
   header.creation_date = oss.str();
   header.archive_date  = header.creation_date;
 

Grid/parallelIO/NerscIO.h

@@ -205,20 +205,11 @@ public:
     std::cout<<GridLogMessage <<"NERSC Configuration "<<file<< " and plaquette, link trace, and checksum agree"<<std::endl;
   }
 
-  // Preferred interface
-  template<class GaugeStats=PeriodicGaugeStatistics>
-  static inline void writeConfiguration(Lattice<vLorentzColourMatrixD > &Umu,
-					std::string file, 
-					std::string ens_label = std::string("DWF"))
-  {
-    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,
+					int bits32)
-					std::string ens_label = std::string("DWF"))
   {
     typedef vLorentzColourMatrixD vobj;
     typedef typename vobj::scalar_object sobj;
@@ -228,8 +219,8 @@ public:
     // Following should become arguments
     ///////////////////////////////////////////
     header.sequence_number = 1;
-    header.ensemble_id     = std::string("UKQCD");
+    header.ensemble_id     = "UKQCD";
-    header.ensemble_label  = ens_label;
+    header.ensemble_label  = "DWF";
 
     typedef LorentzColourMatrixD fobj3D;
     typedef LorentzColour2x3D    fobj2D;
@@ -241,7 +232,7 @@ public:
     GaugeStats Stats; Stats(Umu,header);
     MachineCharacteristics(header);
 
-    uint64_t offset;
+	uint64_t offset;
 
     // Sod it -- always write 3x3 double
     header.floating_point = std::string("IEEE64BIG");

Grid/qcd/action/ActionBase.h

@@ -41,7 +41,7 @@ class Action
 public:
   bool is_smeared = false;
   // Heatbath?
-  virtual void refresh(const GaugeField& U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) = 0; // refresh pseudofermions
+  virtual void refresh(const GaugeField& U, GridParallelRNG& pRNG) = 0; // refresh pseudofermions
   virtual RealD S(const GaugeField& U) = 0;                             // evaluate the action
   virtual void deriv(const GaugeField& U, GaugeField& dSdU) = 0;        // evaluate the action derivative
   virtual std::string action_name()    = 0;                             // return the action name

Grid/qcd/action/fermion/Fermion.h

@@ -291,6 +291,12 @@ typedef ImprovedStaggeredFermion5D<StaggeredImplR> ImprovedStaggeredFermion5DR;
 typedef ImprovedStaggeredFermion5D<StaggeredImplF> ImprovedStaggeredFermion5DF;
 typedef ImprovedStaggeredFermion5D<StaggeredImplD> ImprovedStaggeredFermion5DD;
 
+#ifndef GRID_CUDA
+typedef ImprovedStaggeredFermion5D<StaggeredVec5dImplR> ImprovedStaggeredFermionVec5dR;
+typedef ImprovedStaggeredFermion5D<StaggeredVec5dImplF> ImprovedStaggeredFermionVec5dF;
+typedef ImprovedStaggeredFermion5D<StaggeredVec5dImplD> ImprovedStaggeredFermionVec5dD;
+#endif
+
 NAMESPACE_END(Grid);
 
 ////////////////////

Grid/qcd/action/fermion/FermionOperatorImpl.h

@@ -153,8 +153,8 @@ public:
   typedef typename Impl::StencilImpl             StencilImpl;		\
   typedef typename Impl::ImplParams               ImplParams;	        \
   typedef typename Impl::StencilImpl::View_type  StencilView;		\
-  typedef const typename ViewMap<FermionField>::Type      FermionFieldView;	\
+  typedef typename ViewMap<FermionField>::Type      FermionFieldView;	\
-  typedef const typename ViewMap<DoubledGaugeField>::Type DoubledGaugeFieldView;
+  typedef typename ViewMap<DoubledGaugeField>::Type DoubledGaugeFieldView;
 
 #define INHERIT_IMPL_TYPES(Base)		\
   INHERIT_GIMPL_TYPES(Base)			\
@@ -183,8 +183,7 @@ NAMESPACE_CHECK(ImplStaggered);
 /////////////////////////////////////////////////////////////////////////////
 // Single flavour one component spinors with colour index. 5d vec
 /////////////////////////////////////////////////////////////////////////////
-// Deprecate Vec5d
+#include <Grid/qcd/action/fermion/StaggeredVec5dImpl.h> 
-//#include <Grid/qcd/action/fermion/StaggeredVec5dImpl.h> 
+NAMESPACE_CHECK(ImplStaggered5dVec);  
-//NAMESPACE_CHECK(ImplStaggered5dVec);  
 
 

Grid/qcd/action/fermion/GparityWilsonImpl.h

@@ -88,7 +88,7 @@ public:
 					  const _Spinor &chi, 
 					  int mu, 
 					  StencilEntry *SE,
-					  StencilView &St) 
+					  const StencilView &St) 
   {
     int direction = St._directions[mu];
     int distance  = St._distances[mu];

Grid/qcd/action/fermion/StaggeredImpl.h

@@ -72,23 +72,19 @@ public:
     
   StaggeredImpl(const ImplParams &p = ImplParams()) : Params(p){};
       
-  template<class _Spinor>
+  static accelerator_inline void multLink(SiteSpinor &phi,
-  static accelerator_inline void multLink(_Spinor &phi,
 		       const SiteDoubledGaugeField &U,
-		       const _Spinor &chi,
+		       const SiteSpinor &chi,
 		       int mu)
   {
-    auto UU = coalescedRead(U(mu));
+    mult(&phi(), &U(mu), &chi());
-    mult(&phi(), &UU, &chi());
   }
-  template<class _Spinor>
+  static accelerator_inline void multLinkAdd(SiteSpinor &phi,
-  static accelerator_inline void multLinkAdd(_Spinor &phi,
 			  const SiteDoubledGaugeField &U,
-			  const _Spinor &chi,
+			  const SiteSpinor &chi,
 			  int mu)
   {
-    auto UU = coalescedRead(U(mu));
+    mac(&phi(), &U(mu), &chi());
-    mac(&phi(), &UU, &chi());
   }
       
   template <class ref>

Grid/qcd/action/fermion/StaggeredKernels.h

@@ -56,8 +56,12 @@ template<class Impl> class StaggeredKernels : public FermionOperator<Impl> , pub
 		 DoubledGaugeField &U,
 		 const FermionField &in, FermionField &out, int dag, int interior,int exterior);
   
-  void DhopDirKernel(StencilImpl &st, DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU, SiteSpinor * buf,
+  void DhopDirKernel(StencilImpl &st,
-		     int sF, int sU, const FermionFieldView &in, FermionFieldView &out, int dir,int disp);
+		     const DoubledGaugeFieldView &U,
+		     const DoubledGaugeFieldView &UUU, SiteSpinor * buf,
+		     int sF, int sU,
+		     const FermionFieldView &in,
+		     const FermionFieldView &out, int dir,int disp);
  protected:    
 
    ///////////////////////////////////////////////////////////////////////////////////////
@@ -65,53 +69,67 @@ template<class Impl> class StaggeredKernels : public FermionOperator<Impl> , pub
    ///////////////////////////////////////////////////////////////////////////////////////
    template<int Naik> 
    static accelerator_inline
-   void DhopSiteGeneric(StencilView &st, 
+   void DhopSiteGeneric(const StencilView &st, 
-			DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU, 
+			const DoubledGaugeFieldView &U,
+			const DoubledGaugeFieldView &UUU, 
 			SiteSpinor * buf, int LLs, int sU, 
-			const FermionFieldView &in, FermionFieldView &out,int dag);
+			const FermionFieldView &in,
+			const FermionFieldView &out,int dag);
    
    template<int Naik> static accelerator_inline
-   void DhopSiteGenericInt(StencilView &st, 
+   void DhopSiteGenericInt(const StencilView &st, 
-			   DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU, 
+			   const DoubledGaugeFieldView &U,
+			   const DoubledGaugeFieldView &UUU, 
 			   SiteSpinor * buf, int LLs, int sU, 
-			   const FermionFieldView &in, FermionFieldView &out,int dag);
+			   const FermionFieldView &in,
+			   const FermionFieldView &out,int dag);
    
    template<int Naik> static accelerator_inline
-   void DhopSiteGenericExt(StencilView &st, 
+   void DhopSiteGenericExt(const StencilView &st, 
-			   DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU,
+			   const DoubledGaugeFieldView &U,
-			   SiteSpinor * buf, int LLs, int sU, 
+			   const DoubledGaugeFieldView &UUU,
-			   const FermionFieldView &in, FermionFieldView &out,int dag);
+ 			   SiteSpinor * buf, int LLs, int sU, 
+			   const FermionFieldView &in,
+			   const FermionFieldView &out,int dag);
 
    ///////////////////////////////////////////////////////////////////////////////////////
    // Nc=3 specific kernels
    ///////////////////////////////////////////////////////////////////////////////////////
    
    template<int Naik> static accelerator_inline
-   void DhopSiteHand(StencilView &st, 
+   void DhopSiteHand(const StencilView &st, 
-		     DoubledGaugeFieldView &U,DoubledGaugeFieldView &UUU, 
+		     const DoubledGaugeFieldView &U,
-		     SiteSpinor * buf, int LLs, int sU, 
+		     const DoubledGaugeFieldView &UUU, 
-		     const FermionFieldView &in, FermionFieldView &out,int dag);
+ 		     SiteSpinor * buf, int LLs, int sU, 
+		     const FermionFieldView &in,
+		     const FermionFieldView &out,int dag);
    
    template<int Naik> static accelerator_inline
-   void DhopSiteHandInt(StencilView &st, 
+   void DhopSiteHandInt(const StencilView &st, 
-			DoubledGaugeFieldView &U,DoubledGaugeFieldView &UUU, 
+			const DoubledGaugeFieldView &U,
-			SiteSpinor * buf, int LLs, int sU, 
+			const DoubledGaugeFieldView &UUU, 
-			const FermionFieldView &in, FermionFieldView &out,int dag);
+ 			SiteSpinor * buf, int LLs, int sU, 
+			const FermionFieldView &in,
+			const FermionFieldView &out,int dag);
    
    template<int Naik> static accelerator_inline
-   void DhopSiteHandExt(StencilView &st, 
+   void DhopSiteHandExt(const StencilView &st, 
-			DoubledGaugeFieldView &U,DoubledGaugeFieldView &UUU, 
+			const DoubledGaugeFieldView &U,
-			SiteSpinor * buf, int LLs, int sU, 
+			const DoubledGaugeFieldView &UUU, 
-			const FermionFieldView &in, FermionFieldView &out,int dag);
+ 			SiteSpinor * buf, int LLs, int sU, 
+			const FermionFieldView &in,
+			const FermionFieldView &out,int dag);
 
    ///////////////////////////////////////////////////////////////////////////////////////
    // Asm Nc=3 specific kernels
    ///////////////////////////////////////////////////////////////////////////////////////
    
-   void DhopSiteAsm(StencilView &st, 
+   void DhopSiteAsm(const StencilView &st, 
-		    DoubledGaugeFieldView &U,DoubledGaugeFieldView &UUU, 
+		    const DoubledGaugeFieldView &U,
-		    SiteSpinor * buf, int LLs, int sU, 
+		    const DoubledGaugeFieldView &UUU, 
-		    const FermionFieldView &in, FermionFieldView &out,int dag);
+ 		    SiteSpinor * buf, int LLs, int sU, 
+		    const FermionFieldView &in,
+		    const FermionFieldView &out,int dag);
   
 public:
 

Grid/qcd/action/fermion/WilsonCloverFermion.h

@@ -245,7 +245,7 @@ public:
     return out;
   }
 
-protected:
+private:
   // here fixing the 4 dimensions, make it more general?
 
   RealD csw_r;                                               // Clover coefficient - spatial

Grid/qcd/action/fermion/WilsonImpl.h

@@ -72,7 +72,7 @@ public:
   typedef WilsonCompressor<SiteHalfCommSpinor,SiteHalfSpinor, SiteSpinor> Compressor;
   typedef WilsonImplParams ImplParams;
   typedef WilsonStencil<SiteSpinor, SiteHalfSpinor,ImplParams> StencilImpl;
-  typedef const typename StencilImpl::View_type StencilView;
+  typedef typename StencilImpl::View_type StencilView;
     
   ImplParams Params;
 
@@ -95,7 +95,7 @@ public:
 					  const _Spinor &chi,
 					  int mu,
 					  StencilEntry *SE,
-					  StencilView &St) 
+					  const StencilView &St) 
   {
     multLink(phi,U,chi,mu);
   }
@@ -184,22 +184,18 @@ public:
       mat = TraceIndex<SpinIndex>(P); 
     }
       
-    inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds)
+    inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
-    {
       for (int mu = 0; mu < Nd; mu++)
       mat[mu] = PeekIndex<LorentzIndex>(Uds, mu);
     }
 
-  inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu)
+
-  {
+  inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){
-#undef USE_OLD_INSERT_FORCE    
+      
     int Ls=Btilde.Grid()->_fdimensions[0];
-    autoView( mat_v , mat, AcceleratorWrite);
-#ifdef USE_OLD_INSERT_FORCE    
     GaugeLinkField tmp(mat.Grid());
     tmp = Zero();
     {
-      const int Nsimd = SiteSpinor::Nsimd();
       autoView( tmp_v , tmp, AcceleratorWrite);
       autoView( Btilde_v , Btilde, AcceleratorRead);
       autoView( Atilde_v , Atilde, AcceleratorRead);
@@ -212,29 +208,6 @@ public:
 	});
     }
     PokeIndex<LorentzIndex>(mat,tmp,mu);
-#else
-    {
-      const int Nsimd = SiteSpinor::Nsimd();
-      autoView( Btilde_v , Btilde, AcceleratorRead);
-      autoView( Atilde_v , Atilde, AcceleratorRead);
-      accelerator_for(sss,mat.Grid()->oSites(),Nsimd,{
-	  int sU=sss;
-  	  typedef decltype(coalescedRead(mat_v[sU](mu)() )) ColorMatrixType;
-  	  ColorMatrixType sum;
-	  zeroit(sum);  
-	  for(int s=0;s<Ls;s++){
-	    int sF = s+Ls*sU;
-  	    for(int spn=0;spn<Ns;spn++){ //sum over spin
-  	      auto bb = coalescedRead(Btilde_v[sF]()(spn) ); //color vector
-  	      auto aa = coalescedRead(Atilde_v[sF]()(spn) );
-	      auto op = outerProduct(bb,aa);
-  	      sum = sum + op;
-	    }
-	  }
-  	  coalescedWrite(mat_v[sU](mu)(), sum);
-      });
-    }
-#endif    
   }
 };
 

Grid/qcd/action/fermion/WilsonKernels.h

@@ -49,17 +49,10 @@ public:
 
   INHERIT_IMPL_TYPES(Impl);
   typedef FermionOperator<Impl> Base;
-  typedef AcceleratorVector<int,STENCIL_MAX> StencilVector;   
+  typedef AcceleratorVector<int,STENCIL_MAX> StencilVector;
+   
 public:
 
-#ifdef GRID_SYCL
-#define SYCL_HACK
-#endif  
-#ifdef SYCL_HACK
-  static void HandDhopSiteSycl(StencilVector st_perm,StencilEntry *st_p, SiteDoubledGaugeField *U,SiteHalfSpinor  *buf,
-			       int ss,int sU,const SiteSpinor *in, SiteSpinor *out);
-#endif
-  
   static void DhopKernel(int Opt,StencilImpl &st,  DoubledGaugeField &U, SiteHalfSpinor * buf,
 			 int Ls, int Nsite, const FermionField &in, FermionField &out,
 			 int interior=1,int exterior=1) ;
@@ -76,73 +69,87 @@ public:
 
 private:
 
-  static accelerator_inline void DhopDirK(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor * buf,
+  static accelerator_inline void DhopDirK(const StencilView &st, const DoubledGaugeFieldView &U,
-				   int sF, int sU, const FermionFieldView &in, FermionFieldView &out, int dirdisp, int gamma);
+					  SiteHalfSpinor * buf, int sF, int sU,
+					  const FermionFieldView &in,const FermionFieldView &out, int dirdisp, int gamma);
 
-  static accelerator_inline void DhopDirXp(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,const FermionFieldView &in,FermionFieldView &out,int dirdisp);
+  static accelerator_inline void DhopDirXp(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,
-  static accelerator_inline void DhopDirYp(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,const FermionFieldView &in,FermionFieldView &out,int dirdisp);
+					   const FermionFieldView &in, const FermionFieldView &out,int dirdisp);
-  static accelerator_inline void DhopDirZp(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,const FermionFieldView &in,FermionFieldView &out,int dirdisp);
+  static accelerator_inline void DhopDirYp(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,
-  static accelerator_inline void DhopDirTp(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,const FermionFieldView &in,FermionFieldView &out,int dirdisp);
+					   const FermionFieldView &in, const FermionFieldView &out,int dirdisp);
-  static accelerator_inline void DhopDirXm(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,const FermionFieldView &in,FermionFieldView &out,int dirdisp);
+  static accelerator_inline void DhopDirZp(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,
-  static accelerator_inline void DhopDirYm(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,const FermionFieldView &in,FermionFieldView &out,int dirdisp);
+					   const FermionFieldView &in, const FermionFieldView &out,int dirdisp);
-  static accelerator_inline void DhopDirZm(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,const FermionFieldView &in,FermionFieldView &out,int dirdisp);
+  static accelerator_inline void DhopDirTp(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,
-  static accelerator_inline void DhopDirTm(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,const FermionFieldView &in,FermionFieldView &out,int dirdisp);
+					   const FermionFieldView &in, const FermionFieldView &out,int dirdisp);
+  static accelerator_inline void DhopDirXm(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,
+					   const FermionFieldView &in, const FermionFieldView &out,int dirdisp);
+  static accelerator_inline void DhopDirYm(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,
+					   const FermionFieldView &in, const FermionFieldView &out,int dirdisp);
+  static accelerator_inline void DhopDirZm(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,
+					   const FermionFieldView &in, const FermionFieldView &out,int dirdisp);
+  static accelerator_inline void DhopDirTm(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,int sF,int sU,
+					   const FermionFieldView &in, const FermionFieldView &out,int dirdisp);
       
   // Specialised variants
-  static accelerator void GenericDhopSite(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+  static accelerator void GenericDhopSite(const StencilView &st,
-					  int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+					  const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-      
+					  int sF, int sU, const FermionFieldView &in, const FermionFieldView &out);
-  static accelerator void GenericDhopSiteDag(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+       
-						    int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+  static accelerator void GenericDhopSiteDag(const StencilView &st, const  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-  
+					     int sF, int sU, const FermionFieldView &in, const FermionFieldView &out);
-  static accelerator void GenericDhopSiteInt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+   
-						    int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+  static accelerator void GenericDhopSiteInt(const StencilView &st, const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-      
+					     int sF, int sU, const FermionFieldView &in, const FermionFieldView &out);
-  static accelerator void GenericDhopSiteDagInt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+       
-						int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+  static accelerator void GenericDhopSiteDagInt(const StencilView &st, const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-  
+						int sF, int sU, const FermionFieldView &in, const FermionFieldView &out);
-  static accelerator void GenericDhopSiteExt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+   
-					     int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+  static accelerator void GenericDhopSiteExt(const StencilView &st, const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-      
+					     int sF, int sU, const FermionFieldView &in, const FermionFieldView &out);
-  static accelerator void GenericDhopSiteDagExt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+       
-						       int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+  static accelerator void GenericDhopSiteDagExt(const StencilView &st, const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+						int sF, int sU, const FermionFieldView &in, const FermionFieldView &out);
 
-  static void AsmDhopSite(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+// Keep Hand unrolled 
-			  int sF, int sU, int Ls, int Nsite, const FermionFieldView &in,FermionFieldView &out);
+  static accelerator void HandDhopSiteSycl(StencilVector st_perm, StencilEntry *st_p,  SiteDoubledGaugeField *U, SiteHalfSpinor * buf,
-  
+					   int sF, int sU, const SiteSpinor *in, SiteSpinor *out);
-  static void AsmDhopSiteDag(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-			     int sF, int sU, int Ls, int Nsite, const FermionFieldView &in, FermionFieldView &out);
-  
-  static void AsmDhopSiteInt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-			     int sF, int sU, int Ls, int Nsite, const FermionFieldView &in,FermionFieldView &out);
-  
-  static void AsmDhopSiteDagInt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-				int sF, int sU, int Ls, int Nsite, const FermionFieldView &in, FermionFieldView &out);
-  
-  static void AsmDhopSiteExt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-			     int sF, int sU, int Ls, int Nsite, const FermionFieldView &in,FermionFieldView &out);
-  
-  static void AsmDhopSiteDagExt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-				int sF, int sU, int Ls, int Nsite, const FermionFieldView &in, FermionFieldView &out);
 
-// Keep Hand unrolled temporarily  
+  static accelerator void HandDhopSite(const StencilView &st, const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-  static accelerator void HandDhopSite(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+				       int sF, int sU, const FermionFieldView &in,const FermionFieldView &out);
-				       int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+   
+  static accelerator void HandDhopSiteDag(const StencilView &st, const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+					  int sF, int sU, const FermionFieldView &in, const FermionFieldView &out);
+   
+  static accelerator void HandDhopSiteInt(const StencilView &st, const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+					  int sF, int sU, const FermionFieldView &in, const FermionFieldView &out);
   
-  static accelerator void HandDhopSiteDag(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+  static accelerator void HandDhopSiteDagInt(const StencilView &st, const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-					  int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+					     int sF, int sU, const FermionFieldView &in, const FermionFieldView &out);
   
-  static accelerator void HandDhopSiteInt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+  static accelerator void HandDhopSiteExt(const StencilView &st, const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-					  int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+					  int sF, int sU, const FermionFieldView &in, const FermionFieldView &out);
+   
+  static accelerator void HandDhopSiteDagExt(const StencilView &st, const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+					     int sF, int sU, const FermionFieldView &in, const FermionFieldView &out);
+  //AVX 512 ASM
+  static void AsmDhopSite(const StencilView &st,  const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+			  int sF, int sU, int Ls, int Nsite, const FermionFieldView &in,const FermionFieldView &out);
   
-  static accelerator void HandDhopSiteDagInt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+  static void AsmDhopSiteDag(const StencilView &st,  const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-					     int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+			     int sF, int sU, int Ls, int Nsite, const FermionFieldView &in, const FermionFieldView &out);
   
-  static accelerator void HandDhopSiteExt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+  static void AsmDhopSiteInt(const StencilView &st,  const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-					  int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+			     int sF, int sU, int Ls, int Nsite, const FermionFieldView &in,const FermionFieldView &out);
   
-  static accelerator void HandDhopSiteDagExt(StencilView &st,  DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+  static void AsmDhopSiteDagInt(const StencilView &st,  const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
-					     int sF, int sU, const FermionFieldView &in, FermionFieldView &out);
+				int sF, int sU, int Ls, int Nsite, const FermionFieldView &in, const FermionFieldView &out);
+  
+  static void AsmDhopSiteExt(const StencilView &st,  const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+			     int sF, int sU, int Ls, int Nsite, const FermionFieldView &in,const FermionFieldView &out);
+  
+  static void AsmDhopSiteDagExt(const StencilView &st,  const DoubledGaugeFieldView &U, SiteHalfSpinor * buf,
+				int sF, int sU, int Ls, int Nsite, const FermionFieldView &in, const FermionFieldView &out);
+
  public:
  WilsonKernels(const ImplParams &p = ImplParams()) : Base(p){};
 };

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

@@ -880,23 +880,11 @@ void CayleyFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
   }
 
   std::vector<RealD> G_s(Ls,1.0);
-  Integer sign = 1; // sign flip for vector/tadpole
   if ( curr_type == Current::Axial ) {
     for(int s=0;s<Ls/2;s++){
       G_s[s] = -1.0;
     }
   }
-  else if ( curr_type == Current::Tadpole ) {
-    auto b=this->_b;
-    auto c=this->_c;
-    if ( b == 1 && c == 0 ) {
-      sign = -1;    
-    }
-    else {
-      std::cerr << "Error: Tadpole implementation currently unavailable for non-Shamir actions." << std::endl;
-      assert(b==1 && c==0);
-    }
-  }
 
   for(int s=0;s<Ls;s++){
 
@@ -919,7 +907,7 @@ void CayleyFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
 
     tmp    = Cshift(tmp,mu,1);
     Impl::multLinkField(Utmp,this->Umu,tmp,mu);
-    tmp    = sign*G_s[s]*( Utmp*ph - gmu*Utmp*ph ); // Forward hop
+    tmp    = G_s[s]*( Utmp*ph - gmu*Utmp*ph ); // Forward hop
     tmp    = where((lcoor>=tmin),tmp,zz); // Mask the time 
     L_Q    = where((lcoor<=tmax),tmp,zz); // Position of current complicated
 

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

@@ -618,11 +618,13 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 NAMESPACE_BEGIN(Grid);
 
 template <class Impl>
-void StaggeredKernels<Impl>::DhopSiteAsm(StencilView &st,
+void StaggeredKernels<Impl>::DhopSiteAsm(const StencilView &st,
-					 DoubledGaugeFieldView &U,
+					 const DoubledGaugeFieldView &U,
-					 DoubledGaugeFieldView &UUU,
+					 const DoubledGaugeFieldView &UUU,
-					 SiteSpinor *buf, int sF,
+ 					 SiteSpinor *buf, int sF,
-					 int sU, const FermionFieldView &in, FermionFieldView &out,int dag) 
+					 int sU,
+					 const FermionFieldView &in,
+					 const FermionFieldView &out,int dag) 
 {
   assert(0);
 };
@@ -680,15 +682,16 @@ void StaggeredKernels<Impl>::DhopSiteAsm(StencilView &st,
   gauge2 =(uint64_t)&UU[sU]( Z );				\
   gauge3 =(uint64_t)&UU[sU]( T ); 
   
-#undef STAG_VEC5D
+
-#ifdef STAG_VEC5D
   // This is the single precision 5th direction vectorised kernel
 #include <Grid/simd/Intel512single.h>
-template <> void StaggeredKernels<StaggeredVec5dImplF>::DhopSiteAsm(StencilView &st,
+template <> void StaggeredKernels<StaggeredVec5dImplF>::DhopSiteAsm(const StencilView &st,
-								    DoubledGaugeFieldView &U,
+								    const DoubledGaugeFieldView &U,
-								    DoubledGaugeFieldView &UUU,
+								    const DoubledGaugeFieldView &UUU,
-								    SiteSpinor *buf, int sF,
+ 								    SiteSpinor *buf, int sF,
-								    int sU, const FermionFieldView &in, FermionFieldView &out,int dag) 
+								    int sU,
+								    const FermionFieldView &in,
+								    const FermionFieldView &out,int dag) 
 {
 #ifdef AVX512
   uint64_t gauge0,gauge1,gauge2,gauge3;
@@ -739,11 +742,13 @@ template <> void StaggeredKernels<StaggeredVec5dImplF>::DhopSiteAsm(StencilView
 }
 
 #include <Grid/simd/Intel512double.h>
-template <> void StaggeredKernels<StaggeredVec5dImplD>::DhopSiteAsm(StencilView &st, 
+template <> void StaggeredKernels<StaggeredVec5dImplD>::DhopSiteAsm(const StencilView &st, 
-								    DoubledGaugeFieldView &U,
+								    const DoubledGaugeFieldView &U,
-								    DoubledGaugeFieldView &UUU,
+								    const DoubledGaugeFieldView &UUU,
-								    SiteSpinor *buf, int sF,
+ 								    SiteSpinor *buf, int sF,
-								    int sU, const FermionFieldView &in, FermionFieldView &out, int dag) 
+								    int sU,
+								    const FermionFieldView &in,
+								    const FermionFieldView &out, int dag) 
 {
 #ifdef AVX512
   uint64_t gauge0,gauge1,gauge2,gauge3;
@@ -791,7 +796,7 @@ template <> void StaggeredKernels<StaggeredVec5dImplD>::DhopSiteAsm(StencilView
 #endif
 }
    
-#endif   
+   
 
 
 #define PERMUTE_DIR3 __asm__ (	\
@@ -825,11 +830,13 @@ template <> void StaggeredKernels<StaggeredVec5dImplD>::DhopSiteAsm(StencilView
   // This is the single precision 5th direction vectorised kernel
 
 #include <Grid/simd/Intel512single.h>
-template <> void StaggeredKernels<StaggeredImplF>::DhopSiteAsm(StencilView &st, 
+template <> void StaggeredKernels<StaggeredImplF>::DhopSiteAsm(const StencilView &st, 
-							       DoubledGaugeFieldView &U,
+							       const DoubledGaugeFieldView &U,
-							       DoubledGaugeFieldView &UUU,
+							       const DoubledGaugeFieldView &UUU,
-							       SiteSpinor *buf, int sF,
+ 							       SiteSpinor *buf, int sF,
-							       int sU, const FermionFieldView &in, FermionFieldView &out,int dag) 
+							       int sU,
+							       const FermionFieldView &in,
+							       const FermionFieldView &out,int dag) 
 {
 #ifdef AVX512
   uint64_t gauge0,gauge1,gauge2,gauge3;
@@ -894,11 +901,13 @@ template <> void StaggeredKernels<StaggeredImplF>::DhopSiteAsm(StencilView &st,
 }
 
 #include <Grid/simd/Intel512double.h>
-template <> void StaggeredKernels<StaggeredImplD>::DhopSiteAsm(StencilView &st, 
+template <> void StaggeredKernels<StaggeredImplD>::DhopSiteAsm(const StencilView &st, 
-							       DoubledGaugeFieldView &U,
+							       const DoubledGaugeFieldView &U,
-							       DoubledGaugeFieldView &UUU,
+							       const DoubledGaugeFieldView &UUU,
-							       SiteSpinor *buf, int sF,
+ 							       SiteSpinor *buf, int sF,
-							       int sU, const FermionFieldView &in, FermionFieldView &out,int dag) 
+							       int sU,
+							       const FermionFieldView &in,
+							       const FermionFieldView &out,int dag) 
 {
 #ifdef AVX512
   uint64_t gauge0,gauge1,gauge2,gauge3;

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

@@ -32,50 +32,25 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 NAMESPACE_BEGIN(Grid);
 
-#ifdef GRID_SIMT
+#define LOAD_CHI(b)		\
-
-#define LOAD_CHI(ptype,b)			\
-  const SiteSpinor & ref (b[offset]);				\
-  Chi_0=coalescedReadPermute<ptype>(ref()()(0),perm,lane);	\
-  Chi_1=coalescedReadPermute<ptype>(ref()()(1),perm,lane);	\
-  Chi_2=coalescedReadPermute<ptype>(ref()()(2),perm,lane);
-
-#define LOAD_CHI_COMMS(b)		\
   const SiteSpinor & ref (b[offset]);	\
-  Chi_0=coalescedRead(ref()()(0),lane);	\
+    Chi_0=ref()()(0);\
-  Chi_1=coalescedRead(ref()()(1),lane);	\
+    Chi_1=ref()()(1);\
-  Chi_2=coalescedRead(ref()()(2),lane);
+    Chi_2=ref()()(2);
-
-#define PERMUTE_DIR(dir)	;
-#else
-#define LOAD_CHI(ptype,b)      LOAD_CHI_COMMS(b)
-
-#define LOAD_CHI_COMMS(b)		\
-  const SiteSpinor & ref (b[offset]);	\
-  Chi_0=ref()()(0);			\
-  Chi_1=ref()()(1);			\
-  Chi_2=ref()()(2);
-
-#define PERMUTE_DIR(dir)			\
-  permute##dir(Chi_0,Chi_0);			\
-  permute##dir(Chi_1,Chi_1);			\
-  permute##dir(Chi_2,Chi_2);
-
-#endif
 
 
 // To splat or not to splat depends on the implementation
 #define MULT(A,UChi)				\
   auto & ref(U[sU](A));			\
-    U_00=coalescedRead(ref()(0,0),lane);				\
+   Impl::loadLinkElement(U_00,ref()(0,0));      \
-    U_10=coalescedRead(ref()(1,0),lane);				\
+   Impl::loadLinkElement(U_10,ref()(1,0));      \
-    U_20=coalescedRead(ref()(2,0),lane);				\
+   Impl::loadLinkElement(U_20,ref()(2,0));      \
-    U_01=coalescedRead(ref()(0,1),lane);				\
+   Impl::loadLinkElement(U_01,ref()(0,1));      \
-    U_11=coalescedRead(ref()(1,1),lane);				\
+   Impl::loadLinkElement(U_11,ref()(1,1));      \
-    U_21=coalescedRead(ref()(2,1),lane);				\
+   Impl::loadLinkElement(U_21,ref()(2,1));      \
-    U_02=coalescedRead(ref()(0,2),lane);				\
+   Impl::loadLinkElement(U_02,ref()(0,2));     \
-    U_12=coalescedRead(ref()(1,2),lane);				\
+   Impl::loadLinkElement(U_12,ref()(1,2));     \
-    U_22=coalescedRead(ref()(2,2),lane);				\
+   Impl::loadLinkElement(U_22,ref()(2,2));     \
     UChi ## _0  = U_00*Chi_0;	       \
     UChi ## _1  = U_10*Chi_0;\
     UChi ## _2  = U_20*Chi_0;\
@@ -88,15 +63,15 @@ NAMESPACE_BEGIN(Grid);
 
 #define MULT_ADD(U,A,UChi)			\
   auto & ref(U[sU](A));			\
-    U_00=coalescedRead(ref()(0,0),lane);				\
+   Impl::loadLinkElement(U_00,ref()(0,0));      \
-    U_10=coalescedRead(ref()(1,0),lane);				\
+   Impl::loadLinkElement(U_10,ref()(1,0));      \
-    U_20=coalescedRead(ref()(2,0),lane);				\
+   Impl::loadLinkElement(U_20,ref()(2,0));      \
-    U_01=coalescedRead(ref()(0,1),lane);				\
+   Impl::loadLinkElement(U_01,ref()(0,1));      \
-    U_11=coalescedRead(ref()(1,1),lane);				\
+   Impl::loadLinkElement(U_11,ref()(1,1));      \
-    U_21=coalescedRead(ref()(2,1),lane);				\
+   Impl::loadLinkElement(U_21,ref()(2,1));      \
-    U_02=coalescedRead(ref()(0,2),lane);				\
+   Impl::loadLinkElement(U_02,ref()(0,2));     \
-    U_12=coalescedRead(ref()(1,2),lane);				\
+   Impl::loadLinkElement(U_12,ref()(1,2));     \
-    U_22=coalescedRead(ref()(2,2),lane);				\
+   Impl::loadLinkElement(U_22,ref()(2,2));     \
     UChi ## _0 += U_00*Chi_0;	       \
     UChi ## _1 += U_10*Chi_0;\
     UChi ## _2 += U_20*Chi_0;\
@@ -108,18 +83,24 @@ NAMESPACE_BEGIN(Grid);
     UChi ## _2 += U_22*Chi_2;
 
 
+#define PERMUTE_DIR(dir)			\
+  permute##dir(Chi_0,Chi_0);			\
+  permute##dir(Chi_1,Chi_1);			\
+  permute##dir(Chi_2,Chi_2);
+
+
 #define HAND_STENCIL_LEG_BASE(Dir,Perm,skew)	\
   SE=st.GetEntry(ptype,Dir+skew,sF);	\
   offset = SE->_offset;			\
   local  = SE->_is_local;		\
   perm   = SE->_permute;		\
   if ( local ) {						\
-    LOAD_CHI(Perm,in);						\
+    LOAD_CHI(in);					\
     if ( perm) {						\
       PERMUTE_DIR(Perm);					\
     }								\
   } else {							\
-    LOAD_CHI_COMMS(buf);					\
+    LOAD_CHI(buf);						\
   }								
 
 #define HAND_STENCIL_LEG_BEGIN(Dir,Perm,skew,even)		\
@@ -135,18 +116,19 @@ NAMESPACE_BEGIN(Grid);
   }
 
 
+
 #define HAND_STENCIL_LEG_INT(U,Dir,Perm,skew,even)	\
   SE=st.GetEntry(ptype,Dir+skew,sF);			\
   offset = SE->_offset;					\
   local  = SE->_is_local;				\
   perm   = SE->_permute;				\
   if ( local ) {					\
-    LOAD_CHI(Perm,in);					\
+    LOAD_CHI(in);				\
     if ( perm) {					\
       PERMUTE_DIR(Perm);				\
     }							\
   } else if ( st.same_node[Dir] ) {			\
-    LOAD_CHI_COMMS(buf);				\
+    LOAD_CHI(buf);					\
   }							\
   if (local || st.same_node[Dir] ) {		\
     MULT_ADD(U,Dir,even);				\
@@ -158,51 +140,45 @@ NAMESPACE_BEGIN(Grid);
   local  = SE->_is_local;				\
   if ((!local) && (!st.same_node[Dir]) ) {		\
     nmu++;							\
-    { LOAD_CHI_COMMS(buf);	  }				\
+    { LOAD_CHI(buf);	  }					\
     { MULT_ADD(U,Dir,even); }					\
   }								
 
-#define HAND_DECLARATIONS(Simd) \
-  Simd even_0;			\
-  Simd even_1;			\
-  Simd even_2;			\
-  Simd odd_0;			\
-  Simd odd_1;			\
-  Simd odd_2;		        \
-		      		\
-  Simd Chi_0;			\
-  Simd Chi_1;			\
-  Simd Chi_2;			\
-				\
-  Simd U_00;			\
-  Simd U_10;			\
-  Simd U_20;			\
-  Simd U_01;			\
-  Simd U_11;			\
-  Simd U_21;			\
-  Simd U_02;			\
-  Simd U_12;			\
-  Simd U_22;			
-  
 
 template <class Impl>
 template <int Naik> accelerator_inline
-void StaggeredKernels<Impl>::DhopSiteHand(StencilView &st,
+void StaggeredKernels<Impl>::DhopSiteHand(const StencilView &st,
-					  DoubledGaugeFieldView &U,DoubledGaugeFieldView &UUU,
+					  const DoubledGaugeFieldView &U,
-					  SiteSpinor *buf, int sF, int sU, 
+					  const DoubledGaugeFieldView &UUU,
-					  const FermionFieldView &in, FermionFieldView &out,int dag) 
+ 					  SiteSpinor *buf, int sF, int sU, 
+					  const FermionFieldView &in,
+					  const FermionFieldView &out,int dag) 
 {
   typedef typename Simd::scalar_type S;
   typedef typename Simd::vector_type V;
 
+  Simd even_0; // 12 regs on knc
+  Simd even_1;
+  Simd even_2;
+  Simd odd_0; // 12 regs on knc
+  Simd odd_1;
+  Simd odd_2;
 
-  const int Nsimd = SiteHalfSpinor::Nsimd();
+  Simd Chi_0;    // two spinor; 6 regs
-  const int lane=acceleratorSIMTlane(Nsimd);
+  Simd Chi_1;
-  typedef decltype( coalescedRead( in[0]()()(0) )) Simt;
+  Simd Chi_2;
-  HAND_DECLARATIONS(Simt);
+  
+  Simd U_00;  // two rows of U matrix
+  Simd U_10;
+  Simd U_20;  
+  Simd U_01;
+  Simd U_11;
+  Simd U_21;  // 2 reg left.
+  Simd U_02;
+  Simd U_12;
+  Simd U_22; 
 
-  typedef decltype( coalescedRead( in[0] )) calcSiteSpinor;
+  SiteSpinor result;
-  calcSiteSpinor result;
   int offset,local,perm, ptype;
 
   StencilEntry *SE;
@@ -241,28 +217,45 @@ void StaggeredKernels<Impl>::DhopSiteHand(StencilView &st,
       result()()(1) = even_1 + odd_1;
       result()()(2) = even_2 + odd_2;
     }
-    coalescedWrite(out[sF],result);
+    vstream(out[sF],result);
   }
 }
 
 
 template <class Impl>
 template <int Naik> accelerator_inline
-void StaggeredKernels<Impl>::DhopSiteHandInt(StencilView &st, 
+void StaggeredKernels<Impl>::DhopSiteHandInt(const StencilView &st, 
-					     DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU,
+					     const DoubledGaugeFieldView &U,
-					     SiteSpinor *buf, int sF, int sU, 
+					     const DoubledGaugeFieldView &UUU,
-					     const FermionFieldView &in, FermionFieldView &out,int dag) 
+ 					     SiteSpinor *buf, int sF, int sU, 
+					     const FermionFieldView &in,
+					     const FermionFieldView &out,int dag) 
 {
   typedef typename Simd::scalar_type S;
   typedef typename Simd::vector_type V;
 
-  const int Nsimd = SiteHalfSpinor::Nsimd();
+  Simd even_0; // 12 regs on knc
-  const int lane=acceleratorSIMTlane(Nsimd);
+  Simd even_1;
-  typedef decltype( coalescedRead( in[0]()()(0) )) Simt;
+  Simd even_2;
-  HAND_DECLARATIONS(Simt);
+  Simd odd_0; // 12 regs on knc
+  Simd odd_1;
+  Simd odd_2;
 
-  typedef decltype( coalescedRead( in[0] )) calcSiteSpinor;
+  Simd Chi_0;    // two spinor; 6 regs
-  calcSiteSpinor result;
+  Simd Chi_1;
+  Simd Chi_2;
+  
+  Simd U_00;  // two rows of U matrix
+  Simd U_10;
+  Simd U_20;  
+  Simd U_01;
+  Simd U_11;
+  Simd U_21;  // 2 reg left.
+  Simd U_02;
+  Simd U_12;
+  Simd U_22; 
+
+  SiteSpinor result;
   int offset, ptype, local, perm;
 
   StencilEntry *SE;
@@ -272,8 +265,8 @@ void StaggeredKernels<Impl>::DhopSiteHandInt(StencilView &st,
   //    int sF=s+LLs*sU;
   {
 
-    zeroit(even_0);    zeroit(even_1);    zeroit(even_2);
+    even_0 = Zero();    even_1 = Zero();    even_2 = Zero();
-    zeroit(odd_0);    zeroit(odd_1);    zeroit(odd_2);
+     odd_0 = Zero();     odd_1 = Zero();     odd_2 = Zero();
 
     skew = 0;
     HAND_STENCIL_LEG_INT(U,Xp,3,skew,even);  
@@ -305,28 +298,45 @@ void StaggeredKernels<Impl>::DhopSiteHandInt(StencilView &st,
       result()()(1) = even_1 + odd_1;
       result()()(2) = even_2 + odd_2;
     }
-    coalescedWrite(out[sF],result);
+    vstream(out[sF],result);
   }
 }
 
 
 template <class Impl>
 template <int Naik> accelerator_inline
-void StaggeredKernels<Impl>::DhopSiteHandExt(StencilView &st,
+void StaggeredKernels<Impl>::DhopSiteHandExt(const StencilView &st,
-					     DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU,
+					     const DoubledGaugeFieldView &U,
-					     SiteSpinor *buf, int sF, int sU, 
+					     const DoubledGaugeFieldView &UUU,
-					     const FermionFieldView &in, FermionFieldView &out,int dag) 
+ 					     SiteSpinor *buf, int sF, int sU, 
+					     const FermionFieldView &in,
+					     const FermionFieldView &out,int dag) 
 {
   typedef typename Simd::scalar_type S;
   typedef typename Simd::vector_type V;
 
-  const int Nsimd = SiteHalfSpinor::Nsimd();
+  Simd even_0; // 12 regs on knc
-  const int lane=acceleratorSIMTlane(Nsimd);
+  Simd even_1;
-  typedef decltype( coalescedRead( in[0]()()(0) )) Simt;
+  Simd even_2;
-  HAND_DECLARATIONS(Simt);
+  Simd odd_0; // 12 regs on knc
+  Simd odd_1;
+  Simd odd_2;
 
-  typedef decltype( coalescedRead( in[0] )) calcSiteSpinor;
+  Simd Chi_0;    // two spinor; 6 regs
-  calcSiteSpinor result;
+  Simd Chi_1;
+  Simd Chi_2;
+  
+  Simd U_00;  // two rows of U matrix
+  Simd U_10;
+  Simd U_20;  
+  Simd U_01;
+  Simd U_11;
+  Simd U_21;  // 2 reg left.
+  Simd U_02;
+  Simd U_12;
+  Simd U_22; 
+
+  SiteSpinor result;
   int offset, ptype, local;
 
   StencilEntry *SE;
@@ -336,8 +346,8 @@ void StaggeredKernels<Impl>::DhopSiteHandExt(StencilView &st,
   //    int sF=s+LLs*sU;
   {
 
-    zeroit(even_0);    zeroit(even_1);    zeroit(even_2);
+    even_0 = Zero();    even_1 = Zero();    even_2 = Zero();
-    zeroit(odd_0);    zeroit(odd_1);    zeroit(odd_2);
+     odd_0 = Zero();     odd_1 = Zero();     odd_2 = Zero();
     int nmu=0;
     skew = 0;
     HAND_STENCIL_LEG_EXT(U,Xp,3,skew,even);  
@@ -370,7 +380,7 @@ void StaggeredKernels<Impl>::DhopSiteHandExt(StencilView &st,
 	result()()(1) = even_1 + odd_1;
 	result()()(2) = even_2 + odd_2;
       }
-      coalescedWrite(out[sF] , out(sF)+ result);
+      out[sF] = out[sF] + result;
     }
   }
 }
@@ -393,7 +403,6 @@ void StaggeredKernels<Impl>::DhopSiteHandExt(StencilView &st,
 						     const FermionFieldView &in, FermionFieldView &out, int dag); \
 */
 #undef LOAD_CHI
-#undef HAND_DECLARATIONS
 
 NAMESPACE_END(Grid);
 

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

@@ -35,32 +35,39 @@ NAMESPACE_BEGIN(Grid);
 #define GENERIC_STENCIL_LEG(U,Dir,skew,multLink)		\
   SE = st.GetEntry(ptype, Dir+skew, sF);			\
   if (SE->_is_local ) {						\
-    int perm= SE->_permute;						\
+    if (SE->_permute) {						\
-    chi = coalescedReadPermute(in[SE->_offset],ptype,perm,lane);\
+      chi_p = χ						\
+      permute(chi,  in[SE->_offset], ptype);			\
+    } else {							\
+      chi_p = &in[SE->_offset];					\
+    }								\
   } else {							\
-    chi = coalescedRead(buf[SE->_offset],lane);			\
+    chi_p = &buf[SE->_offset];					\
   }								\
-  acceleratorSynchronise();					\
+  multLink(Uchi, U[sU], *chi_p, Dir);			
-  multLink(Uchi, U[sU], chi, Dir);			
 
 #define GENERIC_STENCIL_LEG_INT(U,Dir,skew,multLink)		\
   SE = st.GetEntry(ptype, Dir+skew, sF);			\
   if (SE->_is_local ) {						\
-    int perm= SE->_permute;						\
+    if (SE->_permute) {						\
-    chi = coalescedReadPermute(in[SE->_offset],ptype,perm,lane);\
+      chi_p = χ						\
+      permute(chi,  in[SE->_offset], ptype);			\
+    } else {							\
+      chi_p = &in[SE->_offset];					\
+    }								\
   } else if ( st.same_node[Dir] ) {				\
-    chi = coalescedRead(buf[SE->_offset],lane);                 \
+    chi_p = &buf[SE->_offset];					\
   }								\
   if (SE->_is_local || st.same_node[Dir] ) {			\
-    multLink(Uchi, U[sU], chi, Dir);				\
+    multLink(Uchi, U[sU], *chi_p, Dir);				\
   }
 
 #define GENERIC_STENCIL_LEG_EXT(U,Dir,skew,multLink)		\
   SE = st.GetEntry(ptype, Dir+skew, sF);			\
   if ((!SE->_is_local) && (!st.same_node[Dir]) ) {		\
     nmu++;							\
-    chi = coalescedRead(buf[SE->_offset],lane);			\
+    chi_p = &buf[SE->_offset];					\
-    multLink(Uchi, U[sU], chi, Dir);				\
+    multLink(Uchi, U[sU], *chi_p, Dir);				\
   }
 
 template <class Impl>
@@ -72,19 +79,17 @@ StaggeredKernels<Impl>::StaggeredKernels(const ImplParams &p) : Base(p){};
 ////////////////////////////////////////////////////////////////////////////////////
 template <class Impl>
 template <int Naik> accelerator_inline
-void StaggeredKernels<Impl>::DhopSiteGeneric(StencilView &st, 
+void StaggeredKernels<Impl>::DhopSiteGeneric(const StencilView &st, 
-					     DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU,
+					     const DoubledGaugeFieldView &U, const DoubledGaugeFieldView &UUU,
-					     SiteSpinor *buf, int sF, int sU, 
+ 					     SiteSpinor *buf, int sF, int sU, 
-					     const FermionFieldView &in, FermionFieldView &out, int dag) 
+					     const FermionFieldView &in, const FermionFieldView &out, int dag) 
 {
-  typedef decltype(coalescedRead(in[0])) calcSpinor;
+  const SiteSpinor *chi_p;
-  calcSpinor chi;
+  SiteSpinor chi;
-  calcSpinor Uchi;
+  SiteSpinor Uchi;
   StencilEntry *SE;
   int ptype;
   int skew;
-  const int Nsimd = SiteHalfSpinor::Nsimd();
-  const int lane=acceleratorSIMTlane(Nsimd);
 
   //  for(int s=0;s<LLs;s++){
   //
@@ -113,7 +118,7 @@ void StaggeredKernels<Impl>::DhopSiteGeneric(StencilView &st,
     if ( dag ) { 
       Uchi = - Uchi;
     } 
-    coalescedWrite(out[sF], Uchi,lane);
+    vstream(out[sF], Uchi);
   }
 };
 
@@ -122,19 +127,17 @@ void StaggeredKernels<Impl>::DhopSiteGeneric(StencilView &st,
   ///////////////////////////////////////////////////
 template <class Impl>
 template <int Naik> accelerator_inline
-void StaggeredKernels<Impl>::DhopSiteGenericInt(StencilView &st, 
+void StaggeredKernels<Impl>::DhopSiteGenericInt(const StencilView &st, 
-						DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU,
+						const DoubledGaugeFieldView &U, const DoubledGaugeFieldView &UUU,
-						SiteSpinor *buf, int sF, int sU, 
+ 						SiteSpinor *buf, int sF, int sU, 
-						const FermionFieldView &in, FermionFieldView &out,int dag)
+						const FermionFieldView &in, const FermionFieldView &out,int dag)
 {
-  typedef decltype(coalescedRead(in[0])) calcSpinor;
+  const SiteSpinor *chi_p;
-  calcSpinor chi;
+  SiteSpinor chi;
-  calcSpinor Uchi;
+  SiteSpinor Uchi;
   StencilEntry *SE;
   int ptype;
   int skew ;
-  const int Nsimd = SiteHalfSpinor::Nsimd();
-  const int lane=acceleratorSIMTlane(Nsimd);
 
   //  for(int s=0;s<LLs;s++){
   //    int sF=LLs*sU+s;
@@ -163,7 +166,7 @@ void StaggeredKernels<Impl>::DhopSiteGenericInt(StencilView &st,
     if ( dag ) {
       Uchi = - Uchi;
     }
-    coalescedWrite(out[sF], Uchi,lane);
+    vstream(out[sF], Uchi);
   }
 };
 
@@ -173,20 +176,20 @@ void StaggeredKernels<Impl>::DhopSiteGenericInt(StencilView &st,
   ///////////////////////////////////////////////////
 template <class Impl>
 template <int Naik> accelerator_inline
-void StaggeredKernels<Impl>::DhopSiteGenericExt(StencilView &st, 
+void StaggeredKernels<Impl>::DhopSiteGenericExt(const StencilView &st, 
-						DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU,
+						const DoubledGaugeFieldView &U,
-						SiteSpinor *buf, int sF, int sU,
+						const DoubledGaugeFieldView &UUU,
-						const FermionFieldView &in, FermionFieldView &out,int dag)
+ 						SiteSpinor *buf, int sF, int sU,
+						const FermionFieldView &in,
+						const FermionFieldView &out,int dag)
 {
-  typedef decltype(coalescedRead(in[0])) calcSpinor;
+  const SiteSpinor *chi_p;
-  calcSpinor chi;
+  //  SiteSpinor chi;
-  calcSpinor Uchi;
+  SiteSpinor Uchi;
   StencilEntry *SE;
   int ptype;
   int nmu=0;
   int skew ;
-  const int Nsimd = SiteHalfSpinor::Nsimd();
-  const int lane=acceleratorSIMTlane(Nsimd);
 
   //  for(int s=0;s<LLs;s++){
   //    int sF=LLs*sU+s;
@@ -212,12 +215,11 @@ void StaggeredKernels<Impl>::DhopSiteGenericExt(StencilView &st,
     GENERIC_STENCIL_LEG_EXT(UUU,Zm,skew,Impl::multLinkAdd);
     GENERIC_STENCIL_LEG_EXT(UUU,Tm,skew,Impl::multLinkAdd);
     }
-    if ( nmu ) {
+    if ( nmu ) { 
-      auto _out = coalescedRead(out[sF],lane);
+      if ( dag ) { 
-      if ( dag ) {
+	out[sF] = out[sF] - Uchi;
-	coalescedWrite(out[sF], _out-Uchi,lane);
       } else { 
-	coalescedWrite(out[sF], _out+Uchi,lane);
+	out[sF] = out[sF] + Uchi;
       }
     }
   }
@@ -227,8 +229,13 @@ void StaggeredKernels<Impl>::DhopSiteGenericExt(StencilView &st,
 // Driving / wrapping routine to select right kernel
 ////////////////////////////////////////////////////////////////////////////////////
 template <class Impl> 
-void StaggeredKernels<Impl>::DhopDirKernel(StencilImpl &st, DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU, SiteSpinor * buf,
+void StaggeredKernels<Impl>::DhopDirKernel(StencilImpl &st,
-					   int sF, int sU, const FermionFieldView &in, FermionFieldView &out, int dir,int disp)
+					   const DoubledGaugeFieldView &U,
+					   const DoubledGaugeFieldView &UUU,
+					   SiteSpinor * buf,
+					   int sF, int sU,
+					   const FermionFieldView &in,
+					   const FermionFieldView &out, int dir,int disp)
 {
   // Disp should be either +1,-1,+3,-3
   // What about "dag" ?
@@ -256,15 +263,14 @@ void StaggeredKernels<Impl>::DhopDirKernel(StencilImpl &st, DoubledGaugeFieldVie
   });
 
 template <class Impl> 
-void StaggeredKernels<Impl>::DhopImproved(StencilImpl &st, LebesgueOrder &lo, 
+void StaggeredKernels<Impl>::DhopImproved(StencilImpl &st,
+					  LebesgueOrder &lo, 
 					  DoubledGaugeField &U, DoubledGaugeField &UUU, 
 					  const FermionField &in, FermionField &out, int dag, int interior,int exterior)
 {
   GridBase *FGrid=in.Grid();  
   GridBase *UGrid=U.Grid();  
   typedef StaggeredKernels<Impl> ThisKernel;
-  const int Nsimd = SiteHalfSpinor::Nsimd();
-  const int lane=acceleratorSIMTlane(Nsimd);
   autoView( UUU_v , UUU, AcceleratorRead);
   autoView( U_v   ,   U, AcceleratorRead);
   autoView( in_v  ,  in, AcceleratorRead);
@@ -305,8 +311,6 @@ void StaggeredKernels<Impl>::DhopNaive(StencilImpl &st, LebesgueOrder &lo,
   GridBase *FGrid=in.Grid();  
   GridBase *UGrid=U.Grid();  
   typedef StaggeredKernels<Impl> ThisKernel;
-  const int Nsimd = SiteHalfSpinor::Nsimd();
-  const int lane=acceleratorSIMTlane(Nsimd);
   autoView( UUU_v ,   U, AcceleratorRead);
   autoView( U_v   ,   U, AcceleratorRead);
   autoView( in_v  ,  in, AcceleratorRead);

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

@@ -397,7 +397,6 @@ void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, co
 template <class Impl>
 void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out, int dag)
 {
-  DhopCalls+=2;
   conformable(in.Grid(), _grid);  // verifies full grid
   conformable(in.Grid(), out.Grid());
 
@@ -409,7 +408,6 @@ void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out, int da
 template <class Impl>
 void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out, int dag)
 {
-  DhopCalls++;
   conformable(in.Grid(), _cbgrid);    // verifies half grid
   conformable(in.Grid(), out.Grid());  // drops the cb check
 
@@ -422,7 +420,6 @@ void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out, int
 template <class Impl>
 void WilsonFermion<Impl>::DhopEO(const FermionField &in, FermionField &out,int dag)
 {
-  DhopCalls++;
   conformable(in.Grid(), _cbgrid);    // verifies half grid
   conformable(in.Grid(), out.Grid());  // drops the cb check
 

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

@@ -38,46 +38,46 @@ NAMESPACE_BEGIN(Grid);
 
 ///////////////////////////////////////////////////////////
 // Default to no assembler implementation
-// Will specialise to 
+// Will specialise to AVX512 if available
 ///////////////////////////////////////////////////////////
 template<class Impl> void 
-WilsonKernels<Impl >::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+WilsonKernels<Impl >::AsmDhopSite(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-				  int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+				  int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, const FermionFieldView &out)
 {
   assert(0);
 }
 
 template<class Impl> void 
-WilsonKernels<Impl >::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+WilsonKernels<Impl >::AsmDhopSiteDag(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-				     int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+				     int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, const FermionFieldView &out)
 {
   assert(0);
 }
 
 template<class Impl> void 
-WilsonKernels<Impl >::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+WilsonKernels<Impl >::AsmDhopSiteInt(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-				     int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+				     int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, const FermionFieldView &out)
 {
   assert(0);
 }
 
 template<class Impl> void 
-WilsonKernels<Impl >::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+WilsonKernels<Impl >::AsmDhopSiteDagInt(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-					int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+					int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, const FermionFieldView &out)
 {
   assert(0);
 }
 
 template<class Impl> void 
-WilsonKernels<Impl >::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+WilsonKernels<Impl >::AsmDhopSiteExt(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-				     int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+				     int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, const FermionFieldView &out)
 {
   assert(0);
 }
 
 template<class Impl> void 
-WilsonKernels<Impl >::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+WilsonKernels<Impl >::AsmDhopSiteDagExt(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-					int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+					int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, const FermionFieldView &out)
 {
   assert(0);
 }

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

@@ -646,9 +646,14 @@ NAMESPACE_BEGIN(Grid);
   HAND_RESULT_EXT(ss,F)
 
 #define HAND_SPECIALISE_GPARITY(IMPL)					\
+  template<> accelerator_inline void					\
+  WilsonKernels<IMPL>::HandDhopSiteSycl(StencilVector st_perm, StencilEntry *st_p, \
+					SiteDoubledGaugeField *U, SiteHalfSpinor * buf, \
+					int sF, int sU, const SiteSpinor *in, SiteSpinor *out) {} \
+  									\
   template<> accelerator_inline void						\
-  WilsonKernels<IMPL>::HandDhopSite(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor  *buf, \
+  WilsonKernels<IMPL>::HandDhopSite(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor  *buf, \
-				    int ss,int sU,const FermionFieldView &in, FermionFieldView &out) \
+				    int ss,int sU,const FermionFieldView &in, const FermionFieldView &out) \
   {									\
     typedef IMPL Impl;							\
     typedef typename Simd::scalar_type S;				\
@@ -663,8 +668,8 @@ NAMESPACE_BEGIN(Grid);
   }									\
 									\
   template<> accelerator_inline void						\
-  WilsonKernels<IMPL>::HandDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf, \
+  WilsonKernels<IMPL>::HandDhopSiteDag(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor *buf, \
-				       int ss,int sU,const FermionFieldView &in, FermionFieldView &out) \
+				       int ss,int sU,const FermionFieldView &in, const FermionFieldView &out) \
   {									\
     typedef IMPL Impl;							\
     typedef typename Simd::scalar_type S;				\
@@ -679,8 +684,8 @@ NAMESPACE_BEGIN(Grid);
   }									\
 									\
   template<> accelerator_inline void						\
-  WilsonKernels<IMPL>::HandDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor  *buf, \
+  WilsonKernels<IMPL>::HandDhopSiteInt(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor  *buf, \
-				       int ss,int sU,const FermionFieldView &in, FermionFieldView &out) \
+				       int ss,int sU,const FermionFieldView &in, const FermionFieldView &out) \
   {									\
     typedef IMPL Impl;							\
     typedef typename Simd::scalar_type S;				\
@@ -695,8 +700,8 @@ NAMESPACE_BEGIN(Grid);
   }									\
 									\
   template<> accelerator_inline void						\
-  WilsonKernels<IMPL>::HandDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf, \
+  WilsonKernels<IMPL>::HandDhopSiteDagInt(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor *buf, \
-					  int ss,int sU,const FermionFieldView &in, FermionFieldView &out) \
+					  int ss,int sU,const FermionFieldView &in, const FermionFieldView &out) \
   {									\
     typedef IMPL Impl;							\
     typedef typename Simd::scalar_type S;				\
@@ -711,8 +716,8 @@ NAMESPACE_BEGIN(Grid);
   }									\
 									\
   template<> accelerator_inline void							\
-  WilsonKernels<IMPL>::HandDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor  *buf, \
+  WilsonKernels<IMPL>::HandDhopSiteExt(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor  *buf, \
-				       int ss,int sU,const FermionFieldView &in, FermionFieldView &out) \
+				       int ss,int sU,const FermionFieldView &in, const FermionFieldView &out) \
   {									\
     typedef IMPL Impl;							\
     typedef typename Simd::scalar_type S;				\
@@ -728,8 +733,8 @@ NAMESPACE_BEGIN(Grid);
     HAND_DOP_SITE_EXT(1, LOAD_CHI_GPARITY,LOAD_CHIMU_GPARITY,MULT_2SPIN_GPARITY); \
   }									\
   template<> accelerator_inline void						\
-  WilsonKernels<IMPL>::HandDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf, \
+  WilsonKernels<IMPL>::HandDhopSiteDagExt(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor *buf, \
-					  int ss,int sU,const FermionFieldView &in, FermionFieldView &out) \
+					  int ss,int sU,const FermionFieldView &in, const FermionFieldView &out) \
   {									\
     typedef IMPL Impl;							\
     typedef typename Simd::scalar_type S;				\

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

@@ -76,24 +76,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 #define REGISTER
 
-#ifdef GRID_SIMT
+#define LOAD_CHIMU \
-#define LOAD_CHIMU(ptype)		\
-  {const SiteSpinor & ref (in[offset]);	\
-    Chimu_00=coalescedReadPermute<ptype>(ref()(0)(0),perm,lane);	\
-    Chimu_01=coalescedReadPermute<ptype>(ref()(0)(1),perm,lane);		\
-    Chimu_02=coalescedReadPermute<ptype>(ref()(0)(2),perm,lane);		\
-    Chimu_10=coalescedReadPermute<ptype>(ref()(1)(0),perm,lane);		\
-    Chimu_11=coalescedReadPermute<ptype>(ref()(1)(1),perm,lane);		\
-    Chimu_12=coalescedReadPermute<ptype>(ref()(1)(2),perm,lane);		\
-    Chimu_20=coalescedReadPermute<ptype>(ref()(2)(0),perm,lane);		\
-    Chimu_21=coalescedReadPermute<ptype>(ref()(2)(1),perm,lane);		\
-    Chimu_22=coalescedReadPermute<ptype>(ref()(2)(2),perm,lane);		\
-    Chimu_30=coalescedReadPermute<ptype>(ref()(3)(0),perm,lane);		\
-    Chimu_31=coalescedReadPermute<ptype>(ref()(3)(1),perm,lane);		\
-    Chimu_32=coalescedReadPermute<ptype>(ref()(3)(2),perm,lane);	}
-#define PERMUTE_DIR(dir) ;
-#else
-#define LOAD_CHIMU(ptype)		\
   {const SiteSpinor & ref (in[offset]);	\
     Chimu_00=ref()(0)(0);\
     Chimu_01=ref()(0)(1);\
@@ -108,54 +91,54 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     Chimu_31=ref()(3)(1);\
     Chimu_32=ref()(3)(2);}
 
-#define PERMUTE_DIR(dir)			\
+#define LOAD_CHI\
-  permute##dir(Chi_00,Chi_00);	\
+  {const SiteHalfSpinor &ref(buf[offset]);	\
-      permute##dir(Chi_01,Chi_01);\
+    Chi_00 = ref()(0)(0);\
-      permute##dir(Chi_02,Chi_02);\
+    Chi_01 = ref()(0)(1);\
-      permute##dir(Chi_10,Chi_10);	\
+    Chi_02 = ref()(0)(2);\
-      permute##dir(Chi_11,Chi_11);\
+    Chi_10 = ref()(1)(0);\
-      permute##dir(Chi_12,Chi_12);
+    Chi_11 = ref()(1)(1);\
-
+    Chi_12 = ref()(1)(2);}
-#endif
 
+// To splat or not to splat depends on the implementation
 #define MULT_2SPIN(A)\
-  {auto & ref(U[sU](A));						\
+  {auto & ref(U[sU](A));			\
-    U_00=coalescedRead(ref()(0,0),lane);				\
+   Impl::loadLinkElement(U_00,ref()(0,0));	\
-    U_10=coalescedRead(ref()(1,0),lane);				\
+   Impl::loadLinkElement(U_10,ref()(1,0));	\
-    U_20=coalescedRead(ref()(2,0),lane);				\
+   Impl::loadLinkElement(U_20,ref()(2,0));	\
-    U_01=coalescedRead(ref()(0,1),lane);				\
+   Impl::loadLinkElement(U_01,ref()(0,1));	\
-    U_11=coalescedRead(ref()(1,1),lane);				\
+   Impl::loadLinkElement(U_11,ref()(1,1));	\
-    U_21=coalescedRead(ref()(2,1),lane);				\
+   Impl::loadLinkElement(U_21,ref()(2,1));	\
-    UChi_00 = U_00*Chi_00;						\
+    UChi_00 = U_00*Chi_00;\
-    UChi_10 = U_00*Chi_10;						\
+    UChi_10 = U_00*Chi_10;\
-    UChi_01 = U_10*Chi_00;						\
+    UChi_01 = U_10*Chi_00;\
-    UChi_11 = U_10*Chi_10;						\
+    UChi_11 = U_10*Chi_10;\
-    UChi_02 = U_20*Chi_00;						\
+    UChi_02 = U_20*Chi_00;\
-    UChi_12 = U_20*Chi_10;						\
+    UChi_12 = U_20*Chi_10;\
-    UChi_00+= U_01*Chi_01;						\
+    UChi_00+= U_01*Chi_01;\
-    UChi_10+= U_01*Chi_11;						\
+    UChi_10+= U_01*Chi_11;\
-    UChi_01+= U_11*Chi_01;						\
+    UChi_01+= U_11*Chi_01;\
-    UChi_11+= U_11*Chi_11;						\
+    UChi_11+= U_11*Chi_11;\
-    UChi_02+= U_21*Chi_01;						\
+    UChi_02+= U_21*Chi_01;\
-    UChi_12+= U_21*Chi_11;						\
+    UChi_12+= U_21*Chi_11;\
-    U_00=coalescedRead(ref()(0,2),lane);				\
+    Impl::loadLinkElement(U_00,ref()(0,2));	\
-    U_10=coalescedRead(ref()(1,2),lane);				\
+    Impl::loadLinkElement(U_10,ref()(1,2));	\
-    U_20=coalescedRead(ref()(2,2),lane);				\
+    Impl::loadLinkElement(U_20,ref()(2,2));	\
-    UChi_00+= U_00*Chi_02;						\
+    UChi_00+= U_00*Chi_02;\
-    UChi_10+= U_00*Chi_12;						\
+    UChi_10+= U_00*Chi_12;\
-    UChi_01+= U_10*Chi_02;						\
+    UChi_01+= U_10*Chi_02;\
-    UChi_11+= U_10*Chi_12;						\
+    UChi_11+= U_10*Chi_12;\
-    UChi_02+= U_20*Chi_02;						\
+    UChi_02+= U_20*Chi_02;\
     UChi_12+= U_20*Chi_12;}
 
-#define LOAD_CHI				\
+
-  {const SiteHalfSpinor &ref(buf[offset]);	\
+#define PERMUTE_DIR(dir)			\
-    Chi_00 = coalescedRead(ref()(0)(0),lane);	\
+      permute##dir(Chi_00,Chi_00);\
-    Chi_01 = coalescedRead(ref()(0)(1),lane);	\
+      permute##dir(Chi_01,Chi_01);\
-    Chi_02 = coalescedRead(ref()(0)(2),lane);	\
+      permute##dir(Chi_02,Chi_02);\
-    Chi_10 = coalescedRead(ref()(1)(0),lane);	\
+      permute##dir(Chi_10,Chi_10);\
-    Chi_11 = coalescedRead(ref()(1)(1),lane);	\
+      permute##dir(Chi_11,Chi_11);\
-    Chi_12 = coalescedRead(ref()(1)(2),lane);}
+      permute##dir(Chi_12,Chi_12);
 
 //      hspin(0)=fspin(0)+timesI(fspin(3));
 //      hspin(1)=fspin(1)+timesI(fspin(2));
@@ -370,13 +353,13 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
   result_31-= UChi_11;	\
   result_32-= UChi_12;
 
-#define HAND_STENCIL_LEGB(PROJ,PERM,DIR,RECON)	\
+#define HAND_STENCIL_LEG(PROJ,PERM,DIR,RECON)	\
   SE=st.GetEntry(ptype,DIR,ss);			\
   offset = SE->_offset;				\
   local  = SE->_is_local;			\
   perm   = SE->_permute;			\
   if ( local ) {				\
-    LOAD_CHIMU(PERM);				\
+    LOAD_CHIMU;					\
     PROJ;					\
     if ( perm) {				\
       PERMUTE_DIR(PERM);			\
@@ -384,37 +367,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
   } else {					\
     LOAD_CHI;					\
   }						\
-  acceleratorSynchronise();			\
-  MULT_2SPIN(DIR);				\
-  RECON;					
-
-#define HAND_STENCIL_LEG(PROJ,PERM,DIR,RECON)	\
-  SE=&st_p[DIR+8*ss];				\
-  ptype=st_perm[DIR];				\
-  offset = SE->_offset;				\
-  local  = SE->_is_local;			\
-  perm   = SE->_permute;			\
-  if ( local ) {				\
-    LOAD_CHIMU(PERM);				\
-    PROJ;					\
-    if ( perm) {				\
-      PERMUTE_DIR(PERM);			\
-    }						\
-  } else {					\
-    LOAD_CHI;					\
-  }						\
-  acceleratorSynchronise();			\
-  MULT_2SPIN(DIR);				\
-  RECON;					
-
-#define HAND_STENCIL_LEGA(PROJ,PERM,DIR,RECON)				\
-  SE=&st_p[DIR+8*ss];							\
-  ptype=st_perm[DIR];							\
- /*SE=st.GetEntry(ptype,DIR,ss);*/					\
-  offset = SE->_offset;				\
-  perm   = SE->_permute;			\
-  LOAD_CHIMU(PERM);				\
-  PROJ;						\
   MULT_2SPIN(DIR);				\
   RECON;					
 
@@ -424,7 +376,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
   local  = SE->_is_local;			\
   perm   = SE->_permute;			\
   if ( local ) {				\
-    LOAD_CHIMU(PERM);				\
+    LOAD_CHIMU;					\
     PROJ;					\
     if ( perm) {				\
       PERMUTE_DIR(PERM);			\
@@ -432,12 +384,10 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
   } else if ( st.same_node[DIR] ) {		\
     LOAD_CHI;					\
   }						\
-  acceleratorSynchronise();			\
   if (local || st.same_node[DIR] ) {		\
     MULT_2SPIN(DIR);				\
     RECON;					\
-  }						\
+  }
-  acceleratorSynchronise();			
 
 #define HAND_STENCIL_LEG_EXT(PROJ,PERM,DIR,RECON)	\
   SE=st.GetEntry(ptype,DIR,ss);			\
@@ -447,44 +397,44 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     MULT_2SPIN(DIR);				\
     RECON;					\
     nmu++;					\
-  }						\
+  }
-  acceleratorSynchronise();			
 
 #define HAND_RESULT(ss)				\
   {						\
-    SiteSpinor & ref (out[ss]);			\
+    SiteSpinor & ref (out[ss]);		\
-    coalescedWrite(ref()(0)(0),result_00,lane);		\
+    vstream(ref()(0)(0),result_00);		\
-    coalescedWrite(ref()(0)(1),result_01,lane);		\
+    vstream(ref()(0)(1),result_01);		\
-    coalescedWrite(ref()(0)(2),result_02,lane);		\
+    vstream(ref()(0)(2),result_02);		\
-    coalescedWrite(ref()(1)(0),result_10,lane);		\
+    vstream(ref()(1)(0),result_10);		\
-    coalescedWrite(ref()(1)(1),result_11,lane);		\
+    vstream(ref()(1)(1),result_11);		\
-    coalescedWrite(ref()(1)(2),result_12,lane);		\
+    vstream(ref()(1)(2),result_12);		\
-    coalescedWrite(ref()(2)(0),result_20,lane);		\
+    vstream(ref()(2)(0),result_20);		\
-    coalescedWrite(ref()(2)(1),result_21,lane);		\
+    vstream(ref()(2)(1),result_21);		\
-    coalescedWrite(ref()(2)(2),result_22,lane);		\
+    vstream(ref()(2)(2),result_22);		\
-    coalescedWrite(ref()(3)(0),result_30,lane);		\
+    vstream(ref()(3)(0),result_30);		\
-    coalescedWrite(ref()(3)(1),result_31,lane);		\
+    vstream(ref()(3)(1),result_31);		\
-    coalescedWrite(ref()(3)(2),result_32,lane);		\
+    vstream(ref()(3)(2),result_32);		\
   }
 
-#define HAND_RESULT_EXT(ss)				\
+#define HAND_RESULT_EXT(ss)			\
-  {							\
+  if (nmu){					\
-    SiteSpinor & ref (out[ss]);				\
+    SiteSpinor & ref (out[ss]);		\
-    coalescedWrite(ref()(0)(0),coalescedRead(ref()(0)(0))+result_00,lane);	\
+    ref()(0)(0)+=result_00;		\
-    coalescedWrite(ref()(0)(1),coalescedRead(ref()(0)(1))+result_01,lane);	\
+    ref()(0)(1)+=result_01;		\
-    coalescedWrite(ref()(0)(2),coalescedRead(ref()(0)(2))+result_02,lane);	\
+    ref()(0)(2)+=result_02;		\
-    coalescedWrite(ref()(1)(0),coalescedRead(ref()(1)(0))+result_10,lane);	\
+    ref()(1)(0)+=result_10;		\
-    coalescedWrite(ref()(1)(1),coalescedRead(ref()(1)(1))+result_11,lane);	\
+    ref()(1)(1)+=result_11;		\
-    coalescedWrite(ref()(1)(2),coalescedRead(ref()(1)(2))+result_12,lane);	\
+    ref()(1)(2)+=result_12;		\
-    coalescedWrite(ref()(2)(0),coalescedRead(ref()(2)(0))+result_20,lane);	\
+    ref()(2)(0)+=result_20;		\
-    coalescedWrite(ref()(2)(1),coalescedRead(ref()(2)(1))+result_21,lane);	\
+    ref()(2)(1)+=result_21;		\
-    coalescedWrite(ref()(2)(2),coalescedRead(ref()(2)(2))+result_22,lane);	\
+    ref()(2)(2)+=result_22;		\
-    coalescedWrite(ref()(3)(0),coalescedRead(ref()(3)(0))+result_30,lane);	\
+    ref()(3)(0)+=result_30;		\
-    coalescedWrite(ref()(3)(1),coalescedRead(ref()(3)(1))+result_31,lane);	\
+    ref()(3)(1)+=result_31;		\
-    coalescedWrite(ref()(3)(2),coalescedRead(ref()(3)(2))+result_32,lane);	\
+    ref()(3)(2)+=result_32;		\
   }
 
-#define HAND_DECLARATIONS(Simd)			\
+
+#define HAND_DECLARATIONS(a)			\
   Simd result_00;				\
   Simd result_01;				\
   Simd result_02;				\
@@ -516,19 +466,19 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
   Simd U_11;					\
   Simd U_21;
 
-#define ZERO_RESULT							\
+#define ZERO_RESULT				\
-  zeroit(result_00);							\
+  result_00=Zero();				\
-  zeroit(result_01);							\
+  result_01=Zero();				\
-  zeroit(result_02);							\
+  result_02=Zero();				\
-  zeroit(result_10);							\
+  result_10=Zero();				\
-  zeroit(result_11);							\
+  result_11=Zero();				\
-  zeroit(result_12);							\
+  result_12=Zero();				\
-  zeroit(result_20);							\
+  result_20=Zero();				\
-  zeroit(result_21);							\
+  result_21=Zero();				\
-  zeroit(result_22);							\
+  result_22=Zero();				\
-  zeroit(result_30);							\
+  result_30=Zero();				\
-  zeroit(result_31);							\
+  result_31=Zero();				\
-  zeroit(result_32);			
+  result_32=Zero();			
 
 #define Chimu_00 Chi_00
 #define Chimu_01 Chi_01
@@ -545,53 +495,15 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 NAMESPACE_BEGIN(Grid);
 
-
-#ifdef SYCL_HACK
 template<class Impl> accelerator_inline void 
-WilsonKernels<Impl>::HandDhopSiteSycl(StencilVector st_perm,StencilEntry *st_p, SiteDoubledGaugeField *U,SiteHalfSpinor  *buf,
+WilsonKernels<Impl>::HandDhopSite(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor  *buf,
-				      int ss,int sU,const SiteSpinor *in, SiteSpinor *out)
+				  int ss,int sU,const FermionFieldView &in, const FermionFieldView &out)
 {
 // T==0, Z==1, Y==2, Z==3 expect 1,2,2,2 simd layout etc...
   typedef typename Simd::scalar_type S;
   typedef typename Simd::vector_type V;
-  typedef iSinglet<Simd> vCplx;
-  //  typedef decltype( coalescedRead( vCplx()()() )) Simt;
-  typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt;
 
-  const int Nsimd = SiteHalfSpinor::Nsimd();
+  HAND_DECLARATIONS(ignore);
-  const int lane=acceleratorSIMTlane(Nsimd);
-
-  HAND_DECLARATIONS(Simt);
-
-  int offset,local,perm, ptype;
-  StencilEntry *SE;
-  HAND_STENCIL_LEG(XM_PROJ,3,Xp,XM_RECON);
-  HAND_STENCIL_LEG(YM_PROJ,2,Yp,YM_RECON_ACCUM);
-  HAND_STENCIL_LEG(ZM_PROJ,1,Zp,ZM_RECON_ACCUM);
-  HAND_STENCIL_LEG(TM_PROJ,0,Tp,TM_RECON_ACCUM);
-  HAND_STENCIL_LEG(XP_PROJ,3,Xm,XP_RECON_ACCUM);
-  HAND_STENCIL_LEG(YP_PROJ,2,Ym,YP_RECON_ACCUM);
-  HAND_STENCIL_LEG(ZP_PROJ,1,Zm,ZP_RECON_ACCUM);
-  HAND_STENCIL_LEG(TP_PROJ,0,Tm,TP_RECON_ACCUM);
-  HAND_RESULT(ss);
-}
-#endif
-
-template<class Impl> accelerator_inline void 
-WilsonKernels<Impl>::HandDhopSite(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor  *buf,
-				  int ss,int sU,const FermionFieldView &in, FermionFieldView &out)
-{
-  auto st_p = st._entries_p;						
-  auto st_perm = st._permute_type;					
-// T==0, Z==1, Y==2, Z==3 expect 1,2,2,2 simd layout etc...
-  typedef typename Simd::scalar_type S;
-  typedef typename Simd::vector_type V;
-  typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt;
-
-  const int Nsimd = SiteHalfSpinor::Nsimd();
-  const int lane=acceleratorSIMTlane(Nsimd);
-
-  HAND_DECLARATIONS(Simt);
 
   int offset,local,perm, ptype;
   StencilEntry *SE;
@@ -608,19 +520,13 @@ WilsonKernels<Impl>::HandDhopSite(StencilView &st, DoubledGaugeFieldView &U,Site
 }
 
 template<class Impl>  accelerator_inline
-void WilsonKernels<Impl>::HandDhopSiteDag(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+void WilsonKernels<Impl>::HandDhopSiteDag(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-					  int ss,int sU,const FermionFieldView &in, FermionFieldView &out)
+					  int ss,int sU,const FermionFieldView &in, const FermionFieldView &out)
 {
-  auto st_p = st._entries_p;						
-  auto st_perm = st._permute_type;					
   typedef typename Simd::scalar_type S;
   typedef typename Simd::vector_type V;
-  typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt;
 
-  const int Nsimd = SiteHalfSpinor::Nsimd();
+  HAND_DECLARATIONS(ignore);
-  const int lane=acceleratorSIMTlane(Nsimd);
-
-  HAND_DECLARATIONS(Simt);
 
   StencilEntry *SE;
   int offset,local,perm, ptype;
@@ -637,20 +543,14 @@ void WilsonKernels<Impl>::HandDhopSiteDag(StencilView &st,DoubledGaugeFieldView
 }
 
 template<class Impl>  accelerator_inline void 
-WilsonKernels<Impl>::HandDhopSiteInt(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor  *buf,
+WilsonKernels<Impl>::HandDhopSiteInt(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor  *buf,
-					  int ss,int sU,const FermionFieldView &in, FermionFieldView &out)
+					  int ss,int sU,const FermionFieldView &in, const FermionFieldView &out)
 {
-  auto st_p = st._entries_p;						
-  auto st_perm = st._permute_type;					
 // T==0, Z==1, Y==2, Z==3 expect 1,2,2,2 simd layout etc...
   typedef typename Simd::scalar_type S;
   typedef typename Simd::vector_type V;
-  typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt;
 
-  const int Nsimd = SiteHalfSpinor::Nsimd();
+  HAND_DECLARATIONS(ignore);
-  const int lane=acceleratorSIMTlane(Nsimd);
-
-  HAND_DECLARATIONS(Simt);
 
   int offset,local,perm, ptype;
   StencilEntry *SE;
@@ -667,19 +567,13 @@ WilsonKernels<Impl>::HandDhopSiteInt(StencilView &st,DoubledGaugeFieldView &U,Si
 }
 
 template<class Impl> accelerator_inline
-void WilsonKernels<Impl>::HandDhopSiteDagInt(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+void WilsonKernels<Impl>::HandDhopSiteDagInt(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-						  int ss,int sU,const FermionFieldView &in, FermionFieldView &out)
+						  int ss,int sU,const FermionFieldView &in, const FermionFieldView &out)
 {
-  auto st_p = st._entries_p;						
-  auto st_perm = st._permute_type;					
   typedef typename Simd::scalar_type S;
   typedef typename Simd::vector_type V;
-  typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt;
 
-  const int Nsimd = SiteHalfSpinor::Nsimd();
+  HAND_DECLARATIONS(ignore);
-  const int lane=acceleratorSIMTlane(Nsimd);
-
-  HAND_DECLARATIONS(Simt);
 
   StencilEntry *SE;
   int offset,local,perm, ptype;
@@ -696,20 +590,14 @@ void WilsonKernels<Impl>::HandDhopSiteDagInt(StencilView &st,DoubledGaugeFieldVi
 }
 
 template<class Impl>  accelerator_inline void 
-WilsonKernels<Impl>::HandDhopSiteExt(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor  *buf,
+WilsonKernels<Impl>::HandDhopSiteExt(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor  *buf,
-					  int ss,int sU,const FermionFieldView &in, FermionFieldView &out)
+					  int ss,int sU,const FermionFieldView &in, const FermionFieldView &out)
 {
-  auto st_p = st._entries_p;						
-  auto st_perm = st._permute_type;					
 // T==0, Z==1, Y==2, Z==3 expect 1,2,2,2 simd layout etc...
   typedef typename Simd::scalar_type S;
   typedef typename Simd::vector_type V;
-  typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt;
 
-  const int Nsimd = SiteHalfSpinor::Nsimd();
+  HAND_DECLARATIONS(ignore);
-  const int lane=acceleratorSIMTlane(Nsimd);
-
-  HAND_DECLARATIONS(Simt);
 
   int offset, ptype;
   StencilEntry *SE;
@@ -727,19 +615,13 @@ WilsonKernels<Impl>::HandDhopSiteExt(StencilView &st,DoubledGaugeFieldView &U,Si
 }
 
 template<class Impl>  accelerator_inline
-void WilsonKernels<Impl>::HandDhopSiteDagExt(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+void WilsonKernels<Impl>::HandDhopSiteDagExt(const StencilView &st,const DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-						  int ss,int sU,const FermionFieldView &in, FermionFieldView &out)
+						  int ss,int sU,const FermionFieldView &in, const FermionFieldView &out)
 {
-  auto st_p = st._entries_p;						
-  auto st_perm = st._permute_type;					
   typedef typename Simd::scalar_type S;
   typedef typename Simd::vector_type V;
-  typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt;
 
-  const int Nsimd = SiteHalfSpinor::Nsimd();
+  HAND_DECLARATIONS(ignore);
-  const int lane=acceleratorSIMTlane(Nsimd);
-
-  HAND_DECLARATIONS(Simt);
 
   StencilEntry *SE;
   int offset, ptype;

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

@@ -0,0 +1,598 @@
+   /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/qcd/action/fermion/WilsonKernelsHand.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 */
+
+#pragma once
+
+#include <Grid/qcd/action/fermion/FermionCore.h>
+
+
+#undef LOAD_CHIMU  
+#undef LOAD_CHI 
+#undef MULT_2SPIN
+#undef PERMUTE_DIR
+#undef XP_PROJ  
+#undef YP_PROJ  
+#undef ZP_PROJ  
+#undef TP_PROJ  
+#undef XM_PROJ  
+#undef YM_PROJ  
+#undef ZM_PROJ  
+#undef TM_PROJ  
+#undef XP_RECON 
+#undef XP_RECON_ACCUM 
+#undef XM_RECON 
+#undef XM_RECON_ACCUM 
+#undef YP_RECON_ACCUM 
+#undef YM_RECON_ACCUM 
+#undef ZP_RECON_ACCUM 
+#undef ZM_RECON_ACCUM 
+#undef TP_RECON_ACCUM 
+#undef TM_RECON_ACCUM 
+#undef ZERO_RESULT				 
+#undef Chimu_00
+#undef Chimu_01
+#undef Chimu_02
+#undef Chimu_10
+#undef Chimu_11
+#undef Chimu_12
+#undef Chimu_20
+#undef Chimu_21
+#undef Chimu_22
+#undef Chimu_30
+#undef Chimu_31
+#undef Chimu_32
+#undef HAND_STENCIL_LEG
+#undef HAND_STENCIL_LEG_INT
+#undef HAND_STENCIL_LEG_EXT
+#undef HAND_RESULT
+#undef HAND_RESULT_INT
+#undef HAND_RESULT_EXT
+
+#define REGISTER
+
+#ifdef GRID_SIMT
+#define LOAD_CHIMU(ptype)		\
+  {const SiteSpinor & ref (in[offset]);	\
+    Chimu_00=coalescedReadPermute<ptype>(ref()(0)(0),perm);	\
+    Chimu_01=coalescedReadPermute<ptype>(ref()(0)(1),perm);	\
+    Chimu_02=coalescedReadPermute<ptype>(ref()(0)(2),perm);	\
+    Chimu_10=coalescedReadPermute<ptype>(ref()(1)(0),perm);	\
+    Chimu_11=coalescedReadPermute<ptype>(ref()(1)(1),perm);	\
+    Chimu_12=coalescedReadPermute<ptype>(ref()(1)(2),perm);	\
+    Chimu_20=coalescedReadPermute<ptype>(ref()(2)(0),perm);	\
+    Chimu_21=coalescedReadPermute<ptype>(ref()(2)(1),perm);	\
+    Chimu_22=coalescedReadPermute<ptype>(ref()(2)(2),perm);	\
+    Chimu_30=coalescedReadPermute<ptype>(ref()(3)(0),perm);	\
+    Chimu_31=coalescedReadPermute<ptype>(ref()(3)(1),perm);	\
+    Chimu_32=coalescedReadPermute<ptype>(ref()(3)(2),perm);	}
+
+#define PERMUTE_DIR(dir) ;
+#else
+#define LOAD_CHIMU(ptype)		\
+  {const SiteSpinor & ref (in[offset]);	\
+    Chimu_00=coalescedRead(ref()(0)(0));	\
+    Chimu_01=coalescedRead(ref()(0)(1));	\
+    Chimu_02=coalescedRead(ref()(0)(2));	\
+    Chimu_10=coalescedRead(ref()(1)(0));	\
+    Chimu_11=coalescedRead(ref()(1)(1));	\
+    Chimu_12=coalescedRead(ref()(1)(2));	\
+    Chimu_20=coalescedRead(ref()(2)(0));	\
+    Chimu_21=coalescedRead(ref()(2)(1));	\
+    Chimu_22=coalescedRead(ref()(2)(2));	\
+    Chimu_30=coalescedRead(ref()(3)(0));	\
+    Chimu_31=coalescedRead(ref()(3)(1));	\
+    Chimu_32=coalescedRead(ref()(3)(2));	}
+
+#define PERMUTE_DIR(dir)			\
+  permute##dir(Chi_00,Chi_00);	\
+      permute##dir(Chi_01,Chi_01);\
+      permute##dir(Chi_02,Chi_02);\
+      permute##dir(Chi_10,Chi_10);	\
+      permute##dir(Chi_11,Chi_11);\
+      permute##dir(Chi_12,Chi_12);
+#endif
+
+#define MULT_2SPIN(A)\
+  {auto & ref(U[sU](A));					\
+  U_00=coalescedRead(ref()(0,0));				\
+  U_10=coalescedRead(ref()(1,0));					\
+  U_20=coalescedRead(ref()(2,0));					\
+  U_01=coalescedRead(ref()(0,1));					\
+  U_11=coalescedRead(ref()(1,1));					\
+  U_21=coalescedRead(ref()(2,1));					\
+    UChi_00 = U_00*Chi_00;					\
+    UChi_10 = U_00*Chi_10;					\
+    UChi_01 = U_10*Chi_00;					\
+    UChi_11 = U_10*Chi_10;					\
+    UChi_02 = U_20*Chi_00;					\
+    UChi_12 = U_20*Chi_10;					\
+    UChi_00+= U_01*Chi_01;					\
+    UChi_10+= U_01*Chi_11;					\
+    UChi_01+= U_11*Chi_01;					\
+    UChi_11+= U_11*Chi_11;					\
+    UChi_02+= U_21*Chi_01;					\
+    UChi_12+= U_21*Chi_11;					\
+    U_00=coalescedRead(ref()(0,2));				\
+    U_10=coalescedRead(ref()(1,2));				\
+    U_20=coalescedRead(ref()(2,2));				\
+    UChi_00+= U_00*Chi_02;					\
+    UChi_10+= U_00*Chi_12;					\
+    UChi_01+= U_10*Chi_02;					\
+    UChi_11+= U_10*Chi_12;					\
+    UChi_02+= U_20*Chi_02;					\
+    UChi_12+= U_20*Chi_12;}
+
+#define LOAD_CHI				\
+  {const SiteHalfSpinor &ref(buf[offset]);	\
+    Chi_00 = coalescedRead(ref()(0)(0));	\
+    Chi_01 = coalescedRead(ref()(0)(1));	\
+    Chi_02 = coalescedRead(ref()(0)(2));	\
+    Chi_10 = coalescedRead(ref()(1)(0));	\
+    Chi_11 = coalescedRead(ref()(1)(1));	\
+    Chi_12 = coalescedRead(ref()(1)(2));}
+
+//      hspin(0)=fspin(0)+timesI(fspin(3));
+//      hspin(1)=fspin(1)+timesI(fspin(2));
+#define XP_PROJ \
+    Chi_00 = Chimu_00+timesI(Chimu_30);\
+    Chi_01 = Chimu_01+timesI(Chimu_31);\
+    Chi_02 = Chimu_02+timesI(Chimu_32);\
+    Chi_10 = Chimu_10+timesI(Chimu_20);\
+    Chi_11 = Chimu_11+timesI(Chimu_21);\
+    Chi_12 = Chimu_12+timesI(Chimu_22);
+
+#define YP_PROJ \
+    Chi_00 = Chimu_00-Chimu_30;\
+    Chi_01 = Chimu_01-Chimu_31;\
+    Chi_02 = Chimu_02-Chimu_32;\
+    Chi_10 = Chimu_10+Chimu_20;\
+    Chi_11 = Chimu_11+Chimu_21;\
+    Chi_12 = Chimu_12+Chimu_22;
+
+#define ZP_PROJ \
+  Chi_00 = Chimu_00+timesI(Chimu_20);		\
+  Chi_01 = Chimu_01+timesI(Chimu_21);		\
+  Chi_02 = Chimu_02+timesI(Chimu_22);		\
+  Chi_10 = Chimu_10-timesI(Chimu_30);		\
+  Chi_11 = Chimu_11-timesI(Chimu_31);		\
+  Chi_12 = Chimu_12-timesI(Chimu_32);
+
+#define TP_PROJ \
+  Chi_00 = Chimu_00+Chimu_20;		\
+  Chi_01 = Chimu_01+Chimu_21;		\
+  Chi_02 = Chimu_02+Chimu_22;		\
+  Chi_10 = Chimu_10+Chimu_30;		\
+  Chi_11 = Chimu_11+Chimu_31;		\
+  Chi_12 = Chimu_12+Chimu_32;
+
+
+//      hspin(0)=fspin(0)-timesI(fspin(3));
+//      hspin(1)=fspin(1)-timesI(fspin(2));
+#define XM_PROJ \
+    Chi_00 = Chimu_00-timesI(Chimu_30);\
+    Chi_01 = Chimu_01-timesI(Chimu_31);\
+    Chi_02 = Chimu_02-timesI(Chimu_32);\
+    Chi_10 = Chimu_10-timesI(Chimu_20);\
+    Chi_11 = Chimu_11-timesI(Chimu_21);\
+    Chi_12 = Chimu_12-timesI(Chimu_22);
+
+#define YM_PROJ \
+    Chi_00 = Chimu_00+Chimu_30;\
+    Chi_01 = Chimu_01+Chimu_31;\
+    Chi_02 = Chimu_02+Chimu_32;\
+    Chi_10 = Chimu_10-Chimu_20;\
+    Chi_11 = Chimu_11-Chimu_21;\
+    Chi_12 = Chimu_12-Chimu_22;
+
+#define ZM_PROJ \
+  Chi_00 = Chimu_00-timesI(Chimu_20);		\
+  Chi_01 = Chimu_01-timesI(Chimu_21);		\
+  Chi_02 = Chimu_02-timesI(Chimu_22);		\
+  Chi_10 = Chimu_10+timesI(Chimu_30);		\
+  Chi_11 = Chimu_11+timesI(Chimu_31);		\
+  Chi_12 = Chimu_12+timesI(Chimu_32);
+
+#define TM_PROJ \
+  Chi_00 = Chimu_00-Chimu_20;		\
+  Chi_01 = Chimu_01-Chimu_21;		\
+  Chi_02 = Chimu_02-Chimu_22;		\
+  Chi_10 = Chimu_10-Chimu_30;		\
+  Chi_11 = Chimu_11-Chimu_31;		\
+  Chi_12 = Chimu_12-Chimu_32;
+
+//      fspin(0)=hspin(0);
+//      fspin(1)=hspin(1);
+//      fspin(2)=timesMinusI(hspin(1));
+//      fspin(3)=timesMinusI(hspin(0));
+#define XP_RECON\
+  result_00 = UChi_00;\
+  result_01 = UChi_01;\
+  result_02 = UChi_02;\
+  result_10 = UChi_10;\
+  result_11 = UChi_11;\
+  result_12 = UChi_12;\
+  result_20 = timesMinusI(UChi_10);\
+  result_21 = timesMinusI(UChi_11);\
+  result_22 = timesMinusI(UChi_12);\
+  result_30 = timesMinusI(UChi_00);\
+  result_31 = timesMinusI(UChi_01);\
+  result_32 = timesMinusI(UChi_02);
+
+#define XP_RECON_ACCUM\
+  result_00+=UChi_00;\
+  result_01+=UChi_01;\
+  result_02+=UChi_02;\
+  result_10+=UChi_10;\
+  result_11+=UChi_11;\
+  result_12+=UChi_12;\
+  result_20-=timesI(UChi_10);\
+  result_21-=timesI(UChi_11);\
+  result_22-=timesI(UChi_12);\
+  result_30-=timesI(UChi_00);\
+  result_31-=timesI(UChi_01);\
+  result_32-=timesI(UChi_02);
+
+#define XM_RECON\
+  result_00 = UChi_00;\
+  result_01 = UChi_01;\
+  result_02 = UChi_02;\
+  result_10 = UChi_10;\
+  result_11 = UChi_11;\
+  result_12 = UChi_12;\
+  result_20 = timesI(UChi_10);\
+  result_21 = timesI(UChi_11);\
+  result_22 = timesI(UChi_12);\
+  result_30 = timesI(UChi_00);\
+  result_31 = timesI(UChi_01);\
+  result_32 = timesI(UChi_02);
+
+#define XM_RECON_ACCUM\
+  result_00+= UChi_00;\
+  result_01+= UChi_01;\
+  result_02+= UChi_02;\
+  result_10+= UChi_10;\
+  result_11+= UChi_11;\
+  result_12+= UChi_12;\
+  result_20+= timesI(UChi_10);\
+  result_21+= timesI(UChi_11);\
+  result_22+= timesI(UChi_12);\
+  result_30+= timesI(UChi_00);\
+  result_31+= timesI(UChi_01);\
+  result_32+= timesI(UChi_02);
+
+#define YP_RECON_ACCUM\
+  result_00+= UChi_00;\
+  result_01+= UChi_01;\
+  result_02+= UChi_02;\
+  result_10+= UChi_10;\
+  result_11+= UChi_11;\
+  result_12+= UChi_12;\
+  result_20+= UChi_10;\
+  result_21+= UChi_11;\
+  result_22+= UChi_12;\
+  result_30-= UChi_00;\
+  result_31-= UChi_01;\
+  result_32-= UChi_02;
+
+#define YM_RECON_ACCUM\
+  result_00+= UChi_00;\
+  result_01+= UChi_01;\
+  result_02+= UChi_02;\
+  result_10+= UChi_10;\
+  result_11+= UChi_11;\
+  result_12+= UChi_12;\
+  result_20-= UChi_10;\
+  result_21-= UChi_11;\
+  result_22-= UChi_12;\
+  result_30+= UChi_00;\
+  result_31+= UChi_01;\
+  result_32+= UChi_02;
+
+#define ZP_RECON_ACCUM\
+  result_00+= UChi_00;\
+  result_01+= UChi_01;\
+  result_02+= UChi_02;\
+  result_10+= UChi_10;\
+  result_11+= UChi_11;\
+  result_12+= UChi_12;\
+  result_20-= timesI(UChi_00);			\
+  result_21-= timesI(UChi_01);			\
+  result_22-= timesI(UChi_02);			\
+  result_30+= timesI(UChi_10);			\
+  result_31+= timesI(UChi_11);			\
+  result_32+= timesI(UChi_12);
+
+#define ZM_RECON_ACCUM\
+  result_00+= UChi_00;\
+  result_01+= UChi_01;\
+  result_02+= UChi_02;\
+  result_10+= UChi_10;\
+  result_11+= UChi_11;\
+  result_12+= UChi_12;\
+  result_20+= timesI(UChi_00);			\
+  result_21+= timesI(UChi_01);			\
+  result_22+= timesI(UChi_02);			\
+  result_30-= timesI(UChi_10);			\
+  result_31-= timesI(UChi_11);			\
+  result_32-= timesI(UChi_12);
+
+#define TP_RECON_ACCUM\
+  result_00+= UChi_00;\
+  result_01+= UChi_01;\
+  result_02+= UChi_02;\
+  result_10+= UChi_10;\
+  result_11+= UChi_11;\
+  result_12+= UChi_12;\
+  result_20+= UChi_00;			\
+  result_21+= UChi_01;			\
+  result_22+= UChi_02;			\
+  result_30+= UChi_10;			\
+  result_31+= UChi_11;			\
+  result_32+= UChi_12;
+
+#define TM_RECON_ACCUM\
+  result_00+= UChi_00;\
+  result_01+= UChi_01;\
+  result_02+= UChi_02;\
+  result_10+= UChi_10;\
+  result_11+= UChi_11;\
+  result_12+= UChi_12;\
+  result_20-= UChi_00;	\
+  result_21-= UChi_01;	\
+  result_22-= UChi_02;	\
+  result_30-= UChi_10;	\
+  result_31-= UChi_11;	\
+  result_32-= UChi_12;
+
+#define HAND_STENCIL_LEGA(PROJ,PERM,DIR,RECON)	\
+  SE=&st_p[DIR+8*ss];			\
+  ptype=st_perm[DIR];			\
+  offset = SE->_offset;				\
+  local  = SE->_is_local;			\
+  perm   = SE->_permute;			\
+  if ( local ) {				\
+    LOAD_CHIMU(PERM);				\
+    PROJ;					\
+    if ( perm) {				\
+      PERMUTE_DIR(PERM);			\
+    }						\
+  } else {					\
+    LOAD_CHI;					\
+  }						\
+  MULT_2SPIN(DIR);				\
+  RECON;					
+
+#define HAND_STENCIL_LEG(PROJ,PERM,DIR,RECON)	\
+  SE=&st_p[DIR+8*ss];			\
+  ptype=st_perm[DIR];			\
+  offset = SE->_offset;				\
+  local  = SE->_is_local;			\
+  perm   = SE->_permute;			\
+  LOAD_CHIMU(PERM);				\
+  PROJ;						\
+  MULT_2SPIN(DIR);				\
+  RECON;					
+
+
+#define HAND_STENCIL_LEG_INT(PROJ,PERM,DIR,RECON)	\
+  SE=&st_p[DIR+8*ss];					\
+  ptype=st_perm[DIR];					\
+  offset = SE->_offset;				\
+  local  = SE->_is_local;			\
+  perm   = SE->_permute;			\
+  if ( local ) {				\
+    LOAD_CHIMU;					\
+    PROJ;					\
+    if ( perm) {				\
+      PERMUTE_DIR(PERM);			\
+    }						\
+  } else if ( st.same_node[DIR] ) {		\
+    LOAD_CHI;					\
+  }						\
+  if (local || st.same_node[DIR] ) {		\
+    MULT_2SPIN(DIR);				\
+    RECON;					\
+  }
+
+#define HAND_STENCIL_LEG_EXT(PROJ,PERM,DIR,RECON)	\
+  SE=st.GetEntry(ptype,DIR,ss);			\
+  offset = SE->_offset;				\
+  if((!SE->_is_local)&&(!st.same_node[DIR]) ) {	\
+    LOAD_CHI;					\
+    MULT_2SPIN(DIR);				\
+    RECON;					\
+    nmu++;					\
+  }
+
+#define HAND_RESULT(ss)				\
+  {						\
+    SiteSpinor & ref (out[ss]);			\
+    coalescedWrite(ref()(0)(0),result_00);		\
+    coalescedWrite(ref()(0)(1),result_01);		\
+    coalescedWrite(ref()(0)(2),result_02);		\
+    coalescedWrite(ref()(1)(0),result_10);		\
+    coalescedWrite(ref()(1)(1),result_11);		\
+    coalescedWrite(ref()(1)(2),result_12);		\
+    coalescedWrite(ref()(2)(0),result_20);		\
+    coalescedWrite(ref()(2)(1),result_21);		\
+    coalescedWrite(ref()(2)(2),result_22);		\
+    coalescedWrite(ref()(3)(0),result_30);		\
+    coalescedWrite(ref()(3)(1),result_31);		\
+    coalescedWrite(ref()(3)(2),result_32);		\
+  }
+
+#define HAND_RESULT_EXT(ss)			\
+  if (nmu){					\
+    SiteSpinor & ref (out[ss]);		\
+    ref()(0)(0)+=result_00;		\
+    ref()(0)(1)+=result_01;		\
+    ref()(0)(2)+=result_02;		\
+    ref()(1)(0)+=result_10;		\
+    ref()(1)(1)+=result_11;		\
+    ref()(1)(2)+=result_12;		\
+    ref()(2)(0)+=result_20;		\
+    ref()(2)(1)+=result_21;		\
+    ref()(2)(2)+=result_22;		\
+    ref()(3)(0)+=result_30;		\
+    ref()(3)(1)+=result_31;		\
+    ref()(3)(2)+=result_32;		\
+  }
+
+#define HAND_DECLARATIONS(Simd)			\
+  Simd result_00;				\
+  Simd result_01;				\
+  Simd result_02;				\
+  Simd result_10;				\
+  Simd result_11;				\
+  Simd result_12;				\
+  Simd result_20;				\
+  Simd result_21;				\
+  Simd result_22;				\
+  Simd result_30;				\
+  Simd result_31;				\
+  Simd result_32;				\
+  Simd Chi_00;					\
+  Simd Chi_01;					\
+  Simd Chi_02;					\
+  Simd Chi_10;					\
+  Simd Chi_11;					\
+  Simd Chi_12;					\
+  Simd UChi_00;					\
+  Simd UChi_01;					\
+  Simd UChi_02;					\
+  Simd UChi_10;					\
+  Simd UChi_11;					\
+  Simd UChi_12;					\
+  Simd U_00;					\
+  Simd U_10;					\
+  Simd U_20;					\
+  Simd U_01;					\
+  Simd U_11;					\
+  Simd U_21;
+
+#define ZERO_RESULT				\
+  result_00=Zero();				\
+  result_01=Zero();				\
+  result_02=Zero();				\
+  result_10=Zero();				\
+  result_11=Zero();				\
+  result_12=Zero();				\
+  result_20=Zero();				\
+  result_21=Zero();				\
+  result_22=Zero();				\
+  result_30=Zero();				\
+  result_31=Zero();				\
+  result_32=Zero();			
+
+#define Chimu_00 Chi_00
+#define Chimu_01 Chi_01
+#define Chimu_02 Chi_02
+#define Chimu_10 Chi_10
+#define Chimu_11 Chi_11
+#define Chimu_12 Chi_12
+#define Chimu_20 UChi_00
+#define Chimu_21 UChi_01
+#define Chimu_22 UChi_02
+#define Chimu_30 UChi_10
+#define Chimu_31 UChi_11
+#define Chimu_32 UChi_12
+
+NAMESPACE_BEGIN(Grid);
+
+template<class Impl> accelerator_inline void 
+WilsonKernels<Impl>::HandDhopSiteSycl(StencilVector st_perm,StencilEntry *st_p, SiteDoubledGaugeField *U,SiteHalfSpinor  *buf,
+				      int ss,int sU,const SiteSpinor *in, SiteSpinor *out)
+{
+// T==0, Z==1, Y==2, Z==3 expect 1,2,2,2 simd layout etc...
+  typedef typename Simd::scalar_type S;
+  typedef typename Simd::vector_type V;
+  typedef iSinglet<Simd> vCplx;
+  //  typedef decltype( coalescedRead( vCplx()()() )) Simt;
+  typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt;
+
+  HAND_DECLARATIONS(Simt);
+
+  int offset,local,perm, ptype;
+  StencilEntry *SE;
+  HAND_STENCIL_LEG(XM_PROJ,3,Xp,XM_RECON);
+  HAND_STENCIL_LEG(YM_PROJ,2,Yp,YM_RECON_ACCUM);
+  HAND_STENCIL_LEG(ZM_PROJ,1,Zp,ZM_RECON_ACCUM);
+  HAND_STENCIL_LEG(TM_PROJ,0,Tp,TM_RECON_ACCUM);
+  HAND_STENCIL_LEG(XP_PROJ,3,Xm,XP_RECON_ACCUM);
+  HAND_STENCIL_LEG(YP_PROJ,2,Ym,YP_RECON_ACCUM);
+  HAND_STENCIL_LEG(ZP_PROJ,1,Zm,ZP_RECON_ACCUM);
+  HAND_STENCIL_LEG(TP_PROJ,0,Tm,TP_RECON_ACCUM);
+  HAND_RESULT(ss);
+}
+
+////////////// Wilson ; uses this implementation /////////////////////
+
+NAMESPACE_END(Grid);
+#undef LOAD_CHIMU  
+#undef LOAD_CHI 
+#undef MULT_2SPIN
+#undef PERMUTE_DIR
+#undef XP_PROJ  
+#undef YP_PROJ  
+#undef ZP_PROJ  
+#undef TP_PROJ  
+#undef XM_PROJ  
+#undef YM_PROJ  
+#undef ZM_PROJ  
+#undef TM_PROJ  
+#undef XP_RECON 
+#undef XP_RECON_ACCUM 
+#undef XM_RECON 
+#undef XM_RECON_ACCUM 
+#undef YP_RECON_ACCUM 
+#undef YM_RECON_ACCUM 
+#undef ZP_RECON_ACCUM 
+#undef ZM_RECON_ACCUM 
+#undef TP_RECON_ACCUM 
+#undef TM_RECON_ACCUM 
+#undef ZERO_RESULT				 
+#undef Chimu_00
+#undef Chimu_01
+#undef Chimu_02
+#undef Chimu_10
+#undef Chimu_11
+#undef Chimu_12
+#undef Chimu_20
+#undef Chimu_21
+#undef Chimu_22
+#undef Chimu_30
+#undef Chimu_31
+#undef Chimu_32
+#undef HAND_STENCIL_LEG
+#undef HAND_STENCIL_LEG_INT
+#undef HAND_STENCIL_LEG_EXT
+#undef HAND_RESULT
+#undef HAND_RESULT_INT
+#undef HAND_RESULT_EXT
+#undef HAND_DECLARATIONS

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

@@ -115,9 +115,9 @@ accelerator_inline void get_stencil(StencilEntry * mem, StencilEntry &chip)
   // All legs kernels ; comms then compute
   ////////////////////////////////////////////////////////////////////
 template <class Impl> accelerator_inline
-void WilsonKernels<Impl>::GenericDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U,
+void WilsonKernels<Impl>::GenericDhopSiteDag(const StencilView &st, const DoubledGaugeFieldView &U,
-					     SiteHalfSpinor *buf, int sF,
+ 					     SiteHalfSpinor *buf, int sF,
-					     int sU, const FermionFieldView &in, FermionFieldView &out)
+					     int sU, const FermionFieldView &in, const FermionFieldView &out)
 {
   typedef decltype(coalescedRead(buf[0]))   calcHalfSpinor;
   typedef decltype(coalescedRead(in[0])) calcSpinor;
@@ -141,9 +141,9 @@ void WilsonKernels<Impl>::GenericDhopSiteDag(StencilView &st, DoubledGaugeFieldV
 };
 
 template <class Impl> accelerator_inline
-void WilsonKernels<Impl>::GenericDhopSite(StencilView &st, DoubledGaugeFieldView &U,
+void WilsonKernels<Impl>::GenericDhopSite(const StencilView &st, const DoubledGaugeFieldView &U,
-					  SiteHalfSpinor *buf, int sF,
+ 					  SiteHalfSpinor *buf, int sF,
-					  int sU, const FermionFieldView &in, FermionFieldView &out)
+					  int sU, const FermionFieldView &in, const FermionFieldView &out)
 {
   typedef decltype(coalescedRead(buf[0])) calcHalfSpinor;
   typedef decltype(coalescedRead(in[0]))  calcSpinor;
@@ -170,9 +170,9 @@ void WilsonKernels<Impl>::GenericDhopSite(StencilView &st, DoubledGaugeFieldView
   // Interior kernels
   ////////////////////////////////////////////////////////////////////
 template <class Impl> accelerator_inline
-void WilsonKernels<Impl>::GenericDhopSiteDagInt(StencilView &st,  DoubledGaugeFieldView &U,
+void WilsonKernels<Impl>::GenericDhopSiteDagInt(const StencilView &st, const DoubledGaugeFieldView &U,
-						SiteHalfSpinor *buf, int sF,
+ 						SiteHalfSpinor *buf, int sF,
-						int sU, const FermionFieldView &in, FermionFieldView &out)
+						int sU, const FermionFieldView &in, const FermionFieldView &out)
 {
   typedef decltype(coalescedRead(buf[0])) calcHalfSpinor;
   typedef decltype(coalescedRead(in[0]))  calcSpinor;
@@ -198,9 +198,9 @@ void WilsonKernels<Impl>::GenericDhopSiteDagInt(StencilView &st,  DoubledGaugeFi
 };
 
 template <class Impl> accelerator_inline
-void WilsonKernels<Impl>::GenericDhopSiteInt(StencilView &st,  DoubledGaugeFieldView &U,
+void WilsonKernels<Impl>::GenericDhopSiteInt(const StencilView &st, const DoubledGaugeFieldView &U,
-							 SiteHalfSpinor *buf, int sF,
+					     SiteHalfSpinor *buf, int sF,
-							 int sU, const FermionFieldView &in, FermionFieldView &out)
+					     int sU, const FermionFieldView &in, const FermionFieldView &out)
 {
   typedef decltype(coalescedRead(buf[0])) calcHalfSpinor;
   typedef decltype(coalescedRead(in[0]))  calcSpinor;
@@ -228,9 +228,9 @@ void WilsonKernels<Impl>::GenericDhopSiteInt(StencilView &st,  DoubledGaugeField
 // Exterior kernels
 ////////////////////////////////////////////////////////////////////
 template <class Impl> accelerator_inline
-void WilsonKernels<Impl>::GenericDhopSiteDagExt(StencilView &st,  DoubledGaugeFieldView &U,
+void WilsonKernels<Impl>::GenericDhopSiteDagExt(const StencilView &st, const DoubledGaugeFieldView &U,
-						SiteHalfSpinor *buf, int sF,
+ 						SiteHalfSpinor *buf, int sF,
-						int sU, const FermionFieldView &in, FermionFieldView &out)
+						int sU, const FermionFieldView &in, const FermionFieldView &out)
 {
   typedef decltype(coalescedRead(buf[0])) calcHalfSpinor;
   typedef decltype(coalescedRead(in[0]))  calcSpinor;
@@ -259,9 +259,9 @@ void WilsonKernels<Impl>::GenericDhopSiteDagExt(StencilView &st,  DoubledGaugeFi
 };
 
 template <class Impl> accelerator_inline
-void WilsonKernels<Impl>::GenericDhopSiteExt(StencilView &st,  DoubledGaugeFieldView &U,
+void WilsonKernels<Impl>::GenericDhopSiteExt(const StencilView &st, const DoubledGaugeFieldView &U,
-					     SiteHalfSpinor *buf, int sF,
+ 					     SiteHalfSpinor *buf, int sF,
-					     int sU, const FermionFieldView &in, FermionFieldView &out)
+					     int sU, const FermionFieldView &in, const FermionFieldView &out)
 {
   typedef decltype(coalescedRead(buf[0])) calcHalfSpinor;
   typedef decltype(coalescedRead(in[0]))  calcSpinor;
@@ -291,8 +291,8 @@ void WilsonKernels<Impl>::GenericDhopSiteExt(StencilView &st,  DoubledGaugeField
 
 #define DhopDirMacro(Dir,spProj,spRecon)	\
   template <class Impl> accelerator_inline				\
-  void WilsonKernels<Impl>::DhopDir##Dir(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf, int sF, \
+  void WilsonKernels<Impl>::DhopDir##Dir(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor *buf, int sF, \
-					 int sU, const FermionFieldView &in, FermionFieldView &out, int dir) \
+					 int sU, const FermionFieldView &in, const FermionFieldView &out, int dir) \
   {									\
   typedef decltype(coalescedRead(buf[0])) calcHalfSpinor;		\
   typedef decltype(coalescedRead(in[0]))  calcSpinor;			\
@@ -302,8 +302,8 @@ void WilsonKernels<Impl>::GenericDhopSiteExt(StencilView &st,  DoubledGaugeField
   StencilEntry *SE;							\
   int ptype;								\
   const int Nsimd = SiteHalfSpinor::Nsimd();				\
-  const int lane=acceleratorSIMTlane(Nsimd);					\
+  const int lane=acceleratorSIMTlane(Nsimd);				\
-									\
+  									\
   SE = st.GetEntry(ptype, dir, sF);					\
   GENERIC_DHOPDIR_LEG_BODY(Dir,spProj,spRecon);				\
   coalescedWrite(out[sF], result,lane);					\
@@ -319,8 +319,8 @@ DhopDirMacro(Zm,spProjZm,spReconZm);
 DhopDirMacro(Tm,spProjTm,spReconTm);
 
 template <class Impl> accelerator_inline
-void WilsonKernels<Impl>::DhopDirK( StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf, int sF,
+void WilsonKernels<Impl>::DhopDirK(const StencilView &st, const DoubledGaugeFieldView &U,SiteHalfSpinor *buf, int sF,
-				    int sU, const FermionFieldView &in, FermionFieldView &out, int dir, int gamma)
+				   int sU, const FermionFieldView &in, const FermionFieldView &out, int dir, int gamma)
 {
   typedef decltype(coalescedRead(buf[0])) calcHalfSpinor;
   typedef decltype(coalescedRead(in[0]))  calcSpinor;
@@ -345,8 +345,8 @@ void WilsonKernels<Impl>::DhopDirK( StencilView &st, DoubledGaugeFieldView &U,Si
 }
 
 template <class Impl>
-void WilsonKernels<Impl>::DhopDirAll( StencilImpl &st, DoubledGaugeField &U,SiteHalfSpinor *buf, int Ls,
+void WilsonKernels<Impl>::DhopDirAll(StencilImpl &st, DoubledGaugeField &U,SiteHalfSpinor *buf, int Ls,
-				      int Nsite, const FermionField &in, std::vector<FermionField> &out)
+				     int Nsite, const FermionField &in, std::vector<FermionField> &out)
 {
    autoView(U_v  ,U,AcceleratorRead);
    autoView(in_v ,in,AcceleratorRead);
@@ -416,6 +416,14 @@ void WilsonKernels<Impl>::DhopDirKernel( StencilImpl &st, DoubledGaugeField &U,S
 #undef LoopBody
 }
 
+#define KERNEL_CALLNB(A) \
+  const uint64_t    NN = Nsite*Ls;					\
+  accelerator_forNB( ss, NN, Simd::Nsimd(), {				\
+      int sF = ss;							\
+      int sU = ss/Ls;							\
+      WilsonKernels<Impl>::A(st_v,U_v,buf,sF,sU,in_v,out_v);		\
+  });
+
 #define KERNEL_CALL_TMP(A) \
   const uint64_t    NN = Nsite*Ls;					\
   auto U_p = & U_v[0];							\
@@ -430,14 +438,6 @@ void WilsonKernels<Impl>::DhopDirKernel( StencilImpl &st, DoubledGaugeField &U,S
     });									\
   accelerator_barrier();
 
-#define KERNEL_CALLNB(A)						\
-  const uint64_t    NN = Nsite*Ls;					\
-  accelerator_forNB( ss, NN, Simd::Nsimd(), {				\
-      int sF = ss;							\
-      int sU = ss/Ls;							\
-      WilsonKernels<Impl>::A(st_v,U_v,buf,sF,sU,in_v,out_v);		\
-  });
-
 #define KERNEL_CALL(A) KERNEL_CALLNB(A); accelerator_barrier();
 
 #define ASM_CALL(A)							\
@@ -459,24 +459,21 @@ void WilsonKernels<Impl>::DhopKernel(int Opt,StencilImpl &st,  DoubledGaugeField
 
    if( interior && exterior ) {
      if (Opt == WilsonKernelsStatic::OptGeneric    ) { KERNEL_CALL(GenericDhopSite); return;}
-#ifdef SYCL_HACK     
-     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL_TMP(HandDhopSiteSycl);    return; }
-#else
-     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSite);    return;}
-#endif     
 #ifndef GRID_CUDA
+     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL_TMP(HandDhopSiteSycl);    return; }
+     //     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSite);    return;}
      if (Opt == WilsonKernelsStatic::OptInlineAsm  ) {  ASM_CALL(AsmDhopSite);    return;}
 #endif
    } else if( interior ) {
      if (Opt == WilsonKernelsStatic::OptGeneric    ) { KERNEL_CALLNB(GenericDhopSiteInt); return;}
-     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALLNB(HandDhopSiteInt);    return;}
 #ifndef GRID_CUDA
+     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALLNB(HandDhopSiteInt);    return;}
      if (Opt == WilsonKernelsStatic::OptInlineAsm  ) {  ASM_CALL(AsmDhopSiteInt);    return;}
 #endif
    } else if( exterior ) {
      if (Opt == WilsonKernelsStatic::OptGeneric    ) { KERNEL_CALL(GenericDhopSiteExt); return;}
-     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteExt);    return;}
 #ifndef GRID_CUDA
+     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteExt);    return;}
      if (Opt == WilsonKernelsStatic::OptInlineAsm  ) {  ASM_CALL(AsmDhopSiteExt);    return;}
 #endif
    }
@@ -494,20 +491,20 @@ void WilsonKernels<Impl>::DhopKernel(int Opt,StencilImpl &st,  DoubledGaugeField
 
    if( interior && exterior ) {
      if (Opt == WilsonKernelsStatic::OptGeneric    ) { KERNEL_CALL(GenericDhopSiteDag); return;}
-     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDag);    return;}
 #ifndef GRID_CUDA
+     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDag);    return;}
      if (Opt == WilsonKernelsStatic::OptInlineAsm  ) {  ASM_CALL(AsmDhopSiteDag);     return;}
 #endif
    } else if( interior ) {
      if (Opt == WilsonKernelsStatic::OptGeneric    ) { KERNEL_CALL(GenericDhopSiteDagInt); return;}
-     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDagInt);    return;}
 #ifndef GRID_CUDA
+     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDagInt);    return;}
      if (Opt == WilsonKernelsStatic::OptInlineAsm  ) {  ASM_CALL(AsmDhopSiteDagInt);     return;}
 #endif
    } else if( exterior ) {
      if (Opt == WilsonKernelsStatic::OptGeneric    ) { KERNEL_CALL(GenericDhopSiteDagExt); return;}
-     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDagExt);    return;}
 #ifndef GRID_CUDA
+     if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDagExt);    return;}
      if (Opt == WilsonKernelsStatic::OptInlineAsm  ) {  ASM_CALL(AsmDhopSiteDagExt);     return;}
 #endif
    }

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

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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/WilsonAdjImplD/WilsonKernelsInstantiationWilsonAdjImplD.cc

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

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

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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

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

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

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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/WilsonImplD/WilsonKernelsInstantiationWilsonImplD.cc

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

Grid/qcd/action/fermion/instantiation/WilsonImplDF/WilsonKernelsInstantiationWilsonImplDF.cc

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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/WilsonImplDF/WilsonKernelsInstantiationWilsonImplDF.cc

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

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

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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/WilsonKernelsInstantiationWilsonImplF.cc

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

Grid/qcd/action/fermion/instantiation/WilsonImplFH/WilsonKernelsInstantiationWilsonImplFH.cc

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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/WilsonImplFH/WilsonKernelsInstantiationWilsonImplFH.cc

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

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

@@ -32,6 +32,7 @@ directory
 #include <Grid/qcd/action/fermion/FermionCore.h>
 #include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
 #include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementationSycl.h>
 
 #ifndef AVX512
 #ifndef QPX

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

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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

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

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

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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

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

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

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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

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

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

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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

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

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

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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

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

Grid/qcd/action/fermion/instantiation/ZWilsonImplDF/WilsonKernelsInstantiationZWilsonImplDF.cc

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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/ZWilsonImplDF/WilsonKernelsInstantiationZWilsonImplDF.cc

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

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

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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

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

Grid/qcd/action/fermion/instantiation/ZWilsonImplFH/WilsonKernelsInstantiationZWilsonImplFH.cc

@@ -1,51 +0,0 @@
-/*************************************************************************************
-
-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/ZWilsonImplFH/WilsonKernelsInstantiationZWilsonImplFH.cc

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

Grid/qcd/action/gauge/GaugeImplTypes.h

@@ -96,7 +96,7 @@ public:
   ///////////////////////////////////////////////////////////
   // Move these to another class
   // HMC auxiliary functions
-  static inline void generate_momenta(Field &P, GridSerialRNG & sRNG, GridParallelRNG &pRNG) 
+  static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) 
   {
     // Zbigniew Srocinsky thesis:
     //

Grid/qcd/action/gauge/PlaqPlusRectangleAction.h

@@ -49,7 +49,7 @@ public:
 
   virtual std::string action_name(){return "PlaqPlusRectangleAction";}
       
-  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {}; // noop as no pseudoferms
+  virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {}; // noop as no pseudoferms
       
   virtual std::string LogParameters(){
     std::stringstream sstream;

Grid/qcd/action/gauge/WilsonGaugeAction.h

@@ -54,7 +54,8 @@ public:
     return sstream.str();
   }
 
-  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG){};  // noop as no pseudoferms
+  virtual void refresh(const GaugeField &U,
+                       GridParallelRNG &pRNG){};  // noop as no pseudoferms
 
   virtual RealD S(const GaugeField &U) {
     RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);

Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h

@@ -124,7 +124,7 @@ NAMESPACE_BEGIN(Grid);
       //
       // As a check of rational require \Phi^dag M_{EOFA} \Phi == eta^dag M^-1/2^dag M M^-1/2 eta = eta^dag eta
       //
-      virtual void refresh(const GaugeField& U, GridSerialRNG &sRNG, GridParallelRNG& pRNG)
+      virtual void refresh(const GaugeField& U, GridParallelRNG& pRNG)
       {
         Lop.ImportGauge(U);
         Rop.ImportGauge(U);

Grid/qcd/action/pseudofermion/OneFlavourEvenOddRational.h

@@ -1,3 +1,4 @@
+
 /*************************************************************************************
 
 Grid physics library, www.github.com/paboyle/Grid
@@ -42,7 +43,8 @@ NAMESPACE_BEGIN(Grid);
 //
 
 template <class Impl>
-class OneFlavourEvenOddRationalPseudoFermionAction : public Action<typename Impl::GaugeField> {
+class OneFlavourEvenOddRationalPseudoFermionAction
+  : public Action<typename Impl::GaugeField> {
 public:
   INHERIT_IMPL_TYPES(Impl);
 
@@ -101,7 +103,7 @@ public:
     return sstream.str();
   }
   
-  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG) {
+  virtual void refresh(const GaugeField &U, GridParallelRNG &pRNG) {
     // P(phi) = e^{- phi^dag (MpcdagMpc)^-1/2 phi}
     //        = e^{- phi^dag (MpcdagMpc)^-1/4 (MpcdagMpc)^-1/4 phi}
     // Phi = MpcdagMpc^{1/4} eta
@@ -154,10 +156,7 @@ public:
 
     msCG(Mpc, PhiOdd, Y);
 
-    auto grid = FermOp.FermionGrid();
+    if ( (rand()%param.BoundsCheckFreq)==0 ) { 
-    auto r=rand();
-    grid->Broadcast(0,r);
-    if ( (r%param.BoundsCheckFreq)==0 ) { 
       FermionField gauss(FermOp.FermionRedBlackGrid());
       gauss = PhiOdd;
       HighBoundCheck(Mpc,gauss,param.hi);

Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h

@@ -101,7 +101,7 @@ NAMESPACE_BEGIN(Grid);
       }
       
       
-      virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
+      virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
 
 	// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi       
 	//
@@ -170,10 +170,7 @@ NAMESPACE_BEGIN(Grid);
 	msCG_M(MdagM,X,Y);
 
 	// Randomly apply rational bounds checks.
-	auto grid = NumOp.FermionGrid();
+	if ( (rand()%param.BoundsCheckFreq)==0 ) { 
-        auto r=rand();
-        grid->Broadcast(0,r);
-        if ( (r%param.BoundsCheckFreq)==0 ) { 
 	  FermionField gauss(NumOp.FermionRedBlackGrid());
 	  gauss = PhiOdd;
 	  HighBoundCheck(MdagM,gauss,param.hi);

Grid/qcd/action/pseudofermion/OneFlavourRational.h

@@ -98,7 +98,7 @@ NAMESPACE_BEGIN(Grid);
 
 
       
-      virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
+      virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
 
 	
 	// P(phi) = e^{- phi^dag (MdagM)^-1/2 phi}
@@ -142,10 +142,7 @@ NAMESPACE_BEGIN(Grid);
 
 	msCG(MdagMOp,Phi,Y);
 
-	auto grid = FermOp.FermionGrid();
+	if ( (rand()%param.BoundsCheckFreq)==0 ) { 
-        auto r=rand();
-        grid->Broadcast(0,r);
-        if ( (r%param.BoundsCheckFreq)==0 ) { 
 	  FermionField gauss(FermOp.FermionGrid());
 	  gauss = Phi;
 	  HighBoundCheck(MdagMOp,gauss,param.hi);

Grid/qcd/action/pseudofermion/OneFlavourRationalRatio.h

@@ -95,7 +95,7 @@ NAMESPACE_BEGIN(Grid);
       }
       
 
-      virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
+      virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
 
 	// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi       
 	//
@@ -156,10 +156,7 @@ NAMESPACE_BEGIN(Grid);
 	msCG_M(MdagM,X,Y);
 
 	// Randomly apply rational bounds checks.
-        auto grid = NumOp.FermionGrid();
+	if ( (rand()%param.BoundsCheckFreq)==0 ) { 
-        auto r=rand();
-        grid->Broadcast(0,r);
-        if ( (r%param.BoundsCheckFreq)==0 ) { 	
 	  FermionField gauss(NumOp.FermionGrid());
 	  gauss = Phi;
 	  HighBoundCheck(MdagM,gauss,param.hi);

Grid/qcd/action/pseudofermion/TwoFlavour.h

@@ -73,7 +73,7 @@ public:
   //////////////////////////////////////////////////////////////////////////////////////
   // Push the gauge field in to the dops. Assume any BC's and smearing already applied
   //////////////////////////////////////////////////////////////////////////////////////
-  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG) {
+  virtual void refresh(const GaugeField &U, GridParallelRNG &pRNG) {
     // P(phi) = e^{- phi^dag (MdagM)^-1 phi}
     // Phi = Mdag eta
     // P(eta) = e^{- eta^dag eta}

Grid/qcd/action/pseudofermion/TwoFlavourEvenOdd.h

@@ -77,7 +77,7 @@ public:
   //////////////////////////////////////////////////////////////////////////////////////
   // Push the gauge field in to the dops. Assume any BC's and smearing already applied
   //////////////////////////////////////////////////////////////////////////////////////
-  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
+  virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
     
     // P(phi) = e^{- phi^dag (MpcdagMpc)^-1 phi}
     // Phi = McpDag eta 

Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h

@@ -84,7 +84,7 @@ NAMESPACE_BEGIN(Grid);
       } 
 
       
-      virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
+      virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
 
         // P(phi) = e^{- phi^dag Vpc (MpcdagMpc)^-1 Vpcdag phi}
         //

Grid/qcd/action/pseudofermion/TwoFlavourRatio.h

@@ -64,7 +64,7 @@ public:
     return sstream.str();
   }  
       
-  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
+  virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
 
     // P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi}
     //

Grid/qcd/action/scalar/ScalarAction.h

@@ -55,7 +55,7 @@ public:
   }
   virtual std::string action_name() {return "ScalarAction";}
 
-  virtual void refresh(const Field &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG) {}  // noop as no pseudoferms
+  virtual void refresh(const Field &U, GridParallelRNG &pRNG) {}  // noop as no pseudoferms
 
   virtual RealD S(const Field &p) {
     return (mass_square * 0.5 + Nd) * ScalarObs<Impl>::sumphisquared(p) +

Grid/qcd/action/scalar/ScalarImpl.h

@@ -27,7 +27,7 @@ public:
   typedef Field              FermionField;
   typedef Field              PropagatorField;
     
-  static inline void generate_momenta(Field& P, GridSerialRNG &sRNG, GridParallelRNG& pRNG){
+  static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){
     RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR); // CPS/UKQCD momentum rescaling
     gaussian(pRNG, P);
     P *= scale; 
@@ -151,7 +151,7 @@ public:
       out = one / out;
     }
 
-    static inline void generate_momenta(Field &P, GridSerialRNG & sRNG, GridParallelRNG &pRNG)
+    static inline void generate_momenta(Field &P, GridParallelRNG &pRNG)
     {
       RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR); // CPS/UKQCD momentum rescaling
 #ifndef USE_FFT_ACCELERATION

Grid/qcd/action/scalar/ScalarInteractionAction.h

@@ -77,7 +77,7 @@ public:
 
   virtual std::string action_name() { return "ScalarAction"; }
 
-  virtual void refresh(const Field &U, GridSerialRNG & sRNG, GridParallelRNG &pRNG) {}
+  virtual void refresh(const Field &U, GridParallelRNG &pRNG) {}
 
   virtual RealD S(const Field &p)
   {

Grid/qcd/hmc/HMC.h

@@ -139,7 +139,7 @@ private:
   // Evolution
   /////////////////////////////////////////////////////////
   RealD evolve_hmc_step(Field &U) {
-    TheIntegrator.refresh(U, sRNG, pRNG);  // set U and initialize P and phi's
+    TheIntegrator.refresh(U, pRNG);  // set U and initialize P and phi's
 
     RealD H0 = TheIntegrator.S(U);  // initial state action
 

Grid/qcd/hmc/integrators/Integrator.h

@@ -33,7 +33,6 @@ directory
 #define INTEGRATOR_INCLUDED
 
 #include <memory>
-#include "MomentumFilter.h"
 
 NAMESPACE_BEGIN(Grid);
 
@@ -79,19 +78,8 @@ protected:
   RepresentationPolicy Representations;
   IntegratorParameters Params;
 
-  //Filters allow the user to manipulate the conjugate momentum, for example to freeze links in DDHMC
-  //It is applied whenever the momentum is updated / refreshed
-  //The default filter does nothing
-  MomentumFilterBase<MomentaField> const* MomFilter;
-
   const ActionSet<Field, RepresentationPolicy> as;
 
-  //Get a pointer to a shared static instance of the "do-nothing" momentum filter to serve as a default
-  static MomentumFilterBase<MomentaField> const* getDefaultMomFilter(){ 
-    static MomentumFilterNone<MomentaField> filter;
-    return &filter;
-  }
-
   void update_P(Field& U, int level, double ep) 
   {
     t_P[level] += ep;
@@ -147,8 +135,6 @@ protected:
 
     // Force from the other representations
     as[level].apply(update_P_hireps, Representations, Mom, U, ep);
-
-    MomFilter->applyFilter(Mom);
   }
 
   void update_U(Field& U, double ep) 
@@ -188,23 +174,11 @@ public:
     t_P.resize(levels, 0.0);
     t_U = 0.0;
     // initialization of smearer delegated outside of Integrator
-
-    //Default the momentum filter to "do-nothing"
-    MomFilter = getDefaultMomFilter();
   };
 
   virtual ~Integrator() {}
 
   virtual std::string integrator_name() = 0;
-  
-  //Set the momentum filter allowing for manipulation of the conjugate momentum
-  void setMomentumFilter(const MomentumFilterBase<MomentaField> &filter){
-    MomFilter = &filter;
-  }
-
-  //Access the conjugate momentum
-  const MomentaField & getMomentum() const{ return P; }
-  
 
   void print_parameters()
   {
@@ -236,9 +210,10 @@ public:
   // over the representations
   struct _refresh {
     template <class FieldType, class Repr>
-    void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep, GridSerialRNG & sRNG, GridParallelRNG& pRNG) {
+    void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep,
+                    GridParallelRNG& pRNG) {
       for (int a = 0; a < repr_set.size(); ++a){
-        repr_set.at(a)->refresh(Rep.U, sRNG, pRNG);
+        repr_set.at(a)->refresh(Rep.U, pRNG);
       
 	std::cout << GridLogDebug << "Hirep refreshing pseudofermions" << std::endl;
       }
@@ -246,12 +221,12 @@ public:
   } refresh_hireps{};
 
   // Initialization of momenta and actions
-  void refresh(Field& U,  GridSerialRNG & sRNG, GridParallelRNG& pRNG) 
+  void refresh(Field& U, GridParallelRNG& pRNG) 
   {
     assert(P.Grid() == U.Grid());
     std::cout << GridLogIntegrator << "Integrator refresh\n";
 
-    FieldImplementation::generate_momenta(P, sRNG, pRNG);
+    FieldImplementation::generate_momenta(P, pRNG);
 
     // Update the smeared fields, can be implemented as observer
     // necessary to keep the fields updated even after a reject
@@ -268,14 +243,12 @@ public:
         // get gauge field from the SmearingPolicy and
         // based on the boolean is_smeared in actionID
         Field& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
-        as[level].actions.at(actionID)->refresh(Us, sRNG, pRNG);
+        as[level].actions.at(actionID)->refresh(Us, pRNG);
       }
 
       // Refresh the higher representation actions
-      as[level].apply(refresh_hireps, Representations, sRNG, pRNG);
+      as[level].apply(refresh_hireps, Representations, pRNG);
     }
-
-    MomFilter->applyFilter(P);
   }
 
   // to be used by the actionlevel class to iterate

Grid/qcd/hmc/integrators/MomentumFilter.h

@@ -1,94 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/hmc/integrators/MomentumFilter.h
-
-Copyright (C) 2015
-
-Author: Christopher Kelly <ckelly@bnl.gov>
-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 */
-//--------------------------------------------------------------------
-#ifndef MOMENTUM_FILTER
-#define MOMENTUM_FILTER
-
-NAMESPACE_BEGIN(Grid);
-
-//These filter objects allow the user to manipulate the conjugate momentum as part of the update / refresh
-
-template<typename MomentaField>
-struct MomentumFilterBase{
-  virtual void applyFilter(MomentaField &P) const;
-};
-
-//Do nothing
-template<typename MomentaField>
-struct MomentumFilterNone: public MomentumFilterBase<MomentaField>{
-  void applyFilter(MomentaField &P) const override{}
-};
-
-//Multiply each site/direction by a Lorentz vector complex number field
-//Can be used to implement a mask, zeroing out sites
-template<typename MomentaField>
-struct MomentumFilterApplyPhase: public MomentumFilterBase<MomentaField>{
-  typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type
-  typedef typename MomentaField::scalar_type scalar_type; //scalar complex type
-  typedef iVector<iScalar<iScalar<vector_type> >, Nd > LorentzScalarType; //complex phase for each site/direction
-  typedef Lattice<LorentzScalarType> LatticeLorentzScalarType;
-  
-  LatticeLorentzScalarType phase;
- 
-  MomentumFilterApplyPhase(const LatticeLorentzScalarType _phase): phase(_phase){}
-
-  //Default to uniform field of (1,0)
-  MomentumFilterApplyPhase(GridBase* _grid): phase(_grid){
-    LorentzScalarType one;
-    for(int mu=0;mu<Nd;mu++)
-      one(mu)()() = scalar_type(1.);
-    
-    phase = one;
-  }
-
-  void applyFilter(MomentaField &P) const override{
-    conformable(P,phase);
-    autoView( P_v , P, AcceleratorWrite);
-    autoView( phase_v , phase, AcceleratorRead);
-
-    accelerator_for(ss,P_v.size(),MomentaField::vector_type::Nsimd(),{
-    	auto site_mom = P_v(ss);
-    	auto site_phase = phase_v(ss);
-	for(int mu=0;mu<Nd;mu++)
-	  site_mom(mu) = site_mom(mu) * site_phase(mu);
-    	coalescedWrite(P_v[ss], site_mom);
-      });
-    
-  }
-
-
-};
-
-
-
-
-NAMESPACE_END(Grid);
-
-#endif

Grid/qcd/smearing/StoutSmearing.h

@@ -85,18 +85,21 @@ public:
 
     std::cout << GridLogDebug << "Stout smearing started\n";
 
-    // C contains the staples multiplied by some rho
+    // Smear the configurations
-    u_smr = U ; // set the smeared field to the current gauge field
     SmearBase->smear(C, U);
 
     for (int mu = 0; mu < Nd; mu++) {
-      if( mu == OrthogDim ) continue ;
+      if( mu == OrthogDim )
-      // u_smr = exp(iQ_mu)*U_mu apart from Orthogdim
+        tmp = 1.0;  // Don't smear in the orthogonal direction
-      Umu = peekLorentz(U, mu);
+      else {
-      tmp = peekLorentz(C, mu);
+        tmp = peekLorentz(C, mu);
-      iq_mu = Ta( tmp * adj(Umu));  
+        Umu = peekLorentz(U, mu);
-      exponentiate_iQ(tmp, iq_mu);
+        iq_mu = Ta(
-      pokeLorentz(u_smr, tmp * Umu, mu);
+                   tmp *
+                   adj(Umu));  // iq_mu = Ta(Omega_mu) to match the signs with the paper
+        exponentiate_iQ(tmp, iq_mu);
+      }
+      pokeLorentz(u_smr, tmp * Umu, mu);  // u_smr = exp(iQ_mu)*U_mu
     }
     std::cout << GridLogDebug << "Stout smearing completed\n";
   };

Grid/qcd/utils/Metric.h

@@ -93,13 +93,13 @@ public:
   GeneralisedMomenta(GridBase* grid, Metric<MomentaField>& M): M(M), Mom(grid), AuxMom(grid), AuxField(grid){}
 
   // Correct
-  void MomentaDistribution(GridSerialRNG & sRNG, GridParallelRNG& pRNG){
+  void MomentaDistribution(GridParallelRNG& pRNG){
     // Generate a distribution for
     // P^dag G P
     // where G = M^-1
 
     // Generate gaussian momenta
-    Implementation::generate_momenta(Mom, sRNG, pRNG);
+    Implementation::generate_momenta(Mom, pRNG);
     // Modify the distribution with the metric
     M.MSquareRoot(Mom);
 
@@ -107,8 +107,8 @@ public:
       // Auxiliary momenta
       // do nothing if trivial, so hide in the metric
       MomentaField AuxMomTemp(Mom.Grid());
-      Implementation::generate_momenta(AuxMom, sRNG, pRNG);
+      Implementation::generate_momenta(AuxMom, pRNG);
-      Implementation::generate_momenta(AuxField, sRNG, pRNG);
+      Implementation::generate_momenta(AuxField, pRNG);
       // Modify the distribution with the metric
       // Aux^dag M Aux
       M.MInvSquareRoot(AuxMom);  // AuxMom = M^{-1/2} AuxMomTemp

Grid/simd/Grid_gpu_vec.h

@@ -67,7 +67,6 @@ public:
   accelerator_inline GpuComplex(const GpuComplex &zz) { z = zz.z;};
   accelerator_inline Real real(void) const { return z.x; };
   accelerator_inline Real imag(void) const { return z.y; };
-  accelerator_inline GpuComplex &operator=(const Zero &zz) { z.x = 0; z.y=0; return *this; };
   accelerator_inline GpuComplex &operator*=(const GpuComplex &r) {
     *this = (*this) * r;
     return *this;

Grid/simd/Grid_vector_types.h

@@ -208,8 +208,8 @@ struct RealPart<complex<T> > {
 //////////////////////////////////////
 // type alias used to simplify the syntax of std::enable_if
 template <typename T> using Invoke = typename T::type;
-template <typename Condition, typename ReturnType = void> using EnableIf    = Invoke<std::enable_if<Condition::value, ReturnType> >;
+template <typename Condition, typename ReturnType> using EnableIf    = Invoke<std::enable_if<Condition::value, ReturnType> >;
-template <typename Condition, typename ReturnType = void> using NotEnableIf = Invoke<std::enable_if<!Condition::value, ReturnType> >;
+template <typename Condition, typename ReturnType> using NotEnableIf = Invoke<std::enable_if<!Condition::value, ReturnType> >;
 
 ////////////////////////////////////////////////////////
 // Check for complexity with type traits

Grid/stencil/Stencil.h

@@ -221,7 +221,7 @@ public:
   typedef typename cobj::vector_type vector_type;
   typedef typename cobj::scalar_type scalar_type;
   typedef typename cobj::scalar_object scalar_object;
-  typedef const CartesianStencilView<vobj,cobj,Parameters> View_type;
+  typedef CartesianStencilView<vobj,cobj,Parameters> View_type;
   typedef typename View_type::StencilVector StencilVector;
   ///////////////////////////////////////////
   // Helper structs

Grid/tensors/Tensor_SIMT.h

@@ -65,9 +65,8 @@ void coalescedWriteNonTemporal(vobj & __restrict__ vec,const vobj & __restrict__
 #else
 
 
-//#ifndef GRID_SYCL
+#if 0
-#if 1
+// Use the scalar as our own complex on GPU
-// Use the scalar as our own complex on GPU ... thrust::complex or std::complex
 template<class vsimd,IfSimd<vsimd> = 0> accelerator_inline
 typename vsimd::scalar_type
 coalescedRead(const vsimd & __restrict__ vec,int lane=acceleratorSIMTlane(vsimd::Nsimd()))
@@ -97,8 +96,6 @@ void coalescedWrite(vsimd & __restrict__ vec,
   p[lane]=extracted;
 }
 #else
-// For SyCL have option to use GpuComplex from inside the vector type in SIMT loops
-// Faster for some reason
 template<class vsimd,IfSimd<vsimd> = 0> accelerator_inline
 typename vsimd::vector_type::datum
 coalescedRead(const vsimd & __restrict__ vec,int lane=acceleratorSIMTlane(vsimd::Nsimd()))

Grid/tensors/Tensor_exp.h

@@ -28,7 +28,7 @@ Author: neo <cossu@post.kek.jp>
 #ifndef GRID_MATH_EXP_H
 #define GRID_MATH_EXP_H
 
-#define DEFAULT_MAT_EXP 20
+#define DEFAULT_MAT_EXP 12
 
 NAMESPACE_BEGIN(Grid);
 

Grid/tensors/Tensor_outer.h

@@ -34,16 +34,6 @@ NAMESPACE_BEGIN(Grid);
 // outerProduct Scalar x Scalar -> Scalar
 //              Vector x Vector -> Matrix
 ///////////////////////////////////////////////////////////////////////////////////////
-template<class CC,IfComplex<CC> = 0>
-accelerator_inline CC outerProduct(const CC &l, const CC& r)
-{
-  return l*conj(r);
-}
-template<class RR,IfReal<RR> = 0>
-accelerator_inline RR outerProduct(const RR &l, const RR& r)
-{
-  return l*r;
-}
 
 template<class l,class r,int N> accelerator_inline
 auto outerProduct (const iVector<l,N>& lhs,const iVector<r,N>& rhs) -> iMatrix<decltype(outerProduct(lhs._internal[0],rhs._internal[0])),N>
@@ -67,6 +57,17 @@ auto outerProduct (const iScalar<l>& lhs,const iScalar<r>& rhs) -> iScalar<declt
   return ret;
 }
 
+template<class CC,IfComplex<CC> = 0>
+accelerator_inline CC outerProduct(const CC &l, const CC& r)
+{
+  return l*conj(r);
+}
+template<class RR,IfReal<RR> = 0>
+accelerator_inline RR outerProduct(const RR &l, const RR& r)
+{
+  return l*r;
+}
+
 NAMESPACE_END(Grid);
 
 #endif

Grid/threads/Accelerator.cc

@@ -53,6 +53,7 @@ void acceleratorInit(void)
     prop = gpu_props[i];
     totalDeviceMem = prop.totalGlobalMem;
     if ( world_rank == 0) {
+#ifndef GRID_DEFAULT_GPU
       if ( i==rank ) {
 	printf("AcceleratorCudaInit[%d]: ========================\n",rank);
 	printf("AcceleratorCudaInit[%d]: Device Number    : %d\n", rank,i);
@@ -66,8 +67,8 @@ void acceleratorInit(void)
 	GPU_PROP(warpSize);
 	GPU_PROP(pciBusID);
 	GPU_PROP(pciDeviceID);
- 	printf("AcceleratorCudaInit[%d]: maxGridSize (%d,%d,%d)\n",rank,prop.maxGridSize[0],prop.maxGridSize[1],prop.maxGridSize[2]);
       }
+#endif
       //      GPU_PROP(unifiedAddressing);
       //      GPU_PROP(l2CacheSize);
       //      GPU_PROP(singleToDoublePrecisionPerfRatio);

Grid/threads/Accelerator.h

@@ -104,7 +104,7 @@ extern int acceleratorAbortOnGpuError;
 
 accelerator_inline int acceleratorSIMTlane(int Nsimd) {
 #ifdef GRID_SIMT
-  return threadIdx.x; 
+  return threadIdx.z; 
 #else
   return 0;
 #endif
@@ -112,76 +112,36 @@ accelerator_inline int acceleratorSIMTlane(int Nsimd) {
 
 #define accelerator_for2dNB( iter1, num1, iter2, num2, nsimd, ... )	\
   {									\
-    int nt=acceleratorThreads();					\
     typedef uint64_t Iterator;						\
     auto lambda = [=] accelerator					\
       (Iterator iter1,Iterator iter2,Iterator lane) mutable {		\
       __VA_ARGS__;							\
     };									\
-    dim3 cu_threads(nsimd,acceleratorThreads(),1);			\
+    int nt=acceleratorThreads();					\
+    dim3 cu_threads(acceleratorThreads(),1,nsimd);			\
     dim3 cu_blocks ((num1+nt-1)/nt,num2,1);				\
     LambdaApply<<<cu_blocks,cu_threads>>>(num1,num2,nsimd,lambda);	\
   }
 
-#define accelerator_for6dNB(iter1, num1,				\
-                            iter2, num2,				\
-                            iter3, num3,				\
-                            iter4, num4,				\
-                            iter5, num5,				\
-			    iter6, num6, ... )				\
-  {									\
-    typedef uint64_t Iterator;						\
-    auto lambda = [=] accelerator					\
-      (Iterator iter1,Iterator iter2,					\
-       Iterator iter3,Iterator iter4,					\
-       Iterator iter5,Iterator iter6) mutable {				\
-      __VA_ARGS__;							\
-    };									\
-    dim3 cu_blocks (num1,num2,num3);					\
-    dim3 cu_threads(num4,num5,num6);					\
-    Lambda6Apply<<<cu_blocks,cu_threads>>>(num1,num2,num3,num4,num5,num6,lambda); \
-  }
-
 template<typename lambda>  __global__
 void LambdaApply(uint64_t num1, uint64_t num2, uint64_t num3, lambda Lambda)
 {
-  // Weird permute is to make lane coalesce for large blocks
+  uint64_t x = threadIdx.x + blockDim.x*blockIdx.x;
-  uint64_t x = threadIdx.y + blockDim.y*blockIdx.x;
+  uint64_t y = threadIdx.y + blockDim.y*blockIdx.y;
-  uint64_t y = threadIdx.z + blockDim.z*blockIdx.y;
+  uint64_t z = threadIdx.z;
-  uint64_t z = threadIdx.x;
   if ( (x < num1) && (y<num2) && (z<num3) ) {
     Lambda(x,y,z);
   }
 }
 
-template<typename lambda>  __global__
-void Lambda6Apply(uint64_t num1, uint64_t num2, uint64_t num3,
-		  uint64_t num4, uint64_t num5, uint64_t num6,
-		  lambda Lambda)
-{
-  uint64_t iter1 = blockIdx.x;
-  uint64_t iter2 = blockIdx.y;
-  uint64_t iter3 = blockIdx.z;
-  uint64_t iter4 = threadIdx.x;
-  uint64_t iter5 = threadIdx.y;
-  uint64_t iter6 = threadIdx.z;
-
-  if ( (iter1 < num1) && (iter2<num2) && (iter3<num3)
-    && (iter4 < num4) && (iter5<num5) && (iter6<num6) )
-  {
-    Lambda(iter1,iter2,iter3,iter4,iter5,iter6);
-  }
-}
-
 #define accelerator_barrier(dummy)					\
   {									\
     cudaDeviceSynchronize();						\
     cudaError err = cudaGetLastError();					\
     if ( cudaSuccess != err ) {						\
-      printf("accelerator_barrier(): Cuda error %s \n",			\
+      printf("Cuda error %s \n", cudaGetErrorString( err ));		\
-	     cudaGetErrorString( err ));				\
+      puts(__FILE__);							\
-      printf("File %s Line %d\n",__FILE__,__LINE__);			\
+      printf("Line %d\n",__LINE__);					\
-      fflush(stdout);							\
       if (acceleratorAbortOnGpuError) assert(err==cudaSuccess);		\
     }									\
   }
@@ -457,7 +417,7 @@ accelerator_inline void acceleratorSynchronise(void)
   __syncwarp();
 #endif
 #ifdef GRID_SYCL
-  //cl::sycl::detail::workGroupBarrier();
+  // No barrier call on SYCL??  // Option get __spir:: stuff to do warp barrier
 #endif
 #ifdef GRID_HIP
   __syncthreads();

Grid/util/CompilerCompatible.h

@@ -1,16 +1,5 @@
 #pragma once 
 
-#if defined(__NVCC__)
-
-#if (__CUDACC_VER_MAJOR__ == 11) && (__CUDACC_VER_MINOR__ == 0)
-#error "NVCC version 11.0 breaks on Ampere, see Github issue 346"
-#endif
-#if (__CUDACC_VER_MAJOR__ == 11) && (__CUDACC_VER_MINOR__ == 1)
-#error "NVCC version 11.1 breaks on Ampere, see Github issue 346"
-#endif
-
-#endif
-
 #if defined(__clang__)
 
   #if __clang_major__ < 3

Grid/util/Init.cc

@@ -140,7 +140,7 @@ void GridCmdOptionCSL(std::string str,std::vector<std::string> & vec)
 }
 
 template<class VectorInt>
-void GridCmdOptionIntVector(const std::string &str,VectorInt & vec)
+void GridCmdOptionIntVector(std::string &str,VectorInt & vec)
 {
   vec.resize(0);
   std::stringstream ss(str);
@@ -153,9 +153,6 @@ void GridCmdOptionIntVector(const std::string &str,VectorInt & vec)
   return;
 }
 
-template void GridCmdOptionIntVector(const std::string &str,std::vector<int> & vec);
-template void GridCmdOptionIntVector(const std::string &str,Coordinate & vec);
-
 void GridCmdOptionInt(std::string &str,int & val)
 {
   std::stringstream ss(str);

Grid/util/Init.h

@@ -55,7 +55,7 @@ template<class VectorInt>
 std::string GridCmdVectorIntToString(const VectorInt & vec);
 void GridCmdOptionCSL(std::string str,std::vector<std::string> & vec);
 template<class VectorInt>
-void GridCmdOptionIntVector(const std::string &str,VectorInt & vec);
+void GridCmdOptionIntVector(std::string &str,VectorInt & vec);
 void GridCmdOptionInt(std::string &str,int & val);
 
 

HMC/Mobius2p1fRHMC.cc

@@ -56,12 +56,12 @@ int main(int argc, char **argv) {
   MD.trajL   = 1.0;
 
   HMCparameters HMCparams;
-  HMCparams.StartTrajectory  = 0;
+  HMCparams.StartTrajectory  = 30;
   HMCparams.Trajectories     = 200;
   HMCparams.NoMetropolisUntil=  0;
   // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
-  HMCparams.StartingType     =std::string("ColdStart");
+  //  HMCparams.StartingType     =std::string("ColdStart");
-  //  HMCparams.StartingType     =std::string("CheckpointStart");
+  HMCparams.StartingType     =std::string("CheckpointStart");
   HMCparams.MD = MD;
   HMCWrapper TheHMC(HMCparams);
 

README.md

@@ -1,4 +1,4 @@
-# Grid [![Teamcity status](http://ci.cliath.ph.ed.ac.uk/app/rest/builds/aggregated/strob:(buildType:(affectedProject(id:GridBasedSoftware_Grid)),branch:name:develop)/statusIcon.svg)](http://ci.cliath.ph.ed.ac.uk/project.html?projectId=GridBasedSoftware_Grid&tab=projectOverview) 
+# Grid [![Teamcity status](http://ci.cliath.ph.ed.ac.uk/app/rest/builds/aggregated/strob:(buildType:(affectedProject(id:GridBasedSoftware_Grid)),branch:name:develop)/statusIcon.svg)](http://ci.cliath.ph.ed.ac.uk/project.html?projectId=GridBasedSoftware_Grid&tab=projectOverview) [![Travis status](https://travis-ci.org/paboyle/Grid.svg?branch=develop)](https://travis-ci.org/paboyle/Grid)
 
 **Data parallel C++ mathematical object library.**
 
@@ -149,6 +149,7 @@ If you want to build all the tests at once just use `make tests`.
 - `--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: 32 bytes).
+- `--enable-precision={single|double}`: set the default precision (default: `double`). **Deprecated option**
 - `--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.

TODO

@@ -1,6 +1,3 @@
--- comms threads issue??
--- Part done: Staggered kernel performance on GPU
-
 =========================================================
 General
 =========================================================
@@ -8,18 +5,28 @@ General
 - Make representations code take Gimpl
 - Simplify the HMCand remove modules
 - Lattice_arith - are the mult, mac etc.. still needed after ET engine?
-- Lattice_rng - faster local only loop in init
+- Lattice_rng
-- Audit: accelerate A2Autils -- off critical path for HMC
+- Lattice_transfer.h
+- accelerate A2Autils -- off critical path for HMC
 
 =========================================================
-GPU  work list
+GPU branch code item work list
 =========================================================
 
-* sum_cpu promote to double during summation for increased precision.
+* sum_cpu promote to double during summation for increased precisoin.
 * Introduce sumD & ReduceD 
 * GPU sum is probably better currently.
+
 * Accelerate the cshift & benchmark
 
+* 0) Single GPU
+- 128 bit integer table load in GPU code.
+  - ImprovedStaggered accelerate & measure perf
+  - Gianluca's changes to Cayley into gpu-port
+  - Mobius kernel fusion.                     -- Gianluca?
+  - Lebesque order reintroduction. StencilView should have pointer to it
+  - Lebesgue reorder in all kernels
+
 * 3) Comms/NVlink
 - OpenMP tasks to run comms threads. Experiment with it 
 - Remove explicit openMP in staggered. 
@@ -28,6 +35,14 @@ GPU  work list
 - Stencil gather ??
 - SIMD dirs in stencil
 
+* 4) ET enhancements
+- eval -> scalar ops in ET engine
+- coalescedRead, coalescedWrite in expressions.
+
+* 5) Misc
+- Conserved current clean up.
+- multLinkProp eliminate
+
 8) Merge develop and test HMC
 
 9) Gamma tables on GPU; check this. Appear to work, but no idea why. Are these done on CPU?
@@ -37,7 +52,7 @@ GPU  work list
 -     Audit NAMESPACE CHANGES
 -     Audit changes
 
----------
+-----
 Gianluca's changes
 - Performance impact of construct in aligned allocator???
 ---------
@@ -47,33 +62,6 @@ Gianluca's changes
 -----------------------------
 DONE:
 -----------------------------
-=====
--- Done: Remez X^-1/2 X^-1/2 X = 1 test.
-         Feed in MdagM^2 as a test and take its sqrt.
-         Automated test that MdagM invsqrt(MdagM)invsqrt(MdagM) = 1 in HMC for bounds satisfaction.
-
--- Done: Sycl Kernels into develop. Compare to existing unroll and just use.
--- Done: sRNG into refresh functions
--- Done: Tuned decomposition on CUDA into develop
--- Done: Sycl friend accessor. Const view attempt via typedef??
-
-
-* Done 5) Misc
-- Conserved current clean up.
-- multLinkProp eliminate
-
-* Done 0) Single GPU
-- 128 bit integer table load in GPU code.
-  - ImprovedStaggered accelerate & measure perf
-  - Gianluca's changes to Cayley into gpu-port
-  - Mobius kernel fusion.                     -- Gianluca?
-  - Lebesque order reintroduction. StencilView should have pointer to it
-  - Lebesgue reorder in all kernels
-
-* 4) ET enhancements
-- Done eval -> scalar ops in ET engine
-- Done coalescedRead, coalescedWrite in expressions.
-
 =============================================================================================
 AUDIT ContractWWVV with respect to develop    -- DONE
 - GPU accelerate EOFA                                                  -- DONE
@@ -137,6 +125,23 @@ AUDIT ContractWWVV with respect to develop    -- DONE
 - -      (4) omp parallel for collapse(n)
 - - Only (1) has a natural mirror in accelerator_loop
 - - Nested loop macros get cumbersome made a generic interface for N deep
+- - Don't like thread_region and thread_loop_in_region
+- - Could replace with 
+
+    thread_nested(1, 
+      for {
+
+      }
+    );
+    thread_nested(2,
+      for (){
+        for (){
+
+	}
+      }
+    );
+
+    and same "in_region".
 
 
 -----------------------------

benchmarks/Benchmark_dwf_fp32.cc

@@ -53,7 +53,7 @@ int main (int argc, char ** argv)
   int threads = GridThread::GetThreads();
 
   Coordinate latt4 = GridDefaultLatt();
-  int Ls=16;
+  int Ls=8;
   for(int i=0;i<argc;i++)
     if(std::string(argv[i]) == "-Ls"){
       std::stringstream ss(argv[i+1]); ss >> Ls;

configure.ac

@@ -7,12 +7,7 @@ AM_INIT_AUTOMAKE([subdir-objects 1.13])
 AM_EXTRA_RECURSIVE_TARGETS([tests bench])
 AC_CONFIG_MACRO_DIR([m4])
 AC_CONFIG_SRCDIR([Grid/Grid.h])
-AC_CONFIG_HEADERS([Grid/Config.h],[[$SED_INPLACE -e 's|PACKAGE_|GRID_|' -e 's|[[:space:]]PACKAGE[[:space:]]| GRID_PACKAGE |' -e 's|[[:space:]]VERSION[[:space:]]| GRID_PACKAGE_VERSION |' Grid/Config.h]],
+AC_CONFIG_HEADERS([Grid/Config.h],[sed -i 's|PACKAGE_|GRID_|' Grid/Config.h])
-    [if test x"$host_os" == x"${host_os#darwin}" ; then]
-        [SED_INPLACE="sed -i"]
-    [else]
-        [SED_INPLACE="sed -i .bak"]
-    [fi])
 m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
 
 ################ Get git info
@@ -130,7 +125,7 @@ esac
 
 ############### fermions
 AC_ARG_ENABLE([fermion-reps],
-     [AC_HELP_STRING([--enable-fermion-reps=yes|no], [enable extra fermion representation support])],
+     [AC_HELP_STRING([--fermion-reps=yes|no], [enable extra fermion representation support])],
      [ac_FERMION_REPS=${enable_fermion_reps}], [ac_FERMION_REPS=yes])
 
 AM_CONDITIONAL(BUILD_FERMION_REPS, [ test "${ac_FERMION_REPS}X" == "yesX" ])
@@ -140,23 +135,12 @@ AC_ARG_ENABLE([gparity],
      [ac_GPARITY=${enable_gparity}], [ac_GPARITY=yes])
 
 AM_CONDITIONAL(BUILD_GPARITY, [ test "${ac_GPARITY}X" == "yesX" ])
-
-AC_ARG_ENABLE([zmobius],
-     [AC_HELP_STRING([--enable-zmobius=yes|no], [enable Zmobius support])],
-     [ac_ZMOBIUS=${enable_zmobius}], [ac_ZMOBIUS=yes])
-
-AM_CONDITIONAL(BUILD_ZMOBIUS, [ test "${ac_ZMOBIUS}X" == "yesX" ])
-
-
 case ${ac_FERMION_REPS} in
    yes) AC_DEFINE([ENABLE_FERMION_REPS],[1],[non QCD fermion reps]);;
 esac
 case ${ac_GPARITY} in
    yes) AC_DEFINE([ENABLE_GPARITY],[1],[fermion actions with GPARITY BCs]);;
 esac
-case ${ac_ZMOBIUS} in
-   yes) AC_DEFINE([ENABLE_ZMOBIUS],[1],[Zmobius fermion actions]);;
-esac
 ############### Nc
 AC_ARG_ENABLE([Nc],
     [AC_HELP_STRING([--enable-Nc=2|3|4], [enable number of colours])],
@@ -444,7 +428,7 @@ case ${ax_cv_cxx_compiler_vendor} in
         SIMD_FLAGS='-mavx2 -mfma -mf16c';;
       AVX512)
         AC_DEFINE([AVX512],[1],[AVX512 intrinsics])
-        SIMD_FLAGS='-mavx512f -mavx512cd';;
+        SIMD_FLAGS='-mavx512f -mavx512pf -mavx512er -mavx512cd';;
       SKL)
         AC_DEFINE([AVX512],[1],[AVX512 intrinsics for SkyLake Xeon])
         SIMD_FLAGS='-march=skylake-avx512';;
@@ -497,9 +481,6 @@ case ${ax_cv_cxx_compiler_vendor} in
         AC_DEFINE([AVX2],[1],[AVX2 intrinsics])
         SIMD_FLAGS='-march=core-avx2 -xcore-avx2';;
       AVX512)
-        AC_DEFINE([AVX512],[1],[AVX512 intrinsics])
-        SIMD_FLAGS='-xcommon-avx512';;
-      SKL)
         AC_DEFINE([AVX512],[1],[AVX512 intrinsics])
         SIMD_FLAGS='-xcore-avx512';;
       KNC)

tests/core/Test_where.cc

@@ -40,9 +40,9 @@ int main (int argc, char ** argv)
 
   int N=16;
   
-  std::vector<int> latt_size  ({N,N,N,N});
+  std::vector<int> latt_size  ({N,4,4});
-  std::vector<int> simd_layout({vComplexD::Nsimd(),1,1,1});
+  std::vector<int> simd_layout({vComplexD::Nsimd(),1,1});
-  std::vector<int> mpi_layout ({1,1,1,1});
+  std::vector<int> mpi_layout ({1,1,1});
 
   int vol = 1;
   int nd  = latt_size.size();
@@ -69,7 +69,7 @@ int main (int argc, char ** argv)
     for(int t=0;t<latt_size[mu];t++){
       LatticeCoordinate(coor,mu);
       sl=where(coor==Integer(t),rn,zz);
-      //      std::cout <<GridLogMessage<< " sl " << sl<<std::endl;
+      std::cout <<GridLogMessage<< " sl " << sl<<std::endl;
       std::cout <<GridLogMessage<<" slice "<<t<<" " << norm2(sl)<<std::endl;
       ns=ns+norm2(sl);
     }

tests/core/Test_where_extended.cc

@@ -1,143 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_poisson_fft.cc
-
-    Copyright (C) 2015
-
-Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
-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 */
-#include <Grid/Grid.h>
-
-using namespace Grid;
- ;
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  int threads = GridThread::GetThreads();
-  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
-
-  int N=16;
-  
-  std::vector<int> latt_size  ({N,4,4});
-  std::vector<int> simd_layout({vComplexD::Nsimd(),1,1});
-  std::vector<int> mpi_layout ({1,1,1});
-
-  int vol = 1;
-  int nd  = latt_size.size();
-  for(int d=0;d<nd;d++){
-    vol = vol * latt_size[d];
-  }
-
-  GridCartesian         GRID(latt_size,simd_layout,mpi_layout);
-  GridParallelRNG RNG(&GRID);
-  RNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));  
-
-  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
-  std::cout<<GridLogMessage<<"== LatticeComplex =="<<std::endl;
-  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
-  {
-    LatticeComplexD      zz(&GRID);
-    LatticeInteger     coor(&GRID);
-    LatticeComplexD  rn(&GRID);
-    LatticeComplexD  sl(&GRID);
-
-    zz  = ComplexD(0.0,0.0);
-
-    gaussian(RNG,rn);
-    
-    RealD nn=norm2(rn);
-    for(int mu=0;mu<nd;mu++){
-      RealD ns=0.0;
-      for(int t=0;t<latt_size[mu];t++){
-	LatticeCoordinate(coor,mu);
-	sl=where(coor==Integer(t),rn,zz);
-	std::cout <<GridLogMessage<<" slice "<<t<<" " << norm2(sl)<<std::endl;
-	ns=ns+norm2(sl);
-      }
-      std::cout <<GridLogMessage <<" sliceNorm" <<mu<<" "<< nn <<" "<<ns<<" err " << nn-ns<<std::endl;
-      assert(abs(nn-ns) < 1.0e-10);
-    }
-  }
-
-
-  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
-  std::cout<<GridLogMessage<<"== LatticeFermion =="<<std::endl;
-  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
-  {
-    LatticeFermionD      zz(&GRID);
-    LatticeInteger     coor(&GRID);
-    LatticeFermionD  rn(&GRID);
-    LatticeFermionD  sl(&GRID);
-
-    zz  = ComplexD(0.0,0.0);
-
-    gaussian(RNG,rn);
-    
-    RealD nn=norm2(rn);
-    for(int mu=0;mu<nd;mu++){
-      RealD ns=0.0;
-      for(int t=0;t<latt_size[mu];t++){
-	LatticeCoordinate(coor,mu);
-	sl=where(coor==Integer(t),rn,zz);
-	std::cout <<GridLogMessage<<" slice "<<t<<" " << norm2(sl)<<std::endl;
-	ns=ns+norm2(sl);
-      }
-      std::cout <<GridLogMessage <<" sliceNorm" <<mu<<" "<< nn <<" "<<ns<<" err " << nn-ns<<std::endl;
-      assert(abs(nn-ns) < 1.0e-10);
-    }
-  }
-
-
-
-  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
-  std::cout<<GridLogMessage<<"== LatticePropagator =="<<std::endl;
-  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
-
-  {
-    LatticePropagatorD      zz(&GRID);
-    LatticeInteger     coor(&GRID);
-    LatticePropagatorD  rn(&GRID);
-    LatticePropagatorD  sl(&GRID);
-
-    zz  = ComplexD(0.0,0.0);
-
-    gaussian(RNG,rn);
-    
-    RealD nn=norm2(rn);
-    for(int mu=0;mu<nd;mu++){
-      RealD ns=0.0;
-      for(int t=0;t<latt_size[mu];t++){
-	LatticeCoordinate(coor,mu);
-	sl=where(coor==Integer(t),rn,zz);
-	std::cout <<GridLogMessage<<" slice "<<t<<" " << norm2(sl)<<std::endl;
-	ns=ns+norm2(sl);
-      }
-      std::cout <<GridLogMessage <<" sliceNorm" <<mu<<" "<< nn <<" "<<ns<<" err " << nn-ns<<std::endl;
-      assert(abs(nn-ns) < 1.0e-10);
-    }
-  }
-
-  Grid_finalize();
-}

tests/debug/Test_cayley_mres.cc

@@ -33,14 +33,13 @@ using namespace Grid;
 
 
 template<class What> 
-void  TestConserved(What & Ddwf,
+void  TestConserved(What & Ddwf, What & Ddwfrev, 
 		    LatticeGaugeField &Umu,
 		    GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
 		    GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
 		    RealD mass, RealD M5,
 		    GridParallelRNG *RNG4,
-		    GridParallelRNG *RNG5,
+		    GridParallelRNG *RNG5);
-                    What *Ddwfrev=nullptr);
 
   Gamma::Algebra Gmu [] = {
     Gamma::Algebra::GammaX,
@@ -103,11 +102,10 @@ int main (int argc, char ** argv)
   GridRedBlackCartesian * FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGridF);
 
 
+  std::vector<int> seeds4({1,2,3,4});
   std::vector<int> seeds5({5,6,7,8});
   GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
-  GridParallelRNG          RNG4(UGrid);
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
-  std::vector<int> seeds4({1,2,3,4}); RNG4.SeedFixedIntegers(seeds4);
-  //const std::string seeds4{ "test-gauge-3000" }; RNG4.SeedUniqueString( seeds4 );
 
   LatticeGaugeField Umu(UGrid);
   if( argc > 1 && argv[1][0] != '-' )
@@ -118,8 +116,8 @@ int main (int argc, char ** argv)
   }
   else
   {
-    std::cout<<GridLogMessage <<"Using hot configuration"<<std::endl;
+    std::cout<<GridLogMessage <<"Using cold configuration"<<std::endl;
-    // SU<Nc>::ColdConfiguration(Umu);
+    //SU<Nc>::ColdConfiguration(Umu);
     SU<Nc>::HotConfiguration(RNG4,Umu);
   }
 
@@ -129,7 +127,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage <<"DomainWallFermion test"<<std::endl;
   std::cout<<GridLogMessage <<"======================"<<std::endl;
   DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  TestConserved<DomainWallFermionR>(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestConserved<DomainWallFermionR>(Ddwf,Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   RealD b=1.5;// Scale factor b+c=2, b-c=1
   RealD c=0.5;
@@ -139,13 +137,13 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage <<"MobiusFermion test"<<std::endl;
   std::cout<<GridLogMessage <<"======================"<<std::endl;
   MobiusFermionR Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
-  TestConserved<MobiusFermionR>(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestConserved<MobiusFermionR>(Dmob,Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"ScaledShamirFermion test"<<std::endl;
   std::cout<<GridLogMessage <<"======================"<<std::endl;
   ScaledShamirFermionR Dsham(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,2.0);
-  TestConserved<ScaledShamirFermionR>(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestConserved<ScaledShamirFermionR>(Dsham,Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"ZMobiusFermion test"<<std::endl;
@@ -154,7 +152,8 @@ int main (int argc, char ** argv)
   //  for(int s=0;s<Ls;s++) omegasrev[s]=omegas[Ls-1-s];
   ZMobiusFermionR ZDmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,omegas,b,c);
   ZMobiusFermionR ZDmobrev(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,omegasrev,b,c);
-  TestConserved<ZMobiusFermionR>(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5,&ZDmobrev);
+
+  TestConserved<ZMobiusFermionR>(ZDmob,ZDmobrev,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   Grid_finalize();
 }
@@ -162,17 +161,22 @@ int main (int argc, char ** argv)
 
 
 template<class Action> 
-void  TestConserved(Action & Ddwf,
+void  TestConserved(Action & Ddwf, 
+		    Action & Ddwfrev, 
 		    LatticeGaugeField &Umu,
 		    GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
 		    GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
 		    RealD mass, RealD M5,
 		    GridParallelRNG *RNG4,
-		    GridParallelRNG *RNG5,
+		    GridParallelRNG *RNG5)
-                    Action * Ddwfrev)
 {
-  LatticePropagator phys_src(UGrid);
+  int Ls=Ddwf.Ls;
-  LatticePropagator seqsrc(FGrid);
+
+  LatticePropagator phys_src(UGrid); 
+
+  std::vector<LatticeColourMatrix> U(4,UGrid);
+  
+  LatticePropagator seqsrc(FGrid); 
   LatticePropagator prop5(FGrid); 
   LatticePropagator prop5rev(FGrid); 
   LatticePropagator prop4(UGrid); 
@@ -190,9 +194,9 @@ void  TestConserved(Action & Ddwf,
   phys_src=Zero();
   pokeSite(kronecker,phys_src,coor);
   
+  MdagMLinearOperator<Action,LatticeFermion> HermOp(Ddwf);
+  MdagMLinearOperator<Action,LatticeFermion> HermOprev(Ddwfrev);
   ConjugateGradient<LatticeFermion> CG(1.0e-16,100000);
-  SchurRedBlackDiagTwoSolve<LatticeFermion> schur(CG);
-  ZeroGuesser<LatticeFermion> zpg;
   for(int s=0;s<Nd;s++){
     for(int c=0;c<Nc;c++){
       LatticeFermion src4  (UGrid); 
@@ -202,20 +206,20 @@ void  TestConserved(Action & Ddwf,
       Ddwf.ImportPhysicalFermionSource(src4,src5);
 
       LatticeFermion result5(FGrid); result5=Zero();
-      schur(Ddwf,src5,result5,zpg);
+
-      std::cout<<GridLogMessage<<"spin "<<s<<" color "<<c<<" norm2(sourc5d) "<<norm2(src5)
+      // CGNE
-               <<" norm2(result5d) "<<norm2(result5)<<std::endl;
+      LatticeFermion Mdagsrc5  (FGrid); 
+      Ddwf.Mdag(src5,Mdagsrc5);
+      CG(HermOp,Mdagsrc5,result5);
       FermToProp<Action>(prop5,result5,s,c);
 
       LatticeFermion result4(UGrid);
       Ddwf.ExportPhysicalFermionSolution(result5,result4);
       FermToProp<Action>(prop4,result4,s,c);
 
-      if( Ddwfrev ) {
+      Ddwfrev.ImportPhysicalFermionSource(src4,src5);
-        Ddwfrev->ImportPhysicalFermionSource(src4,src5);
+      Ddwfrev.Mdag(src5,Mdagsrc5);
-        result5 = Zero();
+      CG(HermOprev,Mdagsrc5,result5);
-        schur(*Ddwfrev,src5,result5,zpg);
-      }
       FermToProp<Action>(prop5rev,result5,s,c);
     }
   }
@@ -247,7 +251,11 @@ void  TestConserved(Action & Ddwf,
       PropToFerm<Action>(src5,seqsrc,s,c);
 
       LatticeFermion result5(FGrid); result5=Zero();
-      schur(Ddwf,src5,result5,zpg);
+
+      // CGNE
+      LatticeFermion Mdagsrc5  (FGrid); 
+      Ddwf.Mdag(src5,Mdagsrc5);
+      CG(HermOp,Mdagsrc5,result5);
 
       LatticeFermion result4(UGrid);
       Ddwf.ExportPhysicalFermionSolution(result5,result4);
@@ -268,10 +276,10 @@ void  TestConserved(Action & Ddwf,
   Ddwf.ContractConservedCurrent(prop5rev,prop5,Vector_mu,phys_src,Current::Vector,Tdir);
   Ddwf.ContractJ5q(prop5,PJ5q);
   
-  PA       = trace(g5*Axial_mu);      // Pseudoscalar-Axial conserved current
+  PA       = trace(g5*Axial_mu);
-  SV       = trace(Vector_mu);        // Scalar-Vector conserved current
+  SV       = trace(Vector_mu);
-  VV       = trace(gT*Vector_mu);     // (local) Vector-Vector conserved current
+  VV       = trace(gT*Vector_mu);
-  PP       = trace(adj(prop4)*prop4); // Pseudoscalar density
+  PP       = trace(adj(prop4)*prop4);
   
   // Spatial sum
   sliceSum(PA,sumPA,Tdir);
@@ -280,17 +288,15 @@ void  TestConserved(Action & Ddwf,
   sliceSum(PP,sumPP,Tdir);
   sliceSum(PJ5q,sumPJ5q,Tdir);
 
-  const int Nt{static_cast<int>(sumPA.size())};
+  int Nt=sumPA.size();
-  std::cout<<GridLogMessage<<"Vector Ward identity by timeslice (~ 0)"<<std::endl;
   for(int t=0;t<Nt;t++){
-    std::cout<<GridLogMessage <<" t "<<t<<" SV "<<real(TensorRemove(sumSV[t]))<<" VV "<<real(TensorRemove(sumVV[t]))<<std::endl;
+    std::cout <<" SV "<<real(TensorRemove(sumSV[t]));
+    std::cout <<" VV "<<real(TensorRemove(sumVV[t]))<<std::endl;
   }
-  std::cout<<GridLogMessage<<"Axial Ward identity by timeslice (defect ~ 0)"<<std::endl;
   for(int t=0;t<Nt;t++){
-    const RealD DmuPAmu{real(TensorRemove(sumPA[t]-sumPA[(t-1+Nt)%Nt]))};
+    std::cout <<" PAc "<<real(TensorRemove(sumPA[t]));
-    std::cout<<GridLogMessage<<" t "<<t<<" DmuPAmu "<<DmuPAmu
+    std::cout <<" PJ5q "<<real(TensorRemove(sumPJ5q[t]));
-             <<" PP "<<real(TensorRemove(sumPP[t]))<<" PJ5q "<<real(TensorRemove(sumPJ5q[t]))
+    std::cout <<" Ward Identity defect " <<real(TensorRemove(sumPA[t]-sumPA[(t-1+Nt)%Nt] - 2.0*(Ddwf.mass*sumPP[t] + sumPJ5q[t]) ))<<"\n";
-             <<" Ward Identity defect " <<(DmuPAmu - 2.*real(TensorRemove(Ddwf.mass*sumPP[t] + sumPJ5q[t])))<<std::endl;
   }
   
   ///////////////////////////////

tests/debug/Test_heatbath_dwf_eofa.cc

@@ -66,9 +66,7 @@ int main(int argc, char** argv)
   // Set up RNGs
   std::vector<int> seeds4({1, 2, 3, 4});
   std::vector<int> seeds5({5, 6, 7, 8});
-  GridSerialRNG sRNG;
   GridParallelRNG RNG5(FGrid);
-  sRNG.SeedFixedIntegers(seeds5);
   RNG5.SeedFixedIntegers(seeds5);
   GridParallelRNG RNG4(UGrid);
   RNG4.SeedFixedIntegers(seeds4);
@@ -86,7 +84,7 @@ int main(int argc, char** argv)
     ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
     ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, false);
 
-    Meofa.refresh(Umu,sRNG, RNG5);
+    Meofa.refresh(Umu, RNG5);
     printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
   }
 
@@ -96,7 +94,7 @@ int main(int argc, char** argv)
     ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
     ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, true);
 
-    Meofa.refresh(Umu,sRNG, RNG5);
+    Meofa.refresh(Umu, RNG5);
     printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
   }