Addedd Ta functionality to the tensor types

Merge remote-tracking branch 'upstream/master'

Conflicts:
	configure

Makefile.am

@@ -1,3 +1,5 @@
 # additional include paths necessary to compile the C++ library
 AM_CXXFLAGS = -I$(top_srcdir)/
 SUBDIRS = lib tests benchmarks
+
+filelist: $(SUBDIRS)

benchmarks/Benchmark_comms.cc

@@ -77,11 +77,12 @@ int main (int argc, char ** argv)
       }
       double stop=usecond();
 
-      double xbytes    = Nloop*bytes*2*ncomm;
+      double dbytes    = bytes;
+      double xbytes    = Nloop*dbytes*2.0*ncomm;
       double rbytes    = xbytes;
       double bidibytes = xbytes+rbytes;
 
-      double time = stop-start;
+      double time = stop-start; // microseconds
 
       std::cout << lat<<"\t\t"<<Ls<<"\t\t"<<bytes<<"\t\t"<<xbytes/time<<"\t\t"<<bidibytes/time<<std::endl;
     }
@@ -151,8 +152,9 @@ int main (int argc, char ** argv)
       }
 
       double stop=usecond();
-
+      
-      double xbytes    = Nloop*bytes*2*ncomm;
+      double dbytes    = bytes;
+      double xbytes    = Nloop*dbytes*2.0*ncomm;
       double rbytes    = xbytes;
       double bidibytes = xbytes+rbytes;
 

benchmarks/Benchmark_dwf.cc

@@ -0,0 +1,156 @@
+#include <Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+template<class d>
+struct scal {
+  d internal;
+};
+
+  Gamma::GammaMatrix Gmu [] = {
+    Gamma::GammaX,
+    Gamma::GammaY,
+    Gamma::GammaZ,
+    Gamma::GammaT
+  };
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  std::vector<int> latt4 = GridDefaultLatt();
+  const int Ls=8;
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  std::vector<int> seeds4({1,2,3,4});
+  std::vector<int> seeds5({5,6,7,8});
+
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+  GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
+
+  LatticeFermion src   (FGrid); random(RNG5,src);
+  LatticeFermion result(FGrid); result=zero;
+  LatticeFermion    ref(FGrid);    ref=zero;
+  LatticeFermion    tmp(FGrid);
+  LatticeFermion    err(FGrid);
+
+  ColourMatrix cm = Complex(1.0,0.0);
+
+  LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
+  LatticeGaugeField Umu5d(FGrid); 
+
+  // replicate across fifth dimension
+  for(int ss=0;ss<Umu._grid->oSites();ss++){
+    for(int s=0;s<Ls;s++){
+      Umu5d._odata[Ls*ss+s] = Umu._odata[ss];
+    }
+  }
+
+  ////////////////////////////////////
+  // Naive wilson implementation
+  ////////////////////////////////////
+  std::vector<LatticeColourMatrix> U(4,FGrid);
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = peekIndex<LorentzIndex>(Umu5d,mu);
+  }
+
+  if (1)
+  {
+    ref = zero;
+    for(int mu=0;mu<Nd;mu++){
+
+      tmp = U[mu]*Cshift(src,mu+1,1);
+      ref=ref + tmp + Gamma(Gmu[mu])*tmp;
+
+      tmp =adj(U[mu])*src;
+      tmp =Cshift(tmp,mu+1,-1);
+      ref=ref + tmp - Gamma(Gmu[mu])*tmp;
+    }
+    ref = -0.5*ref;
+  }
+
+  RealD mass=0.1;
+  RealD M5  =1.8;
+  DomainWallFermion Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+  
+  std::cout << "Calling Dw"<<std::endl;
+  int ncall=10;
+  double t0=usecond();
+  for(int i=0;i<ncall;i++){
+    Dw.Dhop(src,result,0);
+  }
+  double t1=usecond();
+
+  double volume=Ls;  for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
+  double flops=1344*volume*ncall;
+  
+  std::cout << "Called Dw"<<std::endl;
+  std::cout << "norm result "<< norm2(result)<<std::endl;
+  std::cout << "norm ref    "<< norm2(ref)<<std::endl;
+  std::cout << "mflop/s =   "<< flops/(t1-t0)<<std::endl;
+  err = ref-result; 
+  std::cout << "norm diff   "<< norm2(err)<<std::endl;
+
+
+  if (1)
+  { // Naive wilson dag implementation
+    ref = zero;
+    for(int mu=0;mu<Nd;mu++){
+
+      //    ref =  src - Gamma(Gamma::GammaX)* src ; // 1+gamma_x
+      tmp = U[mu]*Cshift(src,mu+1,1);
+      for(int i=0;i<ref._odata.size();i++){
+	ref._odata[i]+= tmp._odata[i] - Gamma(Gmu[mu])*tmp._odata[i]; ;
+      }
+
+      tmp =adj(U[mu])*src;
+      tmp =Cshift(tmp,mu+1,-1);
+      for(int i=0;i<ref._odata.size();i++){
+	ref._odata[i]+= tmp._odata[i] + Gamma(Gmu[mu])*tmp._odata[i]; ;
+      }
+    }
+    ref = -0.5*ref;
+  }
+  Dw.Dhop(src,result,1);
+  std::cout << "Called DwDag"<<std::endl;
+  std::cout << "norm result "<< norm2(result)<<std::endl;
+  std::cout << "norm ref    "<< norm2(ref)<<std::endl;
+  err = ref-result; 
+  std::cout << "norm diff   "<< norm2(err)<<std::endl;
+
+  LatticeFermion src_e (FrbGrid);
+  LatticeFermion src_o (FrbGrid);
+  LatticeFermion r_e   (FrbGrid);
+  LatticeFermion r_o   (FrbGrid);
+  LatticeFermion r_eo  (FGrid);
+
+
+  std::cout << "Calling Deo and Doe"<<std::endl;
+  pickCheckerboard(Even,src_e,src);
+  pickCheckerboard(Odd,src_o,src);
+
+  Dw.DhopEO(src_o,r_e,DaggerNo);
+  Dw.DhopOE(src_e,r_o,DaggerNo);
+  Dw.Dhop(src,result,DaggerNo);
+
+  setCheckerboard(r_eo,r_o);
+  setCheckerboard(r_eo,r_e);
+
+  err = r_eo-result; 
+  std::cout << "norm diff   "<< norm2(err)<<std::endl;
+
+  pickCheckerboard(Even,src_e,err);
+  pickCheckerboard(Odd,src_o,err);
+  std::cout << "norm diff even  "<< norm2(src_e)<<std::endl;
+  std::cout << "norm diff odd   "<< norm2(src_o)<<std::endl;
+
+  Grid_finalize();
+}

benchmarks/Benchmark_wilson.cc

@@ -49,6 +49,7 @@ int main (int argc, char ** argv)
   }  
 
   // Only one non-zero (y)
+#if 0
   Umu=zero;
   Complex cone(1.0,0.0);
   for(int nn=0;nn<Nd;nn++){
@@ -59,6 +60,7 @@ int main (int argc, char ** argv)
     }
     pokeIndex<LorentzIndex>(Umu,U[nn],nn);
   }
+#endif
 
   for(int mu=0;mu<Nd;mu++){
     U[mu] = peekIndex<LorentzIndex>(Umu,mu);
@@ -80,9 +82,9 @@ int main (int argc, char ** argv)
       }
     }
   }
-
+  ref = -0.5*ref;
   RealD mass=0.1;
-  WilsonMatrix Dw(Umu,Grid,RBGrid,mass);
+  WilsonFermion Dw(Umu,Grid,RBGrid,mass);
   
   std::cout << "Calling Dw"<<std::endl;
   int ncall=10000;
@@ -91,7 +93,7 @@ int main (int argc, char ** argv)
     Dw.Dhop(src,result,0);
   }
   double t1=usecond();
-  double flops=1320*volume*ncall;
+  double flops=1344*volume*ncall;
   
   std::cout << "Called Dw"<<std::endl;
   std::cout << "norm result "<< norm2(result)<<std::endl;
@@ -129,6 +131,7 @@ int main (int argc, char ** argv)
       }
     }
   }
+  ref = -0.5*ref;
   Dw.Dhop(src,result,1);
   std::cout << "Called DwDag"<<std::endl;
   std::cout << "norm result "<< norm2(result)<<std::endl;

benchmarks/Make.inc

@@ -0,0 +1,23 @@
+
+bin_PROGRAMS = Benchmark_comms Benchmark_dwf Benchmark_memory_bandwidth Benchmark_su3 Benchmark_wilson
+
+
+Benchmark_comms_SOURCES=Benchmark_comms.cc
+Benchmark_comms_LDADD=-lGrid
+
+
+Benchmark_dwf_SOURCES=Benchmark_dwf.cc
+Benchmark_dwf_LDADD=-lGrid
+
+
+Benchmark_memory_bandwidth_SOURCES=Benchmark_memory_bandwidth.cc
+Benchmark_memory_bandwidth_LDADD=-lGrid
+
+
+Benchmark_su3_SOURCES=Benchmark_su3.cc
+Benchmark_su3_LDADD=-lGrid
+
+
+Benchmark_wilson_SOURCES=Benchmark_wilson.cc
+Benchmark_wilson_LDADD=-lGrid
+

benchmarks/Makefile.am

@@ -5,29 +5,4 @@ AM_LDFLAGS = -L$(top_builddir)/lib
 #
 # Test code
 #
-bin_PROGRAMS = Grid_wilson Grid_comms Grid_memory_bandwidth Grid_su3 Grid_wilson_cg_unprec Grid_wilson_evenodd  Grid_wilson_cg_prec Grid_wilson_cg_schur
+include Make.inc
-
-Grid_wilson_SOURCES = Grid_wilson.cc
-Grid_wilson_LDADD = -lGrid
-
-Grid_wilson_evenodd_SOURCES = Grid_wilson_evenodd.cc
-Grid_wilson_evenodd_LDADD = -lGrid
-
-Grid_wilson_cg_unprec_SOURCES = Grid_wilson_cg_unprec.cc
-Grid_wilson_cg_unprec_LDADD = -lGrid
-
-Grid_wilson_cg_prec_SOURCES = Grid_wilson_cg_prec.cc
-Grid_wilson_cg_prec_LDADD = -lGrid
-
-Grid_wilson_cg_schur_SOURCES = Grid_wilson_cg_schur.cc
-Grid_wilson_cg_schur_LDADD = -lGrid
-
-Grid_comms_SOURCES = Grid_comms.cc
-Grid_comms_LDADD = -lGrid
-
-Grid_su3_SOURCES = Grid_su3.cc Grid_su3_test.cc Grid_su3_expr.cc
-Grid_su3_LDADD = -lGrid
-
-Grid_memory_bandwidth_SOURCES = Grid_memory_bandwidth.cc
-Grid_memory_bandwidth_LDADD = -lGrid
-

configure

@@ -3154,7 +3154,7 @@ END
 fi
 
 
-ac_config_headers="$ac_config_headers lib/Grid_config.h"
+ac_config_headers="$ac_config_headers lib/GridConfig.h"
 
 # Check whether --enable-silent-rules was given.
 if test "${enable_silent_rules+set}" = set; then :
@@ -7671,7 +7671,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1
 for ac_config_target in $ac_config_targets
 do
   case $ac_config_target in
-    "lib/Grid_config.h") CONFIG_HEADERS="$CONFIG_HEADERS lib/Grid_config.h" ;;
+    "lib/GridConfig.h") CONFIG_HEADERS="$CONFIG_HEADERS lib/GridConfig.h" ;;
     "depfiles") CONFIG_COMMANDS="$CONFIG_COMMANDS depfiles" ;;
     "docs/doxy.cfg") CONFIG_FILES="$CONFIG_FILES docs/doxy.cfg" ;;
     "Makefile") CONFIG_FILES="$CONFIG_FILES Makefile" ;;

configure.ac

@@ -11,7 +11,7 @@ AC_CANONICAL_SYSTEM
 AM_INIT_AUTOMAKE(subdir-objects)
 AC_CONFIG_MACRO_DIR([m4])
 AC_CONFIG_SRCDIR([lib/Grid.h])
-AC_CONFIG_HEADERS([lib/Grid_config.h])
+AC_CONFIG_HEADERS([lib/GridConfig.h])
 m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
 
 AC_MSG_NOTICE([

lib/Cartesian.h

@@ -0,0 +1,8 @@
+#ifndef GRID_CARTESIAN_H
+#define GRID_CARTESIAN_H
+
+#include <cartesian/Cartesian_base.h>
+#include <cartesian/Cartesian_full.h>
+#include <cartesian/Cartesian_red_black.h> 
+
+#endif

lib/Communicator.h

@@ -1,6 +1,6 @@
 #ifndef GRID_COMMUNICATOR_H
 #define GRID_COMMUNICATOR_H
 
-#include <communicator/Grid_communicator_base.h>
+#include <communicator/Communicator_base.h>
 
 #endif

lib/Comparison.h

@@ -141,7 +141,7 @@ namespace Grid {
     }
 }
 
-#include <lattice/Grid_lattice_comparison.h>
+#include <lattice/Lattice_comparison.h>
-#include <lattice/Grid_lattice_where.h>
+#include <lattice/Lattice_where.h>
 
 #endif

lib/Cshift.h

@@ -1,13 +1,13 @@
 #ifndef _GRID_CSHIFT_H_
 #define _GRID_CSHIFT_H_
 
-#include <cshift/Grid_cshift_common.h>
+#include <cshift/Cshift_common.h>
 
 #ifdef GRID_COMMS_NONE
-#include <cshift/Grid_cshift_none.h>
+#include <cshift/Cshift_none.h>
 #endif
 
 #ifdef GRID_COMMS_MPI
-#include <cshift/Grid_cshift_mpi.h>
+#include <cshift/Cshift_mpi.h>
 #endif 
 #endif

lib/Grid.h

@@ -33,7 +33,7 @@
 
 #define strong_inline __attribute__((always_inline)) inline
 
-#include <Grid_config.h>
+#include <GridConfig.h>
 
 ////////////////////////////////////////////////////////////
 // Tunable header includes
@@ -46,22 +46,22 @@
 #include <malloc.h>
 #endif
 
-#include <Grid_aligned_allocator.h>
+#include <AlignedAllocator.h>
 
-#include <Grid_simd.h>
+#include <Simd.h>
-#include <Grid_threads.h>
+#include <Threads.h>
 
-#include <Grid_cartesian.h> // subdir aggregate
+#include <Communicator.h> // subdir aggregate
-#include <Grid_math.h>      // subdir aggregate
+#include <Cartesian.h> // subdir aggregate
-#include <Grid_lattice.h>   // subdir aggregate
+#include <Tensors.h>   // subdir aggregate
-#include <Grid_comparison.h>
+#include <Lattice.h>   // subdir aggregate
-#include <Grid_cshift.h>    // subdir aggregate
+#include <Comparison.h>
-#include <Grid_stencil.h>   // subdir aggregate
+#include <Cshift.h>    // subdir aggregate
+#include <Stencil.h>   // subdir aggregate
+#include <Algorithms.h>// subdir aggregate
 
-#include <Grid_algorithms.h>// subdir aggregate
+#include <qcd/QCD.h>
-
+#include <parallelIO/NerscIO.h>
-#include <qcd/Grid_qcd.h>
-#include <parallelIO/GridNerscIO.h>
 
 namespace Grid {
 

lib/GridConfig.h

@@ -1,5 +1,5 @@
-/* lib/Grid_config.h.  Generated from Grid_config.h.in by configure.  */
+/* lib/GridConfig.h.  Generated from GridConfig.h.in by configure.  */
-/* lib/Grid_config.h.in.  Generated from configure.ac by autoheader.  */
+/* lib/GridConfig.h.in.  Generated from configure.ac by autoheader.  */
 
 /* AVX */
 /* #undef AVX1 */

lib/GridConfig.h.in

@@ -1,4 +1,4 @@
-/* lib/Grid_config.h.in.  Generated from configure.ac by autoheader.  */
+/* lib/GridConfig.h.in.  Generated from configure.ac by autoheader.  */
 
 /* AVX */
 #undef AVX1

lib/GridInit.cc

@@ -142,7 +142,11 @@ void Grid_init(int *argc,char ***argv)
     Grid_quiesce_nodes();
   }
   if( GridCmdOptionExists(*argv,*argv+*argc,"--dslash-opt") ){
-    WilsonMatrix::HandOptDslash=1;
+    WilsonFermion::HandOptDslash=1;
+    WilsonFermion5D::HandOptDslash=1;
+  }
+  if( GridCmdOptionExists(*argv,*argv+*argc,"--lebesgue") ){
+    LebesgueOrder::UseLebesgueOrder=1;
   }
   GridParseLayout(*argv,*argc,
 		  Grid_default_latt,

lib/Grid_cartesian.h

@@ -1,8 +0,0 @@
-#ifndef GRID_CARTESIAN_H
-#define GRID_CARTESIAN_H
-
-#include <cartesian/Grid_cartesian_base.h>
-#include <cartesian/Grid_cartesian_full.h>
-#include <cartesian/Grid_cartesian_red_black.h> 
-
-#endif

lib/Grid_math.h

@@ -1,16 +0,0 @@
-#ifndef GRID_MATH_H
-#define GRID_MATH_H
-
-#include <math/Grid_math_traits.h>
-#include <math/Grid_math_tensors.h>
-#include <math/Grid_math_arith.h>
-#include <math/Grid_math_inner.h>
-#include <math/Grid_math_outer.h>
-#include <math/Grid_math_transpose.h>
-#include <math/Grid_math_trace.h>
-#include <math/Grid_math_peek.h>
-#include <math/Grid_math_poke.h>
-#include <math/Grid_math_reality.h>
-
-    
-#endif

lib/Lattice.h

@@ -1,6 +1,6 @@
 #ifndef GRID_LATTICE_H
 #define GRID_LATTICE_H
 
-#include <lattice/Grid_lattice_base.h>
+#include <lattice/Lattice_base.h>
 
 #endif

lib/Make.inc

@@ -0,0 +1,4 @@
+
+HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vRealF.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_poke.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_peek.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/SpaceTimeGrid.h ./qcd/LinalgUtils.h ./qcd/TwoSpinor.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Threads.h ./Comparison.h ./Grid.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./stencil/Lebesgue.h
+
+CCFILES=./qcd/SpaceTimeGrid.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/Dirac.cc ./GridInit.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

lib/Makefile.am

@@ -3,100 +3,26 @@ AM_CXXFLAGS = -I$(top_srcdir)/
 
 extra_sources=
 if BUILD_COMMS_MPI
-  extra_sources+=communicator/Grid_communicator_mpi.cc
+  extra_sources+=communicator/Communicator_mpi.cc
 endif
 
 if BUILD_COMMS_NONE
-  extra_sources+=communicator/Grid_communicator_none.cc
+  extra_sources+=communicator/Communicator_none.cc
 endif
 
 #
 # Libraries
 #
-lib_LIBRARIES = libGrid.a
-libGrid_a_SOURCES =				\
-	Grid_init.cc				\
-	stencil/Grid_stencil_common.cc		\
-	qcd/Grid_qcd_dirac.cc			\
-	qcd/Grid_qcd_dhop.cc			\
-	qcd/Grid_qcd_dhop_hand.cc		\
-	qcd/Grid_qcd_wilson_dop.cc		\
-	algorithms/approx/Zolotarev.cc		\
-	algorithms/approx/Remez.cc		\
-	$(extra_sources)
 
+include Make.inc
+
+lib_LIBRARIES = libGrid.a
+libGrid_a_SOURCES = $(CCFILES) $(extra_sources)
+
+
+#	qcd/action/fermion/PartialFractionFermion5D.cc\	\
 #
 # Include files
 #
-nobase_include_HEADERS = algorithms/approx/bigfloat.h		\
+nobase_include_HEADERS=$(HFILES)
-	algorithms/approx/Chebyshev.h				\
-	algorithms/approx/Remez.h				\
-	algorithms/approx/Zolotarev.h				\
-	algorithms/iterative/ConjugateGradient.h		\
-	algorithms/iterative/NormalEquations.h			\
-	algorithms/iterative/SchurRedBlack.h			\
-	algorithms/LinearOperator.h				\
-	algorithms/SparseMatrix.h				\
-	cartesian/Grid_cartesian_base.h				\
-	cartesian/Grid_cartesian_full.h				\
-	cartesian/Grid_cartesian_red_black.h			\
-	communicator/Grid_communicator_base.h			\
-	cshift/Grid_cshift_common.h				\
-	cshift/Grid_cshift_mpi.h				\
-	cshift/Grid_cshift_none.h				\
-	Grid.h							\
-	Grid_algorithms.h					\
-	Grid_aligned_allocator.h				\
-	Grid_cartesian.h					\
-	Grid_communicator.h					\
-	Grid_comparison.h					\
-	Grid_cshift.h						\
-	Grid_extract.h						\
-	Grid_lattice.h						\
-	Grid_math.h						\
-	Grid_simd.h						\
-	Grid_stencil.h						\
-	Grid_threads.h						\
-	lattice/Grid_lattice_arith.h				\
-	lattice/Grid_lattice_base.h				\
-	lattice/Grid_lattice_comparison.h			\
-	lattice/Grid_lattice_conformable.h			\
-	lattice/Grid_lattice_coordinate.h			\
-	lattice/Grid_lattice_ET.h				\
-	lattice/Grid_lattice_local.h				\
-	lattice/Grid_lattice_overload.h				\
-	lattice/Grid_lattice_peekpoke.h				\
-	lattice/Grid_lattice_reality.h				\
-	lattice/Grid_lattice_reduction.h			\
-	lattice/Grid_lattice_rng.h				\
-	lattice/Grid_lattice_trace.h				\
-	lattice/Grid_lattice_transfer.h				\
-	lattice/Grid_lattice_transpose.h			\
-	lattice/Grid_lattice_where.h				\
-	math/Grid_math_arith.h					\
-	math/Grid_math_arith_add.h				\
-	math/Grid_math_arith_mac.h				\
-	math/Grid_math_arith_mul.h				\
-	math/Grid_math_arith_scalar.h				\
-	math/Grid_math_arith_sub.h				\
-	math/Grid_math_inner.h					\
-	math/Grid_math_outer.h					\
-	math/Grid_math_peek.h					\
-	math/Grid_math_poke.h					\
-	math/Grid_math_reality.h				\
-	math/Grid_math_tensors.h				\
-	math/Grid_math_trace.h					\
-	math/Grid_math_traits.h					\
-	math/Grid_math_transpose.h				\
-	parallelIO/GridNerscIO.h				\
-	qcd/Grid_qcd.h						\
-	qcd/Grid_qcd_2spinor.h					\
-	qcd/Grid_qcd_dirac.h					\
-	qcd/Grid_qcd_wilson_dop.h				\
-	simd/Grid_vector_types.h				\
-	simd/Grid_sse4.h					\
-	simd/Grid_avx.h						\
-	simd/Grid_avx512.h					
-
-
 

lib/Stencil.h

@@ -1,6 +1,8 @@
 #ifndef GRID_STENCIL_H
 #define GRID_STENCIL_H
 
+#include <stencil/Lebesgue.h>   // subdir aggregate
+
 //////////////////////////////////////////////////////////////////////////////////////////
 // Must not lose sight that goal is to be able to construct really efficient
 // gather to a point stencil code. CSHIFT is not the best way, so need
@@ -38,29 +40,12 @@
 //////////////////////////////////////////////////////////////////////////////////////////
 
 namespace Grid {
-
+  
-
 
   class CartesianStencil { // Stencil runs along coordinate axes only; NO diagonal fill in.
   public:
 
       typedef uint32_t StencilInteger;
-      
-
-
-      StencilInteger alignup(StencilInteger n){
-	n--;           // 1000 0011 --> 1000 0010
-	n |= n >> 1;   // 1000 0010 | 0100 0001 = 1100 0011
-	n |= n >> 2;   // 1100 0011 | 0011 0000 = 1111 0011
-	n |= n >> 4;   // 1111 0011 | 0000 1111 = 1111 1111
-	n |= n >> 8;   // ... (At this point all bits are 1, so further bitwise-or
-	n |= n >> 16;  //      operations produce no effect.)
-	n++;           // 1111 1111 --> 1 0000 0000
-	return n;
-      };
-      void LebesgueOrder(void);
-
-      std::vector<StencilInteger> _LebesgueReorder;
 
       int                               _checkerboard;
       int                               _npoints; // Move to template param?
@@ -131,8 +116,8 @@ namespace Grid {
 	  // Gather phase
 	  int sshift [2];
 	  if ( comm_dim ) {
-	    sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,0);
+	    sshift[0] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Even);
-	    sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,1);
+	    sshift[1] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Odd);
 	    if ( sshift[0] == sshift[1] ) {
 	      if (splice_dim) {
 		GatherStartCommsSimd(source,dimension,shift,0x3,u_comm_buf,u_comm_offset,compress);
@@ -179,8 +164,8 @@ namespace Grid {
 	  std::vector<cobj,alignedAllocator<cobj> > send_buf(buffer_size); // hmm...
 	  std::vector<cobj,alignedAllocator<cobj> > recv_buf(buffer_size);
 	  
-	  int cb= (cbmask==0x2)? 1 : 0;
+	  int cb= (cbmask==0x2)? Odd : Even;
-	  int sshift= _grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,cb);
+	  int sshift= _grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,cb);
 	  
 	  for(int x=0;x<rd;x++){       
 	    
@@ -266,8 +251,8 @@ namespace Grid {
 	  // Work out what to send where
 	  ///////////////////////////////////////////
 
-	  int cb    = (cbmask==0x2)? 1 : 0;
+	  int cb    = (cbmask==0x2)? Odd : Even;
-	  int sshift= _grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,cb);
+	  int sshift= _grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,cb);
 	  
 	  // loop over outer coord planes orthog to dim
 	  for(int x=0;x<rd;x++){       

lib/Tensors.h

@@ -0,0 +1,17 @@
+#ifndef GRID_MATH_H
+#define GRID_MATH_H
+
+#include <tensors/Tensor_traits.h>
+#include <tensors/Tensor_class.h>
+#include <tensors/Tensor_arith.h>
+#include <tensors/Tensor_inner.h>
+#include <tensors/Tensor_outer.h>
+#include <tensors/Tensor_transpose.h>
+#include <tensors/Tensor_trace.h>
+#include <tensors/Tensor_Ta.h>
+#include <tensors/Tensor_peek.h>
+#include <tensors/Tensor_poke.h>
+#include <tensors/Tensor_reality.h>
+#include <tensors/Tensor_extract_merge.h>
+    
+#endif

lib/Threads.h

@@ -9,7 +9,7 @@
 
 #ifdef GRID_OMP
 #include <omp.h>
-#define PARALLEL_FOR_LOOP _Pragma("omp parallel for")
+#define PARALLEL_FOR_LOOP _Pragma("omp parallel for ")
 #define PARALLEL_NESTED_LOOP2 _Pragma("omp parallel for collapse(2)")
 #else
 #define PARALLEL_FOR_LOOP 

lib/algorithms/LinearOperator.h

@@ -125,39 +125,7 @@ namespace Grid {
      };
     */
 
-    // Chroma interface defining GaugeAction
-    /*
-      template<typename P, typename Q>   class GaugeAction
-  virtual const CreateGaugeState<P,Q>& getCreateState() const = 0;
-  virtual GaugeState<P,Q>* createState(const Q& q) const
-  virtual const GaugeBC<P,Q>& getGaugeBC() const
-  virtual const Set& getSet(void) const = 0;
-  virtual void deriv(P& result, const Handle< GaugeState<P,Q> >& state) const 
-  virtual Double S(const Handle< GaugeState<P,Q> >& state) const = 0;
 
-  class LinearGaugeAction : public GaugeAction< multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
-  typedef multi1d<LatticeColorMatrix>  P;
-  typedef multi1d<LatticeColorMatrix>  Q;
-  virtual void staple(LatticeColorMatrix& result,
-		      const Handle< GaugeState<P,Q> >& state,
-		      int mu, int cb) const = 0;
-    */
-
-    // Chroma interface defining FermionAction
-    /*
-     template<typename T, typename P, typename Q>  class FermAct4D : public FermionAction<T,P,Q>
-     virtual LinearOperator<T>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
-     virtual LinearOperator<T>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
-     virtual LinOpSystemSolver<T>* invLinOp(Handle< FermState<T,P,Q> > state,
-     virtual MdagMSystemSolver<T>* invMdagM(Handle< FermState<T,P,Q> > state,
-     virtual LinOpMultiSystemSolver<T>* mInvLinOp(Handle< FermState<T,P,Q> > state,
-     virtual MdagMMultiSystemSolver<T>* mInvMdagM(Handle< FermState<T,P,Q> > state,
-     virtual MdagMMultiSystemSolverAccumulate<T>* mInvMdagMAcc(Handle< FermState<T,P,Q> > state,
-     virtual SystemSolver<T>* qprop(Handle< FermState<T,P,Q> > state,
-     class DiffFermAct4D : public FermAct4D<T,P,Q>
-     virtual DiffLinearOperator<T,Q,P>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
-     virtual DiffLinearOperator<T,Q,P>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
-    */
 }
 
 #endif

lib/algorithms/SparseMatrix.h

@@ -10,7 +10,7 @@ namespace Grid {
   /////////////////////////////////////////////////////////////////////////////////////////////
     template<class Field> class SparseMatrixBase {
     public:
-      GridBase *_grid;
+      virtual GridBase *Grid(void) =0;
       // Full checkerboar operations
       virtual RealD M    (const Field &in, Field &out)=0;
       virtual RealD Mdag (const Field &in, Field &out)=0;
@@ -19,7 +19,6 @@ namespace Grid {
 	ni=M(in,tmp);
 	no=Mdag(tmp,out);
       }
-      SparseMatrixBase(GridBase *grid) : _grid(grid) {};
     };
 
   /////////////////////////////////////////////////////////////////////////////////////////////
@@ -27,7 +26,7 @@ namespace Grid {
   /////////////////////////////////////////////////////////////////////////////////////////////
     template<class Field> class CheckerBoardedSparseMatrixBase : public SparseMatrixBase<Field> {
     public:
-      GridBase *_cbgrid;
+      virtual GridBase *RedBlackGrid(void)=0;
       // half checkerboard operaions
       virtual  void Meooe    (const Field &in, Field &out)=0;
       virtual  void Mooee    (const Field &in, Field &out)=0;
@@ -62,9 +61,7 @@ namespace Grid {
 	Field tmp(in._grid);
 	ni=Mpc(in,tmp);
 	no=MpcDag(tmp,out);
-	//	std::cout<<"MpcDagMpc "<<ni<<" "<<no<<std::endl;
       }
-      CheckerBoardedSparseMatrixBase(GridBase *grid,GridBase *cbgrid) : SparseMatrixBase<Field>(grid), _cbgrid(cbgrid) {};
     };
 
 }

lib/algorithms/approx/Zolotarev.cc

@@ -58,6 +58,8 @@
 
 /* Compute the partial fraction expansion coefficients (alpha) from the
  * factored form */
+namespace Grid {
+namespace Approx {
 
 static void construct_partfrac(izd *z) {
   int dn = z -> dn, dd = z -> dd, type = z -> type;
@@ -291,7 +293,7 @@ static void sncndnFK(INTERNAL_PRECISION u, INTERNAL_PRECISION k,
  * Set type = 0 for the Zolotarev approximation, which is zero at x = 0, and
  * type = 1 for the approximation which is infinite at x = 0. */
 
-zolotarev_data* bfm_zolotarev(PRECISION epsilon, int n, int type) {
+zolotarev_data* grid_zolotarev(PRECISION epsilon, int n, int type) {
   INTERNAL_PRECISION A, c, cp, kp, ksq, sn, cn, dn, Kp, Kj, z, z0, t, M, F,
     l, invlambda, xi, xisq, *tv, s, opl;
   int m, czero, ts;
@@ -412,7 +414,7 @@ zolotarev_data* bfm_zolotarev(PRECISION epsilon, int n, int type) {
   return zd;
 }
 
-zolotarev_data* bfm_higham(PRECISION epsilon, int n) {
+zolotarev_data* grid_higham(PRECISION epsilon, int n) {
   INTERNAL_PRECISION A, M, c, cp, z, z0, t, epssq;
   int m, czero;
   zolotarev_data *zd;
@@ -502,6 +504,7 @@ zolotarev_data* bfm_higham(PRECISION epsilon, int n) {
   free(d);
   return zd;
 }
+}}
 
 #ifdef TEST
 
@@ -707,4 +710,6 @@ int main(int argc, char** argv) {
 
   return EXIT_SUCCESS;
 }
+
+
 #endif /* TEST */

lib/algorithms/approx/Zolotarev.h

@@ -1,7 +1,8 @@
 /* -*- Mode: C; comment-column: 22; fill-column: 79; -*- */
 
 #ifdef __cplusplus
-extern "C" {
+namespace Grid {
+namespace Approx {
 #endif
 
 #define HVERSION Header Time-stamp: <14-OCT-2004 09:26:51.00 adk@MISSCONTRARY>
@@ -76,10 +77,10 @@ typedef struct {
  * zolotarev_data structure. The arguments must satisfy the constraints that
  * epsilon > 0, n > 0, and type = 0 or 1. */
 
-ZOLOTAREV_DATA* bfm_higham(PRECISION epsilon, int n) ;
+ZOLOTAREV_DATA* grid_higham(PRECISION epsilon, int n) ;
-ZOLOTAREV_DATA* bfm_zolotarev(PRECISION epsilon, int n, int type);
+ZOLOTAREV_DATA* grid_zolotarev(PRECISION epsilon, int n, int type);
 #endif
 
 #ifdef __cplusplus
-}
+}}
 #endif

lib/algorithms/iterative/SchurRedBlack.h

@@ -60,8 +60,8 @@ namespace Grid {
 
       // FIXME CGdiagonalMee not implemented virtual function
       // FIXME use CBfactorise to control schur decomp
-      GridBase *grid = _Matrix._cbgrid;
+      GridBase *grid = _Matrix.RedBlackGrid();
-      GridBase *fgrid= _Matrix._grid;
+      GridBase *fgrid= _Matrix.Grid();
  
       Field src_e(grid);
       Field src_o(grid);

lib/cartesian/Cartesian_base.h

@@ -2,7 +2,6 @@
 #define GRID_CARTESIAN_BASE_H
 
 #include <Grid.h>
-#include <Grid_communicator.h>
 
 namespace Grid{
 
@@ -21,7 +20,7 @@ public:
     // Give Lattice access
     template<class object> friend class Lattice;
 
-    GridBase(std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};
+    GridBase(const std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};
 
 
     // Physics Grid information.
@@ -52,16 +51,13 @@ public:
     ////////////////////////////////////////////////////////////////
     virtual int CheckerBoarded(int dim)=0;
     virtual int CheckerBoard(std::vector<int> site)=0;
-    virtual int CheckerBoardDestination(int source_cb,int shift)=0;
+    virtual int CheckerBoardDestination(int source_cb,int shift,int dim)=0;
     virtual int CheckerBoardShift(int source_cb,int dim,int shift,int osite)=0;
-    inline int  CheckerBoardFromOindex (int Oindex){
+    virtual int CheckerBoardShiftForCB(int source_cb,int dim,int shift,int cb)=0;
+    int  CheckerBoardFromOindex (int Oindex){
       std::vector<int> ocoor;
       oCoorFromOindex(ocoor,Oindex); 
-      int ss=0;
+      return CheckerBoard(ocoor);
-      for(int d=0;d<_ndimension;d++){
-	ss=ss+ocoor[d];
-      }      
-      return ss&0x1;
     }
 
     //////////////////////////////////////////////////////////////////////////////////////////////

lib/cartesian/Cartesian_full.h

@@ -18,15 +18,18 @@ public:
     virtual int CheckerBoard(std::vector<int> site){
         return 0;
     }
-    virtual int CheckerBoardDestination(int cb,int shift){
+    virtual int CheckerBoardDestination(int cb,int shift,int dim){
         return 0;
     }
-    virtual int CheckerBoardShift(int source_cb,int dim,int shift, int osite){
+    virtual int CheckerBoardShiftForCB(int source_cb,int dim,int shift, int ocb){
-        return shift;
+      return shift;
     }
-    GridCartesian(std::vector<int> &dimensions,
+    virtual int CheckerBoardShift(int source_cb,int dim,int shift, int osite){
-		  std::vector<int> &simd_layout,
+      return shift;
-		  std::vector<int> &processor_grid
+    }
+    GridCartesian(const std::vector<int> &dimensions,
+		  const std::vector<int> &simd_layout,
+		  const std::vector<int> &processor_grid
 		  ) : GridBase(processor_grid)
     {
         ///////////////////////

lib/cartesian/Cartesian_red_black.h

@@ -0,0 +1,196 @@
+#ifndef GRID_CARTESIAN_RED_BLACK_H
+#define GRID_CARTESIAN_RED_BLACK_H
+
+
+namespace Grid {
+
+    static const int CbRed  =0;
+    static const int CbBlack=1;
+    static const int Even   =CbRed;
+    static const int Odd    =CbBlack;
+
+    // Perhaps these are misplaced and 
+    // should be in sparse matrix.
+    // Also should make these a named enum type
+    static const int DaggerNo=0;
+    static const int DaggerYes=1;
+
+// Specialise this for red black grids storing half the data like a chess board.
+class GridRedBlackCartesian : public GridBase
+{
+public:
+    std::vector<int> _checker_dim_mask;
+    int              _checker_dim;
+
+    virtual int CheckerBoarded(int dim){
+      if( dim==_checker_dim) return 1;
+      else return 0;
+    }
+    virtual int CheckerBoard(std::vector<int> site){
+      int linear=0;
+      assert(site.size()==_ndimension);
+      for(int d=0;d<_ndimension;d++){ 
+	if(_checker_dim_mask[d])
+	  linear=linear+site[d];
+      }
+      return (linear&0x1);
+    }
+
+
+    // Depending on the cb of site, we toggle source cb.
+    // for block #b, element #e = (b, e)
+    // we need 
+    virtual int CheckerBoardShiftForCB(int source_cb,int dim,int shift,int ocb){
+      if(dim != _checker_dim) return shift;
+
+      int fulldim =_fdimensions[dim];
+      shift = (shift+fulldim)%fulldim;
+
+      // Probably faster with table lookup;
+      // or by looping over x,y,z and multiply rather than computing checkerboard.
+	  
+      if ( (source_cb+ocb)&1 ) {
+
+	return (shift)/2;
+      } else {
+	return (shift+1)/2;
+      }
+    }
+    virtual int CheckerBoardShift(int source_cb,int dim,int shift,int osite){
+
+      if(dim != _checker_dim) return shift;
+
+      int ocb=CheckerBoardFromOindex(osite);
+      
+      return CheckerBoardShiftForCB(source_cb,dim,shift,ocb);
+    }
+    
+    virtual int CheckerBoardDestination(int source_cb,int shift,int dim){
+      if ( _checker_dim_mask[dim]  ) {
+	// If _fdimensions[checker_dim] is odd, then shifting by 1 in other dims
+	// does NOT cause a parity hop.
+	int add=(dim==_checker_dim) ? 0 : _fdimensions[_checker_dim];
+        if ( (shift+add) &0x1) {
+            return 1-source_cb;
+        } else {
+            return source_cb;
+        }
+      } else {
+	return source_cb;
+
+      }
+    };
+
+    GridRedBlackCartesian(const GridBase *base) : GridRedBlackCartesian(base->_fdimensions,base->_simd_layout,base->_processors)  {};
+
+    GridRedBlackCartesian(const std::vector<int> &dimensions,
+			  const std::vector<int> &simd_layout,
+			  const std::vector<int> &processor_grid,
+			  const std::vector<int> &checker_dim_mask,
+			  int checker_dim
+			  ) :  GridBase(processor_grid) 
+    {
+      Init(dimensions,simd_layout,processor_grid,checker_dim_mask,checker_dim);
+    }
+    GridRedBlackCartesian(const std::vector<int> &dimensions,
+			  const std::vector<int> &simd_layout,
+			  const std::vector<int> &processor_grid) : GridBase(processor_grid) 
+    {
+      std::vector<int> checker_dim_mask(dimensions.size(),1);
+      Init(dimensions,simd_layout,processor_grid,checker_dim_mask,0);
+    }
+    void Init(const std::vector<int> &dimensions,
+	      const std::vector<int> &simd_layout,
+	      const std::vector<int> &processor_grid,
+	      const std::vector<int> &checker_dim_mask,
+	      int checker_dim)
+    {
+    ///////////////////////
+    // Grid information
+    ///////////////////////
+      _checker_dim = checker_dim;
+      assert(checker_dim_mask[checker_dim]==1);
+      _ndimension = dimensions.size();
+      assert(checker_dim_mask.size()==_ndimension);
+      assert(processor_grid.size()==_ndimension);
+      assert(simd_layout.size()==_ndimension);
+      
+      _fdimensions.resize(_ndimension);
+      _gdimensions.resize(_ndimension);
+      _ldimensions.resize(_ndimension);
+      _rdimensions.resize(_ndimension);
+      _simd_layout.resize(_ndimension);
+      
+      _ostride.resize(_ndimension);
+      _istride.resize(_ndimension);
+      
+      _fsites = _gsites = _osites = _isites = 1;
+	
+      _checker_dim_mask=checker_dim_mask;
+
+      for(int d=0;d<_ndimension;d++){
+	_fdimensions[d] = dimensions[d];
+	_gdimensions[d] = _fdimensions[d];
+	_fsites = _fsites * _fdimensions[d];
+	_gsites = _gsites * _gdimensions[d];
+        
+	if (d==_checker_dim) {
+	  _gdimensions[d] = _gdimensions[d]/2; // Remove a checkerboard
+	}
+	_ldimensions[d] = _gdimensions[d]/_processors[d];
+
+	// Use a reduced simd grid
+	_simd_layout[d] = simd_layout[d];
+	_rdimensions[d]= _ldimensions[d]/_simd_layout[d];
+	
+	_osites *= _rdimensions[d];
+	_isites *= _simd_layout[d];
+        
+	// Addressing support
+	if ( d==0 ) {
+	  _ostride[d] = 1;
+	  _istride[d] = 1;
+	} else {
+	  _ostride[d] = _ostride[d-1]*_rdimensions[d-1];
+	  _istride[d] = _istride[d-1]*_simd_layout[d-1];
+	}
+      }
+            
+      ////////////////////////////////////////////////////////////////////////////////////////////
+      // subplane information
+      ////////////////////////////////////////////////////////////////////////////////////////////
+      _slice_block.resize(_ndimension);
+      _slice_stride.resize(_ndimension);
+      _slice_nblock.resize(_ndimension);
+        
+      int block =1;
+      int nblock=1;
+      for(int d=0;d<_ndimension;d++) nblock*=_rdimensions[d];
+      
+      for(int d=0;d<_ndimension;d++){
+	nblock/=_rdimensions[d];
+	_slice_block[d] =block;
+	_slice_stride[d]=_ostride[d]*_rdimensions[d];
+	_slice_nblock[d]=nblock;
+	block = block*_rdimensions[d];
+      }
+      
+    };
+protected:
+    virtual int oIndex(std::vector<int> &coor)
+    {
+      int idx=0;
+      for(int d=0;d<_ndimension;d++) {
+	if( d==_checker_dim ) {
+	  idx+=_ostride[d]*((coor[d]/2)%_rdimensions[d]);
+	} else {
+	  idx+=_ostride[d]*(coor[d]%_rdimensions[d]);
+	}
+      }
+        return idx;
+    };
+        
+};
+
+}
+#endif

lib/cartesian/Grid_cartesian_red_black.h

@@ -1,134 +0,0 @@
-#ifndef GRID_CARTESIAN_RED_BLACK_H
-#define GRID_CARTESIAN_RED_BLACK_H
-
-
-namespace Grid {
-
-    static const int CbRed  =0;
-    static const int CbBlack=1;
-    static const int Even   =CbRed;
-    static const int Odd    =CbBlack;
-    static const int DaggerNo=0;
-    static const int DaggerYes=1;
-
-// Specialise this for red black grids storing half the data like a chess board.
-class GridRedBlackCartesian : public GridBase
-{
-public:
-    virtual int CheckerBoarded(int dim){
-      if( dim==0) return 1;
-      else return 0;
-    }
-    virtual int CheckerBoard(std::vector<int> site){
-      return (site[0]+site[1]+site[2]+site[3])&0x1;
-    }
-
-    // Depending on the cb of site, we toggle source cb.
-    // for block #b, element #e = (b, e)
-    // we need 
-    virtual int CheckerBoardShift(int source_cb,int dim,int shift,int osite){
-
-      if(dim != 0) return shift;
-
-      int fulldim =_fdimensions[0];
-      shift = (shift+fulldim)%fulldim;
-
-      // Probably faster with table lookup;
-      // or by looping over x,y,z and multiply rather than computing checkerboard.
-      int ocb=CheckerBoardFromOindex(osite);
-	  
-      if ( (source_cb+ocb)&1 ) {
-	return (shift)/2;
-      } else {
-	return (shift+1)/2;
-      }
-    }
-
-    virtual int CheckerBoardDestination(int source_cb,int shift){
-        if ((shift+_fdimensions[0])&0x1) {
-            return 1-source_cb;
-        } else {
-            return source_cb;
-        }
-    };
-
-    GridRedBlackCartesian(GridBase *base) : GridRedBlackCartesian(base->_fdimensions,base->_simd_layout,base->_processors)  {};
-
-    GridRedBlackCartesian(std::vector<int> &dimensions,
-			  std::vector<int> &simd_layout,
-			  std::vector<int> &processor_grid ) : GridBase(processor_grid)
-    {
-    ///////////////////////
-    // Grid information
-    ///////////////////////
-        _ndimension = dimensions.size();
-        
-        _fdimensions.resize(_ndimension);
-        _gdimensions.resize(_ndimension);
-        _ldimensions.resize(_ndimension);
-        _rdimensions.resize(_ndimension);
-        _simd_layout.resize(_ndimension);
-        
-        _ostride.resize(_ndimension);
-        _istride.resize(_ndimension);
-        
-        _fsites = _gsites = _osites = _isites = 1;
-
-        for(int d=0;d<_ndimension;d++){
-            _fdimensions[d] = dimensions[d];
-            _gdimensions[d] = _fdimensions[d];
-	    _fsites = _fsites * _fdimensions[d];
-	    _gsites = _gsites * _gdimensions[d];
-                
-            if (d==0) _gdimensions[0] = _gdimensions[0]/2; // Remove a checkerboard
-            _ldimensions[d] = _gdimensions[d]/_processors[d];
-
-            // Use a reduced simd grid
-            _simd_layout[d] = simd_layout[d];
-            _rdimensions[d]= _ldimensions[d]/_simd_layout[d];
-
-            _osites *= _rdimensions[d];
-            _isites *= _simd_layout[d];
-                
-            // Addressing support
-            if ( d==0 ) {
-                _ostride[d] = 1;
-                _istride[d] = 1;
-            } else {
-                _ostride[d] = _ostride[d-1]*_rdimensions[d-1];
-                _istride[d] = _istride[d-1]*_simd_layout[d-1];
-            }
-        }
-            
-        ////////////////////////////////////////////////////////////////////////////////////////////
-        // subplane information
-        ////////////////////////////////////////////////////////////////////////////////////////////
-        _slice_block.resize(_ndimension);
-        _slice_stride.resize(_ndimension);
-        _slice_nblock.resize(_ndimension);
-        
-        int block =1;
-        int nblock=1;
-        for(int d=0;d<_ndimension;d++) nblock*=_rdimensions[d];
-            
-        for(int d=0;d<_ndimension;d++){
-            nblock/=_rdimensions[d];
-            _slice_block[d] =block;
-            _slice_stride[d]=_ostride[d]*_rdimensions[d];
-            _slice_nblock[d]=nblock;
-            block = block*_rdimensions[d];
-        }
-            
-    };
-protected:
-    virtual int oIndex(std::vector<int> &coor)
-    {
-        int idx=_ostride[0]*((coor[0]/2)%_rdimensions[0]);
-        for(int d=1;d<_ndimension;d++) idx+=_ostride[d]*(coor[d]%_rdimensions[d]);
-        return idx;
-    };
-        
-};
-
-}
-#endif

lib/communicator/Communicator_base.h

@@ -27,7 +27,7 @@ class CartesianCommunicator {
 #endif
 
     // Constructor
-    CartesianCommunicator(std::vector<int> &pdimensions_in);
+    CartesianCommunicator(const std::vector<int> &pdimensions_in);
 
     // Wraps MPI_Cart routines
     void ShiftedRanks(int dim,int shift,int & source, int & dest);

lib/communicator/Communicator_mpi.cc

@@ -5,7 +5,7 @@ namespace Grid {
 
   // Should error check all MPI calls.
 
-CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
+CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors)
 {
   _ndimension = processors.size();
   std::vector<int> periodic(_ndimension,1);

lib/communicator/Communicator_none.cc

@@ -1,7 +1,7 @@
 #include "Grid.h"
 namespace Grid {
 
-CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
+CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors)
 {
   _processors = processors;
   _ndimension = processors.size();

lib/cshift/Cshift_common.h

@@ -153,7 +153,7 @@ PARALLEL_NESTED_LOOP2
 //////////////////////////////////////////////////////
 // local to node block strided copies
 //////////////////////////////////////////////////////
-template<class vobj> void Copy_plane(Lattice<vobj>& lhs,Lattice<vobj> &rhs, int dimension,int lplane,int rplane,int cbmask)
+template<class vobj> void Copy_plane(Lattice<vobj>& lhs,const Lattice<vobj> &rhs, int dimension,int lplane,int rplane,int cbmask)
 {
   int rd = rhs._grid->_rdimensions[dimension];
 
@@ -175,7 +175,6 @@ PARALLEL_NESTED_LOOP2
       if ( ocb&cbmask ) {
 	//lhs._odata[lo+o]=rhs._odata[ro+o];
 	vstream(lhs._odata[lo+o],rhs._odata[ro+o]);
-	
       }
 
     }
@@ -183,7 +182,7 @@ PARALLEL_NESTED_LOOP2
   
 }
 
-template<class vobj> void Copy_plane_permute(Lattice<vobj>& lhs,Lattice<vobj> &rhs, int dimension,int lplane,int rplane,int cbmask,int permute_type)
+template<class vobj> void Copy_plane_permute(Lattice<vobj>& lhs,const Lattice<vobj> &rhs, int dimension,int lplane,int rplane,int cbmask,int permute_type)
 {
  
   int rd = rhs._grid->_rdimensions[dimension];
@@ -213,12 +212,12 @@ PARALLEL_NESTED_LOOP2
 //////////////////////////////////////////////////////
 // Local to node Cshift
 //////////////////////////////////////////////////////
-template<class vobj> void Cshift_local(Lattice<vobj>& ret,Lattice<vobj> &rhs,int dimension,int shift)
+template<class vobj> void Cshift_local(Lattice<vobj>& ret,const Lattice<vobj> &rhs,int dimension,int shift)
 {
   int sshift[2];
 
-  sshift[0] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,0);
+  sshift[0] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Even);
-  sshift[1] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,1);
+  sshift[1] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Odd);
 
   if ( sshift[0] == sshift[1] ) {
     Cshift_local(ret,rhs,dimension,shift,0x3);
@@ -228,7 +227,7 @@ template<class vobj> void Cshift_local(Lattice<vobj>& ret,Lattice<vobj> &rhs,int
   }
 }
 
-template<class vobj> Lattice<vobj> Cshift_local(Lattice<vobj> &ret,Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
+template<class vobj> Lattice<vobj> Cshift_local(Lattice<vobj> &ret,const Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
 {
   GridBase *grid = rhs._grid;
   int fd = grid->_fdimensions[dimension];
@@ -239,8 +238,7 @@ template<class vobj> Lattice<vobj> Cshift_local(Lattice<vobj> &ret,Lattice<vobj>
   // Map to always positive shift modulo global full dimension.
   shift = (shift+fd)%fd;
 
-  ret.checkerboard = grid->CheckerBoardDestination(rhs.checkerboard,shift);
+  ret.checkerboard = grid->CheckerBoardDestination(rhs.checkerboard,shift,dimension);
-        
   // the permute type
   int permute_dim =grid->PermuteDim(dimension);
   int permute_type=grid->PermuteType(dimension);
@@ -250,11 +248,11 @@ template<class vobj> Lattice<vobj> Cshift_local(Lattice<vobj> &ret,Lattice<vobj>
     int o   = 0;
     int bo  = x * grid->_ostride[dimension];
     
-    int cb= (cbmask==0x2)? 1 : 0;
+    int cb= (cbmask==0x2)? Odd : Even;
 
-    int sshift = grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,cb);
+    int sshift = grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,cb);
     int sx     = (x+sshift)%rd;
-	
+
     int permute_slice=0;
     if(permute_dim){
       int wrap = sshift/rd;

lib/cshift/Cshift_mpi.h

@@ -4,7 +4,7 @@
 
 namespace Grid { 
 
-template<class vobj> Lattice<vobj> Cshift(Lattice<vobj> &rhs,int dimension,int shift)
+template<class vobj> Lattice<vobj> Cshift(const Lattice<vobj> &rhs,int dimension,int shift)
 {
   typedef typename vobj::vector_type vector_type;
   typedef typename vobj::scalar_type scalar_type;
@@ -17,7 +17,7 @@ template<class vobj> Lattice<vobj> Cshift(Lattice<vobj> &rhs,int dimension,int s
   // Map to always positive shift modulo global full dimension.
   shift = (shift+fd)%fd;
 
-  ret.checkerboard = rhs._grid->CheckerBoardDestination(rhs.checkerboard,shift);
+  ret.checkerboard = rhs._grid->CheckerBoardDestination(rhs.checkerboard,shift,dimension);
         
   // the permute type
   int simd_layout     = rhs._grid->_simd_layout[dimension];
@@ -35,12 +35,12 @@ template<class vobj> Lattice<vobj> Cshift(Lattice<vobj> &rhs,int dimension,int s
   return ret;
 }
 
-template<class vobj> void Cshift_comms(Lattice<vobj>& ret,Lattice<vobj> &rhs,int dimension,int shift)
+template<class vobj> void Cshift_comms(Lattice<vobj>& ret,const Lattice<vobj> &rhs,int dimension,int shift)
 {
   int sshift[2];
 
-  sshift[0] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,0);
+  sshift[0] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Even);
-  sshift[1] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,1);
+  sshift[1] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Odd);
 
   if ( sshift[0] == sshift[1] ) {
     Cshift_comms(ret,rhs,dimension,shift,0x3);
@@ -50,12 +50,12 @@ template<class vobj> void Cshift_comms(Lattice<vobj>& ret,Lattice<vobj> &rhs,int
   }
 }
 
-template<class vobj> void Cshift_comms_simd(Lattice<vobj>& ret,Lattice<vobj> &rhs,int dimension,int shift)
+template<class vobj> void Cshift_comms_simd(Lattice<vobj>& ret,const Lattice<vobj> &rhs,int dimension,int shift)
 {
   int sshift[2];
 
-  sshift[0] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,0);
+  sshift[0] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Even);
-  sshift[1] = rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,1);
+  sshift[1] = rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,Odd);
 
   if ( sshift[0] == sshift[1] ) {
     Cshift_comms_simd(ret,rhs,dimension,shift,0x3);
@@ -65,7 +65,7 @@ template<class vobj> void Cshift_comms_simd(Lattice<vobj>& ret,Lattice<vobj> &rh
   }
 }
 
-template<class vobj> void Cshift_comms(Lattice<vobj> &ret,Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
+template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
 {
   typedef typename vobj::vector_type vector_type;
   typedef typename vobj::scalar_type scalar_type;
@@ -87,8 +87,8 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,Lattice<vobj> &rhs,int
   std::vector<vobj,alignedAllocator<vobj> > send_buf(buffer_size);
   std::vector<vobj,alignedAllocator<vobj> > recv_buf(buffer_size);
 
-  int cb= (cbmask==0x2)? 1 : 0;
+  int cb= (cbmask==0x2)? Odd : Even;
-  int sshift= rhs._grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,cb);
+  int sshift= rhs._grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,cb);
 
   for(int x=0;x<rd;x++){       
 
@@ -124,7 +124,7 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,Lattice<vobj> &rhs,int
   }
 }
 
-template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
+template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
 {
   GridBase *grid=rhs._grid;
   const int Nsimd = grid->Nsimd();
@@ -162,8 +162,8 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,Lattice<vobj> &r
   ///////////////////////////////////////////
   // Work out what to send where
   ///////////////////////////////////////////
-  int cb    = (cbmask==0x2)? 1 : 0;
+  int cb    = (cbmask==0x2)? Odd : Even;
-  int sshift= grid->CheckerBoardShift(rhs.checkerboard,dimension,shift,cb);
+  int sshift= grid->CheckerBoardShiftForCB(rhs.checkerboard,dimension,shift,cb);
 
   // loop over outer coord planes orthog to dim
   for(int x=0;x<rd;x++){       

lib/cshift/Cshift_none.h

@@ -1,10 +1,10 @@
 #ifndef _GRID_CSHIFT_NONE_H_
 #define _GRID_CSHIFT_NONE_H_
 namespace Grid {
-template<class vobj> Lattice<vobj> Cshift(Lattice<vobj> &rhs,int dimension,int shift)
+template<class vobj> Lattice<vobj> Cshift(const Lattice<vobj> &rhs,int dimension,int shift)
 {
   Lattice<vobj> ret(rhs._grid);
-  ret.checkerboard = rhs._grid->CheckerBoardDestination(rhs.checkerboard,shift);
+  ret.checkerboard = rhs._grid->CheckerBoardDestination(rhs.checkerboard,shift,dimension);
   Cshift_local(ret,rhs,dimension,shift);
   return ret;
 }

lib/lattice/Lattice_base.h

@@ -176,7 +176,7 @@ PARALLEL_FOR_LOOP
 PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
 #ifdef STREAMING_STORES
-      vobj tmp = eval(tmp,ss,expr);
+      vobj tmp = eval(ss,expr);
       vstream(_odata[ss] ,tmp);
 #else
       _odata[ss]=eval(ss,expr);
@@ -283,24 +283,23 @@ PARALLEL_FOR_LOOP
 
 
 
-#include <lattice/Grid_lattice_conformable.h>
+#include <lattice/Lattice_conformable.h>
 #define GRID_LATTICE_EXPRESSION_TEMPLATES
 #ifdef  GRID_LATTICE_EXPRESSION_TEMPLATES
-#include <lattice/Grid_lattice_ET.h>
+#include <lattice/Lattice_ET.h>
 #else 
-#include <lattice/Grid_lattice_overload.h>
+#include <lattice/Lattice_overload.h>
 #endif
-#include <lattice/Grid_lattice_arith.h>
+#include <lattice/Lattice_arith.h>
-#include <lattice/Grid_lattice_trace.h>
+#include <lattice/Lattice_trace.h>
-#include <lattice/Grid_lattice_transpose.h>
+#include <lattice/Lattice_transpose.h>
-#include <lattice/Grid_lattice_local.h>
+#include <lattice/Lattice_local.h>
-#include <lattice/Grid_lattice_reduction.h>
+#include <lattice/Lattice_reduction.h>
-#include <lattice/Grid_lattice_peekpoke.h>
+#include <lattice/Lattice_peekpoke.h>
-#include <lattice/Grid_lattice_reality.h>
+#include <lattice/Lattice_reality.h>
-#include <Grid_extract.h>
+#include <lattice/Lattice_coordinate.h>
-#include <lattice/Grid_lattice_coordinate.h>
+#include <lattice/Lattice_rng.h>
-#include <lattice/Grid_lattice_rng.h>
+#include <lattice/Lattice_transfer.h>
-#include <lattice/Grid_lattice_transfer.h>
 
 
 

lib/lattice/Lattice_peekpoke.h

@@ -116,7 +116,7 @@ PARALLEL_FOR_LOOP
 
       int Nsimd = grid->Nsimd();
 
-      assert( l.checkerboard== l._grid->CheckerBoard(site));
+      assert( l.checkerboard == l._grid->CheckerBoard(site));
       assert( sizeof(sobj)*Nsimd == sizeof(vobj));
 
       int rank,odx,idx;

lib/lattice/Lattice_transfer.h

@@ -32,7 +32,6 @@ PARALLEL_FOR_LOOP
       cbos=half._grid->CheckerBoard(coor);
       
       if (cbos==cb) {
-	
 	half._odata[ssh] = full._odata[ss];
 	ssh++;
       }
@@ -45,7 +44,7 @@ PARALLEL_FOR_LOOP
     for(int ss=0;ss<full._grid->oSites();ss++){
       std::vector<int> coor;
       int cbos;
-      
+
       full._grid->oCoorFromOindex(coor,ss);
       cbos=half._grid->CheckerBoard(coor);
       

lib/math/Grid_math_arith.h

@@ -1,11 +0,0 @@
-#ifndef GRID_MATH_ARITH_H
-#define GRID_MATH_ARITH_H
-
-#include <math/Grid_math_arith_add.h>
-#include <math/Grid_math_arith_sub.h>
-#include <math/Grid_math_arith_mac.h>
-#include <math/Grid_math_arith_mul.h>
-#include <math/Grid_math_arith_scalar.h>
-
-#endif
-

lib/qcd/Grid_qcd_wilson_dop.cc

@@ -1,217 +0,0 @@
-#include <Grid.h>
-
-namespace Grid {
-namespace QCD {
-
-const std::vector<int> WilsonMatrix::directions   ({0,1,2,3, 0, 1, 2, 3});
-const std::vector<int> WilsonMatrix::displacements({1,1,1,1,-1,-1,-1,-1});
-
-  int WilsonMatrix::HandOptDslash;
-
-  class WilsonCompressor {
-  public:
-    int mu;
-    int dag;
-
-    WilsonCompressor(int _dag){
-      mu=0;
-      dag=_dag;
-      assert((dag==0)||(dag==1));
-    }
-    void Point(int p) { 
-      mu=p;
-    };
-
-    vHalfSpinColourVector operator () (const vSpinColourVector &in)
-    {
-      vHalfSpinColourVector ret;
-      int mudag=mu;
-      if (dag) {
-	mudag=(mu+Nd)%(2*Nd);
-      }
-      switch(mudag) {
-      case Xp:
-	spProjXp(ret,in);
-	break;
-      case Yp:
-	spProjYp(ret,in);
-	break;
-      case Zp:
-	spProjZp(ret,in);
-	break;
-      case Tp:
-	spProjTp(ret,in);
-	break;
-      case Xm:
-	spProjXm(ret,in);
-	break;
-      case Ym:
-	spProjYm(ret,in);
-	break;
-      case Zm:
-	spProjZm(ret,in);
-	break;
-      case Tm:
-	spProjTm(ret,in);
-	break;
-      default: 
-	assert(0);
-	break;
-      }
-      return ret;
-    }
-  };
-
-  WilsonMatrix::WilsonMatrix(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid, double _mass)  : 
-    CheckerBoardedSparseMatrixBase<LatticeFermion>(&Fgrid,&Hgrid),
-    Stencil    (  _grid,npoint,Even,directions,displacements),
-    StencilEven(_cbgrid,npoint,Even,directions,displacements), // source is Even
-    StencilOdd (_cbgrid,npoint,Odd ,directions,displacements), // source is Odd
-    mass(_mass),
-    Umu(_grid),
-    UmuEven(_cbgrid),
-    UmuOdd (_cbgrid)
-  {
-    // Allocate the required comms buffer
-    comm_buf.resize(Stencil._unified_buffer_size); // this is always big enough to contain EO
-    
-    DoubleStore(Umu,_Umu);
-    pickCheckerboard(Even,UmuEven,Umu);
-    pickCheckerboard(Odd ,UmuOdd,Umu);
-  }
-      
-void WilsonMatrix::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu)
-{
-  LatticeColourMatrix U(_grid);
-
-  for(int mu=0;mu<Nd;mu++){
-    U = peekIndex<LorentzIndex>(Umu,mu);
-    pokeIndex<LorentzIndex>(Uds,U,mu);
-    U = adj(Cshift(U,mu,-1));
-    pokeIndex<LorentzIndex>(Uds,U,mu+4);
-  }
-}
-
-RealD WilsonMatrix::M(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard=in.checkerboard;
-  Dhop(in,out,DaggerNo);
-  out = (4+mass)*in - 0.5*out  ; // FIXME : axpby_norm! fusion fun
-  return norm2(out);
-}
-RealD WilsonMatrix::Mdag(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard=in.checkerboard;
-  Dhop(in,out,DaggerYes);
-  out = (4+mass)*in - 0.5*out  ; // FIXME : axpby_norm! fusion fun
-  return norm2(out);
-}
-
-void WilsonMatrix::Meooe(const LatticeFermion &in, LatticeFermion &out)
-{
-  if ( in.checkerboard == Odd ) {
-    DhopEO(in,out,DaggerNo);
-  } else {
-    DhopOE(in,out,DaggerNo);
-  }
-  out = (-0.5)*out; // FIXME : scale factor in Dhop
-}
-void WilsonMatrix::MeooeDag(const LatticeFermion &in, LatticeFermion &out)
-{
-  if ( in.checkerboard == Odd ) {
-    DhopEO(in,out,DaggerYes);
-  } else {
-    DhopOE(in,out,DaggerYes);
-  }
-  out = (-0.5)*out; // FIXME : scale factor in Dhop
-}
-void WilsonMatrix::Mooee(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard = in.checkerboard;
-  out = (4.0+mass)*in;
-  return ;
-}
-void WilsonMatrix::MooeeDag(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard = in.checkerboard;
-  Mooee(in,out);
-}
-void WilsonMatrix::MooeeInv(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard = in.checkerboard;
-  out = (1.0/(4.0+mass))*in;
-  return ;
-}
-void WilsonMatrix::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard = in.checkerboard;
-  MooeeInv(in,out);
-}
-
-void WilsonMatrix::DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField & U,
-				const LatticeFermion &in, LatticeFermion &out,int dag)
-{
-  assert((dag==DaggerNo) ||(dag==DaggerYes));
-  WilsonCompressor compressor(dag);
-  st.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
-
-  if ( dag == DaggerYes ) {
-    if( HandOptDslash ) {
-PARALLEL_FOR_LOOP
-      for(int sss=0;sss<in._grid->oSites();sss++){
-        DiracOptHand::DhopSiteDag(st,U,comm_buf,sss,in,out);
-      }
-    } else { 
-PARALLEL_FOR_LOOP
-      for(int sss=0;sss<in._grid->oSites();sss++){
-        DiracOpt::DhopSiteDag(st,U,comm_buf,sss,in,out);
-      }
-    }
-  } else {
-    if( HandOptDslash ) {
-PARALLEL_FOR_LOOP
-      for(int sss=0;sss<in._grid->oSites();sss++){
-        DiracOptHand::DhopSite(st,U,comm_buf,sss,in,out);
-      }
-    } else { 
-PARALLEL_FOR_LOOP
-      for(int sss=0;sss<in._grid->oSites();sss++){
-        DiracOpt::DhopSite(st,U,comm_buf,sss,in,out);
-      }
-    }
-  }
-}
-void WilsonMatrix::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag)
-{
-  conformable(in._grid,_cbgrid);    // verifies half grid
-  conformable(in._grid,out._grid); // drops the cb check
-
-  assert(in.checkerboard==Even);
-  out.checkerboard = Odd;
-
-  DhopInternal(StencilEven,UmuOdd,in,out,dag);
-}
-void WilsonMatrix::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag)
-{
-  conformable(in._grid,_cbgrid);    // verifies half grid
-  conformable(in._grid,out._grid); // drops the cb check
-
-  assert(in.checkerboard==Odd);
-  out.checkerboard = Even;
-
-  DhopInternal(StencilOdd,UmuEven,in,out,dag);
-}
-void WilsonMatrix::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag)
-{
-  conformable(in._grid,_grid); // verifies full grid
-  conformable(in._grid,out._grid);
-
-  out.checkerboard = in.checkerboard;
-
-  DhopInternal(Stencil,Umu,in,out,dag);
-}
-
-}}
-
-
-

lib/qcd/Grid_qcd_wilson_dop.h

@@ -1,105 +0,0 @@
-#ifndef  GRID_QCD_WILSON_DOP_H
-#define  GRID_QCD_WILSON_DOP_H
-
-
-namespace Grid {
-
-  namespace QCD {
-
-  // Should be in header?
-    const int Xp = 0;
-    const int Yp = 1;
-    const int Zp = 2;
-    const int Tp = 3;
-    const int Xm = 4;
-    const int Ym = 5;
-    const int Zm = 6;
-    const int Tm = 7;
-
-    class WilsonMatrix : public CheckerBoardedSparseMatrixBase<LatticeFermion>
-    {
-      //NB r=1;
-    public:
-      static int HandOptDslash;
-
-      double                        mass;
-      //      GridBase                     *    grid; // Inherited
-      //      GridBase                     *  cbgrid;
-
-      //Defines the stencils for even and odd
-      CartesianStencil Stencil; 
-      CartesianStencil StencilEven; 
-      CartesianStencil StencilOdd; 
-
-      // Copy of the gauge field , with even and odd subsets
-      LatticeDoubledGaugeField Umu;
-      LatticeDoubledGaugeField UmuEven;
-      LatticeDoubledGaugeField UmuOdd;
-
-      static const int npoint=8;
-      static const std::vector<int> directions   ;
-      static const std::vector<int> displacements;
-      static const int Xp,Xm,Yp,Ym,Zp,Zm,Tp,Tm;
-
-      // Comms buffer
-      std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  comm_buf;
-
-      // Constructor
-      WilsonMatrix(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,double _mass);
-
-      // DoubleStore
-      void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
-
-      // override multiply
-      virtual RealD  M    (const LatticeFermion &in, LatticeFermion &out);
-      virtual RealD  Mdag (const LatticeFermion &in, LatticeFermion &out);
-
-      // half checkerboard operaions
-      virtual void   Meooe       (const LatticeFermion &in, LatticeFermion &out);
-      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out);
-      virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out);
-      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out);
-      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
-      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
-
-      // non-hermitian hopping term; half cb or both
-      void Dhop  (const LatticeFermion &in, LatticeFermion &out,int dag);
-      void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag);
-      void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag);
-      void DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField &U,
-			const LatticeFermion &in, LatticeFermion &out,int dag);
-
-      typedef iScalar<iMatrix<vComplex, Nc> > matrix;
-
-      
-    };
-
-
-    class DiracOpt {
-    public:
-      // These ones will need to be package intelligently. WilsonType base class
-      // for use by DWF etc..
-      static void DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
-		    std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
-		    int ss,const LatticeFermion &in, LatticeFermion &out);
-      static void DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
-		       std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
-		       int ss,const LatticeFermion &in, LatticeFermion &out);
-
-    };
-    class DiracOptHand {
-    public:
-      // These ones will need to be package intelligently. WilsonType base class
-      // for use by DWF etc..
-      static void DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
-		    std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
-		    int ss,const LatticeFermion &in, LatticeFermion &out);
-      static void DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
-		       std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
-		       int ss,const LatticeFermion &in, LatticeFermion &out);
-
-    };
-
-  }
-}
-#endif

lib/qcd/LinalgUtils.h

@@ -0,0 +1,113 @@
+#ifndef GRID_QCD_LINALG_UTILS_H
+#define GRID_QCD_LINALG_UTILS_H
+
+namespace Grid{
+namespace QCD{
+////////////////////////////////////////////////////////////////////////
+//This file brings additional linear combination assist that is helpful
+//to QCD such as chiral projectors and spin matrices applied to one of the inputs.
+//These routines support five-D chiral fermions and contain s-subslice indexing 
+//on the 5d (rb4d) checkerboarded lattices
+////////////////////////////////////////////////////////////////////////
+template<class vobj> 
+void axpby_ssp(Lattice<vobj> &z, RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp = a*x._odata[ss+s]+b*y._odata[ss+sp];
+    vstream(z._odata[ss+s],tmp);
+  }
+}
+
+template<class vobj> 
+void ag5xpby_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp;
+    multGamma5(tmp(),a*x._odata[ss+s]());
+    tmp = tmp + b*y._odata[ss+sp];
+    vstream(z._odata[ss+s],tmp);
+  }
+}
+
+template<class vobj> 
+void axpbg5y_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp;
+    multGamma5(tmp(),b*y._odata[ss+sp]());
+    tmp = tmp + a*x._odata[ss+s];
+    vstream(z._odata[ss+s],tmp);
+  }
+}
+
+template<class vobj> 
+void ag5xpbg5y_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp1;
+    vobj tmp2;
+    tmp1 = a*x._odata[ss+s]+b*y._odata[ss+sp];
+    multGamma5(tmp2(),tmp1());
+    vstream(z._odata[ss+s],tmp2);
+  }
+}
+
+template<class vobj> 
+void axpby_ssp_pminus(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp;
+    spProj5m(tmp,y._odata[ss+sp]);
+    tmp = a*x._odata[ss+s]+b*tmp;
+    vstream(z._odata[ss+s],tmp);
+  }
+}
+
+template<class vobj> 
+void axpby_ssp_pplus(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp;
+    spProj5p(tmp,y._odata[ss+sp]);
+    tmp = a*x._odata[ss+s]+b*tmp;
+    vstream(z._odata[ss+s],tmp);
+  }
+}
+}}
+#endif 

lib/qcd/QCD.h

@@ -4,6 +4,16 @@ namespace Grid{
 
 namespace QCD {
 
+
+    static const int Xp = 0;
+    static const int Yp = 1;
+    static const int Zp = 2;
+    static const int Tp = 3;
+    static const int Xm = 4;
+    static const int Ym = 5;
+    static const int Zm = 6;
+    static const int Tm = 7;
+
     static const int Nc=3;
     static const int Ns=4;
     static const int Nd=4;
@@ -297,9 +307,10 @@ namespace QCD {
 }   //namespace QCD
 } // Grid
 
-#include <qcd/Grid_qcd_dirac.h>
+#include <qcd/SpaceTimeGrid.h>
-#include <qcd/Grid_qcd_2spinor.h>
+#include <qcd/Dirac.h>
-//#include <qcd/Grid_qcd_pauli.h>
+#include <qcd/TwoSpinor.h>
-#include <qcd/Grid_qcd_wilson_dop.h>
+#include <qcd/LinalgUtils.h>
+#include <qcd/action/Actions.h>
 
 #endif

lib/qcd/SpaceTimeGrid.cc

@@ -0,0 +1,52 @@
+#include <Grid.h>
+
+namespace Grid { 
+  namespace QCD {
+
+/////////////////////////////////////////////////////////////////
+// Public interface
+/////////////////////////////////////////////////////////////////
+GridCartesian *SpaceTimeGrid::makeFourDimGrid(const std::vector<int> & latt,const std::vector<int> &simd,const std::vector<int> &mpi)
+{
+  return new GridCartesian(latt,simd,mpi); 
+}
+GridRedBlackCartesian *SpaceTimeGrid::makeFourDimRedBlackGrid(const GridCartesian *FourDimGrid)
+{
+  return new GridRedBlackCartesian(FourDimGrid); 
+}
+
+GridCartesian         *SpaceTimeGrid::makeFiveDimGrid(int Ls,const GridCartesian *FourDimGrid)
+{
+  int N4=FourDimGrid->_ndimension;
+
+  std::vector<int> latt5(1,Ls);
+  std::vector<int> simd5(1,1);
+  std::vector<int>  mpi5(1,1);
+  
+  for(int d=0;d<N4;d++){
+    latt5.push_back(FourDimGrid->_fdimensions[d]);
+    simd5.push_back(FourDimGrid->_simd_layout[d]);
+     mpi5.push_back(FourDimGrid->_processors[d]);
+  }
+  return new GridCartesian(latt5,simd5,mpi5); 
+}
+
+GridRedBlackCartesian *SpaceTimeGrid::makeFiveDimRedBlackGrid(int Ls,const GridCartesian *FourDimGrid)
+{
+  int N4=FourDimGrid->_ndimension;
+  int cbd=1;
+  std::vector<int> latt5(1,Ls);
+  std::vector<int> simd5(1,1);
+  std::vector<int>  mpi5(1,1);
+  std::vector<int>   cb5(1,0);
+    
+  for(int d=0;d<N4;d++){
+    latt5.push_back(FourDimGrid->_fdimensions[d]);
+    simd5.push_back(FourDimGrid->_simd_layout[d]);
+     mpi5.push_back(FourDimGrid->_processors[d]);
+      cb5.push_back(  1);
+    }
+  return new GridRedBlackCartesian(latt5,simd5,mpi5,cb5,cbd); 
+}
+
+}}

lib/qcd/SpaceTimeGrid.h

@@ -0,0 +1,18 @@
+#ifndef GRID_QCD_SPACE_TIME_GRID_H
+#define GRID_QCD_SPACE_TIME_GRID_H
+namespace Grid {
+namespace QCD {
+
+class SpaceTimeGrid {
+ public:
+
+  static GridCartesian         *makeFourDimGrid(const std::vector<int> & latt,const std::vector<int> &simd,const std::vector<int> &mpi);
+  static GridRedBlackCartesian *makeFourDimRedBlackGrid       (const GridCartesian *FourDimGrid);
+  static GridCartesian         *makeFiveDimGrid        (int Ls,const GridCartesian *FourDimGrid);
+  static GridRedBlackCartesian *makeFiveDimRedBlackGrid(int Ls,const GridCartesian *FourDimGrid);
+
+};
+
+}}
+
+#endif

lib/qcd/action/Actions.h

@@ -0,0 +1,98 @@
+#ifndef GRID_QCD_ACTIONS_H
+#define GRID_QCD_ACTIONS_H
+
+
+// Some reorganisation likely required as both Chroma and IroIro
+// are separating the concept of the operator from that of action.
+//
+// The FermAction contains methods to create 
+//
+// * Linear operators             (Hermitian and non-hermitian)  .. my LinearOperator
+// * System solvers               (Hermitian and non-hermitian)  .. my OperatorFunction
+// * MultiShift System solvers    (Hermitian and non-hermitian)  .. my OperatorFunction
+
+
+////////////////////////////////////////////
+// Abstract base interface
+////////////////////////////////////////////
+#include <qcd/action/fermion/FermionOperator.h>
+
+////////////////////////////////////////////
+// Utility functions
+////////////////////////////////////////////
+#include <qcd/action/fermion/WilsonCompressor.h>     //used by all wilson type fermions
+#include <qcd/action/fermion/WilsonKernels.h>        //used by all wilson type fermions
+
+////////////////////////////////////////////
+// 4D formulations
+////////////////////////////////////////////
+#include <qcd/action/fermion/WilsonFermion.h>
+//#include <qcd/action/fermion/CloverFermion.h>
+
+////////////////////////////////////////////
+// 5D formulations...
+////////////////////////////////////////////
+
+#include <qcd/action/fermion/WilsonFermion5D.h> // used by all 5d overlap types
+
+//////////
+// Cayley
+//////////
+#include <qcd/action/fermion/CayleyFermion5D.h>
+
+#include <qcd/action/fermion/DomainWallFermion.h>
+#include <qcd/action/fermion/DomainWallFermion.h>
+
+#include <qcd/action/fermion/MobiusFermion.h>
+#include <qcd/action/fermion/ScaledShamirFermion.h>
+#include <qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h>
+
+#include <qcd/action/fermion/MobiusZolotarevFermion.h>
+#include <qcd/action/fermion/ShamirZolotarevFermion.h>
+#include <qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h>
+
+//////////////////////
+// Continued fraction
+//////////////////////
+#include <qcd/action/fermion/ContinuedFractionFermion5D.h>
+#include <qcd/action/fermion/OverlapWilsonContfracTanhFermion.h>
+#include <qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h>
+
+//////////////////////
+// Partial fraction
+//////////////////////
+#include <qcd/action/fermion/PartialFractionFermion5D.h>
+
+
+    // Chroma interface defining FermionAction
+    /*
+     template<typename T, typename P, typename Q>  class FermAct4D : public FermionAction<T,P,Q>
+     virtual LinearOperator<T>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
+     virtual LinearOperator<T>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
+     virtual LinOpSystemSolver<T>* invLinOp(Handle< FermState<T,P,Q> > state,
+     virtual MdagMSystemSolver<T>* invMdagM(Handle< FermState<T,P,Q> > state,
+     virtual LinOpMultiSystemSolver<T>* mInvLinOp(Handle< FermState<T,P,Q> > state,
+     virtual MdagMMultiSystemSolver<T>* mInvMdagM(Handle< FermState<T,P,Q> > state,
+     virtual MdagMMultiSystemSolverAccumulate<T>* mInvMdagMAcc(Handle< FermState<T,P,Q> > state,
+     virtual SystemSolver<T>* qprop(Handle< FermState<T,P,Q> > state,
+     class DiffFermAct4D : public FermAct4D<T,P,Q>
+     virtual DiffLinearOperator<T,Q,P>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
+     virtual DiffLinearOperator<T,Q,P>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
+    */
+
+
+    // Chroma interface defining GaugeAction
+    /*
+      template<typename P, typename Q>   class GaugeAction
+  virtual const CreateGaugeState<P,Q>& getCreateState() const = 0;
+  virtual GaugeState<P,Q>* createState(const Q& q) const
+  virtual const GaugeBC<P,Q>& getGaugeBC() const
+  virtual const Set& getSet(void) const = 0;
+  virtual void deriv(P& result, const Handle< GaugeState<P,Q> >& state) const 
+  virtual Double S(const Handle< GaugeState<P,Q> >& state) const = 0;
+
+  class LinearGaugeAction : public GaugeAction< multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
+  typedef multi1d<LatticeColorMatrix>  P;
+    */
+
+#endif

lib/qcd/action/fermion/CayleyFermion5D.cc

@@ -0,0 +1,346 @@
+#include <Grid.h>
+namespace Grid {
+namespace QCD {
+
+ CayleyFermion5D::CayleyFermion5D(LatticeGaugeField &_Umu,
+				  GridCartesian         &FiveDimGrid,
+				  GridRedBlackCartesian &FiveDimRedBlackGrid,
+				  GridCartesian         &FourDimGrid,
+				  GridRedBlackCartesian &FourDimRedBlackGrid,
+				  RealD _mass,RealD _M5) :
+   WilsonFermion5D(_Umu,
+		   FiveDimGrid,
+		   FiveDimRedBlackGrid,
+		   FourDimGrid,
+		   FourDimRedBlackGrid,_M5),
+   mass(_mass)
+ {
+ }
+
+  // override multiply
+  RealD CayleyFermion5D::M    (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    LatticeFermion Din(psi._grid);
+
+    // Assemble Din
+    for(int s=0;s<Ls;s++){
+      if ( s==0 ) {
+	//	Din = bs psi[s] + cs[s] psi[s+1}
+	axpby_ssp_pminus(Din,bs[s],psi,cs[s],psi,s,s+1);
+	//      Din+= -mass*cs[s] psi[s+1}
+	axpby_ssp_pplus (Din,1.0,Din,-mass*cs[s],psi,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pminus(Din,bs[s],psi,-mass*cs[s],psi,s,0);
+	axpby_ssp_pplus (Din,1.0,Din,cs[s],psi,s,s-1);
+      } else {
+	axpby_ssp_pminus(Din,bs[s],psi,cs[s],psi,s,s+1);
+	axpby_ssp_pplus(Din,1.0,Din,cs[s],psi,s,s-1);
+      }
+    }
+
+    DW(Din,chi,DaggerNo);
+    // ((b D_W + D_w hop terms +1) on s-diag
+    axpby(chi,1.0,1.0,chi,psi); 
+
+    for(int s=0;s<Ls;s++){
+      if ( s==0 ){
+	axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s+1);
+	axpby_ssp_pplus (chi,1.0,chi,mass,psi,s,Ls-1);
+      } else if ( s==(Ls-1)) {
+	axpby_ssp_pminus(chi,1.0,chi,mass,psi,s,0);
+	axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s-1);
+      } else {
+	axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s+1);
+	axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s-1);
+      }
+    }
+    return norm2(chi);
+  }
+
+  RealD CayleyFermion5D::Mdag (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    // Under adjoint
+    //D1+        D1- P-    ->   D1+^dag   P+ D2-^dag
+    //D2- P+     D2+            P-D1-^dag D2+dag
+
+    LatticeFermion Din(psi._grid);
+    // Apply Dw
+    DW(psi,Din,DaggerYes); 
+
+    for(int s=0;s<Ls;s++){
+      // Collect the terms in DW
+      //	Chi = bs Din[s] + cs[s] Din[s+1}
+      //    Chi+= -mass*cs[s] psi[s+1}
+      if ( s==0 ) {
+	axpby_ssp_pplus (chi,bs[s],Din,cs[s+1],Din,s,s+1);
+	axpby_ssp_pminus(chi,1.0,chi,-mass*cs[Ls-1],Din,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pplus (chi,bs[s],Din,-mass*cs[0],Din,s,0);
+	axpby_ssp_pminus(chi,1.0,chi,cs[s-1],Din,s,s-1);
+      } else {
+	axpby_ssp_pplus (chi,bs[s],Din,cs[s+1],Din,s,s+1);
+	axpby_ssp_pminus(chi,1.0,chi,cs[s-1],Din,s,s-1);
+      }
+      // Collect the terms indept of DW
+      if ( s==0 ){
+	axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s+1);
+	axpby_ssp_pminus(chi,1.0,chi,mass,psi,s,Ls-1);
+      } else if ( s==(Ls-1)) {
+	axpby_ssp_pplus (chi,1.0,chi,mass,psi,s,0);
+	axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s-1);
+      } else {
+	axpby_ssp_pplus(chi,1.0,chi,-1.0,psi,s,s+1);
+	axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s-1);
+      }
+    }
+    // ((b D_W + D_w hop terms +1) on s-diag
+    axpby (chi,1.0,1.0,chi,psi); 
+    return norm2(chi);
+  }
+
+  // half checkerboard operations
+  void CayleyFermion5D::Meooe       (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    LatticeFermion tmp(psi._grid);
+    // Assemble the 5d matrix
+    for(int s=0;s<Ls;s++){
+      if ( s==0 ) {
+	//	tmp = bs psi[s] + cs[s] psi[s+1}
+	//      tmp+= -mass*cs[s] psi[s+1}
+	axpby_ssp_pminus(tmp,beo[s],psi,-ceo[s],psi ,s, s+1);
+	axpby_ssp_pplus(tmp,1.0,tmp,mass*ceo[s],psi,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pminus(tmp,beo[s],psi,mass*ceo[s],psi,s,0);
+	axpby_ssp_pplus(tmp,1.0,tmp,-ceo[s],psi,s,s-1);
+      } else {
+	axpby_ssp_pminus(tmp,beo[s],psi,-ceo[s],psi,s,s+1);
+	axpby_ssp_pplus (tmp,1.0,tmp,-ceo[s],psi,s,s-1);
+      }
+    }
+    // Apply 4d dslash
+    if ( psi.checkerboard == Odd ) {
+      DhopEO(tmp,chi,DaggerNo);
+    } else {
+      DhopOE(tmp,chi,DaggerNo);
+    }
+  }
+
+  void CayleyFermion5D::MeooeDag    (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    LatticeFermion tmp(psi._grid);
+    // Apply 4d dslash
+    if ( psi.checkerboard == Odd ) {
+      DhopEO(psi,tmp,DaggerYes);
+    } else {
+      DhopOE(psi,tmp,DaggerYes);
+    }
+    // Assemble the 5d matrix
+    for(int s=0;s<Ls;s++){
+      if ( s==0 ) {
+	axpby_ssp_pplus(chi,beo[s],tmp,   -ceo[s+1]  ,tmp,s,s+1);
+	axpby_ssp_pminus(chi,   1.0,chi,mass*ceo[Ls-1],tmp,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pplus(chi,beo[s],tmp,mass*ceo[0],tmp,s,0);
+	axpby_ssp_pminus(chi,1.0,chi,-ceo[s-1],tmp,s,s-1);
+      } else {
+	axpby_ssp_pplus(chi,beo[s],tmp,-ceo[s+1],tmp,s,s+1);
+	axpby_ssp_pminus(chi,1.0   ,chi,-ceo[s-1],tmp,s,s-1);
+      }
+    }
+  }
+
+  void CayleyFermion5D::Mooee       (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    for (int s=0;s<Ls;s++){
+      if ( s==0 ) {
+	axpby_ssp_pminus(chi,bee[s],psi ,-cee[s],psi,s,s+1);
+	axpby_ssp_pplus (chi,1.0,chi,mass*cee[s],psi,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pminus(chi,bee[s],psi,mass*cee[s],psi,s,0);
+	axpby_ssp_pplus (chi,1.0,chi,-cee[s],psi,s,s-1);
+      } else {
+	axpby_ssp_pminus(chi,bee[s],psi,-cee[s],psi,s,s+1);
+	axpby_ssp_pplus (chi,1.0,chi,-cee[s],psi,s,s-1);
+      }
+    }
+  }
+
+  void CayleyFermion5D::MooeeDag    (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    for (int s=0;s<Ls;s++){
+      // Assemble the 5d matrix
+      if ( s==0 ) {
+	axpby_ssp_pplus(chi,bee[s],psi,-cee[s+1]  ,psi,s,s+1);
+	axpby_ssp_pminus(chi,1.0,chi,mass*cee[Ls-1],psi,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pplus(chi,bee[s],psi,mass*cee[0],psi,s,0);
+	axpby_ssp_pminus(chi,1.0,chi,-cee[s-1],psi,s,s-1);
+      } else {
+	axpby_ssp_pplus(chi,bee[s],psi,-cee[s+1],psi,s,s+1);
+	axpby_ssp_pminus(chi,1.0   ,chi,-cee[s-1],psi,s,s-1);
+      }
+    }
+  }
+
+  void CayleyFermion5D::MooeeInv    (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    // Apply (L^{\prime})^{-1}
+    axpby_ssp (chi,1.0,psi,     0.0,psi,0,0);      // chi[0]=psi[0]
+    for (int s=1;s<Ls;s++){
+      axpby_ssp_pplus(chi,1.0,psi,-lee[s-1],chi,s,s-1);// recursion Psi[s] -lee P_+ chi[s-1]
+    }
+    // L_m^{-1} 
+    for (int s=0;s<Ls-1;s++){ // Chi[ee] = 1 - sum[s<Ls-1] -leem[s]P_- chi
+      axpby_ssp_pminus(chi,1.0,chi,-leem[s],chi,Ls-1,s);
+    }
+    // U_m^{-1} D^{-1}
+    for (int s=0;s<Ls-1;s++){
+      // Chi[s] + 1/d chi[s] 
+      axpby_ssp_pplus(chi,1.0/dee[s],chi,-ueem[s]/dee[Ls-1],chi,s,Ls-1);
+    }	
+    axpby_ssp(chi,1.0/dee[Ls-1],chi,0.0,chi,Ls-1,Ls-1); // Modest avoidable 
+    
+    // Apply U^{-1}
+    for (int s=Ls-2;s>=0;s--){
+      axpby_ssp_pminus (chi,1.0,chi,-uee[s],chi,s,s+1);  // chi[Ls]
+    }
+  }
+
+  void CayleyFermion5D::MooeeInvDag (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    // Apply (U^{\prime})^{-dagger}
+    axpby_ssp (chi,1.0,psi,     0.0,psi,0,0);      // chi[0]=psi[0]
+    for (int s=1;s<Ls;s++){
+      axpby_ssp_pminus(chi,1.0,psi,-uee[s-1],chi,s,s-1);
+    }
+    // U_m^{-\dagger} 
+    for (int s=0;s<Ls-1;s++){
+      axpby_ssp_pplus(chi,1.0,chi,-ueem[s],chi,Ls-1,s);
+    }
+    // L_m^{-\dagger} D^{-dagger}
+    for (int s=0;s<Ls-1;s++){
+      axpby_ssp_pminus(chi,1.0/dee[s],chi,-leem[s]/dee[Ls-1],chi,s,Ls-1);
+    }	
+    axpby_ssp(chi,1.0/dee[Ls-1],chi,0.0,chi,Ls-1,Ls-1); // Modest avoidable 
+    
+    // Apply L^{-dagger}
+    for (int s=Ls-2;s>=0;s--){
+      axpby_ssp_pplus (chi,1.0,chi,-lee[s],chi,s,s+1);  // chi[Ls]
+    }
+  }
+  
+  // Tanh
+  void CayleyFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c)
+  {
+    SetCoefficientsZolotarev(1.0,zdata,b,c);
+
+  }
+  //Zolo
+  void CayleyFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c)
+  {
+
+    ///////////////////////////////////////////////////////////
+    // The Cayley coeffs (unprec)
+    ///////////////////////////////////////////////////////////
+    omega.resize(Ls);
+    bs.resize(Ls);
+    cs.resize(Ls);
+    as.resize(Ls);
+    
+    // 
+    // Ts = (    [bs+cs]Dw        )^-1 (    (bs+cs) Dw         )
+    //     -(g5  -------       -1 )    ( g5 ---------     + 1  )
+    //      (   {2+(bs-cs)Dw}     )    (    2+(bs-cs) Dw       )
+    //
+    //  bs = 1/2( (1/omega_s + 1)*b + (1/omega - 1)*c ) = 1/2(  1/omega(b+c) + (b-c) )
+    //  cs = 1/2( (1/omega_s - 1)*b + (1/omega + 1)*c ) = 1/2(  1/omega(b+c) - (b-c) )
+    //
+    // bs+cs = 0.5*( 1/omega(b+c) + (b-c) + 1/omega(b+c) - (b-c) ) = 1/omega(b+c)
+    // bs-cs = 0.5*( 1/omega(b+c) + (b-c) - 1/omega(b+c) + (b-c) ) = b-c
+    //
+    // So 
+    //
+    // Ts = (    [b+c]Dw/omega_s    )^-1 (    (b+c) Dw /omega_s        )
+    //     -(g5  -------         -1 )    ( g5 ---------           + 1  )
+    //      (   {2+(b-c)Dw}         )    (    2+(b-c) Dw               )
+    //
+    // Ts = (    [b+c]Dw            )^-1 (    (b+c) Dw                 )
+    //     -(g5  -------    -omega_s)    ( g5 ---------      + omega_s )
+    //      (   {2+(b-c)Dw}         )    (    2+(b-c) Dw               )
+    // 
+    
+    double bpc = b+c;
+    double bmc = b-c;
+    for(int i=0; i < Ls; i++){
+      as[i] = 1.0;
+      omega[i] = ((double)zdata->gamma[i])*zolo_hi; //NB reciprocal relative to Chroma NEF code
+      bs[i] = 0.5*(bpc/omega[i] + bmc);
+      cs[i] = 0.5*(bpc/omega[i] - bmc);
+    }
+
+    ////////////////////////////////////////////////////////
+    // Constants for the preconditioned matrix Cayley form
+    ////////////////////////////////////////////////////////
+    bee.resize(Ls);
+    cee.resize(Ls);
+    beo.resize(Ls);
+    ceo.resize(Ls);
+    
+    for(int i=0;i<Ls;i++){
+      bee[i]=as[i]*(bs[i]*(4.0-M5) +1.0);
+      cee[i]=as[i]*(1.0-cs[i]*(4.0-M5));
+      beo[i]=as[i]*bs[i];
+      ceo[i]=-as[i]*cs[i];
+    }
+
+    aee.resize(Ls);
+    aeo.resize(Ls);
+    for(int i=0;i<Ls;i++){
+      aee[i]=cee[i];
+      aeo[i]=ceo[i];
+    }
+
+    //////////////////////////////////////////
+    // LDU decomposition of eeoo
+    //////////////////////////////////////////
+    dee.resize(Ls);
+    lee.resize(Ls);
+    leem.resize(Ls);
+    uee.resize(Ls);
+    ueem.resize(Ls);
+    
+    for(int i=0;i<Ls;i++){
+      
+      dee[i] = bee[i];
+      
+      if ( i < Ls-1 ) {
+	
+	lee[i] =-cee[i+1]/bee[i]; // sub-diag entry on the ith column
+	    
+	leem[i]=mass*cee[Ls-1]/bee[0];
+	for(int j=0;j<i;j++)  leem[i]*= aee[j]/bee[j+1];
+	
+	uee[i] =-aee[i]/bee[i];   // up-diag entry on the ith row
+	
+	ueem[i]=mass;
+	for(int j=1;j<=i;j++) ueem[i]*= cee[j]/bee[j];
+	ueem[i]*= aee[0]/bee[0];
+	    
+      } else { 
+	lee[i] =0.0;
+	leem[i]=0.0;
+	uee[i] =0.0;
+	ueem[i]=0.0;
+      }
+    }
+	
+    { 
+      double delta_d=mass*cee[Ls-1];
+      for(int j=0;j<Ls-1;j++) delta_d *= cee[j]/bee[j];
+      dee[Ls-1] += delta_d;
+    }
+  }
+
+}}
+
+

lib/qcd/action/fermion/CayleyFermion5D.h

@@ -0,0 +1,62 @@
+#ifndef  GRID_QCD_CAYLEY_FERMION_H
+#define  GRID_QCD_CAYLEY_FERMION_H
+
+namespace Grid {
+
+  namespace QCD {
+
+    class CayleyFermion5D : public WilsonFermion5D
+    {
+    public:
+
+      // override multiply
+      virtual RealD  M    (const LatticeFermion &in, LatticeFermion &out);
+      virtual RealD  Mdag (const LatticeFermion &in, LatticeFermion &out);
+
+      // half checkerboard operations
+      virtual void   Meooe       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   Instantiatable(void)=0;
+      //    protected:
+      RealD mass;
+
+      // Cayley form Moebius (tanh and zolotarev)
+      std::vector<RealD> omega; 
+      std::vector<RealD> bs;    // S dependent coeffs
+      std::vector<RealD> cs;    
+      std::vector<RealD> as;    
+      // For preconditioning Cayley form
+      std::vector<RealD> bee;    
+      std::vector<RealD> cee;    
+      std::vector<RealD> aee;    
+      std::vector<RealD> beo;    
+      std::vector<RealD> ceo;    
+      std::vector<RealD> aeo;    
+      // LDU factorisation of the eeoo matrix
+      std::vector<RealD> lee;    
+      std::vector<RealD> leem;    
+      std::vector<RealD> uee;    
+      std::vector<RealD> ueem;    
+      std::vector<RealD> dee;    
+
+      // Constructors
+      CayleyFermion5D(LatticeGaugeField &_Umu,
+		      GridCartesian         &FiveDimGrid,
+		      GridRedBlackCartesian &FiveDimRedBlackGrid,
+		      GridCartesian         &FourDimGrid,
+		      GridRedBlackCartesian &FourDimRedBlackGrid,
+		      RealD _mass,RealD _M5);
+
+    protected:
+      void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c);
+      void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c);
+    };
+
+  }
+}
+
+#endif

lib/qcd/action/fermion/ContinuedFractionFermion5D.cc

@@ -0,0 +1,171 @@
+#include <Grid.h>
+
+namespace Grid {
+  namespace QCD {
+
+    void ContinuedFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale)
+    {
+      SetCoefficientsZolotarev(1.0/scale,zdata);
+    }
+    void ContinuedFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata)
+    {
+      R=(1+this->mass)/(1-this->mass);
+
+      Beta.resize(Ls);
+      cc.resize(Ls);
+      cc_d.resize(Ls);
+      sqrt_cc.resize(Ls);
+      for(int i=0; i < Ls ; i++){
+	Beta[i] = zdata -> beta[i];
+	cc[i] = 1.0/Beta[i];
+	cc_d[i]=sqrt(cc[i]);
+      }
+    
+      cc_d[Ls-1]=1.0;
+      for(int i=0; i < Ls-1 ; i++){
+	sqrt_cc[i]= sqrt(cc[i]*cc[i+1]);
+      }    
+      sqrt_cc[Ls-2]=sqrt(cc[Ls-2]);
+
+
+      ZoloHiInv =1.0/zolo_hi;
+      double dw_diag = (4.0-M5)*ZoloHiInv;
+    
+      See.resize(Ls);
+      Aee.resize(Ls);
+      int sign=1;
+      for(int s=0;s<Ls;s++){
+	Aee[s] = sign * Beta[s] * dw_diag;
+	sign   = - sign;
+      }
+      Aee[Ls-1] += R;
+    
+      See[0] = Aee[0];
+      for(int s=1;s<Ls;s++){
+	See[s] = Aee[s] - 1.0/See[s-1];
+      }
+      for(int s=0;s<Ls;s++){
+	std::cout <<"s = "<<s<<" Beta "<<Beta[s]<<" Aee "<<Aee[s] <<" See "<<See[s] <<std::endl;
+      }
+    }
+
+
+
+    RealD  ContinuedFractionFermion5D::M           (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      LatticeFermion D(psi._grid);
+
+      DW(psi,D,DaggerNo); 
+
+      int sign=1;
+      for(int s=0;s<Ls;s++){
+	if ( s==0 ) {
+	  ag5xpby_ssp(chi,cc[0]*Beta[0]*sign*ZoloHiInv,D,sqrt_cc[0],psi,s,s+1); // Multiplies Dw by G5 so Hw
+	} else if ( s==(Ls-1) ){
+	  RealD R=(1.0+mass)/(1.0-mass);
+	  ag5xpby_ssp(chi,Beta[s]*ZoloHiInv,D,sqrt_cc[s-1],psi,s,s-1);
+	  ag5xpby_ssp(chi,R,psi,1.0,chi,s,s);
+	} else {
+	  ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*ZoloHiInv,D,sqrt_cc[s],psi,s,s+1);
+  	  axpby_ssp(chi,1.0,chi,sqrt_cc[s-1],psi,s,s-1);
+	}
+	sign=-sign; 
+      }
+      return norm2(chi);
+    }
+    RealD  ContinuedFractionFermion5D::Mdag        (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      // This matrix is already hermitian. (g5 Dw) = Dw dag g5 = (g5 Dw)dag
+      // The rest of matrix is symmetric.
+      // Can ignore "dag"
+      return M(psi,chi);
+    }
+    void   ContinuedFractionFermion5D::Meooe       (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      // Apply 4d dslash
+      if ( psi.checkerboard == Odd ) {
+	DhopEO(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian
+      } else {
+	DhopOE(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian
+      }
+      
+      int sign=1;
+      for(int s=0;s<Ls;s++){
+	if ( s==(Ls-1) ){
+	  ag5xpby_ssp(chi,Beta[s]*ZoloHiInv,chi,0.0,chi,s,s);
+	} else {
+	  ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*ZoloHiInv,chi,0.0,chi,s,s);
+	}
+	sign=-sign; 
+      }
+    }
+    void   ContinuedFractionFermion5D::MeooeDag    (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      Meooe(psi,chi);
+    }
+    void   ContinuedFractionFermion5D::Mooee       (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      double dw_diag = (4.0-M5)*ZoloHiInv;
+    
+      int sign=1;
+      for(int s=0;s<Ls;s++){
+	if ( s==0 ) {
+	  ag5xpby_ssp(chi,cc[0]*Beta[0]*sign*dw_diag,psi,sqrt_cc[0],psi,s,s+1); // Multiplies Dw by G5 so Hw
+	} else if ( s==(Ls-1) ){
+	  // Drop the CC here.
+	  double R=(1+mass)/(1-mass);
+	  ag5xpby_ssp(chi,Beta[s]*dw_diag,psi,sqrt_cc[s-1],psi,s,s-1);
+	  ag5xpby_ssp(chi,R,psi,1.0,chi,s,s);
+	} else {
+	  ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*dw_diag,psi,sqrt_cc[s],psi,s,s+1);
+	  axpby_ssp(chi,1.0,chi,sqrt_cc[s-1],psi,s,s-1);
+	}
+	sign=-sign; 
+      }
+    }
+
+    void   ContinuedFractionFermion5D::MooeeDag    (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      Mooee(psi,chi);
+    }
+    void   ContinuedFractionFermion5D::MooeeInv    (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      // Apply Linv
+      axpby_ssp(chi,1.0/cc_d[0],psi,0.0,psi,0,0); 
+      for(int s=1;s<Ls;s++){
+	axpbg5y_ssp(chi,1.0/cc_d[s],psi,-1.0/See[s-1],chi,s,s-1);
+      }
+      // Apply Dinv
+      for(int s=0;s<Ls;s++){
+	ag5xpby_ssp(chi,1.0/See[s],chi,0.0,chi,s,s); //only appearance of See[0]
+      }
+      // Apply Uinv = (Linv)^T
+      axpby_ssp(chi,1.0/cc_d[Ls-1],chi,0.0,chi,Ls-1,Ls-1);
+      for(int s=Ls-2;s>=0;s--){
+	axpbg5y_ssp(chi,1.0/cc_d[s],chi,-1.0*cc_d[s+1]/See[s]/cc_d[s],chi,s,s+1);
+      }
+    }
+    void   ContinuedFractionFermion5D::MooeeInvDag (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      MooeeInv(psi,chi);
+    }
+    
+    // Constructors
+    ContinuedFractionFermion5D::ContinuedFractionFermion5D(
+							   LatticeGaugeField &_Umu,
+							   GridCartesian         &FiveDimGrid,
+							   GridRedBlackCartesian &FiveDimRedBlackGrid,
+							   GridCartesian         &FourDimGrid,
+							   GridRedBlackCartesian &FourDimRedBlackGrid,
+							   RealD _mass,RealD M5) :
+      WilsonFermion5D(_Umu,
+		      FiveDimGrid, FiveDimRedBlackGrid,
+		      FourDimGrid, FourDimRedBlackGrid,M5),
+      mass(_mass)
+    {
+      assert((Ls&0x1)==1); // Odd Ls required
+    }
+
+  }
+}
+

lib/qcd/action/fermion/ContinuedFractionFermion5D.h

@@ -0,0 +1,59 @@
+#ifndef  GRID_QCD_CONTINUED_FRACTION_H
+#define  GRID_QCD_CONTINUED_FRACTION_H
+
+namespace Grid {
+
+  namespace QCD {
+
+    class ContinuedFractionFermion5D : public WilsonFermion5D
+    {
+    public:
+
+      // override multiply
+      virtual RealD  M    (const LatticeFermion &in, LatticeFermion &out);
+      virtual RealD  Mdag (const LatticeFermion &in, LatticeFermion &out);
+
+      // half checkerboard operaions
+      virtual void   Meooe       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
+
+      //      virtual void   Instantiatable(void)=0;
+      virtual void   Instantiatable(void) =0;
+
+      // Constructors
+      ContinuedFractionFermion5D(LatticeGaugeField &_Umu,
+				 GridCartesian         &FiveDimGrid,
+				 GridRedBlackCartesian &FiveDimRedBlackGrid,
+				 GridCartesian         &FourDimGrid,
+				 GridRedBlackCartesian &FourDimRedBlackGrid,
+				 RealD _mass,RealD M5);
+
+    protected:
+
+      void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale);
+      void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata);;
+
+      Approx::zolotarev_data *zdata;
+
+      // Cont frac
+      RealD mass;
+      RealD R;
+      RealD ZoloHiInv;
+      std::vector<double> Beta;
+      std::vector<double> cc;;
+      std::vector<double> cc_d;;
+      std::vector<double> sqrt_cc;
+      std::vector<double> See;
+      std::vector<double> Aee;
+
+    };
+
+
+  }
+}
+
+#endif

lib/qcd/action/fermion/DomainWallFermion.h

@@ -0,0 +1,46 @@
+#ifndef  GRID_QCD_DOMAIN_WALL_FERMION_H
+#define  GRID_QCD_DOMAIN_WALL_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class DomainWallFermion : public CayleyFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void) {};
+      // Constructors
+      DomainWallFermion(LatticeGaugeField &_Umu,
+			GridCartesian         &FiveDimGrid,
+			GridRedBlackCartesian &FiveDimRedBlackGrid,
+			GridCartesian         &FourDimGrid,
+			GridRedBlackCartesian &FourDimRedBlackGrid,
+			RealD _mass,RealD _M5) : 
+
+      CayleyFermion5D(_Umu,
+		      FiveDimGrid,
+		      FiveDimRedBlackGrid,
+		      FourDimGrid,
+		      FourDimRedBlackGrid,_mass,_M5)
+
+      {
+	RealD eps = 1.0;
+
+	Approx::zolotarev_data *zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham
+	assert(zdata->n==this->Ls);
+	
+	std::cout << "DomainWallFermion with Ls="<<Ls<<std::endl;
+	// Call base setter
+	this->CayleyFermion5D::SetCoefficientsTanh(zdata,1.0,0.0);
+ 
+      }
+
+    };
+
+  }
+}
+
+#endif

lib/qcd/action/fermion/FermionOperator.h

@@ -0,0 +1,48 @@
+#ifndef  GRID_QCD_FERMION_OPERATOR_H
+#define  GRID_QCD_FERMION_OPERATOR_H
+
+namespace Grid {
+
+  namespace QCD {
+
+    //////////////////////////////////////////////////////////////////////////////
+    // Four component fermions
+    // Should type template the vector and gauge types
+    // Think about multiple representations
+    //////////////////////////////////////////////////////////////////////////////
+    template<class FermionField,class GaugeField>
+    class FermionOperator : public CheckerBoardedSparseMatrixBase<FermionField>
+    {
+    public:
+
+      GridBase * Grid(void)   { return FermionGrid(); };   // this is all the linalg routines need to know
+      GridBase * RedBlackGrid(void) { return FermionRedBlackGrid(); };
+
+      virtual GridBase *FermionGrid(void)         =0;
+      virtual GridBase *FermionRedBlackGrid(void) =0;
+      virtual GridBase *GaugeGrid(void)           =0;
+      virtual GridBase *GaugeRedBlackGrid(void)   =0;
+
+      // override multiply
+      virtual RealD  M    (const FermionField &in, FermionField &out)=0;
+      virtual RealD  Mdag (const FermionField &in, FermionField &out)=0;
+
+      // half checkerboard operaions
+      virtual void   Meooe       (const FermionField &in, FermionField &out)=0;
+      virtual void   MeooeDag    (const FermionField &in, FermionField &out)=0;
+      virtual void   Mooee       (const FermionField &in, FermionField &out)=0;
+      virtual void   MooeeDag    (const FermionField &in, FermionField &out)=0;
+      virtual void   MooeeInv    (const FermionField &in, FermionField &out)=0;
+      virtual void   MooeeInvDag (const FermionField &in, FermionField &out)=0;
+
+      // non-hermitian hopping term; half cb or both
+      virtual void Dhop  (const FermionField &in, FermionField &out,int dag)=0;
+      virtual void DhopOE(const FermionField &in, FermionField &out,int dag)=0;
+      virtual void DhopEO(const FermionField &in, FermionField &out,int dag)=0;
+
+
+    };
+
+  }
+}
+#endif

lib/qcd/action/fermion/MobiusFermion.h

@@ -0,0 +1,47 @@
+#ifndef  GRID_QCD_MOBIUS_FERMION_H
+#define  GRID_QCD_MOBIUS_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class MobiusFermion : public CayleyFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void) {};
+      // Constructors
+      MobiusFermion(LatticeGaugeField &_Umu,
+		    GridCartesian         &FiveDimGrid,
+		    GridRedBlackCartesian &FiveDimRedBlackGrid,
+		    GridCartesian         &FourDimGrid,
+		    GridRedBlackCartesian &FourDimRedBlackGrid,
+		    RealD _mass,RealD _M5,
+		    RealD b, RealD c) : 
+      
+      CayleyFermion5D(_Umu,
+		      FiveDimGrid,
+		      FiveDimRedBlackGrid,
+		      FourDimGrid,
+		      FourDimRedBlackGrid,_mass,_M5)
+
+      {
+	RealD eps = 1.0;
+
+	std::cout << "MobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Tanh approx"<<std::endl;
+	Approx::zolotarev_data *zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham
+	assert(zdata->n==this->Ls);
+	
+	// Call base setter
+	this->CayleyFermion5D::SetCoefficientsTanh(zdata,b,c);
+ 
+      }
+
+    };
+
+  }
+}
+
+#endif

lib/qcd/action/fermion/MobiusZolotarevFermion.h

@@ -0,0 +1,49 @@
+#ifndef  GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H
+#define  GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class MobiusZolotarevFermion : public CayleyFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void) {};
+      // Constructors
+       MobiusZolotarevFermion(LatticeGaugeField &_Umu,
+			      GridCartesian         &FiveDimGrid,
+			      GridRedBlackCartesian &FiveDimRedBlackGrid,
+			      GridCartesian         &FourDimGrid,
+			      GridRedBlackCartesian &FourDimRedBlackGrid,
+			      RealD _mass,RealD _M5,
+			      RealD b, RealD c,
+			      RealD lo, RealD hi) : 
+      
+      CayleyFermion5D(_Umu,
+		      FiveDimGrid,
+		      FiveDimRedBlackGrid,
+		      FourDimGrid,
+		      FourDimRedBlackGrid,_mass,_M5)
+
+      {
+	RealD eps = lo/hi;
+
+	Approx::zolotarev_data *zdata = Approx::grid_zolotarev(eps,this->Ls,0);// eps is ignored for higham
+	assert(zdata->n==this->Ls);
+
+	std::cout << "MobiusZolotarevFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Zolotarev range ["<<lo<<","<<hi<<"]"<<std::endl;
+	
+	// Call base setter
+	this->CayleyFermion5D::SetCoefficientsZolotarev(hi,zdata,b,c);
+ 
+      }
+
+    };
+
+  }
+}
+
+#endif

lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h

@@ -0,0 +1,34 @@
+#ifndef OVERLAP_WILSON_CAYLEY_TANH_FERMION_H
+#define OVERLAP_WILSON_CAYLEY_TANH_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class OverlapWilsonCayleyTanhFermion : public MobiusFermion
+    {
+    public:
+
+      // Constructors
+    OverlapWilsonCayleyTanhFermion(LatticeGaugeField &_Umu,
+				   GridCartesian         &FiveDimGrid,
+				   GridRedBlackCartesian &FiveDimRedBlackGrid,
+				   GridCartesian         &FourDimGrid,
+				   GridRedBlackCartesian &FourDimRedBlackGrid,
+				   RealD _mass,RealD _M5,
+				   RealD scale) :
+      
+      // b+c=scale, b-c = 0 <=> b =c = scale/2
+      MobiusFermion(_Umu,
+		    FiveDimGrid,
+		    FiveDimRedBlackGrid,
+		    FourDimGrid,
+		    FourDimRedBlackGrid,_mass,_M5,0.5*scale,0.5*scale)
+	{
+	}
+    };
+  }
+}
+#endif

lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h

@@ -0,0 +1,37 @@
+#ifndef  OVERLAP_WILSON_CAYLEY_ZOLOTAREV_FERMION_H
+#define  OVERLAP_WILSON_CAYLEY_ZOLOTAREV_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class OverlapWilsonCayleyZolotarevFermion : public MobiusZolotarevFermion
+    {
+    public:
+
+      // Constructors
+
+    OverlapWilsonCayleyZolotarevFermion(LatticeGaugeField &_Umu,
+					GridCartesian         &FiveDimGrid,
+					GridRedBlackCartesian &FiveDimRedBlackGrid,
+					GridCartesian         &FourDimGrid,
+					GridRedBlackCartesian &FourDimRedBlackGrid,
+					RealD _mass,RealD _M5,
+					RealD lo, RealD hi) : 
+      // b+c=1.0, b-c = 0 <=> b =c = 1/2
+      MobiusZolotarevFermion(_Umu,
+			     FiveDimGrid,
+			     FiveDimRedBlackGrid,
+			     FourDimGrid,
+			     FourDimRedBlackGrid,_mass,_M5,0.5,0.5,lo,hi)
+
+      {}
+
+    };
+
+  }
+}
+
+#endif

lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h

@@ -0,0 +1,39 @@
+#ifndef OVERLAP_WILSON_CONTFRAC_TANH_FERMION_H
+#define OVERLAP_WILSON_CONTFRAC_TANH_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class OverlapWilsonContFracTanhFermion : public ContinuedFractionFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void){};
+      // Constructors
+    OverlapWilsonContFracTanhFermion(LatticeGaugeField &_Umu,
+				     GridCartesian         &FiveDimGrid,
+				     GridRedBlackCartesian &FiveDimRedBlackGrid,
+				     GridCartesian         &FourDimGrid,
+				     GridRedBlackCartesian &FourDimRedBlackGrid,
+				     RealD _mass,RealD _M5,
+				     RealD scale) :
+      
+      // b+c=scale, b-c = 0 <=> b =c = scale/2
+      ContinuedFractionFermion5D(_Umu,
+				 FiveDimGrid,
+				 FiveDimRedBlackGrid,
+				 FourDimGrid,
+				 FourDimRedBlackGrid,_mass,_M5)
+	{
+	  assert((Ls&0x1)==1); // Odd Ls required
+	  int nrational=Ls-1;// Even rational order
+	  zdata = Approx::grid_higham(1.0,nrational);// eps is ignored for higham
+	  SetCoefficientsTanh(zdata,scale);
+	}
+    };
+  }
+}
+#endif

lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h

@@ -0,0 +1,44 @@
+#ifndef OVERLAP_WILSON_CONTFRAC_ZOLOTAREV_FERMION_H
+#define OVERLAP_WILSON_CONTFRAC_ZOLOTAREV_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class OverlapWilsonContFracZolotarevFermion : public ContinuedFractionFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void){};
+      // Constructors
+    OverlapWilsonContFracZolotarevFermion(LatticeGaugeField &_Umu,
+					  GridCartesian         &FiveDimGrid,
+					  GridRedBlackCartesian &FiveDimRedBlackGrid,
+					  GridCartesian         &FourDimGrid,
+					  GridRedBlackCartesian &FourDimRedBlackGrid,
+					  RealD _mass,RealD _M5,
+					  RealD lo,RealD hi):
+      
+      // b+c=scale, b-c = 0 <=> b =c = scale/2
+      ContinuedFractionFermion5D(_Umu,
+				 FiveDimGrid,
+				 FiveDimRedBlackGrid,
+				 FourDimGrid,
+				 FourDimRedBlackGrid,_mass,_M5)
+	{
+	  assert((Ls&0x1)==1); // Odd Ls required
+
+	  int nrational=Ls-1;// Even rational order
+	  RealD eps = lo/hi;
+
+	  Approx::zolotarev_data *zdata = Approx::grid_zolotarev(eps,nrational,0);
+
+	  SetCoefficientsZolotarev(hi,zdata);
+
+	}
+    };
+  }
+}
+#endif

lib/qcd/action/fermion/PartialFractionFermion5D.cc

@@ -0,0 +1 @@
+

lib/qcd/action/fermion/PartialFractionFermion5D.h

@@ -0,0 +1,49 @@
+#ifndef  GRID_QCD_PARTIAL_FRACTION_H
+#define  GRID_QCD_PARTIAL_FRACTION_H
+
+namespace Grid {
+
+  namespace QCD {
+
+    class PartialFractionFermion5D : public WilsonFermion5D
+    {
+    public:
+
+      // override multiply
+      virtual RealD  M    (const LatticeFermion &in, LatticeFermion &out);
+      virtual RealD  Mdag (const LatticeFermion &in, LatticeFermion &out);
+
+      // half checkerboard operaions
+      virtual void   Meooe       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
+
+    private:
+
+      virtual void PartialFractionCoefficients(void);
+
+      Approx::zolotarev_data *zdata;
+
+      // Part frac
+      double R;
+      std::vector<double> p; 
+      std::vector<double> q;
+
+      // Constructors
+      PartialFractionFermion5D(LatticeGaugeField &_Umu,
+				    GridCartesian         &FiveDimGrid,
+				    GridRedBlackCartesian &FiveDimRedBlackGrid,
+				    GridCartesian         &FourDimGrid,
+				    GridRedBlackCartesian &FourDimRedBlackGrid,
+				    RealD _mass,RealD M5);
+
+    };
+
+
+  }
+}
+
+#endif

lib/qcd/action/fermion/ScaledShamirFermion.h

@@ -0,0 +1,37 @@
+#ifndef  GRID_QCD_SCALED_SHAMIR_FERMION_H
+#define  GRID_QCD_SCALED_SHAMIR_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class ScaledShamirFermion : public MobiusFermion
+    {
+    public:
+
+      // Constructors
+    ScaledShamirFermion(LatticeGaugeField &_Umu,
+			GridCartesian         &FiveDimGrid,
+			GridRedBlackCartesian &FiveDimRedBlackGrid,
+			GridCartesian         &FourDimGrid,
+			GridRedBlackCartesian &FourDimRedBlackGrid,
+			RealD _mass,RealD _M5,
+			RealD scale) :
+      
+      // b+c=scale, b-c = 1 <=> 2b = scale+1; 2c = scale-1
+      MobiusFermion(_Umu,
+		    FiveDimGrid,
+		    FiveDimRedBlackGrid,
+		    FourDimGrid,
+		    FourDimRedBlackGrid,_mass,_M5,0.5*(scale+1.0),0.5*(scale-1.0))
+      {
+      }
+
+    };
+
+  }
+}
+
+#endif

lib/qcd/action/fermion/ShamirZolotarevFermion.h

@@ -0,0 +1,39 @@
+#ifndef  GRID_QCD_SHAMIR_ZOLOTAREV_FERMION_H
+#define  GRID_QCD_SHAMIR_ZOLOTAREV_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class ShamirZolotarevFermion : public MobiusZolotarevFermion
+    {
+    public:
+
+      // Constructors
+
+
+    ShamirZolotarevFermion(LatticeGaugeField &_Umu,
+			   GridCartesian         &FiveDimGrid,
+			   GridRedBlackCartesian &FiveDimRedBlackGrid,
+			   GridCartesian         &FourDimGrid,
+			   GridRedBlackCartesian &FourDimRedBlackGrid,
+			   RealD _mass,RealD _M5,
+			   RealD lo, RealD hi) : 
+      
+      // b+c = 1; b-c = 1 => b=1, c=0
+      MobiusZolotarevFermion(_Umu,
+			     FiveDimGrid,
+			     FiveDimRedBlackGrid,
+			     FourDimGrid,
+			     FourDimRedBlackGrid,_mass,_M5,1.0,0.0,lo,hi)
+      
+      {}
+
+    };
+
+  }
+}
+
+#endif

lib/qcd/action/fermion/WilsonCompressor.h

@@ -0,0 +1,61 @@
+#ifndef  GRID_QCD_WILSON_COMPRESSOR_H
+#define  GRID_QCD_WILSON_COMPRESSOR_H
+
+namespace Grid {
+namespace QCD {
+
+  class WilsonCompressor {
+  public:
+    int mu;
+    int dag;
+
+    WilsonCompressor(int _dag){
+      mu=0;
+      dag=_dag;
+      assert((dag==0)||(dag==1));
+    }
+    void Point(int p) { 
+      mu=p;
+    };
+
+    vHalfSpinColourVector operator () (const vSpinColourVector &in)
+    {
+      vHalfSpinColourVector ret;
+      int mudag=mu;
+      if (dag) {
+	mudag=(mu+Nd)%(2*Nd);
+      }
+      switch(mudag) {
+      case Xp:
+	spProjXp(ret,in);
+	break;
+      case Yp:
+	spProjYp(ret,in);
+	break;
+      case Zp:
+	spProjZp(ret,in);
+	break;
+      case Tp:
+	spProjTp(ret,in);
+	break;
+      case Xm:
+	spProjXm(ret,in);
+	break;
+      case Ym:
+	spProjYm(ret,in);
+	break;
+      case Zm:
+	spProjZm(ret,in);
+	break;
+      case Tm:
+	spProjTm(ret,in);
+	break;
+      default: 
+	assert(0);
+	break;
+      }
+      return ret;
+    }
+  };
+}} // namespace close
+#endif

lib/qcd/action/fermion/WilsonFermion.cc

@@ -0,0 +1,163 @@
+#include <Grid.h>
+
+namespace Grid {
+namespace QCD {
+
+const std::vector<int> WilsonFermion::directions   ({0,1,2,3, 0, 1, 2, 3});
+const std::vector<int> WilsonFermion::displacements({1,1,1,1,-1,-1,-1,-1});
+
+int WilsonFermion::HandOptDslash;
+
+WilsonFermion::WilsonFermion(LatticeGaugeField &_Umu,
+			     GridCartesian         &Fgrid,
+			     GridRedBlackCartesian &Hgrid, 
+			     RealD _mass) :
+  _grid(&Fgrid),
+  _cbgrid(&Hgrid),
+  Stencil    (&Fgrid,npoint,Even,directions,displacements),
+  StencilEven(&Hgrid,npoint,Even,directions,displacements), // source is Even
+  StencilOdd (&Hgrid,npoint,Odd ,directions,displacements), // source is Odd
+  mass(_mass),
+  Umu(&Fgrid),
+  UmuEven(&Hgrid),
+  UmuOdd (&Hgrid)
+{
+  // Allocate the required comms buffer
+  comm_buf.resize(Stencil._unified_buffer_size); // this is always big enough to contain EO
+  DoubleStore(Umu,_Umu);
+  pickCheckerboard(Even,UmuEven,Umu);
+  pickCheckerboard(Odd ,UmuOdd,Umu);
+}
+      
+void WilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu)
+{
+  conformable(Uds._grid,GaugeGrid());
+  conformable(Umu._grid,GaugeGrid());
+  LatticeColourMatrix U(GaugeGrid());
+  for(int mu=0;mu<Nd;mu++){
+    U = peekIndex<LorentzIndex>(Umu,mu);
+    pokeIndex<LorentzIndex>(Uds,U,mu);
+    U = adj(Cshift(U,mu,-1));
+    pokeIndex<LorentzIndex>(Uds,U,mu+4);
+  }
+}
+
+RealD WilsonFermion::M(const LatticeFermion &in, LatticeFermion &out)
+{
+  out.checkerboard=in.checkerboard;
+  Dhop(in,out,DaggerNo);
+  return axpy_norm(out,4+mass,in,out);
+}
+RealD WilsonFermion::Mdag(const LatticeFermion &in, LatticeFermion &out)
+{
+  out.checkerboard=in.checkerboard;
+  Dhop(in,out,DaggerYes);
+  return axpy_norm(out,4+mass,in,out);
+}
+
+void WilsonFermion::Meooe(const LatticeFermion &in, LatticeFermion &out)
+{
+  if ( in.checkerboard == Odd ) {
+    DhopEO(in,out,DaggerNo);
+  } else {
+    DhopOE(in,out,DaggerNo);
+  }
+}
+void WilsonFermion::MeooeDag(const LatticeFermion &in, LatticeFermion &out)
+{
+  if ( in.checkerboard == Odd ) {
+    DhopEO(in,out,DaggerYes);
+  } else {
+    DhopOE(in,out,DaggerYes);
+  }
+}
+void WilsonFermion::Mooee(const LatticeFermion &in, LatticeFermion &out)
+{
+  out.checkerboard = in.checkerboard;
+  out = (4.0+mass)*in;
+  return ;
+}
+void WilsonFermion::MooeeDag(const LatticeFermion &in, LatticeFermion &out)
+{
+  out.checkerboard = in.checkerboard;
+  Mooee(in,out);
+}
+void WilsonFermion::MooeeInv(const LatticeFermion &in, LatticeFermion &out)
+{
+  out.checkerboard = in.checkerboard;
+  out = (1.0/(4.0+mass))*in;
+  return ;
+}
+void WilsonFermion::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out)
+{
+  out.checkerboard = in.checkerboard;
+  MooeeInv(in,out);
+}
+
+void WilsonFermion::DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField & U,
+				const LatticeFermion &in, LatticeFermion &out,int dag)
+{
+  assert((dag==DaggerNo) ||(dag==DaggerYes));
+  WilsonCompressor compressor(dag);
+  st.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
+
+  if ( dag == DaggerYes ) {
+    if( HandOptDslash ) {
+PARALLEL_FOR_LOOP
+      for(int sss=0;sss<in._grid->oSites();sss++){
+        DiracOptHand::DhopSiteDag(st,U,comm_buf,sss,sss,in,out);
+      }
+    } else { 
+PARALLEL_FOR_LOOP
+      for(int sss=0;sss<in._grid->oSites();sss++){
+        DiracOpt::DhopSiteDag(st,U,comm_buf,sss,sss,in,out);
+      }
+    }
+  } else {
+    if( HandOptDslash ) {
+PARALLEL_FOR_LOOP
+      for(int sss=0;sss<in._grid->oSites();sss++){
+        DiracOptHand::DhopSite(st,U,comm_buf,sss,sss,in,out);
+      }
+    } else { 
+PARALLEL_FOR_LOOP
+      for(int sss=0;sss<in._grid->oSites();sss++){
+        DiracOpt::DhopSite(st,U,comm_buf,sss,sss,in,out);
+      }
+    }
+  }
+}
+void WilsonFermion::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag)
+{
+  conformable(in._grid,_cbgrid);    // verifies half grid
+  conformable(in._grid,out._grid); // drops the cb check
+
+  assert(in.checkerboard==Even);
+  out.checkerboard = Odd;
+
+  DhopInternal(StencilEven,UmuOdd,in,out,dag);
+}
+void WilsonFermion::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag)
+{
+  conformable(in._grid,_cbgrid);    // verifies half grid
+  conformable(in._grid,out._grid); // drops the cb check
+
+  assert(in.checkerboard==Odd);
+  out.checkerboard = Even;
+
+  DhopInternal(StencilOdd,UmuEven,in,out,dag);
+}
+void WilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag)
+{
+  conformable(in._grid,_grid); // verifies full grid
+  conformable(in._grid,out._grid);
+
+  out.checkerboard = in.checkerboard;
+
+  DhopInternal(Stencil,Umu,in,out,dag);
+}
+
+}}
+
+
+

lib/qcd/action/fermion/WilsonFermion.h

@@ -0,0 +1,87 @@
+#ifndef  GRID_QCD_WILSON_FERMION_H
+#define  GRID_QCD_WILSON_FERMION_H
+
+namespace Grid {
+
+  namespace QCD {
+
+    class WilsonFermion : public FermionOperator<LatticeFermion,LatticeGaugeField>
+    {
+    public:
+
+      ///////////////////////////////////////////////////////////////
+      // Implement the abstract base
+      ///////////////////////////////////////////////////////////////
+      GridBase *GaugeGrid(void)              { return _grid ;}
+      GridBase *GaugeRedBlackGrid(void)      { return _cbgrid ;}
+      GridBase *FermionGrid(void)            { return _grid;}
+      GridBase *FermionRedBlackGrid(void)    { return _cbgrid;}
+
+      // override multiply
+      virtual RealD  M    (const LatticeFermion &in, LatticeFermion &out);
+      virtual RealD  Mdag (const LatticeFermion &in, LatticeFermion &out);
+
+      // half checkerboard operaions
+      void   Meooe       (const LatticeFermion &in, LatticeFermion &out);
+      void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out); // remain virtual so we 
+      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out); // can derive Clover
+      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out); // from Wilson bas
+      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
+
+      // non-hermitian hopping term; half cb or both
+      void Dhop  (const LatticeFermion &in, LatticeFermion &out,int dag);
+      void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag);
+      void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag);
+
+      ///////////////////////////////////////////////////////////////
+      // Extra methods added by derived
+      ///////////////////////////////////////////////////////////////
+      void DhopInternal(CartesianStencil & st,
+			LatticeDoubledGaugeField &U,
+			const LatticeFermion &in, 
+			LatticeFermion &out,
+			int dag);
+
+      // Constructor
+      WilsonFermion(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,RealD _mass);
+
+      // DoubleStore
+      void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
+
+      ///////////////////////////////////////////////////////////////
+      // Data members require to support the functionality
+      ///////////////////////////////////////////////////////////////
+      static int HandOptDslash; // these are a temporary hack
+      static int MortonOrder;
+
+    protected:
+
+      RealD                        mass;
+
+      GridBase                     *    _grid; 
+      GridBase                     *  _cbgrid;
+
+      static const int npoint=8;
+      static const std::vector<int> directions   ;
+      static const std::vector<int> displacements;
+
+      //Defines the stencils for even and odd
+      CartesianStencil Stencil; 
+      CartesianStencil StencilEven; 
+      CartesianStencil StencilOdd; 
+
+      // Copy of the gauge field , with even and odd subsets
+      LatticeDoubledGaugeField Umu;
+      LatticeDoubledGaugeField UmuEven;
+      LatticeDoubledGaugeField UmuOdd;
+
+      // Comms buffer
+      std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  comm_buf;
+
+      
+    };
+
+  }
+}
+#endif

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -0,0 +1,186 @@
+#include <Grid.h>
+
+namespace Grid {
+namespace QCD {
+  
+  // S-direction is INNERMOST and takes no part in the parity.
+  const std::vector<int> WilsonFermion5D::directions   ({1,2,3,4, 1, 2, 3, 4});
+  const std::vector<int> WilsonFermion5D::displacements({1,1,1,1,-1,-1,-1,-1});
+
+  int WilsonFermion5D::HandOptDslash;
+
+  // 5d lattice for DWF.
+  WilsonFermion5D::WilsonFermion5D(LatticeGaugeField &_Umu,
+					   GridCartesian         &FiveDimGrid,
+					   GridRedBlackCartesian &FiveDimRedBlackGrid,
+					   GridCartesian         &FourDimGrid,
+					   GridRedBlackCartesian &FourDimRedBlackGrid,
+					   RealD _M5) :
+  _FiveDimGrid(&FiveDimGrid),
+  _FiveDimRedBlackGrid(&FiveDimRedBlackGrid),
+  _FourDimGrid(&FourDimGrid),
+  _FourDimRedBlackGrid(&FourDimRedBlackGrid),
+  Stencil    (_FiveDimGrid,npoint,Even,directions,displacements),
+  StencilEven(_FiveDimRedBlackGrid,npoint,Even,directions,displacements), // source is Even
+  StencilOdd (_FiveDimRedBlackGrid,npoint,Odd ,directions,displacements), // source is Odd
+  M5(_M5),
+  Umu(_FourDimGrid),
+  UmuEven(_FourDimRedBlackGrid),
+  UmuOdd (_FourDimRedBlackGrid),
+  Lebesgue(_FourDimGrid),
+  LebesgueEvenOdd(_FourDimRedBlackGrid)
+{
+  // some assertions
+  assert(FiveDimGrid._ndimension==5);
+  assert(FourDimGrid._ndimension==4);
+  
+  assert(FiveDimRedBlackGrid._ndimension==5);
+  assert(FourDimRedBlackGrid._ndimension==4);
+
+  assert(FiveDimRedBlackGrid._checker_dim==1);
+
+  // Dimension zero of the five-d is the Ls direction
+  Ls=FiveDimGrid._fdimensions[0];
+  assert(FiveDimRedBlackGrid._fdimensions[0]==Ls);
+  assert(FiveDimRedBlackGrid._processors[0] ==1);
+  assert(FiveDimRedBlackGrid._simd_layout[0]==1);
+  assert(FiveDimGrid._processors[0]         ==1);
+  assert(FiveDimGrid._simd_layout[0]        ==1);
+
+  // Other dimensions must match the decomposition of the four-D fields 
+  for(int d=0;d<4;d++){
+    assert(FourDimRedBlackGrid._fdimensions[d]  ==FourDimGrid._fdimensions[d]);
+    assert(FiveDimRedBlackGrid._fdimensions[d+1]==FourDimGrid._fdimensions[d]);
+
+    assert(FourDimRedBlackGrid._processors[d]   ==FourDimGrid._processors[d]);
+    assert(FiveDimRedBlackGrid._processors[d+1] ==FourDimGrid._processors[d]);
+
+    assert(FourDimRedBlackGrid._simd_layout[d]  ==FourDimGrid._simd_layout[d]);
+    assert(FiveDimRedBlackGrid._simd_layout[d+1]==FourDimGrid._simd_layout[d]);
+
+    assert(FiveDimGrid._fdimensions[d+1]        ==FourDimGrid._fdimensions[d]);
+    assert(FiveDimGrid._processors[d+1]         ==FourDimGrid._processors[d]);
+    assert(FiveDimGrid._simd_layout[d+1]        ==FourDimGrid._simd_layout[d]);
+  }
+
+  // Allocate the required comms buffer
+  comm_buf.resize(Stencil._unified_buffer_size); // this is always big enough to contain EO
+  
+  DoubleStore(Umu,_Umu);
+  pickCheckerboard(Even,UmuEven,Umu);
+  pickCheckerboard(Odd ,UmuOdd,Umu);
+}
+void WilsonFermion5D::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu)
+{
+  conformable(Uds._grid,GaugeGrid());
+  conformable(Umu._grid,GaugeGrid());
+  LatticeColourMatrix U(GaugeGrid());
+  for(int mu=0;mu<Nd;mu++){
+    U = peekIndex<LorentzIndex>(Umu,mu);
+    pokeIndex<LorentzIndex>(Uds,U,mu);
+    U = adj(Cshift(U,mu,-1));
+    pokeIndex<LorentzIndex>(Uds,U,mu+4);
+  }
+}
+void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
+				   LatticeDoubledGaugeField & U,
+				   const LatticeFermion &in, LatticeFermion &out,int dag)
+{
+  assert((dag==DaggerNo) ||(dag==DaggerYes));
+
+  WilsonCompressor compressor(dag);
+
+  st.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
+  
+  // Dhop takes the 4d grid from U, and makes a 5d index for fermion
+  // Not loop ordering and data layout.
+  // Designed to create 
+  // - per thread reuse in L1 cache for U
+  // - 8 linear access unit stride streams per thread for Fermion for hw prefetchable.
+  if ( dag == DaggerYes ) {
+    if( HandOptDslash ) {
+PARALLEL_FOR_LOOP
+      for(int ss=0;ss<U._grid->oSites();ss++){
+	for(int s=0;s<Ls;s++){
+	  //int sU=lo.Reorder(ss);
+	  int sU=ss;
+	  int sF = s+Ls*sU;
+	  DiracOptHand::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
+	}
+      }
+    } else { 
+PARALLEL_FOR_LOOP
+      for(int ss=0;ss<U._grid->oSites();ss++){
+	for(int s=0;s<Ls;s++){
+	  //	  int sU=lo.Reorder(ss);
+	  int sU=ss;
+	  int sF = s+Ls*sU;
+	  DiracOpt::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
+	}
+      }
+    }
+  } else {
+    if( HandOptDslash ) {
+PARALLEL_FOR_LOOP
+      for(int ss=0;ss<U._grid->oSites();ss++){
+	for(int s=0;s<Ls;s++){
+	  //	  int sU=lo.Reorder(ss);
+	  int sU=ss;
+	  int sF = s+Ls*sU;
+	  DiracOptHand::DhopSite(st,U,comm_buf,sF,sU,in,out);
+	}
+      }
+
+    } else { 
+PARALLEL_FOR_LOOP
+      for(int ss=0;ss<U._grid->oSites();ss++){
+	for(int s=0;s<Ls;s++){
+	  //	  int sU=lo.Reorder(ss);
+	  int sU=ss;
+	  int sF = s+Ls*sU; 
+	  DiracOpt::DhopSite(st,U,comm_buf,sF,sU,in,out);
+	}
+      }
+    }
+  }
+}
+void WilsonFermion5D::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag)
+{
+  conformable(in._grid,FermionRedBlackGrid());    // verifies half grid
+  conformable(in._grid,out._grid); // drops the cb check
+
+  assert(in.checkerboard==Even);
+  out.checkerboard = Odd;
+
+  DhopInternal(StencilEven,LebesgueEvenOdd,UmuOdd,in,out,dag);
+}
+void WilsonFermion5D::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag)
+{
+  conformable(in._grid,FermionRedBlackGrid());    // verifies half grid
+  conformable(in._grid,out._grid); // drops the cb check
+
+  assert(in.checkerboard==Odd);
+  out.checkerboard = Even;
+
+  DhopInternal(StencilOdd,LebesgueEvenOdd,UmuEven,in,out,dag);
+}
+void WilsonFermion5D::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag)
+{
+  conformable(in._grid,FermionGrid()); // verifies full grid
+  conformable(in._grid,out._grid);
+
+  out.checkerboard = in.checkerboard;
+
+  DhopInternal(Stencil,Lebesgue,Umu,in,out,dag);
+}
+void WilsonFermion5D::DW(const LatticeFermion &in, LatticeFermion &out,int dag)
+{
+  out.checkerboard=in.checkerboard;
+  Dhop(in,out,dag); // -0.5 is included
+  axpy(out,4.0-M5,in,out);
+}
+}
+}
+
+
+

lib/qcd/action/fermion/WilsonFermion5D.h

@@ -0,0 +1,121 @@
+#ifndef  GRID_QCD_DWF_H
+#define  GRID_QCD_DWF_H
+
+namespace Grid {
+
+  namespace QCD {
+
+    ////////////////////////////////////////////////////////////////////////////////
+    // This is the 4d red black case appropriate to support
+    //
+    // parity = (x+y+z+t)|2;
+    // generalised five dim fermions like mobius, zolotarev etc..	
+    //
+    // i.e. even even contains fifth dim hopping term.
+    //
+    // [DIFFERS from original CPS red black implementation parity = (x+y+z+t+s)|2 ]
+    ////////////////////////////
+    //ContFrac:
+    //  Ls always odd. Rational poly deg is either Ls or Ls-1
+    //PartFrac 
+    //  Ls always odd. Rational poly deg is either Ls or Ls-1
+    //
+    //Cayley: Ls always even, Rational poly deg is Ls
+    // 
+    // Just set nrational as Ls. Forget about Ls-1 cases.
+    //
+    // Require odd Ls for cont and part frac
+    ////////////////////////////
+    ////////////////////////////////////////////////////////////////////////////////
+    class WilsonFermion5D : public FermionOperator<LatticeFermion,LatticeGaugeField>
+    {
+    public:
+      ///////////////////////////////////////////////////////////////
+      // Implement the abstract base
+      ///////////////////////////////////////////////////////////////
+      GridBase *GaugeGrid(void)              { return _FourDimGrid ;}
+      GridBase *GaugeRedBlackGrid(void)      { return _FourDimRedBlackGrid ;}
+      GridBase *FermionGrid(void)            { return _FiveDimGrid;}
+      GridBase *FermionRedBlackGrid(void)    { return _FiveDimRedBlackGrid;}
+
+      // full checkerboard operations; leave unimplemented as abstract for now
+      //virtual RealD  M    (const LatticeFermion &in, LatticeFermion &out)=0;
+      //virtual RealD  Mdag (const LatticeFermion &in, LatticeFermion &out)=0;
+
+      // half checkerboard operations; leave unimplemented as abstract for now
+      //      virtual void   Meooe       (const LatticeFermion &in, LatticeFermion &out)=0;
+      //      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out)=0;
+      //      virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out)=0;
+      //      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out)=0;
+      //      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out)=0;
+      //      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out)=0;
+
+      // Implement hopping term non-hermitian hopping term; half cb or both
+      // Implement s-diagonal DW
+      void DW    (const LatticeFermion &in, LatticeFermion &out,int dag);
+      void Dhop  (const LatticeFermion &in, LatticeFermion &out,int dag);
+      void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag);
+      void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag);
+
+      ///////////////////////////////////////////////////////////////
+      // New methods added 
+      ///////////////////////////////////////////////////////////////
+      void DhopInternal(CartesianStencil & st,
+			LebesgueOrder &lo,
+			LatticeDoubledGaugeField &U,
+			const LatticeFermion &in, 
+			LatticeFermion &out,
+			int dag);
+
+      // Constructors
+      WilsonFermion5D(LatticeGaugeField &_Umu,
+			  GridCartesian         &FiveDimGrid,
+			  GridRedBlackCartesian &FiveDimRedBlackGrid,
+			  GridCartesian         &FourDimGrid,
+			  GridRedBlackCartesian &FourDimRedBlackGrid,
+			  double _M5);
+
+      // DoubleStore
+      void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
+
+      ///////////////////////////////////////////////////////////////
+      // Data members require to support the functionality
+      ///////////////////////////////////////////////////////////////
+      static int HandOptDslash; // these are a temporary hack
+
+    protected:
+
+      // Add these to the support from Wilson
+      GridBase *_FourDimGrid;
+      GridBase *_FourDimRedBlackGrid;
+      GridBase *_FiveDimGrid;
+      GridBase *_FiveDimRedBlackGrid;
+
+      static const int npoint=8;
+      static const std::vector<int> directions   ;
+      static const std::vector<int> displacements;
+
+      double                        M5;
+      int Ls;
+
+      //Defines the stencils for even and odd
+      CartesianStencil Stencil; 
+      CartesianStencil StencilEven; 
+      CartesianStencil StencilOdd; 
+
+      // Copy of the gauge field , with even and odd subsets
+      LatticeDoubledGaugeField Umu;
+      LatticeDoubledGaugeField UmuEven;
+      LatticeDoubledGaugeField UmuOdd;
+
+      LebesgueOrder Lebesgue;
+      LebesgueOrder LebesgueEvenOdd;
+
+      // Comms buffer
+      std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  comm_buf;
+      
+    };
+  }
+}
+
+#endif

lib/qcd/action/fermion/WilsonKernels.cc

@@ -4,8 +4,8 @@ namespace Grid {
 namespace QCD {
 
 void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
-			    std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
+			std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
-			    int ss,const LatticeFermion &in, LatticeFermion &out)
+			int sF,int sU,const LatticeFermion &in, LatticeFermion &out)
 {
     vHalfSpinColourVector  tmp;    
     vHalfSpinColourVector  chi;    
@@ -16,6 +16,7 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     //#define VERBOSE( A)  if ( ss<10 ) { std::cout << "site " <<ss << " " #A " neigh " << offset << " perm "<< perm <<std::endl;}    
 
     // Xp
+    int ss = sF;
     offset = st._offsets [Xp][ss];
     local  = st._is_local[Xp][ss];
     perm   = st._permute[Xp][ss];
@@ -29,7 +30,7 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Xp),&chi());
+    mult(&Uchi(),&U._odata[sU](Xp),&chi());
     spReconXp(result,Uchi);
 
     //    std::cout << "XP_RECON"<<std::endl;
@@ -51,7 +52,7 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Yp),&chi());
+    mult(&Uchi(),&U._odata[sU](Yp),&chi());
     accumReconYp(result,Uchi);
 
     // Zp
@@ -67,7 +68,7 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Zp),&chi());
+    mult(&Uchi(),&U._odata[sU](Zp),&chi());
     accumReconZp(result,Uchi);
 
     // Tp
@@ -83,7 +84,7 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Tp),&chi());
+    mult(&Uchi(),&U._odata[sU](Tp),&chi());
     accumReconTp(result,Uchi);
 
     // Xm
@@ -101,7 +102,7 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Xm),&chi());
+    mult(&Uchi(),&U._odata[sU](Xm),&chi());
     accumReconXm(result,Uchi);
     //  std::cout << "XM_RECON_ACCUM"<<std::endl;
     //    std::cout << result()(0)(0) <<" "<<result()(0)(1) <<" "<<result()(0)(2) <<std::endl;
@@ -124,7 +125,7 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Ym),&chi());
+    mult(&Uchi(),&U._odata[sU](Ym),&chi());
     accumReconYm(result,Uchi);
 
     // Zm
@@ -140,7 +141,7 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Zm),&chi());
+    mult(&Uchi(),&U._odata[sU](Zm),&chi());
     accumReconZm(result,Uchi);
 
     // Tm
@@ -156,15 +157,15 @@ void DiracOpt::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Tm),&chi());
+    mult(&Uchi(),&U._odata[sU](Tm),&chi());
     accumReconTm(result,Uchi);
 
-    vstream(out._odata[ss],result);
+    vstream(out._odata[ss],result*(-0.5));
 }
 
 void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
-			       std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
+			   std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
-			       int ss,const LatticeFermion &in, LatticeFermion &out)
+			   int sF,int sU,const LatticeFermion &in, LatticeFermion &out)
 {
     vHalfSpinColourVector  tmp;    
     vHalfSpinColourVector  chi;    
@@ -173,6 +174,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
     int offset,local,perm, ptype;
 
     // Xp
+    int ss=sF;
     offset = st._offsets [Xm][ss];
     local  = st._is_local[Xm][ss];
     perm   = st._permute[Xm][ss];
@@ -186,7 +188,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Xm),&chi());
+    mult(&Uchi(),&U._odata[sU](Xm),&chi());
     spReconXp(result,Uchi);
 
     // Yp
@@ -202,7 +204,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Ym),&chi());
+    mult(&Uchi(),&U._odata[sU](Ym),&chi());
     accumReconYp(result,Uchi);
 
     // Zp
@@ -218,7 +220,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Zm),&chi());
+    mult(&Uchi(),&U._odata[sU](Zm),&chi());
     accumReconZp(result,Uchi);
 
     // Tp
@@ -234,7 +236,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Tm),&chi());
+    mult(&Uchi(),&U._odata[sU](Tm),&chi());
     accumReconTp(result,Uchi);
 
     // Xm
@@ -252,7 +254,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Xp),&chi());
+    mult(&Uchi(),&U._odata[sU](Xp),&chi());
     accumReconXm(result,Uchi);
 
     // Ym
@@ -269,7 +271,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Yp),&chi());
+    mult(&Uchi(),&U._odata[sU](Yp),&chi());
     accumReconYm(result,Uchi);
 
     // Zm
@@ -285,7 +287,7 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Zp),&chi());
+    mult(&Uchi(),&U._odata[sU](Zp),&chi());
     accumReconZm(result,Uchi);
 
     // Tm
@@ -301,9 +303,9 @@ void DiracOpt::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
     } else { 
       chi=buf[offset];
     }
-    mult(&Uchi(),&U._odata[ss](Tp),&chi());
+    mult(&Uchi(),&U._odata[sU](Tp),&chi());
     accumReconTm(result,Uchi);
 
-    vstream(out._odata[ss],result);
+    vstream(out._odata[ss],result*(-0.5));
 }
 }}

lib/qcd/action/fermion/WilsonKernels.h

@@ -0,0 +1,42 @@
+#ifndef  GRID_QCD_DHOP_H
+#define  GRID_QCD_DHOP_H
+
+namespace Grid {
+
+  namespace QCD {
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    // Helper classes that implement Wilson stencil for a single site.
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    
+    // Generic version works for any Nc and with extra flavour indices
+    class DiracOpt {
+    public:
+      // These ones will need to be package intelligently. WilsonType base class
+      // for use by DWF etc..
+      static void DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
+			   std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
+			   int sF,int sU,const LatticeFermion &in, LatticeFermion &out);
+      static void DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
+			      std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
+			      int sF,int sU,const LatticeFermion &in, LatticeFermion &out);
+
+    };
+
+    // Hand unrolled for Nc=3, one flavour
+    class DiracOptHand {
+    public:
+      // These ones will need to be package intelligently. WilsonType base class
+      // for use by DWF etc..
+      static void DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
+			   std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
+			   int sF,int sU,const LatticeFermion &in, LatticeFermion &out);
+      static void DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
+			      std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
+			      int sF,int sU,const LatticeFermion &in, LatticeFermion &out);
+
+    };
+
+  }
+}
+#endif

lib/qcd/action/fermion/WilsonKernelsHand.cc

@@ -27,7 +27,7 @@
     Chi_12 = ref()(1)(2);
 
 #define MULT_2SPIN(A)\
-   auto & ref(U._odata[ss](A));	\
+   auto & ref(U._odata[sU](A));	\
     U_00 = ref()(0,0);\
     U_10 = ref()(1,0);\
     U_20 = ref()(2,0);\
@@ -282,7 +282,7 @@ namespace QCD {
 
 void DiracOptHand::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
 			    std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
-			    int ss,const LatticeFermion &in, LatticeFermion &out)
+			    int sF,int sU,const LatticeFermion &in, LatticeFermion &out)
 {
   REGISTER vComplex result_00; // 12 regs on knc
   REGISTER vComplex result_01;
@@ -338,7 +338,8 @@ void DiracOptHand::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
 
 
   int offset,local,perm, ptype;
-
+  int ss=sF;
+  
   // Xp
   offset = st._offsets [Xp][ss];
   local  = st._is_local[Xp][ss];
@@ -514,24 +515,24 @@ void DiracOptHand::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
 
   {
     vSpinColourVector & ref (out._odata[ss]);
-    vstream(ref()(0)(0),result_00);
+    vstream(ref()(0)(0),result_00*(-0.5));
-    vstream(ref()(0)(1),result_01);
+    vstream(ref()(0)(1),result_01*(-0.5));
-    vstream(ref()(0)(2),result_02);
+    vstream(ref()(0)(2),result_02*(-0.5));
-    vstream(ref()(1)(0),result_10);
+    vstream(ref()(1)(0),result_10*(-0.5));
-    vstream(ref()(1)(1),result_11);
+    vstream(ref()(1)(1),result_11*(-0.5));
-    vstream(ref()(1)(2),result_12);
+    vstream(ref()(1)(2),result_12*(-0.5));
-    vstream(ref()(2)(0),result_20);
+    vstream(ref()(2)(0),result_20*(-0.5));
-    vstream(ref()(2)(1),result_21);
+    vstream(ref()(2)(1),result_21*(-0.5));
-    vstream(ref()(2)(2),result_22);
+    vstream(ref()(2)(2),result_22*(-0.5));
-    vstream(ref()(3)(0),result_30);
+    vstream(ref()(3)(0),result_30*(-0.5));
-    vstream(ref()(3)(1),result_31);
+    vstream(ref()(3)(1),result_31*(-0.5));
-    vstream(ref()(3)(2),result_32);
+    vstream(ref()(3)(2),result_32*(-0.5));
   }
 }
 
 void DiracOptHand::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
 			       std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> >  &buf,
-			       int ss,const LatticeFermion &in, LatticeFermion &out)
+			       int ss,int sU,const LatticeFermion &in, LatticeFermion &out)
 {
   REGISTER vComplex result_00; // 12 regs on knc
   REGISTER vComplex result_01;
@@ -752,18 +753,18 @@ void DiracOptHand::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
 
   {
     vSpinColourVector & ref (out._odata[ss]);
-    vstream(ref()(0)(0),result_00);
+    vstream(ref()(0)(0),result_00*(-0.5));
-    vstream(ref()(0)(1),result_01);
+    vstream(ref()(0)(1),result_01*(-0.5));
-    vstream(ref()(0)(2),result_02);
+    vstream(ref()(0)(2),result_02*(-0.5));
-    vstream(ref()(1)(0),result_10);
+    vstream(ref()(1)(0),result_10*(-0.5));
-    vstream(ref()(1)(1),result_11);
+    vstream(ref()(1)(1),result_11*(-0.5));
-    vstream(ref()(1)(2),result_12);
+    vstream(ref()(1)(2),result_12*(-0.5));
-    vstream(ref()(2)(0),result_20);
+    vstream(ref()(2)(0),result_20*(-0.5));
-    vstream(ref()(2)(1),result_21);
+    vstream(ref()(2)(1),result_21*(-0.5));
-    vstream(ref()(2)(2),result_22);
+    vstream(ref()(2)(2),result_22*(-0.5));
-    vstream(ref()(3)(0),result_30);
+    vstream(ref()(3)(0),result_30*(-0.5));
-    vstream(ref()(3)(1),result_31);
+    vstream(ref()(3)(1),result_31*(-0.5));
-    vstream(ref()(3)(2),result_32);
+    vstream(ref()(3)(2),result_32*(-0.5));
   }
 }
 }}