finer timers in Benchmark_IO

2026-10-29 08:53:19 +00:00 · 2021-06-17 11:57:02 +01:00
209 changed files with 1818 additions and 14515 deletions
@@ -88,7 +88,6 @@ Thumbs.db
 # build directory #
 ###################
 build*/*
-Documentation/_build

 # IDE related files #
 #####################
@@ -34,9 +34,6 @@ directory

 #if defined __GNUC__ && __GNUC__>=6
 #pragma GCC diagnostic ignored "-Wignored-attributes"
-#endif
-#if defined __GNUC__ 
-#pragma GCC diagnostic ignored "-Wpsabi"
 #endif

 //disables and intel compiler specific warning (in json.hpp)
@@ -36,7 +36,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/GridCore.h>
 #include <Grid/qcd/QCD.h>
 #include <Grid/qcd/spin/Spin.h>
-#include <Grid/qcd/gparity/Gparity.h>
 #include <Grid/qcd/utils/Utils.h>
 #include <Grid/qcd/representations/Representations.h>
 NAMESPACE_CHECK(GridQCDCore);
@@ -54,7 +54,6 @@ NAMESPACE_CHECK(BiCGSTAB);
 #include <Grid/algorithms/iterative/SchurRedBlack.h>
 #include <Grid/algorithms/iterative/ConjugateGradientMultiShift.h>
 #include <Grid/algorithms/iterative/ConjugateGradientMixedPrec.h>
-#include <Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h>
 #include <Grid/algorithms/iterative/BiCGSTABMixedPrec.h>
 #include <Grid/algorithms/iterative/BlockConjugateGradient.h>
 #include <Grid/algorithms/iterative/ConjugateGradientReliableUpdate.h>
@@ -442,8 +442,6 @@ public:
    for(int p=0; p<geom.npoint; p++)
      points[p] = geom.points_dagger[p];

-    auto points_p = &points[0];
-
    RealD* dag_factor_p = &dag_factor[0];

    accelerator_for(sss, Grid()->oSites()*nbasis, Nsimd, {
@@ -455,7 +453,7 @@ public:
      StencilEntry *SE;

      for(int p=0;p<geom_v.npoint;p++){
-        int point = points_p[p];
+        int point = points[p];

 	SE=Stencil_v.GetEntry(ptype,point,ss);

@@ -710,8 +708,6 @@ public:
    for(int p=0; p<npoint; p++)
      points[p] = (dag && !hermitian) ? geom.points_dagger[p] : p;

-    auto points_p = &points[0];
-
    Vector<Aview> AcceleratorViewContainer;
    for(int p=0;p<npoint;p++) AcceleratorViewContainer.push_back(a[p].View(AcceleratorRead));
    Aview *Aview_p = & AcceleratorViewContainer[0];
@@ -732,7 +728,7 @@ public:
        StencilEntry *SE;

        for(int p=0;p<npoint;p++){
-          int point = points_p[p];
+          int point = points[p];
          SE=st_v.GetEntry(ptype,point,ss);

          if(SE->_is_local) {
@@ -758,7 +754,7 @@ public:
        StencilEntry *SE;

        for(int p=0;p<npoint;p++){
-          int point = points_p[p];
+          int point = points[p];
          SE=st_v.GetEntry(ptype,point,ss);

          if(SE->_is_local) {
@@ -136,7 +136,7 @@ public:
    flops=0;
    usec =0;
    Coordinate layout(Nd,1);
-    sgrid = new GridCartesian(dimensions,layout,processors,*grid);
+    sgrid = new GridCartesian(dimensions,layout,processors);
  };
    
  ~FFT ( void)  {
@@ -182,7 +182,7 @@ public:
    pencil_gd[dim] = G*processors[dim];
      
    // Pencil global vol LxLxGxLxL per node
-    GridCartesian pencil_g(pencil_gd,layout,processors,*vgrid);
+    GridCartesian pencil_g(pencil_gd,layout,processors);
      
    // Construct pencils
    typedef typename vobj::scalar_object sobj;
@@ -530,16 +530,6 @@ public:
 template<class Field> class LinearFunction {
 public:
  virtual void operator() (const Field &in, Field &out) = 0;
-
-  virtual void operator() (const std::vector<Field> &in, std::vector<Field> &out)
-  {
-    assert(in.size() == out.size());
-
-    for (unsigned int i = 0; i < in.size(); ++i)
-    {
-      (*this)(in[i], out[i]);
-    }
-  }
 };

 template<class Field> class IdentityLinearFunction : public LinearFunction<Field> {
@@ -292,7 +292,6 @@ public:
 template<class Field>
 class ChebyshevLanczos : public Chebyshev<Field> {
 private:
-
  std::vector<RealD> Coeffs;
  int order;
  RealD alpha;
@@ -48,7 +48,6 @@ NAMESPACE_BEGIN(Grid);
    Integer TotalInnerIterations; //Number of inner CG iterations
    Integer TotalOuterIterations; //Number of restarts
    Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step
-    RealD TrueResidual;

    //Option to speed up *inner single precision* solves using a LinearFunction that produces a guess
    LinearFunction<FieldF> *guesser;
@@ -68,7 +67,6 @@ NAMESPACE_BEGIN(Grid);
    }
  
  void operator() (const FieldD &src_d_in, FieldD &sol_d){
-    std::cout << GridLogMessage << "MixedPrecisionConjugateGradient: Starting mixed precision CG with outer tolerance " << Tolerance << " and inner tolerance " << InnerTolerance << std::endl;
    TotalInnerIterations = 0;
 	
    GridStopWatch TotalTimer;
@@ -81,11 +79,6 @@ NAMESPACE_BEGIN(Grid);
    RealD stop = src_norm * Tolerance*Tolerance;

    GridBase* DoublePrecGrid = src_d_in.Grid();
-
-    //Generate precision change workspaces
-    precisionChangeWorkspace wk_dp_from_sp(DoublePrecGrid, SinglePrecGrid);
-    precisionChangeWorkspace wk_sp_from_dp(SinglePrecGrid, DoublePrecGrid);
-
    FieldD tmp_d(DoublePrecGrid);
    tmp_d.Checkerboard() = cb;
    
@@ -103,7 +96,6 @@ NAMESPACE_BEGIN(Grid);
    FieldF sol_f(SinglePrecGrid);
    sol_f.Checkerboard() = cb;
    
-    std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Starting initial inner CG with tolerance " << inner_tol << std::endl;
    ConjugateGradient<FieldF> CG_f(inner_tol, MaxInnerIterations);
    CG_f.ErrorOnNoConverge = false;

@@ -127,7 +119,7 @@ NAMESPACE_BEGIN(Grid);
      while(norm * inner_tol * inner_tol < stop) inner_tol *= 2;  // inner_tol = sqrt(stop/norm) ??

      PrecChangeTimer.Start();
-      precisionChange(src_f, src_d, wk_sp_from_dp);
+      precisionChange(src_f, src_d);
      PrecChangeTimer.Stop();
      
      sol_f = Zero();
@@ -137,7 +129,6 @@ NAMESPACE_BEGIN(Grid);
 	(*guesser)(src_f, sol_f);

      //Inner CG
-      std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Outer iteration " << outer_iter << " starting inner CG with tolerance " << inner_tol << std::endl;
      CG_f.Tolerance = inner_tol;
      InnerCGtimer.Start();
      CG_f(Linop_f, src_f, sol_f);
@@ -146,7 +137,7 @@ NAMESPACE_BEGIN(Grid);
      
      //Convert sol back to double and add to double prec solution
      PrecChangeTimer.Start();
-      precisionChange(tmp_d, sol_f, wk_dp_from_sp);
+      precisionChange(tmp_d, sol_f);
      PrecChangeTimer.Stop();
      
      axpy(sol_d, 1.0, tmp_d, sol_d);
@@ -158,7 +149,6 @@ NAMESPACE_BEGIN(Grid);
    ConjugateGradient<FieldD> CG_d(Tolerance, MaxInnerIterations);
    CG_d(Linop_d, src_d_in, sol_d);
    TotalFinalStepIterations = CG_d.IterationsToComplete;
-    TrueResidual = CG_d.TrueResidual;

    TotalTimer.Stop();
    std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Inner CG iterations " << TotalInnerIterations << " Restarts " << TotalOuterIterations << " Final CG iterations " << TotalFinalStepIterations << std::endl;
@@ -52,7 +52,7 @@ public:
  MultiShiftFunction shifts;
  std::vector<RealD> TrueResidualShift;

-  ConjugateGradientMultiShift(Integer maxit, const MultiShiftFunction &_shifts) : 
+  ConjugateGradientMultiShift(Integer maxit,MultiShiftFunction &_shifts) : 
    MaxIterations(maxit),
    shifts(_shifts)
  { 
@@ -182,9 +182,6 @@ public:
    for(int s=0;s<nshift;s++) {
      axpby(psi[s],0.,-bs[s]*alpha[s],src,src);
    }
-
-    std::cout << GridLogIterative << "ConjugateGradientMultiShift: initial rn (|src|^2) =" << rn << " qq (|MdagM src|^2) =" << qq << " d ( dot(src, [MdagM + m_0]src) ) =" << d << " c=" << c << std::endl;
-    
  
  ///////////////////////////////////////
  // Timers
@@ -1,411 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/algorithms/iterative/ConjugateGradientMultiShift.h
-
-    Copyright (C) 2015
-
-Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Christopher Kelly <ckelly@bnl.gov>
-
-    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_CONJUGATE_GRADIENT_MULTI_SHIFT_MIXEDPREC_H
-#define GRID_CONJUGATE_GRADIENT_MULTI_SHIFT_MIXEDPREC_H
-
-NAMESPACE_BEGIN(Grid);
-
-//CK 2020: A variant of the multi-shift conjugate gradient with the matrix multiplication in single precision. 
-//The residual is stored in single precision, but the search directions and solution are stored in double precision. 
-//Every update_freq iterations the residual is corrected in double precision. 
-    
-//For safety the a final regular CG is applied to clean up if necessary
-
-//Linop to add shift to input linop, used in cleanup CG
-namespace ConjugateGradientMultiShiftMixedPrecSupport{
-template<typename Field>
-class ShiftedLinop: public LinearOperatorBase<Field>{
-public:
-  LinearOperatorBase<Field> &linop_base;
-  RealD shift;
-
-  ShiftedLinop(LinearOperatorBase<Field> &_linop_base, RealD _shift): linop_base(_linop_base), shift(_shift){}
-
-  void OpDiag (const Field &in, Field &out){ assert(0); }
-  void OpDir  (const Field &in, Field &out,int dir,int disp){ assert(0); }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){ assert(0); }
-  
-  void Op     (const Field &in, Field &out){ assert(0); }
-  void AdjOp  (const Field &in, Field &out){ assert(0); }
-
-  void HermOp(const Field &in, Field &out){
-    linop_base.HermOp(in, out);
-    axpy(out, shift, in, out);
-  }    
-
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-    HermOp(in,out);
-    ComplexD dot = innerProduct(in,out);
-    n1=real(dot);
-    n2=norm2(out);
-  }
-};
-};
-
-
-template<class FieldD, class FieldF,
-	 typename std::enable_if< getPrecision<FieldD>::value == 2, int>::type = 0,
-	 typename std::enable_if< getPrecision<FieldF>::value == 1, int>::type = 0> 
-class ConjugateGradientMultiShiftMixedPrec : public OperatorMultiFunction<FieldD>,
-					     public OperatorFunction<FieldD>
-{
-public:                                                
-
-  using OperatorFunction<FieldD>::operator();
-
-  RealD   Tolerance;
-  Integer MaxIterations;
-  Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
-  std::vector<int> IterationsToCompleteShift;  // Iterations for this shift
-  int verbose;
-  MultiShiftFunction shifts;
-  std::vector<RealD> TrueResidualShift;
-
-  int ReliableUpdateFreq; //number of iterations between reliable updates
-
-  GridBase* SinglePrecGrid; //Grid for single-precision fields
-  LinearOperatorBase<FieldF> &Linop_f; //single precision
-
-  ConjugateGradientMultiShiftMixedPrec(Integer maxit, const MultiShiftFunction &_shifts,
-				       GridBase* _SinglePrecGrid, LinearOperatorBase<FieldF> &_Linop_f,
-				       int _ReliableUpdateFreq
-				       ) : 
-    MaxIterations(maxit),  shifts(_shifts), SinglePrecGrid(_SinglePrecGrid), Linop_f(_Linop_f), ReliableUpdateFreq(_ReliableUpdateFreq)
-  { 
-    verbose=1;
-    IterationsToCompleteShift.resize(_shifts.order);
-    TrueResidualShift.resize(_shifts.order);
-  }
-
-  void operator() (LinearOperatorBase<FieldD> &Linop, const FieldD &src, FieldD &psi)
-  {
-    GridBase *grid = src.Grid();
-    int nshift = shifts.order;
-    std::vector<FieldD> results(nshift,grid);
-    (*this)(Linop,src,results,psi);
-  }
-  void operator() (LinearOperatorBase<FieldD> &Linop, const FieldD &src, std::vector<FieldD> &results, FieldD &psi)
-  {
-    int nshift = shifts.order;
-
-    (*this)(Linop,src,results);
-  
-    psi = shifts.norm*src;
-    for(int i=0;i<nshift;i++){
-      psi = psi + shifts.residues[i]*results[i];
-    }
-
-    return;
-  }
-
-  void operator() (LinearOperatorBase<FieldD> &Linop_d, const FieldD &src_d, std::vector<FieldD> &psi_d)
-  { 
-    GridBase *DoublePrecGrid = src_d.Grid();
-    precisionChangeWorkspace wk_f_from_d(SinglePrecGrid, DoublePrecGrid);
-    precisionChangeWorkspace wk_d_from_f(DoublePrecGrid, SinglePrecGrid);
-
-    ////////////////////////////////////////////////////////////////////////
-    // Convenience references to the info stored in "MultiShiftFunction"
-    ////////////////////////////////////////////////////////////////////////
-    int nshift = shifts.order;
-
-    std::vector<RealD> &mass(shifts.poles); // Make references to array in "shifts"
-    std::vector<RealD> &mresidual(shifts.tolerances);
-    std::vector<RealD> alpha(nshift,1.0);
-
-    //Double precision search directions
-    FieldD p_d(DoublePrecGrid);
-    std::vector<FieldD> ps_d(nshift, DoublePrecGrid);// Search directions (double precision)
-
-    FieldD tmp_d(DoublePrecGrid);
-    FieldD r_d(DoublePrecGrid);
-    FieldD mmp_d(DoublePrecGrid);
-
-    assert(psi_d.size()==nshift);
-    assert(mass.size()==nshift);
-    assert(mresidual.size()==nshift);
-  
-    // dynamic sized arrays on stack; 2d is a pain with vector
-    RealD  bs[nshift];
-    RealD  rsq[nshift];
-    RealD  z[nshift][2];
-    int     converged[nshift];
-  
-    const int       primary =0;
-  
-    //Primary shift fields CG iteration
-    RealD a,b,c,d;
-    RealD cp,bp,qq; //prev
-  
-    // Matrix mult fields
-    FieldF r_f(SinglePrecGrid);
-    FieldF p_f(SinglePrecGrid);
-    FieldF tmp_f(SinglePrecGrid);
-    FieldF mmp_f(SinglePrecGrid);
-    FieldF src_f(SinglePrecGrid);
-    precisionChange(src_f, src_d, wk_f_from_d);
-
-    // Check lightest mass
-    for(int s=0;s<nshift;s++){
-      assert( mass[s]>= mass[primary] );
-      converged[s]=0;
-    }
-  
-    // Wire guess to zero
-    // Residuals "r" are src
-    // First search direction "p" is also src
-    cp = norm2(src_d);
-
-    // Handle trivial case of zero src.
-    if( cp == 0. ){
-      for(int s=0;s<nshift;s++){
-	psi_d[s] = Zero();
-	IterationsToCompleteShift[s] = 1;
-	TrueResidualShift[s] = 0.;
-      }
-      return;
-    }
-
-    for(int s=0;s<nshift;s++){
-      rsq[s] = cp * mresidual[s] * mresidual[s];
-      std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec: shift "<< s <<" target resid "<<rsq[s]<<std::endl;
-      ps_d[s] = src_d;
-    }
-    // r and p for primary
-    r_f=src_f; //residual maintained in single
-    p_f=src_f;
-    p_d = src_d; //primary copy --- make this a reference to ps_d to save axpys
-  
-    //MdagM+m[0]
-    Linop_f.HermOpAndNorm(p_f,mmp_f,d,qq); // mmp = MdagM p        d=real(dot(p, mmp)),  qq=norm2(mmp)
-    axpy(mmp_f,mass[0],p_f,mmp_f);
-    RealD rn = norm2(p_f);
-    d += rn*mass[0];
-
-    b = -cp /d;
-  
-    // Set up the various shift variables
-    int       iz=0;
-    z[0][1-iz] = 1.0;
-    z[0][iz]   = 1.0;
-    bs[0]      = b;
-    for(int s=1;s<nshift;s++){
-      z[s][1-iz] = 1.0;
-      z[s][iz]   = 1.0/( 1.0 - b*(mass[s]-mass[0]));
-      bs[s]      = b*z[s][iz]; 
-    }
-  
-    // r += b[0] A.p[0]
-    // c= norm(r)
-    c=axpy_norm(r_f,b,mmp_f,r_f);
-  
-    for(int s=0;s<nshift;s++) {
-      axpby(psi_d[s],0.,-bs[s]*alpha[s],src_d,src_d);
-    }
-  
-    ///////////////////////////////////////
-    // Timers
-    ///////////////////////////////////////
-    GridStopWatch AXPYTimer, ShiftTimer, QRTimer, MatrixTimer, SolverTimer, PrecChangeTimer, CleanupTimer;
-
-    SolverTimer.Start();
-  
-    // Iteration loop
-    int k;
-  
-    for (k=1;k<=MaxIterations;k++){    
-      a = c /cp;
-
-      //Update double precision search direction by residual
-      PrecChangeTimer.Start();
-      precisionChange(r_d, r_f, wk_d_from_f);
-      PrecChangeTimer.Stop();
-
-      AXPYTimer.Start();
-      axpy(p_d,a,p_d,r_d); 
-
-      for(int s=0;s<nshift;s++){
-	if ( ! converged[s] ) { 
-	  if (s==0){
-	    axpy(ps_d[s],a,ps_d[s],r_d);
-	  } else{
-	    RealD as =a *z[s][iz]*bs[s] /(z[s][1-iz]*b);
-	    axpby(ps_d[s],z[s][iz],as,r_d,ps_d[s]);
-	  }
-	}
-      }
-      AXPYTimer.Stop();
-
-      PrecChangeTimer.Start();
-      precisionChange(p_f, p_d, wk_f_from_d); //get back single prec search direction for linop
-      PrecChangeTimer.Stop();
-
-      cp=c;
-      MatrixTimer.Start();  
-      Linop_f.HermOp(p_f,mmp_f); 
-      d=real(innerProduct(p_f,mmp_f));    
-      MatrixTimer.Stop();  
-
-      AXPYTimer.Start();
-      axpy(mmp_f,mass[0],p_f,mmp_f);
-      AXPYTimer.Stop();
-      RealD rn = norm2(p_f);
-      d += rn*mass[0];
-    
-      bp=b;
-      b=-cp/d;
-    
-      // Toggle the recurrence history
-      bs[0] = b;
-      iz = 1-iz;
-      ShiftTimer.Start();
-      for(int s=1;s<nshift;s++){
-	if((!converged[s])){
-	  RealD z0 = z[s][1-iz];
-	  RealD z1 = z[s][iz];
-	  z[s][iz] = z0*z1*bp
-	    / (b*a*(z1-z0) + z1*bp*(1- (mass[s]-mass[0])*b)); 
-	  bs[s] = b*z[s][iz]/z0; // NB sign  rel to Mike
-	}
-      }
-      ShiftTimer.Stop();
-
-      //Update double precision solutions
-      AXPYTimer.Start();
-      for(int s=0;s<nshift;s++){
-	int ss = s;
-	if( (!converged[s]) ) { 
-	  axpy(psi_d[ss],-bs[s]*alpha[s],ps_d[s],psi_d[ss]);
-	}
-      }
-
-      //Perform reliable update if necessary; otherwise update residual from single-prec mmp
-      RealD c_f = axpy_norm(r_f,b,mmp_f,r_f);
-      AXPYTimer.Stop();
-
-      c = c_f;
-
-      if(k % ReliableUpdateFreq == 0){
-	//Replace r with true residual
-	MatrixTimer.Start();  
-	Linop_d.HermOp(psi_d[0],mmp_d); 
-	MatrixTimer.Stop();  
-
-	AXPYTimer.Start();
-	axpy(mmp_d,mass[0],psi_d[0],mmp_d);
-
-	RealD c_d = axpy_norm(r_d, -1.0, mmp_d, src_d);
-	AXPYTimer.Stop();
-
-	std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec k="<<k<< ", replaced |r|^2 = "<<c_f <<" with |r|^2 = "<<c_d<<std::endl;
-	
-	PrecChangeTimer.Start();
-	precisionChange(r_f, r_d, wk_f_from_d);
-	PrecChangeTimer.Stop();
-	c = c_d;
-      }
-    
-      // Convergence checks
-      int all_converged = 1;
-      for(int s=0;s<nshift;s++){
-      
-	if ( (!converged[s]) ){
-	  IterationsToCompleteShift[s] = k;
-	
-	  RealD css  = c * z[s][iz]* z[s][iz];
-	
-	  if(css<rsq[s]){
-	    if ( ! converged[s] )
-	      std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec k="<<k<<" Shift "<<s<<" has converged"<<std::endl;
-	    converged[s]=1;
-	  } else {
-	    all_converged=0;
-	  }
-
-	}
-      }
-
-      if ( all_converged ){
-
-	SolverTimer.Stop();
-	std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrec: All shifts have converged iteration "<<k<<std::endl;
-	std::cout<<GridLogMessage<< "ConjugateGradientMultiShiftMixedPrec: Checking solutions"<<std::endl;
-      
-	// Check answers 
-	for(int s=0; s < nshift; s++) { 
-	  Linop_d.HermOpAndNorm(psi_d[s],mmp_d,d,qq);
-	  axpy(tmp_d,mass[s],psi_d[s],mmp_d);
-	  axpy(r_d,-alpha[s],src_d,tmp_d);
-	  RealD rn = norm2(r_d);
-	  RealD cn = norm2(src_d);
-	  TrueResidualShift[s] = std::sqrt(rn/cn);
-	  std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec: shift["<<s<<"] true residual "<< TrueResidualShift[s] << " target " << mresidual[s] << std::endl;
-
-	  //If we have not reached the desired tolerance, do a (mixed precision) CG cleanup
-	  if(rn >= rsq[s]){
-	    CleanupTimer.Start();
-	    std::cout<<GridLogMessage<<"ConjugateGradientMultiShiftMixedPrec: performing cleanup step for shift " << s << std::endl;
-
-	    //Setup linear operators for final cleanup
-	    ConjugateGradientMultiShiftMixedPrecSupport::ShiftedLinop<FieldD> Linop_shift_d(Linop_d, mass[s]);
-	    ConjugateGradientMultiShiftMixedPrecSupport::ShiftedLinop<FieldF> Linop_shift_f(Linop_f, mass[s]);
-					       
-	    MixedPrecisionConjugateGradient<FieldD,FieldF> cg(mresidual[s], MaxIterations, MaxIterations, SinglePrecGrid, Linop_shift_f, Linop_shift_d); 
-	    cg(src_d, psi_d[s]);
-	    
-	    TrueResidualShift[s] = cg.TrueResidual;
-	    CleanupTimer.Stop();
-	  }
-	}
-
-	std::cout << GridLogMessage << "ConjugateGradientMultiShiftMixedPrec: Time Breakdown for body"<<std::endl;
-	std::cout << GridLogMessage << "\tSolver    " << SolverTimer.Elapsed()     <<std::endl;
-	std::cout << GridLogMessage << "\t\tAXPY    " << AXPYTimer.Elapsed()     <<std::endl;
-	std::cout << GridLogMessage << "\t\tMatrix    " << MatrixTimer.Elapsed()     <<std::endl;
-	std::cout << GridLogMessage << "\t\tShift    " << ShiftTimer.Elapsed()     <<std::endl;
-	std::cout << GridLogMessage << "\t\tPrecision Change " << PrecChangeTimer.Elapsed()     <<std::endl;
-	std::cout << GridLogMessage << "\tFinal Cleanup " << CleanupTimer.Elapsed()     <<std::endl;
-	std::cout << GridLogMessage << "\tSolver+Cleanup " << SolverTimer.Elapsed() + CleanupTimer.Elapsed() << std::endl;
-
-	IterationsToComplete = k;	
-
-	return;
-      }
-
-   
-    }
-    // ugly hack
-    std::cout<<GridLogMessage<<"CG multi shift did not converge"<<std::endl;
-    //  assert(0);
-  }
-
-};
-NAMESPACE_END(Grid);
-#endif
@@ -44,7 +44,6 @@ public:
 				  int, MinRes);    // Must restart
 };

-//This class is the input parameter class for some testing programs
 struct LocalCoherenceLanczosParams : Serializable {
 public:
  GRID_SERIALIZABLE_CLASS_MEMBERS(LocalCoherenceLanczosParams,
@@ -154,7 +153,6 @@ public:
      _coarse_relax_tol(coarse_relax_tol)  
  {    };

-  //evalMaxApprox: approximation of largest eval of the fine Chebyshev operator (suitably wrapped by block projection)
  int TestConvergence(int j,RealD eresid,CoarseField &B, RealD &eval,RealD evalMaxApprox)
  {
    CoarseField v(B);
@@ -181,16 +179,8 @@ public:
    if( (vv<eresid*eresid) ) conv = 1;
    return conv;
  }
-
-  //This function is called at the end of the coarse grid Lanczos. It promotes the coarse eigenvector 'B' to the fine grid,
-  //applies a smoother to the result then computes the computes the *fine grid* eigenvalue (output as 'eval').
-
-  //evalMaxApprox should be the approximation of the largest eval of the fine Hermop. However when this function is called by IRL it actually passes the largest eval of the *Chebyshev* operator (as this is the max approx used for the TestConvergence above)
-  //As the largest eval of the Chebyshev is typically several orders of magnitude larger this makes the convergence test pass even when it should not.
-  //We therefore ignore evalMaxApprox here and use a value of 1.0 (note this value is already used by TestCoarse)
-  int ReconstructEval(int j,RealD eresid,CoarseField &B, RealD &eval,RealD evalMaxApprox)  
+  int ReconstructEval(int j,RealD eresid,CoarseField &B, RealD &eval,RealD evalMaxApprox)
  {
-    evalMaxApprox = 1.0; //cf above
    GridBase *FineGrid = _subspace[0].Grid();    
    int checkerboard   = _subspace[0].Checkerboard();
    FineField fB(FineGrid);fB.Checkerboard() =checkerboard;
@@ -209,13 +199,13 @@ public:
    eval   = vnum/vden;
    fv -= eval*fB;
    RealD vv = norm2(fv) / ::pow(evalMaxApprox,2.0);
-    if ( j > nbasis ) eresid = eresid*_coarse_relax_tol;
-    
+
    std::cout.precision(13);
    std::cout<<GridLogIRL  << "[" << std::setw(3)<<j<<"] "
 	     <<"eval = "<<std::setw(25)<< eval << " (" << eval_poly << ")"
-	     <<" |H B[i] - eval[i]B[i]|^2 / evalMaxApprox^2 " << std::setw(25) << vv << " target " << eresid*eresid
+	     <<" |H B[i] - eval[i]B[i]|^2 / evalMaxApprox^2 " << std::setw(25) << vv
 	     <<std::endl;
+    if ( j > nbasis ) eresid = eresid*_coarse_relax_tol;
    if( (vv<eresid*eresid) ) return 1;
    return 0;
  }
@@ -293,10 +283,6 @@ public:
    evals_coarse.resize(0);
  };

-  //The block inner product is the inner product on the fine grid locally summed over the blocks
-  //to give a Lattice<Scalar> on the coarse grid. This function orthnormalizes the fine-grid subspace
-  //vectors under the block inner product. This step must be performed after computing the fine grid
-  //eigenvectors and before computing the coarse grid eigenvectors.    
  void Orthogonalise(void ) {
    CoarseScalar InnerProd(_CoarseGrid);
    std::cout << GridLogMessage <<" Gramm-Schmidt pass 1"<<std::endl;
@@ -340,8 +326,6 @@ public:
    }
  }

-  //While this method serves to check the coarse eigenvectors, it also recomputes the eigenvalues from the smoothed reconstructed eigenvectors
-  //hence the smoother can be tuned after running the coarse Lanczos by using a different smoother here
  void testCoarse(RealD resid,ChebyParams cheby_smooth,RealD relax) 
  {
    assert(evals_fine.size() == nbasis);
@@ -390,31 +374,25 @@ public:
    evals_fine.resize(nbasis);
    subspace.resize(nbasis,_FineGrid);
  }
-
-
-  //cheby_op: Parameters of the fine grid Chebyshev polynomial used for the Lanczos acceleration
-  //cheby_smooth: Parameters of a separate Chebyshev polynomial used after the Lanczos has completed to smooth out high frequency noise in the reconstructed fine grid eigenvectors prior to computing the eigenvalue
-  //relax: Reconstructed eigenvectors (post smoothing) are naturally not as precise as true eigenvectors. This factor acts as a multiplier on the stopping condition when determining whether the results satisfy the user provided stopping condition
  void calcCoarse(ChebyParams cheby_op,ChebyParams cheby_smooth,RealD relax,
 		  int Nstop, int Nk, int Nm,RealD resid, 
 		  RealD MaxIt, RealD betastp, int MinRes)
  {
-    Chebyshev<FineField>                          Cheby(cheby_op); //Chebyshev of fine operator on fine grid
-    ProjectedHermOp<Fobj,CComplex,nbasis>         Op(_FineOp,subspace); //Fine operator on coarse grid with intermediate fine grid conversion
-    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (Cheby,_FineOp,subspace); //Chebyshev of fine operator on coarse grid with intermediate fine grid conversion
+    Chebyshev<FineField>                          Cheby(cheby_op);
+    ProjectedHermOp<Fobj,CComplex,nbasis>         Op(_FineOp,subspace);
+    ProjectedFunctionHermOp<Fobj,CComplex,nbasis> ChebyOp (Cheby,_FineOp,subspace);
    //////////////////////////////////////////////////////////////////////////////////////////////////
    // create a smoother and see if we can get a cheap convergence test and smooth inside the IRL
    //////////////////////////////////////////////////////////////////////////////////////////////////

-    Chebyshev<FineField>                                           ChebySmooth(cheby_smooth); //lower order Chebyshev of fine operator on fine grid used to smooth regenerated eigenvectors
-    ImplicitlyRestartedLanczosSmoothedTester<Fobj,CComplex,nbasis> ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax); 
+    Chebyshev<FineField>                                           ChebySmooth(cheby_smooth);
+    ImplicitlyRestartedLanczosSmoothedTester<Fobj,CComplex,nbasis> ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax);

    evals_coarse.resize(Nm);
    evec_coarse.resize(Nm,_CoarseGrid);

    CoarseField src(_CoarseGrid);     src=1.0; 

-    //Note the "tester" here is also responsible for generating the fine grid eigenvalues which are output into the "evals_coarse" array
    ImplicitlyRestartedLanczos<CoarseField> IRL(ChebyOp,ChebyOp,ChebySmoothTester,Nstop,Nk,Nm,resid,MaxIt,betastp,MinRes);
    int Nconv=0;
    IRL.calc(evals_coarse,evec_coarse,src,Nconv,false);
@@ -425,14 +403,6 @@ public:
      std::cout << i << " Coarse eval = " << evals_coarse[i]  << std::endl;
    }
  }
-
-  //Get the fine eigenvector 'i' by reconstruction
-  void getFineEvecEval(FineField &evec, RealD &eval, const int i) const{
-    blockPromote(evec_coarse[i],evec,subspace);  
-    eval = evals_coarse[i];
-  }
-    
-    
 };

 NAMESPACE_END(Grid);
@@ -29,8 +29,6 @@ template<class Field> class PowerMethod
      RealD vnum = real(innerProduct(src_n,tmp)); // HermOp. 
      RealD vden = norm2(src_n); 
      RealD na = vnum/vden; 
-
-      std::cout << GridLogIterative << "PowerMethod: Current approximation of largest eigenvalue " << na << std::endl;
      
      if ( (fabs(evalMaxApprox/na - 1.0) < 0.001) || (i==_MAX_ITER_EST_-1) ) { 
 	evalMaxApprox = na; 
@@ -132,31 +132,6 @@ namespace Grid {
      (*this)(_Matrix,in,out,guess);
    }

-    void RedBlackSource(Matrix &_Matrix, const std::vector<Field> &in, std::vector<Field> &src_o) 
-    {
-      GridBase *grid = _Matrix.RedBlackGrid();
-      Field tmp(grid);
-      int nblock = in.size();
-      for(int b=0;b<nblock;b++){
-	RedBlackSource(_Matrix,in[b],tmp,src_o[b]);
-      }
-    }
-    // James can write his own deflated guesser
-    // with optimised code for the inner products
-    //    RedBlackSolveSplitGrid();
-    //    RedBlackSolve(_Matrix,src_o,sol_o); 
-
-    void RedBlackSolution(Matrix &_Matrix, const std::vector<Field> &in, const std::vector<Field> &sol_o, std::vector<Field> &out)
-    {
-      GridBase *grid = _Matrix.RedBlackGrid();
-      Field tmp(grid);
-      int nblock = in.size();
-      for(int b=0;b<nblock;b++) {
-	pickCheckerboard(Even,tmp,in[b]);
-	RedBlackSolution(_Matrix,sol_o[b],tmp,out[b]);
-      }
-    }
-
    template<class Guesser>
    void operator()(Matrix &_Matrix, const std::vector<Field> &in, std::vector<Field> &out,Guesser &guess) 
    {
@@ -175,29 +150,24 @@ namespace Grid {
      ////////////////////////////////////////////////
      // Prepare RedBlack source
      ////////////////////////////////////////////////
-      RedBlackSource(_Matrix,in,src_o);
-	//      for(int b=0;b<nblock;b++){
-	//	RedBlackSource(_Matrix,in[b],tmp,src_o[b]);
-	//      }
-      
+      for(int b=0;b<nblock;b++){
+	RedBlackSource(_Matrix,in[b],tmp,src_o[b]);
+      }
      ////////////////////////////////////////////////
      // Make the guesses
      ////////////////////////////////////////////////
      if ( subGuess ) guess_save.resize(nblock,grid);

-      
-      if(useSolnAsInitGuess) {
-        for(int b=0;b<nblock;b++){
+      for(int b=0;b<nblock;b++){
+        if(useSolnAsInitGuess) {
          pickCheckerboard(Odd, sol_o[b], out[b]);
+        } else {
+          guess(src_o[b],sol_o[b]); 
        }
-      } else {
-        guess(src_o, sol_o); 
-      }

-	    if ( subGuess ) { 
-        for(int b=0;b<nblock;b++){
-          guess_save[b] = sol_o[b];
-        }
+	if ( subGuess ) { 
+	  guess_save[b] = sol_o[b];
+	}
      }
      //////////////////////////////////////////////////////////////
      // Call the block solver
@@ -33,8 +33,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>

 NAMESPACE_BEGIN(Grid);

-bool Stencil_force_mpi = true;
-
 ///////////////////////////////////////////////////////////////
 // Info that is setup once and indept of cartesian layout
 ///////////////////////////////////////////////////////////////
@@ -35,8 +35,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>

 NAMESPACE_BEGIN(Grid);

-extern bool Stencil_force_mpi ;
-
 class CartesianCommunicator : public SharedMemory {

 public:    
@@ -370,7 +370,7 @@ double CartesianCommunicator::StencilSendToRecvFromBegin(std::vector<CommsReques
  double off_node_bytes=0.0;
  int tag;

-  if ( (gfrom ==MPI_UNDEFINED) || Stencil_force_mpi ) {
+  if ( gfrom ==MPI_UNDEFINED) {
    tag= dir+from*32;
    ierr=MPI_Irecv(recv, bytes, MPI_CHAR,from,tag,communicator_halo[commdir],&rrq);
    assert(ierr==0);
@@ -378,18 +378,12 @@ double CartesianCommunicator::StencilSendToRecvFromBegin(std::vector<CommsReques
    off_node_bytes+=bytes;
  }

-  if ( (gdest == MPI_UNDEFINED) || Stencil_force_mpi ) {
+  if ( gdest == MPI_UNDEFINED ) {
    tag= dir+_processor*32;
    ierr =MPI_Isend(xmit, bytes, MPI_CHAR,dest,tag,communicator_halo[commdir],&xrq);
    assert(ierr==0);
    list.push_back(xrq);
    off_node_bytes+=bytes;
-  } else {
-    // TODO : make a OMP loop on CPU, call threaded bcopy
-    void *shm = (void *) this->ShmBufferTranslate(dest,recv);
-    assert(shm!=NULL);
-    acceleratorCopyDeviceToDeviceAsynch(xmit,shm,bytes);
-    acceleratorCopySynchronise(); // MPI prob slower
  }

  if ( CommunicatorPolicy == CommunicatorPolicySequential ) {
@@ -35,9 +35,6 @@ Author: Christoph Lehner <christoph@lhnr.de>
 #endif
 #ifdef GRID_HIP
 #include <hip/hip_runtime_api.h>
-#endif
-#ifdef GRID_SYCl
-
 #endif

 NAMESPACE_BEGIN(Grid); 
@@ -73,7 +70,6 @@ void GlobalSharedMemory::Init(Grid_MPI_Comm comm)
  WorldNodes = WorldSize/WorldShmSize;
  assert( (WorldNodes * WorldShmSize) == WorldSize );

-
  // FIXME: Check all WorldShmSize are the same ?

  /////////////////////////////////////////////////////////////////////
@@ -450,47 +446,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 ////////////////////////////////////////////////////////////////////////////////////////////
 // Hugetlbfs mapping intended
 ////////////////////////////////////////////////////////////////////////////////////////////
-#if defined(GRID_CUDA) ||defined(GRID_HIP)  || defined(GRID_SYCL)
-
-//if defined(GRID_SYCL)
-#if 0
-void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
-{
-  void * ShmCommBuf ; 
-  assert(_ShmSetup==1);
-  assert(_ShmAlloc==0);
-
-  //////////////////////////////////////////////////////////////////////////////////////////////////////////
-  // allocate the pointer array for shared windows for our group
-  //////////////////////////////////////////////////////////////////////////////////////////////////////////
-  MPI_Barrier(WorldShmComm);
-  WorldShmCommBufs.resize(WorldShmSize);
-
-  ///////////////////////////////////////////////////////////////////////////////////////////////////////////
-  // Each MPI rank should allocate our own buffer
-  ///////////////////////////////////////////////////////////////////////////////////////////////////////////
-  ShmCommBuf = acceleratorAllocDevice(bytes);
-
-  if (ShmCommBuf == (void *)NULL ) {
-    std::cerr << " SharedMemoryMPI.cc acceleratorAllocDevice failed NULL pointer for " << bytes<<" bytes " << std::endl;
-    exit(EXIT_FAILURE);  
-  }
-
-  std::cout << WorldRank << header " SharedMemoryMPI.cc acceleratorAllocDevice "<< bytes 
-	    << "bytes at "<< std::hex<< ShmCommBuf <<std::dec<<" for comms buffers " <<std::endl;
-
-  SharedMemoryZero(ShmCommBuf,bytes);
-
-  assert(WorldShmSize == 1);
-  for(int r=0;r<WorldShmSize;r++){
-    WorldShmCommBufs[r] = ShmCommBuf;
-  }
-  _ShmAllocBytes=bytes;
-  _ShmAlloc=1;
-}
-#endif
-
-#if defined(GRID_CUDA) ||defined(GRID_HIP) ||defined(GRID_SYCL)  
+#if defined(GRID_CUDA) ||defined(GRID_HIP)
 void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 {
  void * ShmCommBuf ; 
@@ -514,16 +470,18 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
  // Each MPI rank should allocate our own buffer
  ///////////////////////////////////////////////////////////////////////////////////////////////////////////
  ShmCommBuf = acceleratorAllocDevice(bytes);
+
  if (ShmCommBuf == (void *)NULL ) {
    std::cerr << " SharedMemoryMPI.cc acceleratorAllocDevice failed NULL pointer for " << bytes<<" bytes " << std::endl;
    exit(EXIT_FAILURE);  
  }
-  if ( WorldRank == 0 ){
+  //  if ( WorldRank == 0 ){
+  if ( 1 ){
    std::cout << WorldRank << header " SharedMemoryMPI.cc acceleratorAllocDevice "<< bytes 
 	      << "bytes at "<< std::hex<< ShmCommBuf <<std::dec<<" for comms buffers " <<std::endl;
  }
  SharedMemoryZero(ShmCommBuf,bytes);
-  std::cout<< "Setting up IPC"<<std::endl;
+
  ///////////////////////////////////////////////////////////////////////////////////////////////////////////
  // Loop over ranks/gpu's on our node
  ///////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -533,29 +491,6 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
    //////////////////////////////////////////////////
    // If it is me, pass around the IPC access key
    //////////////////////////////////////////////////
-    void * thisBuf = ShmCommBuf;
-    if(!Stencil_force_mpi) {
-#ifdef GRID_SYCL_LEVEL_ZERO_IPC
-    typedef struct { int fd; pid_t pid ; } clone_mem_t;
-
-    auto zeDevice    = cl::sycl::get_native<cl::sycl::backend::level_zero>(theGridAccelerator->get_device());
-    auto zeContext   = cl::sycl::get_native<cl::sycl::backend::level_zero>(theGridAccelerator->get_context());
-      
-    ze_ipc_mem_handle_t ihandle;
-    clone_mem_t handle;
-
-    if ( r==WorldShmRank ) { 
-      auto err = zeMemGetIpcHandle(zeContext,ShmCommBuf,&ihandle);
-      if ( err != ZE_RESULT_SUCCESS ) {
-	std::cout << "SharedMemoryMPI.cc zeMemGetIpcHandle failed for rank "<<r<<" "<<std::hex<<err<<std::dec<<std::endl;
-	exit(EXIT_FAILURE);
-      } else {
-	std::cout << "SharedMemoryMPI.cc zeMemGetIpcHandle succeeded for rank "<<r<<" "<<std::hex<<err<<std::dec<<std::endl;
-      }
-      memcpy((void *)&handle.fd,(void *)&ihandle,sizeof(int));
-      handle.pid = getpid();
-    }
-#endif
 #ifdef GRID_CUDA
    cudaIpcMemHandle_t handle;
    if ( r==WorldShmRank ) { 
@@ -576,7 +511,6 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
      }
    }
 #endif
-
    //////////////////////////////////////////////////
    // Share this IPC handle across the Shm Comm
    //////////////////////////////////////////////////
@@ -592,35 +526,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
    ///////////////////////////////////////////////////////////////
    // If I am not the source, overwrite thisBuf with remote buffer
    ///////////////////////////////////////////////////////////////
-
-#ifdef GRID_SYCL_LEVEL_ZERO_IPC
-    if ( r!=WorldShmRank ) {
-      thisBuf = nullptr;
-      std::cout<<"mapping seeking remote pid/fd "
-	       <<handle.pid<<"/"
-	       <<handle.fd<<std::endl;
-
-      int pidfd = syscall(SYS_pidfd_open,handle.pid,0);
-      std::cout<<"Using IpcHandle pidfd "<<pidfd<<"\n";
-      //      int myfd  = syscall(SYS_pidfd_getfd,pidfd,handle.fd,0);
-      int myfd  = syscall(438,pidfd,handle.fd,0);
-
-      std::cout<<"Using IpcHandle myfd "<<myfd<<"\n";
-      
-      memcpy((void *)&ihandle,(void *)&myfd,sizeof(int));
-
-      auto err = zeMemOpenIpcHandle(zeContext,zeDevice,ihandle,0,&thisBuf);
-      if ( err != ZE_RESULT_SUCCESS ) {
-	std::cout << "SharedMemoryMPI.cc "<<zeContext<<" "<<zeDevice<<std::endl;
-	std::cout << "SharedMemoryMPI.cc zeMemOpenIpcHandle failed for rank "<<r<<" "<<std::hex<<err<<std::dec<<std::endl; 
-	exit(EXIT_FAILURE);
-      } else {
-	std::cout << "SharedMemoryMPI.cc zeMemOpenIpcHandle succeeded for rank "<<r<<std::endl;
-	std::cout << "SharedMemoryMPI.cc zeMemOpenIpcHandle pointer is "<<std::hex<<thisBuf<<std::dec<<std::endl;
-      }
-      assert(thisBuf!=nullptr);
-    }
-#endif
+    void * thisBuf = ShmCommBuf;
 #ifdef GRID_CUDA
    if ( r!=WorldShmRank ) { 
      auto err = cudaIpcOpenMemHandle(&thisBuf,handle,cudaIpcMemLazyEnablePeerAccess);
@@ -642,7 +548,6 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
    ///////////////////////////////////////////////////////////////
    // Save a copy of the device buffers
    ///////////////////////////////////////////////////////////////
-    }
    WorldShmCommBufs[r] = thisBuf;
 #else
    WorldShmCommBufs[r] = ShmCommBuf;
@@ -652,8 +557,6 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
  _ShmAllocBytes=bytes;
  _ShmAlloc=1;
 }
-#endif
-
 #else 
 #ifdef GRID_MPI3_SHMMMAP
 void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
@@ -824,16 +727,16 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 /////////////////////////////////////////////////////////////////////////
 void GlobalSharedMemory::SharedMemoryZero(void *dest,size_t bytes)
 {
-#if defined(GRID_CUDA) || defined(GRID_HIP) || defined(GRID_SYCL)
-  acceleratorMemSet(dest,0,bytes);
+#ifdef GRID_CUDA
+  cudaMemset(dest,0,bytes);
 #else
  bzero(dest,bytes);
 #endif
 }
 void GlobalSharedMemory::SharedMemoryCopy(void *dest,void *src,size_t bytes)
 {
-#if defined(GRID_CUDA) || defined(GRID_HIP) || defined(GRID_SYCL)
-  acceleratorCopyToDevice(src,dest,bytes);
+#ifdef GRID_CUDA
+  cudaMemcpy(dest,src,bytes,cudaMemcpyDefault);
 #else   
  bcopy(src,dest,bytes);
 #endif
@@ -897,7 +800,7 @@ void SharedMemory::SetCommunicator(Grid_MPI_Comm comm)
  }
 #endif

-  //SharedMemoryTest();
+  SharedMemoryTest();
 }
 //////////////////////////////////////////////////////////////////
 // On node barrier
@@ -46,4 +46,3 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/lattice/Lattice_unary.h>
 #include <Grid/lattice/Lattice_transfer.h>
 #include <Grid/lattice/Lattice_basis.h>
-#include <Grid/lattice/Lattice_crc.h>
@@ -225,7 +225,7 @@ void axpy(Lattice<vobj> &ret,sobj a,const Lattice<vobj> &x,const Lattice<vobj> &
  autoView( x_v , x, AcceleratorRead);
  autoView( y_v , y, AcceleratorRead);
  accelerator_for(ss,x_v.size(),vobj::Nsimd(),{
-    auto tmp = a*coalescedRead(x_v[ss])+coalescedRead(y_v[ss]);
+    auto tmp = a*x_v(ss)+y_v(ss);
    coalescedWrite(ret_v[ss],tmp);
  });
 }
@@ -125,7 +125,7 @@ void basisRotate(VField &basis,Matrix& Qt,int j0, int j1, int k0,int k1,int Nm)

 	for(int k=k0; k<k1; ++k){
 	  auto tmp = coalescedRead(Bp[ss*nrot+j]);
-	  coalescedWrite(Bp[ss*nrot+j],tmp+ Qt_p[jj*Nm+k] * coalescedRead(basis_vp[k][sss]));
+	  coalescedWrite(Bp[ss*nrot+j],tmp+ Qt_p[jj*Nm+k] * coalescedRead(basis_v[k][sss]));
 	}
      });

@@ -134,7 +134,7 @@ void basisRotate(VField &basis,Matrix& Qt,int j0, int j1, int k0,int k1,int Nm)
 	int jj  =j0+j;
 	int ss =sj/nrot;
 	int sss=ss+s;
-	coalescedWrite(basis_vp[jj][sss],coalescedRead(Bp[ss*nrot+j]));
+	coalescedWrite(basis_v[jj][sss],coalescedRead(Bp[ss*nrot+j]));
      });
  }
 #endif
@@ -1,42 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/lattice/Lattice_crc.h
-
-    Copyright (C) 2021
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#pragma once
-
-NAMESPACE_BEGIN(Grid);
-
-template<class vobj> uint32_t crc(Lattice<vobj> & buf)
-{
-  autoView( buf_v , buf, CpuRead);
-  return ::crc32(0L,(unsigned char *)&buf_v[0],(size_t)sizeof(vobj)*buf.oSites());
-}
-
-#define CRC(U) std::cout << "FingerPrint "<<__FILE__ <<" "<< __LINE__ <<" "<< #U <<" "<<crc(U)<<std::endl;
-
-NAMESPACE_END(Grid);
-
-
@@ -361,7 +361,6 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
  // But easily avoided by using double precision fields
  ///////////////////////////////////////////////////////
  typedef typename vobj::scalar_object sobj;
-  typedef typename vobj::scalar_object::scalar_type scalar_type;
  GridBase  *grid = Data.Grid();
  assert(grid!=NULL);

@@ -420,19 +419,20 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
  }
  
  // sum over nodes.
+  sobj gsum;
  for(int t=0;t<fd;t++){
    int pt = t/ld; // processor plane
    int lt = t%ld;
    if ( pt == grid->_processor_coor[orthogdim] ) {
-      result[t]=lsSum[lt];
+      gsum=lsSum[lt];
    } else {
-      result[t]=Zero();
+      gsum=Zero();
    }

+    grid->GlobalSum(gsum);
+
+    result[t]=gsum;
  }
-  scalar_type * ptr = (scalar_type *) &result[0];
-  int words = fd*sizeof(sobj)/sizeof(scalar_type);
-  grid->GlobalSumVector(ptr, words);
 }

 template<class vobj>
@@ -32,9 +32,8 @@
 #include <random>

 #ifdef RNG_SITMO
-#include <Grid/random/sitmo_prng_engine.hpp>
+#include <Grid/sitmo_rng/sitmo_prng_engine.hpp>
 #endif 
-#include <Grid/random/gaussian.h>

 #if defined(RNG_SITMO)
 #define RNG_FAST_DISCARD
@@ -143,8 +142,8 @@ public:

  std::vector<RngEngine>                             _generators;
  std::vector<std::uniform_real_distribution<RealD> > _uniform;
-  std::vector<Grid::gaussian_distribution<RealD> >       _gaussian;
-  //  std::vector<std::discrete_distribution<int32_t> >   _bernoulli;
+  std::vector<std::normal_distribution<RealD> >       _gaussian;
+  std::vector<std::discrete_distribution<int32_t> >   _bernoulli;
  std::vector<std::uniform_int_distribution<uint32_t> > _uid;

  ///////////////////////
@@ -244,8 +243,8 @@ public:
  GridSerialRNG() : GridRNGbase() {
    _generators.resize(1);
    _uniform.resize(1,std::uniform_real_distribution<RealD>{0,1});
-    _gaussian.resize(1,gaussian_distribution<RealD>(0.0,1.0) );
-    //    _bernoulli.resize(1,std::discrete_distribution<int32_t>{1,1});
+    _gaussian.resize(1,std::normal_distribution<RealD>(0.0,1.0) );
+    _bernoulli.resize(1,std::discrete_distribution<int32_t>{1,1});
    _uid.resize(1,std::uniform_int_distribution<uint32_t>() );
  }

@@ -358,8 +357,8 @@ public:

    _generators.resize(_vol);
    _uniform.resize(_vol,std::uniform_real_distribution<RealD>{0,1});
-    _gaussian.resize(_vol,gaussian_distribution<RealD>(0.0,1.0) );
-    //    _bernoulli.resize(_vol,std::discrete_distribution<int32_t>{1,1});
+    _gaussian.resize(_vol,std::normal_distribution<RealD>(0.0,1.0) );
+    _bernoulli.resize(_vol,std::discrete_distribution<int32_t>{1,1});
    _uid.resize(_vol,std::uniform_int_distribution<uint32_t>() );
  }

@@ -516,11 +515,11 @@ public:

 template <class vobj> inline void random(GridParallelRNG &rng,Lattice<vobj> &l)   { rng.fill(l,rng._uniform);  }
 template <class vobj> inline void gaussian(GridParallelRNG &rng,Lattice<vobj> &l) { rng.fill(l,rng._gaussian); }
-//template <class vobj> inline void bernoulli(GridParallelRNG &rng,Lattice<vobj> &l){ rng.fill(l,rng._bernoulli);}
+template <class vobj> inline void bernoulli(GridParallelRNG &rng,Lattice<vobj> &l){ rng.fill(l,rng._bernoulli);}

 template <class sobj> inline void random(GridSerialRNG &rng,sobj &l)   { rng.fill(l,rng._uniform  ); }
 template <class sobj> inline void gaussian(GridSerialRNG &rng,sobj &l) { rng.fill(l,rng._gaussian ); }
-//template <class sobj> inline void bernoulli(GridSerialRNG &rng,sobj &l){ rng.fill(l,rng._bernoulli); }
+template <class sobj> inline void bernoulli(GridSerialRNG &rng,sobj &l){ rng.fill(l,rng._bernoulli); }

 NAMESPACE_END(Grid);
 #endif
@@ -364,22 +364,16 @@ inline void blockSum(Lattice<vobj> &coarseData,const Lattice<vobj> &fineData)
  autoView( coarseData_ , coarseData, AcceleratorWrite);
  autoView( fineData_   , fineData, AcceleratorRead);

-  auto coarseData_p = &coarseData_[0];
-  auto fineData_p = &fineData_[0];
-  
  Coordinate fine_rdimensions = fine->_rdimensions;
  Coordinate coarse_rdimensions = coarse->_rdimensions;
-
-  vobj zz = Zero();
  
  accelerator_for(sc,coarse->oSites(),1,{

      // One thread per sub block
      Coordinate coor_c(_ndimension);
      Lexicographic::CoorFromIndex(coor_c,sc,coarse_rdimensions);  // Block coordinate
+      coarseData_[sc]=Zero();

-      vobj cd = zz;
-      
      for(int sb=0;sb<blockVol;sb++){

 	int sf;
@@ -389,11 +383,9 @@ inline void blockSum(Lattice<vobj> &coarseData,const Lattice<vobj> &fineData)
 	for(int d=0;d<_ndimension;d++) coor_f[d]=coor_c[d]*block_r[d] + coor_b[d];
 	Lexicographic::IndexFromCoor(coor_f,sf,fine_rdimensions);

-	cd=cd+fineData_p[sf];
+	coarseData_[sc]=coarseData_[sc]+fineData_[sf];
      }

-      coarseData_p[sc] = cd;
-
    });
  return;
 }
@@ -785,7 +777,7 @@ void ExtractSliceLocal(Lattice<vobj> &lowDim,const Lattice<vobj> & higherDim,int


 template<class vobj>
-void Replicate(const Lattice<vobj> &coarse,Lattice<vobj> & fine)
+void Replicate(Lattice<vobj> &coarse,Lattice<vobj> & fine)
 {
  typedef typename vobj::scalar_object sobj;

@@ -1010,95 +1002,53 @@ vectorizeFromRevLexOrdArray( std::vector<sobj> &in, Lattice<vobj> &out)
  });
 }

-//The workspace for a precision change operation allowing for the reuse of the mapping to save time on subsequent calls
-class precisionChangeWorkspace{
-  std::pair<Integer,Integer>* fmap_device; //device pointer
-public:
-  precisionChangeWorkspace(GridBase *out_grid, GridBase *in_grid){
-    //Build a map between the sites and lanes of the output field and the input field as we cannot use the Grids on the device
-    assert(out_grid->Nd() == in_grid->Nd());
-    for(int d=0;d<out_grid->Nd();d++){
-      assert(out_grid->FullDimensions()[d] == in_grid->FullDimensions()[d]);
-    }
-    int Nsimd_out = out_grid->Nsimd();
-
-    std::vector<Coordinate> out_icorrs(out_grid->Nsimd()); //reuse these
-    for(int lane=0; lane < out_grid->Nsimd(); lane++)
-      out_grid->iCoorFromIindex(out_icorrs[lane], lane);
-  
-    std::vector<std::pair<Integer,Integer> > fmap_host(out_grid->lSites()); //lsites = osites*Nsimd
-    thread_for(out_oidx,out_grid->oSites(),{
-	Coordinate out_ocorr; 
-	out_grid->oCoorFromOindex(out_ocorr, out_oidx);
-      
-	Coordinate lcorr; //the local coordinate (common to both in and out as full coordinate)
-	for(int out_lane=0; out_lane < Nsimd_out; out_lane++){
-	  out_grid->InOutCoorToLocalCoor(out_ocorr, out_icorrs[out_lane], lcorr);
-	
-	  //int in_oidx = in_grid->oIndex(lcorr), in_lane = in_grid->iIndex(lcorr);
-	  //Note oIndex and OcorrFromOindex (and same for iIndex) are not inverse for checkerboarded lattice, the former coordinates being defined on the full lattice and the latter on the reduced lattice
-	  //Until this is fixed we need to circumvent the problem locally. Here I will use the coordinates defined on the reduced lattice for simplicity
-	  int in_oidx = 0, in_lane = 0;
-	  for(int d=0;d<in_grid->_ndimension;d++){
-	    in_oidx += in_grid->_ostride[d] * ( lcorr[d] % in_grid->_rdimensions[d] );
-	    in_lane += in_grid->_istride[d] * ( lcorr[d] / in_grid->_rdimensions[d] );
-	  }
-	  fmap_host[out_lane + Nsimd_out*out_oidx] = std::pair<Integer,Integer>( in_oidx, in_lane );
-	}
-      });
-
-    //Copy the map to the device (if we had a way to tell if an accelerator is in use we could avoid this copy for CPU-only machines)
-    size_t fmap_bytes = out_grid->lSites() * sizeof(std::pair<Integer,Integer>);
-    fmap_device = (std::pair<Integer,Integer>*)acceleratorAllocDevice(fmap_bytes);
-    acceleratorCopyToDevice(fmap_host.data(), fmap_device, fmap_bytes); 
-  }
-
-  //Prevent moving or copying
-  precisionChangeWorkspace(const precisionChangeWorkspace &r) = delete;
-  precisionChangeWorkspace(precisionChangeWorkspace &&r) = delete;
-  precisionChangeWorkspace &operator=(const precisionChangeWorkspace &r) = delete;
-  precisionChangeWorkspace &operator=(precisionChangeWorkspace &&r) = delete;
-  
-  std::pair<Integer,Integer> const* getMap() const{ return fmap_device; }
-
-  ~precisionChangeWorkspace(){
-    acceleratorFreeDevice(fmap_device);
-  }
-};
-
-
-//Convert a lattice of one precision to another. The input workspace contains the mapping data.
-template<class VobjOut, class VobjIn>
-void precisionChange(Lattice<VobjOut> &out, const Lattice<VobjIn> &in, const precisionChangeWorkspace &workspace){
-  static_assert( std::is_same<typename VobjOut::DoublePrecision, typename VobjIn::DoublePrecision>::value == 1, "copyLane: tensor types must be the same" ); //if tensor types are same the DoublePrecision type must be the same
-
-  out.Checkerboard() = in.Checkerboard();
-  constexpr int Nsimd_out = VobjOut::Nsimd();
-
-  std::pair<Integer,Integer> const* fmap_device = workspace.getMap();
-
-  //Do the copy/precision change
-  autoView( out_v , out, AcceleratorWrite);
-  autoView( in_v , in, AcceleratorRead);
-
-  accelerator_for(out_oidx, out.Grid()->oSites(), 1,{
-      std::pair<Integer,Integer> const* fmap_osite = fmap_device + out_oidx*Nsimd_out;
-      for(int out_lane=0; out_lane < Nsimd_out; out_lane++){      
-	int in_oidx = fmap_osite[out_lane].first;
-	int in_lane = fmap_osite[out_lane].second;
-	copyLane(out_v[out_oidx], out_lane, in_v[in_oidx], in_lane);
-      }
-    });
-}
-
 //Convert a Lattice from one precision to another
-//Generate the workspace in place; if multiple calls with the same mapping are performed, consider pregenerating the workspace and reusing
 template<class VobjOut, class VobjIn>
-void precisionChange(Lattice<VobjOut> &out, const Lattice<VobjIn> &in){
-  precisionChangeWorkspace workspace(out.Grid(), in.Grid());
-  precisionChange(out, in, workspace);
-}
+void precisionChange(Lattice<VobjOut> &out, const Lattice<VobjIn> &in)
+{
+  assert(out.Grid()->Nd() == in.Grid()->Nd());
+  for(int d=0;d<out.Grid()->Nd();d++){
+    assert(out.Grid()->FullDimensions()[d] == in.Grid()->FullDimensions()[d]);
+  }
+  out.Checkerboard() = in.Checkerboard();
+  GridBase *in_grid=in.Grid();
+  GridBase *out_grid = out.Grid();

+  typedef typename VobjOut::scalar_object SobjOut;
+  typedef typename VobjIn::scalar_object SobjIn;
+
+  int ndim = out.Grid()->Nd();
+  int out_nsimd = out_grid->Nsimd();
+    
+  std::vector<Coordinate > out_icoor(out_nsimd);
+      
+  for(int lane=0; lane < out_nsimd; lane++){
+    out_icoor[lane].resize(ndim);
+    out_grid->iCoorFromIindex(out_icoor[lane], lane);
+  }
+        
+  std::vector<SobjOut> in_slex_conv(in_grid->lSites());
+  unvectorizeToLexOrdArray(in_slex_conv, in);
+    
+  autoView( out_v , out, CpuWrite);
+  thread_for(out_oidx,out_grid->oSites(),{
+    Coordinate out_ocoor(ndim);
+    out_grid->oCoorFromOindex(out_ocoor, out_oidx);
+
+    ExtractPointerArray<SobjOut> ptrs(out_nsimd);      
+
+    Coordinate lcoor(out_grid->Nd());
+      
+    for(int lane=0; lane < out_nsimd; lane++){
+      for(int mu=0;mu<ndim;mu++)
+	lcoor[mu] = out_ocoor[mu] + out_grid->_rdimensions[mu]*out_icoor[lane][mu];
+	
+      int llex; Lexicographic::IndexFromCoor(lcoor, llex, out_grid->_ldimensions);
+      ptrs[lane] = &in_slex_conv[llex];
+    }
+    merge(out_v[out_oidx], ptrs, 0);
+  });
+}

 ////////////////////////////////////////////////////////////////////////////////
 // Communicate between grids
@@ -69,7 +69,6 @@ GridLogger GridLogDebug  (1, "Debug", GridLogColours, "PURPLE");
 GridLogger GridLogPerformance(1, "Performance", GridLogColours, "GREEN");
 GridLogger GridLogIterative  (1, "Iterative", GridLogColours, "BLUE");
 GridLogger GridLogIntegrator (1, "Integrator", GridLogColours, "BLUE");
-GridLogger GridLogHMC (1, "HMC", GridLogColours, "BLUE");

 void GridLogConfigure(std::vector<std::string> &logstreams) {
  GridLogError.Active(0);
@@ -80,7 +79,6 @@ void GridLogConfigure(std::vector<std::string> &logstreams) {
  GridLogPerformance.Active(0);
  GridLogIntegrator.Active(1);
  GridLogColours.Active(0);
-  GridLogHMC.Active(1);

  for (int i = 0; i < logstreams.size(); i++) {
    if (logstreams[i] == std::string("Error"))       GridLogError.Active(1);
@@ -89,8 +87,7 @@ void GridLogConfigure(std::vector<std::string> &logstreams) {
    if (logstreams[i] == std::string("Iterative"))   GridLogIterative.Active(1);
    if (logstreams[i] == std::string("Debug"))       GridLogDebug.Active(1);
    if (logstreams[i] == std::string("Performance")) GridLogPerformance.Active(1);
-    if (logstreams[i] == std::string("NoIntegrator"))  GridLogIntegrator.Active(0);
-    if (logstreams[i] == std::string("NoHMC"))         GridLogHMC.Active(0);
+    if (logstreams[i] == std::string("Integrator"))  GridLogIntegrator.Active(1);
    if (logstreams[i] == std::string("Colours"))     GridLogColours.Active(1);
  }
 }
@@ -182,7 +182,6 @@ extern GridLogger GridLogDebug  ;
 extern GridLogger GridLogPerformance;
 extern GridLogger GridLogIterative  ;
 extern GridLogger GridLogIntegrator  ;
-extern GridLogger GridLogHMC;
 extern Colours    GridLogColours;

 std::string demangle(const char* name) ;
@@ -271,7 +271,7 @@ class BinaryIO {
 			      uint32_t &scidac_csumb)
  {
    grid->Barrier();
-    GridStopWatch timer; 
+    GridStopWatch timer, insideTimer; 
    GridStopWatch bstimer;
    
    nersc_csum=0;
@@ -363,7 +363,10 @@ class BinaryIO {
 	std::cout<< GridLogMessage<<"IOobject: MPI read I/O "<< file<< std::endl;
 	ierr=MPI_File_open(grid->communicator,(char *) file.c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);    assert(ierr==0);
 	ierr=MPI_File_set_view(fh, disp, mpiObject, fileArray, "native", MPI_INFO_NULL);    assert(ierr==0);
-	ierr=MPI_File_read_all(fh, &iodata[0], 1, localArray, &status);    assert(ierr==0);
+  insideTimer.Start();
+	ierr=MPI_File_read_all(fh, &iodata[0], 1, localArray, &status);
+  insideTimer.Stop();
+  assert(ierr==0);
 	MPI_File_close(&fh);
 	MPI_Type_free(&fileArray);
 	MPI_Type_free(&localArray);
@@ -438,7 +441,9 @@ class BinaryIO {
        assert(ierr == 0);

        std::cout << GridLogDebug << "MPI write I/O write all " << file << std::endl;
+        insideTimer.Start();
        ierr = MPI_File_write_all(fh, &iodata[0], 1, localArray, &status);
+        insideTimer.Stop();
        assert(ierr == 0);

        MPI_Offset os;
@@ -516,8 +521,13 @@ class BinaryIO {
    if ( control & BINARYIO_READ) std::cout << " read  ";
    else                          std::cout << " write ";
    uint64_t bytes = sizeof(fobj)*iodata.size()*nrank;
-    std::cout<< lastPerf.size <<" bytes in "<< timer.Elapsed() <<" "
+    std::cout<< lastPerf.size <<"bytes in "<< timer.Elapsed() <<" "
 	     << lastPerf.mbytesPerSecond <<" MB/s "<<std::endl;
+    std::cout << GridLogMessage << "IOobject: pure MPI IO call " 
+              << lastPerf.size <<" bytes in " 
+              << insideTimer.Elapsed() << " "
+              << lastPerf.size/1024./1024./(insideTimer.useconds()/1.0e6)
+              <<" MB/s "<<std::endl;

    std::cout<<GridLogMessage<<"IOobject: endian and checksum overhead "<<bstimer.Elapsed()  <<std::endl;

@@ -39,10 +39,8 @@ using namespace Grid;
 ////////////////////////////////////////////////////////////////////////////////
 class NerscIO : public BinaryIO { 
 public:
-  typedef Lattice<vLorentzColourMatrixD> GaugeField;

-  // Enable/disable exiting if the plaquette in the header does not match the value computed (default true)
-  static bool & exitOnReadPlaquetteMismatch(){ static bool v=true; return v; }
+  typedef Lattice<vLorentzColourMatrixD> GaugeField;

  static inline void truncate(std::string file){
    std::ofstream fout(file,std::ios::out);
@@ -200,7 +198,7 @@ public:
      std::cerr << " nersc_csum  " <<std::hex<< nersc_csum << " " << header.checksum<< std::dec<< std::endl;
      exit(0);
    }
-    if(exitOnReadPlaquetteMismatch()) assert(fabs(clone.plaquette -header.plaquette ) < 1.0e-5 );
+    assert(fabs(clone.plaquette -header.plaquette ) < 1.0e-5 );
    assert(fabs(clone.link_trace-header.link_trace) < 1.0e-6 );
    assert(nersc_csum == header.checksum );
      
@@ -63,7 +63,6 @@ static constexpr int Ngp=2; // gparity index range
 #define ColourIndex  (2)
 #define SpinIndex    (1)
 #define LorentzIndex (0)
-#define GparityFlavourIndex (0)

 // Also should make these a named enum type
 static constexpr int DaggerNo=0;
@@ -88,8 +87,6 @@ template<typename T> struct isCoarsened {
 template <typename T> using IfCoarsened    = Invoke<std::enable_if< isCoarsened<T>::value,int> > ;
 template <typename T> using IfNotCoarsened = Invoke<std::enable_if<!isCoarsened<T>::value,int> > ;

-const int GparityFlavourTensorIndex = 3; //TensorLevel counts from the bottom!
-
 // ChrisK very keen to add extra space for Gparity doubling.
 //
 // Also add domain wall index, in a way where Wilson operator 
@@ -113,10 +110,8 @@ template<typename vtype> using iHalfSpinColourVector      = iScalar<iVector<iVec
    template<typename vtype> using iSpinColourSpinColourMatrix  = iScalar<iMatrix<iMatrix<iMatrix<iMatrix<vtype, Nc>, Ns>, Nc>, Ns> >;


-template<typename vtype> using iGparityFlavourVector                = iVector<iScalar<iScalar<vtype> >, Ngp>;
 template<typename vtype> using iGparitySpinColourVector       = iVector<iVector<iVector<vtype, Nc>, Ns>, Ngp >;
 template<typename vtype> using iGparityHalfSpinColourVector   = iVector<iVector<iVector<vtype, Nc>, Nhs>, Ngp >;
-template<typename vtype> using iGparityFlavourMatrix = iMatrix<iScalar<iScalar<vtype> >, Ngp>;

 // Spin matrix
 typedef iSpinMatrix<Complex  >          SpinMatrix;
@@ -181,16 +176,6 @@ typedef iDoubleStoredColourMatrix<vComplex > vDoubleStoredColourMatrix;
 typedef iDoubleStoredColourMatrix<vComplexF> vDoubleStoredColourMatrixF;
 typedef iDoubleStoredColourMatrix<vComplexD> vDoubleStoredColourMatrixD;

-//G-parity flavour matrix
-typedef iGparityFlavourMatrix<Complex> GparityFlavourMatrix;
-typedef iGparityFlavourMatrix<ComplexF> GparityFlavourMatrixF;
-typedef iGparityFlavourMatrix<ComplexD> GparityFlavourMatrixD;
-
-typedef iGparityFlavourMatrix<vComplex> vGparityFlavourMatrix;
-typedef iGparityFlavourMatrix<vComplexF> vGparityFlavourMatrixF;
-typedef iGparityFlavourMatrix<vComplexD> vGparityFlavourMatrixD;
-
-
 // Spin vector
 typedef iSpinVector<Complex >           SpinVector;
 typedef iSpinVector<ComplexF>           SpinVectorF;
@@ -235,16 +220,6 @@ typedef iHalfSpinColourVector<ComplexD> HalfSpinColourVectorD;
 typedef iHalfSpinColourVector<vComplex > vHalfSpinColourVector;
 typedef iHalfSpinColourVector<vComplexF> vHalfSpinColourVectorF;
 typedef iHalfSpinColourVector<vComplexD> vHalfSpinColourVectorD;
-
-//G-parity flavour vector
-typedef iGparityFlavourVector<Complex >         GparityFlavourVector;
-typedef iGparityFlavourVector<ComplexF>         GparityFlavourVectorF;
-typedef iGparityFlavourVector<ComplexD>         GparityFlavourVectorD;
-
-typedef iGparityFlavourVector<vComplex >         vGparityFlavourVector;
-typedef iGparityFlavourVector<vComplexF>         vGparityFlavourVectorF;
-typedef iGparityFlavourVector<vComplexD>         vGparityFlavourVectorD;
-
    
 // singlets
 typedef iSinglet<Complex >         TComplex;     // FIXME This is painful. Tensor singlet complex type.
@@ -36,8 +36,7 @@ NAMESPACE_BEGIN(Grid);

 // These can move into a params header and be given MacroMagic serialisation
 struct GparityWilsonImplParams {
-  Coordinate twists; //Here the first Nd-1 directions are treated as "spatial", and a twist value of 1 indicates G-parity BCs in that direction. 
-                     //mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs
+  Coordinate twists;
  GparityWilsonImplParams() : twists(Nd, 0) {};
 };
  
@@ -66,8 +65,7 @@ struct StaggeredImplParams {
 				    RealD, tolerance, 
 				    int,   degree, 
 				    int,   precision,
-				    int,   BoundsCheckFreq,
-				    RealD, BoundsCheckTol);
+				    int,   BoundsCheckFreq);
    
  // MaxIter and tolerance, vectors??
    
@@ -78,61 +76,15 @@ struct StaggeredImplParams {
 				RealD tol      = 1.0e-8, 
                           	int _degree    = 10,
 				int _precision = 64,
-				int _BoundsCheckFreq=20,
-				double _BoundsCheckTol=1e-6)
+				int _BoundsCheckFreq=20)
      : lo(_lo),
 	hi(_hi),
 	MaxIter(_maxit),
 	tolerance(tol),
 	degree(_degree),
        precision(_precision),
-        BoundsCheckFreq(_BoundsCheckFreq),
-        BoundsCheckTol(_BoundsCheckTol){};
+        BoundsCheckFreq(_BoundsCheckFreq){};
  };
-
-
-  /*Action parameters for the generalized rational action
-    The approximation is for (M^dag M)^{1/inv_pow}
-    where inv_pow is the denominator of the fractional power.
-    Default inv_pow=2 for square root, making this equivalent to 
-    the OneFlavourRational action
-  */
-    struct RationalActionParams : Serializable {
-    GRID_SERIALIZABLE_CLASS_MEMBERS(RationalActionParams, 
-				    int, inv_pow, 
-				    RealD, lo, //low eigenvalue bound of rational approx
-				    RealD, hi, //high eigenvalue bound of rational approx
-				    int,   MaxIter,  //maximum iterations in msCG
-				    RealD, action_tolerance,  //msCG tolerance in action evaluation
-				    int,   action_degree, //rational approx tolerance in action evaluation
-				    RealD, md_tolerance,  //msCG tolerance in MD integration
-				    int,   md_degree, //rational approx tolerance in MD integration
-				    int,   precision, //precision of floating point arithmetic
-				    int,   BoundsCheckFreq); //frequency the approximation is tested (with Metropolis degree/tolerance); 0 disables the check
-  // constructor 
-  RationalActionParams(int _inv_pow = 2,
-		       RealD _lo      = 0.0, 
-		       RealD _hi      = 1.0, 
-		       int _maxit     = 1000,
-		       RealD _action_tolerance      = 1.0e-8, 
-		       int _action_degree    = 10,
-		       RealD _md_tolerance      = 1.0e-8, 
-		       int _md_degree    = 10,
-		       int _precision = 64,
-		       int _BoundsCheckFreq=20)
-    : inv_pow(_inv_pow), 
-      lo(_lo),
-      hi(_hi),
-      MaxIter(_maxit),
-      action_tolerance(_action_tolerance),
-      action_degree(_action_degree),
-      md_tolerance(_md_tolerance),
-      md_degree(_md_degree),
-      precision(_precision),
-      BoundsCheckFreq(_BoundsCheckFreq){};
-  };
-
-
  
 NAMESPACE_END(Grid);

@@ -115,9 +115,9 @@ typedef WilsonFermion<WilsonImplR> WilsonFermionR;
 typedef WilsonFermion<WilsonImplF> WilsonFermionF;
 typedef WilsonFermion<WilsonImplD> WilsonFermionD;

-//typedef WilsonFermion<WilsonImplRL> WilsonFermionRL;
-//typedef WilsonFermion<WilsonImplFH> WilsonFermionFH;
-//typedef WilsonFermion<WilsonImplDF> WilsonFermionDF;
+typedef WilsonFermion<WilsonImplRL> WilsonFermionRL;
+typedef WilsonFermion<WilsonImplFH> WilsonFermionFH;
+typedef WilsonFermion<WilsonImplDF> WilsonFermionDF;

 typedef WilsonFermion<WilsonAdjImplR> WilsonAdjFermionR;
 typedef WilsonFermion<WilsonAdjImplF> WilsonAdjFermionF;
@@ -158,41 +158,41 @@ typedef DomainWallFermion<WilsonImplR> DomainWallFermionR;
 typedef DomainWallFermion<WilsonImplF> DomainWallFermionF;
 typedef DomainWallFermion<WilsonImplD> DomainWallFermionD;

-//typedef DomainWallFermion<WilsonImplRL> DomainWallFermionRL;
-//typedef DomainWallFermion<WilsonImplFH> DomainWallFermionFH;
-//typedef DomainWallFermion<WilsonImplDF> DomainWallFermionDF;
+typedef DomainWallFermion<WilsonImplRL> DomainWallFermionRL;
+typedef DomainWallFermion<WilsonImplFH> DomainWallFermionFH;
+typedef DomainWallFermion<WilsonImplDF> DomainWallFermionDF;

 typedef DomainWallEOFAFermion<WilsonImplR> DomainWallEOFAFermionR;
 typedef DomainWallEOFAFermion<WilsonImplF> DomainWallEOFAFermionF;
 typedef DomainWallEOFAFermion<WilsonImplD> DomainWallEOFAFermionD;

-//typedef DomainWallEOFAFermion<WilsonImplRL> DomainWallEOFAFermionRL;
-//typedef DomainWallEOFAFermion<WilsonImplFH> DomainWallEOFAFermionFH;
-//typedef DomainWallEOFAFermion<WilsonImplDF> DomainWallEOFAFermionDF;
+typedef DomainWallEOFAFermion<WilsonImplRL> DomainWallEOFAFermionRL;
+typedef DomainWallEOFAFermion<WilsonImplFH> DomainWallEOFAFermionFH;
+typedef DomainWallEOFAFermion<WilsonImplDF> DomainWallEOFAFermionDF;

 typedef MobiusFermion<WilsonImplR> MobiusFermionR;
 typedef MobiusFermion<WilsonImplF> MobiusFermionF;
 typedef MobiusFermion<WilsonImplD> MobiusFermionD;

-//typedef MobiusFermion<WilsonImplRL> MobiusFermionRL;
-//typedef MobiusFermion<WilsonImplFH> MobiusFermionFH;
-//typedef MobiusFermion<WilsonImplDF> MobiusFermionDF;
+typedef MobiusFermion<WilsonImplRL> MobiusFermionRL;
+typedef MobiusFermion<WilsonImplFH> MobiusFermionFH;
+typedef MobiusFermion<WilsonImplDF> MobiusFermionDF;

 typedef MobiusEOFAFermion<WilsonImplR> MobiusEOFAFermionR;
 typedef MobiusEOFAFermion<WilsonImplF> MobiusEOFAFermionF;
 typedef MobiusEOFAFermion<WilsonImplD> MobiusEOFAFermionD;

-//typedef MobiusEOFAFermion<WilsonImplRL> MobiusEOFAFermionRL;
-//typedef MobiusEOFAFermion<WilsonImplFH> MobiusEOFAFermionFH;
-//typedef MobiusEOFAFermion<WilsonImplDF> MobiusEOFAFermionDF;
+typedef MobiusEOFAFermion<WilsonImplRL> MobiusEOFAFermionRL;
+typedef MobiusEOFAFermion<WilsonImplFH> MobiusEOFAFermionFH;
+typedef MobiusEOFAFermion<WilsonImplDF> MobiusEOFAFermionDF;

 typedef ZMobiusFermion<ZWilsonImplR> ZMobiusFermionR;
 typedef ZMobiusFermion<ZWilsonImplF> ZMobiusFermionF;
 typedef ZMobiusFermion<ZWilsonImplD> ZMobiusFermionD;

-//typedef ZMobiusFermion<ZWilsonImplRL> ZMobiusFermionRL;
-//typedef ZMobiusFermion<ZWilsonImplFH> ZMobiusFermionFH;
-//typedef ZMobiusFermion<ZWilsonImplDF> ZMobiusFermionDF;
+typedef ZMobiusFermion<ZWilsonImplRL> ZMobiusFermionRL;
+typedef ZMobiusFermion<ZWilsonImplFH> ZMobiusFermionFH;
+typedef ZMobiusFermion<ZWilsonImplDF> ZMobiusFermionDF;

 // Ls vectorised
 typedef ScaledShamirFermion<WilsonImplR> ScaledShamirFermionR;
@@ -235,49 +235,49 @@ typedef WilsonFermion<GparityWilsonImplR>     GparityWilsonFermionR;
 typedef WilsonFermion<GparityWilsonImplF>     GparityWilsonFermionF;
 typedef WilsonFermion<GparityWilsonImplD>     GparityWilsonFermionD;

-//typedef WilsonFermion<GparityWilsonImplRL>     GparityWilsonFermionRL;
-//typedef WilsonFermion<GparityWilsonImplFH>     GparityWilsonFermionFH;
-//typedef WilsonFermion<GparityWilsonImplDF>     GparityWilsonFermionDF;
+typedef WilsonFermion<GparityWilsonImplRL>     GparityWilsonFermionRL;
+typedef WilsonFermion<GparityWilsonImplFH>     GparityWilsonFermionFH;
+typedef WilsonFermion<GparityWilsonImplDF>     GparityWilsonFermionDF;

 typedef DomainWallFermion<GparityWilsonImplR> GparityDomainWallFermionR;
 typedef DomainWallFermion<GparityWilsonImplF> GparityDomainWallFermionF;
 typedef DomainWallFermion<GparityWilsonImplD> GparityDomainWallFermionD;

-//typedef DomainWallFermion<GparityWilsonImplRL> GparityDomainWallFermionRL;
-//typedef DomainWallFermion<GparityWilsonImplFH> GparityDomainWallFermionFH;
-//typedef DomainWallFermion<GparityWilsonImplDF> GparityDomainWallFermionDF;
+typedef DomainWallFermion<GparityWilsonImplRL> GparityDomainWallFermionRL;
+typedef DomainWallFermion<GparityWilsonImplFH> GparityDomainWallFermionFH;
+typedef DomainWallFermion<GparityWilsonImplDF> GparityDomainWallFermionDF;

 typedef DomainWallEOFAFermion<GparityWilsonImplR> GparityDomainWallEOFAFermionR;
 typedef DomainWallEOFAFermion<GparityWilsonImplF> GparityDomainWallEOFAFermionF;
 typedef DomainWallEOFAFermion<GparityWilsonImplD> GparityDomainWallEOFAFermionD;

-//typedef DomainWallEOFAFermion<GparityWilsonImplRL> GparityDomainWallEOFAFermionRL;
-//typedef DomainWallEOFAFermion<GparityWilsonImplFH> GparityDomainWallEOFAFermionFH;
-//typedef DomainWallEOFAFermion<GparityWilsonImplDF> GparityDomainWallEOFAFermionDF;
+typedef DomainWallEOFAFermion<GparityWilsonImplRL> GparityDomainWallEOFAFermionRL;
+typedef DomainWallEOFAFermion<GparityWilsonImplFH> GparityDomainWallEOFAFermionFH;
+typedef DomainWallEOFAFermion<GparityWilsonImplDF> GparityDomainWallEOFAFermionDF;

 typedef WilsonTMFermion<GparityWilsonImplR> GparityWilsonTMFermionR;
 typedef WilsonTMFermion<GparityWilsonImplF> GparityWilsonTMFermionF;
 typedef WilsonTMFermion<GparityWilsonImplD> GparityWilsonTMFermionD;

-//typedef WilsonTMFermion<GparityWilsonImplRL> GparityWilsonTMFermionRL;
-//typedef WilsonTMFermion<GparityWilsonImplFH> GparityWilsonTMFermionFH;
-//typedef WilsonTMFermion<GparityWilsonImplDF> GparityWilsonTMFermionDF;
+typedef WilsonTMFermion<GparityWilsonImplRL> GparityWilsonTMFermionRL;
+typedef WilsonTMFermion<GparityWilsonImplFH> GparityWilsonTMFermionFH;
+typedef WilsonTMFermion<GparityWilsonImplDF> GparityWilsonTMFermionDF;

 typedef MobiusFermion<GparityWilsonImplR> GparityMobiusFermionR;
 typedef MobiusFermion<GparityWilsonImplF> GparityMobiusFermionF;
 typedef MobiusFermion<GparityWilsonImplD> GparityMobiusFermionD;

-//typedef MobiusFermion<GparityWilsonImplRL> GparityMobiusFermionRL;
-//typedef MobiusFermion<GparityWilsonImplFH> GparityMobiusFermionFH;
-//typedef MobiusFermion<GparityWilsonImplDF> GparityMobiusFermionDF;
+typedef MobiusFermion<GparityWilsonImplRL> GparityMobiusFermionRL;
+typedef MobiusFermion<GparityWilsonImplFH> GparityMobiusFermionFH;
+typedef MobiusFermion<GparityWilsonImplDF> GparityMobiusFermionDF;

 typedef MobiusEOFAFermion<GparityWilsonImplR> GparityMobiusEOFAFermionR;
 typedef MobiusEOFAFermion<GparityWilsonImplF> GparityMobiusEOFAFermionF;
 typedef MobiusEOFAFermion<GparityWilsonImplD> GparityMobiusEOFAFermionD;

-//typedef MobiusEOFAFermion<GparityWilsonImplRL> GparityMobiusEOFAFermionRL;
-//typedef MobiusEOFAFermion<GparityWilsonImplFH> GparityMobiusEOFAFermionFH;
-//typedef MobiusEOFAFermion<GparityWilsonImplDF> GparityMobiusEOFAFermionDF;
+typedef MobiusEOFAFermion<GparityWilsonImplRL> GparityMobiusEOFAFermionRL;
+typedef MobiusEOFAFermion<GparityWilsonImplFH> GparityMobiusEOFAFermionFH;
+typedef MobiusEOFAFermion<GparityWilsonImplDF> GparityMobiusEOFAFermionDF;

 typedef ImprovedStaggeredFermion<StaggeredImplR> ImprovedStaggeredFermionR;
 typedef ImprovedStaggeredFermion<StaggeredImplF> ImprovedStaggeredFermionF;
@@ -30,18 +30,6 @@ directory

 NAMESPACE_BEGIN(Grid);

-/*
-  Policy implementation for G-parity boundary conditions
-
-  Rather than treating the gauge field as a flavored field, the Grid implementation of G-parity treats the gauge field as a regular
-  field with complex conjugate boundary conditions. In order to ensure the second flavor interacts with the conjugate links and the first
-  with the regular links we overload the functionality of doubleStore, whose purpose is to store the gauge field and the barrel-shifted gauge field
-  to avoid communicating links when applying the Dirac operator, such that the double-stored field contains also a flavor index which maps to
-  either the link or the conjugate link. This flavored field is then used by multLink to apply the correct link to a spinor.
-
-  Here the first Nd-1 directions are treated as "spatial", and a twist value of 1 indicates G-parity BCs in that direction. 
-  mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs
- */
 template <class S, class Representation = FundamentalRepresentation, class Options=CoeffReal>
 class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Representation::Dimension> > {
 public:
@@ -125,7 +113,7 @@ public:
    || ((distance== 1)&&(icoor[direction]==1))
    || ((distance==-1)&&(icoor[direction]==0));

-    permute_lane = permute_lane && SE->_around_the_world && St.parameters.twists[mmu] && mmu < Nd-1; //only if we are going around the world in a spatial direction
+    permute_lane = permute_lane && SE->_around_the_world && St.parameters.twists[mmu]; //only if we are going around the world

    //Apply the links
    int f_upper = permute_lane ? 1 : 0;
@@ -151,10 +139,10 @@ public:
    assert((distance == 1) || (distance == -1));  // nearest neighbour stencil hard code
    assert((sl == 1) || (sl == 2));

-    //If this site is an global boundary site, perform the G-parity flavor twist
-    if ( mmu < Nd-1 && SE->_around_the_world && St.parameters.twists[mmu] ) {
+    if ( SE->_around_the_world && St.parameters.twists[mmu] ) {
+
      if ( sl == 2 ) {
-	//Only do the twist for lanes on the edge of the physical node
+       
 	ExtractBuffer<sobj> vals(Nsimd);

 	extract(chi,vals);
@@ -209,19 +197,6 @@ public:
    reg = memory;
  }

-
-  //Poke 'poke_f0' onto flavor 0 and 'poke_f1' onto flavor 1 in direction mu of the doubled gauge field Uds
-  inline void pokeGparityDoubledGaugeField(DoubledGaugeField &Uds, const GaugeLinkField &poke_f0, const GaugeLinkField &poke_f1, const int mu){
-    autoView(poke_f0_v, poke_f0, CpuRead);
-    autoView(poke_f1_v, poke_f1, CpuRead);
-    autoView(Uds_v, Uds, CpuWrite);
-    thread_foreach(ss,poke_f0_v,{
-	Uds_v[ss](0)(mu) = poke_f0_v[ss]();
-	Uds_v[ss](1)(mu) = poke_f1_v[ss]();
-      });
-  }
-    
-
  inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
  {
    conformable(Uds.Grid(),GaugeGrid);
@@ -232,19 +207,14 @@ public:
    GaugeLinkField Uconj(GaugeGrid);
   
    Lattice<iScalar<vInteger> > coor(GaugeGrid);
-
-    //Here the first Nd-1 directions are treated as "spatial", and a twist value of 1 indicates G-parity BCs in that direction. 
-    //mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs        
-    for(int mu=0;mu<Nd-1;mu++){
-
-      if( Params.twists[mu] ){
-	LatticeCoordinate(coor,mu);
-      }
+        
+    for(int mu=0;mu<Nd;mu++){
+          
+      LatticeCoordinate(coor,mu);
          
      U     = PeekIndex<LorentzIndex>(Umu,mu);
      Uconj = conjugate(U);
     
-      // Implement the isospin rotation sign on the boundary between f=1 and f=0
      // This phase could come from a simple bc 1,1,-1,1 ..
      int neglink = GaugeGrid->GlobalDimensions()[mu]-1;
      if ( Params.twists[mu] ) { 
@@ -259,7 +229,7 @@ public:
 	thread_foreach(ss,U_v,{
 	    Uds_v[ss](0)(mu) = U_v[ss]();
 	    Uds_v[ss](1)(mu) = Uconj_v[ss]();
-	});
+	  });
      }
          
      U     = adj(Cshift(U    ,mu,-1));      // correct except for spanning the boundary
@@ -290,38 +260,6 @@ public:
        });
      }
    }
-
-    { //periodic / antiperiodic temporal BCs
-      int mu = Nd-1;
-      int L   = GaugeGrid->GlobalDimensions()[mu];
-      int Lmu = L - 1;
-
-      LatticeCoordinate(coor, mu);
-
-      U = PeekIndex<LorentzIndex>(Umu, mu); //Get t-directed links
-      
-      GaugeLinkField *Upoke = &U;
-
-      if(Params.twists[mu]){ //antiperiodic
-	Utmp =  where(coor == Lmu, -U, U);
-	Upoke = &Utmp;
-      }
-    
-      Uconj = conjugate(*Upoke); //second flavor interacts with conjugate links      
-      pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu);
-
-      //Get the barrel-shifted field
-      Utmp = adj(Cshift(U, mu, -1)); //is a forward shift!
-      Upoke = &Utmp;
-
-      if(Params.twists[mu]){
-	U = where(coor == 0, -Utmp, Utmp);  //boundary phase
-	Upoke = &U;
-      }
-      
-      Uconj = conjugate(*Upoke);
-      pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu + 4);
-    }
  }
      
  inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A, int mu) {
@@ -360,48 +298,28 @@ public:
  inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
    assert(0);
  }
- 
-  inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde, int mu) {
-    int Ls=Btilde.Grid()->_fdimensions[0];
-    
-    {
-      GridBase *GaugeGrid = mat.Grid();
-      Lattice<iScalar<vInteger> > coor(GaugeGrid);
-
-      if( Params.twists[mu] ){
-	LatticeCoordinate(coor,mu);
-      }
-
-      autoView( mat_v , mat, AcceleratorWrite);
-      autoView( Btilde_v , Btilde, AcceleratorRead);
-      autoView( Atilde_v , Atilde, AcceleratorRead);
-      accelerator_for(sss,mat.Grid()->oSites(), FermionField::vector_type::Nsimd(),{	  
-  	  int sU=sss;
-  	  typedef decltype(coalescedRead(mat_v[sU](mu)() )) ColorMatrixType;
-  	  ColorMatrixType sum;
-  	  zeroit(sum);
-  	  for(int s=0;s<Ls;s++){
-  	    int sF = s+Ls*sU;
-  	    for(int spn=0;spn<Ns;spn++){ //sum over spin
-	      //Flavor 0
-  	      auto bb = coalescedRead(Btilde_v[sF](0)(spn) ); //color vector
-  	      auto aa = coalescedRead(Atilde_v[sF](0)(spn) );
-  	      sum = sum + outerProduct(bb,aa);
-
-  	      //Flavor 1
-  	      bb = coalescedRead(Btilde_v[sF](1)(spn) );
-  	      aa = coalescedRead(Atilde_v[sF](1)(spn) );
-  	      sum = sum + conjugate(outerProduct(bb,aa));
-  	    }
-  	  }	    
-  	  coalescedWrite(mat_v[sU](mu)(), sum);
-  	});
-    }
-  }
-
-
  
+  inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde, int mu) {

+    int Ls = Btilde.Grid()->_fdimensions[0];
+        
+    GaugeLinkField tmp(mat.Grid());
+    tmp = Zero();
+    {
+      autoView( tmp_v , tmp, CpuWrite);
+      autoView( Atilde_v , Atilde, CpuRead);
+      autoView( Btilde_v , Btilde, CpuRead);
+      thread_for(ss,tmp.Grid()->oSites(),{
+	  for (int s = 0; s < Ls; s++) {
+	    int sF = s + Ls * ss;
+	    auto ttmp = traceIndex<SpinIndex>(outerProduct(Btilde_v[sF], Atilde_v[sF]));
+	    tmp_v[ss]() = tmp_v[ss]() + ttmp(0, 0) + conjugate(ttmp(1, 1));
+	  }
+	});
+    }
+    PokeIndex<LorentzIndex>(mat, tmp, mu);
+    return;
+  }
  
 };

@@ -409,8 +327,8 @@ typedef GparityWilsonImpl<vComplex , FundamentalRepresentation,CoeffReal> Gparit
 typedef GparityWilsonImpl<vComplexF, FundamentalRepresentation,CoeffReal> GparityWilsonImplF;  // Float
 typedef GparityWilsonImpl<vComplexD, FundamentalRepresentation,CoeffReal> GparityWilsonImplD;  // Double
 
-//typedef GparityWilsonImpl<vComplex , FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplRL;  // Real.. whichever prec
-//typedef GparityWilsonImpl<vComplexF, FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplFH;  // Float
-//typedef GparityWilsonImpl<vComplexD, FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplDF;  // Double
+typedef GparityWilsonImpl<vComplex , FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplRL;  // Real.. whichever prec
+typedef GparityWilsonImpl<vComplexF, FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplFH;  // Float
+typedef GparityWilsonImpl<vComplexD, FundamentalRepresentation,CoeffRealHalfComms> GparityWilsonImplDF;  // Double

 NAMESPACE_END(Grid);
@@ -68,12 +68,11 @@ public:
  /*****************************************************/
  /* Compress includes precision change if mpi data is not same */
  /*****************************************************/
-  accelerator_inline void Compress(SiteHalfSpinor &buf,const SiteSpinor &in) const {
-    typedef decltype(coalescedRead(buf)) sobj;
-    sobj sp;
-    auto sin = coalescedRead(in);
-    projector::Proj(sp,sin,mu,dag);
-    coalescedWrite(buf,sp);
+  template<class _SiteHalfSpinor, class _SiteSpinor>
+  accelerator_inline void Compress(_SiteHalfSpinor *buf,Integer o,const _SiteSpinor &in) const {
+    _SiteHalfSpinor tmp;
+    projector::Proj(tmp,in,mu,dag);
+    vstream(buf[o],tmp);
  }

  /*****************************************************/
@@ -83,18 +82,13 @@ public:
 				   const SiteHalfSpinor * __restrict__ vp0,
 				   const SiteHalfSpinor * __restrict__ vp1,
 				   Integer type,Integer o) const {
-#ifdef GRID_SIMT
-    exchangeSIMT(mp[2*o],mp[2*o+1],vp0[o],vp1[o],type);
-#else
    SiteHalfSpinor tmp1;
    SiteHalfSpinor tmp2;
    exchange(tmp1,tmp2,vp0[o],vp1[o],type);
    vstream(mp[2*o  ],tmp1);
    vstream(mp[2*o+1],tmp2);
-#endif
  }

-
  /*****************************************************/
  /* Have a decompression step if mpi data is not same */
  /*****************************************************/
@@ -111,28 +105,6 @@ public:
 					   const SiteSpinor * __restrict__ in,
 					   Integer j,Integer k, Integer m,Integer type) const
  {
-#ifdef GRID_SIMT
-    typedef SiteSpinor vobj;
-    typedef SiteHalfSpinor hvobj;
-    typedef decltype(coalescedRead(*in))    sobj;
-    typedef decltype(coalescedRead(*out0)) hsobj;
-
-    unsigned int Nsimd = vobj::Nsimd();
-    unsigned int mask = Nsimd >> (type + 1);
-    int lane = acceleratorSIMTlane(Nsimd);
-    int j0 = lane &(~mask); // inner coor zero
-    int j1 = lane |(mask) ; // inner coor one
-    const vobj *vp0 = &in[k];
-    const vobj *vp1 = &in[m];
-    const vobj *vp = (lane&mask) ? vp1:vp0;
-    auto sa = coalescedRead(*vp,j0);
-    auto sb = coalescedRead(*vp,j1);
-    hsobj psa, psb;
-    projector::Proj(psa,sa,mu,dag);
-    projector::Proj(psb,sb,mu,dag);
-    coalescedWrite(out0[j],psa);
-    coalescedWrite(out1[j],psb);
-#else
    SiteHalfSpinor temp1, temp2;
    SiteHalfSpinor temp3, temp4;
    projector::Proj(temp1,in[k],mu,dag);
@@ -140,7 +112,6 @@ public:
    exchange(temp3,temp4,temp1,temp2,type);
    vstream(out0[j],temp3);
    vstream(out1[j],temp4);
-#endif
  }

  /*****************************************************/
@@ -150,7 +121,6 @@ public:

 };

-#if 0
 template<class _HCspinor,class _Hspinor,class _Spinor, class projector>
 class WilsonCompressorTemplate< _HCspinor, _Hspinor, _Spinor, projector,
 				typename std::enable_if<!std::is_same<_HCspinor,_Hspinor>::value>::type >
@@ -179,23 +149,13 @@ public:
  /*****************************************************/
  /* Compress includes precision change if mpi data is not same */
  /*****************************************************/
-  accelerator_inline void Compress(SiteHalfSpinor &buf,const SiteSpinor &in) const {
-    SiteHalfSpinor hsp;
+  template<class _SiteHalfSpinor, class _SiteSpinor>
+  accelerator_inline void Compress(_SiteHalfSpinor *buf,Integer o,const _SiteSpinor &in) const {
+    _SiteHalfSpinor hsp;
    SiteHalfCommSpinor *hbuf = (SiteHalfCommSpinor *)buf;
    projector::Proj(hsp,in,mu,dag);
    precisionChange((vComplexLow *)&hbuf[o],(vComplexHigh *)&hsp,Nw);
  }
-  accelerator_inline void Compress(SiteHalfSpinor &buf,const SiteSpinor &in) const {
-#ifdef GRID_SIMT
-    typedef decltype(coalescedRead(buf)) sobj;
-    sobj sp;
-    auto sin = coalescedRead(in);
-    projector::Proj(sp,sin,mu,dag);
-    coalescedWrite(buf,sp);
-#else
-    projector::Proj(buf,in,mu,dag);
-#endif
-  }

  /*****************************************************/
  /* Exchange includes precision change if mpi data is not same */
@@ -243,7 +203,6 @@ public:
  accelerator_inline bool DecompressionStep(void) const { return true; }

 };
-#endif

 #define DECLARE_PROJ(Projector,Compressor,spProj)			\
  class Projector {							\
@@ -294,8 +253,33 @@ public:
  typedef typename Base::View_type View_type;
  typedef typename Base::StencilVector StencilVector;

-  void ZeroCountersi(void)  {  }
-  void Reporti(int calls)  {  }
+  double timer0;
+  double timer1;
+  double timer2;
+  double timer3;
+  double timer4;
+  double timer5;
+  double timer6;
+  uint64_t callsi;
+  void ZeroCountersi(void)
+  {
+    timer0=0;
+    timer1=0;
+    timer2=0;
+    timer3=0;
+    timer4=0;
+    timer5=0;
+    timer6=0;
+    callsi=0;
+  }
+  void Reporti(int calls)
+  {
+    if ( timer0 ) std::cout << GridLogMessage << " timer0 (HaloGatherOpt) " <<timer0/calls <<std::endl;
+    if ( timer1 ) std::cout << GridLogMessage << " timer1 (Communicate)   " <<timer1/calls <<std::endl;
+    if ( timer2 ) std::cout << GridLogMessage << " timer2 (CommsMerge )   " <<timer2/calls <<std::endl;
+    if ( timer3 ) std::cout << GridLogMessage << " timer3 (commsMergeShm) " <<timer3/calls <<std::endl;
+    if ( timer4 ) std::cout << GridLogMessage << " timer4 " <<timer4 <<std::endl;
+  }

  std::vector<int> surface_list;

@@ -337,18 +321,26 @@ public:
  {
    std::vector<std::vector<CommsRequest_t> > reqs;
    this->HaloExchangeOptGather(source,compress);
+    double t1=usecond();
    // Asynchronous MPI calls multidirectional, Isend etc...
    // Non-overlapped directions within a thread. Asynchronous calls except MPI3, threaded up to comm threads ways.
    this->Communicate();
+    double t2=usecond(); timer1 += t2-t1;
    this->CommsMerge(compress);
+    double t3=usecond(); timer2 += t3-t2;
    this->CommsMergeSHM(compress);
+    double t4=usecond(); timer3 += t4-t3;
  }
  
  template <class compressor>
  void HaloExchangeOptGather(const Lattice<vobj> &source,compressor &compress) 
  {
    this->Prepare();
+    double t0=usecond();
    this->HaloGatherOpt(source,compress);
+    double t1=usecond();
+    timer0 += t1-t0;
+    callsi++;
  }

  template <class compressor>
@@ -360,9 +352,12 @@ public:
    typedef typename compressor::SiteHalfSpinor     SiteHalfSpinor;
    typedef typename compressor::SiteHalfCommSpinor SiteHalfCommSpinor;

+    this->mpi3synctime_g-=usecond();
    this->_grid->StencilBarrier();
+    this->mpi3synctime_g+=usecond();

    assert(source.Grid()==this->_grid);
+    this->halogtime-=usecond();
    
    this->u_comm_offset=0;
      
@@ -398,6 +393,7 @@ public:
    }
    this->face_table_computed=1;
    assert(this->u_comm_offset==this->_unified_buffer_size);
+    this->halogtime+=usecond();
    accelerator_barrier();
  }

@@ -243,17 +243,17 @@ typedef WilsonImpl<vComplex,  FundamentalRepresentation, CoeffReal > WilsonImplR
 typedef WilsonImpl<vComplexF, FundamentalRepresentation, CoeffReal > WilsonImplF;  // Float
 typedef WilsonImpl<vComplexD, FundamentalRepresentation, CoeffReal > WilsonImplD;  // Double

-//typedef WilsonImpl<vComplex,  FundamentalRepresentation, CoeffRealHalfComms > WilsonImplRL;  // Real.. whichever prec
-//typedef WilsonImpl<vComplexF, FundamentalRepresentation, CoeffRealHalfComms > WilsonImplFH;  // Float
-//typedef WilsonImpl<vComplexD, FundamentalRepresentation, CoeffRealHalfComms > WilsonImplDF;  // Double
+typedef WilsonImpl<vComplex,  FundamentalRepresentation, CoeffRealHalfComms > WilsonImplRL;  // Real.. whichever prec
+typedef WilsonImpl<vComplexF, FundamentalRepresentation, CoeffRealHalfComms > WilsonImplFH;  // Float
+typedef WilsonImpl<vComplexD, FundamentalRepresentation, CoeffRealHalfComms > WilsonImplDF;  // Double

 typedef WilsonImpl<vComplex,  FundamentalRepresentation, CoeffComplex > ZWilsonImplR; // Real.. whichever prec
 typedef WilsonImpl<vComplexF, FundamentalRepresentation, CoeffComplex > ZWilsonImplF; // Float
 typedef WilsonImpl<vComplexD, FundamentalRepresentation, CoeffComplex > ZWilsonImplD; // Double

-//typedef WilsonImpl<vComplex,  FundamentalRepresentation, CoeffComplexHalfComms > ZWilsonImplRL; // Real.. whichever prec
-//typedef WilsonImpl<vComplexF, FundamentalRepresentation, CoeffComplexHalfComms > ZWilsonImplFH; // Float
-//typedef WilsonImpl<vComplexD, FundamentalRepresentation, CoeffComplexHalfComms > ZWilsonImplDF; // Double
+typedef WilsonImpl<vComplex,  FundamentalRepresentation, CoeffComplexHalfComms > ZWilsonImplRL; // Real.. whichever prec
+typedef WilsonImpl<vComplexF, FundamentalRepresentation, CoeffComplexHalfComms > ZWilsonImplFH; // Float
+typedef WilsonImpl<vComplexD, FundamentalRepresentation, CoeffComplexHalfComms > ZWilsonImplDF; // Double
 
 typedef WilsonImpl<vComplex,  AdjointRepresentation, CoeffReal > WilsonAdjImplR;   // Real.. whichever prec
 typedef WilsonImpl<vComplexF, AdjointRepresentation, CoeffReal > WilsonAdjImplF;  // Float
@@ -880,7 +880,7 @@ void CayleyFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
  }

  std::vector<RealD> G_s(Ls,1.0);
-  RealD sign = 1; // sign flip for vector/tadpole
+  Integer sign = 1; // sign flip for vector/tadpole
  if ( curr_type == Current::Axial ) {
    for(int s=0;s<Ls/2;s++){
      G_s[s] = -1.0;
@@ -901,8 +901,8 @@ void CayleyFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
  for(int s=0;s<Ls;s++){

    int sp = (s+1)%Ls;
-    //    int sr = Ls-1-s;
-    //    int srp= (sr+1)%Ls;
+    int sr = Ls-1-s;
+    int srp= (sr+1)%Ls;

    // Mobius parameters
    auto b=this->bs[s];
@@ -73,17 +73,17 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>


 #undef INTERIOR_AND_EXTERIOR
@@ -102,17 +102,17 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>


 #undef INTERIOR_AND_EXTERIOR
@@ -131,17 +131,17 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>



@@ -165,17 +165,17 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>


 #undef INTERIOR_AND_EXTERIOR
@@ -194,17 +194,17 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFiel
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>


 #undef INTERIOR_AND_EXTERIOR
@@ -223,17 +223,17 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFiel
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-//#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>



@@ -280,17 +280,17 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<WilsonImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<ZWilsonImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>


 #undef INTERIOR_AND_EXTERIOR
@@ -309,17 +309,17 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<WilsonImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<ZWilsonImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>


 #undef INTERIOR_AND_EXTERIOR
@@ -338,17 +338,17 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<WilsonImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<ZWilsonImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>


 /////////////////////////////////////////////////////////////////
@@ -371,17 +371,17 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<WilsonImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>


 #undef INTERIOR_AND_EXTERIOR
@@ -400,17 +400,17 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFiel
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<WilsonImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>


 #undef INTERIOR_AND_EXTERIOR
@@ -429,17 +429,17 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFiel
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<WilsonImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<WilsonImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>

-// #pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
-// template<> void
-// WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-// 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-// #include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>
+#pragma GCC optimize ("-O3", "-fno-schedule-insns", "-fno-schedule-insns2")
+template<> void
+WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBodyA64FX.h>



@@ -74,15 +74,15 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+
+template<> void 
+WilsonKernels<ZWilsonImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #define INTERIOR
@@ -97,15 +97,15 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+
+template<> void 
+WilsonKernels<ZWilsonImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>


 #undef INTERIOR_AND_EXTERIOR
@@ -121,15 +121,15 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+
+template<> void 
+WilsonKernels<ZWilsonImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
      
 /////////////////////////////////////////////////////////////////
 // XYZT vectorised, dag Kernel, single
@@ -148,15 +148,15 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+
+template<> void 
+WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #define INTERIOR
@@ -171,15 +171,15 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFiel
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+
+template<> void 
+WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #undef INTERIOR
@@ -194,15 +194,15 @@ WilsonKernels<ZWilsonImplF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFiel
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
 				    
-//template<> void
-//WilsonKernels<WilsonImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//
-//template<> void
-//WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+				    
+template<> void 
+WilsonKernels<ZWilsonImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
 				    
 #undef MAYBEPERM
 #undef MULT_2SPIN
@@ -228,14 +228,14 @@ WilsonKernels<ZDomainWallVec5dImplF>::AsmDhopSite(StencilView &st, DoubledGaugeF
 							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #define INTERIOR
@@ -249,14 +249,14 @@ WilsonKernels<ZDomainWallVec5dImplF>::AsmDhopSiteInt(StencilView &st, DoubledGau
 							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #undef INTERIOR
@@ -273,15 +273,15 @@ WilsonKernels<ZDomainWallVec5dImplF>::AsmDhopSiteExt(StencilView &st, DoubledGau
 							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
 				    
-//template<> void
-//WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+				    
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
 				    
 /////////////////////////////////////////////////////////////////
 // Ls vectorised, dag Kernel, single
@@ -299,14 +299,14 @@ WilsonKernels<ZDomainWallVec5dImplF>::AsmDhopSiteDag(StencilView &st, DoubledGau
 							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #define INTERIOR
@@ -320,14 +320,14 @@ WilsonKernels<ZDomainWallVec5dImplF>::AsmDhopSiteDagInt(StencilView &st, Doubled
 							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #undef INTERIOR
@@ -341,14 +341,14 @@ WilsonKernels<ZDomainWallVec5dImplF>::AsmDhopSiteDagExt(StencilView &st, Doubled
 							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplFH>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #endif  // VEC 5D

@@ -392,14 +392,14 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<WilsonImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZWilsonImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZWilsonImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #define INTERIOR
@@ -413,14 +413,14 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<WilsonImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZWilsonImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZWilsonImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #undef INTERIOR
@@ -434,14 +434,14 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
      
-//template<> void
-//WilsonKernels<WilsonImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZWilsonImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZWilsonImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
      
 /////////////////////////////////////////////////////////////////
 // XYZT vectorised, dag Kernel, single
@@ -459,14 +459,14 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldVi
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<WilsonImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #define INTERIOR
@@ -480,14 +480,14 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFiel
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<WilsonImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #undef INTERIOR
@@ -501,14 +501,14 @@ WilsonKernels<ZWilsonImplD>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFiel
 						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
 				    
-//template<> void
-//WilsonKernels<WilsonImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<WilsonImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZWilsonImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+						int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
 				    
 #undef MAYBEPERM
 #undef MULT_2SPIN
@@ -533,14 +533,14 @@ WilsonKernels<ZDomainWallVec5dImplD>::AsmDhopSite(StencilView &st, DoubledGaugeF
 							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSite(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #define INTERIOR
@@ -554,14 +554,14 @@ WilsonKernels<ZDomainWallVec5dImplD>::AsmDhopSiteInt(StencilView &st, DoubledGau
 							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSiteInt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #undef INTERIOR
@@ -577,14 +577,14 @@ WilsonKernels<ZDomainWallVec5dImplD>::AsmDhopSiteExt(StencilView &st, DoubledGau
 							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
 				    
-//template<> void
-//WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
-//							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSiteExt(StencilView &st, DoubledGaugeFieldView &U, SiteHalfSpinor *buf,
+							 int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
 				    
 /////////////////////////////////////////////////////////////////
 // Ls vectorised, dag Kernel, single
@@ -602,14 +602,14 @@ WilsonKernels<ZDomainWallVec5dImplD>::AsmDhopSiteDag(StencilView &st, DoubledGau
 							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSiteDag(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #define INTERIOR
@@ -623,14 +623,14 @@ WilsonKernels<ZDomainWallVec5dImplD>::AsmDhopSiteDagInt(StencilView &st, Doubled
 							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSiteDagInt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #undef INTERIOR_AND_EXTERIOR
 #undef INTERIOR
@@ -645,14 +645,14 @@ WilsonKernels<ZDomainWallVec5dImplD>::AsmDhopSiteDagExt(StencilView &st, Doubled
 							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
 #include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

-//template<> void
-//WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
-//template<> void
-//WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
-//							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
-//#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<DomainWallVec5dImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>
+template<> void 
+WilsonKernels<ZDomainWallVec5dImplDF>::AsmDhopSiteDagExt(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf,
+							    int ss,int ssU,int Ls,int Ns,const FermionFieldView &in, FermionFieldView &out)
+#include <qcd/action/fermion/implementation/WilsonKernelsAsmBody.h>

 #endif  // VEC 5D

@@ -0,0 +1 @@
+../CayleyFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../ContinuedFractionFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../DomainWallEOFAFermionInstantiation.cc.master
@@ -0,0 +1 @@
+../MobiusEOFAFermionInstantiation.cc.master
@@ -0,0 +1 @@
+../PartialFractionFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../WilsonCloverFermionInstantiation.cc.master
@@ -0,0 +1 @@
+../WilsonFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../WilsonFermionInstantiation.cc.master
@@ -0,0 +1 @@
+../WilsonKernelsInstantiationGparity.cc.master
@@ -0,0 +1 @@
+../WilsonTMFermionInstantiation.cc.master
@@ -0,0 +1 @@
+#define IMPLEMENTATION GparityWilsonImplDF
@@ -0,0 +1 @@
+../CayleyFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../ContinuedFractionFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../DomainWallEOFAFermionInstantiation.cc.master
@@ -0,0 +1 @@
+../MobiusEOFAFermionInstantiation.cc.master
@@ -0,0 +1 @@
+../PartialFractionFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../WilsonCloverFermionInstantiation.cc.master
@@ -0,0 +1 @@
+../WilsonFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../WilsonFermionInstantiation.cc.master
@@ -0,0 +1 @@
+../WilsonKernelsInstantiationGparity.cc.master
@@ -0,0 +1 @@
+../WilsonTMFermionInstantiation.cc.master
@@ -0,0 +1 @@
+#define IMPLEMENTATION GparityWilsonImplFH
@@ -0,0 +1 @@
+../CayleyFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../ContinuedFractionFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../DomainWallEOFAFermionInstantiation.cc.master
@@ -0,0 +1 @@
+../MobiusEOFAFermionInstantiation.cc.master
@@ -0,0 +1 @@
+../PartialFractionFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../WilsonCloverFermionInstantiation.cc.master
@@ -0,0 +1 @@
+../WilsonFermion5DInstantiation.cc.master
@@ -0,0 +1 @@
+../WilsonFermionInstantiation.cc.master
@@ -2,11 +2,14 @@

 Grid physics library, www.github.com/paboyle/Grid

-Source file: ./lib/serialisation/BaseIO.h
+Source file: ./lib/qcd/action/fermion/WilsonKernels.cc

-Copyright (C) 2015
+Copyright (C) 2015, 2020

-Author: Michael Marshall <michael.marshall@ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Nils Meyer <nils.meyer@ur.de> Regensburg University

 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -22,14 +25,27 @@ 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
+See the full license in the file "LICENSE" in the top level distribution
+directory
 *************************************************************************************/
 /*  END LEGAL */
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>

-#include <Grid/GridCore.h>
+#ifndef AVX512
+#ifndef QPX
+#ifndef A64FX
+#ifndef A64FXFIXEDSIZE
+#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
+#endif
+#endif
+#endif
+#endif

-NAMESPACE_BEGIN(Grid)
+NAMESPACE_BEGIN(Grid);

-std::uint64_t EigenIO::EigenResizeCounter(0);
+#include "impl.h"
+template class WilsonKernels<IMPLEMENTATION>;

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

 DWF_IMPL_LIST=" \
 	   WilsonImplF \
 	   WilsonImplD \
+	   WilsonImplFH \
+	   WilsonImplDF \
 	   ZWilsonImplF \
-	   ZWilsonImplD "
+	   ZWilsonImplD \
+	   ZWilsonImplFH \
+	   ZWilsonImplDF "

 GDWF_IMPL_LIST=" \
 	   GparityWilsonImplF \
-	   GparityWilsonImplD "
+	   GparityWilsonImplD \
+	   GparityWilsonImplFH \
+	   GparityWilsonImplDF"
+

 IMPL_LIST="$STAG_IMPL_LIST  $WILSON_IMPL_LIST $DWF_IMPL_LIST $GDWF_IMPL_LIST"

@@ -69,11 +69,6 @@ public:
    return PeriodicBC::ShiftStaple(Link,mu);
  }

-  //Same as Cshift for periodic BCs
-  static inline GaugeLinkField CshiftLink(const GaugeLinkField &Link, int mu, int shift){
-    return PeriodicBC::CshiftLink(Link,mu,shift);
-  }
-
  static inline bool isPeriodicGaugeField(void) { return true; }
 };

@@ -115,11 +110,6 @@ public:
      return PeriodicBC::CovShiftBackward(Link, mu, field);
  }

-  //If mu is a conjugate BC direction
-  //Out(x) = U^dag_\mu(x-mu)  | x_\mu != 0
-  //       = U^T_\mu(L-1)  | x_\mu == 0
-  //else
-  //Out(x) = U^dag_\mu(x-mu mod L)
  static inline GaugeLinkField
  CovShiftIdentityBackward(const GaugeLinkField &Link, int mu)
  {
@@ -139,13 +129,6 @@ public:
      return PeriodicBC::CovShiftIdentityForward(Link,mu);
  }

-
-  //If mu is a conjugate BC direction
-  //Out(x) = S_\mu(x+mu)  | x_\mu != L-1
-  //       = S*_\mu(x+mu)  | x_\mu == L-1
-  //else
-  //Out(x) = S_\mu(x+mu mod L)
-  //Note: While this is used for Staples it is also applicable for shifting gauge links or gauge transformation matrices
  static inline GaugeLinkField ShiftStaple(const GaugeLinkField &Link, int mu)
  {
    assert(_conjDirs.size() == Nd);
@@ -155,27 +138,6 @@ public:
      return PeriodicBC::ShiftStaple(Link,mu);
  }

-  //Boundary-aware C-shift of gauge links / gauge transformation matrices
-  //For conjugate BC direction
-  //shift = 1
-  //Out(x) = U_\mu(x+\hat\mu)  | x_\mu != L-1
-  //       = U*_\mu(0)  | x_\mu == L-1
-  //shift = -1
-  //Out(x) = U_\mu(x-mu)  | x_\mu != 0
-  //       = U*_\mu(L-1)  | x_\mu == 0
-  //else
-  //shift = 1
-  //Out(x) = U_\mu(x+\hat\mu mod L)
-  //shift = -1
-  //Out(x) = U_\mu(x-\hat\mu mod L)
-  static inline GaugeLinkField CshiftLink(const GaugeLinkField &Link, int mu, int shift){
-    assert(_conjDirs.size() == Nd);
-    if(_conjDirs[mu]) 
-      return ConjugateBC::CshiftLink(Link,mu,shift);
-    else     
-      return PeriodicBC::CshiftLink(Link,mu,shift);
-  }
-
  static inline void       setDirections(std::vector<int> &conjDirs) { _conjDirs=conjDirs; }
  static inline std::vector<int> getDirections(void) { return _conjDirs; }
  static inline bool isPeriodicGaugeField(void) { return false; }
				`@@ -0,0 +1 @@`
				`../ContinuedFractionFermion5DInstantiation.cc.master`
				`@@ -0,0 +1 @@`
				`../DomainWallEOFAFermionInstantiation.cc.master`
				`@@ -0,0 +1 @@`
				`../MobiusEOFAFermionInstantiation.cc.master`
				`@@ -0,0 +1 @@`
				`../PartialFractionFermion5DInstantiation.cc.master`
				`@@ -0,0 +1 @@`
				`../WilsonCloverFermionInstantiation.cc.master`
				`@@ -0,0 +1 @@`
				`../WilsonKernelsInstantiationGparity.cc.master`