Switching Block lanczos precision to explicitly single

Adding sample run script and input file

.gitignore

@@ -10,6 +10,8 @@
 *~
 *#
 *.sublime-*
+.ctags
+tags
 
 # Precompiled Headers #
 #######################

.travis.yml

@@ -1,56 +0,0 @@
-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/algorithms/approx/Chebyshev.h

@@ -240,12 +240,14 @@ public:
     Field T0(grid); T0 = in;  
     Field T1(grid); 
     Field T2(grid);
+    Field Tout(grid);
     Field y(grid);
       
     Field *Tnm = &T0;
     Field *Tn  = &T1;
     Field *Tnp = &T2;
 
+    std::cout << GridLogMessage << "Chebyshev() starts"<<std::endl;
     // Tn=T1 = (xscale M + mscale)in
     RealD xscale = 2.0/(hi-lo);
     RealD mscale = -(hi+lo)/(hi-lo);
@@ -254,7 +256,7 @@ public:
 
     // sum = .5 c[0] T0 + c[1] T1
     //    out = ()*T0 + Coeffs[1]*T1;
-    axpby(out,0.5*Coeffs[0],Coeffs[1],T0,T1);
+    axpby(Tout,0.5*Coeffs[0],Coeffs[1],T0,T1);
     for(int n=2;n<order;n++){
 
       Linop.HermOp(*Tn,y);
@@ -275,7 +277,7 @@ public:
       axpby(y,xscale,mscale,y,(*Tn));
       axpby(*Tnp,2.0,-1.0,y,(*Tnm));
       if ( Coeffs[n] != 0.0) {
-	axpy(out,Coeffs[n],*Tnp,out);
+	axpy(Tout,Coeffs[n],*Tnp,Tout);
       }
 #endif
       // Cycle pointers to avoid copies
@@ -285,6 +287,8 @@ public:
       Tnp    =swizzle;
 	  
     }
+    out = Tout;
+    std::cout << GridLogMessage << "Chebyshev() ends"<<std::endl;
   }
 };
 
@@ -377,24 +381,26 @@ public:
     Field T0(grid); T0 = in;  
     Field T1(grid); 
     Field T2(grid);
+    Field Tout(grid);
     Field  y(grid);
       
     Field *Tnm = &T0;
     Field *Tn  = &T1;
     Field *Tnp = &T2;
 
+    std::cout << GridLogMessage << "ChebyshevLanczos() starts"<<std::endl;
     // Tn=T1 = (xscale M )*in
     AminusMuSq(Linop,T0,T1);
 
     // sum = .5 c[0] T0 + c[1] T1
-    out = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1;
+    Tout = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1;
     for(int n=2;n<order;n++){
 	
       AminusMuSq(Linop,*Tn,y);
 
       *Tnp=2.0*y-(*Tnm);
 
-      out=out+Coeffs[n]* (*Tnp);
+      Tout=Tout+Coeffs[n]* (*Tnp);
 
       // Cycle pointers to avoid copies
       Field *swizzle = Tnm;
@@ -403,6 +409,8 @@ public:
       Tnp    =swizzle;
 	  
     }
+    out=Tout;
+    std::cout << GridLogMessage << "ChebyshevLanczos() ends"<<std::endl;
   }
 };
 NAMESPACE_END(Grid);

Grid/allocator/MemoryManager.cc

@@ -44,7 +44,7 @@ void  MemoryManager::AcceleratorFree    (void *ptr,size_t bytes)
   if ( __freeme ) {
     acceleratorFreeDevice(__freeme);
     total_device-=bytes;
-    //    PrintBytes();
+//       PrintBytes();
   }
 }
 void *MemoryManager::SharedAllocate(size_t bytes)
@@ -53,8 +53,8 @@ void *MemoryManager::SharedAllocate(size_t bytes)
   if ( ptr == (void *) NULL ) {
     ptr = (void *) acceleratorAllocShared(bytes);
     total_shared+=bytes;
-    //    std::cout <<"AcceleratorAllocate: allocated Shared pointer "<<std::hex<<ptr<<std::dec<<std::endl;
-    //    PrintBytes();
+        std::cout <<GridLogMessage<<"AcceleratorAllocate: allocated Shared pointer "<<std::hex<<ptr<<std::dec<<std::endl;
+//        PrintBytes();
   }
   return ptr;
 }
@@ -74,6 +74,7 @@ void *MemoryManager::CpuAllocate(size_t bytes)
   if ( ptr == (void *) NULL ) {
     ptr = (void *) acceleratorAllocShared(bytes);
     total_host+=bytes;
+//    std::cout << GridLogMessage<< "MemoryManager:: CpuAllocate  total_host= "<<total_host<<" "<< ptr << std::endl;
   }
   return ptr;
 }
@@ -83,6 +84,7 @@ void  MemoryManager::CpuFree    (void *_ptr,size_t bytes)
   void *__freeme = Insert(_ptr,bytes,Cpu);
   if ( __freeme ) { 
     acceleratorFreeShared(__freeme);
+//    std::cout << GridLogMessage<< "MemoryManager:: CpuFree  total_host= "<<total_host<<" "<< __freeme << std::endl;
     total_host-=bytes;
   }
 }

Grid/communicator/SharedMemoryMPI.cc

@@ -35,6 +35,9 @@ Author: Christoph Lehner <christoph@lhnr.de>
 #endif
 #ifdef GRID_HIP
 #include <hip/hip_runtime_api.h>
+#endif
+#ifdef GRID_SYCl
+
 #endif
 
 NAMESPACE_BEGIN(Grid); 
@@ -70,6 +73,7 @@ void GlobalSharedMemory::Init(Grid_MPI_Comm comm)
   WorldNodes = WorldSize/WorldShmSize;
   assert( (WorldNodes * WorldShmSize) == WorldSize );
 
+
   // FIXME: Check all WorldShmSize are the same ?
 
   /////////////////////////////////////////////////////////////////////
@@ -446,7 +450,47 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 ////////////////////////////////////////////////////////////////////////////////////////////
 // Hugetlbfs mapping intended
 ////////////////////////////////////////////////////////////////////////////////////////////
-#if defined(GRID_CUDA) ||defined(GRID_HIP)
+#if defined(GRID_CUDA) ||defined(GRID_HIP)  || defined(GRID_SYCL)
+
+//if defined(GRID_SYCL)
+#if 0
+void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
+{
+  void * ShmCommBuf ; 
+  assert(_ShmSetup==1);
+  assert(_ShmAlloc==0);
+
+  //////////////////////////////////////////////////////////////////////////////////////////////////////////
+  // allocate the pointer array for shared windows for our group
+  //////////////////////////////////////////////////////////////////////////////////////////////////////////
+  MPI_Barrier(WorldShmComm);
+  WorldShmCommBufs.resize(WorldShmSize);
+
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////////
+  // Each MPI rank should allocate our own buffer
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////////
+  ShmCommBuf = acceleratorAllocDevice(bytes);
+
+  if (ShmCommBuf == (void *)NULL ) {
+    std::cerr << " SharedMemoryMPI.cc acceleratorAllocDevice failed NULL pointer for " << bytes<<" bytes " << std::endl;
+    exit(EXIT_FAILURE);  
+  }
+
+  std::cout << WorldRank << header " SharedMemoryMPI.cc acceleratorAllocDevice "<< bytes 
+	    << "bytes at "<< std::hex<< ShmCommBuf <<std::dec<<" for comms buffers " <<std::endl;
+
+  SharedMemoryZero(ShmCommBuf,bytes);
+
+  assert(WorldShmSize == 1);
+  for(int r=0;r<WorldShmSize;r++){
+    WorldShmCommBufs[r] = ShmCommBuf;
+  }
+  _ShmAllocBytes=bytes;
+  _ShmAlloc=1;
+}
+#endif
+
+#if defined(GRID_CUDA) ||defined(GRID_HIP) ||defined(GRID_SYCL)  
 void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 {
   void * ShmCommBuf ; 
@@ -469,8 +513,16 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
   ///////////////////////////////////////////////////////////////////////////////////////////////////////////
   // Each MPI rank should allocate our own buffer
   ///////////////////////////////////////////////////////////////////////////////////////////////////////////
+#ifdef GRID_SYCL_LEVEL_ZERO_IPC
+  auto zeDevice = cl::sycl::get_native<cl::sycl::backend::level_zero>(theGridAccelerator->get_device());
+  auto zeContext= cl::sycl::get_native<cl::sycl::backend::level_zero>(theGridAccelerator->get_context());
+  ze_device_mem_alloc_desc_t zeDesc = {};
+  zeMemAllocDevice(zeContext,&zeDesc,bytes,2*1024*1024,zeDevice,&ShmCommBuf);
+  std::cout << WorldRank << header " SharedMemoryMPI.cc zeMemAllocDevice "<< bytes 
+	      << "bytes at "<< std::hex<< ShmCommBuf <<std::dec<<" for comms buffers " <<std::endl;
+#else  
   ShmCommBuf = acceleratorAllocDevice(bytes);
-
+#endif  
   if (ShmCommBuf == (void *)NULL ) {
     std::cerr << " SharedMemoryMPI.cc acceleratorAllocDevice failed NULL pointer for " << bytes<<" bytes " << std::endl;
     exit(EXIT_FAILURE);  
@@ -480,8 +532,8 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
     std::cout << WorldRank << header " SharedMemoryMPI.cc acceleratorAllocDevice "<< bytes 
 	      << "bytes at "<< std::hex<< ShmCommBuf <<std::dec<<" for comms buffers " <<std::endl;
   }
-  SharedMemoryZero(ShmCommBuf,bytes);
-
+  //  SharedMemoryZero(ShmCommBuf,bytes);
+  std::cout<< "Setting up IPC"<<std::endl;
   ///////////////////////////////////////////////////////////////////////////////////////////////////////////
   // Loop over ranks/gpu's on our node
   ///////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -491,6 +543,23 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
     //////////////////////////////////////////////////
     // If it is me, pass around the IPC access key
     //////////////////////////////////////////////////
+#ifdef GRID_SYCL_LEVEL_ZERO_IPC
+    ze_ipc_mem_handle_t handle;
+    if ( r==WorldShmRank ) { 
+      auto err = zeMemGetIpcHandle(zeContext,ShmCommBuf,&handle);
+      if ( err != ZE_RESULT_SUCCESS ) {
+	std::cerr << "SharedMemoryMPI.cc zeMemGetIpcHandle failed for rank "<<r<<" "<<std::hex<<err<<std::dec<<std::endl;
+	exit(EXIT_FAILURE);
+      } else {
+	std::cerr << "SharedMemoryMPI.cc zeMemGetIpcHandle succeeded for rank "<<r<<" "<<std::hex<<err<<std::dec<<std::endl;
+      }
+      std::cerr<<"Allocated IpcHandle rank "<<r<<" (hex) ";
+      for(int c=0;c<ZE_MAX_IPC_HANDLE_SIZE;c++){
+	std::cerr<<std::hex<<(uint32_t)((uint8_t)handle.data[c])<<std::dec;
+      }
+      std::cerr<<std::endl;
+    }
+#endif
 #ifdef GRID_CUDA
     cudaIpcMemHandle_t handle;
     if ( r==WorldShmRank ) { 
@@ -527,6 +596,25 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
     // If I am not the source, overwrite thisBuf with remote buffer
     ///////////////////////////////////////////////////////////////
     void * thisBuf = ShmCommBuf;
+#ifdef GRID_SYCL_LEVEL_ZERO_IPC
+    if ( r!=WorldShmRank ) {
+      thisBuf = nullptr;
+      std::cerr<<"Using IpcHandle rank "<<r<<" ";
+      for(int c=0;c<ZE_MAX_IPC_HANDLE_SIZE;c++){
+	std::cerr<<std::hex<<(uint32_t)((uint8_t)handle.data[c])<<std::dec;
+      }
+      std::cerr<<std::endl;
+      auto err = zeMemOpenIpcHandle(zeContext,zeDevice,handle,0,&thisBuf);
+      if ( err != ZE_RESULT_SUCCESS ) {
+	std::cerr << "SharedMemoryMPI.cc "<<zeContext<<" "<<zeDevice<<std::endl;
+	std::cerr << "SharedMemoryMPI.cc zeMemOpenIpcHandle failed for rank "<<r<<" "<<std::hex<<err<<std::dec<<std::endl; 
+	exit(EXIT_FAILURE);
+      } else {
+	std::cerr << "SharedMemoryMPI.cc zeMemOpenIpcHandle succeeded for rank "<<r<<" "<<std::hex<<err<<std::dec<<std::endl;
+      }
+      assert(thisBuf!=nullptr);
+    }
+#endif
 #ifdef GRID_CUDA
     if ( r!=WorldShmRank ) { 
       auto err = cudaIpcOpenMemHandle(&thisBuf,handle,cudaIpcMemLazyEnablePeerAccess);
@@ -557,6 +645,8 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
   _ShmAllocBytes=bytes;
   _ShmAlloc=1;
 }
+#endif
+
 #else 
 #ifdef GRID_MPI3_SHMMMAP
 void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
@@ -727,16 +817,16 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 /////////////////////////////////////////////////////////////////////////
 void GlobalSharedMemory::SharedMemoryZero(void *dest,size_t bytes)
 {
-#ifdef GRID_CUDA
-  cudaMemset(dest,0,bytes);
+#if defined(GRID_CUDA) || defined(GRID_HIP) || defined(GRID_SYCL)
+  acceleratorMemSet(dest,0,bytes);
 #else
   bzero(dest,bytes);
 #endif
 }
 void GlobalSharedMemory::SharedMemoryCopy(void *dest,void *src,size_t bytes)
 {
-#ifdef GRID_CUDA
-  cudaMemcpy(dest,src,bytes,cudaMemcpyDefault);
+#if defined(GRID_CUDA) || defined(GRID_HIP) || defined(GRID_SYCL)
+  acceleratorCopyToDevice(src,dest,bytes);
 #else   
   bcopy(src,dest,bytes);
 #endif
@@ -800,7 +890,7 @@ void SharedMemory::SetCommunicator(Grid_MPI_Comm comm)
   }
 #endif
 
-  SharedMemoryTest();
+  //SharedMemoryTest();
 }
 //////////////////////////////////////////////////////////////////
 // On node barrier

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];
-  cshiftVector<vobj> send_buf(buffer_size);
-  cshiftVector<vobj> recv_buf(buffer_size);
+  static cshiftVector<vobj> send_buf; send_buf.resize(buffer_size);
+  static cshiftVector<vobj> recv_buf; recv_buf.resize(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);
 
-  std::vector<cshiftVector<scalar_object> >  send_buf_extract(Nsimd);
-  std::vector<cshiftVector<scalar_object> >  recv_buf_extract(Nsimd);
+  static std::vector<cshiftVector<scalar_object> >  send_buf_extract; send_buf_extract.resize(Nsimd);
+  static std::vector<cshiftVector<scalar_object> >  recv_buf_extract; recv_buf_extract.resize(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];
-  cshiftVector<vobj> send_buf_v(buffer_size);
-  cshiftVector<vobj> recv_buf_v(buffer_size);
+  static cshiftVector<vobj> send_buf_v; send_buf_v.resize(buffer_size);
+  static cshiftVector<vobj> recv_buf_v; recv_buf_v.resize(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);
 
-  std::vector<cshiftVector<scalar_object> >  send_buf_extract(Nsimd);
-  std::vector<cshiftVector<scalar_object> >  recv_buf_extract(Nsimd);
+  static std::vector<cshiftVector<scalar_object> >  send_buf_extract; send_buf_extract.resize(Nsimd);
+  static std::vector<cshiftVector<scalar_object> >  recv_buf_extract; recv_buf_extract.resize(Nsimd);
   scalar_object *  recv_buf_extract_mpi;
   scalar_object *  send_buf_extract_mpi;
   {

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

@@ -198,18 +198,27 @@ public:
       std::cerr << " nersc_csum  " <<std::hex<< nersc_csum << " " << header.checksum<< std::dec<< std::endl;
       exit(0);
     }
-    assert(fabs(clone.plaquette -header.plaquette ) < 1.0e-5 );
+    assert(fabs(clone.plaquette -header.plaquette ) < 1.0e-1 );
     assert(fabs(clone.link_trace-header.link_trace) < 1.0e-6 );
     assert(nersc_csum == header.checksum );
       
     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;
@@ -219,8 +228,8 @@ public:
     // Following should become arguments
     ///////////////////////////////////////////
     header.sequence_number = 1;
-    header.ensemble_id     = "UKQCD";
-    header.ensemble_label  = "DWF";
+    header.ensemble_id     = std::string("UKQCD");
+    header.ensemble_label  = ens_label;
 
     typedef LorentzColourMatrixD fobj3D;
     typedef LorentzColour2x3D    fobj2D;
@@ -232,7 +241,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/fermion/Fermion.h

@@ -291,12 +291,6 @@ 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

@@ -183,7 +183,8 @@ NAMESPACE_CHECK(ImplStaggered);
 /////////////////////////////////////////////////////////////////////////////
 // Single flavour one component spinors with colour index. 5d vec
 /////////////////////////////////////////////////////////////////////////////
-#include <Grid/qcd/action/fermion/StaggeredVec5dImpl.h> 
-NAMESPACE_CHECK(ImplStaggered5dVec);  
+// Deprecate Vec5d
+//#include <Grid/qcd/action/fermion/StaggeredVec5dImpl.h> 
+//NAMESPACE_CHECK(ImplStaggered5dVec);  
 
 

Grid/qcd/action/fermion/StaggeredImpl.h

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

Grid/qcd/action/fermion/WilsonImpl.h

@@ -184,18 +184,22 @@ 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);
@@ -208,6 +212,29 @@ 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/implementation/CayleyFermion5DImplementation.h

@@ -880,11 +880,23 @@ 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++){
 
@@ -907,7 +919,7 @@ void CayleyFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
 
     tmp    = Cshift(tmp,mu,1);
     Impl::multLinkField(Utmp,this->Umu,tmp,mu);
-    tmp    = G_s[s]*( Utmp*ph - gmu*Utmp*ph ); // Forward hop
+    tmp    = sign*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

@@ -680,7 +680,8 @@ 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,
@@ -790,7 +791,7 @@ template <> void StaggeredKernels<StaggeredVec5dImplD>::DhopSiteAsm(StencilView
 #endif
 }
    
-   
+#endif   
 
 
 #define PERMUTE_DIR3 __asm__ (	\

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

@@ -32,25 +32,50 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 NAMESPACE_BEGIN(Grid);
 
-#define LOAD_CHI(b)		\
+#ifdef GRID_SIMT
+
+#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=ref()()(0);\
-    Chi_1=ref()()(1);\
-    Chi_2=ref()()(2);
+  Chi_0=coalescedRead(ref()()(0),lane);	\
+  Chi_1=coalescedRead(ref()()(1),lane);	\
+  Chi_2=coalescedRead(ref()()(2),lane);
+
+#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));			\
-   Impl::loadLinkElement(U_00,ref()(0,0));      \
-   Impl::loadLinkElement(U_10,ref()(1,0));      \
-   Impl::loadLinkElement(U_20,ref()(2,0));      \
-   Impl::loadLinkElement(U_01,ref()(0,1));      \
-   Impl::loadLinkElement(U_11,ref()(1,1));      \
-   Impl::loadLinkElement(U_21,ref()(2,1));      \
-   Impl::loadLinkElement(U_02,ref()(0,2));     \
-   Impl::loadLinkElement(U_12,ref()(1,2));     \
-   Impl::loadLinkElement(U_22,ref()(2,2));     \
+    U_00=coalescedRead(ref()(0,0),lane);				\
+    U_10=coalescedRead(ref()(1,0),lane);				\
+    U_20=coalescedRead(ref()(2,0),lane);				\
+    U_01=coalescedRead(ref()(0,1),lane);				\
+    U_11=coalescedRead(ref()(1,1),lane);				\
+    U_21=coalescedRead(ref()(2,1),lane);				\
+    U_02=coalescedRead(ref()(0,2),lane);				\
+    U_12=coalescedRead(ref()(1,2),lane);				\
+    U_22=coalescedRead(ref()(2,2),lane);				\
     UChi ## _0  = U_00*Chi_0;	       \
     UChi ## _1  = U_10*Chi_0;\
     UChi ## _2  = U_20*Chi_0;\
@@ -63,15 +88,15 @@ NAMESPACE_BEGIN(Grid);
 
 #define MULT_ADD(U,A,UChi)			\
   auto & ref(U[sU](A));			\
-   Impl::loadLinkElement(U_00,ref()(0,0));      \
-   Impl::loadLinkElement(U_10,ref()(1,0));      \
-   Impl::loadLinkElement(U_20,ref()(2,0));      \
-   Impl::loadLinkElement(U_01,ref()(0,1));      \
-   Impl::loadLinkElement(U_11,ref()(1,1));      \
-   Impl::loadLinkElement(U_21,ref()(2,1));      \
-   Impl::loadLinkElement(U_02,ref()(0,2));     \
-   Impl::loadLinkElement(U_12,ref()(1,2));     \
-   Impl::loadLinkElement(U_22,ref()(2,2));     \
+    U_00=coalescedRead(ref()(0,0),lane);				\
+    U_10=coalescedRead(ref()(1,0),lane);				\
+    U_20=coalescedRead(ref()(2,0),lane);				\
+    U_01=coalescedRead(ref()(0,1),lane);				\
+    U_11=coalescedRead(ref()(1,1),lane);				\
+    U_21=coalescedRead(ref()(2,1),lane);				\
+    U_02=coalescedRead(ref()(0,2),lane);				\
+    U_12=coalescedRead(ref()(1,2),lane);				\
+    U_22=coalescedRead(ref()(2,2),lane);				\
     UChi ## _0 += U_00*Chi_0;	       \
     UChi ## _1 += U_10*Chi_0;\
     UChi ## _2 += U_20*Chi_0;\
@@ -83,24 +108,18 @@ 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(in);					\
+    LOAD_CHI(Perm,in);						\
     if ( perm) {						\
       PERMUTE_DIR(Perm);					\
     }								\
   } else {							\
-    LOAD_CHI(buf);						\
+    LOAD_CHI_COMMS(buf);					\
   }								
 
 #define HAND_STENCIL_LEG_BEGIN(Dir,Perm,skew,even)		\
@@ -116,19 +135,18 @@ 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(in);				\
+    LOAD_CHI(Perm,in);					\
     if ( perm) {					\
       PERMUTE_DIR(Perm);				\
     }							\
   } else if ( st.same_node[Dir] ) {			\
-    LOAD_CHI(buf);					\
+    LOAD_CHI_COMMS(buf);				\
   }							\
   if (local || st.same_node[Dir] ) {		\
     MULT_ADD(U,Dir,even);				\
@@ -140,10 +158,32 @@ NAMESPACE_BEGIN(Grid);
   local  = SE->_is_local;				\
   if ((!local) && (!st.same_node[Dir]) ) {		\
     nmu++;							\
-    { LOAD_CHI(buf);	  }					\
+    { LOAD_CHI_COMMS(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
@@ -155,28 +195,14 @@ void StaggeredKernels<Impl>::DhopSiteHand(StencilView &st,
   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;
 
-  Simd Chi_0;    // two spinor; 6 regs
-  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; 
+  const int Nsimd = SiteHalfSpinor::Nsimd();
+  const int lane=acceleratorSIMTlane(Nsimd);
+  typedef decltype( coalescedRead( in[0]()()(0) )) Simt;
+  HAND_DECLARATIONS(Simt);
 
-  SiteSpinor result;
+  typedef decltype( coalescedRead( in[0] )) calcSiteSpinor;
+  calcSiteSpinor result;
   int offset,local,perm, ptype;
 
   StencilEntry *SE;
@@ -215,7 +241,7 @@ void StaggeredKernels<Impl>::DhopSiteHand(StencilView &st,
       result()()(1) = even_1 + odd_1;
       result()()(2) = even_2 + odd_2;
     }
-    vstream(out[sF],result);
+    coalescedWrite(out[sF],result);
   }
 }
 
@@ -230,28 +256,13 @@ void StaggeredKernels<Impl>::DhopSiteHandInt(StencilView &st,
   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();
+  const int lane=acceleratorSIMTlane(Nsimd);
+  typedef decltype( coalescedRead( in[0]()()(0) )) Simt;
+  HAND_DECLARATIONS(Simt);
 
-  Simd Chi_0;    // two spinor; 6 regs
-  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;
+  typedef decltype( coalescedRead( in[0] )) calcSiteSpinor;
+  calcSiteSpinor result;
   int offset, ptype, local, perm;
 
   StencilEntry *SE;
@@ -261,8 +272,8 @@ void StaggeredKernels<Impl>::DhopSiteHandInt(StencilView &st,
   //    int sF=s+LLs*sU;
   {
 
-    even_0 = Zero();    even_1 = Zero();    even_2 = Zero();
-     odd_0 = Zero();     odd_1 = Zero();     odd_2 = Zero();
+    zeroit(even_0);    zeroit(even_1);    zeroit(even_2);
+    zeroit(odd_0);    zeroit(odd_1);    zeroit(odd_2);
 
     skew = 0;
     HAND_STENCIL_LEG_INT(U,Xp,3,skew,even);  
@@ -294,7 +305,7 @@ void StaggeredKernels<Impl>::DhopSiteHandInt(StencilView &st,
       result()()(1) = even_1 + odd_1;
       result()()(2) = even_2 + odd_2;
     }
-    vstream(out[sF],result);
+    coalescedWrite(out[sF],result);
   }
 }
 
@@ -309,28 +320,13 @@ void StaggeredKernels<Impl>::DhopSiteHandExt(StencilView &st,
   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();
+  const int lane=acceleratorSIMTlane(Nsimd);
+  typedef decltype( coalescedRead( in[0]()()(0) )) Simt;
+  HAND_DECLARATIONS(Simt);
 
-  Simd Chi_0;    // two spinor; 6 regs
-  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;
+  typedef decltype( coalescedRead( in[0] )) calcSiteSpinor;
+  calcSiteSpinor result;
   int offset, ptype, local;
 
   StencilEntry *SE;
@@ -340,8 +336,8 @@ void StaggeredKernels<Impl>::DhopSiteHandExt(StencilView &st,
   //    int sF=s+LLs*sU;
   {
 
-    even_0 = Zero();    even_1 = Zero();    even_2 = Zero();
-     odd_0 = Zero();     odd_1 = Zero();     odd_2 = Zero();
+    zeroit(even_0);    zeroit(even_1);    zeroit(even_2);
+    zeroit(odd_0);    zeroit(odd_1);    zeroit(odd_2);
     int nmu=0;
     skew = 0;
     HAND_STENCIL_LEG_EXT(U,Xp,3,skew,even);  
@@ -374,7 +370,7 @@ void StaggeredKernels<Impl>::DhopSiteHandExt(StencilView &st,
 	result()()(1) = even_1 + odd_1;
 	result()()(2) = even_2 + odd_2;
       }
-      out[sF] = out[sF] + result;
+      coalescedWrite(out[sF] , out(sF)+ result);
     }
   }
 }
@@ -397,6 +393,7 @@ 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,39 +35,32 @@ NAMESPACE_BEGIN(Grid);
 #define GENERIC_STENCIL_LEG(U,Dir,skew,multLink)		\
   SE = st.GetEntry(ptype, Dir+skew, sF);			\
   if (SE->_is_local ) {						\
-    if (SE->_permute) {						\
-      chi_p = χ						\
-      permute(chi,  in[SE->_offset], ptype);			\
-    } else {							\
-      chi_p = &in[SE->_offset];					\
-    }								\
+    int perm= SE->_permute;						\
+    chi = coalescedReadPermute(in[SE->_offset],ptype,perm,lane);\
   } else {							\
-    chi_p = &buf[SE->_offset];					\
+    chi = coalescedRead(buf[SE->_offset],lane);			\
   }								\
-  multLink(Uchi, U[sU], *chi_p, Dir);			
+  acceleratorSynchronise();					\
+  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 ) {						\
-    if (SE->_permute) {						\
-      chi_p = χ						\
-      permute(chi,  in[SE->_offset], ptype);			\
-    } else {							\
-      chi_p = &in[SE->_offset];					\
-    }								\
+    int perm= SE->_permute;						\
+    chi = coalescedReadPermute(in[SE->_offset],ptype,perm,lane);\
   } else if ( st.same_node[Dir] ) {				\
-    chi_p = &buf[SE->_offset];					\
+    chi = coalescedRead(buf[SE->_offset],lane);                 \
   }								\
   if (SE->_is_local || st.same_node[Dir] ) {			\
-    multLink(Uchi, U[sU], *chi_p, Dir);				\
+    multLink(Uchi, U[sU], chi, 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_p = &buf[SE->_offset];					\
-    multLink(Uchi, U[sU], *chi_p, Dir);				\
+    chi = coalescedRead(buf[SE->_offset],lane);			\
+    multLink(Uchi, U[sU], chi, Dir);				\
   }
 
 template <class Impl>
@@ -84,12 +77,14 @@ void StaggeredKernels<Impl>::DhopSiteGeneric(StencilView &st,
 					     SiteSpinor *buf, int sF, int sU, 
 					     const FermionFieldView &in, FermionFieldView &out, int dag) 
 {
-  const SiteSpinor *chi_p;
-  SiteSpinor chi;
-  SiteSpinor Uchi;
+  typedef decltype(coalescedRead(in[0])) calcSpinor;
+  calcSpinor chi;
+  calcSpinor Uchi;
   StencilEntry *SE;
   int ptype;
   int skew;
+  const int Nsimd = SiteHalfSpinor::Nsimd();
+  const int lane=acceleratorSIMTlane(Nsimd);
 
   //  for(int s=0;s<LLs;s++){
   //
@@ -118,7 +113,7 @@ void StaggeredKernels<Impl>::DhopSiteGeneric(StencilView &st,
     if ( dag ) { 
       Uchi = - Uchi;
     } 
-    vstream(out[sF], Uchi);
+    coalescedWrite(out[sF], Uchi,lane);
   }
 };
 
@@ -130,13 +125,16 @@ template <int Naik> accelerator_inline
 void StaggeredKernels<Impl>::DhopSiteGenericInt(StencilView &st, 
 						DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU,
 						SiteSpinor *buf, int sF, int sU, 
-						const FermionFieldView &in, FermionFieldView &out,int dag) {
-  const SiteSpinor *chi_p;
-  SiteSpinor chi;
-  SiteSpinor Uchi;
+						const FermionFieldView &in, FermionFieldView &out,int dag)
+{
+  typedef decltype(coalescedRead(in[0])) calcSpinor;
+  calcSpinor chi;
+  calcSpinor 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;
@@ -165,7 +163,7 @@ void StaggeredKernels<Impl>::DhopSiteGenericInt(StencilView &st,
     if ( dag ) {
       Uchi = - Uchi;
     }
-    vstream(out[sF], Uchi);
+    coalescedWrite(out[sF], Uchi,lane);
   }
 };
 
@@ -178,14 +176,17 @@ template <int Naik> accelerator_inline
 void StaggeredKernels<Impl>::DhopSiteGenericExt(StencilView &st, 
 						DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU,
 						SiteSpinor *buf, int sF, int sU,
-						const FermionFieldView &in, FermionFieldView &out,int dag) {
-  const SiteSpinor *chi_p;
-  //  SiteSpinor chi;
-  SiteSpinor Uchi;
+						const FermionFieldView &in, FermionFieldView &out,int dag)
+{
+  typedef decltype(coalescedRead(in[0])) calcSpinor;
+  calcSpinor chi;
+  calcSpinor 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;
@@ -211,11 +212,12 @@ 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 ( dag ) { 
-	out[sF] = out[sF] - Uchi;
+    if ( nmu ) {
+      auto _out = coalescedRead(out[sF],lane);
+      if ( dag ) {
+	coalescedWrite(out[sF], _out-Uchi,lane);
       } else { 
-	out[sF] = out[sF] + Uchi;
+	coalescedWrite(out[sF], _out+Uchi,lane);
       }
     }
   }
@@ -261,6 +263,8 @@ void StaggeredKernels<Impl>::DhopImproved(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 , UUU, AcceleratorRead);
   autoView( U_v   ,   U, AcceleratorRead);
   autoView( in_v  ,  in, AcceleratorRead);
@@ -301,6 +305,8 @@ 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/smearing/StoutSmearing.h

@@ -85,21 +85,18 @@ public:
 
     std::cout << GridLogDebug << "Stout smearing started\n";
 
-    // Smear the configurations
+    // C contains the staples multiplied by some rho
+    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 )
-        tmp = 1.0;  // Don't smear in the orthogonal direction
-      else {
-        tmp = peekLorentz(C, mu);
-        Umu = peekLorentz(U, mu);
-        iq_mu = Ta(
-                   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
+      if( mu == OrthogDim ) continue ;
+      // u_smr = exp(iQ_mu)*U_mu apart from Orthogdim
+      Umu = peekLorentz(U, mu);
+      tmp = peekLorentz(C, mu);
+      iq_mu = Ta( tmp * adj(Umu));  
+      exponentiate_iQ(tmp, iq_mu);
+      pokeLorentz(u_smr, tmp * Umu, mu);
     }
     std::cout << GridLogDebug << "Stout smearing completed\n";
   };

Grid/tensors/Tensor_SIMT.h

@@ -65,7 +65,8 @@ void coalescedWriteNonTemporal(vobj & __restrict__ vec,const vobj & __restrict__
 #else
 
 
-#ifndef GRID_SYCL
+//#ifndef GRID_SYCL
+#if 1
 // 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

Grid/tensors/Tensor_extract_merge.h

@@ -1,5 +1,5 @@
 /*************************************************************************************
-
+n
     Grid physics library, www.github.com/paboyle/Grid 
 
     Source file: ./lib/tensors/Tensor_extract_merge.h
@@ -153,7 +153,7 @@ void insertLane(int lane, vobj & __restrict__ vec,const typename vobj::scalar_ob
 // Extract to a bunch of scalar object pointers of different scalar type, with offset. Useful for precision change
 ////////////////////////////////////////////////////////////////////////
 template<class vobj, class sobj> accelerator
-void extract(const vobj &vec,ExtractPointerArray<sobj> &extracted, int offset)
+void extract(const vobj &vec,const ExtractPointerArray<sobj> &extracted, int offset)
 {
   typedef typename GridTypeMapper<sobj>::scalar_type sobj_scalar_type;
   typedef typename GridTypeMapper<vobj>::scalar_type scalar_type;
@@ -181,7 +181,7 @@ void extract(const vobj &vec,ExtractPointerArray<sobj> &extracted, int offset)
 // Merge bunch of scalar object pointers of different scalar type, with offset. Useful for precision change
 ////////////////////////////////////////////////////////////////////////
 template<class vobj, class sobj> accelerator
-void merge(vobj &vec,ExtractPointerArray<sobj> &extracted, int offset)
+void merge(vobj &vec,const ExtractPointerArray<sobj> &extracted, int offset)
 {
   typedef typename GridTypeMapper<sobj>::scalar_type sobj_scalar_type;
   typedef typename GridTypeMapper<vobj>::scalar_type scalar_type;

Grid/tensors/Tensor_outer.h

@@ -34,6 +34,16 @@ 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>
@@ -57,17 +67,6 @@ 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

@@ -171,7 +171,6 @@ void acceleratorInit(void)
 #ifdef GRID_SYCL
 
 cl::sycl::queue *theGridAccelerator;
-
 void acceleratorInit(void)
 {
   int nDevices = 1;
@@ -179,6 +178,10 @@ void acceleratorInit(void)
   cl::sycl::device selectedDevice { selector };
   theGridAccelerator = new sycl::queue (selectedDevice);
 
+#ifdef GRID_SYCL_LEVEL_ZERO_IPC
+  zeInit(0);
+#endif
+  
   char * localRankStr = NULL;
   int rank = 0, world_rank=0; 
 #define ENV_LOCAL_RANK_OMPI    "OMPI_COMM_WORLD_LOCAL_RANK"

Grid/threads/Accelerator.h

@@ -39,6 +39,10 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #ifdef HAVE_MM_MALLOC_H
 #include <mm_malloc.h>
 #endif
+#ifdef __APPLE__
+// no memalign
+inline void *memalign(size_t align, size_t bytes) { return malloc(bytes); }
+#endif
 
 NAMESPACE_BEGIN(Grid);
 
@@ -194,6 +198,9 @@ inline void *acceleratorAllocShared(size_t bytes)
     ptr = (void *) NULL;
     printf(" cudaMallocManaged failed for %d %s \n",bytes,cudaGetErrorString(err));
   }
+  size_t free,total;
+  cudaMemGetInfo(&free,&total);
+  std::cout<<"cudaMemGetInfo  "<<free<<" / "<<total<<std::endl;
   return ptr;
 };
 inline void *acceleratorAllocDevice(size_t bytes)
@@ -204,6 +211,9 @@ inline void *acceleratorAllocDevice(size_t bytes)
     ptr = (void *) NULL;
     printf(" cudaMalloc failed for %d %s \n",bytes,cudaGetErrorString(err));
   }
+  size_t free,total;
+  cudaMemGetInfo(&free,&total);
+  std::cout<<"cudaMemGetInfo  "<<free<<" / "<<total<<std::endl;
   return ptr;
 };
 inline void acceleratorFreeShared(void *ptr){ cudaFree(ptr);};
@@ -233,6 +243,13 @@ inline int  acceleratorIsCommunicable(void *ptr)
 NAMESPACE_END(Grid);
 #include <CL/sycl.hpp>
 #include <CL/sycl/usm.hpp>
+
+#define GRID_SYCL_LEVEL_ZERO_IPC
+
+#ifdef GRID_SYCL_LEVEL_ZERO_IPC
+#include <level_zero/ze_api.h>
+#include <CL/sycl/backend/level_zero.hpp>
+#endif
 NAMESPACE_BEGIN(Grid);
 
 extern cl::sycl::queue *theGridAccelerator;
@@ -257,11 +274,14 @@ accelerator_inline int acceleratorSIMTlane(int Nsimd) {
       unsigned long nt=acceleratorThreads();				\
       unsigned long unum1 = num1;					\
       unsigned long unum2 = num2;					\
+      if(nt < 8)nt=8;							\
       cl::sycl::range<3> local {nt,1,nsimd};				\
       cl::sycl::range<3> global{unum1,unum2,nsimd};			\
       cgh.parallel_for<class dslash>(					\
       cl::sycl::nd_range<3>(global,local), \
-      [=] (cl::sycl::nd_item<3> item) /*mutable*/ {   \
+      [=] (cl::sycl::nd_item<3> item) /*mutable*/     \
+      [[intel::reqd_sub_group_size(8)]]	      \
+      {						      \
       auto iter1    = item.get_global_id(0);	      \
       auto iter2    = item.get_global_id(1);	      \
       auto lane     = item.get_global_id(2);	      \
@@ -409,6 +429,8 @@ inline void acceleratorMemSet(void *base,int value,size_t bytes) { hipMemset(bas
 
 #undef GRID_SIMT
 
+
+
 #define accelerator 
 #define accelerator_inline strong_inline
 #define accelerator_for(iterator,num,nsimd, ... )   thread_for(iterator, num, { __VA_ARGS__ });
@@ -457,7 +479,7 @@ accelerator_inline void acceleratorSynchronise(void)
   __syncwarp();
 #endif
 #ifdef GRID_SYCL
-  cl::sycl::detail::workGroupBarrier();
+  //cl::sycl::detail::workGroupBarrier();
 #endif
 #ifdef GRID_HIP
   __syncthreads();

Grid/util/Init.cc

@@ -56,6 +56,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 static int
 feenableexcept (unsigned int excepts)
 {
+#if 0
+  // Fails on Apple M1
   static fenv_t fenv;
   unsigned int new_excepts = excepts & FE_ALL_EXCEPT;
   unsigned int old_excepts;  // previous masks
@@ -70,6 +72,8 @@ feenableexcept (unsigned int excepts)
 
   iold_excepts  = (int) old_excepts;
   return ( fesetenv (&fenv) ? -1 : iold_excepts );
+#endif
+  return 0;
 }
 #endif
 
@@ -163,6 +167,13 @@ void GridCmdOptionInt(std::string &str,int & val)
   return;
 }
 
+// ypj [add]
+void GridCmdOptionFloat(std::string &str,double & val)
+{
+  std::stringstream ss(str);
+  ss>>val;
+  return;
+}
 
 void GridParseLayout(char **argv,int argc,
 		     Coordinate &latt_c,

Grid/util/Init.h

@@ -57,7 +57,8 @@ void GridCmdOptionCSL(std::string str,std::vector<std::string> & vec);
 template<class VectorInt>
 void GridCmdOptionIntVector(const std::string &str,VectorInt & vec);
 void GridCmdOptionInt(std::string &str,int & val);
-
+// ypj [add]
+void GridCmdOptionFloat(std::string &str,double & val);
 
 void GridParseLayout(char **argv,int argc,
 		     std::vector<int> &latt,

Makefile.am

@@ -1,5 +1,5 @@
 # additional include paths necessary to compile the C++ library
-SUBDIRS = Grid HMC benchmarks tests
+SUBDIRS = Grid HMC benchmarks tests examples
 
 include $(top_srcdir)/doxygen.inc
 

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) [![Travis status](https://travis-ci.org/paboyle/Grid.svg?branch=develop)](https://travis-ci.org/paboyle/Grid)
+# 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) 
 
 **Data parallel C++ mathematical object library.**
 

benchmarks/Benchmark_ITT.cc

@@ -133,34 +133,30 @@ public:
 
 	std::vector<HalfSpinColourVectorD *> xbuf(8);
 	std::vector<HalfSpinColourVectorD *> rbuf(8);
-	Grid.ShmBufferFreeAll();
+	//Grid.ShmBufferFreeAll();
+	uint64_t bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD);
 	for(int d=0;d<8;d++){
-	  xbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
-	  rbuf[d] = (HalfSpinColourVectorD *)Grid.ShmBufferMalloc(lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
+	  xbuf[d] = (HalfSpinColourVectorD *)acceleratorAllocDevice(bytes);
+	  rbuf[d] = (HalfSpinColourVectorD *)acceleratorAllocDevice(bytes);
 	  //	  bzero((void *)xbuf[d],lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
 	  //	  bzero((void *)rbuf[d],lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
 	}
 
-	int bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD);
 	int ncomm;
 	double dbytes;
-	std::vector<double> times(Nloop);
-	for(int i=0;i<Nloop;i++){
 
-	  double start=usecond();
+        for(int dir=0;dir<8;dir++) {
+	  int mu =dir % 4;
+	  if (mpi_layout[mu]>1 ) {
 
-	  dbytes=0;
-	  ncomm=0;
+	    std::vector<double> times(Nloop);
+	    for(int i=0;i<Nloop;i++){
 
-	  thread_for(dir,8,{
-		     
-	    double tbytes;
-	    int mu =dir % 4;
-
-	    if (mpi_layout[mu]>1 ) {
-	        
+	      dbytes=0;	        
+	      double start=usecond();
 	      int xmit_to_rank;
 	      int recv_from_rank;
+
 	      if ( dir == mu ) { 
 		int comm_proc=1;
 		Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
@@ -168,40 +164,38 @@ public:
 		int comm_proc = mpi_layout[mu]-1;
 		Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
 	      }
-	      tbytes= Grid.StencilSendToRecvFrom((void *)&xbuf[dir][0], xmit_to_rank,
-						 (void *)&rbuf[dir][0], recv_from_rank,
-						 bytes,dir);
-	      thread_critical {
-		ncomm++;
-		dbytes+=tbytes;
-	      }
+	      Grid.SendToRecvFrom((void *)&xbuf[dir][0], xmit_to_rank,
+				  (void *)&rbuf[dir][0], recv_from_rank,
+				  bytes);
+	      dbytes+=bytes;
+	     
+	      double stop=usecond();
+	      t_time[i] = stop-start; // microseconds
+
 	    }
-          });
-	  Grid.Barrier();
-	  double stop=usecond();
-	  t_time[i] = stop-start; // microseconds
+	    timestat.statistics(t_time);
+	  
+	    dbytes=dbytes*ppn;
+	    double xbytes    = dbytes*0.5;
+	    double bidibytes = dbytes;
+	  
+	    std::cout<<GridLogMessage << lat<<"\t"<<Ls<<"\t "
+		     << bytes << " \t "
+		     <<xbytes/timestat.mean<<" \t "<< xbytes*timestat.err/(timestat.mean*timestat.mean)<< " \t "
+		     <<xbytes/timestat.max <<" "<< xbytes/timestat.min  
+		     << "\t\t"<< bidibytes/timestat.mean<< "  " << bidibytes*timestat.err/(timestat.mean*timestat.mean) << " "
+		     << bidibytes/timestat.max << " " << bidibytes/timestat.min << std::endl;
+	  }
 	}
-
-	timestat.statistics(t_time);
-
-	dbytes=dbytes*ppn;
-	double xbytes    = dbytes*0.5;
-	//	double rbytes    = dbytes*0.5;
-	double bidibytes = dbytes;
-
-	std::cout<<GridLogMessage << lat<<"\t"<<Ls<<"\t "
-		 << bytes << " \t "
-		 <<xbytes/timestat.mean<<" \t "<< xbytes*timestat.err/(timestat.mean*timestat.mean)<< " \t "
-		 <<xbytes/timestat.max <<" "<< xbytes/timestat.min  
-		 << "\t\t"<< bidibytes/timestat.mean<< "  " << bidibytes*timestat.err/(timestat.mean*timestat.mean) << " "
-		 << bidibytes/timestat.max << " " << bidibytes/timestat.min << std::endl;
-	
-	    }
-    }    
-
+	for(int d=0;d<8;d++){
+	  acceleratorFreeDevice(xbuf[d]);
+	  acceleratorFreeDevice(rbuf[d]);
+	}
+      }
+    } 
     return;
   }
-
+  
 
   
   static void Memory(void)
@@ -281,7 +275,6 @@ public:
 
 
     uint64_t lmax=32;
-#define NLOOP (1000*lmax*lmax*lmax*lmax/lat/lat/lat/lat)
 
     GridSerialRNG          sRNG;      sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
     for(int lat=8;lat<=lmax;lat+=8){
@@ -445,7 +438,7 @@ public:
 	// 1344= 3*(2*8+6)*2*8 + 8*3*2*2 + 3*4*2*8
 	// 1344 = Nc* (6+(Nc-1)*8)*2*Nd + Nd*Nc*2*2  + Nd*Nc*Ns*2
 	//	double flops=(1344.0*volume)/2;
-#if 1
+#if 0
 	double fps = Nc* (6+(Nc-1)*8)*Ns*Nd + Nd*Nc*Ns  + Nd*Nc*Ns*2;
 #else
 	double fps = Nc* (6+(Nc-1)*8)*Ns*Nd + 2*Nd*Nc*Ns  + 2*Nd*Nc*Ns*2;
@@ -716,12 +709,12 @@ int main (int argc, char ** argv)
 
   if ( do_su4 ) {
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << " Memory benchmark " <<std::endl;
+    std::cout<<GridLogMessage << " SU(4) benchmark " <<std::endl;
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
     Benchmark::SU4();
   }
   
-  if ( do_comms && (NN>1) ) {
+  if ( do_comms ) {
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
     std::cout<<GridLogMessage << " Communications benchmark " <<std::endl;
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;

configure.ac

@@ -815,6 +815,7 @@ AC_CONFIG_FILES(tests/smearing/Makefile)
 AC_CONFIG_FILES(tests/qdpxx/Makefile)
 AC_CONFIG_FILES(tests/testu01/Makefile)
 AC_CONFIG_FILES(benchmarks/Makefile)
+AC_CONFIG_FILES(examples/Makefile)
 AC_OUTPUT
 
 echo ""

examples/Example_Laplacian.cc

@@ -0,0 +1,396 @@
+#include <Grid/Grid.h>
+using namespace Grid;
+
+/*
+/////////////////////////////////////////////////////////////////////////////////////////////
+// Grid/algorithms/SparseMatrix.h: Interface defining what I expect of a general sparse matrix, such as a Fermion action
+/////////////////////////////////////////////////////////////////////////////////////////////
+template<class Field> class SparseMatrixBase {
+public:
+  virtual GridBase *Grid(void) =0;
+
+  virtual void  M    (const Field &in, Field &out)=0;
+  virtual void  Mdag (const Field &in, Field &out)=0;
+  virtual void  MdagM(const Field &in, Field &out) {
+    Field tmp (in.Grid());
+    M(in,tmp);
+    Mdag(tmp,out);
+  }
+  virtual  void Mdiag    (const Field &in, Field &out)=0;
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp)=0;
+  virtual  void MdirAll  (const Field &in, std::vector<Field> &out)=0;
+};
+*/
+
+const std::vector<int> directions   ({Xdir,Ydir,Zdir,Xdir,Ydir,Zdir});
+const std::vector<int> displacements({1,1,1,-1,-1,-1});
+
+template<class Field> class FreeLaplacianCshift : public SparseMatrixBase<Field>
+{
+public:
+  GridBase *grid;
+  FreeLaplacianCshift(GridBase *_grid)
+  {
+    grid=_grid;
+  };
+  virtual GridBase *Grid(void) { return grid; };
+
+  virtual void  M    (const Field &in, Field &out)
+  {
+    out = Zero();
+    for(int mu=0;mu<Nd-1;mu++) {
+      out = out + Cshift(in,mu,1) + Cshift(in,mu,-1) - 2.0 * in;
+    }
+  };
+  virtual void  Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
+  virtual  void Mdiag    (const Field &in, Field &out)                  {assert(0);}; // Unimplemented need only for multigrid
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
+  virtual  void MdirAll  (const Field &in, std::vector<Field> &out)     {assert(0);}; // Unimplemented need only for multigrid
+};
+
+template<class Gimpl,class Field> class CovariantLaplacianCshift : public SparseMatrixBase<Field>
+{
+public:
+  INHERIT_GIMPL_TYPES(Gimpl);
+
+  GridBase *grid;
+  GaugeField U;
+  
+  CovariantLaplacianCshift(GaugeField &_U)    :
+    grid(_U.Grid()),
+    U(_U) {  };
+
+  virtual GridBase *Grid(void) { return grid; };
+
+  virtual void  M    (const Field &in, Field &out)
+  {
+    out=Zero();
+    for(int mu=0;mu<Nd-1;mu++) {
+      GaugeLinkField Umu = PeekIndex<LorentzIndex>(U, mu); // NB: Inefficent
+      out = out + Gimpl::CovShiftForward(Umu,mu,in);    
+      out = out + Gimpl::CovShiftBackward(Umu,mu,in);    
+      out = out - 2.0*in;
+    }
+  };
+  virtual void  Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
+  virtual  void Mdiag    (const Field &in, Field &out)                  {assert(0);}; // Unimplemented need only for multigrid
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
+  virtual  void MdirAll  (const Field &in, std::vector<Field> &out)     {assert(0);}; // Unimplemented need only for multigrid
+};
+
+
+#define LEG_LOAD(Dir)						 \
+  SE = st.GetEntry(ptype, Dir, ss);				 \
+  if (SE->_is_local ) {						 \
+    int perm= SE->_permute;					 \
+    chi = coalescedReadPermute(in[SE->_offset],ptype,perm,lane); \
+  } else {							 \
+    chi = coalescedRead(buf[SE->_offset],lane);			 \
+  }								 \
+  acceleratorSynchronise();
+
+template<class Field> class FreeLaplacianStencil : public SparseMatrixBase<Field>
+{
+public:
+  typedef typename Field::vector_object siteObject;
+  typedef CartesianStencil<siteObject, siteObject, int> StencilImpl;
+
+  GridBase *grid;
+  StencilImpl Stencil;
+  SimpleCompressor<siteObject> Compressor;
+  
+  FreeLaplacianStencil(GridBase *_grid)
+    : Stencil    (_grid,6,Even,directions,displacements,0), grid(_grid)
+  {  };
+  
+  virtual GridBase *Grid(void) { return grid; };
+
+  virtual void  M    (const Field &_in, Field &_out)
+  {
+
+    ///////////////////////////////////////////////
+    // Halo exchange for this geometry of stencil
+    ///////////////////////////////////////////////
+    Stencil.HaloExchange(_in, Compressor);
+
+    ///////////////////////////////////
+    // Arithmetic expressions
+    ///////////////////////////////////
+
+    // Views; device friendly/accessible pointers
+    auto st = Stencil.View(AcceleratorRead);
+    auto buf = st.CommBuf();
+    autoView( in     , _in    , AcceleratorRead);
+    autoView( out    , _out   , AcceleratorWrite);
+
+    typedef typename Field::vector_object        vobj;
+    typedef decltype(coalescedRead(in[0])) calcObj;
+
+    const int      Nsimd = vobj::Nsimd();
+    const uint64_t NN = grid->oSites();
+
+    accelerator_for( ss, NN, Nsimd, {
+
+	StencilEntry *SE;
+	
+	const int lane=acceleratorSIMTlane(Nsimd);
+
+	calcObj chi;
+	calcObj res;
+	int ptype;
+
+	res                 = coalescedRead(in[ss])*(-6.0);
+	LEG_LOAD(0);	res = res + chi;
+	LEG_LOAD(1);	res = res + chi;
+	LEG_LOAD(2);	res = res + chi;
+	LEG_LOAD(3);	res = res + chi;
+	LEG_LOAD(4);	res = res + chi;
+	LEG_LOAD(5);	res = res + chi;
+
+	coalescedWrite(out[ss], res,lane);
+	
+    });
+  };
+  virtual void  Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
+  virtual  void Mdiag    (const Field &in, Field &out)                  {assert(0);}; // Unimplemented need only for multigrid
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
+  virtual  void MdirAll  (const Field &in, std::vector<Field> &out)     {assert(0);}; // Unimplemented need only for multigrid
+};
+
+template<class Gimpl,class Field> class CovariantLaplacianStencil : public SparseMatrixBase<Field>
+{
+public:
+  INHERIT_GIMPL_TYPES(Gimpl);
+
+  typedef typename Field::vector_object siteObject;
+
+  template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Nc> >, Nds>;
+  typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
+  typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
+
+  typedef CartesianStencil<siteObject, siteObject, int> StencilImpl;
+
+  GridBase *grid;
+  StencilImpl Stencil;
+  SimpleCompressor<siteObject> Compressor;
+  DoubledGaugeField Uds;
+  CovariantLaplacianStencil(GaugeField &Umu)
+    :
+      grid(Umu.Grid()),
+      Stencil    (grid,6,Even,directions,displacements,0),
+      Uds(grid)
+  {
+    for (int mu = 0; mu < Nd; mu++) {
+      auto U = PeekIndex<LorentzIndex>(Umu, mu);
+      PokeIndex<LorentzIndex>(Uds, U, mu );
+      U = adj(Cshift(U, mu, -1));
+      PokeIndex<LorentzIndex>(Uds, U, mu + 4);
+    }
+  };
+  
+  virtual GridBase *Grid(void) { return grid; };
+
+  virtual void  M    (const Field &_in, Field &_out)
+  {
+    ///////////////////////////////////////////////
+    // Halo exchange for this geometry of stencil
+    ///////////////////////////////////////////////
+    Stencil.HaloExchange(_in, Compressor);
+
+    ///////////////////////////////////
+    // Arithmetic expressions
+    ///////////////////////////////////
+    auto st = Stencil.View(AcceleratorRead);
+    auto buf = st.CommBuf();
+
+    autoView( in     , _in    , AcceleratorRead);
+    autoView( out    , _out   , AcceleratorWrite);
+    autoView( U     , Uds    , AcceleratorRead);
+
+    typedef typename Field::vector_object        vobj;
+    typedef decltype(coalescedRead(in[0]))    calcObj;
+    typedef decltype(coalescedRead(U[0](0))) calcLink;
+
+    const int      Nsimd = vobj::Nsimd();
+    const uint64_t NN = grid->oSites();
+
+    accelerator_for( ss, NN, Nsimd, {
+
+	StencilEntry *SE;
+	
+	const int lane=acceleratorSIMTlane(Nsimd);
+
+	calcObj chi;
+	calcObj res;
+	calcObj Uchi;
+	calcLink UU;
+	int ptype;
+
+	res                 = coalescedRead(in[ss])*(-6.0);
+
+#define LEG_LOAD_MULT(leg,polarisation)			\
+	UU = coalescedRead(U[ss](polarisation));	\
+	LEG_LOAD(leg);					\
+	mult(&Uchi(), &UU, &chi());			\
+	res = res + Uchi;
+	
+	LEG_LOAD_MULT(0,Xp);
+	LEG_LOAD_MULT(1,Yp);
+	LEG_LOAD_MULT(2,Zp);
+	LEG_LOAD_MULT(3,Xm);
+	LEG_LOAD_MULT(4,Ym);
+	LEG_LOAD_MULT(5,Zm);
+
+	coalescedWrite(out[ss], res,lane);
+    });
+  };
+  virtual void  Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
+  virtual  void Mdiag    (const Field &in, Field &out)                  {assert(0);}; // Unimplemented need only for multigrid
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
+  virtual  void MdirAll  (const Field &in, std::vector<Field> &out)     {assert(0);}; // Unimplemented need only for multigrid
+};
+
+#undef LEG_LOAD_MULT
+#undef LEG_LOAD
+	
+int main(int argc, char ** argv)
+{
+  Grid_init(&argc, &argv);
+
+  typedef LatticeColourVector Field;
+
+  auto latt_size   = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  auto mpi_layout  = GridDefaultMpi();
+
+  GridCartesian    Grid(latt_size,simd_layout,mpi_layout);
+  GridParallelRNG  RNG(&Grid);  RNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+
+  FreeLaplacianCshift<Field>  FLcs(&Grid);
+  FreeLaplacianStencil<Field> FLst(&Grid);
+
+  LatticeGaugeField U(&Grid);
+
+  SU<Nc>::ColdConfiguration(RNG,U);
+
+  std::cout << " Gauge field has norm " <<norm2(U)<<std::endl;
+
+  CovariantLaplacianCshift <PeriodicGimplR,Field> CLcs(U);
+  CovariantLaplacianStencil<PeriodicGimplR,Field> CLst(U);
+
+  Field in(&Grid); gaussian(RNG,in);
+  Field out_FLcs(&Grid);
+  Field out_FLst(&Grid);
+  Field out_CLcs(&Grid);
+  Field out_CLst(&Grid);
+  Field diff(&Grid);
+
+  ////////////////////////////////////////////////////////
+  // First test: in free field these should all agree
+  ////////////////////////////////////////////////////////
+  FLcs.M(in,out_FLcs);
+  FLst.M(in,out_FLst);
+  CLcs.M(in,out_CLcs);
+  CLst.M(in,out_CLst);
+
+  std:: cout << "******************************************************************" <<std::endl;
+  std:: cout << " Test A: consistency of four different Laplacian implementations " <<std::endl;
+  std:: cout << "******************************************************************" <<std::endl;
+  std:: cout << " Input test vector " <<norm2(in)<<std::endl;
+  std:: cout << "--------------------------------------------------------" <<std::endl;
+  std:: cout << " Free cshift  output vector " <<norm2(out_FLcs)<<std::endl;
+  std:: cout << " Free stencil output vector " <<norm2(out_FLst)<<std::endl;
+  std:: cout << " Cov  cshift  output vector " <<norm2(out_CLcs)<<std::endl;
+  std:: cout << " Cov  stencil output vector " <<norm2(out_CLst)<<std::endl;
+  std:: cout << "--------------------------------------------------------" <<std::endl;
+
+  diff = out_FLcs - out_FLst;
+  std:: cout << " Difference between free Cshift Laplacian and free Stencil Laplacian      = " <<norm2(diff)<<std::endl;
+
+  diff = out_FLcs - out_CLcs;
+  std:: cout << " Difference between free Cshift Laplacian and covariant Cshift Laplacian  = " <<norm2(diff)<<std::endl;
+
+  diff = out_FLcs - out_CLst;
+  std:: cout << " Difference between free Cshift Laplacian and covariant Stencil Laplacian = " <<norm2(diff)<<std::endl;
+  std:: cout << "--------------------------------------------------------" <<std::endl;
+  
+
+  std:: cout << "******************************************************************" <<std::endl;
+  std:: cout << " Test B: gauge covariance  " <<std::endl;
+  std:: cout << "******************************************************************" <<std::endl;
+
+  LatticeGaugeField     U_GT(&Grid); // Gauge transformed field
+  LatticeColourMatrix   g(&Grid);   // local Gauge xform matrix
+
+  U_GT = U;
+  // Make a random xform to teh gauge field
+  SU<Nc>::RandomGaugeTransform(RNG,U_GT,g); // Unit gauge
+
+  Field in_GT(&Grid); 
+  Field out_GT(&Grid);
+
+  Field out_CLcs_GT(&Grid);
+  Field out_CLst_GT(&Grid);
+
+  CovariantLaplacianCshift <PeriodicGimplR,Field> CLcs_GT(U_GT);
+  CovariantLaplacianStencil<PeriodicGimplR,Field> CLst_GT(U_GT);
+
+  in_GT  = g*in;
+  out_GT = g*out_FLcs;
+
+  // Check M^GT_xy in_GT = g(x) M_xy g^dag(y) g(y) in = g(x) out(x)
+  CLcs_GT.M(in_GT,out_CLcs_GT);
+  CLst_GT.M(in_GT,out_CLst_GT);
+
+  diff = out_CLcs_GT - out_GT;
+  std:: cout << " Difference between Gauge xformed result and covariant Cshift Laplacian in xformed gauge  = " <<norm2(diff)<<std::endl;
+
+  diff = out_CLst_GT - out_GT;
+  std:: cout << " Difference between Gauge xformed result and covariant Stencil Laplacian in xformed gauge  = " <<norm2(diff)<<std::endl;
+  std:: cout << "--------------------------------------------------------" <<std::endl;
+
+
+  std:: cout << "******************************************************************" <<std::endl;
+  std:: cout << " Test C: compare in free Field to \"Feynman rule\"  " <<std::endl;
+  std:: cout << "******************************************************************" <<std::endl;
+
+  std::vector<int> dim_mask({1,1,1,0}); // 3d FFT
+  FFT theFFT(&Grid);
+  Field out(&Grid);
+  Field F_out(&Grid);
+  Field F_in(&Grid);
+
+  // FFT the random input vector
+  theFFT.FFT_dim_mask(F_in,in,dim_mask,FFT::forward);
+
+  // Convolution theorem: multiply by Fourier representation of (discrete) Laplacian to apply diff op
+  LatticeComplexD    lap(&Grid); lap = Zero();
+  LatticeComplexD    kmu(&Grid); 
+  ComplexD ci(0.0,1.0);
+  for(int mu=0;mu<3;mu++) {
+
+    RealD TwoPiL =  M_PI * 2.0/ latt_size[mu];
+
+    LatticeCoordinate(kmu,mu);
+    kmu = TwoPiL * kmu;
+
+    // (e^ik_mu + e^-ik_mu - 2) = 2( cos kmu - 1) ~ 2 (1 - k_mu^2/2 -1 ) =  - k_mu^2 + O(k^4)
+    lap = lap + 2.0*cos(kmu) - 2.0;
+    
+  }  
+  F_out = lap * F_in;
+  
+  // Inverse FFT the result
+  theFFT.FFT_dim_mask(out,F_out,dim_mask,FFT::backward);
+  
+  std::cout<<"Fourier xformed (in)             "<<norm2(F_in)<<std::endl;
+  std::cout<<"Fourier xformed Laplacian x (in) "<<norm2(F_out)<<std::endl;
+  
+  std::cout<<"Momentum space Laplacian application  "<< norm2(out)<<std::endl;
+  std::cout<<"Stencil Laplacian application         "<< norm2(out_CLcs)<<std::endl;
+    
+  diff = out_CLcs - out;
+  std::cout<<"diff "<< norm2(diff)<<std::endl;
+
+  Grid_finalize();
+}

examples/Example_Laplacian_smearing.cc

@@ -0,0 +1,127 @@
+#include <Grid/Grid.h>
+using namespace Grid;
+
+// Function used for Chebyshev smearing
+// 
+Real MomentumSmearing(Real p2)
+{
+  return (1 - 4.0*p2) * exp(-p2/4);
+}
+Real DistillationSmearing(Real p2)
+{
+  if ( p2 > 0.5 ) return 0.0;
+  else return 1.0;
+}
+
+// Flip sign to make prop to p^2, not -p^2 relative to last example
+template<class Gimpl,class Field> class CovariantLaplacianCshift : public SparseMatrixBase<Field>
+{
+public:
+  INHERIT_GIMPL_TYPES(Gimpl);
+
+  GridBase *grid;
+  GaugeField U;
+  
+  CovariantLaplacianCshift(GaugeField &_U)    :
+    grid(_U.Grid()),
+    U(_U) {  };
+
+  virtual GridBase *Grid(void) { return grid; };
+
+  virtual void  M    (const Field &in, Field &out)
+  {
+    out=Zero();
+    for(int mu=0;mu<Nd-1;mu++) {
+      GaugeLinkField Umu = PeekIndex<LorentzIndex>(U, mu); // NB: Inefficent
+      out = out - Gimpl::CovShiftForward(Umu,mu,in);    
+      out = out - Gimpl::CovShiftBackward(Umu,mu,in);    
+      out = out + 2.0*in;
+    }
+  };
+  virtual void  Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
+  virtual  void Mdiag    (const Field &in, Field &out)                  {assert(0);}; // Unimplemented need only for multigrid
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
+  virtual  void MdirAll  (const Field &in, std::vector<Field> &out)     {assert(0);}; // Unimplemented need only for multigrid
+};
+
+
+
+int main(int argc, char ** argv)
+{
+  Grid_init(&argc, &argv);
+
+  typedef LatticeColourVector Field;
+
+  auto latt_size   = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  auto mpi_layout  = GridDefaultMpi();
+
+  GridCartesian    Grid(latt_size,simd_layout,mpi_layout);
+  GridParallelRNG  RNG(&Grid);  RNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+
+
+  LatticeGaugeField U(&Grid);
+
+  SU<Nc>::ColdConfiguration(RNG,U);
+
+  typedef CovariantLaplacianCshift <PeriodicGimplR,Field> Laplacian_t;
+  Laplacian_t Laplacian(U);
+
+  
+  ColourVector ColourKronecker;
+  ColourKronecker = Zero();
+  ColourKronecker()()(0) = 1.0;
+
+  Coordinate site({latt_size[0]/2,
+                   latt_size[1]/2,
+                   latt_size[2]/2,
+		   0});
+
+  Field kronecker(&Grid);
+  kronecker = Zero();
+  pokeSite(ColourKronecker,kronecker,site);
+
+
+  Field psi(&Grid), chi(&Grid);
+
+  //////////////////////////////////////
+  // Classic Wuppertal smearing
+  //////////////////////////////////////
+
+  Integer Iterations = 80;
+  Real width = 2.0;
+  Real coeff = (width*width) / Real(4*Iterations);
+
+  chi=kronecker;
+  //  chi = (1-p^2/2N)^N kronecker
+  for(int n = 0; n < Iterations; ++n) {
+    Laplacian.M(chi,psi);
+    chi = chi - coeff*psi;
+  }
+
+  std::cout << " Wuppertal smeared operator is chi = \n" << chi <<std::endl;
+
+  /////////////////////////////////////
+  // Chebyshev smearing
+  /////////////////////////////////////
+  RealD lo = 0.0;
+  RealD hi = 12.0; // Analytic free field bound
+  HermitianLinearOperator<Laplacian_t,Field> HermOp(Laplacian);
+
+  std::cout << " Checking spectral range of our POSITIVE definite operator \n";
+  PowerMethod<Field> PM;
+  PM(HermOp,kronecker);
+  
+  //  Chebyshev<Field> ChebySmear(lo,hi,20,DistillationSmearing);
+  Chebyshev<Field> ChebySmear(lo,hi,20,MomentumSmearing);
+  {
+    std::ofstream of("chebysmear");
+    ChebySmear.csv(of);
+  }
+
+  ChebySmear(HermOp,kronecker,chi);
+  
+  std::cout << " Chebyshev smeared operator is chi = \n" << chi <<std::endl;
+
+  Grid_finalize();
+}

examples/Example_Laplacian_solver.cc

@@ -0,0 +1,129 @@
+#include <Grid/Grid.h>
+using namespace Grid;
+
+template<class Field>
+void SimpleConjugateGradient(LinearOperatorBase<Field> &HPDop,const Field &b, Field &x)
+{
+    RealD cp, c, alpha, d, beta, ssq, qq;
+    RealD Tolerance=1.0e-10;
+    int MaxIterations=10000;
+    
+    Field p(b), mmp(b), r(b);
+
+    HPDop.HermOpAndNorm(x, mmp, d, beta);
+    
+    r = b - mmp;
+    p = r;
+
+    cp = alpha = norm2(p);
+    ssq = norm2(b);
+
+    RealD rsq = Tolerance * Tolerance * ssq;
+
+    for (int k = 1; k <= MaxIterations; k++) {
+      c = cp;
+
+      HPDop.HermOp(p, mmp);
+
+      d = real(innerProduct(p,mmp));
+
+      alpha = c / d;
+
+      r = r - alpha *mmp;
+      cp = norm2(r);
+      beta = cp / c;
+
+      x   = x   + alpha* p ;
+      p   = r   + beta* p ;
+
+      std::cout << "iteration "<<k<<" cp " <<std::sqrt(cp/ssq) << std::endl;
+      if (cp <= rsq) {
+        return;
+      }
+    }
+    assert(0);
+}
+
+
+
+// Flip sign to make prop to p^2, not -p^2 relative to last example
+template<class Gimpl,class Field> class CovariantLaplacianCshift : public SparseMatrixBase<Field>
+{
+public:
+  INHERIT_GIMPL_TYPES(Gimpl);
+
+  GridBase *grid;
+  GaugeField U;
+  RealD m2=1.0e-2;
+  CovariantLaplacianCshift(GaugeField &_U)    :
+    grid(_U.Grid()),
+    U(_U) {  };
+
+  virtual GridBase *Grid(void) { return grid; };
+
+  virtual void  M    (const Field &in, Field &out)
+  {
+    out=Zero();
+    for(int mu=0;mu<Nd-1;mu++) {
+      GaugeLinkField Umu = PeekIndex<LorentzIndex>(U, mu); // NB: Inefficent
+      out = out - Gimpl::CovShiftForward(Umu,mu,in);    
+      out = out - Gimpl::CovShiftBackward(Umu,mu,in);    
+      out = out + 2.0*in + m2*in;
+    }
+  };
+  virtual void  Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
+  virtual  void Mdiag    (const Field &in, Field &out)                  {assert(0);}; // Unimplemented need only for multigrid
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
+  virtual  void MdirAll  (const Field &in, std::vector<Field> &out)     {assert(0);}; // Unimplemented need only for multigrid
+};
+
+
+
+int main(int argc, char ** argv)
+{
+  Grid_init(&argc, &argv);
+
+  typedef LatticeColourVector Field;
+
+  auto latt_size   = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  auto mpi_layout  = GridDefaultMpi();
+
+  GridCartesian    Grid(latt_size,simd_layout,mpi_layout);
+  GridParallelRNG  RNG(&Grid);  RNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+
+
+  LatticeGaugeField U(&Grid);
+
+  SU<Nc>::ColdConfiguration(RNG,U);
+
+  typedef CovariantLaplacianCshift <PeriodicGimplR,Field> Laplacian_t;
+  Laplacian_t Laplacian(U);
+
+  
+  ColourVector ColourKronecker;
+  ColourKronecker = Zero();
+  ColourKronecker()()(0) = 1.0;
+
+  Coordinate site({0,0,0,0}); // Point source at origin
+
+  Field kronecker(&Grid);
+  kronecker = Zero();
+  pokeSite(ColourKronecker,kronecker,site);
+
+  Field psi(&Grid); psi=Zero();
+
+  HermitianLinearOperator<Laplacian_t,Field> HermOp(Laplacian);
+  SimpleConjugateGradient(HermOp, kronecker,psi);
+
+  Field r(&Grid);
+  Laplacian.M(psi,r);
+  r=kronecker-r;
+  
+  std::cout << "True residual "<< norm2(r) <<std::endl;
+
+  // Optionally print the result vector
+  // std::cout << psi<<std::endl;
+
+  Grid_finalize();
+}

examples/Example_Mobius_spectrum.cc

@@ -0,0 +1,314 @@
+/*
+ * Warning: This code illustrative only: not well tested, and not meant for production use
+ * without regression / tests being applied
+ */
+
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+template<class Gimpl,class Field> class CovariantLaplacianCshift : public SparseMatrixBase<Field>
+{
+public:
+  INHERIT_GIMPL_TYPES(Gimpl);
+
+  GridBase *grid;
+  GaugeField U;
+  
+  CovariantLaplacianCshift(GaugeField &_U)    :
+    grid(_U.Grid()),
+    U(_U) {  };
+
+  virtual GridBase *Grid(void) { return grid; };
+
+  virtual void  M    (const Field &in, Field &out)
+  {
+    out=Zero();
+    for(int mu=0;mu<Nd-1;mu++) {
+      GaugeLinkField Umu = PeekIndex<LorentzIndex>(U, mu); // NB: Inefficent
+      out = out - Gimpl::CovShiftForward(Umu,mu,in);    
+      out = out - Gimpl::CovShiftBackward(Umu,mu,in);    
+      out = out + 2.0*in;
+    }
+  };
+  virtual void  Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
+  virtual  void Mdiag    (const Field &in, Field &out)                  {assert(0);}; // Unimplemented need only for multigrid
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
+  virtual  void MdirAll  (const Field &in, std::vector<Field> &out)     {assert(0);}; // Unimplemented need only for multigrid
+};
+
+void MakePhase(Coordinate mom,LatticeComplex &phase)
+{
+  GridBase *grid = phase.Grid();
+  auto latt_size = grid->GlobalDimensions();
+  ComplexD ci(0.0,1.0);
+  phase=Zero();
+
+  LatticeComplex coor(phase.Grid());
+  for(int mu=0;mu<Nd;mu++){
+    RealD TwoPiL =  M_PI * 2.0/ latt_size[mu];
+    LatticeCoordinate(coor,mu);
+    phase = phase + (TwoPiL * mom[mu]) * coor;
+  }
+  phase = exp(phase*ci);
+}
+void PointSource(Coordinate &coor,LatticePropagator &source)
+{
+  //  Coordinate coor({0,0,0,0});
+  source=Zero();
+  SpinColourMatrix kronecker; kronecker=1.0;
+  pokeSite(kronecker,source,coor);
+}
+void Z2WallSource(GridParallelRNG &RNG,int tslice,LatticePropagator &source)
+{
+  GridBase *grid = source.Grid();
+  LatticeComplex noise(grid);
+  LatticeComplex zz(grid); zz=Zero();
+  LatticeInteger t(grid);
+
+  RealD nrm=1.0/sqrt(2);
+  bernoulli(RNG, noise); // 0,1 50:50
+
+  noise = (2.*noise - Complex(1,1))*nrm;
+
+  LatticeCoordinate(t,Tdir);
+  noise = where(t==Integer(tslice), noise, zz);
+
+  source = 1.0;
+  source = source*noise;
+  std::cout << " Z2 wall " << norm2(source) << std::endl;
+}
+template<class Field>
+void GaussianSmear(LatticeGaugeField &U,Field &unsmeared,Field &smeared)
+{
+  typedef CovariantLaplacianCshift <PeriodicGimplR,Field> Laplacian_t;
+  Laplacian_t Laplacian(U);
+
+  Integer Iterations = 40;
+  Real width = 2.0;
+  Real coeff = (width*width) / Real(4*Iterations);
+
+  Field tmp(U.Grid());
+  smeared=unsmeared;
+  //  chi = (1-p^2/2N)^N kronecker
+  for(int n = 0; n < Iterations; ++n) {
+    Laplacian.M(smeared,tmp);
+    smeared = smeared - coeff*tmp;
+    std::cout << " smear iter " << n<<" " <<norm2(smeared)<<std::endl;
+  }
+}
+void GaussianSource(Coordinate &site,LatticeGaugeField &U,LatticePropagator &source)
+{
+  LatticePropagator tmp(source.Grid());
+  PointSource(site,source);
+  std::cout << " GaussianSource Kronecker "<< norm2(source)<<std::endl;
+  tmp = source;
+  GaussianSmear(U,tmp,source);
+  std::cout << " GaussianSource Smeared "<< norm2(source)<<std::endl;
+}
+void GaussianWallSource(GridParallelRNG &RNG,int tslice,LatticeGaugeField &U,LatticePropagator &source)
+{
+  Z2WallSource(RNG,tslice,source);
+  auto tmp = source;
+  GaussianSmear(U,tmp,source);
+}
+void SequentialSource(int tslice,Coordinate &mom,LatticePropagator &spectator,LatticePropagator &source)
+{
+  assert(mom.size()==Nd);
+  assert(mom[Tdir] == 0);
+
+  GridBase * grid = spectator.Grid();
+
+
+  LatticeInteger ts(grid);
+  LatticeCoordinate(ts,Tdir);
+  source = Zero();
+  source = where(ts==Integer(tslice),spectator,source); // Stick in a slice of the spectator, zero everywhere else
+
+  LatticeComplex phase(grid);
+  MakePhase(mom,phase);
+
+  source = source *phase;
+}
+template<class Action>
+void Solve(Action &D,LatticePropagator &source,LatticePropagator &propagator)
+{
+  GridBase *UGrid = D.GaugeGrid();
+  GridBase *FGrid = D.FermionGrid();
+
+  LatticeFermion src4  (UGrid); 
+  LatticeFermion src5  (FGrid); 
+  LatticeFermion result5(FGrid);
+  LatticeFermion result4(UGrid);
+  
+  ConjugateGradient<LatticeFermion> CG(1.0e-8,100000);
+  SchurRedBlackDiagMooeeSolve<LatticeFermion> schur(CG);
+  ZeroGuesser<LatticeFermion> ZG; // Could be a DeflatedGuesser if have eigenvectors
+  for(int s=0;s<Nd;s++){
+    for(int c=0;c<Nc;c++){
+      PropToFerm<Action>(src4,source,s,c);
+
+      D.ImportPhysicalFermionSource(src4,src5);
+
+      result5=Zero();
+      schur(D,src5,result5,ZG);
+      std::cout<<GridLogMessage
+	       <<"spin "<<s<<" color "<<c
+	       <<" norm2(src5d) "   <<norm2(src5)
+               <<" norm2(result5d) "<<norm2(result5)<<std::endl;
+
+      D.ExportPhysicalFermionSolution(result5,result4);
+
+      FermToProp<Action>(propagator,result4,s,c);
+    }
+  }
+}
+
+class MesonFile: Serializable {
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(MesonFile, std::vector<std::vector<Complex> >, data);
+};
+
+void MesonTrace(std::string file,LatticePropagator &q1,LatticePropagator &q2,LatticeComplex &phase)
+{
+  const int nchannel=4;
+  Gamma::Algebra Gammas[nchannel][2] = {
+    {Gamma::Algebra::Gamma5      ,Gamma::Algebra::Gamma5},
+    {Gamma::Algebra::GammaTGamma5,Gamma::Algebra::GammaTGamma5},
+    {Gamma::Algebra::GammaTGamma5,Gamma::Algebra::Gamma5},
+    {Gamma::Algebra::Gamma5      ,Gamma::Algebra::GammaTGamma5}
+  };
+
+  Gamma G5(Gamma::Algebra::Gamma5);
+
+  LatticeComplex meson_CF(q1.Grid());
+  MesonFile MF;
+
+  for(int ch=0;ch<nchannel;ch++){
+
+    Gamma Gsrc(Gammas[ch][0]);
+    Gamma Gsnk(Gammas[ch][1]);
+
+    meson_CF = trace(G5*adj(q1)*G5*Gsnk*q2*adj(Gsrc));
+
+    std::vector<TComplex> meson_T;
+    sliceSum(meson_CF,meson_T, Tdir);
+
+    int nt=meson_T.size();
+
+    std::vector<Complex> corr(nt);
+    for(int t=0;t<nt;t++){
+      corr[t] = TensorRemove(meson_T[t]); // Yes this is ugly, not figured a work around
+      std::cout << " channel "<<ch<<" t "<<t<<" " <<corr[t]<<std::endl;
+    }
+    MF.data.push_back(corr);
+  }
+
+  {
+    XmlWriter WR(file);
+    write(WR,"MesonFile",MF);
+  }
+}
+
+int main (int argc, char ** argv)
+{
+  const int Ls=8;
+
+  Grid_init(&argc,&argv);
+
+  // Double precision grids
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), 
+								   GridDefaultSimd(Nd,vComplex::Nsimd()),
+								   GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  //////////////////////////////////////////////////////////////////////
+  // You can manage seeds however you like.
+  // Recommend SeedUniqueString.
+  //////////////////////////////////////////////////////////////////////
+  std::vector<int> seeds4({1,2,3,4}); 
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+
+  LatticeGaugeField Umu(UGrid);
+  std::string config;
+  if( argc > 1 && argv[1][0] != '-' )
+  {
+    std::cout<<GridLogMessage <<"Loading configuration from "<<argv[1]<<std::endl;
+    FieldMetaData header;
+    NerscIO::readConfiguration(Umu, header, argv[1]);
+    config=argv[1];
+  }
+  else
+  {
+    std::cout<<GridLogMessage <<"Using hot configuration"<<std::endl;
+    SU<Nc>::ColdConfiguration(Umu);
+    //    SU<Nc>::HotConfiguration(RNG4,Umu);
+    config="HotConfig";
+  }
+
+  std::vector<RealD> masses({ 0.03,0.04,0.45} ); // u/d, s, c ??
+
+  int nmass = masses.size();
+
+  std::vector<MobiusFermionR *> FermActs;
+  
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
+  std::cout<<GridLogMessage <<"MobiusFermion action as Scaled Shamir kernel"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
+
+  for(auto mass: masses) {
+
+    RealD M5=1.0;
+    RealD b=1.5;// Scale factor b+c=2, b-c=1
+    RealD c=0.5;
+    
+    FermActs.push_back(new MobiusFermionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c));
+   
+  }
+
+  LatticePropagator point_source(UGrid);
+  LatticePropagator wall_source(UGrid);
+  LatticePropagator gaussian_source(UGrid);
+
+  Coordinate Origin({0,0,0,0});
+  PointSource   (Origin,point_source);
+  Z2WallSource  (RNG4,0,wall_source);
+  GaussianSource(Origin,Umu,gaussian_source);
+  
+  std::vector<LatticePropagator> PointProps(nmass,UGrid);
+  std::vector<LatticePropagator> GaussProps(nmass,UGrid);
+  std::vector<LatticePropagator> Z2Props   (nmass,UGrid);
+
+  for(int m=0;m<nmass;m++) {
+    
+    Solve(*FermActs[m],point_source   ,PointProps[m]);
+    Solve(*FermActs[m],gaussian_source,GaussProps[m]);
+    Solve(*FermActs[m],wall_source    ,Z2Props[m]);
+  
+  }
+
+  LatticeComplex phase(UGrid);
+  Coordinate mom({0,0,0,0});
+  MakePhase(mom,phase);
+  
+  for(int m1=0 ;m1<nmass;m1++) {
+  for(int m2=m1;m2<nmass;m2++) {
+    std::stringstream ssp,ssg,ssz;
+
+    ssp<<config<< "_m" << m1 << "_m"<< m2 << "_point_meson.xml";
+    ssg<<config<< "_m" << m1 << "_m"<< m2 << "_smeared_meson.xml";
+    ssz<<config<< "_m" << m1 << "_m"<< m2 << "_wall_meson.xml";
+
+    MesonTrace(ssp.str(),PointProps[m1],PointProps[m2],phase);
+    MesonTrace(ssg.str(),GaussProps[m1],GaussProps[m2],phase);
+    MesonTrace(ssz.str(),Z2Props[m1],Z2Props[m2],phase);
+  }}
+
+  Grid_finalize();
+}
+
+
+

examples/Makefile.am

@@ -0,0 +1,6 @@
+SUBDIRS = . 
+
+include Make.inc
+
+
+

scripts/filelist

@@ -82,3 +82,18 @@ for f in $TESTS; do
     echo >> Make.inc
 done
 cd ..
+
+# examples Make.inc
+cd $home/examples/
+echo> Make.inc
+TESTS=`ls *.cc`
+TESTLIST=`echo ${TESTS} | sed s/.cc//g `
+echo bin_PROGRAMS = ${TESTLIST} > Make.inc
+echo >> Make.inc
+for f in $TESTS; do
+    BNAME=`basename $f .cc`
+    echo ${BNAME}_SOURCES=$f  >> Make.inc
+    echo ${BNAME}_LDADD='$(top_builddir)/Grid/libGrid.a'>> Make.inc
+    echo >> Make.inc
+done
+cd ..

tests/debug/Test_heatbath_dwf_eofa.cc

@@ -66,7 +66,9 @@ 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);
@@ -84,7 +86,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, RNG5);
+    Meofa.refresh(Umu,sRNG, RNG5);
     printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
   }
 
@@ -94,7 +96,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, RNG5);
+    Meofa.refresh(Umu,sRNG, RNG5);
     printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
   }
 

tests/debug/Test_heatbath_dwf_eofa_gparity.cc

@@ -74,6 +74,9 @@ int main(int argc, char** argv)
   RNG5.SeedFixedIntegers(seeds5);
   GridParallelRNG RNG4(UGrid);
   RNG4.SeedFixedIntegers(seeds4);
+  GridSerialRNG sRNG;
+  RNG4.SeedFixedIntegers(seeds4);
+  sRNG.SeedFixedIntegers(seeds5);
 
   // Random gauge field
   LatticeGaugeField Umu(UGrid);
@@ -90,7 +93,7 @@ int main(int argc, char** argv)
     ConjugateGradient<FermionField> CG(1.0e-12, 5000);
     ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, false);
 
-    Meofa.refresh(Umu, RNG5);
+    Meofa.refresh(Umu,sRNG, RNG5);
     printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
   }
 
@@ -100,7 +103,7 @@ int main(int argc, char** argv)
     ConjugateGradient<FermionField> CG(1.0e-12, 5000);
     ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, true);
 
-    Meofa.refresh(Umu, RNG5);
+    Meofa.refresh(Umu,sRNG, RNG5);
     printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
   }
 

tests/debug/Test_heatbath_mobius_eofa.cc

@@ -68,8 +68,10 @@ 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);
   RNG5.SeedFixedIntegers(seeds5);
+  sRNG.SeedFixedIntegers(seeds5);
   GridParallelRNG RNG4(UGrid);
   RNG4.SeedFixedIntegers(seeds4);
 
@@ -86,7 +88,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, RNG5);
+    Meofa.refresh(Umu, sRNG,RNG5);
     printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
   }
 
@@ -96,7 +98,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, RNG5);
+    Meofa.refresh(Umu, sRNG,RNG5);
     printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
   }
 

tests/debug/Test_heatbath_mobius_eofa_gparity.cc

@@ -73,7 +73,9 @@ int main(int argc, char** argv)
   std::vector<int> seeds4({1, 2, 3, 4});
   std::vector<int> seeds5({5, 6, 7, 8});
   GridParallelRNG RNG5(FGrid);
+  GridSerialRNG   sRNG;
   RNG5.SeedFixedIntegers(seeds5);
+  sRNG.SeedFixedIntegers(seeds5);
   GridParallelRNG RNG4(UGrid);
   RNG4.SeedFixedIntegers(seeds4);
 
@@ -91,7 +93,7 @@ int main(int argc, char** argv)
     ConjugateGradient<FermionField> CG(1.0e-12, 5000);
     ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, false);
 
-    Meofa.refresh(Umu, RNG5);
+    Meofa.refresh(Umu, sRNG, RNG5);
     printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
   }
 
@@ -101,7 +103,7 @@ int main(int argc, char** argv)
     ConjugateGradient<FermionField> CG(1.0e-12, 5000);
     ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, true);
 
-    Meofa.refresh(Umu, RNG5);
+    Meofa.refresh(Umu, sRNG, RNG5);
     printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
   }
 

tests/forces/Test_gp_rect_force.cc

@@ -29,7 +29,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 using namespace std;
 using namespace Grid;
- ;
 
  
 
@@ -59,6 +58,10 @@ int main (int argc, char ** argv)
   double beta = 1.0;
   double c1   = 0.331;
 
+  const int nu = 1;
+  std::vector<int> twists(Nd,0);
+  twists[nu] = 1;
+  ConjugateGimplD::setDirections(twists);
   ConjugatePlaqPlusRectangleActionR Action(beta,c1);
   //ConjugateWilsonGaugeActionR Action(beta);
   //WilsonGaugeActionR Action(beta);

tests/forces/Test_momentum_filter.cc

@@ -61,7 +61,9 @@ int main (int argc, char ** argv)
   std::vector<int> seeds({1,2,3,4});
 
   GridParallelRNG          pRNG(&Grid);
+  GridSerialRNG            sRNG;
   pRNG.SeedFixedIntegers(seeds);
+  sRNG.SeedFixedIntegers(seeds);
 
   typedef PeriodicGimplR Gimpl;
   typedef WilsonGaugeAction<Gimpl> GaugeAction;
@@ -115,7 +117,7 @@ int main (int argc, char ** argv)
   
   integrator.setMomentumFilter(filter);
 
-  integrator.refresh(U, pRNG); //doesn't actually change the gauge field
+  integrator.refresh(U, sRNG, pRNG); //doesn't actually change the gauge field
 
   //Check the momentum is zero on the boundary
   const auto &P = integrator.getMomentum();

tests/lanczos/Test_dwf_block_lanczos.README

@@ -0,0 +1,73 @@
+#Example script 
+DIR=/gpfs/alpine/phy157/proj-shared/phy157dwf/chulwoo/Grid/BL/build/tests/lanczos
+BIN=${DIR}/Test_dwf_block_lanczos
+
+VOL='--grid 16.16.16.32 '
+GRID='--mpi 1.1.1.4 ' 
+CONF='--gconf ckpoint_lat.IEEE64BIG.2000 '
+OPT='--mass 0.01 --M5 1.8 --phase in.params --omega in.params --shm 4096' 
+#BL='--rbl 16.1024.128.1000.10 --split 1.1.4.4 --check_int 100 --resid 1.0e-5 --cheby_l 0.007 --cheby_u 7 --cheby_n 51'
+BL='--rbl 4.128.16.100.10 --split 1.1.1.4 --check_int 25 --resid 1.0e-5 --cheby_l 0.007 --cheby_u 7 --cheby_n 51'
+
+ARGS=${CONF}" "${OPT}" "${BL}" "${VOL}" "${GRID}
+export APP="${BIN}  ${ARGS}"
+echo APP=${APP}
+#export JS="jsrun --nrs 32 -a4 -g4 -c42 -dpacked -b packed:7 --smpiargs="-gpu" "
+export JS="jsrun --nrs 1 -a4 -g4 -c42 -dpacked -b  packed:10  --smpiargs="-gpu" "
+$JS  $APP
+
+#sample in.param 
+
+boundary_phase 0 1 0
+boundary_phase 1 1 0
+boundary_phase 2 1 0
+boundary_phase 3 -1 0
+
+omega 0 0.5 0
+omega 1 0.5 0
+omega 2 0.5 0
+omega 3 0.5 0
+omega 4 0.5 0
+omega 5 0.5 0
+omega 6 0.5 0
+omega 7 0.5 0
+omega 8 0.5 0
+omega 9 0.5 0
+omega 10 0.5 0
+omega 11 0.5 0
+
+
+#output 
+
+Grid : Message : 1.717474 s :  Gauge Configuration ckpoint_lat.IEEE64BIG.2000
+Grid : Message : 1.717478 s :  boundary_phase[0] = (1,0)
+Grid : Message : 1.717497 s :  boundary_phase[1] = (1,0)
+Grid : Message : 1.717500 s :  boundary_phase[2] = (1,0)
+Grid : Message : 1.717503 s :  boundary_phase[3] = (-1,0)
+Grid : Message : 1.717506 s :  Ls 12
+Grid : Message : 1.717507 s :  mass 0.01
+Grid : Message : 1.717510 s :  M5 1.8
+Grid : Message : 1.717512 s :  mob_b 1.5
+Grid : Message : 1.717514 s :  omega[0] = (0.5,0)
+Grid : Message : 1.717517 s :  omega[1] = (0.5,0)
+Grid : Message : 1.717520 s :  omega[2] = (0.5,0)
+Grid : Message : 1.717523 s :  omega[3] = (0.5,0)
+Grid : Message : 1.717526 s :  omega[4] = (0.5,0)
+Grid : Message : 1.717529 s :  omega[5] = (0.5,0)
+Grid : Message : 1.717532 s :  omega[6] = (0.5,0)
+Grid : Message : 1.717535 s :  omega[7] = (0.5,0)
+Grid : Message : 1.717538 s :  omega[8] = (0.5,0)
+Grid : Message : 1.717541 s :  omega[9] = (0.5,0)
+Grid : Message : 1.717544 s :  omega[10] = (0.5,0)
+Grid : Message : 1.717547 s :  omega[11] = (0.5,0)
+Grid : Message : 1.717550 s :  Nu 4
+Grid : Message : 1.717551 s :  Nk 128
+Grid : Message : 1.717552 s :  Np 16
+Grid : Message : 1.717553 s :  Nm 288
+Grid : Message : 1.717554 s :  Nstop 100
+Grid : Message : 1.717555 s :  Ntest 25
+Grid : Message : 1.717557 s :  MaxIter 10
+Grid : Message : 1.717558 s :  resid 1e-05
+Grid : Message : 1.717560 s :  Cheby Poly 0.007,7,51
+
+

tests/lanczos/Test_dwf_block_lanczos.cc

@@ -0,0 +1,403 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_dwf_block_lanczos.cc
+
+    Copyright (C) 2015
+
+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>
+#include <Grid/util/Init.h>
+#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h>
+
+using namespace std;
+using namespace Grid;
+//using namespace Grid::QCD;
+
+//typedef typename GparityDomainWallFermionR::FermionField FermionField;
+typedef typename ZMobiusFermionF::FermionField FermionField;
+
+RealD AllZero(RealD x){ return 0.;}
+
+class CmdJobParams 
+{
+  public:
+    std::string gaugefile;
+
+    int Ls;
+    double mass;
+    double M5;
+    double mob_b;
+    std::vector<ComplexD> omega;
+    std::vector<Complex> boundary_phase;
+    std::vector<int> mpi_split;
+    
+    LanczosType Impl;
+    int Nu;
+    int Nk;
+    int Np;
+    int Nm;
+    int Nstop;
+    int Ntest;
+    int MaxIter;
+    double resid;
+    
+    double low;
+    double high;
+    int order;
+
+    CmdJobParams()
+      : gaugefile("Hot"),
+        Ls(8), mass(0.01), M5(1.8), mob_b(1.5),
+        Impl(LanczosType::irbl),mpi_split(4,1),
+        Nu(4), Nk(200), Np(200), Nstop(100), Ntest(1), MaxIter(10), resid(1.0e-8), 
+        low(0.2), high(5.5), order(11)
+    {Nm=Nk+Np;};
+    
+    void Parse(char **argv, int argc);
+};
+
+
+void CmdJobParams::Parse(char **argv,int argc)
+{
+  std::string arg;
+  std::vector<int> vi;
+  double re,im;
+  int expect, idx;
+  std::string vstr;
+  std::ifstream pfile;
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--gconf") ){
+    gaugefile = GridCmdOptionPayload(argv,argv+argc,"--gconf");
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--phase") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--phase");
+    pfile.open(arg);
+    assert(pfile);
+    expect = 0;
+    while( pfile >> vstr ) {
+      if ( vstr.compare("boundary_phase") == 0 ) {
+        pfile >> vstr;
+        GridCmdOptionInt(vstr,idx);
+        assert(expect==idx);
+        pfile >> vstr;
+        GridCmdOptionFloat(vstr,re);
+        pfile >> vstr;
+        GridCmdOptionFloat(vstr,im);
+        boundary_phase.push_back({re,im});
+        expect++;
+      }
+    }
+    pfile.close();
+  } else {
+    for (int i=0; i<4; ++i) boundary_phase.push_back({1.,0.});
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--omega") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--omega");
+    pfile.open(arg);
+    assert(pfile);
+    Ls = 0;
+    while( pfile >> vstr ) {
+      if ( vstr.compare("omega") == 0 ) {
+        pfile >> vstr;
+        GridCmdOptionInt(vstr,idx);
+        assert(Ls==idx);
+        pfile >> vstr;
+        GridCmdOptionFloat(vstr,re);
+        pfile >> vstr;
+        GridCmdOptionFloat(vstr,im);
+        omega.push_back({re,im});
+        Ls++;
+      }
+    }
+    pfile.close();
+  } else {
+    if( GridCmdOptionExists(argv,argv+argc,"--Ls") ){
+      arg = GridCmdOptionPayload(argv,argv+argc,"--Ls");
+      GridCmdOptionInt(arg,Ls);
+    }
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--mass") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--mass");
+    GridCmdOptionFloat(arg,mass);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--M5") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--M5");
+    GridCmdOptionFloat(arg,M5);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--mob_b") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--mob_b");
+    GridCmdOptionFloat(arg,mob_b);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--irbl") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--irbl");
+    GridCmdOptionIntVector(arg,vi);
+    Nu = vi[0];
+    Nk = vi[1];
+    Np = vi[2];
+    Nstop = vi[3];
+    MaxIter = vi[4];
+    // ypj[fixme] mode overriding message is needed.
+    Impl = LanczosType::irbl;
+    Nm = Nk+Np;
+  }
+  
+  // block Lanczos with explicit extension of its dimensions
+  if( GridCmdOptionExists(argv,argv+argc,"--rbl") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--rbl");
+    GridCmdOptionIntVector(arg,vi);
+    Nu = vi[0];
+    Nk = vi[1];
+    Np = vi[2]; // vector space is enlarged by adding Np vectors
+    Nstop = vi[3];
+    MaxIter = vi[4];
+    // ypj[fixme] mode overriding message is needed.
+    Impl = LanczosType::rbl;
+    Nm = Nk+Np*MaxIter;
+  }
+  
+#if 1
+  // block Lanczos with explicit extension of its dimensions
+  if( GridCmdOptionExists(argv,argv+argc,"--split") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--split");
+    GridCmdOptionIntVector(arg,vi);
+    for(int i=0;i<mpi_split.size();i++)
+    mpi_split[i] = vi[i];
+  }
+#endif
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--check_int") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--check_int");
+    GridCmdOptionInt(arg,Ntest);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--resid") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--resid");
+    GridCmdOptionFloat(arg,resid);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--cheby_l") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_l");
+    GridCmdOptionFloat(arg,low);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--cheby_u") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_u");
+    GridCmdOptionFloat(arg,high);
+  }
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--cheby_n") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--cheby_n");
+    GridCmdOptionInt(arg,order);
+  }
+  
+  if ( CartesianCommunicator::RankWorld() == 0 ) {
+    std::streamsize ss = std::cout.precision();
+    std::cout << GridLogMessage <<" Gauge Configuration "<< gaugefile << '\n';
+    std::cout.precision(15);
+    for ( int i=0; i<4; ++i ) std::cout << GridLogMessage <<" boundary_phase["<< i << "] = " << boundary_phase[i] << '\n';
+    std::cout.precision(ss);
+    std::cout << GridLogMessage <<" Ls "<< Ls << '\n';
+    std::cout << GridLogMessage <<" mass "<< mass << '\n';
+    std::cout << GridLogMessage <<" M5 "<< M5 << '\n';
+    std::cout << GridLogMessage <<" mob_b "<< mob_b << '\n';
+    std::cout.precision(15);
+    for ( int i=0; i<Ls; ++i ) std::cout << GridLogMessage <<" omega["<< i << "] = " << omega[i] << '\n';
+    std::cout.precision(ss);
+    std::cout << GridLogMessage <<" Nu "<< Nu << '\n'; 
+    std::cout << GridLogMessage <<" Nk "<< Nk << '\n'; 
+    std::cout << GridLogMessage <<" Np "<< Np << '\n'; 
+    std::cout << GridLogMessage <<" Nm "<< Nm << '\n'; 
+    std::cout << GridLogMessage <<" Nstop "<< Nstop << '\n'; 
+    std::cout << GridLogMessage <<" Ntest "<< Ntest << '\n'; 
+    std::cout << GridLogMessage <<" MaxIter "<< MaxIter << '\n'; 
+    std::cout << GridLogMessage <<" resid "<< resid << '\n'; 
+    std::cout << GridLogMessage <<" Cheby Poly "<< low << "," << high << "," << order << std::endl; 
+  }
+}
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+  
+  CmdJobParams JP;
+  JP.Parse(argv,argc);
+
+  GridCartesian         * UGridD   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(JP.Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(JP.Ls,UGrid);
+
+  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);  RNG4.SeedFixedIntegers(seeds4);
+  // ypj [note] why seed RNG5 again? bug? In this case, run with a default seed().
+  GridParallelRNG          RNG5rb(FrbGrid);  RNG5rb.SeedFixedIntegers(seeds5);
+
+  LatticeGaugeField UmuD(UGridD); 
+  LatticeGaugeFieldF Umu(UGrid); 
+  std::vector<LatticeColourMatrixF> U(4,UGrid);
+  
+  if ( JP.gaugefile.compare("Hot") == 0 ) {
+    SU3::HotConfiguration(RNG4, UmuD);
+  } else {
+    FieldMetaData header;
+    NerscIO::readConfiguration(UmuD,header,JP.gaugefile);
+    // ypj [fixme] additional checks for the loaded configuration?
+  }
+  precisionChange(Umu,UmuD);
+  
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+  }
+  
+  RealD mass = JP.mass;
+  RealD M5 = JP.M5;
+
+// ypj [fixme] flexible support for a various Fermions
+//  RealD mob_b = JP.mob_b;      // Gparity
+//  std::vector<ComplexD> omega; // ZMobius
+  
+//  GparityMobiusFermionD ::ImplParams params;
+//  std::vector<int> twists({1,1,1,0});
+//  params.twists = twists;
+//  GparityMobiusFermionR  Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params);
+//  SchurDiagTwoOperator<GparityMobiusFermionR,FermionField> HermOp(Ddwf);
+
+
+//  int mrhs = JP.Nu;
+  int Ndir=4;
+  auto mpi_layout  = GridDefaultMpi();
+  std::vector<int> mpi_split (Ndir,1);
+#if 0
+    int tmp=mrhs, dir=0;
+    std::cout << GridLogMessage  << "dir= "<<dir <<"tmp= "<<tmp<<"mpi_split= "<<mpi_split[dir]<<"mpi_layout= "<<mpi_split[dir]<<std::endl;
+    while ( tmp> 1) {
+    if ((mpi_split[dir]*2) <= mpi_layout[dir]){
+        mpi_split[dir] *=2;
+        tmp = tmp/2;
+    }
+    std::cout << GridLogMessage  << "dir= "<<dir <<"tmp= "<<tmp<<"mpi_split= "<<mpi_split[dir]<<"mpi_layout= "<<mpi_layout[dir]<<std::endl;
+        dir = (dir+1)%Ndir;
+    }
+#endif
+    int mrhs=1;
+    for(int i =0;i<Ndir;i++){
+      mpi_split[i] = mpi_layout[i] / JP.mpi_split[i] ;
+      mrhs *= JP.mpi_split[i];
+    }
+    std::cout << GridLogMessage  << "mpi_layout= " << mpi_layout << std::endl;
+    std::cout << GridLogMessage  << "mpi_split= " << mpi_split << std::endl;
+    std::cout << GridLogMessage  << "mrhs= " << mrhs << std::endl;
+//    assert(JP.Nu==tmp);
+
+  /////////////////////////////////////////////
+  // Split into 1^4 mpi communicators
+  /////////////////////////////////////////////
+  GridCartesian         * SGrid = new GridCartesian(GridDefaultLatt(),
+                                                    GridDefaultSimd(Nd,vComplexF::Nsimd()),
+                                                    mpi_split,
+                                                    *UGrid);
+
+  GridCartesian         * SFGrid   = SpaceTimeGrid::makeFiveDimGrid(JP.Ls,SGrid);
+  GridRedBlackCartesian * SrbGrid  = SpaceTimeGrid::makeFourDimRedBlackGrid(SGrid);
+  GridRedBlackCartesian * SFrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(JP.Ls,SGrid);
+
+  LatticeGaugeFieldF s_Umu(SGrid);
+  Grid_split  (Umu,s_Umu);
+
+  //WilsonFermionR::ImplParams params;
+  ZMobiusFermionF::ImplParams params;
+  params.overlapCommsCompute = true;
+  params.boundary_phases = JP.boundary_phase;
+  ZMobiusFermionF  Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,JP.omega,1.,0.,params);
+//  SchurDiagTwoOperator<ZMobiusFermionR,FermionField> HermOp(Ddwf);
+  SchurDiagOneOperator<ZMobiusFermionF,FermionField> HermOp(Ddwf);
+  ZMobiusFermionF  Dsplit(s_Umu,*SFGrid,*SFrbGrid,*SGrid,*SrbGrid,mass,M5,JP.omega,1.,0.,params);
+//  SchurDiagTwoOperator<ZMobiusFermionR,FermionField> SHermOp(Dsplit);
+  SchurDiagOneOperator<ZMobiusFermionF,FermionField> SHermOp(Dsplit);
+
+  //std::vector<double> Coeffs { 0.,-1.}; 
+  // ypj [note] this may not be supported by some compilers
+  std::vector<double> Coeffs({ 0.,-1.}); 
+  Polynomial<FermionField> PolyX(Coeffs);
+  //Chebyshev<FermionField> Cheb(0.2,5.5,11);
+  Chebyshev<FermionField> Cheb(JP.low,JP.high,JP.order);
+//  Cheb.csv(std::cout);
+  ImplicitlyRestartedBlockLanczos<FermionField> IRBL(HermOp, SHermOp,
+						     FrbGrid,SFrbGrid,mrhs,
+                                                     Cheb,
+                                                     JP.Nstop, JP.Ntest,
+                                                     JP.Nu, JP.Nk, JP.Nm,
+                                                     JP.resid,
+                                                     JP.MaxIter,
+						     IRBLdiagonaliseWithEigen);
+//						     IRBLdiagonaliseWithLAPACK);
+  IRBL.split_test=0;
+  
+  std::vector<RealD> eval(JP.Nm);
+  
+  std::vector<FermionField> src(JP.Nu,FrbGrid);
+if (0)
+{
+// in case RNG is too slow
+  std::cout << GridLogMessage << "Using RNG5"<<std::endl;
+  FermionField src_tmp(FGrid);
+  for ( int i=0; i<JP.Nu; ++i ){
+//    gaussian(RNG5,src_tmp);
+     ComplexD rnd;
+     RealD re;
+     fillScalar(re,RNG5._gaussian[0],RNG5._generators[0]);
+    std::cout << i <<" / "<< JP.Nm  <<" re "<< re  << std::endl;
+// printf("%d / %d re %e\n",i,FGrid->_processor,re);
+    src_tmp=re;
+    pickCheckerboard(Odd,src[i],src_tmp);
+  }
+  RNG5.Report();
+} else {
+  std::cout << GridLogMessage << "Using RNG5rb"<<std::endl;
+  for ( int i=0; i<JP.Nu; ++i )
+    gaussian(RNG5rb,src[i]);
+  RNG5rb.Report();
+
+}
+  
+  std::vector<FermionField> evec(JP.Nm,FrbGrid);
+  for(int i=0;i<1;++i){
+    std::cout << GridLogMessage << i <<" / "<< JP.Nm <<" grid pointer "<< evec[i].Grid() << std::endl;
+  };
+
+  int Nconv;
+  IRBL.calc(eval,evec,src,Nconv,JP.Impl);
+
+
+  Grid_finalize();
+}

tests/lanczos/dwf.batch

@@ -0,0 +1,27 @@
+# run script on Cori
+#!/bin/bash
+#SBATCH  -t 30
+#SBATCH  -N 16
+#SBATCH  -C knl
+#SBATCH  -A mp13
+#SBATCH  --ntasks-per-node=1
+
+export OMP_PROC_BIND=true
+export OMP_PLACES=threads
+export OMP_NUM_THREADS=64
+
+#BIN=benchmarks/Benchmark_ITT
+#OPT='--cpu-bind=cores' 
+BIN=./Test_dwf_block_lanczos
+rundir=.
+CONF='--gconf '${rundir}'/ckpoint_lat.IEEE64BIG.5000'
+VOL='--grid 16.16.16.32'
+GRID='--mpi 1.2.2.4'
+OPT='--mass 0.00054 --M5 1.8 --phase in.params --omega in.params --shm 4096'
+BL='--rbl 4.32.32.30.7 --split 1.2.2.1 --check_int 8 --resid 1.0e-5 --cheby_l 0.0027 --cheby_u 7 --cheby_n 51'
+
+#echo srun $OPT $BIN $gridoptions $jobparams
+#srun $OPT $BIN $gridoptions $jobparams
+
+echo srun $BIN $CONF $OPT $BL $VOL $GRID
+srun $BIN $CONF $OPT $BL $VOL $GRID

tests/lanczos/in.params

@@ -0,0 +1,17 @@
+boundary_phase 0 1 0
+boundary_phase 1 1 0
+boundary_phase 2 1 0
+boundary_phase 3 -1 0
+
+omega 0 0.375 0
+omega 1 0.375 0
+omega 2 0.375 0
+omega 3 0.375 0
+omega 4 0.375 0
+omega 5 0.375 0
+omega 6 0.375 0
+omega 7 0.375 0
+omega 8 0.375 0
+omega 9 0.375 0
+omega 10 0.375 0
+omega 11 0.375 0