* Finished the template/policy style introduction of gparity, except the gparity force terms.

So valence sector looks ok.

FermionOperatorImpl.h provides the policy classes.

Expect HMC will introduce a smearing policy and a fermion representation change policy template
param. Will also probably need multi-precision work.

* HMC is running even-odd and non-checkerboarded (checked 4^4 wilson fermion/wilson gauge).

There appears to be a bug in the multi-level integrator -- <e-dH> passes with single level but
not with multi-level.

In any case there looks to be quite a bit to clean up.

This is the "const det" style implementation that is not appropriate  yet for clover since
it assumes that Mee is indept of the gauge fields. Easily fixed in future.

lib/Make.inc

@@ -1,4 +1,4 @@
 
-HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Config.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./MacroMagic.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/NerscIO.h ./qcd/action/ActionBase.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionImplTypedefs.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h
+HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Config.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./MacroMagic.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/NerscIO.h ./qcd/action/ActionBase.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeMap.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h
 
 CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./Init.cc ./Log.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/hmc/HMC.cc ./qcd/hmc/integrators/Integrator.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

lib/qcd/action/Actions.h

@@ -27,6 +27,8 @@
 // Utility functions
 ////////////////////////////////////////////
 #include <qcd/action/fermion/WilsonCompressor.h>     //used by all wilson type fermions
+#include <qcd/action/fermion/FermionOperatorImpl.h>
+#include <qcd/action/fermion/FermionOperator.h>
 #include <qcd/action/fermion/WilsonKernels.h>        //used by all wilson type fermions
 
 
@@ -54,8 +56,6 @@
 ////////////////////////////////////////////
 // Fermion operators / actions
 ////////////////////////////////////////////
-#include <qcd/action/fermion/FermionOperatorImpl.h>
-#include <qcd/action/fermion/FermionOperator.h>
 
 #include <qcd/action/fermion/WilsonFermion.h>       // 4d wilson like
 //#include <qcd/action/fermion/CloverFermion.h>

lib/qcd/action/fermion/CayleyFermion5D.cc

@@ -8,12 +8,12 @@ namespace QCD {
 					GridRedBlackCartesian &FiveDimRedBlackGrid,
 					GridCartesian         &FourDimGrid,
 					GridRedBlackCartesian &FourDimRedBlackGrid,
-					RealD _mass,RealD _M5) :
+					RealD _mass,RealD _M5,const ImplParams &p) :
    WilsonFermion5D<Impl>(_Umu,
 		   FiveDimGrid,
 		   FiveDimRedBlackGrid,
 		   FourDimGrid,
-		   FourDimRedBlackGrid,_M5),
+ 	 	   FourDimRedBlackGrid,_M5,p),
    mass(_mass)
  {
  }

lib/qcd/action/fermion/CayleyFermion5D.h

@@ -9,7 +9,7 @@ namespace Grid {
     class CayleyFermion5D : public WilsonFermion5D<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       // override multiply
@@ -64,7 +64,7 @@ namespace Grid {
 		      GridRedBlackCartesian &FiveDimRedBlackGrid,
 		      GridCartesian         &FourDimGrid,
 		      GridRedBlackCartesian &FourDimRedBlackGrid,
-		      RealD _mass,RealD _M5);
+		      RealD _mass,RealD _M5,const ImplParams &p= ImplParams());
 
     protected:
       void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c);

lib/qcd/action/fermion/ContinuedFractionFermion5D.cc

@@ -256,10 +256,10 @@ namespace Grid {
 							   GridRedBlackCartesian &FiveDimRedBlackGrid,
 							   GridCartesian         &FourDimGrid,
 							   GridRedBlackCartesian &FourDimRedBlackGrid,
-							   RealD _mass,RealD M5) :
+							   RealD _mass,RealD M5,const ImplParams &p) :
       WilsonFermion5D<Impl>(_Umu,
-		      FiveDimGrid, FiveDimRedBlackGrid,
+			    FiveDimGrid, FiveDimRedBlackGrid,
-		      FourDimGrid, FourDimRedBlackGrid,M5),
+			    FourDimGrid, FourDimRedBlackGrid,M5,p),
       mass(_mass)
     {
       int Ls = this->Ls;

lib/qcd/action/fermion/ContinuedFractionFermion5D.h

@@ -9,7 +9,7 @@ namespace Grid {
     class ContinuedFractionFermion5D : public WilsonFermion5D<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       // override multiply
@@ -41,7 +41,7 @@ namespace Grid {
 				 GridRedBlackCartesian &FiveDimRedBlackGrid,
 				 GridCartesian         &FourDimGrid,
 				 GridRedBlackCartesian &FourDimRedBlackGrid,
-				 RealD _mass,RealD M5);
+				 RealD _mass,RealD M5,const ImplParams &p= ImplParams());
 
     protected:
 

lib/qcd/action/fermion/DomainWallFermion.h

@@ -11,7 +11,7 @@ namespace Grid {
     class DomainWallFermion : public CayleyFermion5D<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       virtual void   Instantiatable(void) {};
@@ -21,13 +21,13 @@ namespace Grid {
 			GridRedBlackCartesian &FiveDimRedBlackGrid,
 			GridCartesian         &FourDimGrid,
 			GridRedBlackCartesian &FourDimRedBlackGrid,
-			RealD _mass,RealD _M5) : 
+			RealD _mass,RealD _M5,const ImplParams &p= ImplParams()) : 
-
+      
       CayleyFermion5D<Impl>(_Umu,
 			    FiveDimGrid,
 			    FiveDimRedBlackGrid,
 			    FourDimGrid,
-			    FourDimRedBlackGrid,_mass,_M5)
+			    FourDimRedBlackGrid,_mass,_M5,p)
 
       {
 	RealD eps = 1.0;

lib/qcd/action/fermion/FermionImplTypedefs.h

@@ -1,11 +0,0 @@
-    public:
-      typedef typename Impl::Simd                           Simd;
-      typedef typename Impl::FermionField           FermionField;
-      typedef typename Impl::GaugeLinkField       GaugeLinkField;
-      typedef typename Impl::GaugeField               GaugeField;
-      typedef typename Impl::DoubledGaugeField DoubledGaugeField;
-      typedef typename Impl::SiteSpinor               SiteSpinor;
-      typedef typename Impl::SiteHalfSpinor       SiteHalfSpinor;
-      typedef typename Impl::Compressor               Compressor;
-      typedef WilsonKernels<Impl> Kernels;
-

lib/qcd/action/fermion/FermionOperator.h

@@ -21,7 +21,10 @@ namespace Grid {
     class FermionOperator : public CheckerBoardedSparseMatrixBase<typename Impl::FermionField>, public Impl
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+
+      INHERIT_IMPL_TYPES(Impl);
+
+      FermionOperator(const ImplParams &p= ImplParams()) : Impl(p) {};
 
     public:
 

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -5,35 +5,136 @@ namespace Grid {
 
   namespace QCD {
 
+    //////////////////////////////////////////////
+    // Template parameter class constructs to package
+    // externally control Fermion implementations
+    // in orthogonal directions
+    //
+    // Ultimately need Impl to always define types where XXX is opaque
+    //
+    //    typedef typename XXX               Simd;
+    //    typedef typename XXX     GaugeLinkField;	
+    //    typedef typename XXX         GaugeField;
+    //    typedef typename XXX      GaugeActField;
+    //    typedef typename XXX       FermionField;
+    //    typedef typename XXX  DoubledGaugeField;
+    //    typedef typename XXX         SiteSpinor;
+    //    typedef typename XXX     SiteHalfSpinor;	
+    //    typedef typename XXX         Compressor;	
+    //
+    // and Methods:
+    //    void ImportGauge(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
+    //    void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
+    //    void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,CartesianStencil &St)
+    //    void InsertForce4D(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu)
+    //    void InsertForce5D(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu)
+    //
+    //
+    // To acquire the typedefs from "Base" (either a base class or template param) use:
+    //
+    // INHERIT_GIMPL_TYPES(Base)
+    // INHERIT_FIMPL_TYPES(Base)
+    // INHERIT_IMPL_TYPES(Base)
+    //
+    // The Fermion operators will do the following:
+    //
+    // struct MyOpParams { 
+    //   RealD mass;
+    // };
+    //
+    //
+    // template<class Impl>
+    // class MyOp : pubic<Impl> { 
+    // public:
+    //
+    //    INHERIT_ALL_IMPL_TYPES(Impl);
+    //
+    //    MyOp(MyOpParams Myparm, ImplParams &ImplParam) :  Impl(ImplParam)
+    //    {
+    //
+    //    };
+    //    
+    //  }
+    //////////////////////////////////////////////
+
+    
+    ////////////////////////////////////////////////////////////////////////
+    // Implementation dependent gauge types
+    ////////////////////////////////////////////////////////////////////////
+
+#define INHERIT_IMPL_TYPES(Base) \
+    INHERIT_GIMPL_TYPES(Base)\
+    INHERIT_FIMPL_TYPES(Base)
+
+#define INHERIT_GIMPL_TYPES(GImpl) \
+    typedef typename GImpl::Simd                           Simd;\
+    typedef typename GImpl::GaugeLinkField       GaugeLinkField;\
+    typedef typename GImpl::GaugeField               GaugeField;	
+    
+    // Composition with smeared link, bc's etc.. probably need multiple inheritance
     // Variable precision "S" and variable Nc
     template<class S,int Nrepresentation=Nc>
-    class WilsonImpl { 
+    class ImplGauge { 
+    public:
+    
+      typedef S Simd;
+    
+      template<typename vtype> using iImplGaugeLink          = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
+      template<typename vtype> using iImplGaugeField         = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd  >;
+    
+      typedef iImplGaugeLink    <Simd>           SiteGaugeLink;
+      typedef iImplGaugeField   <Simd>           SiteGaugeField;
+    
+      typedef Lattice<SiteGaugeLink>                GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly
+      typedef Lattice<SiteGaugeField>                   GaugeField;
+
+    };
+
+    ////////////////////////////////////////////////////////////////////////
+    // Implementation dependent fermion types
+    ////////////////////////////////////////////////////////////////////////
+
+#define INHERIT_FIMPL_TYPES(Impl)\
+    typedef typename Impl::FermionField           FermionField;		\
+    typedef typename Impl::DoubledGaugeField DoubledGaugeField;		\
+    typedef typename Impl::SiteSpinor               SiteSpinor;		\
+    typedef typename Impl::SiteHalfSpinor       SiteHalfSpinor;		\
+    typedef typename Impl::Compressor               Compressor;		\
+    typedef typename Impl::ImplParams ImplParams;
+
+    ///////
+    // Single flavour four spinors with colour index
+    ///////
+    template<class S,int Nrepresentation=Nc>
+    class WilsonImpl :  public ImplGauge<S,Nrepresentation> { 
     public:
 
-      typedef S Simd;
+      typedef ImplGauge<S,Nrepresentation> Gimpl;
+
+
+      INHERIT_GIMPL_TYPES(Gimpl);
 
       template<typename vtype> using iImplSpinor             = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
       template<typename vtype> using iImplHalfSpinor         = iScalar<iVector<iVector<vtype, Nrepresentation>, Nhs> >;
-      template<typename vtype> using iImplGaugeLink          = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
-      template<typename vtype> using iImplGaugeField         = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd  >;
       template<typename vtype> using iImplDoubledGaugeField  = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds >;
     
       typedef iImplSpinor    <Simd>           SiteSpinor;
       typedef iImplHalfSpinor<Simd>           SiteHalfSpinor;
-      typedef iImplGaugeLink <Simd>           SiteGaugeLink;
-      typedef iImplGaugeField<Simd>           SiteGaugeField;
       typedef iImplDoubledGaugeField<Simd>    SiteDoubledGaugeField;
 
       typedef Lattice<SiteSpinor>                 FermionField;
-      typedef Lattice<SiteGaugeLink>            GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly
-      typedef Lattice<SiteGaugeField>               GaugeField;
       typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
 
       typedef WilsonCompressor<SiteHalfSpinor,SiteSpinor> Compressor;
 
-      static inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,CartesianStencil &St){
+      typedef struct WilsonImplParams { } ImplParams;
+      ImplParams Params;
+      WilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {}; 
+
+      inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,CartesianStencil &St){
         mult(&phi(),&U(mu),&chi());
       }
+
       inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
       {
         conformable(Uds._grid,GaugeGrid);
@@ -46,52 +147,69 @@ namespace Grid {
 	  PokeIndex<LorentzIndex>(Uds,U,mu+4);
 	}
       }
-
+      
-      inline void InsertForce(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu){
+      inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){
 	GaugeLinkField link(mat._grid);
 	link = TraceIndex<SpinIndex>(outerProduct(Btilde,A)); 
 	PokeIndex<LorentzIndex>(mat,link,mu);
       }   
 
+      inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){
+
+	int Ls=Btilde._grid->_fdimensions[0];
+
+	GaugeLinkField tmp(mat._grid);
+	tmp = zero;
+PARALLEL_FOR_LOOP
+	for(int sss=0;sss<tmp._grid->oSites();sss++){
+	  int sU=sss;
+	  for(int s=0;s<Ls;s++){
+	    int sF = s+Ls*sU;
+	    tmp[sU] = tmp[sU]+ traceIndex<SpinIndex>(outerProduct(Btilde[sF],Atilde[sF])); // ordering here
+	  }
+	}
+	PokeIndex<LorentzIndex>(mat,tmp,mu);
+	
+      }
+
     };
 
-    template<class S,int Nrepresentation=Nc>
+    ////////////////////////////////////////////////////////////////////////////////////////
-    class GparityWilsonImpl { 
+    // Flavour doubled spinors; is Gparity the only? what about C*?
+    ////////////////////////////////////////////////////////////////////////////////////////
+    template<class S,int Nrepresentation>
+    class GparityWilsonImpl : public ImplGauge<S,Nrepresentation> { 
     public:
 
-      typedef S Simd;
+      typedef ImplGauge<S,Nrepresentation> Gimpl;
+
+      INHERIT_GIMPL_TYPES(Gimpl);
 
       template<typename vtype> using iImplSpinor             = iVector<iVector<iVector<vtype, Nrepresentation>, Ns>, Ngp >;
       template<typename vtype> using iImplHalfSpinor         = iVector<iVector<iVector<vtype, Nrepresentation>, Nhs>, Ngp >;
-      template<typename vtype> using iImplGaugeField         = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd  >;
-
-      template<typename vtype> using iImplGaugeLink          = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
       template<typename vtype> using iImplDoubledGaugeField  = iVector<iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds >, Ngp >;
     
       typedef iImplSpinor    <Simd>           SiteSpinor;
       typedef iImplHalfSpinor<Simd>           SiteHalfSpinor;
-      typedef iImplGaugeLink <Simd>           SiteGaugeLink;
-      typedef iImplGaugeField<Simd>           SiteGaugeField;
       typedef iImplDoubledGaugeField<Simd>    SiteDoubledGaugeField;
 
       typedef Lattice<SiteSpinor>                 FermionField;
-      typedef Lattice<SiteGaugeLink>            GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly
-      typedef Lattice<SiteGaugeField>               GaugeField;
       typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
 
-      //      typedef GparityWilsonCompressor<SiteHalfSpinor,SiteSpinor> Compressor;
       typedef WilsonCompressor<SiteHalfSpinor,SiteSpinor> Compressor;
 
-      // provide the multiply by link that is differentiated between Gparity (with flavour index) and 
+      typedef struct GparityWilsonImplParams {std::vector<int> twists; } ImplParams;
-      // non-Gparity
+      ImplParams Params;
-      static inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,CartesianStencil &St){
+      GparityWilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {}; 
+      
+      // provide the multiply by link that is differentiated between Gparity (with flavour index) and non-Gparity
+      inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,CartesianStencil &St){
 
 	typedef SiteHalfSpinor vobj;
 	typedef typename SiteHalfSpinor::scalar_object sobj;
 
 	vobj vtmp;
 	sobj stmp;
-	std::vector<int> gpbc({0,0,0,1,0,0,0,1});
 	
 	GridBase *grid = St._grid;
       
@@ -99,125 +217,131 @@ namespace Grid {
 	
 	int direction    = St._directions[mu];
 	int distance     = St._distances[mu];
-	int ptype     = St._permute_type[mu]; 
+	int ptype        = St._permute_type[mu]; 
-	int sl        = St._grid->_simd_layout[direction];
+	int sl           = St._grid->_simd_layout[direction];
-	
+
+	// Fixme X.Y.Z.T hardcode in stencil
+	int mmu          = mu % Nd;
+
 	// assert our assumptions
 	assert((distance==1)||(distance==-1)); // nearest neighbour stencil hard code
 	assert((sl==1)||(sl==2));
 	
 	std::vector<int> icoor;
       
-	if ( SE->_around_the_world && gpbc[mu] ) {
+	if ( SE->_around_the_world && Params.twists[mmu] ) {
 
 	  if ( sl == 2 ) {
 
-	  std::vector<sobj> vals(Nsimd);
+	    std::vector<sobj> vals(Nsimd);
 
-	  extract(chi,vals);
+	    extract(chi,vals);
-	  for(int s=0;s<Nsimd;s++){
+	    for(int s=0;s<Nsimd;s++){
 
-	    grid->iCoorFromIindex(icoor,s);
+	      grid->iCoorFromIindex(icoor,s);
-
+	      
-	    assert((icoor[direction]==0)||(icoor[direction]==1));
+	      assert((icoor[direction]==0)||(icoor[direction]==1));
-
+	      
-	    int permute_lane;
+	      int permute_lane;
-	    if ( distance == 1) {
+	      if ( distance == 1) {
-	      permute_lane = icoor[direction]?1:0;
+		permute_lane = icoor[direction]?1:0;
-	    } else {
+	      } else {
-	      permute_lane = icoor[direction]?0:1;
+		permute_lane = icoor[direction]?0:1;
+	      }
+	      
+	      if ( permute_lane ) { 
+		stmp(0) = vals[s](1);
+		stmp(1) = vals[s](0);
+		vals[s] = stmp;
+	      }
 	    }
+	    merge(vtmp,vals);
 
-            if ( permute_lane ) { 
+	  } else { 
-	      stmp(0) = vals[s](1);
+	    vtmp(0) = chi(1);
-	      stmp(1) = vals[s](0);
+	    vtmp(1) = chi(0);
-	      vals[s] = stmp;
-	    }
 	  }
-	  merge(vtmp,vals);
+	  mult(&phi(0),&U(0)(mu),&vtmp(0));
-
+	  mult(&phi(1),&U(1)(mu),&vtmp(1));
+	  
 	} else { 
-	  vtmp(0) = chi(1);
+	  mult(&phi(0),&U(0)(mu),&chi(0));
-	  vtmp(1) = chi(0);
+	  mult(&phi(1),&U(1)(mu),&chi(1));
 	}
-	mult(&phi(0),&U(0)(mu),&vtmp(0));
-	mult(&phi(1),&U(1)(mu),&vtmp(1));
 	
-      } else { 
-	mult(&phi(0),&U(0)(mu),&chi(0));
-	mult(&phi(1),&U(1)(mu),&chi(1));
       }
 
-    }
+      inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
-
+      {
-    static inline void InsertForce(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu){
-      // Fixme
-      return;
-    }
-    static inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
-    {
-
-      conformable(Uds._grid,GaugeGrid);
-      conformable(Umu._grid,GaugeGrid);
-
-      GaugeLinkField Utmp(GaugeGrid);
-      GaugeLinkField U(GaugeGrid);
-      GaugeLinkField Uconj(GaugeGrid);
-
-      Lattice<iScalar<vInteger> > coor(GaugeGrid);
-      
-      std::vector<int> gpdirs({0,0,0,1});
-      
-      for(int mu=0;mu<Nd;mu++){
 	
-	LatticeCoordinate(coor,mu);
+	conformable(Uds._grid,GaugeGrid);
+	conformable(Umu._grid,GaugeGrid);
 	
-	U     = PeekIndex<LorentzIndex>(Umu,mu);
+	GaugeLinkField Utmp(GaugeGrid);
-	Uconj = conjugate(U);
+	GaugeLinkField U(GaugeGrid);
-
+	GaugeLinkField Uconj(GaugeGrid);
-	int neglink = GaugeGrid->GlobalDimensions()[mu]-1;
 	
-	if ( gpdirs[mu] ) { 
+	Lattice<iScalar<vInteger> > coor(GaugeGrid);
-	  Uconj = where(coor==neglink,-Uconj,Uconj);
-	}
 
+	
+	for(int mu=0;mu<Nd;mu++){
+	  
+	  LatticeCoordinate(coor,mu);
+	  
+	  U     = PeekIndex<LorentzIndex>(Umu,mu);
+	  Uconj = conjugate(U);
+
+	  // This phase could come from a simple bc 1,1,-1,1 ..
+	  int neglink = GaugeGrid->GlobalDimensions()[mu]-1;
+	  if ( Params.twists[mu] ) { 
+	    Uconj = where(coor==neglink,-Uconj,Uconj);
+	  }
+
+	  
 PARALLEL_FOR_LOOP
-        for(auto ss=U.begin();ss<U.end();ss++){
+	  for(auto ss=U.begin();ss<U.end();ss++){
-	  Uds[ss](0)(mu) = U[ss]();
+	    Uds[ss](0)(mu) = U[ss]();
-	  Uds[ss](1)(mu) = Uconj[ss]();
+	    Uds[ss](1)(mu) = Uconj[ss]();
-	}
+	  }
- 
+	  
-        U     = adj(Cshift(U    ,mu,-1));      // correct except for spanning the boundary
+	  U     = adj(Cshift(U    ,mu,-1));      // correct except for spanning the boundary
-        Uconj = adj(Cshift(Uconj,mu,-1));
+	  Uconj = adj(Cshift(Uconj,mu,-1));
-
+	  
-	Utmp = U;
+	  Utmp = U;
-	if ( gpdirs[mu] ) { 
+	  if ( Params.twists[mu] ) { 
-	  Utmp = where(coor==0,Uconj,Utmp);
+	    Utmp = where(coor==0,Uconj,Utmp);
-	}
+	  }
-
+	  
 PARALLEL_FOR_LOOP
-        for(auto ss=U.begin();ss<U.end();ss++){
+	  for(auto ss=U.begin();ss<U.end();ss++){
-	  Uds[ss](0)(mu+4) = Utmp[ss]();
+	    Uds[ss](0)(mu+4) = Utmp[ss]();
-	}
+	  }
-
+	  
-        Utmp = Uconj;
+	  Utmp = Uconj;
-	if ( gpdirs[mu] ) { 
+	  if ( Params.twists[mu] ) { 
-	  Utmp = where(coor==0,U,Utmp);
+	    Utmp = where(coor==0,U,Utmp);
-	}
+	  }
-
+	  
 PARALLEL_FOR_LOOP
-	for(auto ss=U.begin();ss<U.end();ss++){
+	  for(auto ss=U.begin();ss<U.end();ss++){
-	  Uds[ss](1)(mu+4) = Utmp[ss]();
+	    Uds[ss](1)(mu+4) = Utmp[ss]();
+	  }
+	  
 	}
-
       }
-    }
 
+      inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){
+	assert(0);
+	// Fixme
+	return;
+      }
+      inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){
+	assert(0);
+	// Fixme
+	return;
+      }
     };
 
-    // Could QCD, SU2, QED etc...
-    // Fund, adj, etc...
-    // Composition with smeared link, bc's etc.. probably need multiple inheritance
     typedef WilsonImpl<vComplex ,Nc> WilsonImplR; // Real.. whichever prec
     typedef WilsonImpl<vComplexF,Nc> WilsonImplF; // Float
     typedef WilsonImpl<vComplexD,Nc> WilsonImplD; // Double
@@ -226,7 +350,6 @@ PARALLEL_FOR_LOOP
     typedef GparityWilsonImpl<vComplexF,Nc> GparityWilsonImplF; // Float
     typedef GparityWilsonImpl<vComplexD,Nc> GparityWilsonImplD; // Double
 
-
   }
 }
 #endif

lib/qcd/action/fermion/MobiusFermion.h

@@ -11,7 +11,7 @@ namespace Grid {
     class MobiusFermion : public CayleyFermion5D<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       virtual void   Instantiatable(void) {};
@@ -22,13 +22,13 @@ namespace Grid {
 		    GridCartesian         &FourDimGrid,
 		    GridRedBlackCartesian &FourDimRedBlackGrid,
 		    RealD _mass,RealD _M5,
-		    RealD b, RealD c) : 
+		    RealD b, RealD c,const ImplParams &p= ImplParams()) : 
       
       CayleyFermion5D<Impl>(_Umu,
-		      FiveDimGrid,
+			    FiveDimGrid,
-		      FiveDimRedBlackGrid,
+			    FiveDimRedBlackGrid,
-		      FourDimGrid,
+			    FourDimGrid,
-		      FourDimRedBlackGrid,_mass,_M5)
+			    FourDimRedBlackGrid,_mass,_M5,p)
 
       {
 	RealD eps = 1.0;

lib/qcd/action/fermion/MobiusZolotarevFermion.h

@@ -11,7 +11,7 @@ namespace Grid {
     class MobiusZolotarevFermion : public CayleyFermion5D<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       virtual void   Instantiatable(void) {};
@@ -23,13 +23,13 @@ namespace Grid {
 			      GridRedBlackCartesian &FourDimRedBlackGrid,
 			      RealD _mass,RealD _M5,
 			      RealD b, RealD c,
-			      RealD lo, RealD hi) : 
+			      RealD lo, RealD hi,const ImplParams &p= ImplParams()) : 
       
       CayleyFermion5D<Impl>(_Umu,
-		      FiveDimGrid,
+			    FiveDimGrid,
-		      FiveDimRedBlackGrid,
+			    FiveDimRedBlackGrid,
-		      FourDimGrid,
+			    FourDimGrid,
-		      FourDimRedBlackGrid,_mass,_M5)
+			    FourDimRedBlackGrid,_mass,_M5,p)
 
       {
 	RealD eps = lo/hi;

lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h

@@ -11,7 +11,7 @@ namespace Grid {
     class OverlapWilsonCayleyTanhFermion : public MobiusFermion<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       // Constructors
@@ -21,14 +21,14 @@ namespace Grid {
 				   GridCartesian         &FourDimGrid,
 				   GridRedBlackCartesian &FourDimRedBlackGrid,
 				   RealD _mass,RealD _M5,
-				   RealD scale) :
+				   RealD scale,const ImplParams &p= ImplParams()) :
       
       // b+c=scale, b-c = 0 <=> b =c = scale/2
       MobiusFermion<Impl>(_Umu,
 			  FiveDimGrid,
 			  FiveDimRedBlackGrid,
 			  FourDimGrid,
-			  FourDimRedBlackGrid,_mass,_M5,0.5*scale,0.5*scale)
+			  FourDimRedBlackGrid,_mass,_M5,0.5*scale,0.5*scale,p)
 	{
 	}
     };

lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h

@@ -11,7 +11,7 @@ namespace Grid {
     class OverlapWilsonCayleyZolotarevFermion : public MobiusZolotarevFermion<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       // Constructors
@@ -22,13 +22,13 @@ namespace Grid {
 					GridCartesian         &FourDimGrid,
 					GridRedBlackCartesian &FourDimRedBlackGrid,
 					RealD _mass,RealD _M5,
-					RealD lo, RealD hi) : 
+					RealD lo, RealD hi,const ImplParams &p= ImplParams()) : 
       // b+c=1.0, b-c = 0 <=> b =c = 1/2
       MobiusZolotarevFermion<Impl>(_Umu,
-			     FiveDimGrid,
+				   FiveDimGrid,
-			     FiveDimRedBlackGrid,
+				   FiveDimRedBlackGrid,
-			     FourDimGrid,
+				   FourDimGrid,
-			     FourDimRedBlackGrid,_mass,_M5,0.5,0.5,lo,hi)
+				   FourDimRedBlackGrid,_mass,_M5,0.5,0.5,lo,hi,p)
 
       {}
 

lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h

@@ -11,7 +11,7 @@ namespace Grid {
     class OverlapWilsonContFracTanhFermion : public ContinuedFractionFermion5D<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       virtual void   Instantiatable(void){};
@@ -22,14 +22,14 @@ namespace Grid {
 				     GridCartesian         &FourDimGrid,
 				     GridRedBlackCartesian &FourDimRedBlackGrid,
 				     RealD _mass,RealD _M5,
-				     RealD scale) :
+				     RealD scale,const ImplParams &p= ImplParams()) :
       
       // b+c=scale, b-c = 0 <=> b =c = scale/2
       ContinuedFractionFermion5D<Impl>(_Umu,
-				 FiveDimGrid,
+				       FiveDimGrid,
-				 FiveDimRedBlackGrid,
+				       FiveDimRedBlackGrid,
-				 FourDimGrid,
+				       FourDimGrid,
-				 FourDimRedBlackGrid,_mass,_M5)
+				       FourDimRedBlackGrid,_mass,_M5,p)
 	{
 	  assert((this->Ls&0x1)==1); // Odd Ls required
 	  int nrational=this->Ls-1;// Even rational order

lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h

@@ -11,7 +11,7 @@ namespace Grid {
     class OverlapWilsonContFracZolotarevFermion : public ContinuedFractionFermion5D<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       virtual void   Instantiatable(void){};
@@ -22,14 +22,14 @@ namespace Grid {
 					  GridCartesian         &FourDimGrid,
 					  GridRedBlackCartesian &FourDimRedBlackGrid,
 					  RealD _mass,RealD _M5,
-					  RealD lo,RealD hi):
+					  RealD lo,RealD hi,const ImplParams &p= ImplParams()):
       
       // b+c=scale, b-c = 0 <=> b =c = scale/2
       ContinuedFractionFermion5D<Impl>(_Umu,
-				 FiveDimGrid,
+				       FiveDimGrid,
-				 FiveDimRedBlackGrid,
+				       FiveDimRedBlackGrid,
-				 FourDimGrid,
+				       FourDimGrid,
-				 FourDimRedBlackGrid,_mass,_M5)
+				       FourDimRedBlackGrid,_mass,_M5,p)
 	{
 	  assert((this->Ls&0x1)==1); // Odd Ls required
 

lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h

@@ -11,7 +11,7 @@ namespace Grid {
     class OverlapWilsonPartialFractionTanhFermion : public PartialFractionFermion5D<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       virtual void   Instantiatable(void){};
@@ -22,14 +22,14 @@ namespace Grid {
 					    GridCartesian         &FourDimGrid,
 					    GridRedBlackCartesian &FourDimRedBlackGrid,
 					    RealD _mass,RealD _M5,
-					    RealD scale) :
+					    RealD scale,const ImplParams &p= ImplParams()) :
       
       // b+c=scale, b-c = 0 <=> b =c = scale/2
       PartialFractionFermion5D<Impl>(_Umu,
 				     FiveDimGrid,
 				     FiveDimRedBlackGrid,
 				     FourDimGrid,
-				     FourDimRedBlackGrid,_mass,_M5)
+				     FourDimRedBlackGrid,_mass,_M5,p)
 	{
 	  assert((this->Ls&0x1)==1); // Odd Ls required
 	  int nrational=this->Ls-1;// Even rational order

lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h

@@ -11,25 +11,25 @@ namespace Grid {
     class OverlapWilsonPartialFractionZolotarevFermion : public PartialFractionFermion5D<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       virtual void   Instantiatable(void){};
       // Constructors
     OverlapWilsonPartialFractionZolotarevFermion(GaugeField &_Umu,
-					  GridCartesian         &FiveDimGrid,
+						 GridCartesian         &FiveDimGrid,
-					  GridRedBlackCartesian &FiveDimRedBlackGrid,
+						 GridRedBlackCartesian &FiveDimRedBlackGrid,
-					  GridCartesian         &FourDimGrid,
+						 GridCartesian         &FourDimGrid,
-					  GridRedBlackCartesian &FourDimRedBlackGrid,
+						 GridRedBlackCartesian &FourDimRedBlackGrid,
-					  RealD _mass,RealD _M5,
+						 RealD _mass,RealD _M5,
-					  RealD lo,RealD hi):
+						 RealD lo,RealD hi,const ImplParams &p= ImplParams()):
       
       // b+c=scale, b-c = 0 <=> b =c = scale/2
       PartialFractionFermion5D<Impl>(_Umu,
 				     FiveDimGrid,
 				     FiveDimRedBlackGrid,
 				     FourDimGrid,
-				     FourDimRedBlackGrid,_mass,_M5)
+				     FourDimRedBlackGrid,_mass,_M5,p)
 	{
 	  assert((this->Ls&0x1)==1); // Odd Ls required
 

lib/qcd/action/fermion/PartialFractionFermion5D.cc

@@ -369,14 +369,15 @@ namespace Grid {
       // Constructors
     template<class Impl>
     PartialFractionFermion5D<Impl>::PartialFractionFermion5D(GaugeField &_Umu,
-						       GridCartesian         &FiveDimGrid,
+							     GridCartesian         &FiveDimGrid,
-						       GridRedBlackCartesian &FiveDimRedBlackGrid,
+							     GridRedBlackCartesian &FiveDimRedBlackGrid,
-						       GridCartesian         &FourDimGrid,
+							     GridCartesian         &FourDimGrid,
-						       GridRedBlackCartesian &FourDimRedBlackGrid,
+							     GridRedBlackCartesian &FourDimRedBlackGrid,
-						       RealD _mass,RealD M5) :
+							     RealD _mass,RealD M5,
+							     const ImplParams &p) :
       WilsonFermion5D<Impl>(_Umu,
-		      FiveDimGrid, FiveDimRedBlackGrid,
+			    FiveDimGrid, FiveDimRedBlackGrid,
-		      FourDimGrid, FourDimRedBlackGrid,M5),
+			    FourDimGrid, FourDimRedBlackGrid,M5,p),
       mass(_mass)
 
     {

lib/qcd/action/fermion/PartialFractionFermion5D.h

@@ -9,7 +9,7 @@ namespace Grid {
     class PartialFractionFermion5D : public WilsonFermion5D<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       const int part_frac_chroma_convention=1;
@@ -47,7 +47,7 @@ namespace Grid {
 			       GridRedBlackCartesian &FiveDimRedBlackGrid,
 			       GridCartesian         &FourDimGrid,
 			       GridRedBlackCartesian &FourDimRedBlackGrid,
-			       RealD _mass,RealD M5);
+			       RealD _mass,RealD M5,const ImplParams &p= ImplParams());
 
     protected:
 

lib/qcd/action/fermion/ScaledShamirFermion.h

@@ -11,7 +11,7 @@ namespace Grid {
     class ScaledShamirFermion : public MobiusFermion<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       // Constructors

lib/qcd/action/fermion/ShamirZolotarevFermion.h

@@ -11,7 +11,7 @@ namespace Grid {
     class ShamirZolotarevFermion : public MobiusZolotarevFermion<Impl>
     {
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
     public:
 
       // Constructors
@@ -23,14 +23,14 @@ namespace Grid {
 			   GridCartesian         &FourDimGrid,
 			   GridRedBlackCartesian &FourDimRedBlackGrid,
 			   RealD _mass,RealD _M5,
-			   RealD lo, RealD hi) : 
+			   RealD lo, RealD hi,const ImplParams &p= ImplParams()) : 
       
       // b+c = 1; b-c = 1 => b=1, c=0
       MobiusZolotarevFermion<Impl>(_Umu,
-			     FiveDimGrid,
+				   FiveDimGrid,
-			     FiveDimRedBlackGrid,
+				   FiveDimRedBlackGrid,
-			     FourDimGrid,
+				   FourDimGrid,
-			     FourDimRedBlackGrid,_mass,_M5,1.0,0.0,lo,hi)
+				   FourDimRedBlackGrid,_mass,_M5,1.0,0.0,lo,hi,p)
       
       {}
 

lib/qcd/action/fermion/WilsonFermion.cc

@@ -15,7 +15,8 @@ namespace QCD {
   WilsonFermion<Impl>::WilsonFermion(GaugeField &_Umu,
 				     GridCartesian         &Fgrid,
 				     GridRedBlackCartesian &Hgrid, 
-				     RealD _mass) :
+				     RealD _mass,const ImplParams &p) :
+        Kernels(p),
         _grid(&Fgrid),
 	_cbgrid(&Hgrid),
 	Stencil    (&Fgrid,npoint,Even,directions,displacements),
@@ -118,7 +119,9 @@ namespace QCD {
     Compressor compressor(dag);
     
     FermionField Btilde(B._grid);
-    
+    FermionField Atilde(B._grid);
+    Atilde = A;
+
     st.HaloExchange<SiteSpinor,SiteHalfSpinor,Compressor>(B,comm_buf,compressor);
     
     for(int mu=0;mu<Nd;mu++){
@@ -140,7 +143,7 @@ PARALLEL_FOR_LOOP
       //////////////////////////////////////////////////
       // spin trace outer product
       //////////////////////////////////////////////////
-      Impl::InsertForce(mat,Btilde,A,mu);
+      Impl::InsertForce4D(mat,Btilde,Atilde,mu);
 
     }
   }

lib/qcd/action/fermion/WilsonFermion.h

@@ -15,11 +15,11 @@ namespace Grid {
     };
 
     template<class Impl>
-    class WilsonFermion : public FermionOperator<Impl>, public WilsonFermionStatic
+    class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic
     {
-#include <qcd/action/fermion/FermionImplTypedefs.h>
-
     public:
+    INHERIT_IMPL_TYPES(Impl);
+    typedef WilsonKernels<Impl> Kernels;
 
       ///////////////////////////////////////////////////////////////
       // Implement the abstract base
@@ -88,7 +88,9 @@ namespace Grid {
       WilsonFermion(GaugeField &_Umu,
 		    GridCartesian         &Fgrid,
 		    GridRedBlackCartesian &Hgrid, 
-		    RealD _mass) ;
+		    RealD _mass,
+		    const ImplParams &p= ImplParams()
+		    ) ;
 
       // DoubleStore impl dependent
       void ImportGauge(const GaugeField &_Umu);

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -15,7 +15,8 @@ WilsonFermion5D<Impl>::WilsonFermion5D(GaugeField &_Umu,
 				       GridRedBlackCartesian &FiveDimRedBlackGrid,
 				       GridCartesian         &FourDimGrid,
 				       GridRedBlackCartesian &FourDimRedBlackGrid,
-				       RealD _M5) :
+				       RealD _M5,const ImplParams &p) :
+  Kernels(p),
   _FiveDimGrid(&FiveDimGrid),
   _FiveDimRedBlackGrid(&FiveDimRedBlackGrid),
   _FourDimGrid(&FourDimGrid),
@@ -105,11 +106,11 @@ PARALLEL_FOR_LOOP
 
 template<class Impl>
 void WilsonFermion5D<Impl>::DerivInternal(CartesianStencil & st,
-				    DoubledGaugeField & U,
+					  DoubledGaugeField & U,
-				    GaugeField &mat,
+					  GaugeField &mat,
-				    const FermionField &A,
+					  const FermionField &A,
-				    const FermionField &B,
+					  const FermionField &B,
-				    int dag)
+					  int dag)
 {
   assert((dag==DaggerNo) ||(dag==DaggerYes));
 
@@ -118,7 +119,6 @@ void WilsonFermion5D<Impl>::DerivInternal(CartesianStencil & st,
 
   Compressor compressor(dag);
   
-  GaugeLinkField tmp(mat._grid);
   FermionField Btilde(B._grid);
   FermionField Atilde(B._grid);
 
@@ -157,14 +157,7 @@ PARALLEL_FOR_LOOP
 
     }
 
-    Impl::InsertForce(mat,Btilde,A,mu);
+    Impl::InsertForce5D(mat,Btilde,Atilde,mu);
-    /*
-    tmp = zero;
-    for(int sss=0;sss<U._grid->oSites();sss++){
-      tmp[sU] = tmp[sU]+ traceIndex<SpinIndex>(outerProduct(Btilde[sF],Atilde[sF])); // ordering here
-    }
-    PokeIndex<LorentzIndex>(mat,tmp,mu);
-    */
 
   }
 }

lib/qcd/action/fermion/WilsonFermion5D.h

@@ -26,9 +26,11 @@ namespace Grid {
     };
 
     template<class Impl>
-    class WilsonFermion5D : public FermionOperator<Impl>, public WilsonFermion5DStatic
+    class WilsonFermion5D : public WilsonKernels<Impl>, public WilsonFermion5DStatic
     {
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
+     typedef WilsonKernels<Impl> Kernels;
+
     public:
       ///////////////////////////////////////////////////////////////
       // Implement the abstract base
@@ -90,7 +92,7 @@ namespace Grid {
 		      GridRedBlackCartesian &FiveDimRedBlackGrid,
 		      GridCartesian         &FourDimGrid,
 		      GridRedBlackCartesian &FourDimRedBlackGrid,
-		      double _M5);
+		      double _M5,const ImplParams &p= ImplParams());
 
       // DoubleStore
       void ImportGauge(const GaugeField &_Umu);

lib/qcd/action/fermion/WilsonKernels.h

@@ -7,43 +7,41 @@ namespace Grid {
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     // Helper routines that implement Wilson stencil for a single site.
+    // Common to both the WilsonFermion and WilsonFermion5D
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-    template<class Impl> 
+    template<class Impl> class WilsonKernels : public FermionOperator<Impl> { 
-    class WilsonKernels { 
     public:
 
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+     INHERIT_IMPL_TYPES(Impl);
-
+     typedef FermionOperator<Impl> Base;
+     
     public:
-      static void 
+     void DiracOptDhopSite(CartesianStencil &st,DoubledGaugeField &U,
-	DiracOptDhopSite(CartesianStencil &st,DoubledGaugeField &U,
+			   std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-			 std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
+			   int sF,int sU,const FermionField &in, FermionField &out);
-			 int sF,int sU,const FermionField &in, FermionField &out);
       
-      static void 
+     void DiracOptDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U,
-	DiracOptDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U,
+			      std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-			    std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
+			      int sF,int sU,const FermionField &in,FermionField &out);
-			    int sF,int sU,const FermionField &in,FermionField &out);
 
-      static void 
+     void DiracOptDhopDir(CartesianStencil &st,DoubledGaugeField &U,
-	DiracOptDhopDir(CartesianStencil &st,DoubledGaugeField &U,
+			  std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-			std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
+			  int sF,int sU,const FermionField &in, FermionField &out,int dirdisp,int gamma);
-			int sF,int sU,const FermionField &in, FermionField &out,int dirdisp,int gamma);
 
-      static void 
+     void DiracOptHandDhopSite(CartesianStencil &st,DoubledGaugeField &U,
-	DiracOptHandDhopSite(CartesianStencil &st,DoubledGaugeField &U,
+			       std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-			     std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
+			       int sF,int sU,const FermionField &in, FermionField &out){
-			     int sF,int sU,const FermionField &in, FermionField &out){
+       DiracOptDhopSite(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3
-	DiracOptDhopSite(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3
+     }
-      }
 
-      static void 
+     void DiracOptHandDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U,
-	DiracOptHandDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U,
+				  std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-				std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
+				  int sF,int sU,const FermionField &in, FermionField &out){
-				int sF,int sU,const FermionField &in, FermionField &out){
+       DiracOptDhopSiteDag(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3
-	DiracOptDhopSiteDag(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3
+     }
-      }
+
+     WilsonKernels(const ImplParams &p= ImplParams()) : Base(p) {};
 
     };
 

lib/qcd/action/pseudofermion/TwoFlavour.h

@@ -97,9 +97,10 @@ namespace Grid{
     ////////////////////////////////////////////////////////////////////////
     template<class Impl>
     class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
+    public:
+      INHERIT_IMPL_TYPES(Impl);
 
     private:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
       
       FermionOperator<Impl> & FermOp;// the basic operator
 

lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h

@@ -4,14 +4,21 @@
 namespace Grid{
   namespace QCD{
 
+    // Base even odd HMC on the normal Mee based schur decomposition.
+    //
+    //     M = (Mee Meo) =  (1             0 )   (Mee   0               )  (1 Mee^{-1} Meo)
+    //         (Moe Moo)    (Moe Mee^-1    1 )   (0   Moo-Moe Mee^-1 Meo)  (0   1         )
+    //
+    // Determinant is det of middle factor
+    // This assumes Mee is indept of U.
+    //
     template<class Impl>
     class SchurDifferentiableOperator :  public SchurDiagMooeeOperator<FermionOperator<Impl>,typename Impl::FermionField> 
       {
       public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+      INHERIT_IMPL_TYPES(Impl);
-	
+
-    public:
+ 	typedef FermionOperator<Impl> Matrix;
-	typedef FermionOperator<Impl> Matrix;
 
 	SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator<Matrix,FermionField>(Mat) {};
 
@@ -19,10 +26,11 @@ namespace Grid{
 	
 	  GridBase *fgrid   = this->_Mat.FermionGrid();
 	  GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid();
-	  GridBase *ugrid = this->_Mat.GaugeGrid();
+	  GridBase *ugrid   = this->_Mat.GaugeGrid();
 	  GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid();
 
 	  Real coeff = 1.0;
+
 	  FermionField tmp1(fcbgrid);
 	  FermionField tmp2(fcbgrid);
 
@@ -31,7 +39,7 @@ namespace Grid{
 
 	  // Assert the checkerboard?? or code for either
 	  assert(U.checkerboard==Odd);
-	  assert(V.checkerboard==V.checkerboard);
+	  assert(V.checkerboard==U.checkerboard);
 
 	  GaugeField ForceO(ucbgrid);
 	  GaugeField ForceE(ucbgrid);
@@ -47,6 +55,9 @@ namespace Grid{
 	  this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even 
 	  this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerNo);
 	  
+	  assert(ForceE.checkerboard==Even);
+	  assert(ForceO.checkerboard==Odd);
+
 	  setCheckerboard(Force,ForceE); 
 	  setCheckerboard(Force,ForceO);
 	  Force=-Force;
@@ -61,6 +72,7 @@ namespace Grid{
 	  GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid();
 
 	  Real coeff = 1.0;
+
 	  FermionField tmp1(fcbgrid);
 	  FermionField tmp2(fcbgrid);
 
@@ -85,6 +97,9 @@ namespace Grid{
 	  this->_Mat.MooeeInv(tmp1,tmp2); // even->even 
 	  this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerYes);
 
+	  assert(ForceE.checkerboard==Even);
+	  assert(ForceO.checkerboard==Odd);
+
 	  setCheckerboard(Force,ForceE); 
 	  setCheckerboard(Force,ForceO);
 	  Force=-Force;
@@ -100,7 +115,7 @@ namespace Grid{
     class TwoFlavourEvenOddPseudoFermionAction : public Action<typename Impl::GaugeField> {
 
     public:
-#include <qcd/action/fermion/FermionImplTypedefs.h>
+      INHERIT_IMPL_TYPES(Impl);
 
     private:
       
@@ -180,8 +195,8 @@ namespace Grid{
 
 	// The EE factorised block; normally can replace with zero if det is constant (gauge field indept)
 	// Only really clover term that creates this.
-	FermOp.MooeeInvDag(PhiEven,Y);
+	//	FermOp.MooeeInvDag(PhiEven,Y);
-	action = action + norm2(Y);
+	//	action = action + norm2(Y);
 
 	std::cout << GridLogMessage << "Pseudofermion EO action "<<action<<std::endl;
 	return action;
@@ -217,10 +232,10 @@ namespace Grid{
 
 	// Treat the EE case. (MdagM)^-1 = Minv Minvdag
 	// Deriv defaults to zero.
-        FermOp.MooeeInvDag(PhiOdd,Y);
+	//        FermOp.MooeeInvDag(PhiOdd,Y);
-        FermOp.MooeeInv(Y,X);
+	//      FermOp.MooeeInv(Y,X);
-  	FermOp.MeeDeriv(tmp , Y, X,DaggerNo );    dSdU=tmp;
+	//	FermOp.MeeDeriv(tmp , Y, X,DaggerNo );    dSdU=tmp;
-	FermOp.MeeDeriv(tmp , X, Y,DaggerYes);  dSdU=dSdU+tmp;
+	//  FermOp.MeeDeriv(tmp , X, Y,DaggerYes);  dSdU=dSdU+tmp;
 	
 	dSdU = Ta(dSdU);
 

lib/qcd/hmc/integrators/Integrator_base.h

@@ -20,7 +20,7 @@ namespace Grid{
   namespace QCD{
 
     typedef Action<LatticeLorentzColourMatrix>*  ActPtr; // now force the same size as the rest of the code
-    typedef std::vector<ActPtr> ActionLevel;
+    typedef std::vector<ActPtr>      ActionLevel;
     typedef std::vector<ActionLevel> ActionSet;
     typedef std::vector<Observer*> ObserverList;
     
@@ -45,10 +45,6 @@ namespace Grid{
       void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&);
     }
 
-    
-    
-
-
     template< class IntegratorPolicy >
     class Integrator{
     private:

scripts/hmc.sh

@@ -0,0 +1,6 @@
+#!/bin/bash
+
+LOG=$1
+SWEEPS=$2
+grep Plaq $LOG | tail -n $SWEEPS | awk '{ S=S+$10} END { print S/NR} ' 
+grep dH log_cb | tail -n $SWEEPS | awk '{ S=S+exp(-$10)} END { print S/NR} '

scripts/substitute

@@ -10,7 +10,7 @@ shift
 while (( "$#" )); do
 
 echo $1
-sed -e "s/${WAS}/${IS}/g" -i .bak $1
+sed -e "s:${WAS}:${IS}:g" -i .bak $1
 
 shift
 

tests/Test_hmc_EOWilsonFermionGauge.cc

@@ -35,16 +35,14 @@ int main (int argc, char ** argv)
   //Collect actions
   ActionLevel Level1;
   Level1.push_back(&WilsonNf2);
-  ActionLevel Level2(3);
+  Level1.push_back(&Waction);
-  Level2.push_back(&Waction);
   ActionSet FullSet;
   FullSet.push_back(Level1);
-  FullSet.push_back(Level2);
 
   // Create integrator
-  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
+  //  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
-  //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
+  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
-  IntegratorParameters MDpar(12,20,1.0);
+  IntegratorParameters MDpar(12,40,1.0);
   std::vector<int> rel ={1};
   Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 

tests/Test_hmc_WilsonFermionGauge.cc

@@ -35,17 +35,15 @@ int main (int argc, char ** argv)
   //Collect actions
   ActionLevel Level1;
   Level1.push_back(&WilsonNf2);
-  ActionLevel Level2(3);
+  Level1.push_back(&Waction);
-  Level2.push_back(&Waction);
   ActionSet FullSet;
   FullSet.push_back(Level1);
-  FullSet.push_back(Level2);
+  //  FullSet.push_back(Level2);
 
   // Create integrator
   typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
   //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
   IntegratorParameters MDpar(12,20,1.0);
-  std::vector<int> rel ={1};
   Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC