Staggered

lib/qcd/action/ActionParams.h

@@ -45,6 +45,10 @@ namespace QCD {
       WilsonImplParams() : overlapCommsCompute(false) {};
     };
 
+    struct StaggeredImplParams {
+      StaggeredImplParams()  {};
+    };
+
     struct OneFlavourRationalParams { 
       RealD  lo;
       RealD  hi;

lib/qcd/action/Actions.h

@@ -53,6 +53,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/qcd/action/fermion/FermionOperatorImpl.h>
 #include <Grid/qcd/action/fermion/FermionOperator.h>
 #include <Grid/qcd/action/fermion/WilsonKernels.h>        //used by all wilson type fermions
+#include <Grid/qcd/action/fermion/StaggeredKernels.h>        //used by all wilson type fermions
 
 ////////////////////////////////////////////
 // Gauge Actions
@@ -108,6 +109,10 @@ typedef SymanzikGaugeAction<ConjugateGimplD>        ConjugateSymanzikGaugeAction
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 
 
+#define FermOpStaggeredTemplateInstantiate(A) \
+  template class A<StaggeredImplF>;		\
+  template class A<StaggeredImplD>;		
+
 #define FermOp4dVecTemplateInstantiate(A) \
   template class A<WilsonImplF>;		\
   template class A<WilsonImplD>;		\
@@ -147,6 +152,8 @@ typedef SymanzikGaugeAction<ConjugateGimplD>        ConjugateSymanzikGaugeAction
 
 //#include <Grid/qcd/action/fermion/CloverFermion.h>
 
+#include <Grid/qcd/action/fermion/ImprovedStaggeredFermion.h>
+
 #include <Grid/qcd/action/fermion/CayleyFermion5D.h>     // Cayley types
 #include <Grid/qcd/action/fermion/DomainWallFermion.h>
 #include <Grid/qcd/action/fermion/DomainWallFermion.h>

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -499,6 +499,97 @@ PARALLEL_FOR_LOOP
 
 };
 
+
+  /////////////////////////////////////////////////////////////////////////////
+  // Single flavour one component spinors with colour index
+  /////////////////////////////////////////////////////////////////////////////
+  template <class S, class Representation = FundamentalRepresentation >
+  class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation::Dimension > > {
+
+    public:
+
+    typedef RealD  _Coeff_t ;
+    static const int Dimension = Representation::Dimension;
+    typedef PeriodicGaugeImpl<GaugeImplTypes<S, Dimension > > Gimpl;
+      
+    //Necessary?
+    constexpr bool is_fundamental() const{return Dimension == Nc ? 1 : 0;}
+    
+    const bool LsVectorised=false;
+    typedef _Coeff_t Coeff_t;
+
+    INHERIT_GIMPL_TYPES(Gimpl);
+      
+    template <typename vtype> using iImplSpinor            = iScalar<iScalar<iVector<vtype, Dimension> > >;
+    template <typename vtype> using iImplHalfSpinor        = iVector<iScalar<iVector<vtype, Dimension> >, Ngp>;
+    template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>;
+    
+    typedef iImplSpinor<Simd>            SiteSpinor;
+    typedef iImplHalfSpinor<Simd>        SiteHalfSpinor;
+    typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
+    
+    typedef Lattice<SiteSpinor>            FermionField;
+    typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
+    
+    typedef SimpleCompressor<SiteSpinor> Compressor;
+    typedef StaggeredImplParams ImplParams;
+    typedef CartesianStencil<SiteSpinor, SiteSpinor> StencilImpl;
+    
+    ImplParams Params;
+    
+    StaggeredImpl(const ImplParams &p = ImplParams()) : Params(p){};
+      
+    inline void multLink(SiteSpinor &phi,
+			 const SiteDoubledGaugeField &U,
+			 const SiteSpinor &chi,
+			 int mu){
+      mult(&phi(), &U(mu), &chi());
+    }
+    inline void multLinkAdd(SiteSpinor &phi,
+			    const SiteDoubledGaugeField &U,
+			    const SiteSpinor &chi,
+			    int mu){
+      mac(&phi(), &U(mu), &chi());
+    }
+      
+    template <class ref>
+    inline void loadLinkElement(Simd &reg, ref &memory) {
+      reg = memory;
+    }
+      
+    inline void DoubleStore(GridBase *GaugeGrid,
+			    DoubledGaugeField &Uds,
+			    DoubledGaugeField &UUUds, // for Naik term
+			    const GaugeField &Umu) {
+      conformable(Uds._grid, GaugeGrid);
+      conformable(Umu._grid, GaugeGrid);
+      GaugeLinkField U(GaugeGrid);
+      for (int mu = 0; mu < Nd; mu++) {
+	U = PeekIndex<LorentzIndex>(Umu, mu);
+	PokeIndex<LorentzIndex>(Uds, U, mu);
+	PokeIndex<LorentzIndex>(UUUds, U, mu);
+	std::cout << GridLogMessage << " NOT created the treble links for staggered yet" <<std::endl;
+	std::cout << GridLogMessage << " Must do this and also apply the staggered phases which requires understanding the action conventions" <<std::endl;
+	U = adj(Cshift(U, mu, -1));
+	PokeIndex<LorentzIndex>(Uds, U, mu + 4);
+	PokeIndex<LorentzIndex>(UUUds, U, mu+4);
+      }
+    }
+
+    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){
+      assert (0); 
+      // Must never hit
+    }
+  };
+
+
+
  typedef WilsonImpl<vComplex,  FundamentalRepresentation > WilsonImplR;   // Real.. whichever prec
  typedef WilsonImpl<vComplexF, FundamentalRepresentation > WilsonImplF;  // Float
  typedef WilsonImpl<vComplexD, FundamentalRepresentation > WilsonImplD;  // Double
@@ -527,6 +618,10 @@ PARALLEL_FOR_LOOP
  typedef GparityWilsonImpl<vComplexF, Nc> GparityWilsonImplF;  // Float
  typedef GparityWilsonImpl<vComplexD, Nc> GparityWilsonImplD;  // Double
 
+ typedef StaggeredImpl<vComplex,  FundamentalRepresentation > StaggeredImplR;   // Real.. whichever prec
+ typedef StaggeredImpl<vComplexF, FundamentalRepresentation > StaggeredImplF;  // Float
+ typedef StaggeredImpl<vComplexD, FundamentalRepresentation > StaggeredImplD;  // Double
+
 }}
 
 #endif

lib/qcd/action/fermion/ImprovedStaggeredFermion.cc

@@ -0,0 +1,309 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/ImprovedStaggeredFermion.cc
+
+Copyright (C) 2015
+
+Author: Azusa Yamaguchi, Peter Boyle
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid.h>
+
+namespace Grid {
+namespace QCD {
+
+const std::vector<int> 
+ImprovedStaggeredFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3});
+const std::vector<int> 
+ImprovedStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3, 3, 3, -3, -3, -3, -3});
+
+/////////////////////////////////
+// Constructor and gauge import
+/////////////////////////////////
+
+template <class Impl>
+ImprovedStaggeredFermion<Impl>::ImprovedStaggeredFermion(GaugeField &_Umu, GridCartesian &Fgrid,
+							 GridRedBlackCartesian &Hgrid, RealD _mass,
+							 const ImplParams &p)
+    : Kernels(p),
+      _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),
+      Lebesgue(_grid),
+      LebesgueEvenOdd(_cbgrid),
+      Umu(&Fgrid),
+      UmuEven(&Hgrid),
+      UmuOdd(&Hgrid),
+      UUUmu(&Fgrid),
+      UUUmuEven(&Hgrid),
+      UUUmuOdd(&Hgrid) {
+  // Allocate the required comms buffer
+  ImportGauge(_Umu);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::ImportGauge(const GaugeField &_Umu) {
+  GaugeField HUmu(_Umu._grid);
+  HUmu = _Umu * (-0.5);
+  Impl::DoubleStore(GaugeGrid(), Umu, UUUmu, HUmu);
+  pickCheckerboard(Even, UmuEven, Umu);
+  pickCheckerboard(Odd, UmuOdd, Umu);
+  pickCheckerboard(Even, UUUmuEven, UUUmu);
+  pickCheckerboard(Odd, UUUmuOdd, UUUmu);
+}
+
+/////////////////////////////
+// Implement the interface
+/////////////////////////////
+
+template <class Impl>
+RealD ImprovedStaggeredFermion<Impl>::M(const FermionField &in, FermionField &out) {
+  out.checkerboard = in.checkerboard;
+  Dhop(in, out, DaggerNo);
+  return axpy_norm(out, mass, in, out);
+}
+
+template <class Impl>
+RealD ImprovedStaggeredFermion<Impl>::Mdag(const FermionField &in, FermionField &out) {
+  out.checkerboard = in.checkerboard;
+  Dhop(in, out, DaggerYes);
+  return axpy_norm(out, mass, in, out);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::Meooe(const FermionField &in, FermionField &out) {
+  if (in.checkerboard == Odd) {
+    DhopEO(in, out, DaggerNo);
+  } else {
+    DhopOE(in, out, DaggerNo);
+  }
+}
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
+  if (in.checkerboard == Odd) {
+    DhopEO(in, out, DaggerYes);
+  } else {
+    DhopOE(in, out, DaggerYes);
+  }
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::Mooee(const FermionField &in, FermionField &out) {
+  out.checkerboard = in.checkerboard;
+  typename FermionField::scalar_type scal(mass);
+  out = scal * in;
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) {
+  out.checkerboard = in.checkerboard;
+  Mooee(in, out);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) {
+  out.checkerboard = in.checkerboard;
+  out = (1.0 / (mass)) * in;
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::MooeeInvDag(const FermionField &in,
+                                      FermionField &out) {
+  out.checkerboard = in.checkerboard;
+  MooeeInv(in, out);
+}
+
+///////////////////////////////////
+// Internal
+///////////////////////////////////
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U, DoubledGaugeField &UUU, 
+						   GaugeField & mat,
+						   const FermionField &A, const FermionField &B, int dag) {
+  assert((dag == DaggerNo) || (dag == DaggerYes));
+
+  Compressor compressor;
+
+  FermionField Btilde(B._grid);
+  FermionField Atilde(B._grid);
+  Atilde = A;
+
+  st.HaloExchange(B, compressor);
+
+  for (int mu = 0; mu < Nd; mu++) {
+
+    ////////////////////////
+    // Call the single hop
+    ////////////////////////
+    PARALLEL_FOR_LOOP
+    for (int sss = 0; sss < B._grid->oSites(); sss++) {
+      Kernels::DhopDir(st, U, UUU, st.CommBuf(), sss, sss, B, Btilde, mu,1);
+    }
+
+    // Force in three link terms
+    //
+    //    Impl::InsertForce4D(mat, Btilde, Atilde, mu);
+    //
+    // dU_ac(x)/dt = i p_ab U_bc(x)
+    //
+    // => dS_f/dt = dS_f/dU_ac(x) . dU_ac(x)/dt =  i p_ab U_bc(x) dS_f/dU_ac(x) 
+    //
+    // One link: form fragments S_f = A U B 
+    //
+    //         write Btilde = U(x) B(x+mu)
+    //
+    // mat+= TraceIndex<SpinIndex>(outerProduct(Btilde,A)); 
+    // 
+    // Three link: form fragments S_f = A UUU B 
+    //
+    // mat+= outer ( A, UUUB) <-- Best take DhopDeriv with one linke or identity matrix
+    // mat+= outer ( AU, UUB) <-- and then use covariant cshift?
+    // mat+= outer ( AUU, UB) <-- Returned from call to DhopDir
+
+    assert(0);// need to figure out the force interface with a blasted three link term.
+    
+  }
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
+
+  conformable(U._grid, _grid);
+  conformable(U._grid, V._grid);
+  conformable(U._grid, mat._grid);
+
+  mat.checkerboard = U.checkerboard;
+
+  DerivInternal(Stencil, Umu, UUUmu, mat, U, V, dag);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
+
+  conformable(U._grid, _cbgrid);
+  conformable(U._grid, V._grid);
+  conformable(U._grid, mat._grid);
+
+  assert(V.checkerboard == Even);
+  assert(U.checkerboard == Odd);
+  mat.checkerboard = Odd;
+
+  DerivInternal(StencilEven, UmuOdd, UUUmuOdd, mat, U, V, dag);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
+
+  conformable(U._grid, _cbgrid);
+  conformable(U._grid, V._grid);
+  conformable(U._grid, mat._grid);
+
+  assert(V.checkerboard == Odd);
+  assert(U.checkerboard == Even);
+  mat.checkerboard = Even;
+
+  DerivInternal(StencilOdd, UmuEven, UUUmuEven, mat, U, V, dag);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::Dhop(const FermionField &in, FermionField &out, int dag) {
+  conformable(in._grid, _grid);  // verifies full grid
+  conformable(in._grid, out._grid);
+
+  out.checkerboard = in.checkerboard;
+
+  DhopInternal(Stencil, Lebesgue, Umu, UUUmu, in, out, dag);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::DhopOE(const FermionField &in, FermionField &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, LebesgueEvenOdd, UmuOdd, UUUmuOdd, in, out, dag);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::DhopEO(const FermionField &in, FermionField &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, LebesgueEvenOdd, UmuEven, UUUmuEven, in, out, dag);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::Mdir(const FermionField &in, FermionField &out, int dir, int disp) {
+  DhopDir(in, out, dir, disp);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::DhopDir(const FermionField &in, FermionField &out, int dir, int disp) {
+
+  Compressor compressor;
+  Stencil.HaloExchange(in, compressor);
+
+PARALLEL_FOR_LOOP
+  for (int sss = 0; sss < in._grid->oSites(); sss++) {
+    Kernels::DhopDir(Stencil, Umu, UUUmu, Stencil.CommBuf(), sss, sss, in, out, dir, disp);
+  }
+};
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder &lo,
+						  DoubledGaugeField &U,
+						  DoubledGaugeField &UUU,
+						  const FermionField &in,
+						  FermionField &out, int dag) {
+  assert((dag == DaggerNo) || (dag == DaggerYes));
+
+  Compressor compressor;
+  st.HaloExchange(in, compressor);
+
+  if (dag == DaggerYes) {
+    PARALLEL_FOR_LOOP
+    for (int sss = 0; sss < in._grid->oSites(); sss++) {
+      Kernels::DhopSiteDag(st, lo, U, UUU, st.CommBuf(), sss, sss, in, out);
+    }
+  } else {
+    PARALLEL_FOR_LOOP
+    for (int sss = 0; sss < in._grid->oSites(); sss++) {
+      Kernels::DhopSite(st, lo, U, UUU, st.CommBuf(), sss, sss, in, out);
+    }
+  }
+};
+
+FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion);
+
+  //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion);
+  //TwoIndexFermOpTemplateInstantiate(ImprovedStaggeredFermion);
+
+}}

lib/qcd/action/fermion/ImprovedStaggeredFermion.h

@@ -0,0 +1,152 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/ImprovedStaggered.h
+
+Copyright (C) 2015
+
+Author: Azusa Yamaguchi, Peter Boyle
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef GRID_QCD_IMPR_STAG_FERMION_H
+#define GRID_QCD_IMPR_STAG_FERMION_H
+
+namespace Grid {
+
+namespace QCD {
+
+class ImprovedStaggeredFermionStatic {
+ public:
+  static const std::vector<int> directions;
+  static const std::vector<int> displacements;
+  static const int npoint = 16;
+};
+
+template <class Impl>
+class ImprovedStaggeredFermion : public StaggeredKernels<Impl>, public ImprovedStaggeredFermionStatic {
+ public:
+  INHERIT_IMPL_TYPES(Impl);
+  typedef StaggeredKernels<Impl> Kernels;
+
+  ///////////////////////////////////////////////////////////////
+  // 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; cut number routines if pass dagger argument
+  // and also make interface more uniformly consistent
+  //////////////////////////////////////////////////////////////////
+  RealD M(const FermionField &in, FermionField &out);
+  RealD Mdag(const FermionField &in, FermionField &out);
+
+  /////////////////////////////////////////////////////////
+  // half checkerboard operations
+  /////////////////////////////////////////////////////////
+  void Meooe(const FermionField &in, FermionField &out);
+  void MeooeDag(const FermionField &in, FermionField &out);
+
+  // allow override for twisted mass and clover
+  virtual void Mooee(const FermionField &in, FermionField &out);
+  virtual void MooeeDag(const FermionField &in, FermionField &out);
+  virtual void MooeeInv(const FermionField &in, FermionField &out);
+  virtual void MooeeInvDag(const FermionField &in, FermionField &out);
+
+  ////////////////////////
+  // Derivative interface
+  ////////////////////////
+  // Interface calls an internal routine
+  void DhopDeriv  (GaugeField &mat, const FermionField &U, const FermionField &V, int dag);
+  void DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag);
+  void DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag);
+
+  ///////////////////////////////////////////////////////////////
+  // non-hermitian hopping term; half cb or both
+  ///////////////////////////////////////////////////////////////
+  void Dhop  (const FermionField &in, FermionField &out, int dag);
+  void DhopOE(const FermionField &in, FermionField &out, int dag);
+  void DhopEO(const FermionField &in, FermionField &out, int dag);
+
+  ///////////////////////////////////////////////////////////////
+  // Multigrid assistance; force term uses too
+  ///////////////////////////////////////////////////////////////
+  void Mdir(const FermionField &in, FermionField &out, int dir, int disp);
+  void DhopDir(const FermionField &in, FermionField &out, int dir, int disp);
+
+  ///////////////////////////////////////////////////////////////
+  // Extra methods added by derived
+  ///////////////////////////////////////////////////////////////
+  void DerivInternal(StencilImpl &st, 
+		     DoubledGaugeField &U,DoubledGaugeField &UUU,
+		     GaugeField &mat, 
+		     const FermionField &A, const FermionField &B, int dag);
+
+  void DhopInternal(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,DoubledGaugeField &UUU,
+                    const FermionField &in, FermionField &out, int dag);
+
+  // Constructor
+  ImprovedStaggeredFermion(GaugeField &_Umu, GridCartesian &Fgrid,
+			   GridRedBlackCartesian &Hgrid, RealD _mass,
+			   const ImplParams &p = ImplParams());
+
+  // DoubleStore impl dependent
+  void ImportGauge(const GaugeField &_Umu);
+
+  ///////////////////////////////////////////////////////////////
+  // Data members require to support the functionality
+  ///////////////////////////////////////////////////////////////
+
+  //    protected:
+ public:
+  // any other parameters of action ???
+
+  RealD mass;
+
+  GridBase *_grid;
+  GridBase *_cbgrid;
+
+  // Defines the stencils for even and odd
+  StencilImpl Stencil;
+  StencilImpl StencilEven;
+  StencilImpl StencilOdd;
+
+  // Copy of the gauge field , with even and odd subsets
+  DoubledGaugeField Umu;
+  DoubledGaugeField UmuEven;
+  DoubledGaugeField UmuOdd;
+
+  DoubledGaugeField UUUmu;
+  DoubledGaugeField UUUmuEven;
+  DoubledGaugeField UUUmuOdd;
+
+  LebesgueOrder Lebesgue;
+  LebesgueOrder LebesgueEvenOdd;
+};
+
+typedef ImprovedStaggeredFermion<StaggeredImplF> ImprovedStaggeredFermionF;
+typedef ImprovedStaggeredFermion<StaggeredImplD> ImprovedStaggeredFermionD;
+
+}
+}
+#endif

lib/qcd/action/fermion/StaggeredKernels.cc

@@ -0,0 +1,223 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
+
+Copyright (C) 2015
+
+Author: Azusa Yamaguchi, Peter Boyle
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid.h>
+namespace Grid {
+namespace QCD {
+
+template <class Impl>
+StaggeredKernels<Impl>::StaggeredKernels(const ImplParams &p) : Base(p){};
+
+////////////////////////////////////////////
+// Generic implementation; move to different file?
+////////////////////////////////////////////
+
+template <class Impl>
+void StaggeredKernels<Impl>::DhopSiteDepth(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+					   SiteSpinor *buf, int sF,
+					   int sU, const FermionField &in, SiteSpinor &out,int threeLink) {
+  const SiteSpinor *chi_p;
+  SiteSpinor chi;
+  SiteSpinor Uchi;
+  StencilEntry *SE;
+  int ptype;
+  int skew = 0;
+  if (threeLink) skew=8;
+  ///////////////////////////
+  // Xp
+  ///////////////////////////
+
+  SE = st.GetEntry(ptype, Xp+skew, sF);
+  if (SE->_is_local) {
+    if (SE->_permute) {
+      chi_p = &chi;
+      permute(chi,  in._odata[SE->_offset], ptype);
+    } else {
+      chi_p = &in._odata[SE->_offset];
+    }
+  } else {
+    chi_p = &buf[SE->_offset];
+  }
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Xp);
+
+  ///////////////////////////
+  // Yp
+  ///////////////////////////
+  SE = st.GetEntry(ptype, Yp+skew, sF);
+  if (SE->_is_local) {
+    if (SE->_permute) {
+      chi_p = &chi;
+      permute(chi,  in._odata[SE->_offset], ptype);
+    } else {
+      chi_p = &in._odata[SE->_offset];
+    }
+  } else {
+    chi_p = &buf[SE->_offset];
+  }
+  Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Yp);
+
+  ///////////////////////////
+  // Zp
+  ///////////////////////////
+  SE = st.GetEntry(ptype, Zp+skew, sF);
+  if (SE->_is_local) {
+    if (SE->_permute) {
+      chi_p = &chi;
+      permute(chi,  in._odata[SE->_offset], ptype);
+    } else {
+      chi_p = &in._odata[SE->_offset];
+    }
+  } else {
+    chi_p = &buf[SE->_offset];
+  }
+  Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Zp);
+
+  ///////////////////////////
+  // Tp
+  ///////////////////////////
+  SE = st.GetEntry(ptype, Tp+skew, sF);
+  if (SE->_is_local) {
+    if (SE->_permute) {
+      chi_p = &chi;
+      permute(chi,  in._odata[SE->_offset], ptype);
+    } else {
+      chi_p = &in._odata[SE->_offset];
+    }
+  } else {
+    chi_p = &buf[SE->_offset];
+  }
+  Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Tp);
+
+  ///////////////////////////
+  // Xm
+  ///////////////////////////
+  SE = st.GetEntry(ptype, Xm+skew, sF);
+  if (SE->_is_local) {
+    if (SE->_permute) {
+      chi_p = &chi;
+      permute(chi,  in._odata[SE->_offset], ptype);
+    } else {
+      chi_p = &in._odata[SE->_offset];
+    }
+  } else {
+    chi_p = &buf[SE->_offset];
+  }
+  Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Xm);
+
+  ///////////////////////////
+  // Ym
+  ///////////////////////////
+  SE = st.GetEntry(ptype, Ym+skew, sF);
+  if (SE->_is_local) {
+    if (SE->_permute) {
+      chi_p = &chi;
+      permute(chi,  in._odata[SE->_offset], ptype);
+    } else {
+      chi_p = &in._odata[SE->_offset];
+    }
+  } else {
+    chi_p = &buf[SE->_offset];
+  }
+  Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Ym);
+
+  ///////////////////////////
+  // Zm
+  ///////////////////////////
+  SE = st.GetEntry(ptype, Zm+skew, sF);
+  if (SE->_is_local) {
+    if (SE->_permute) {
+      chi_p = &chi;
+      permute(chi,  in._odata[SE->_offset], ptype);
+    } else {
+      chi_p = &in._odata[SE->_offset];
+    }
+  } else {
+    chi_p = &buf[SE->_offset];
+  }
+  Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Zm);
+
+  ///////////////////////////
+  // Tm
+  ///////////////////////////
+  SE = st.GetEntry(ptype, Tm+skew, sF);
+  if (SE->_is_local) {
+    if (SE->_permute) {
+      chi_p = &chi;
+      permute(chi,  in._odata[SE->_offset], ptype);
+    } else {
+      chi_p = &in._odata[SE->_offset];
+    }
+  } else {
+    chi_p = &buf[SE->_offset];
+  }
+  Impl::multLinkAdd(Uchi, U._odata[sU], *chi_p, Tm);
+
+  vstream(out, Uchi);
+};
+
+// Need controls to do interior, exterior, or both
+template <class Impl>
+void StaggeredKernels<Impl>::DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, DoubledGaugeField &UUU,
+						  SiteSpinor *buf, int sF,
+						  int sU, const FermionField &in, FermionField &out) {
+  SiteSpinor naik;
+  SiteSpinor naive;
+  int oneLink  =0;
+  int threeLink=1;
+  DhopSiteDepth(st,lo,U,buf,sF,sU,in,naive,oneLink);
+  DhopSiteDepth(st,lo,UUU,buf,sF,sU,in,naik,threeLink);
+  out._odata[sF] =naive+naik;
+};
+template <class Impl>
+void StaggeredKernels<Impl>::DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, DoubledGaugeField &UUU,
+				      SiteSpinor *buf, int sF,
+				      int sU, const FermionField &in, FermionField &out) {
+  SiteSpinor naik;
+  SiteSpinor naive;
+  int oneLink  =0;
+  int threeLink=1;
+  DhopSiteDepth(st,lo,U,buf,sF,sU,in,naive,oneLink);
+  DhopSiteDepth(st,lo,UUU,buf,sF,sU,in,naik,threeLink);
+  out._odata[sF] =-naive-naik;
+};
+
+template <class Impl>
+void StaggeredKernels<Impl>::DhopDir( StencilImpl &st, DoubledGaugeField &U,  DoubledGaugeField &UUU, SiteSpinor *buf, int sF,
+				      int sU, const FermionField &in, FermionField &out, int dir, int disp) 
+{
+  // Disp should be either +1,-1,+3,-3
+  // What about "dag" ?
+  // Because we work out pU . dS/dU 
+  // U
+  assert(0);
+}
+
+FermOpStaggeredTemplateInstantiate(StaggeredKernels);
+
+}}
+

lib/qcd/action/fermion/StaggeredKernels.h

@@ -0,0 +1,70 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/StaggeredKernels.h
+
+Copyright (C) 2015
+
+Author: Azusa Yamaguchi, Peter Boyle
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef GRID_QCD_STAGGERED_KERNELS_H
+#define GRID_QCD_STAGGERED_KERNELS_H
+
+namespace Grid {
+namespace QCD {
+
+  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+  // Helper routines that implement Staggered stencil for a single site.
+  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+class StaggeredKernelsStatic { 
+ public:
+};
+ 
+template<class Impl> class StaggeredKernels : public FermionOperator<Impl> , public StaggeredKernelsStatic { 
+ public:
+   
+  INHERIT_IMPL_TYPES(Impl);
+  typedef FermionOperator<Impl> Base;
+   
+public:
+    
+   void DhopDir(StencilImpl &st, DoubledGaugeField &U, DoubledGaugeField &UUU, SiteSpinor * buf,
+		      int sF, int sU, const FermionField &in, FermionField &out, int dir,int disp);
+
+   void DhopSiteDepth(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteSpinor * buf,
+		     int sF, int sU, const FermionField &in, SiteSpinor &out,int threeLink);
+      
+   void DhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, DoubledGaugeField &UUU, SiteSpinor * buf,
+		int sF, int sU, const FermionField &in, FermionField &out);
+
+   void DhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, DoubledGaugeField &UUU, SiteSpinor * buf,
+		   int sF, int sU, const FermionField &in, FermionField &out);
+  
+public:
+
+  StaggeredKernels(const ImplParams &p = ImplParams());
+
+};
+    
+}}
+
+#endif

lib/qcd/action/fermion/WilsonFermion.cc

@@ -34,10 +34,9 @@ directory
 namespace Grid {
 namespace QCD {
 
-const std::vector<int> WilsonFermionStatic::directions({0, 1, 2, 3, 0, 1, 2,
-                                                        3});
-const std::vector<int> WilsonFermionStatic::displacements({1, 1, 1, 1, -1, -1,
-                                                           -1, -1});
+const std::vector<int> WilsonFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3});
+const std::vector<int> WilsonFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1});
+
 int WilsonFermionStatic::HandOptDslash;
 
 /////////////////////////////////
@@ -166,8 +165,7 @@ void WilsonFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U,
     ////////////////////////
     PARALLEL_FOR_LOOP
     for (int sss = 0; sss < B._grid->oSites(); sss++) {
-      Kernels::DiracOptDhopDir(st, U, st.CommBuf(), sss, sss, B, Btilde, mu,
-                               gamma);
+      Kernels::DiracOptDhopDir(st, U, st.CommBuf(), sss, sss, B, Btilde, mu, gamma);
     }
 
     //////////////////////////////////////////////////
@@ -277,8 +275,7 @@ void WilsonFermion<Impl>::DhopDirDisp(const FermionField &in, FermionField &out,
 
   PARALLEL_FOR_LOOP
   for (int sss = 0; sss < in._grid->oSites(); sss++) {
-    Kernels::DiracOptDhopDir(Stencil, Umu, Stencil.CommBuf(), sss, sss, in, out,
-                             dirdisp, gamma);
+    Kernels::DiracOptDhopDir(Stencil, Umu, Stencil.CommBuf(), sss, sss, in, out, dirdisp, gamma);
   }
 };
 
@@ -295,14 +292,12 @@ void WilsonFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder &lo,
   if (dag == DaggerYes) {
     PARALLEL_FOR_LOOP
     for (int sss = 0; sss < in._grid->oSites(); sss++) {
-      Kernels::DiracOptDhopSiteDag(st, lo, U, st.CommBuf(), sss, sss, 1, 1, in,
-                                   out);
+      Kernels::DiracOptDhopSiteDag(st, lo, U, st.CommBuf(), sss, sss, 1, 1, in, out);
     }
   } else {
     PARALLEL_FOR_LOOP
     for (int sss = 0; sss < in._grid->oSites(); sss++) {
-      Kernels::DiracOptDhopSite(st, lo, U, st.CommBuf(), sss, sss, 1, 1, in,
-                                out);
+      Kernels::DiracOptDhopSite(st, lo, U, st.CommBuf(), sss, sss, 1, 1, in, out);
     }
   }
 };

lib/qcd/action/pseudofermion/TwoFlavour.h

@@ -63,8 +63,7 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
         Phi(Op.FermionGrid()){};
 
   //////////////////////////////////////////////////////////////////////////////////////
-  // Push the gauge field in to the dops. Assume any BC's and smearing already
-  // applied
+  // Push the gauge field in to the dops. Assume any BC's and smearing already applied
   //////////////////////////////////////////////////////////////////////////////////////
   virtual void refresh(const GaugeField &U, GridParallelRNG &pRNG) {
     // P(phi) = e^{- phi^dag (MdagM)^-1 phi}
@@ -107,8 +106,7 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
     MdagMOp.Op(X, Y);
 
     RealD action = norm2(Y);
-    std::cout << GridLogMessage << "Pseudofermion action " << action
-              << std::endl;
+    std::cout << GridLogMessage << "Pseudofermion action " << action << std::endl;
     return action;
   };
 
@@ -119,6 +117,7 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
   //
   //       = - Ydag dM X  - Xdag dMdag Y
   //
+  // 
   //////////////////////////////////////////////////////
   virtual void deriv(const GaugeField &U, GaugeField &dSdU) {
     FermOp.ImportGauge(U);
@@ -133,8 +132,7 @@ class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
     DerivativeSolver(MdagMOp, Phi, X); // X = (MdagM)^-1 phi    
     MdagMOp.Op(X, Y);                  // Y = M X = (Mdag)^-1 phi
 
-    // Our conventions really make this UdSdU; We do not differentiate wrt Udag
-    // here.
+    // Our conventions really make this UdSdU; We do not differentiate wrt Udag here.
     // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt.
 
     FermOp.MDeriv(tmp, Y, X, DaggerNo);