Merge pull request #396 from fjosw/fix/readd_config.h

fix: readded Config.h and Version.h to HFILEs in Grid/Makefile.am

Grid/lattice/Lattice_reduction.h

@@ -142,6 +142,15 @@ inline typename vobj::scalar_objectD sumD(const vobj *arg, Integer osites)
   return sumD_cpu(arg,osites);
 #endif  
 }
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_large(const vobj *arg, Integer osites)
+{
+#if defined(GRID_CUDA)||defined(GRID_HIP)
+  return sumD_gpu_large(arg,osites);
+#else
+  return sumD_cpu(arg,osites);
+#endif  
+}
 
 template<class vobj>
 inline typename vobj::scalar_object sum(const Lattice<vobj> &arg)
@@ -159,6 +168,22 @@ inline typename vobj::scalar_object sum(const Lattice<vobj> &arg)
   return ssum;
 }
 
+template<class vobj>
+inline typename vobj::scalar_object sum_large(const Lattice<vobj> &arg)
+{
+#if defined(GRID_CUDA)||defined(GRID_HIP)
+  autoView( arg_v, arg, AcceleratorRead);
+  Integer osites = arg.Grid()->oSites();
+  auto ssum= sum_gpu_large(&arg_v[0],osites);
+#else
+  autoView(arg_v, arg, CpuRead);
+  Integer osites = arg.Grid()->oSites();
+  auto ssum= sum_cpu(&arg_v[0],osites);
+#endif
+  arg.Grid()->GlobalSum(ssum);
+  return ssum;
+}
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Deterministic Reduction operations
 ////////////////////////////////////////////////////////////////////////////////////////////////////

Grid/lattice/Lattice_reduction_gpu.h

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

Grid/parallelIO/IldgIO.h

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

Grid/parallelIO/MetaData.h

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

Grid/parallelIO/NerscIO.h

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

Grid/qcd/action/fermion/CayleyFermion5D.h

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

Grid/qcd/action/fermion/CloverHelpers.h

@@ -0,0 +1,432 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./lib/qcd/action/fermion/WilsonCloverFermionImplementation.h
+
+    Copyright (C) 2017 - 2022
+
+    Author: paboyle <paboyle@ph.ed.ac.uk>
+    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
+    Author: Mattia Bruno <mattia.bruno@cern.ch>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+/*  END LEGAL */
+
+#pragma once
+
+#include <Grid/Grid.h>
+#include <Grid/qcd/spin/Dirac.h>
+#include <Grid/qcd/action/fermion/WilsonCloverHelpers.h>
+
+////////////////////////////////////////////
+// Standard Clover
+//   (4+m0) + csw * clover_term
+// Exp Clover
+//   (4+m0) * exp(csw/(4+m0) clover_term)
+//   = (4+m0) + csw * clover_term + ...
+////////////////////////////////////////////
+
+NAMESPACE_BEGIN(Grid);
+
+
+//////////////////////////////////
+// Generic Standard Clover
+//////////////////////////////////
+
+template<class Impl>
+class CloverHelpers: public WilsonCloverHelpers<Impl> {
+public:
+
+  INHERIT_IMPL_TYPES(Impl);
+  INHERIT_CLOVER_TYPES(Impl);
+
+  typedef WilsonCloverHelpers<Impl> Helpers;
+
+  static void Instantiate(CloverField& CloverTerm, CloverField& CloverTermInv, RealD csw_t, RealD diag_mass) {
+    GridBase *grid = CloverTerm.Grid();
+    CloverTerm += diag_mass;
+
+    int lvol = grid->lSites();
+    int DimRep = Impl::Dimension;
+    {
+      autoView(CTv,CloverTerm,CpuRead);
+      autoView(CTIv,CloverTermInv,CpuWrite);
+      thread_for(site, lvol, {
+        Coordinate lcoor;
+        grid->LocalIndexToLocalCoor(site, lcoor);
+        Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
+        Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
+        typename SiteClover::scalar_object Qx = Zero(), Qxinv = Zero();
+        peekLocalSite(Qx, CTv, lcoor);
+
+        for (int j = 0; j < Ns; j++)
+          for (int k = 0; k < Ns; k++)
+            for (int a = 0; a < DimRep; a++)
+              for (int b = 0; b < DimRep; b++){
+                auto zz =  Qx()(j, k)(a, b);
+                EigenCloverOp(a + j * DimRep, b + k * DimRep) = std::complex<double>(zz);
+              }
+
+        EigenInvCloverOp = EigenCloverOp.inverse();
+        for (int j = 0; j < Ns; j++)
+          for (int k = 0; k < Ns; k++)
+            for (int a = 0; a < DimRep; a++)
+              for (int b = 0; b < DimRep; b++)
+                Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep);
+               pokeLocalSite(Qxinv, CTIv, lcoor);
+      });
+    }
+  }
+
+  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu) {
+    return Helpers::Cmunu(U, lambda, mu, nu);
+  }
+
+};
+
+
+//////////////////////////////////
+// Generic Exp Clover
+//////////////////////////////////
+
+template<class Impl>
+class ExpCloverHelpers: public WilsonCloverHelpers<Impl> {
+public:
+
+  INHERIT_IMPL_TYPES(Impl);
+  INHERIT_CLOVER_TYPES(Impl);
+
+  template <typename vtype> using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
+  typedef WilsonCloverHelpers<Impl> Helpers;
+
+  // Can this be avoided?
+  static void IdentityTimesC(const CloverField& in, RealD c) {
+    int DimRep = Impl::Dimension;
+
+    autoView(in_v, in, AcceleratorWrite);
+
+    accelerator_for(ss, in.Grid()->oSites(), 1, {
+      for (int sa=0; sa<Ns; sa++)
+        for (int ca=0; ca<DimRep; ca++)
+          in_v[ss]()(sa,sa)(ca,ca) = c;
+    });
+  }
+
+  static int getNMAX(RealD prec, RealD R) {
+    /* compute stop condition for exponential */
+    int NMAX=1;
+    RealD cond=R*R/2.;
+
+    while (cond*std::exp(R)>prec) {
+      NMAX++;
+      cond*=R/(double)(NMAX+1);
+    }
+    return NMAX;
+  }
+
+  static int getNMAX(Lattice<iImplClover<vComplexD>> &t, RealD R) {return getNMAX(1e-12,R);}
+  static int getNMAX(Lattice<iImplClover<vComplexF>> &t, RealD R) {return getNMAX(1e-6,R);}
+
+  static void Instantiate(CloverField& Clover, CloverField& CloverInv, RealD csw_t, RealD diag_mass) {
+    GridBase* grid = Clover.Grid();
+    CloverField ExpClover(grid);
+
+    int NMAX = getNMAX(Clover, 3.*csw_t/diag_mass);
+
+    Clover *= (1.0/diag_mass);
+
+    // Taylor expansion, slow but generic
+    // Horner scheme: a0 + a1 x + a2 x^2 + .. = a0 + x (a1 + x(...))
+    // qN = cN
+    // qn = cn + qn+1 X
+    std::vector<RealD> cn(NMAX+1);
+    cn[0] = 1.0;
+    for (int i=1; i<=NMAX; i++)
+      cn[i] = cn[i-1] / RealD(i);
+
+    ExpClover = Zero();
+    IdentityTimesC(ExpClover, cn[NMAX]);
+    for (int i=NMAX-1; i>=0; i--)
+      ExpClover = ExpClover * Clover + cn[i];
+
+    // prepare inverse
+    CloverInv = (-1.0)*Clover;
+
+    Clover = ExpClover * diag_mass;
+
+    ExpClover = Zero();
+    IdentityTimesC(ExpClover, cn[NMAX]);
+    for (int i=NMAX-1; i>=0; i--)
+      ExpClover = ExpClover * CloverInv + cn[i];
+
+    CloverInv = ExpClover * (1.0/diag_mass);
+
+  }
+
+  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu) {
+    assert(0);
+  }
+
+};
+
+
+//////////////////////////////////
+// Compact Standard Clover
+//////////////////////////////////
+
+
+template<class Impl>
+class CompactCloverHelpers: public CompactWilsonCloverHelpers<Impl>,
+                            public WilsonCloverHelpers<Impl> {
+public:
+
+  INHERIT_IMPL_TYPES(Impl);
+  INHERIT_CLOVER_TYPES(Impl);
+  INHERIT_COMPACT_CLOVER_TYPES(Impl);
+
+  typedef WilsonCloverHelpers<Impl> Helpers;
+  typedef CompactWilsonCloverHelpers<Impl> CompactHelpers;
+
+  static void MassTerm(CloverField& Clover, RealD diag_mass) {
+    Clover += diag_mass;
+  }
+
+  static void Exponentiate_Clover(CloverDiagonalField& Diagonal,
+                          CloverTriangleField& Triangle,
+                          RealD csw_t, RealD diag_mass) {
+
+    // Do nothing
+  }
+
+  // TODO: implement Cmunu for better performances with compact layout, but don't do it
+  // here, but rather in WilsonCloverHelpers.h -> CompactWilsonCloverHelpers
+  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu) {
+    return Helpers::Cmunu(U, lambda, mu, nu);
+  }
+};
+
+//////////////////////////////////
+// Compact Exp Clover
+//////////////////////////////////
+
+template<class Impl>
+class CompactExpCloverHelpers: public CompactWilsonCloverHelpers<Impl> {
+public:
+
+  INHERIT_IMPL_TYPES(Impl);
+  INHERIT_CLOVER_TYPES(Impl);
+  INHERIT_COMPACT_CLOVER_TYPES(Impl);
+
+  template <typename vtype> using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
+  typedef CompactWilsonCloverHelpers<Impl> CompactHelpers;
+
+  static void MassTerm(CloverField& Clover, RealD diag_mass) {
+    // do nothing!
+    // mass term is multiplied to exp(Clover) below
+  }
+
+  static int getNMAX(RealD prec, RealD R) {
+    /* compute stop condition for exponential */
+    int NMAX=1;
+    RealD cond=R*R/2.;
+
+    while (cond*std::exp(R)>prec) {
+      NMAX++;
+      cond*=R/(double)(NMAX+1);
+    }
+    return NMAX;
+  }
+
+  static int getNMAX(Lattice<iImplCloverDiagonal<vComplexD>> &t, RealD R) {return getNMAX(1e-12,R);}
+  static int getNMAX(Lattice<iImplCloverDiagonal<vComplexF>> &t, RealD R) {return getNMAX(1e-6,R);}
+
+  static void ExponentiateHermitean6by6(const iMatrix<ComplexD,6> &arg, const RealD& alpha, const std::vector<RealD>& cN, const int Niter, iMatrix<ComplexD,6>& dest){
+
+  	  typedef iMatrix<ComplexD,6> mat;
+
+  	  RealD qn[6];
+  	  RealD qnold[6];
+  	  RealD p[5];
+  	  RealD trA2, trA3, trA4;
+
+  	  mat A2, A3, A4, A5;
+  	  A2 = alpha * alpha * arg * arg;
+  	  A3 = alpha * arg * A2;
+  	  A4 = A2 * A2;
+  	  A5 = A2 * A3;
+
+  	  trA2 = toReal( trace(A2) );
+  	  trA3 = toReal( trace(A3) );
+  	  trA4 = toReal( trace(A4));
+
+  	  p[0] = toReal( trace(A3 * A3)) / 6.0 - 0.125 * trA4 * trA2 - trA3 * trA3 / 18.0 + trA2 * trA2 * trA2/ 48.0;
+  	  p[1] = toReal( trace(A5)) / 5.0 - trA3 * trA2 / 6.0;
+  	  p[2] = toReal( trace(A4)) / 4.0 - 0.125 * trA2 * trA2;
+  	  p[3] = trA3 / 3.0;
+  	  p[4] = 0.5 * trA2;
+
+  	  qnold[0] = cN[Niter];
+  	  qnold[1] = 0.0;
+  	  qnold[2] = 0.0;
+  	  qnold[3] = 0.0;
+  	  qnold[4] = 0.0;
+  	  qnold[5] = 0.0;
+
+  	  for(int i = Niter-1; i >= 0; i--)
+  	  {
+  	   qn[0] = p[0] * qnold[5] + cN[i];
+  	   qn[1] = p[1] * qnold[5] + qnold[0];
+  	   qn[2] = p[2] * qnold[5] + qnold[1];
+  	   qn[3] = p[3] * qnold[5] + qnold[2];
+  	   qn[4] = p[4] * qnold[5] + qnold[3];
+  	   qn[5] = qnold[4];
+
+  	   qnold[0] = qn[0];
+  	   qnold[1] = qn[1];
+  	   qnold[2] = qn[2];
+  	   qnold[3] = qn[3];
+  	   qnold[4] = qn[4];
+  	   qnold[5] = qn[5];
+  	  }
+
+  	  mat unit(1.0);
+
+  	  dest = (qn[0] * unit + qn[1] * alpha * arg + qn[2] * A2 + qn[3] * A3 + qn[4] * A4 + qn[5] * A5);
+
+    }
+
+  static void Exponentiate_Clover(CloverDiagonalField& Diagonal, CloverTriangleField& Triangle, RealD csw_t, RealD diag_mass) {
+
+    GridBase* grid = Diagonal.Grid();
+    int NMAX = getNMAX(Diagonal, 3.*csw_t/diag_mass);
+
+    //
+    // Implementation completely in Daniel's layout
+    //
+
+    // Taylor expansion with Cayley-Hamilton recursion
+    // underlying Horner scheme as above
+    std::vector<RealD> cn(NMAX+1);
+    cn[0] = 1.0;
+    for (int i=1; i<=NMAX; i++){
+      cn[i] = cn[i-1] / RealD(i);
+    }
+
+      // Taken over from Daniel's implementation
+      conformable(Diagonal, Triangle);
+
+      long lsites = grid->lSites();
+
+      typedef typename SiteCloverDiagonal::scalar_object scalar_object_diagonal;
+      typedef typename SiteCloverTriangle::scalar_object scalar_object_triangle;
+      typedef iMatrix<ComplexD,6> mat;
+
+      autoView(diagonal_v,  Diagonal,  CpuRead);
+      autoView(triangle_v,  Triangle,  CpuRead);
+      autoView(diagonalExp_v, Diagonal, CpuWrite);
+      autoView(triangleExp_v, Triangle, CpuWrite);
+
+      thread_for(site, lsites, { // NOTE: Not on GPU because of (peek/poke)LocalSite
+
+    	  mat srcCloverOpUL(0.0); // upper left block
+    	  mat srcCloverOpLR(0.0); // lower right block
+    	  mat ExpCloverOp;
+
+        scalar_object_diagonal diagonal_tmp     = Zero();
+        scalar_object_diagonal diagonal_exp_tmp = Zero();
+        scalar_object_triangle triangle_tmp     = Zero();
+        scalar_object_triangle triangle_exp_tmp = Zero();
+
+        Coordinate lcoor;
+        grid->LocalIndexToLocalCoor(site, lcoor);
+
+        peekLocalSite(diagonal_tmp, diagonal_v, lcoor);
+        peekLocalSite(triangle_tmp, triangle_v, lcoor);
+
+        int block;
+        block = 0;
+        for(int i = 0; i < 6; i++){
+        	for(int j = 0; j < 6; j++){
+        		if (i == j){
+        			srcCloverOpUL(i,j) = static_cast<ComplexD>(TensorRemove(diagonal_tmp()(block)(i)));
+        		}
+        		else{
+        			srcCloverOpUL(i,j) = static_cast<ComplexD>(TensorRemove(CompactHelpers::triangle_elem(triangle_tmp, block, i, j)));
+        		}
+        	}
+        }
+        block = 1;
+        for(int i = 0; i < 6; i++){
+          	for(int j = 0; j < 6; j++){
+           		if (i == j){
+           			srcCloverOpLR(i,j) = static_cast<ComplexD>(TensorRemove(diagonal_tmp()(block)(i)));
+           		}
+           		else{
+           			srcCloverOpLR(i,j) = static_cast<ComplexD>(TensorRemove(CompactHelpers::triangle_elem(triangle_tmp, block, i, j)));
+           		}
+            }
+        }
+
+        // exp(Clover)
+
+        ExponentiateHermitean6by6(srcCloverOpUL,1.0/diag_mass,cn,NMAX,ExpCloverOp);
+
+        block = 0;
+        for(int i = 0; i < 6; i++){
+        	for(int j = 0; j < 6; j++){
+            	if (i == j){
+            		diagonal_exp_tmp()(block)(i) = ExpCloverOp(i,j);
+            	}
+            	else if(i < j){
+            		triangle_exp_tmp()(block)(CompactHelpers::triangle_index(i, j)) = ExpCloverOp(i,j);
+            	}
+           	}
+        }
+
+        ExponentiateHermitean6by6(srcCloverOpLR,1.0/diag_mass,cn,NMAX,ExpCloverOp);
+
+        block = 1;
+        for(int i = 0; i < 6; i++){
+        	for(int j = 0; j < 6; j++){
+              	if (i == j){
+              		diagonal_exp_tmp()(block)(i) = ExpCloverOp(i,j);
+               	}
+               	else if(i < j){
+               		triangle_exp_tmp()(block)(CompactHelpers::triangle_index(i, j)) = ExpCloverOp(i,j);
+               	}
+            }
+        }
+
+        pokeLocalSite(diagonal_exp_tmp, diagonalExp_v, lcoor);
+        pokeLocalSite(triangle_exp_tmp, triangleExp_v, lcoor);
+      });
+
+    Diagonal = Diagonal * diag_mass;
+    Triangle = Triangle * diag_mass;
+  }
+
+
+  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu) {
+    assert(0);
+  }
+
+};
+
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/CompactWilsonCloverFermion.h

@@ -31,6 +31,7 @@
 
 #include <Grid/qcd/action/fermion/WilsonCloverTypes.h>
 #include <Grid/qcd/action/fermion/WilsonCloverHelpers.h>
+#include <Grid/qcd/action/fermion/CloverHelpers.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -85,7 +86,7 @@ NAMESPACE_BEGIN(Grid);
 //           + (2 * 1 + 4 * 1/2) triangle parts = 4 triangle parts =  60 complex words per site
 //                                                                 =  84 complex words per site
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
 class CompactWilsonCloverFermion : public WilsonFermion<Impl>,
                                    public WilsonCloverHelpers<Impl>,
                                    public CompactWilsonCloverHelpers<Impl> {

Grid/qcd/action/fermion/Fermion.h

@@ -138,38 +138,52 @@ typedef WilsonTMFermion<WilsonImplF> WilsonTMFermionF;
 typedef WilsonTMFermion<WilsonImplD> WilsonTMFermionD;
 
 // Clover fermions
-typedef WilsonCloverFermion<WilsonImplR> WilsonCloverFermionR;
+template <typename WImpl> using WilsonClover = WilsonCloverFermion<WImpl, CloverHelpers<WImpl>>;
-typedef WilsonCloverFermion<WilsonImplF> WilsonCloverFermionF;
+template <typename WImpl> using WilsonExpClover = WilsonCloverFermion<WImpl, ExpCloverHelpers<WImpl>>;
-typedef WilsonCloverFermion<WilsonImplD> WilsonCloverFermionD;
 
-typedef WilsonCloverFermion<WilsonAdjImplR> WilsonCloverAdjFermionR;
+typedef WilsonClover<WilsonImplR> WilsonCloverFermionR;
-typedef WilsonCloverFermion<WilsonAdjImplF> WilsonCloverAdjFermionF;
+typedef WilsonClover<WilsonImplF> WilsonCloverFermionF;
-typedef WilsonCloverFermion<WilsonAdjImplD> WilsonCloverAdjFermionD;
+typedef WilsonClover<WilsonImplD> WilsonCloverFermionD;
 
-typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplR> WilsonCloverTwoIndexSymmetricFermionR;
+typedef WilsonExpClover<WilsonImplR> WilsonExpCloverFermionR;
-typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplF> WilsonCloverTwoIndexSymmetricFermionF;
+typedef WilsonExpClover<WilsonImplF> WilsonExpCloverFermionF;
-typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplD> WilsonCloverTwoIndexSymmetricFermionD;
+typedef WilsonExpClover<WilsonImplD> WilsonExpCloverFermionD;
 
-typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplR> WilsonCloverTwoIndexAntiSymmetricFermionR;
+typedef WilsonClover<WilsonAdjImplR> WilsonCloverAdjFermionR;
-typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplF> WilsonCloverTwoIndexAntiSymmetricFermionF;
+typedef WilsonClover<WilsonAdjImplF> WilsonCloverAdjFermionF;
-typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplD> WilsonCloverTwoIndexAntiSymmetricFermionD;
+typedef WilsonClover<WilsonAdjImplD> WilsonCloverAdjFermionD;
+
+typedef WilsonClover<WilsonTwoIndexSymmetricImplR> WilsonCloverTwoIndexSymmetricFermionR;
+typedef WilsonClover<WilsonTwoIndexSymmetricImplF> WilsonCloverTwoIndexSymmetricFermionF;
+typedef WilsonClover<WilsonTwoIndexSymmetricImplD> WilsonCloverTwoIndexSymmetricFermionD;
+
+typedef WilsonClover<WilsonTwoIndexAntiSymmetricImplR> WilsonCloverTwoIndexAntiSymmetricFermionR;
+typedef WilsonClover<WilsonTwoIndexAntiSymmetricImplF> WilsonCloverTwoIndexAntiSymmetricFermionF;
+typedef WilsonClover<WilsonTwoIndexAntiSymmetricImplD> WilsonCloverTwoIndexAntiSymmetricFermionD;
 
 // Compact Clover fermions
-typedef CompactWilsonCloverFermion<WilsonImplR> CompactWilsonCloverFermionR;
+template <typename WImpl> using CompactWilsonClover = CompactWilsonCloverFermion<WImpl, CompactCloverHelpers<WImpl>>;
-typedef CompactWilsonCloverFermion<WilsonImplF> CompactWilsonCloverFermionF;
+template <typename WImpl> using CompactWilsonExpClover = CompactWilsonCloverFermion<WImpl, CompactExpCloverHelpers<WImpl>>;
-typedef CompactWilsonCloverFermion<WilsonImplD> CompactWilsonCloverFermionD;
 
-typedef CompactWilsonCloverFermion<WilsonAdjImplR> CompactWilsonCloverAdjFermionR;
+typedef CompactWilsonClover<WilsonImplR> CompactWilsonCloverFermionR;
-typedef CompactWilsonCloverFermion<WilsonAdjImplF> CompactWilsonCloverAdjFermionF;
+typedef CompactWilsonClover<WilsonImplF> CompactWilsonCloverFermionF;
-typedef CompactWilsonCloverFermion<WilsonAdjImplD> CompactWilsonCloverAdjFermionD;
+typedef CompactWilsonClover<WilsonImplD> CompactWilsonCloverFermionD;
 
-typedef CompactWilsonCloverFermion<WilsonTwoIndexSymmetricImplR> CompactWilsonCloverTwoIndexSymmetricFermionR;
+typedef CompactWilsonExpClover<WilsonImplR> CompactWilsonExpCloverFermionR;
-typedef CompactWilsonCloverFermion<WilsonTwoIndexSymmetricImplF> CompactWilsonCloverTwoIndexSymmetricFermionF;
+typedef CompactWilsonExpClover<WilsonImplF> CompactWilsonExpCloverFermionF;
-typedef CompactWilsonCloverFermion<WilsonTwoIndexSymmetricImplD> CompactWilsonCloverTwoIndexSymmetricFermionD;
+typedef CompactWilsonExpClover<WilsonImplD> CompactWilsonExpCloverFermionD;
 
-typedef CompactWilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplR> CompactWilsonCloverTwoIndexAntiSymmetricFermionR;
+typedef CompactWilsonClover<WilsonAdjImplR> CompactWilsonCloverAdjFermionR;
-typedef CompactWilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplF> CompactWilsonCloverTwoIndexAntiSymmetricFermionF;
+typedef CompactWilsonClover<WilsonAdjImplF> CompactWilsonCloverAdjFermionF;
-typedef CompactWilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplD> CompactWilsonCloverTwoIndexAntiSymmetricFermionD;
+typedef CompactWilsonClover<WilsonAdjImplD> CompactWilsonCloverAdjFermionD;
+
+typedef CompactWilsonClover<WilsonTwoIndexSymmetricImplR> CompactWilsonCloverTwoIndexSymmetricFermionR;
+typedef CompactWilsonClover<WilsonTwoIndexSymmetricImplF> CompactWilsonCloverTwoIndexSymmetricFermionF;
+typedef CompactWilsonClover<WilsonTwoIndexSymmetricImplD> CompactWilsonCloverTwoIndexSymmetricFermionD;
+
+typedef CompactWilsonClover<WilsonTwoIndexAntiSymmetricImplR> CompactWilsonCloverTwoIndexAntiSymmetricFermionR;
+typedef CompactWilsonClover<WilsonTwoIndexAntiSymmetricImplF> CompactWilsonCloverTwoIndexAntiSymmetricFermionF;
+typedef CompactWilsonClover<WilsonTwoIndexAntiSymmetricImplD> CompactWilsonCloverTwoIndexAntiSymmetricFermionD;
 
 // Domain Wall fermions
 typedef DomainWallFermion<WilsonImplR> DomainWallFermionR;

Grid/qcd/action/fermion/WilsonCloverFermion.h

@@ -32,6 +32,7 @@
 
 #include <Grid/qcd/action/fermion/WilsonCloverTypes.h>
 #include <Grid/qcd/action/fermion/WilsonCloverHelpers.h>
+#include <Grid/qcd/action/fermion/CloverHelpers.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -51,7 +52,7 @@ NAMESPACE_BEGIN(Grid);
 // csw_r = csw_t to recover the isotropic version
 //////////////////////////////////////////////////////////////////
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
 class WilsonCloverFermion : public WilsonFermion<Impl>,
                             public WilsonCloverHelpers<Impl>
 {

Grid/qcd/action/fermion/WilsonCloverHelpers.h

@@ -209,6 +209,8 @@ public:
 };
 
 
+////////////////////////////////////////////////////////
+
 template<class Impl> class CompactWilsonCloverHelpers {
 public:
 

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

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

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

@@ -32,17 +32,18 @@
 #include <Grid/qcd/spin/Dirac.h>
 #include <Grid/qcd/action/fermion/CompactWilsonCloverFermion.h>
 
+
 NAMESPACE_BEGIN(Grid);
-template<class Impl>
+template<class Impl, class CloverHelpers>
-CompactWilsonCloverFermion<Impl>::CompactWilsonCloverFermion(GaugeField& _Umu,
+CompactWilsonCloverFermion<Impl, CloverHelpers>::CompactWilsonCloverFermion(GaugeField& _Umu,
-                                                             GridCartesian& Fgrid,
+                                                                            GridCartesian& Fgrid,
-                                                             GridRedBlackCartesian& Hgrid,
+                                                                            GridRedBlackCartesian& Hgrid,
-                                                             const RealD _mass,
+                                                                            const RealD _mass,
-                                                             const RealD _csw_r,
+                                                                            const RealD _csw_r,
-                                                             const RealD _csw_t,
+                                                                            const RealD _csw_t,
-                                                             const RealD _cF,
+                                                                            const RealD _cF,
-                                                             const WilsonAnisotropyCoefficients& clover_anisotropy,
+                                                                            const WilsonAnisotropyCoefficients& clover_anisotropy,
-                                                             const ImplParams& impl_p)
+                                                                            const ImplParams& impl_p)
   : WilsonBase(_Umu, Fgrid, Hgrid, _mass, impl_p, clover_anisotropy)
   , csw_r(_csw_r)
   , csw_t(_csw_t)
@@ -64,44 +65,47 @@ CompactWilsonCloverFermion<Impl>::CompactWilsonCloverFermion(GaugeField& _Umu,
     csw_r /= clover_anisotropy.xi_0;
 
   ImportGauge(_Umu);
-  if (open_boundaries)
+  if (open_boundaries) {
+    this->BoundaryMaskEven.Checkerboard() = Even;
+    this->BoundaryMaskOdd.Checkerboard() = Odd;
     CompactHelpers::SetupMasks(this->BoundaryMask, this->BoundaryMaskEven, this->BoundaryMaskOdd);
+  }
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::Dhop(const FermionField& in, FermionField& out, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::Dhop(const FermionField& in, FermionField& out, int dag) {
   WilsonBase::Dhop(in, out, dag);
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::DhopOE(const FermionField& in, FermionField& out, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopOE(const FermionField& in, FermionField& out, int dag) {
   WilsonBase::DhopOE(in, out, dag);
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::DhopEO(const FermionField& in, FermionField& out, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopEO(const FermionField& in, FermionField& out, int dag) {
   WilsonBase::DhopEO(in, out, dag);
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::DhopDir(const FermionField& in, FermionField& out, int dir, int disp) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopDir(const FermionField& in, FermionField& out, int dir, int disp) {
   WilsonBase::DhopDir(in, out, dir, disp);
   if(this->open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::DhopDirAll(const FermionField& in, std::vector<FermionField>& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopDirAll(const FermionField& in, std::vector<FermionField>& out) {
   WilsonBase::DhopDirAll(in, out);
   if(this->open_boundaries) {
     for(auto& o : out) ApplyBoundaryMask(o);
   }
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::M(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::M(const FermionField& in, FermionField& out) {
   out.Checkerboard() = in.Checkerboard();
   WilsonBase::Dhop(in, out, DaggerNo); // call base to save applying bc
   Mooee(in, Tmp);
@@ -109,8 +113,8 @@ void CompactWilsonCloverFermion<Impl>::M(const FermionField& in, FermionField& o
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::Mdag(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::Mdag(const FermionField& in, FermionField& out) {
   out.Checkerboard() = in.Checkerboard();
   WilsonBase::Dhop(in, out, DaggerYes);  // call base to save applying bc
   MooeeDag(in, Tmp);
@@ -118,20 +122,20 @@ void CompactWilsonCloverFermion<Impl>::Mdag(const FermionField& in, FermionField
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::Meooe(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::Meooe(const FermionField& in, FermionField& out) {
   WilsonBase::Meooe(in, out);
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MeooeDag(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MeooeDag(const FermionField& in, FermionField& out) {
   WilsonBase::MeooeDag(in, out);
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::Mooee(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::Mooee(const FermionField& in, FermionField& out) {
   if(in.Grid()->_isCheckerBoarded) {
     if(in.Checkerboard() == Odd) {
       MooeeInternal(in, out, DiagonalOdd, TriangleOdd);
@@ -144,13 +148,13 @@ void CompactWilsonCloverFermion<Impl>::Mooee(const FermionField& in, FermionFiel
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MooeeDag(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeDag(const FermionField& in, FermionField& out) {
   Mooee(in, out); // blocks are hermitian
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MooeeInv(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeInv(const FermionField& in, FermionField& out) {
   if(in.Grid()->_isCheckerBoarded) {
     if(in.Checkerboard() == Odd) {
       MooeeInternal(in, out, DiagonalInvOdd, TriangleInvOdd);
@@ -163,23 +167,23 @@ void CompactWilsonCloverFermion<Impl>::MooeeInv(const FermionField& in, FermionF
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MooeeInvDag(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeInvDag(const FermionField& in, FermionField& out) {
   MooeeInv(in, out); // blocks are hermitian
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::Mdir(const FermionField& in, FermionField& out, int dir, int disp) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::Mdir(const FermionField& in, FermionField& out, int dir, int disp) {
   DhopDir(in, out, dir, disp);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MdirAll(const FermionField& in, std::vector<FermionField>& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MdirAll(const FermionField& in, std::vector<FermionField>& out) {
   DhopDirAll(in, out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MDeriv(GaugeField& force, const FermionField& X, const FermionField& Y, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MDeriv(GaugeField& force, const FermionField& X, const FermionField& Y, int dag) {
   assert(!open_boundaries); // TODO check for changes required for open bc
 
   // NOTE: code copied from original clover term
@@ -251,7 +255,7 @@ void CompactWilsonCloverFermion<Impl>::MDeriv(GaugeField& force, const FermionFi
       }
       PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
       Impl::TraceSpinImpl(lambda, Slambda);                   // traceSpin ok
-      force_mu -= factor*Helpers::Cmunu(U, lambda, mu, nu);   // checked
+      force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu);   // checked
       count++;
     }
 
@@ -261,18 +265,18 @@ void CompactWilsonCloverFermion<Impl>::MDeriv(GaugeField& force, const FermionFi
   force += clover_force;
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MooDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
   assert(0);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MeeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MeeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
   assert(0);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MooeeInternal(const FermionField&        in,
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeInternal(const FermionField&        in,
                     FermionField&              out,
                     const CloverDiagonalField& diagonal,
                     const CloverTriangleField& triangle) {
@@ -285,8 +289,8 @@ void CompactWilsonCloverFermion<Impl>::MooeeInternal(const FermionField&
   CompactHelpers::MooeeKernel(diagonal.oSites(), 1, in, out, diagonal, triangle);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::ImportGauge(const GaugeField& _Umu) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::ImportGauge(const GaugeField& _Umu) {
   // NOTE: parts copied from original implementation
 
   // Import gauge into base class
@@ -318,22 +322,27 @@ void CompactWilsonCloverFermion<Impl>::ImportGauge(const GaugeField& _Umu) {
   TmpOriginal += Helpers::fillCloverXT(Ex) * csw_t;
   TmpOriginal += Helpers::fillCloverYT(Ey) * csw_t;
   TmpOriginal += Helpers::fillCloverZT(Ez) * csw_t;
-  TmpOriginal += this->diag_mass;
+  // Handle mass term based on clover policy
-
+  CloverHelpers::MassTerm(TmpOriginal, this->diag_mass);
+  
   // Convert the data layout of the clover term
   double t4 = usecond();
   CompactHelpers::ConvertLayout(TmpOriginal, Diagonal, Triangle);
 
-  // Possible modify the boundary values
+  // Exponentiate the clover (nothing happens in case of the standard clover)
   double t5 = usecond();
+  CloverHelpers::Exponentiate_Clover(Diagonal, Triangle, csw_t, this->diag_mass);
+
+  // Possible modify the boundary values
+  double t6 = usecond();
   if(open_boundaries) CompactHelpers::ModifyBoundaries(Diagonal, Triangle, csw_t, cF, this->diag_mass);
 
-  // Invert the clover term in the improved layout
+  // Invert the Clover term (explicit inversion needed for the improvement in case of open boundary conditions)
-  double t6 = usecond();
+  double t7 = usecond();
   CompactHelpers::Invert(Diagonal, Triangle, DiagonalInv, TriangleInv);
 
   // Fill the remaining clover fields
-  double t7 = usecond();
+  double t8 = usecond();
   pickCheckerboard(Even, DiagonalEven,    Diagonal);
   pickCheckerboard(Even, TriangleEven,    Triangle);
   pickCheckerboard(Odd,  DiagonalOdd,     Diagonal);
@@ -344,20 +353,19 @@ void CompactWilsonCloverFermion<Impl>::ImportGauge(const GaugeField& _Umu) {
   pickCheckerboard(Odd,  TriangleInvOdd,  TriangleInv);
 
   // Report timings
-  double t8 = usecond();
+  double t9 = usecond();
-#if 0
+
-  std::cout << GridLogMessage << "CompactWilsonCloverFermion::ImportGauge timings:"
+  std::cout << GridLogDebug << "CompactWilsonCloverFermion::ImportGauge timings:" << std::endl;
-            << " WilsonFermion::Importgauge = " << (t1 - t0) / 1e6
+  std::cout << GridLogDebug << "WilsonFermion::Importgauge = " << (t1 - t0) / 1e6 << std::endl;
-            << ", allocations = "               << (t2 - t1) / 1e6
+  std::cout << GridLogDebug << "allocations =                " << (t2 - t1) / 1e6 << std::endl;
-            << ", field strength = "            << (t3 - t2) / 1e6
+  std::cout << GridLogDebug << "field strength =             " << (t3 - t2) / 1e6 << std::endl;
-            << ", fill clover = "               << (t4 - t3) / 1e6
+  std::cout << GridLogDebug << "fill clover =                " << (t4 - t3) / 1e6 << std::endl;
-            << ", convert = "                   << (t5 - t4) / 1e6
+  std::cout << GridLogDebug << "convert =                    " << (t5 - t4) / 1e6 << std::endl;
-            << ", boundaries = "                << (t6 - t5) / 1e6
+  std::cout << GridLogDebug << "exponentiation =             " << (t6 - t5) / 1e6 << std::endl;
-            << ", inversions = "                << (t7 - t6) / 1e6
+  std::cout << GridLogDebug << "boundaries =                 " << (t7 - t6) / 1e6 << std::endl;
-            << ", pick cbs = "                  << (t8 - t7) / 1e6
+  std::cout << GridLogDebug << "inversions =                 " << (t8 - t7) / 1e6 << std::endl;
-            << ", total = "                     << (t8 - t0) / 1e6
+  std::cout << GridLogDebug << "pick cbs =                   " << (t9 - t8) / 1e6 << std::endl;
-            << std::endl;
+  std::cout << GridLogDebug << "total =                      " << (t9 - t0) / 1e6 << std::endl;
-#endif
 }
 
 NAMESPACE_END(Grid);

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

@@ -34,8 +34,8 @@
 
 NAMESPACE_BEGIN(Grid);
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-WilsonCloverFermion<Impl>::WilsonCloverFermion(GaugeField&                         _Umu,
+WilsonCloverFermion<Impl, CloverHelpers>::WilsonCloverFermion(GaugeField&                         _Umu,
                                                GridCartesian&                      Fgrid,
                                                GridRedBlackCartesian&              Hgrid,
                                                const RealD                         _mass,
@@ -74,8 +74,8 @@ WilsonCloverFermion<Impl>::WilsonCloverFermion(GaugeField&
 }
 
 // *NOT* EO
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::M(const FermionField &in, FermionField &out)
 {
   FermionField temp(out.Grid());
 
@@ -89,8 +89,8 @@ void WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
   out += temp;
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::Mdag(const FermionField &in, FermionField &out)
 {
   FermionField temp(out.Grid());
 
@@ -104,8 +104,8 @@ void WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
   out += temp;
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
+void WilsonCloverFermion<Impl, CloverHelpers>::ImportGauge(const GaugeField &_Umu)
 {
   double t0 = usecond();
   WilsonFermion<Impl>::ImportGauge(_Umu);
@@ -131,47 +131,11 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
   CloverTerm += Helpers::fillCloverXT(Ex) * csw_t;
   CloverTerm += Helpers::fillCloverYT(Ey) * csw_t;
   CloverTerm += Helpers::fillCloverZT(Ez) * csw_t;
-  CloverTerm += diag_mass;
+   
-
   double t4 = usecond();
-  int lvol = _Umu.Grid()->lSites();
+  CloverHelpers::Instantiate(CloverTerm, CloverTermInv, csw_t, this->diag_mass);
-  int DimRep = Impl::Dimension;
 
   double t5 = usecond();
-  {
-    autoView(CTv,CloverTerm,CpuRead);
-    autoView(CTIv,CloverTermInv,CpuWrite);
-    thread_for(site, lvol, {
-      Coordinate lcoor;
-      grid->LocalIndexToLocalCoor(site, lcoor);
-      Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
-      Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
-      typename SiteClover::scalar_object Qx = Zero(), Qxinv = Zero();
-      peekLocalSite(Qx, CTv, lcoor);
-      //if (csw!=0){
-      for (int j = 0; j < Ns; j++)
-	for (int k = 0; k < Ns; k++)
-	  for (int a = 0; a < DimRep; a++)
-	    for (int b = 0; b < DimRep; b++){
-	      auto zz =  Qx()(j, k)(a, b);
-	      EigenCloverOp(a + j * DimRep, b + k * DimRep) = std::complex<double>(zz);
-	    }
-      //   if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl;
-      
-      EigenInvCloverOp = EigenCloverOp.inverse();
-      //std::cout << EigenInvCloverOp << std::endl;
-      for (int j = 0; j < Ns; j++)
-	for (int k = 0; k < Ns; k++)
-	  for (int a = 0; a < DimRep; a++)
-	    for (int b = 0; b < DimRep; b++)
-	      Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep);
-      //    if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl;
-      //  }
-      pokeLocalSite(Qxinv, CTIv, lcoor);
-    });
-  }
-
-  double t6 = usecond();
   // Separate the even and odd parts
   pickCheckerboard(Even, CloverTermEven, CloverTerm);
   pickCheckerboard(Odd, CloverTermOdd, CloverTerm);
@@ -184,48 +148,44 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
 
   pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv));
   pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv));
-  double t7 = usecond();
+  double t6 = usecond();
 
-#if 0
+  std::cout << GridLogDebug << "WilsonCloverFermion::ImportGauge timings:" << std::endl;
-  std::cout << GridLogMessage << "WilsonCloverFermion::ImportGauge timings:"
+  std::cout << GridLogDebug << "WilsonFermion::Importgauge = " << (t1 - t0) / 1e6 << std::endl;
-            << " WilsonFermion::Importgauge = " << (t1 - t0) / 1e6
+  std::cout << GridLogDebug << "allocations =                " << (t2 - t1) / 1e6 << std::endl;
-            << ", allocations = "               << (t2 - t1) / 1e6
+  std::cout << GridLogDebug << "field strength =             " << (t3 - t2) / 1e6 << std::endl;
-            << ", field strength = "            << (t3 - t2) / 1e6
+  std::cout << GridLogDebug << "fill clover =                " << (t4 - t3) / 1e6 << std::endl;
-            << ", fill clover = "               << (t4 - t3) / 1e6
+  std::cout << GridLogDebug << "instantiation =              " << (t5 - t4) / 1e6 << std::endl;
-            << ", misc = "                      << (t5 - t4) / 1e6
+  std::cout << GridLogDebug << "pick cbs =                   " << (t6 - t5) / 1e6 << std::endl;
-            << ", inversions = "                << (t6 - t5) / 1e6
+  std::cout << GridLogDebug << "total =                      " << (t6 - t0) / 1e6 << std::endl;
-            << ", pick cbs = "                  << (t7 - t6) / 1e6
-            << ", total = "                     << (t7 - t0) / 1e6
-            << std::endl;
-#endif
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::Mooee(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::Mooee(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerNo, InverseNo);
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::MooeeDag(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerYes, InverseNo);
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::MooeeInv(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerNo, InverseYes);
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::MooeeInvDag(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerYes, InverseYes);
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv)
+void WilsonCloverFermion<Impl, CloverHelpers>::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv)
 {
   out.Checkerboard() = in.Checkerboard();
   CloverField *Clover;
@@ -278,8 +238,8 @@ void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionFie
 } // MooeeInternal
 
 // Derivative parts unpreconditioned pseudofermions
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag)
+void WilsonCloverFermion<Impl, CloverHelpers>::MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag)
 {
   conformable(X.Grid(), Y.Grid());
   conformable(X.Grid(), force.Grid());
@@ -349,7 +309,7 @@ void WilsonCloverFermion<Impl>::MDeriv(GaugeField &force, const FermionField &X,
       }
       PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
       Impl::TraceSpinImpl(lambda, Slambda);                   // traceSpin ok
-      force_mu -= factor*Helpers::Cmunu(U, lambda, mu, nu);                   // checked
+      force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu);                   // checked
       count++;
     }
 
@@ -360,15 +320,15 @@ void WilsonCloverFermion<Impl>::MDeriv(GaugeField &force, const FermionField &X,
 }
 
 // Derivative parts
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag)
+void WilsonCloverFermion<Impl, CloverHelpers>::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag)
 {
   assert(0);
 }
 
 // Derivative parts
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
+void WilsonCloverFermion<Impl, CloverHelpers>::MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
 {
   assert(0); // not implemented yet
 }

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

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

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

@@ -9,6 +9,7 @@
     Author: paboyle <paboyle@ph.ed.ac.uk>
     Author: Guido Cossu <guido.cossu@ed.ac.uk>
     Author: Daniel Richtmann <daniel.richtmann@gmail.com>
+    Author: Mattia Bruno <mattia.bruno@cern.ch>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -32,10 +33,12 @@
 #include <Grid/qcd/spin/Dirac.h>
 #include <Grid/qcd/action/fermion/CompactWilsonCloverFermion.h>
 #include <Grid/qcd/action/fermion/implementation/CompactWilsonCloverFermionImplementation.h>
+#include <Grid/qcd/action/fermion/CloverHelpers.h>
 
 NAMESPACE_BEGIN(Grid);
 
 #include "impl.h"
-template class CompactWilsonCloverFermion<IMPLEMENTATION>; 
+template class CompactWilsonCloverFermion<IMPLEMENTATION, CompactCloverHelpers<IMPLEMENTATION>>; 
+template class CompactWilsonCloverFermion<IMPLEMENTATION, CompactExpCloverHelpers<IMPLEMENTATION>>; 
 
 NAMESPACE_END(Grid);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Grid/qcd/utils/GaugeFix.h

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

Grid/qcd/utils/WilsonLoops.h

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

benchmarks/Benchmark_mooee.cc

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

configure.ac

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

tests/IO/Test_nersc_io.cc

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

tests/core/Test_compact_wilson_clover_speedup.cc

@@ -117,8 +117,8 @@ void runBenchmark(int* argc, char*** argv) {
 
   // type definitions
   typedef WilsonImpl<vCoeff_t, FundamentalRepresentation, CoeffReal> WImpl;
-  typedef WilsonCloverFermion<WImpl> WilsonCloverOperator;
+  typedef WilsonCloverFermion<WImpl, CloverHelpers<WImpl>> WilsonCloverOperator;
-  typedef CompactWilsonCloverFermion<WImpl> CompactWilsonCloverOperator;
+  typedef CompactWilsonCloverFermion<WImpl, CompactCloverHelpers<WImpl>> CompactWilsonCloverOperator;
   typedef typename WilsonCloverOperator::FermionField Fermion;
   typedef typename WilsonCloverOperator::GaugeField Gauge;
 

tests/core/Test_lie_generators.cc

@@ -9,6 +9,7 @@ Copyright (C) 2015
 Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: Guido Cossu <guido.cossu@ed.ac.uk>
+Author: Jamie Hudspith <renwick.james.hudspth@gmail.com>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -42,14 +43,14 @@ directory
 
 using namespace std;
 using namespace Grid;
- ;
+;
 
 int main(int argc, char** argv) {
   Grid_init(&argc, &argv);
 
   std::vector<int> latt({4, 4, 4, 8});
   GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid(
-      latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
+						       latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
 
   GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid);
 
@@ -60,15 +61,19 @@ int main(int argc, char** argv) {
             << std::endl;
   SU2::printGenerators();
   std::cout << "Dimension of adjoint representation: "<< SU2Adjoint::Dimension << std::endl;
+
+  // guard as this code fails to compile for Nc != 3
+#if (Nc == 3)
+    
   SU2Adjoint::printGenerators();
   SU2::testGenerators();
   SU2Adjoint::testGenerators();
-
+    
   std::cout << GridLogMessage << "*********************************************"
-            << std::endl;
+	    << std::endl;
   std::cout << GridLogMessage << "* Generators for SU(Nc" << std::endl;
   std::cout << GridLogMessage << "*********************************************"
-            << std::endl;
+	    << std::endl;
   SU3::printGenerators();
   std::cout << "Dimension of adjoint representation: "<< SU3Adjoint::Dimension << std::endl;
   SU3Adjoint::printGenerators();
@@ -111,12 +116,10 @@ int main(int argc, char** argv) {
 
   // AdjointRepresentation has the predefined number of colours Nc
   //  Representations<FundamentalRepresentation, AdjointRepresentation, TwoIndexSymmetricRepresentation> RepresentationTypes(grid);  
-
-  
   LatticeGaugeField U(grid), V(grid);
   SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U);
   SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V);
-
+    
   // Adjoint representation
   // Test group structure
   // (U_f * V_f)_r = U_r * V_r
@@ -127,17 +130,17 @@ int main(int argc, char** argv) {
     SU3::LatticeMatrix Vmu = peekLorentz(V,mu);
     pokeLorentz(UV,Umu*Vmu, mu);
   }
-
+    
   AdjRep.update_representation(UV);
   typename AdjointRep<Nc>::LatticeField UVr = AdjRep.U;  // (U_f * V_f)_r 
-
+    
-
+    
   AdjRep.update_representation(U);
   typename AdjointRep<Nc>::LatticeField Ur = AdjRep.U;  // U_r
-
+    
   AdjRep.update_representation(V);
   typename AdjointRep<Nc>::LatticeField Vr = AdjRep.U;  // V_r
-
+    
   typename AdjointRep<Nc>::LatticeField UrVr(grid);
   UrVr = Zero();
   for (int mu = 0; mu < Nd; mu++) {
@@ -145,10 +148,10 @@ int main(int argc, char** argv) {
     typename AdjointRep<Nc>::LatticeMatrix Vrmu = peekLorentz(Vr,mu);
     pokeLorentz(UrVr,Urmu*Vrmu, mu);
   }
-
+    
   typename AdjointRep<Nc>::LatticeField Diff_check = UVr - UrVr;
   std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Adjoint representation) : " << norm2(Diff_check) << std::endl;
-
+    
   // Check correspondence of algebra and group transformations
   // Create a random vector
   SU3::LatticeAlgebraVector h_adj(grid);
@@ -156,32 +159,31 @@ int main(int argc, char** argv) {
   random(gridRNG,h_adj);
   h_adj = real(h_adj);
   SU_Adjoint<Nc>::AdjointLieAlgebraMatrix(h_adj,Ar);
-
+    
   // Re-extract h_adj
   SU3::LatticeAlgebraVector h_adj2(grid);
   SU_Adjoint<Nc>::projectOnAlgebra(h_adj2, Ar);
   SU3::LatticeAlgebraVector h_diff = h_adj - h_adj2;
   std::cout << GridLogMessage << "Projections structure check vector difference (Adjoint representation) : " << norm2(h_diff) << std::endl;
-
+    
   // Exponentiate
   typename AdjointRep<Nc>::LatticeMatrix Uadj(grid);
   Uadj  = expMat(Ar, 1.0, 16);
-
+    
   typename AdjointRep<Nc>::LatticeMatrix uno(grid);
   uno = 1.0;
   // Check matrix Uadj, must be real orthogonal
   typename AdjointRep<Nc>::LatticeMatrix Ucheck = Uadj - conjugate(Uadj);
   std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck)
-            << std::endl;
+	    << std::endl;
-
+    
   Ucheck = Uadj * adj(Uadj) - uno;
   std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck)
-            << std::endl;
+	    << std::endl;
   Ucheck = adj(Uadj) * Uadj - uno;
   std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck)
-            << std::endl;
+	    << std::endl;
-      
+    
-
   // Construct the fundamental matrix in the group
   SU3::LatticeMatrix Af(grid);
   SU3::FundamentalLieAlgebraMatrix(h_adj,Af);
@@ -193,72 +195,65 @@ int main(int argc, char** argv) {
   SU3::LatticeMatrix UnitCheck(grid);
   UnitCheck = Ufund * adj(Ufund) - uno_f;
   std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck)
-            << std::endl;
+	    << std::endl;
   UnitCheck = adj(Ufund) * Ufund - uno_f;
   std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck)
-            << std::endl;
+	    << std::endl;
-
 
   // Tranform to the adjoint representation
   U = Zero(); // fill this with only one direction
   pokeLorentz(U,Ufund,0); // the representation transf acts on full gauge fields
-
+    
   AdjRep.update_representation(U);
   Ur = AdjRep.U;  // U_r  
   typename AdjointRep<Nc>::LatticeMatrix Ur0 = peekLorentz(Ur,0); // this should be the same as Uadj
-
+    
   typename AdjointRep<Nc>::LatticeMatrix Diff_check_mat = Ur0 - Uadj;
   std::cout << GridLogMessage << "Projections structure check group difference : " << norm2(Diff_check_mat) << std::endl;
 
-
-
-
   // TwoIndexRep tests
-
   std::cout << GridLogMessage << "*********************************************"
-            << std::endl;
+	    << std::endl;
   std::cout << GridLogMessage << "*********************************************"
-            << std::endl;
+	    << std::endl;
 
 
   
   std::cout << GridLogMessage << "* eS^{ij} base for SU(2)" << std::endl;
   std::cout << GridLogMessage << "*********************************************"
-            << std::endl;
+	    << std::endl;
   std::cout << GridLogMessage << "Dimension of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl;
   SU2TwoIndexSymm::printBase();
-      std::cout << GridLogMessage << "*********************************************"
+  std::cout << GridLogMessage << "*********************************************"
-            << std::endl;
+	    << std::endl;
-        std::cout << GridLogMessage << "Generators of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl;
+  std::cout << GridLogMessage << "Generators of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl;
   SU2TwoIndexSymm::printGenerators();
-        std::cout << GridLogMessage << "Test of Two Index Symmetric Generators: "<< SU2TwoIndexSymm::Dimension << std::endl;
+  std::cout << GridLogMessage << "Test of Two Index Symmetric Generators: "<< SU2TwoIndexSymm::Dimension << std::endl;
   SU2TwoIndexSymm::testGenerators();
   std::cout << GridLogMessage << "*********************************************"
-            << std::endl;
+	    << std::endl;
 
 
   
   std::cout << GridLogMessage << "*********************************************"
-            << std::endl;
+	    << std::endl;
   std::cout << GridLogMessage << "* eAS^{ij} base for SU(2)" << std::endl;
   
   std::cout << GridLogMessage << "*********************************************"
-            << std::endl;
+	    << std::endl;
   std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl;
   SU2TwoIndexAntiSymm::printBase();
-      std::cout << GridLogMessage << "*********************************************"
+  std::cout << GridLogMessage << "*********************************************"
-            << std::endl;
+	    << std::endl;
-        std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl;
+  std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl;
   SU2TwoIndexAntiSymm::printGenerators();
   std::cout << GridLogMessage << "Test of Two Index anti-Symmetric Generators: "<< SU2TwoIndexAntiSymm::Dimension << std::endl;
   SU2TwoIndexAntiSymm::testGenerators();
   
-  
-  
   std::cout << GridLogMessage << "*********************************************"
-      << std::endl;
+	    << std::endl;
   std::cout << GridLogMessage << "Test for the Two Index Symmetric projectors"
-      << std::endl;
+	    << std::endl;
   // Projectors 
   SU3TwoIndexSymm::LatticeTwoIndexMatrix Gauss2(grid);
   random(gridRNG,Gauss2);
@@ -276,13 +271,13 @@ int main(int argc, char** argv) {
   SU3::LatticeAlgebraVector diff2 = ha - hb;
   std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl;
   std::cout << GridLogMessage << "*********************************************"
-      << std::endl;
+	    << std::endl;
 
   
-    std::cout << GridLogMessage << "*********************************************"
+  std::cout << GridLogMessage << "*********************************************"
-      << std::endl;
+	    << std::endl;
   std::cout << GridLogMessage << "Test for the Two index anti-Symmetric projectors"
-      << std::endl;
+	    << std::endl;
   // Projectors
   SU3TwoIndexAntiSymm::LatticeTwoIndexMatrix Gauss2a(grid);
   random(gridRNG,Gauss2a);
@@ -300,11 +295,11 @@ int main(int argc, char** argv) {
   SU3::LatticeAlgebraVector diff2a = ha - hb;
   std::cout << GridLogMessage << "Difference: " << norm2(diff2a) << std::endl;
   std::cout << GridLogMessage << "*********************************************"
-      << std::endl;
+	    << std::endl;
     
   
   std::cout << GridLogMessage << "Two index Symmetric: Checking Group Structure"
-      << std::endl;
+	    << std::endl;
   // Testing HMC representation classes
   TwoIndexRep< Nc, Symmetric > TIndexRep(grid);
 
@@ -313,7 +308,7 @@ int main(int argc, char** argv) {
   LatticeGaugeField U2(grid), V2(grid);
   SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2);
   SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2);
-  
+    
   LatticeGaugeField UV2(grid);
   UV2 = Zero();
   for (int mu = 0; mu < Nd; mu++) {
@@ -321,16 +316,16 @@ int main(int argc, char** argv) {
     SU3::LatticeMatrix Vmu2 = peekLorentz(V2,mu);
     pokeLorentz(UV2,Umu2*Vmu2, mu);
   }
-  
+    
   TIndexRep.update_representation(UV2);
   typename TwoIndexRep< Nc, Symmetric >::LatticeField UVr2 = TIndexRep.U;  // (U_f * V_f)_r
-  
+    
   TIndexRep.update_representation(U2);
   typename TwoIndexRep< Nc, Symmetric >::LatticeField Ur2 = TIndexRep.U;  // U_r
-  
+    
   TIndexRep.update_representation(V2);
   typename TwoIndexRep< Nc, Symmetric >::LatticeField Vr2 = TIndexRep.U;  // V_r
-  
+    
   typename TwoIndexRep< Nc, Symmetric >::LatticeField Ur2Vr2(grid);
   Ur2Vr2 = Zero();
   for (int mu = 0; mu < Nd; mu++) {
@@ -338,11 +333,11 @@ int main(int argc, char** argv) {
     typename TwoIndexRep< Nc, Symmetric >::LatticeMatrix Vrmu2 = peekLorentz(Vr2,mu);
     pokeLorentz(Ur2Vr2,Urmu2*Vrmu2, mu);
   }
-  
+    
   typename TwoIndexRep< Nc, Symmetric >::LatticeField Diff_check2 = UVr2 - Ur2Vr2;
   std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Two Index Symmetric): " << norm2(Diff_check2) << std::endl;
-
+    
-  
+    
   // Check correspondence of algebra and group transformations
   // Create a random vector
   SU3::LatticeAlgebraVector h_sym(grid);
@@ -350,34 +345,31 @@ int main(int argc, char** argv) {
   random(gridRNG,h_sym);
   h_sym = real(h_sym);
   SU_TwoIndex<Nc,Symmetric>::TwoIndexLieAlgebraMatrix(h_sym,Ar_sym);
-  
+    
   // Re-extract h_sym
   SU3::LatticeAlgebraVector h_sym2(grid);
   SU_TwoIndex< Nc, Symmetric>::projectOnAlgebra(h_sym2, Ar_sym);
   SU3::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2;
   std::cout << GridLogMessage << "Projections structure check vector difference (Two Index Symmetric): " << norm2(h_diff_sym) << std::endl;
-
+    
-  
   // Exponentiate
   typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix U2iS(grid);
   U2iS  = expMat(Ar_sym, 1.0, 16);
-  
+    
   typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix uno2iS(grid);
   uno2iS = 1.0;
   // Check matrix U2iS, must be real orthogonal
   typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ucheck2iS = U2iS - conjugate(U2iS);
   std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iS)
-      << std::endl;
+	    << std::endl;
-  
+    
   Ucheck2iS = U2iS * adj(U2iS) - uno2iS;
   std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iS)
-      << std::endl;
+	    << std::endl;
   Ucheck2iS = adj(U2iS) * U2iS - uno2iS;
   std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iS)
-      << std::endl;
+	    << std::endl;
-  
+    
-  
-  
   // Construct the fundamental matrix in the group
   SU3::LatticeMatrix Af_sym(grid);
   SU3::FundamentalLieAlgebraMatrix(h_sym,Af_sym);
@@ -386,147 +378,137 @@ int main(int argc, char** argv) {
   SU3::LatticeMatrix UnitCheck2(grid);
   UnitCheck2 = Ufund2 * adj(Ufund2) - uno_f;
   std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2)
-      << std::endl;
+	    << std::endl;
   UnitCheck2 = adj(Ufund2) * Ufund2 - uno_f;
   std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2)
-      << std::endl;
+	    << std::endl;
-  
+    
-
+    
   // Tranform to the 2Index Sym representation
   U = Zero(); // fill this with only one direction
   pokeLorentz(U,Ufund2,0); // the representation transf acts on full gauge fields
-  
+    
   TIndexRep.update_representation(U);
   Ur2 = TIndexRep.U;  // U_r  
   typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ur02 = peekLorentz(Ur2,0); // this should be the same as U2iS
-  
+    
   typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Diff_check_mat2 = Ur02 - U2iS;
   std::cout << GridLogMessage << "Projections structure check group difference (Two Index Symmetric): " << norm2(Diff_check_mat2) << std::endl;
-  
+    
-
-
-
   if (TwoIndexRep<Nc, AntiSymmetric >::Dimension != 1){
 
-  std::cout << GridLogMessage << "*********************************************"
+    std::cout << GridLogMessage << "*********************************************"
-      << std::endl;
+	      << std::endl;
     
   
-  std::cout << GridLogMessage << "Two Index anti-Symmetric: Check Group Structure"
+    std::cout << GridLogMessage << "Two Index anti-Symmetric: Check Group Structure"
-      << std::endl;
+	      << std::endl;
-  // Testing HMC representation classes
+    // Testing HMC representation classes
-  TwoIndexRep< Nc, AntiSymmetric > TIndexRepA(grid);
+    TwoIndexRep< Nc, AntiSymmetric > TIndexRepA(grid);
 
 
-  // Test group structure
+    // Test group structure
-  // (U_f * V_f)_r = U_r * V_r
+    // (U_f * V_f)_r = U_r * V_r
-  LatticeGaugeField U2A(grid), V2A(grid);
+    LatticeGaugeField U2A(grid), V2A(grid);
-  SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2A);
+    SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2A);
-  SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2A);
+    SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2A);
   
-  LatticeGaugeField UV2A(grid);
+    LatticeGaugeField UV2A(grid);
-  UV2A = Zero();
+    UV2A = Zero();
-  for (int mu = 0; mu < Nd; mu++) {
+    for (int mu = 0; mu < Nd; mu++) {
-    SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu);
+      SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu);
-    SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu);
+      SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu);
-    pokeLorentz(UV2A,Umu2A*Vmu2A, mu);
+      pokeLorentz(UV2A,Umu2A*Vmu2A, mu);
+    }
+  
+    TIndexRep.update_representation(UV2A);
+    typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField UVr2A = TIndexRepA.U;  // (U_f * V_f)_r
+  
+    TIndexRep.update_representation(U2A);
+    typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2A = TIndexRepA.U;  // U_r
+  
+    TIndexRep.update_representation(V2A);
+    typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Vr2A = TIndexRepA.U;  // V_r
+  
+    typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2Vr2A(grid);
+    Ur2Vr2A = Zero();
+    for (int mu = 0; mu < Nd; mu++) {
+      typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Urmu2A = peekLorentz(Ur2A,mu);
+      typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Vrmu2A = peekLorentz(Vr2A,mu);
+      pokeLorentz(Ur2Vr2A,Urmu2A*Vrmu2A, mu);
+    }
+  
+    typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Diff_check2A = UVr2A - Ur2Vr2A;
+    std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Two Index anti-Symmetric): " << norm2(Diff_check2A) << std::endl;
+
+  
+    // Check correspondence of algebra and group transformations
+    // Create a random vector
+    SU3::LatticeAlgebraVector h_Asym(grid);
+    typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ar_Asym(grid);
+    random(gridRNG,h_Asym);
+    h_Asym = real(h_Asym);
+    SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym);
+  
+    // Re-extract h_sym
+    SU3::LatticeAlgebraVector h_Asym2(grid);
+    SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym);
+    SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2;
+    std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl;
+
+  
+    // Exponentiate
+    typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix U2iAS(grid);
+    U2iAS  = expMat(Ar_Asym, 1.0, 16);
+  
+    typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix uno2iAS(grid);
+    uno2iAS = 1.0;
+    // Check matrix U2iS, must be real orthogonal
+    typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ucheck2iAS = U2iAS - conjugate(U2iAS);
+    std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iAS)
+	      << std::endl;
+  
+    Ucheck2iAS = U2iAS * adj(U2iAS) - uno2iAS;
+    std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iAS)
+	      << std::endl;
+    Ucheck2iAS = adj(U2iAS) * U2iAS - uno2iAS;
+    std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iAS)
+	      << std::endl;
+  
+  
+  
+    // Construct the fundamental matrix in the group
+    SU3::LatticeMatrix Af_Asym(grid);
+    SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym);
+    SU3::LatticeMatrix Ufund2A(grid);
+    Ufund2A  = expMat(Af_Asym, 1.0, 16);
+    SU3::LatticeMatrix UnitCheck2A(grid);
+    UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f;
+    std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A)
+	      << std::endl;
+    UnitCheck2A = adj(Ufund2A) * Ufund2A - uno_f;
+    std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2A)
+	      << std::endl;
+  
+
+    // Tranform to the 2Index Sym representation
+    U = Zero(); // fill this with only one direction
+    pokeLorentz(U,Ufund2A,0); // the representation transf acts on full gauge fields
+  
+    TIndexRepA.update_representation(U);
+    Ur2A = TIndexRepA.U;  // U_r  
+    typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ur02A = peekLorentz(Ur2A,0); // this should be the same as U2iS
+  
+    typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Diff_check_mat2A = Ur02A - U2iAS;
+    std::cout << GridLogMessage << "Projections structure check group difference (Two Index anti-Symmetric): " << norm2(Diff_check_mat2A) << std::endl;
+  
+  } else  {
+    std::cout << GridLogMessage << "Skipping Two Index anti-Symmetric tests "
+      "because representation is trivial (dim = 1)"
+	      << std::endl;
   }
-  
-  TIndexRep.update_representation(UV2A);
-  typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField UVr2A = TIndexRepA.U;  // (U_f * V_f)_r
-  
-  TIndexRep.update_representation(U2A);
-  typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2A = TIndexRepA.U;  // U_r
-  
-  TIndexRep.update_representation(V2A);
-  typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Vr2A = TIndexRepA.U;  // V_r
-  
-  typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2Vr2A(grid);
-  Ur2Vr2A = Zero();
-  for (int mu = 0; mu < Nd; mu++) {
-    typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Urmu2A = peekLorentz(Ur2A,mu);
-    typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Vrmu2A = peekLorentz(Vr2A,mu);
-    pokeLorentz(Ur2Vr2A,Urmu2A*Vrmu2A, mu);
-  }
-  
-  typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Diff_check2A = UVr2A - Ur2Vr2A;
-  std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Two Index anti-Symmetric): " << norm2(Diff_check2A) << std::endl;
-
-  
-  // Check correspondence of algebra and group transformations
-  // Create a random vector
-  SU3::LatticeAlgebraVector h_Asym(grid);
-  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ar_Asym(grid);
-  random(gridRNG,h_Asym);
-  h_Asym = real(h_Asym);
-  SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym);
-  
-  // Re-extract h_sym
-  SU3::LatticeAlgebraVector h_Asym2(grid);
-  SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym);
-  SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2;
-  std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl;
-
-  
-  // Exponentiate
-  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix U2iAS(grid);
-  U2iAS  = expMat(Ar_Asym, 1.0, 16);
-  
-  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix uno2iAS(grid);
-  uno2iAS = 1.0;
-  // Check matrix U2iS, must be real orthogonal
-  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ucheck2iAS = U2iAS - conjugate(U2iAS);
-  std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iAS)
-      << std::endl;
-  
-  Ucheck2iAS = U2iAS * adj(U2iAS) - uno2iAS;
-  std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iAS)
-      << std::endl;
-  Ucheck2iAS = adj(U2iAS) * U2iAS - uno2iAS;
-  std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iAS)
-      << std::endl;
-  
-  
-  
-  // Construct the fundamental matrix in the group
-  SU3::LatticeMatrix Af_Asym(grid);
-  SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym);
-  SU3::LatticeMatrix Ufund2A(grid);
-  Ufund2A  = expMat(Af_Asym, 1.0, 16);
-  SU3::LatticeMatrix UnitCheck2A(grid);
-  UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f;
-  std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A)
-      << std::endl;
-  UnitCheck2A = adj(Ufund2A) * Ufund2A - uno_f;
-  std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2A)
-      << std::endl;
-  
-
-  // Tranform to the 2Index Sym representation
-  U = Zero(); // fill this with only one direction
-  pokeLorentz(U,Ufund2A,0); // the representation transf acts on full gauge fields
-  
-  TIndexRepA.update_representation(U);
-  Ur2A = TIndexRepA.U;  // U_r  
-  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ur02A = peekLorentz(Ur2A,0); // this should be the same as U2iS
-  
-  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Diff_check_mat2A = Ur02A - U2iAS;
-  std::cout << GridLogMessage << "Projections structure check group difference (Two Index anti-Symmetric): " << norm2(Diff_check_mat2A) << std::endl;
-  
-} else  {
-  std::cout << GridLogMessage << "Skipping Two Index anti-Symmetric tests "
-                                 "because representation is trivial (dim = 1)"
-            << std::endl;
-}
-
-
-
-
-
-
-
-
 
+#endif
   
   Grid_finalize();
 }

tests/core/Test_reunitarise.cc

@@ -122,14 +122,15 @@ int main (int argc, char ** argv)
   std::cout << "Determinant defect before projection " <<norm2(detU)<<std::endl;
   tmp = U*adj(U) - ident;
   std::cout << "Unitarity check before projection    " << norm2(tmp)<<std::endl; 
-  
+#if (Nc == 3)
   ProjectSU3(U);
   detU= Determinant(U) ;
   detU= detU -1.0;
   std::cout << "Determinant ProjectSU3 defect " <<norm2(detU)<<std::endl;
   tmp = U*adj(U) - ident;
   std::cout << "Unitarity check after projection    " << norm2(tmp)<<std::endl; 
-
+#endif
+  
   ProjectSUn(UU);
   detUU= Determinant(UU);
   detUU= detUU -1.0;

tests/core/Test_wilson_clover.cc

@@ -2,11 +2,12 @@
 
     Grid physics library, www.github.com/paboyle/Grid 
 
-    Source file: ./benchmarks/Benchmark_wilson.cc
+    Source file: ./tests/core/Test_wilson_clover.cc
 
     Copyright (C) 2015
 
     Author: Guido Cossu <guido.cossu@ed.ac.uk>
+            Fabian Joswig <fabian.joswig@ed.ac.uk>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -67,8 +68,6 @@ int main(int argc, char **argv)
   tmp = Zero();
   FermionField err(&Grid);
   err = Zero();
-  FermionField err2(&Grid);
-  err2 = Zero();
   FermionField phi(&Grid);
   random(pRNG, phi);
   FermionField chi(&Grid);
@@ -77,6 +76,8 @@ int main(int argc, char **argv)
   SU<Nc>::HotConfiguration(pRNG, Umu);
   std::vector<LatticeColourMatrix> U(4, &Grid);
 
+  double tolerance = 1e-4;
+
   double volume = 1;
   for (int mu = 0; mu < Nd; mu++)
   {
@@ -88,7 +89,7 @@ int main(int argc, char **argv)
   RealD csw_t = 1.0;
 
   WilsonCloverFermionR Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
-  //Dwc.ImportGauge(Umu); // not necessary, included in the constructor
+  CompactWilsonCloverFermionR Dwc_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
 
   std::cout << GridLogMessage << "==========================================================" << std::endl;
   std::cout << GridLogMessage << "= Testing that Deo + Doe = Dunprec " << std::endl;
@@ -112,7 +113,24 @@ int main(int argc, char **argv)
   setCheckerboard(r_eo, r_e);
 
   err = ref - r_eo;
-  std::cout << GridLogMessage << "EO norm diff   " << norm2(err) << " " << norm2(ref) << " " << norm2(r_eo) << std::endl;
+  std::cout << GridLogMessage << "EO norm diff\t" << norm2(err) << " (" << norm2(ref) << " - " << norm2(r_eo) << ")" << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+
+
+  Dwc_compact.Meooe(src_e, r_o);
+  std::cout << GridLogMessage << "Applied Meo" << std::endl;
+  Dwc_compact.Meooe(src_o, r_e);
+  std::cout << GridLogMessage << "Applied Moe" << std::endl;
+  Dwc_compact.Dhop(src, ref, DaggerNo);
+
+  setCheckerboard(r_eo, r_o);
+  setCheckerboard(r_eo, r_e);
+
+  err = ref - r_eo;
+  std::cout << GridLogMessage << "EO norm diff compact\t" << norm2(err) << " (" << norm2(ref) << " - " << norm2(r_eo) << ")" << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
 
   std::cout << GridLogMessage << "==============================================================" << std::endl;
   std::cout << GridLogMessage << "= Test Ddagger is the dagger of D by requiring                " << std::endl;
@@ -152,6 +170,22 @@ int main(int argc, char **argv)
   std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDce - conj(cDpo) << std::endl;
   std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDco - conj(cDpe) << std::endl;
 
+  Dwc_compact.Meooe(chi_e, dchi_o);
+  Dwc_compact.Meooe(chi_o, dchi_e);
+  Dwc_compact.MeooeDag(phi_e, dphi_o);
+  Dwc_compact.MeooeDag(phi_o, dphi_e);
+
+  pDce = innerProduct(phi_e, dchi_e);
+  pDco = innerProduct(phi_o, dchi_o);
+  cDpe = innerProduct(chi_e, dphi_e);
+  cDpo = innerProduct(chi_o, dphi_o);
+
+  std::cout << GridLogMessage << "e compact " << pDce << " " << cDpe << std::endl;
+  std::cout << GridLogMessage << "o compact " << pDco << " " << cDpo << std::endl;
+
+  std::cout << GridLogMessage << "pDce - conj(cDpo) compact " << pDce - conj(cDpo) << std::endl;
+  std::cout << GridLogMessage << "pDco - conj(cDpe) compact " << pDco - conj(cDpe) << std::endl;
+
   std::cout << GridLogMessage << "==============================================================" << std::endl;
   std::cout << GridLogMessage << "= Test MeeInv Mee = 1   (if csw!=0)                           " << std::endl;
   std::cout << GridLogMessage << "==============================================================" << std::endl;
@@ -169,7 +203,21 @@ int main(int argc, char **argv)
   setCheckerboard(phi, phi_o);
 
   err = phi - chi;
-  std::cout << GridLogMessage << "norm diff   " << norm2(err) << std::endl;
+  std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  Dwc_compact.Mooee(chi_e, src_e);
+  Dwc_compact.MooeeInv(src_e, phi_e);
+
+  Dwc_compact.Mooee(chi_o, src_o);
+  Dwc_compact.MooeeInv(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
 
   std::cout << GridLogMessage << "==============================================================" << std::endl;
   std::cout << GridLogMessage << "= Test MeeDag MeeInvDag = 1    (if csw!=0)                    " << std::endl;
@@ -188,7 +236,21 @@ int main(int argc, char **argv)
   setCheckerboard(phi, phi_o);
 
   err = phi - chi;
-  std::cout << GridLogMessage << "norm diff   " << norm2(err) << std::endl;
+  std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  Dwc_compact.MooeeDag(chi_e, src_e);
+  Dwc_compact.MooeeInvDag(src_e, phi_e);
+
+  Dwc_compact.MooeeDag(chi_o, src_o);
+  Dwc_compact.MooeeInvDag(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
 
   std::cout << GridLogMessage << "==============================================================" << std::endl;
   std::cout << GridLogMessage << "= Test MeeInv MeeDag = 1      (if csw!=0)                     " << std::endl;
@@ -207,7 +269,21 @@ int main(int argc, char **argv)
   setCheckerboard(phi, phi_o);
 
   err = phi - chi;
-  std::cout << GridLogMessage << "norm diff   " << norm2(err) << std::endl;
+  std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  Dwc_compact.MooeeDag(chi_e, src_e);
+  Dwc_compact.MooeeInv(src_e, phi_e);
+
+  Dwc_compact.MooeeDag(chi_o, src_o);
+  Dwc_compact.MooeeInv(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
 
   std::cout << GridLogMessage << "================================================================" << std::endl;
   std::cout << GridLogMessage << "= Testing gauge covariance Clover term with EO preconditioning  " << std::endl;
@@ -249,7 +325,7 @@ int main(int argc, char **argv)
   /////////////////
 
   WilsonCloverFermionR Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
-  Dwc_prime.ImportGauge(U_prime);
+  CompactWilsonCloverFermionR Dwc_compact_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
 
   tmp = Omega * src;
   pickCheckerboard(Even, src_e, tmp);
@@ -262,7 +338,37 @@ int main(int argc, char **argv)
   setCheckerboard(phi, phi_o);
 
   err = chi - adj(Omega) * phi;
-  std::cout << GridLogMessage << "norm diff   " << norm2(err) << std::endl;
+  std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  chi = Zero();
+  phi = Zero();
+  tmp = Zero();
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+
+  Dwc_compact.Mooee(src_e, chi_e);
+  Dwc_compact.Mooee(src_o, chi_o);
+  setCheckerboard(chi, chi_e);
+  setCheckerboard(chi, chi_o);
+  setCheckerboard(src, src_e);
+  setCheckerboard(src, src_o);
+
+  tmp = Omega * src;
+  pickCheckerboard(Even, src_e, tmp);
+  pickCheckerboard(Odd, src_o, tmp);
+
+  Dwc_compact_prime.Mooee(src_e, phi_e);
+  Dwc_compact_prime.Mooee(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = chi - adj(Omega) * phi;
+  std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
 
   std::cout << GridLogMessage << "=================================================================" << std::endl;
   std::cout << GridLogMessage << "= Testing gauge covariance Clover term w/o EO preconditioning  " << std::endl;
@@ -272,7 +378,6 @@ int main(int argc, char **argv)
   phi = Zero();
 
   WilsonFermionR Dw(Umu, Grid, RBGrid, mass, params);
-  Dw.ImportGauge(Umu);
 
   Dw.M(src, result);
   Dwc.M(src, chi);
@@ -280,13 +385,24 @@ int main(int argc, char **argv)
   Dwc_prime.M(Omega * src, phi);
 
   WilsonFermionR Dw_prime(U_prime, Grid, RBGrid, mass, params);
-  Dw_prime.ImportGauge(U_prime);
   Dw_prime.M(Omega * src, result2);
 
+  err = result - adj(Omega) * result2;
+  std::cout << GridLogMessage << "norm diff Wilson              " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
   err = chi - adj(Omega) * phi;
-  err2 = result - adj(Omega) * result2;
+  std::cout << GridLogMessage << "norm diff WilsonClover        " << norm2(err) << std::endl;
-  std::cout << GridLogMessage << "norm diff Wilson   " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
-  std::cout << GridLogMessage << "norm diff WilsonClover  " << norm2(err2) << std::endl;
+
+  chi = Zero();
+  phi = Zero();
+
+  Dwc_compact.M(src, chi);
+  Dwc_compact_prime.M(Omega * src, phi);
+
+  err = chi - adj(Omega) * phi;
+  std::cout << GridLogMessage << "norm diff CompactWilsonClover " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
 
   std::cout << GridLogMessage << "==========================================================" << std::endl;
   std::cout << GridLogMessage << "= Testing Mooee(csw=0) Clover to reproduce Mooee Wilson   " << std::endl;
@@ -296,7 +412,6 @@ int main(int argc, char **argv)
   phi = Zero();
   err = Zero();
   WilsonCloverFermionR Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); //  <-- Notice: csw=0
-  Dwc_csw0.ImportGauge(Umu);
 
   pickCheckerboard(Even, phi_e, phi);
   pickCheckerboard(Odd, phi_o, phi);
@@ -316,7 +431,34 @@ int main(int argc, char **argv)
   setCheckerboard(src, src_o);
 
   err = chi - phi;
-  std::cout << GridLogMessage << "norm diff  " << norm2(err) << std::endl;
+  std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  chi = Zero();
+  phi = Zero();
+  err = Zero();
+  CompactWilsonCloverFermionR Dwc_compact_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, 1.0, anis, params); //  <-- Notice: csw=0
+
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  Dw.Mooee(src_e, chi_e);
+  Dw.Mooee(src_o, chi_o);
+  Dwc_compact_csw0.Mooee(src_e, phi_e);
+  Dwc_compact_csw0.Mooee(src_o, phi_o);
+
+  setCheckerboard(chi, chi_e);
+  setCheckerboard(chi, chi_o);
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+  setCheckerboard(src, src_e);
+  setCheckerboard(src, src_o);
+
+  err = chi - phi;
+  std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
 
   std::cout << GridLogMessage << "==========================================================" << std::endl;
   std::cout << GridLogMessage << "= Testing EO operator is equal to the unprec              " << std::endl;
@@ -348,9 +490,41 @@ int main(int argc, char **argv)
   setCheckerboard(phi, phi_o);
 
   err = ref - phi;
-  std::cout << GridLogMessage << "ref (unpreconditioned operator) diff  :" << norm2(ref) << std::endl;
+  std::cout << GridLogMessage << "ref (unpreconditioned operator) diff         : " << norm2(ref) << std::endl;
-  std::cout << GridLogMessage << "phi (EO decomposition)          diff  :" << norm2(phi) << std::endl;
+  std::cout << GridLogMessage << "phi (EO decomposition)          diff         : " << norm2(phi) << std::endl;
-  std::cout << GridLogMessage << "norm diff                             :" << norm2(err) << std::endl;
+  std::cout << GridLogMessage << "norm diff                                    : " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  chi = Zero();
+  phi = Zero();
+  err = Zero();
+
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  // M phi = (Mooee src_e + Meooe src_o , Meooe src_e + Mooee src_o)
+
+  Dwc_compact.M(src, ref); // Reference result from the unpreconditioned operator
+
+  // EO matrix
+  Dwc_compact.Mooee(src_e, chi_e); 
+  Dwc_compact.Mooee(src_o, chi_o);
+  Dwc_compact.Meooe(src_o, phi_e);
+  Dwc_compact.Meooe(src_e, phi_o);
+
+  phi_o += chi_o;
+  phi_e += chi_e;
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = ref - phi;
+  std::cout << GridLogMessage << "ref (unpreconditioned operator) diff compact : " << norm2(ref) << std::endl;
+  std::cout << GridLogMessage << "phi (EO decomposition)          diff compact : " << norm2(phi) << std::endl;
+  std::cout << GridLogMessage << "norm diff compact                            : " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
 
   Grid_finalize();
 }

tests/core/Test_wilson_conserved_current.cc

@@ -0,0 +1,253 @@
+/*************************************************************************************
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_cayley_cg.cc
+
+Copyright (C) 2022
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Fabian Joswig <fabian.joswig@ed.ac.uk>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+
+template<class What>
+void  TestConserved(What & Dw,
+		    LatticeGaugeField &Umu,
+		    GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+		    GridParallelRNG *RNG4);
+
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT,
+    Gamma::Algebra::Gamma5
+  };
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+								   GridDefaultSimd(Nd,vComplex::Nsimd()),
+								   GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+  std::vector<int> seeds5({5,6,7,8});
+  GridParallelRNG          RNG4(UGrid);
+  std::vector<int> seeds4({1,2,3,4}); RNG4.SeedFixedIntegers(seeds4);
+
+  LatticeGaugeField Umu(UGrid);
+  if( argc > 1 && argv[1][0] != '-' )
+  {
+    std::cout<<GridLogMessage <<"Loading configuration from "<<argv[1]<<std::endl;
+    FieldMetaData header;
+    NerscIO::readConfiguration(Umu, header, argv[1]);
+  }
+  else
+  {
+    std::cout<<GridLogMessage <<"Using hot configuration"<<std::endl;
+    SU<Nc>::HotConfiguration(RNG4,Umu);
+  }
+
+  typename WilsonCloverFermionR::ImplParams params;
+  WilsonAnisotropyCoefficients anis;
+  RealD mass = 0.1;
+  RealD csw_r = 1.0;
+  RealD csw_t = 1.0;
+
+  std::cout<<GridLogMessage <<"=================================="<<std::endl;
+  std::cout<<GridLogMessage <<"WilsonFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"=================================="<<std::endl;
+  WilsonFermionR Dw(Umu,*UGrid,*UrbGrid,mass,params);
+  TestConserved<WilsonFermionR>(Dw,Umu,UGrid,UrbGrid,&RNG4);
+
+  std::cout<<GridLogMessage <<"=================================="<<std::endl;
+  std::cout<<GridLogMessage <<"WilsonCloverFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"=================================="<<std::endl;
+  WilsonCloverFermionR Dwc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, anis, params);
+  TestConserved<WilsonCloverFermionR>(Dwc,Umu,UGrid,UrbGrid,&RNG4);
+
+  std::cout<<GridLogMessage <<"=================================="<<std::endl;
+  std::cout<<GridLogMessage <<"CompactWilsonCloverFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"=================================="<<std::endl;
+  CompactWilsonCloverFermionR Dwcc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, 1.0, anis, params);
+  TestConserved<CompactWilsonCloverFermionR>(Dwcc,Umu,UGrid,UrbGrid,&RNG4);
+
+  std::cout<<GridLogMessage <<"=================================="<<std::endl;
+  std::cout<<GridLogMessage <<"WilsonExpCloverFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"=================================="<<std::endl;
+  WilsonExpCloverFermionR Dewc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, anis, params);
+  TestConserved<WilsonExpCloverFermionR>(Dewc,Umu,UGrid,UrbGrid,&RNG4);
+
+  std::cout<<GridLogMessage <<"=================================="<<std::endl;
+  std::cout<<GridLogMessage <<"CompactWilsonExpCloverFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"=================================="<<std::endl;
+  CompactWilsonExpCloverFermionR Dewcc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, 1.0, anis, params);
+  TestConserved<CompactWilsonExpCloverFermionR>(Dewcc,Umu,UGrid,UrbGrid,&RNG4);
+
+  Grid_finalize();
+}
+
+
+
+template<class Action>
+void  TestConserved(Action & Dw,
+		    LatticeGaugeField &Umu,
+		    GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+		    GridParallelRNG *RNG4)
+{
+  LatticePropagator phys_src(UGrid);
+  LatticePropagator seqsrc(UGrid);
+  LatticePropagator prop4(UGrid);
+  LatticePropagator Vector_mu(UGrid);
+  LatticeComplex    SV (UGrid);
+  LatticeComplex    VV (UGrid);
+  LatticePropagator seqprop(UGrid);
+
+  SpinColourMatrix kronecker; kronecker=1.0;
+  Coordinate coor({0,0,0,0});
+  phys_src=Zero();
+  pokeSite(kronecker,phys_src,coor);
+
+  ConjugateGradient<LatticeFermion> CG(1.0e-16,100000);
+  SchurRedBlackDiagTwoSolve<LatticeFermion> schur(CG);
+  ZeroGuesser<LatticeFermion> zpg;
+  for(int s=0;s<Nd;s++){
+    for(int c=0;c<Nc;c++){
+      LatticeFermion src4  (UGrid);
+      PropToFerm<Action>(src4,phys_src,s,c);
+
+      LatticeFermion result4(UGrid); result4=Zero();
+      schur(Dw,src4,result4,zpg);
+      std::cout<<GridLogMessage<<"spin "<<s<<" color "<<c<<" norm2(sourc4d) "<<norm2(src4)
+               <<" norm2(result4d) "<<norm2(result4)<<std::endl;
+      FermToProp<Action>(prop4,result4,s,c);
+    }
+  }
+
+  auto curr = Current::Vector;
+  const int mu_J=0;
+  const int t_J=0;
+
+  LatticeComplex    ph (UGrid); ph=1.0;
+
+  Dw.SeqConservedCurrent(prop4,
+			   seqsrc,
+			   phys_src,
+			   curr,
+			   mu_J,
+			   t_J,
+			   t_J,// whole lattice
+			   ph);
+
+  for(int s=0;s<Nd;s++){
+    for(int c=0;c<Nc;c++){
+
+      LatticeFermion src4  (UGrid);
+      PropToFerm<Action>(src4,seqsrc,s,c);
+
+      LatticeFermion result4(UGrid); result4=Zero();
+      schur(Dw,src4,result4,zpg);
+
+      FermToProp<Action>(seqprop,result4,s,c);
+    }
+  }
+
+  Gamma g5(Gamma::Algebra::Gamma5);
+  Gamma gT(Gamma::Algebra::GammaT);
+
+  std::vector<TComplex> sumSV;
+  std::vector<TComplex> sumVV;
+
+  Dw.ContractConservedCurrent(prop4,prop4,Vector_mu,phys_src,Current::Vector,Tdir);
+
+  SV       = trace(Vector_mu);        // Scalar-Vector conserved current
+  VV       = trace(gT*Vector_mu);     // (local) Vector-Vector conserved current
+
+  // Spatial sum
+  sliceSum(SV,sumSV,Tdir);
+  sliceSum(VV,sumVV,Tdir);
+
+  const int Nt{static_cast<int>(sumSV.size())};
+
+  std::cout<<GridLogMessage<<"Vector Ward identity by timeslice (~ 0)"<<std::endl;
+  for(int t=0;t<Nt;t++){
+    std::cout<<GridLogMessage <<" t "<<t<<" SV "<<real(TensorRemove(sumSV[t]))<<" VV "<<real(TensorRemove(sumVV[t]))<<std::endl;
+    assert(abs(real(TensorRemove(sumSV[t]))) < 1e-10);
+    assert(abs(real(TensorRemove(sumVV[t]))) < 1e-2);
+  }
+
+  ///////////////////////////////
+  // 3pt vs 2pt check
+  ///////////////////////////////
+  {
+    Gamma::Algebra        gA = Gamma::Algebra::Identity;
+    Gamma                 g(gA);
+
+    LatticePropagator cur(UGrid);
+    LatticePropagator tmp(UGrid);
+    LatticeComplex c(UGrid);
+    SpinColourMatrix qSite;
+    peekSite(qSite, seqprop, coor);
+
+    Complex               test_S, test_V, check_S, check_V;
+
+    std::vector<TComplex> check_buf;
+
+    test_S = trace(qSite*g);
+    test_V = trace(qSite*g*Gamma::gmu[mu_J]);
+
+    Dw.ContractConservedCurrent(prop4,prop4,cur,phys_src,curr,mu_J);
+
+    c = trace(cur*g);
+    sliceSum(c, check_buf, Tp);
+    check_S = TensorRemove(check_buf[t_J]);
+
+    auto gmu=Gamma::gmu[mu_J];
+    c = trace(cur*g*gmu);
+    sliceSum(c, check_buf, Tp);
+    check_V = TensorRemove(check_buf[t_J]);
+
+
+    std::cout<<GridLogMessage << std::setprecision(14)<<"Test S  = " << abs(test_S)   << std::endl;
+    std::cout<<GridLogMessage << "Test V  = " << abs(test_V) << std::endl;
+    std::cout<<GridLogMessage << "Check S = " << abs(check_S) << std::endl;
+    std::cout<<GridLogMessage << "Check V = " << abs(check_V) << std::endl;
+
+    // Check difference = 0
+    check_S = check_S - test_S;
+    check_V = check_V - test_V;
+
+    std::cout<<GridLogMessage << "Consistency check for sequential conserved " <<std::endl;
+    std::cout<<GridLogMessage << "Diff S  = " << abs(check_S) << std::endl;
+    assert(abs(check_S) < 1e-8);
+    std::cout<<GridLogMessage << "Diff V  = " << abs(check_V) << std::endl;
+    assert(abs(check_V) < 1e-8);
+  }
+
+}

tests/core/Test_wilson_exp_clover.cc

@@ -0,0 +1,530 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/core/Test_wilson_exp_clover.cc
+
+    Copyright (C) 2022
+
+    Author: Guido Cossu <guido.cossu@ed.ac.uk>
+            Fabian Joswig <fabian.joswig@ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  auto latt_size = GridDefaultLatt();
+  auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  auto mpi_layout = GridDefaultMpi();
+  GridCartesian Grid(latt_size, simd_layout, mpi_layout);
+  GridRedBlackCartesian RBGrid(&Grid);
+
+  int threads = GridThread::GetThreads();
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+  std::cout << GridLogMessage << "Grid floating point word size is REALF" << sizeof(RealF) << std::endl;
+  std::cout << GridLogMessage << "Grid floating point word size is REALD" << sizeof(RealD) << std::endl;
+  std::cout << GridLogMessage << "Grid floating point word size is REAL" << sizeof(Real) << std::endl;
+
+  std::vector<int> seeds({1, 2, 3, 4});
+  GridParallelRNG pRNG(&Grid);
+  pRNG.SeedFixedIntegers(seeds);
+  //  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});
+
+  typedef typename WilsonExpCloverFermionR::FermionField FermionField;
+  typename WilsonExpCloverFermionR::ImplParams params;
+  WilsonAnisotropyCoefficients anis;
+
+  FermionField src(&Grid);
+  random(pRNG, src);
+  FermionField result(&Grid);
+  result = Zero();
+  FermionField result2(&Grid);
+  result2 = Zero();
+  FermionField ref(&Grid);
+  ref = Zero();
+  FermionField tmp(&Grid);
+  tmp = Zero();
+  FermionField err(&Grid);
+  err = Zero();
+  FermionField phi(&Grid);
+  random(pRNG, phi);
+  FermionField chi(&Grid);
+  random(pRNG, chi);
+  LatticeGaugeField Umu(&Grid);
+  SU<Nc>::HotConfiguration(pRNG, Umu);
+  std::vector<LatticeColourMatrix> U(4, &Grid);
+
+  double tolerance = 1e-4;
+
+  double volume = 1;
+  for (int mu = 0; mu < Nd; mu++)
+  {
+    volume = volume * latt_size[mu];
+  }
+
+  RealD mass = 0.1;
+  RealD csw_r = 1.0;
+  RealD csw_t = 1.0;
+
+  WilsonExpCloverFermionR Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
+  CompactWilsonExpCloverFermionR Dwc_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
+
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+  std::cout << GridLogMessage << "= Testing that Deo + Doe = Dunprec " << std::endl;
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+
+  FermionField src_e(&RBGrid);
+  FermionField src_o(&RBGrid);
+  FermionField r_e(&RBGrid);
+  FermionField r_o(&RBGrid);
+  FermionField r_eo(&Grid);
+  pickCheckerboard(Even, src_e, src);
+  pickCheckerboard(Odd, src_o, src);
+
+  Dwc.Meooe(src_e, r_o);
+  std::cout << GridLogMessage << "Applied Meo" << std::endl;
+  Dwc.Meooe(src_o, r_e);
+  std::cout << GridLogMessage << "Applied Moe" << std::endl;
+  Dwc.Dhop(src, ref, DaggerNo);
+
+  setCheckerboard(r_eo, r_o);
+  setCheckerboard(r_eo, r_e);
+
+  err = ref - r_eo;
+  std::cout << GridLogMessage << "EO norm diff\t" << norm2(err) << " (" << norm2(ref) << " - " << norm2(r_eo) << ")" << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+
+
+  Dwc_compact.Meooe(src_e, r_o);
+  std::cout << GridLogMessage << "Applied Meo" << std::endl;
+  Dwc_compact.Meooe(src_o, r_e);
+  std::cout << GridLogMessage << "Applied Moe" << std::endl;
+  Dwc_compact.Dhop(src, ref, DaggerNo);
+
+  setCheckerboard(r_eo, r_o);
+  setCheckerboard(r_eo, r_e);
+
+  err = ref - r_eo;
+  std::cout << GridLogMessage << "EO norm diff compact\t" << norm2(err) << " (" << norm2(ref) << " - " << norm2(r_eo) << ")" << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+  std::cout << GridLogMessage << "= Test Ddagger is the dagger of D by requiring                " << std::endl;
+  std::cout << GridLogMessage << "=  < phi | Deo | chi > * = < chi | Deo^dag| phi>  " << std::endl;
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+
+  FermionField chi_e(&RBGrid);
+  FermionField chi_o(&RBGrid);
+
+  FermionField dchi_e(&RBGrid);
+  FermionField dchi_o(&RBGrid);
+
+  FermionField phi_e(&RBGrid);
+  FermionField phi_o(&RBGrid);
+
+  FermionField dphi_e(&RBGrid);
+  FermionField dphi_o(&RBGrid);
+
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+
+  Dwc.Meooe(chi_e, dchi_o);
+  Dwc.Meooe(chi_o, dchi_e);
+  Dwc.MeooeDag(phi_e, dphi_o);
+  Dwc.MeooeDag(phi_o, dphi_e);
+
+  ComplexD pDce = innerProduct(phi_e, dchi_e);
+  ComplexD pDco = innerProduct(phi_o, dchi_o);
+  ComplexD cDpe = innerProduct(chi_e, dphi_e);
+  ComplexD cDpo = innerProduct(chi_o, dphi_o);
+
+  std::cout << GridLogMessage << "e " << pDce << " " << cDpe << std::endl;
+  std::cout << GridLogMessage << "o " << pDco << " " << cDpo << std::endl;
+
+  std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDce - conj(cDpo) << std::endl;
+  std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDco - conj(cDpe) << std::endl;
+
+  Dwc_compact.Meooe(chi_e, dchi_o);
+  Dwc_compact.Meooe(chi_o, dchi_e);
+  Dwc_compact.MeooeDag(phi_e, dphi_o);
+  Dwc_compact.MeooeDag(phi_o, dphi_e);
+
+  pDce = innerProduct(phi_e, dchi_e);
+  pDco = innerProduct(phi_o, dchi_o);
+  cDpe = innerProduct(chi_e, dphi_e);
+  cDpo = innerProduct(chi_o, dphi_o);
+
+  std::cout << GridLogMessage << "e compact " << pDce << " " << cDpe << std::endl;
+  std::cout << GridLogMessage << "o compact " << pDco << " " << cDpo << std::endl;
+
+  std::cout << GridLogMessage << "pDce - conj(cDpo) compact " << pDce - conj(cDpo) << std::endl;
+  std::cout << GridLogMessage << "pDco - conj(cDpe) compact " << pDco - conj(cDpe) << std::endl;
+
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+  std::cout << GridLogMessage << "= Test MeeInv Mee = 1   (if csw!=0)                           " << std::endl;
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  Dwc.Mooee(chi_e, src_e);
+  Dwc.MooeeInv(src_e, phi_e);
+
+  Dwc.Mooee(chi_o, src_o);
+  Dwc.MooeeInv(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  Dwc_compact.Mooee(chi_e, src_e);
+  Dwc_compact.MooeeInv(src_e, phi_e);
+
+  Dwc_compact.Mooee(chi_o, src_o);
+  Dwc_compact.MooeeInv(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+  std::cout << GridLogMessage << "= Test MeeDag MeeInvDag = 1    (if csw!=0)                    " << std::endl;
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  Dwc.MooeeDag(chi_e, src_e);
+  Dwc.MooeeInvDag(src_e, phi_e);
+
+  Dwc.MooeeDag(chi_o, src_o);
+  Dwc.MooeeInvDag(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  Dwc_compact.MooeeDag(chi_e, src_e);
+  Dwc_compact.MooeeInvDag(src_e, phi_e);
+
+  Dwc_compact.MooeeDag(chi_o, src_o);
+  Dwc_compact.MooeeInvDag(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+  std::cout << GridLogMessage << "= Test MeeInv MeeDag = 1      (if csw!=0)                     " << std::endl;
+  std::cout << GridLogMessage << "==============================================================" << std::endl;
+
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  Dwc.MooeeDag(chi_e, src_e);
+  Dwc.MooeeInv(src_e, phi_e);
+
+  Dwc.MooeeDag(chi_o, src_o);
+  Dwc.MooeeInv(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  Dwc_compact.MooeeDag(chi_e, src_e);
+  Dwc_compact.MooeeInv(src_e, phi_e);
+
+  Dwc_compact.MooeeDag(chi_o, src_o);
+  Dwc_compact.MooeeInv(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = phi - chi;
+  std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  std::cout << GridLogMessage << "================================================================" << std::endl;
+  std::cout << GridLogMessage << "= Testing gauge covariance Clover term with EO preconditioning  " << std::endl;
+  std::cout << GridLogMessage << "================================================================" << std::endl;
+
+  chi = Zero();
+  phi = Zero();
+  tmp = Zero();
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+
+  Dwc.Mooee(src_e, chi_e);
+  Dwc.Mooee(src_o, chi_o);
+  setCheckerboard(chi, chi_e);
+  setCheckerboard(chi, chi_o);
+  setCheckerboard(src, src_e);
+  setCheckerboard(src, src_o);
+
+  ////////////////////// Gauge Transformation
+  std::vector<int> seeds2({5, 6, 7, 8});
+  GridParallelRNG pRNG2(&Grid);
+  pRNG2.SeedFixedIntegers(seeds2);
+  LatticeColourMatrix Omega(&Grid);
+  LatticeColourMatrix ShiftedOmega(&Grid);
+  LatticeGaugeField U_prime(&Grid);
+  U_prime = Zero();
+  LatticeColourMatrix U_prime_mu(&Grid);
+  U_prime_mu = Zero();
+  SU<Nc>::LieRandomize(pRNG2, Omega, 1.0);
+  for (int mu = 0; mu < Nd; mu++)
+  {
+    U[mu] = peekLorentz(Umu, mu);
+    ShiftedOmega = Cshift(Omega, mu, 1);
+    U_prime_mu = Omega * U[mu] * adj(ShiftedOmega);
+    pokeLorentz(U_prime, U_prime_mu, mu);
+  }
+  /////////////////
+
+  WilsonExpCloverFermionR Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
+  CompactWilsonExpCloverFermionR Dwc_compact_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
+
+  tmp = Omega * src;
+  pickCheckerboard(Even, src_e, tmp);
+  pickCheckerboard(Odd, src_o, tmp);
+
+  Dwc_prime.Mooee(src_e, phi_e);
+  Dwc_prime.Mooee(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = chi - adj(Omega) * phi;
+  std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  chi = Zero();
+  phi = Zero();
+  tmp = Zero();
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+
+  Dwc_compact.Mooee(src_e, chi_e);
+  Dwc_compact.Mooee(src_o, chi_o);
+  setCheckerboard(chi, chi_e);
+  setCheckerboard(chi, chi_o);
+  setCheckerboard(src, src_e);
+  setCheckerboard(src, src_o);
+
+  tmp = Omega * src;
+  pickCheckerboard(Even, src_e, tmp);
+  pickCheckerboard(Odd, src_o, tmp);
+
+  Dwc_compact_prime.Mooee(src_e, phi_e);
+  Dwc_compact_prime.Mooee(src_o, phi_o);
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = chi - adj(Omega) * phi;
+  std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  std::cout << GridLogMessage << "=================================================================" << std::endl;
+  std::cout << GridLogMessage << "= Testing gauge covariance Clover term w/o EO preconditioning  " << std::endl;
+  std::cout << GridLogMessage << "================================================================" << std::endl;
+
+  chi = Zero();
+  phi = Zero();
+
+  WilsonFermionR Dw(Umu, Grid, RBGrid, mass, params);
+
+  Dw.M(src, result);
+  Dwc.M(src, chi);
+
+  Dwc_prime.M(Omega * src, phi);
+
+  WilsonFermionR Dw_prime(U_prime, Grid, RBGrid, mass, params);
+  Dw_prime.M(Omega * src, result2);
+
+  err = result - adj(Omega) * result2;
+  std::cout << GridLogMessage << "norm diff Wilson                 " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+  err = chi - adj(Omega) * phi;
+  std::cout << GridLogMessage << "norm diff WilsonExpClover        " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  chi = Zero();
+  phi = Zero();
+
+  Dwc_compact.M(src, chi);
+  Dwc_compact_prime.M(Omega * src, phi);
+
+  err = chi - adj(Omega) * phi;
+  std::cout << GridLogMessage << "norm diff CompactWilsonExpClover " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+  std::cout << GridLogMessage << "= Testing Mooee(csw=0) Clover to reproduce Mooee Wilson   " << std::endl;
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+
+  chi = Zero();
+  phi = Zero();
+  err = Zero();
+  WilsonExpCloverFermionR Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); //  <-- Notice: csw=0
+
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  Dw.Mooee(src_e, chi_e);
+  Dw.Mooee(src_o, chi_o);
+  Dwc_csw0.Mooee(src_e, phi_e);
+  Dwc_csw0.Mooee(src_o, phi_o);
+
+  setCheckerboard(chi, chi_e);
+  setCheckerboard(chi, chi_o);
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+  setCheckerboard(src, src_e);
+  setCheckerboard(src, src_o);
+
+  err = chi - phi;
+  std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  chi = Zero();
+  phi = Zero();
+  err = Zero();
+  CompactWilsonExpCloverFermionR Dwc_compact_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, 1.0, anis, params); //  <-- Notice: csw=0
+
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  Dw.Mooee(src_e, chi_e);
+  Dw.Mooee(src_o, chi_o);
+  Dwc_compact_csw0.Mooee(src_e, phi_e);
+  Dwc_compact_csw0.Mooee(src_o, phi_o);
+
+  setCheckerboard(chi, chi_e);
+  setCheckerboard(chi, chi_o);
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+  setCheckerboard(src, src_e);
+  setCheckerboard(src, src_o);
+
+  err = chi - phi;
+  std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+  std::cout << GridLogMessage << "= Testing EO operator is equal to the unprec              " << std::endl;
+  std::cout << GridLogMessage << "==========================================================" << std::endl;
+
+  chi = Zero();
+  phi = Zero();
+  err = Zero();
+
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  // M phi = (Mooee src_e + Meooe src_o , Meooe src_e + Mooee src_o)
+
+  Dwc.M(src, ref); // Reference result from the unpreconditioned operator
+
+  // EO matrix
+  Dwc.Mooee(src_e, chi_e); 
+  Dwc.Mooee(src_o, chi_o);
+  Dwc.Meooe(src_o, phi_e);
+  Dwc.Meooe(src_e, phi_o);
+
+  phi_o += chi_o;
+  phi_e += chi_e;
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = ref - phi;
+  std::cout << GridLogMessage << "ref (unpreconditioned operator) diff         : " << norm2(ref) << std::endl;
+  std::cout << GridLogMessage << "phi (EO decomposition)          diff         : " << norm2(phi) << std::endl;
+  std::cout << GridLogMessage << "norm diff                                    : " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  chi = Zero();
+  phi = Zero();
+  err = Zero();
+
+  pickCheckerboard(Even, phi_e, phi);
+  pickCheckerboard(Odd, phi_o, phi);
+  pickCheckerboard(Even, chi_e, chi);
+  pickCheckerboard(Odd, chi_o, chi);
+
+  // M phi = (Mooee src_e + Meooe src_o , Meooe src_e + Mooee src_o)
+
+  Dwc_compact.M(src, ref); // Reference result from the unpreconditioned operator
+
+  // EO matrix
+  Dwc_compact.Mooee(src_e, chi_e); 
+  Dwc_compact.Mooee(src_o, chi_o);
+  Dwc_compact.Meooe(src_o, phi_e);
+  Dwc_compact.Meooe(src_e, phi_o);
+
+  phi_o += chi_o;
+  phi_e += chi_e;
+
+  setCheckerboard(phi, phi_e);
+  setCheckerboard(phi, phi_o);
+
+  err = ref - phi;
+  std::cout << GridLogMessage << "ref (unpreconditioned operator) diff compact : " << norm2(ref) << std::endl;
+  std::cout << GridLogMessage << "phi (EO decomposition)          diff compact : " << norm2(phi) << std::endl;
+  std::cout << GridLogMessage << "norm diff compact                            : " << norm2(err) << std::endl;
+  assert(fabs(norm2(err)) < tolerance);
+
+  Grid_finalize();
+}

tests/hmc/Test_hmc_ScalarActionNxN.cc

@@ -132,8 +132,8 @@ int main(int argc, char **argv) {
 
   // Checkpointer definition
   CheckpointerParameters CPparams(Reader);
-  //TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
+  TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
-  TheHMC.Resources.LoadScidacCheckpointer(CPparams, SPar);
+  //TheHMC.Resources.LoadScidacCheckpointer(CPparams, SPar); this breaks for compilation without lime
 
   RNGModuleParameters RNGpar(Reader);
   TheHMC.Resources.SetRNGSeeds(RNGpar);

tests/hmc/Test_hmc_WG_Production.cc

@@ -74,10 +74,10 @@ int main(int argc, char **argv) {
 
   // Checkpointer definition
   CheckpointerParameters CPparams(Reader);
-  //TheHMC.Resources.LoadNerscCheckpointer(CPparams);
+  TheHMC.Resources.LoadNerscCheckpointer(CPparams);
 
-  // Store metadata in the Scidac checkpointer
+  // Store metadata in the Scidac checkpointer - obviously breaks without LIME
-  TheHMC.Resources.LoadScidacCheckpointer(CPparams, WilsonPar);
+  //TheHMC.Resources.LoadScidacCheckpointer(CPparams, WilsonPar);
 
   RNGModuleParameters RNGpar(Reader);
   TheHMC.Resources.SetRNGSeeds(RNGpar);

tests/lanczos/Test_compressed_lanczos.cc

@@ -37,6 +37,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 using namespace std;
 using namespace Grid;
 
+#ifdef HAVE_LIME
+
 template<class Fobj,class CComplex,int nbasis>
 class LocalCoherenceLanczosScidac : public LocalCoherenceLanczos<Fobj,CComplex,nbasis>
 { 
@@ -249,3 +251,11 @@ int main (int argc, char ** argv) {
   Grid_finalize();
 }
 
+#else
+
+int main( void )
+{
+  return 0 ;
+}
+
+#endif // HAVE_LIME_H

tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc

@@ -36,7 +36,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 using namespace std;
 using namespace Grid;
- ;
+
+#ifdef HAVE_LIME
 
 template<class Fobj,class CComplex,int nbasis>
 class LocalCoherenceLanczosScidac : public LocalCoherenceLanczos<Fobj,CComplex,nbasis>
@@ -250,3 +251,11 @@ int main (int argc, char ** argv) {
   Grid_finalize();
 }
 
+#else
+
+int main( void )
+{
+  return 0 ;
+}
+
+#endif

tests/solver/Test_coarse_even_odd.cc

@@ -93,8 +93,16 @@ int main(int argc, char** argv) {
   //                    Setup of Dirac Matrix and Operator                   //
   /////////////////////////////////////////////////////////////////////////////
 
-  LatticeGaugeField Umu(Grid_f); SU3::HotConfiguration(pRNG_f, Umu);
+  LatticeGaugeField Umu(Grid_f);
-
+#if (Nc==2)
+  SU2::HotConfiguration(pRNG_f, Umu);
+#elif (defined Nc==3)
+  SU3::HotConfiguration(pRNG_f, Umu);
+#elif (defined Nc==4)
+  SU4::HotConfiguration(pRNG_f, Umu);
+#elif (defined Nc==5)
+  SU5::HotConfiguration(pRNG_f, Umu);
+#endif
   RealD checkTolerance = (getPrecision<LatticeFermion>::value == 1) ? 1e-7 : 1e-15;
 
   RealD mass = -0.30;

tests/solver/Test_dwf_mrhs_cg.cc

@@ -34,6 +34,7 @@ using namespace Grid;
 
 int main (int argc, char ** argv)
 {
+#ifdef HAVE_LIME 
   typedef typename DomainWallFermionR::FermionField FermionField; 
   typedef typename DomainWallFermionR::ComplexField ComplexField; 
   typename DomainWallFermionR::ImplParams params; 
@@ -237,4 +238,5 @@ int main (int argc, char ** argv)
   }
 
   Grid_finalize();
+#endif // HAVE_LIME
 }

tests/solver/Test_wilsonclover_cg_prec.cc

@@ -71,7 +71,12 @@ int main (int argc, char ** argv)
   RealD mass = -0.1;
   RealD csw_r = 1.0;
   RealD csw_t = 1.0;
+  RealD cF = 1.0;
   WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t);
+  CompactWilsonCloverFermionR Dw_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0);
+  WilsonExpCloverFermionR Dwe(Umu, Grid, RBGrid, mass, csw_r, csw_t);
+  CompactWilsonExpCloverFermionR Dwe_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0);
+
 
   //  HermitianOperator<WilsonFermion,LatticeFermion> HermOp(Dw);
   //  ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
@@ -80,12 +85,28 @@ int main (int argc, char ** argv)
   LatticeFermion    src_o(&RBGrid);
   LatticeFermion result_o(&RBGrid);
   pickCheckerboard(Odd,src_o,src);
-  result_o=Zero();
 
-  SchurDiagMooeeOperator<WilsonCloverFermionR,LatticeFermion> HermOpEO(Dw);
   ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
+
+  std::cout << GridLogMessage << "Testing Wilson Clover" << std::endl;
+  SchurDiagMooeeOperator<WilsonCloverFermionR,LatticeFermion> HermOpEO(Dw);
+  result_o=Zero();
   CG(HermOpEO,src_o,result_o);
-  
+
+  std::cout << GridLogMessage << "Testing Compact Wilson Clover" << std::endl;
+  SchurDiagMooeeOperator<CompactWilsonCloverFermionR,LatticeFermion> HermOpEO_compact(Dw_compact);
+  result_o=Zero();
+  CG(HermOpEO_compact,src_o,result_o);
+
+  std::cout << GridLogMessage << "Testing Wilson Exp Clover" << std::endl;
+  SchurDiagMooeeOperator<WilsonExpCloverFermionR,LatticeFermion> HermOpEO_exp(Dwe);
+  result_o=Zero();
+  CG(HermOpEO_exp,src_o,result_o);
+
+  std::cout << GridLogMessage << "Testing Compact Wilson Exp Clover" << std::endl;
+  SchurDiagMooeeOperator<CompactWilsonExpCloverFermionR,LatticeFermion> HermOpEO_exp_compact(Dwe_compact);
+  result_o=Zero();
+  CG(HermOpEO_exp_compact,src_o,result_o);
 
   Grid_finalize();
 }

tests/solver/Test_wilsonclover_cg_schur.cc

@@ -60,18 +60,36 @@ int main (int argc, char ** argv)
   LatticeGaugeField Umu(&Grid); SU<Nc>::HotConfiguration(pRNG,Umu);
 
   LatticeFermion    src(&Grid); random(pRNG,src);
-  LatticeFermion result(&Grid); result=Zero();
+  LatticeFermion result(&Grid);
-  LatticeFermion resid(&Grid); 
+  LatticeFermion resid(&Grid);
-
-  RealD mass = -0.1;
-  RealD csw_r = 1.0;
-  RealD csw_t = 1.0;
-  WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t);
 
   ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
   SchurRedBlackDiagMooeeSolve<LatticeFermion> SchurSolver(CG);
 
+  RealD mass = -0.1;
+  RealD csw_r = 1.0;
+  RealD csw_t = 1.0;
+  RealD cF = 1.0;
+
+  std::cout << GridLogMessage << "Testing Wilson Clover" << std::endl;
+  WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t);
+  result=Zero();
   SchurSolver(Dw,src,result);
+
+  std::cout << GridLogMessage << "Testing Compact Wilson Clover" << std::endl;
+  CompactWilsonCloverFermionR Dw_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0);
+  result=Zero();
+  SchurSolver(Dw_compact,src,result);
+
+  std::cout << GridLogMessage << "Testing Wilson Exp Clover" << std::endl;
+  WilsonExpCloverFermionR Dwe(Umu, Grid, RBGrid, mass, csw_r, csw_t);
+  result=Zero();
+  SchurSolver(Dwe,src,result);
+
+  std::cout << GridLogMessage << "Testing Compact Wilson Exp Clover" << std::endl;
+  CompactWilsonExpCloverFermionR Dwe_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0);
+  result=Zero();
+  SchurSolver(Dwe_compact,src,result);
   
   Grid_finalize();
 }

tests/solver/Test_wilsonclover_cg_unprec.cc

@@ -59,7 +59,7 @@ int main (int argc, char ** argv)
 
   LatticeFermion src(&Grid); random(pRNG,src);
   RealD nrm = norm2(src);
-  LatticeFermion result(&Grid); result=Zero();
+  LatticeFermion result(&Grid);
   LatticeGaugeField Umu(&Grid); SU<Nc>::HotConfiguration(pRNG,Umu);
 
   double volume=1;
@@ -70,11 +70,34 @@ int main (int argc, char ** argv)
   RealD mass = -0.1;
   RealD csw_r = 1.0;
   RealD csw_t = 1.0;
+  RealD cF = 1.0;
   WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t);
+  CompactWilsonCloverFermionR Dw_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0);
+  WilsonExpCloverFermionR Dwe(Umu, Grid, RBGrid, mass, csw_r, csw_t);
+  CompactWilsonExpCloverFermionR Dwe_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0);
 
-  MdagMLinearOperator<WilsonFermionR,LatticeFermion> HermOp(Dw);
   ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
+
+  std::cout << GridLogMessage << "Testing Wilson Clover" << std::endl;
+  MdagMLinearOperator<WilsonCloverFermionR,LatticeFermion> HermOp(Dw);
+  result=Zero();
   CG(HermOp,src,result);
 
+  std::cout << GridLogMessage << "Testing Compact Wilson Clover" << std::endl;
+  MdagMLinearOperator<CompactWilsonCloverFermionR,LatticeFermion> HermOp_compact(Dw_compact);
+  result=Zero();
+  CG(HermOp_compact,src,result);
+
+
+  std::cout << GridLogMessage << "Testing Wilson Exp Clover" << std::endl;
+  MdagMLinearOperator<WilsonExpCloverFermionR,LatticeFermion> HermOp_exp(Dwe);
+  result=Zero();
+  CG(HermOp_exp,src,result);
+
+  std::cout << GridLogMessage << "Testing Compact Wilson Exp Clover" << std::endl;
+  MdagMLinearOperator<CompactWilsonExpCloverFermionR,LatticeFermion> HermOp_exp_compact(Dwe_compact);
+  result=Zero();
+  CG(HermOp_exp_compact,src,result);
+
   Grid_finalize();
 }