Modified the Dirac Kernel class to compile with different number of colours

Added the general push_back functionality to accomodate for all defined representations Compiles, not tested
2025-11-04 05:54:32 +00:00 · 2016-07-18 16:36:28 +01:00
parent 9c77bb69a5
commit b93e18ed50
12 changed files with 579 additions and 425 deletions
--- a/lib/qcd/action/ActionBase.h
+++ b/lib/qcd/action/ActionBase.h
@@ -66,6 +66,20 @@ Phi(&_Grid), pRNG(_pRNG) {
 };
 */

+// Indexing of tuple types
+template <class T, class Tuple>
+struct Index;
+
+template <class T, class... Types>
+struct Index<T, std::tuple<T, Types...>> {
+  static const std::size_t value = 0;
+};
+
+template <class T, class U, class... Types>
+struct Index<T, std::tuple<U, Types...>> {
+  static const std::size_t value = 1 + Index<T, std::tuple<Types...>>::value;
+};
+
 template <class GaugeField>
 struct ActionLevel {
 public:
@@ -99,38 +113,39 @@ struct ActionLevelHirep {
  // representation fields
  typedef typename AccessTypes<Action, Repr>::VectorCollection action_collection;
  action_collection actions_hirep;
-  typedef typename  AccessTypes<Action, Repr>::ClassCollection actions_hirep_ptrs_type;
+  typedef typename  AccessTypes<Action, Repr>::FieldTypeCollection action_hirep_types;

  std::vector<ActPtr>& actions;

  // Temporary conversion between ActionLevel and ActionLevelHirep
  ActionLevelHirep(ActionLevel<GaugeField>& AL ):actions(AL.actions), multiplier(AL.multiplier){}

-
-
  ActionLevelHirep(unsigned int mul = 1) : actions(std::get<0>(actions_hirep)), multiplier(mul) {
    // initialize the hirep vectors to zero.
    //apply(this->resize, actions_hirep, 0); //need a working resize
    assert(mul >= 1);
  };

-  void push_back(ActPtr ptr) { actions.push_back(ptr); }
+  //void push_back(ActPtr ptr) { actions.push_back(ptr); }

-// SFINAE construct, check
-  template <class actionpointer, size_t N>
-  void push_back(actionpointer ptr, decltype(std::tuple_element<N, actions_hirep_ptrs_type>::value)* = 0) {
-    //insert only in the correct vector
-    std::get<N>(actions_hirep).push_back(ptr);
+
+
+  template < class Field >
+  void push_back(Action<Field>* ptr) {
+    // insert only in the correct vector
+    std::get< Index < Field, action_hirep_types>::value >(actions_hirep).push_back(ptr);
  };

+ 
+
  template < class ActPtr>
  static void resize(ActPtr ap, unsigned int n){
    ap->resize(n);

  }

-  template <std::size_t I>
-  auto getRepresentation(Repr& R)->decltype(std::get<I>(R).U)  {return std::get<I>(R).U;}
+  //template <std::size_t I>
+  //auto getRepresentation(Repr& R)->decltype(std::get<I>(R).U)  {return std::get<I>(R).U;}

  // Loop on tuple for a callable function
  template <std::size_t I = 1, typename Callable, typename ...Args>
--- a/lib/qcd/action/Actions.h
+++ b/lib/qcd/action/Actions.h
@@ -113,6 +113,10 @@ typedef SymanzikGaugeAction<ConjugateGimplD>        ConjugateSymanzikGaugeAction
  template class A<GparityWilsonImplF>;		\
  template class A<GparityWilsonImplD>;		

+#define AdjointFermOpTemplateInstantiate(A) \
+  template class A<WilsonAdjImplF>; \
+  template class A<WilsonAdjImplD>; 
+
 #define GparityFermOpTemplateInstantiate(A) 

 ////////////////////////////////////////////
@@ -157,6 +161,10 @@ typedef WilsonFermion<WilsonImplR> WilsonFermionR;
 typedef WilsonFermion<WilsonImplF> WilsonFermionF;
 typedef WilsonFermion<WilsonImplD> WilsonFermionD;

+typedef WilsonFermion<WilsonAdjImplR> WilsonAdjFermionR;
+typedef WilsonFermion<WilsonAdjImplF> WilsonAdjFermionF;
+typedef WilsonFermion<WilsonAdjImplD> WilsonAdjFermionD;
+
 typedef WilsonTMFermion<WilsonImplR> WilsonTMFermionR;
 typedef WilsonTMFermion<WilsonImplF> WilsonTMFermionF;
 typedef WilsonTMFermion<WilsonImplD> WilsonTMFermionD;
--- a/lib/qcd/action/fermion/FermionOperatorImpl.h
+++ b/lib/qcd/action/fermion/FermionOperatorImpl.h
@@ -115,22 +115,23 @@ template <class S, class Representation = FundamentalRepresentation >
 class WilsonImpl
    : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation::Dimension > > {
 public:
-  static const int Nrepresentation = Representation::Dimension;
-  typedef PeriodicGaugeImpl<GaugeImplTypes<S, Representation::Dimension > > Gimpl;
+  static const int Dimension = Representation::Dimension;
+//  static const int Nrepresentation = Representation::Dimension;
+  typedef PeriodicGaugeImpl<GaugeImplTypes<S, Dimension > > Gimpl;
  
  //Necessary?
-  constexpr bool is_fundamental() const{return Representation::Dimension == Nc ? 1 : 0;}
+  constexpr bool is_fundamental() const{return Dimension == Nc ? 1 : 0;}

  INHERIT_GIMPL_TYPES(Gimpl);

  template <typename vtype>
-  using iImplSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
+  using iImplSpinor = iScalar<iVector<iVector<vtype, Dimension>, Ns> >;
  template <typename vtype>
  using iImplHalfSpinor =
-      iScalar<iVector<iVector<vtype, Nrepresentation>, Nhs> >;
+      iScalar<iVector<iVector<vtype, Dimension>, Nhs> >;
  template <typename vtype>
  using iImplDoubledGaugeField =
-      iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>;
+      iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>;

  typedef iImplSpinor<Simd> SiteSpinor;
  typedef iImplHalfSpinor<Simd> SiteHalfSpinor;
@@ -214,6 +215,8 @@ template <class S, int Nrepresentation = Nc>
 class DomainWallRedBlack5dImpl
    : public PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > {
 public:
+  static const int Dimension = Nrepresentation;
+
  typedef PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > Gimpl;

  INHERIT_GIMPL_TYPES(Gimpl);
@@ -318,6 +321,7 @@ template <class S, int Nrepresentation>
 class GparityWilsonImpl
    : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresentation> > {
 public:
+  static const int Dimension = Nrepresentation;
  typedef ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresentation> > Gimpl;

  INHERIT_GIMPL_TYPES(Gimpl);
--- a/lib/qcd/action/fermion/WilsonFermion.cc
+++ b/lib/qcd/action/fermion/WilsonFermion.cc
@@ -308,6 +308,7 @@ void WilsonFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder &lo,
 };

 FermOpTemplateInstantiate(WilsonFermion);
+AdjointFermOpTemplateInstantiate(WilsonFermion);
 GparityFermOpTemplateInstantiate(WilsonFermion);
 }
 }
--- a/lib/qcd/action/fermion/WilsonFermion.h
+++ b/lib/qcd/action/fermion/WilsonFermion.h
@@ -148,6 +148,8 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {

 typedef WilsonFermion<WilsonImplF> WilsonFermionF;
 typedef WilsonFermion<WilsonImplD> WilsonFermionD;
+
+
 }
 }
 #endif
--- a/lib/qcd/action/fermion/WilsonKernels.cc
+++ b/lib/qcd/action/fermion/WilsonKernels.cc
@@ -1,72 +1,94 @@
-    /*************************************************************************************
+/*************************************************************************************

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

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

-    Copyright (C) 2015
+Copyright (C) 2015

 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
 Author: paboyle <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 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.
+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.
+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 */
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
 #include <Grid.h>
 namespace Grid {
 namespace QCD {

-  int WilsonKernelsStatic::HandOpt;
-  int WilsonKernelsStatic::AsmOpt;
+int WilsonKernelsStatic::HandOpt;
+int WilsonKernelsStatic::AsmOpt;

-template<class Impl> 
-WilsonKernels<Impl>::WilsonKernels(const ImplParams &p): Base(p) {};
+template <class Impl>
+WilsonKernels<Impl>::WilsonKernels(const ImplParams &p) : Base(p){};

-template<class Impl> 
-void WilsonKernels<Impl>::DiracOptDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
-						  std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-						  int sF,int sU,int Ls, int Ns, const FermionField &in, FermionField &out)
-{
+/*
+template <class Impl>
+typename std::enable_if<Impl::Dimension == 3>::type WilsonKernels<Impl>::DiracOptDhopSite(
+    StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+    std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf, int sF,
+    int sU, int Ls, int Ns, const FermionField &in, FermionField &out) {
 #ifdef AVX512
-  if ( AsmOpt ) {
-
-    WilsonKernels<Impl>::DiracOptAsmDhopSite(st,lo,U,buf,sF,sU,Ls,Ns,in,out);
+  if (AsmOpt) {
+    WilsonKernels<Impl>::DiracOptAsmDhopSite(st, lo, U, buf, sF, sU, Ls, Ns, in,
+                                             out);

  } else {
 #else
-  {  
+  {
 #endif
-    for(int site=0;site<Ns;site++) {
-      for(int s=0;s<Ls;s++) {
-	if (HandOpt) WilsonKernels<Impl>::DiracOptHandDhopSite(st,lo,U,buf,sF,sU,in,out);
-	else         WilsonKernels<Impl>::DiracOptGenericDhopSite(st,lo,U,buf,sF,sU,in,out);
-	sF++;
+    for (int site = 0; site < Ns; site++) {
+      for (int s = 0; s < Ls; s++) {
+        if (HandOpt)
+          WilsonKernels<Impl>::DiracOptHandDhopSite(st, lo, U, buf, sF, sU, in,
+                                                    out);
+        else
+          WilsonKernels<Impl>::DiracOptGenericDhopSite(st, lo, U, buf, sF, sU,
+                                                       in, out);
+        sF++;
      }
      sU++;
    }
-
  }
 }

+template <class Impl>
+typename std::enable_if<Impl::Dimension != 3>::type WilsonKernels<Impl>::DiracOptDhopSite(
+    StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+    std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf, int sF,
+    int sU, int Ls, int Ns, const FermionField &in, FermionField &out) {
+  for (int site = 0; site < Ns; site++) {
+    for (int s = 0; s < Ls; s++) {
+      WilsonKernels<Impl>::DiracOptGenericDhopSite(st, lo, U, buf, sF, sU, in,
+                                                   out);
+      sF++;
+    }
+    sU++;
+  }
+}
+
+
 template<class Impl> 
 void WilsonKernels<Impl>::DiracOptDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
-					   std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-					   int sF,int sU,int Ls, int Ns, const FermionField &in, FermionField &out)
+             std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
+             int sF,int sU,int Ls, int Ns, const FermionField &in, FermionField &out,
+             typename std::enable_if<Impl::Dimension == 3, int>::type = 0)
 {
  // No asm implementation yet.
  //  if ( AsmOpt )     WilsonKernels<Impl>::DiracOptAsmDhopSiteDag(st,lo,U,buf,sF,sU,in,out);
@@ -82,17 +104,35 @@ void WilsonKernels<Impl>::DiracOptDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,
 }


-  ////////////////////////////////////////////
-  // Generic implementation; move to different file?
-  ////////////////////////////////////////////

-template<class Impl> 
-void WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
-					   std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-					   int sF,int sU,const FermionField &in, FermionField &out)
-{
-  SiteHalfSpinor  tmp;    
-  SiteHalfSpinor  chi;    
+
+template <class Impl>
+void WilsonKernels<Impl>::DiracOptDhopSiteDag(
+    StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+    std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf, int sF,
+    int sU, int Ls, int Ns, const FermionField &in, FermionField &out,
+    typename std::enable_if<Impl::Dimension != 3, int>::type = 0) {
+  for (int site = 0; site < Ns; site++) {
+    for (int s = 0; s < Ls; s++) {
+      WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(st, lo, U, buf, sF, sU,
+                                                      in, out);
+      sF++;
+    }
+    sU++;
+  }
+}
+*/
+////////////////////////////////////////////
+// Generic implementation; move to different file?
+////////////////////////////////////////////
+
+template <class Impl>
+void WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(
+    StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+    std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf, int sF,
+    int sU, const FermionField &in, FermionField &out) {
+  SiteHalfSpinor tmp;
+  SiteHalfSpinor chi;
  SiteHalfSpinor *chi_p;
  SiteHalfSpinor Uchi;
  SiteSpinor result;
@@ -102,176 +142,175 @@ void WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(StencilImpl &st,LebesgueOrd
  ///////////////////////////
  // Xp
  ///////////////////////////
-  SE=st.GetEntry(ptype,Xp,sF);
+  SE = st.GetEntry(ptype, Xp, sF);

-  if (SE->_is_local ) { 
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjXp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjXp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjXp(chi,in._odata[SE->_offset]);
+      spProjXp(chi, in._odata[SE->_offset]);
    }
-  } else { 
-    chi_p=&buf[SE->_offset];
+  } else {
+    chi_p = &buf[SE->_offset];
  }
-  
-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Xp,SE,st);
-  spReconXp(result,Uchi);
-    
+
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Xp, SE, st);
+  spReconXp(result, Uchi);
+
  ///////////////////////////
  // Yp
  ///////////////////////////
-  SE=st.GetEntry(ptype,Yp,sF);
+  SE = st.GetEntry(ptype, Yp, sF);

-  if ( SE->_is_local ) { 
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjYp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjYp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjYp(chi,in._odata[SE->_offset]);
+      spProjYp(chi, in._odata[SE->_offset]);
    }
-  } else { 
-    chi_p=&buf[SE->_offset];
+  } else {
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Yp,SE,st);
-  accumReconYp(result,Uchi);
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Yp, SE, st);
+  accumReconYp(result, Uchi);

  ///////////////////////////
  // Zp
  ///////////////////////////
-  SE=st.GetEntry(ptype,Zp,sF);
+  SE = st.GetEntry(ptype, Zp, sF);

-  if ( SE->_is_local ) { 
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjZp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjZp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjZp(chi,in._odata[SE->_offset]);
+      spProjZp(chi, in._odata[SE->_offset]);
    }
-  } else { 
-    chi_p=&buf[SE->_offset];
+  } else {
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Zp,SE,st);
-  accumReconZp(result,Uchi);
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Zp, SE, st);
+  accumReconZp(result, Uchi);

  ///////////////////////////
  // Tp
  ///////////////////////////
-  SE=st.GetEntry(ptype,Tp,sF);
+  SE = st.GetEntry(ptype, Tp, sF);

-  if ( SE->_is_local ) {
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjTp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjTp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjTp(chi,in._odata[SE->_offset]);
+      spProjTp(chi, in._odata[SE->_offset]);
    }
  } else {
-    chi_p=&buf[SE->_offset];
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Tp,SE,st);
-  accumReconTp(result,Uchi);
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Tp, SE, st);
+  accumReconTp(result, Uchi);

  ///////////////////////////
  // Xm
  ///////////////////////////
-  SE=st.GetEntry(ptype,Xm,sF);
+  SE = st.GetEntry(ptype, Xm, sF);

-  if ( SE->_is_local ) {
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjXm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjXm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjXm(chi,in._odata[SE->_offset]);
+      spProjXm(chi, in._odata[SE->_offset]);
    }
  } else {
-    chi_p=&buf[SE->_offset];
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Xm,SE,st);
-  accumReconXm(result,Uchi);
-  
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Xm, SE, st);
+  accumReconXm(result, Uchi);
+
  ///////////////////////////
  // Ym
  ///////////////////////////
-  SE=st.GetEntry(ptype,Ym,sF);
+  SE = st.GetEntry(ptype, Ym, sF);

-  if ( SE->_is_local ) {
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjYm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjYm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjYm(chi,in._odata[SE->_offset]);
+      spProjYm(chi, in._odata[SE->_offset]);
    }
  } else {
-    chi_p=&buf[SE->_offset];
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Ym,SE,st);
-  accumReconYm(result,Uchi);
-  
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Ym, SE, st);
+  accumReconYm(result, Uchi);
+
  ///////////////////////////
  // Zm
  ///////////////////////////
-  SE=st.GetEntry(ptype,Zm,sF);
+  SE = st.GetEntry(ptype, Zm, sF);

-  if ( SE->_is_local ) {
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjZm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjZm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjZm(chi,in._odata[SE->_offset]);
+      spProjZm(chi, in._odata[SE->_offset]);
    }
  } else {
-    chi_p=&buf[SE->_offset];
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Zm,SE,st);
-  accumReconZm(result,Uchi);
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Zm, SE, st);
+  accumReconZm(result, Uchi);

  ///////////////////////////
  // Tm
  ///////////////////////////
-  SE=st.GetEntry(ptype,Tm,sF);
+  SE = st.GetEntry(ptype, Tm, sF);

-  if ( SE->_is_local ) {
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjTm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
-    } else { 
-      spProjTm(chi,in._odata[SE->_offset]);
+    if (SE->_permute) {
+      spProjTm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
+    } else {
+      spProjTm(chi, in._odata[SE->_offset]);
    }
  } else {
-    chi_p=&buf[SE->_offset];
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Tm,SE,st);
-  accumReconTm(result,Uchi);
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Tm, SE, st);
+  accumReconTm(result, Uchi);

-  vstream(out._odata[sF],result);
+  vstream(out._odata[sF], result);
 };

-
-  // Need controls to do interior, exterior, or both
-template<class Impl> 
-void WilsonKernels<Impl>::DiracOptGenericDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
-						  std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-						  int sF,int sU,const FermionField &in, FermionField &out)
-{
-  SiteHalfSpinor  tmp;    
-  SiteHalfSpinor  chi;    
-  SiteHalfSpinor *chi_p;    
+// Need controls to do interior, exterior, or both
+template <class Impl>
+void WilsonKernels<Impl>::DiracOptGenericDhopSite(
+    StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+    std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf, int sF,
+    int sU, const FermionField &in, FermionField &out) {
+  SiteHalfSpinor tmp;
+  SiteHalfSpinor chi;
+  SiteHalfSpinor *chi_p;
  SiteHalfSpinor Uchi;
  SiteSpinor result;
  StencilEntry *SE;
@@ -280,299 +319,299 @@ void WilsonKernels<Impl>::DiracOptGenericDhopSite(StencilImpl &st,LebesgueOrder
  ///////////////////////////
  // Xp
  ///////////////////////////
-  SE=st.GetEntry(ptype,Xm,sF);
+  SE = st.GetEntry(ptype, Xm, sF);

-  if ( SE->_is_local ) { 
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjXp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjXp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjXp(chi,in._odata[SE->_offset]);
+      spProjXp(chi, in._odata[SE->_offset]);
    }
-  } else { 
-    chi_p=&buf[SE->_offset];
+  } else {
+    chi_p = &buf[SE->_offset];
  }
-  
-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Xm,SE,st);
-  spReconXp(result,Uchi);
-    
+
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Xm, SE, st);
+  spReconXp(result, Uchi);
+
  ///////////////////////////
  // Yp
  ///////////////////////////
-  SE=st.GetEntry(ptype,Ym,sF);
+  SE = st.GetEntry(ptype, Ym, sF);

-  if ( SE->_is_local ) { 
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjYp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjYp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjYp(chi,in._odata[SE->_offset]);
+      spProjYp(chi, in._odata[SE->_offset]);
    }
-  } else { 
-    chi_p=&buf[SE->_offset];
+  } else {
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Ym,SE,st);
-  accumReconYp(result,Uchi);
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Ym, SE, st);
+  accumReconYp(result, Uchi);

  ///////////////////////////
  // Zp
  ///////////////////////////
-  SE=st.GetEntry(ptype,Zm,sF);
+  SE = st.GetEntry(ptype, Zm, sF);

-  if ( SE->_is_local ) { 
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjZp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjZp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjZp(chi,in._odata[SE->_offset]);
+      spProjZp(chi, in._odata[SE->_offset]);
    }
-  } else { 
-    chi_p=&buf[SE->_offset];
+  } else {
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Zm,SE,st);
-  accumReconZp(result,Uchi);
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Zm, SE, st);
+  accumReconZp(result, Uchi);

  ///////////////////////////
  // Tp
  ///////////////////////////
-  SE=st.GetEntry(ptype,Tm,sF);
+  SE = st.GetEntry(ptype, Tm, sF);

-  if ( SE->_is_local ) {
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjTp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjTp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjTp(chi,in._odata[SE->_offset]);
+      spProjTp(chi, in._odata[SE->_offset]);
    }
  } else {
-    chi_p=&buf[SE->_offset];
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Tm,SE,st);
-  accumReconTp(result,Uchi);
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Tm, SE, st);
+  accumReconTp(result, Uchi);

  ///////////////////////////
  // Xm
  ///////////////////////////
-  SE=st.GetEntry(ptype,Xp,sF);
+  SE = st.GetEntry(ptype, Xp, sF);

-  if ( SE->_is_local ) {
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjXm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjXm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjXm(chi,in._odata[SE->_offset]);
+      spProjXm(chi, in._odata[SE->_offset]);
    }
  } else {
-    chi_p=&buf[SE->_offset];
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Xp,SE,st);
-  accumReconXm(result,Uchi);
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Xp, SE, st);
+  accumReconXm(result, Uchi);

  ///////////////////////////
  // Ym
  ///////////////////////////
-  SE=st.GetEntry(ptype,Yp,sF);
+  SE = st.GetEntry(ptype, Yp, sF);

-  if ( SE->_is_local ) {
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjYm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjYm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjYm(chi,in._odata[SE->_offset]);
+      spProjYm(chi, in._odata[SE->_offset]);
    }
  } else {
-    chi_p=&buf[SE->_offset];
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Yp,SE,st);
-  accumReconYm(result,Uchi);
-  
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Yp, SE, st);
+  accumReconYm(result, Uchi);
+
  ///////////////////////////
  // Zm
  ///////////////////////////
-  SE=st.GetEntry(ptype,Zp,sF);
+  SE = st.GetEntry(ptype, Zp, sF);

-  if ( SE->_is_local ) {
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjZm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
+    if (SE->_permute) {
+      spProjZm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
    } else {
-      spProjZm(chi,in._odata[SE->_offset]);
+      spProjZm(chi, in._odata[SE->_offset]);
    }
  } else {
-    chi_p=&buf[SE->_offset];
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Zp,SE,st);
-  accumReconZm(result,Uchi);
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Zp, SE, st);
+  accumReconZm(result, Uchi);

  ///////////////////////////
  // Tm
  ///////////////////////////
-  SE=st.GetEntry(ptype,Tp,sF);
+  SE = st.GetEntry(ptype, Tp, sF);

-  if ( SE->_is_local ) {
+  if (SE->_is_local) {
    chi_p = &chi;
-    if ( SE->_permute ) {
-      spProjTm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
-    } else { 
-      spProjTm(chi,in._odata[SE->_offset]);
+    if (SE->_permute) {
+      spProjTm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
+    } else {
+      spProjTm(chi, in._odata[SE->_offset]);
    }
  } else {
-    chi_p=&buf[SE->_offset];
+    chi_p = &buf[SE->_offset];
  }

-  Impl::multLink(Uchi,U._odata[sU],*chi_p,Tp,SE,st);
-  accumReconTm(result,Uchi);
+  Impl::multLink(Uchi, U._odata[sU], *chi_p, Tp, SE, st);
+  accumReconTm(result, Uchi);

-  vstream(out._odata[sF],result);
+  vstream(out._odata[sF], result);
 };

-template<class Impl> 
-void WilsonKernels<Impl>::DiracOptDhopDir(StencilImpl &st,DoubledGaugeField &U,
-					  std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-					  int sF,int sU,const FermionField &in, FermionField &out,int dir,int gamma)
-{
-  SiteHalfSpinor  tmp;    
-  SiteHalfSpinor  chi;    
-  SiteSpinor   result;
+template <class Impl>
+void WilsonKernels<Impl>::DiracOptDhopDir(
+    StencilImpl &st, DoubledGaugeField &U,
+    std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf, int sF,
+    int sU, const FermionField &in, FermionField &out, int dir, int gamma) {
+  SiteHalfSpinor tmp;
+  SiteHalfSpinor chi;
+  SiteSpinor result;
  SiteHalfSpinor Uchi;
  StencilEntry *SE;
  int ptype;

-  SE=st.GetEntry(ptype,dir,sF);
+  SE = st.GetEntry(ptype, dir, sF);

  // Xp
-  if(gamma==Xp){
-    if (  SE->_is_local && SE->_permute ) {
-      spProjXp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
-    } else if ( SE->_is_local ) {
-      spProjXp(chi,in._odata[SE->_offset]);
-    } else { 
-      chi=buf[SE->_offset];
+  if (gamma == Xp) {
+    if (SE->_is_local && SE->_permute) {
+      spProjXp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
+    } else if (SE->_is_local) {
+      spProjXp(chi, in._odata[SE->_offset]);
+    } else {
+      chi = buf[SE->_offset];
    }
-    Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
-    spReconXp(result,Uchi);
+    Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
+    spReconXp(result, Uchi);
  }

  // Yp
-  if ( gamma==Yp ){
-    if (  SE->_is_local && SE->_permute ) {
-      spProjYp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
-    } else if ( SE->_is_local ) {
-      spProjYp(chi,in._odata[SE->_offset]);
-    } else { 
-      chi=buf[SE->_offset];
+  if (gamma == Yp) {
+    if (SE->_is_local && SE->_permute) {
+      spProjYp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
+    } else if (SE->_is_local) {
+      spProjYp(chi, in._odata[SE->_offset]);
+    } else {
+      chi = buf[SE->_offset];
    }
-    Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
-    spReconYp(result,Uchi);
+    Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
+    spReconYp(result, Uchi);
  }
-  
+
  // Zp
-  if ( gamma ==Zp ){
-    if (  SE->_is_local && SE->_permute ) {
-      spProjZp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
-    } else if ( SE->_is_local ) {
-      spProjZp(chi,in._odata[SE->_offset]);
-    } else { 
-      chi=buf[SE->_offset];
+  if (gamma == Zp) {
+    if (SE->_is_local && SE->_permute) {
+      spProjZp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
+    } else if (SE->_is_local) {
+      spProjZp(chi, in._odata[SE->_offset]);
+    } else {
+      chi = buf[SE->_offset];
    }
-    Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
-    spReconZp(result,Uchi);
+    Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
+    spReconZp(result, Uchi);
  }
-  
+
  // Tp
-  if ( gamma ==Tp ){
-    if (  SE->_is_local && SE->_permute ) {
-      spProjTp(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
-    } else if ( SE->_is_local ) {
-      spProjTp(chi,in._odata[SE->_offset]);
-    } else { 
-      chi=buf[SE->_offset];
+  if (gamma == Tp) {
+    if (SE->_is_local && SE->_permute) {
+      spProjTp(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
+    } else if (SE->_is_local) {
+      spProjTp(chi, in._odata[SE->_offset]);
+    } else {
+      chi = buf[SE->_offset];
    }
-    Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
-    spReconTp(result,Uchi);
+    Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
+    spReconTp(result, Uchi);
  }

  // Xm
-  if ( gamma==Xm ){
-    if (  SE->_is_local && SE->_permute ) {
-      spProjXm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
-    } else if ( SE->_is_local ) {
-      spProjXm(chi,in._odata[SE->_offset]);
-    } else { 
-      chi=buf[SE->_offset];
+  if (gamma == Xm) {
+    if (SE->_is_local && SE->_permute) {
+      spProjXm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
+    } else if (SE->_is_local) {
+      spProjXm(chi, in._odata[SE->_offset]);
+    } else {
+      chi = buf[SE->_offset];
    }
-    Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
-    spReconXm(result,Uchi);
+    Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
+    spReconXm(result, Uchi);
  }

  // Ym
-  if ( gamma == Ym ){
-    if (  SE->_is_local && SE->_permute ) {
-      spProjYm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
-    } else if ( SE->_is_local ) {
-      spProjYm(chi,in._odata[SE->_offset]);
-    } else { 
-      chi=buf[SE->_offset];
+  if (gamma == Ym) {
+    if (SE->_is_local && SE->_permute) {
+      spProjYm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
+    } else if (SE->_is_local) {
+      spProjYm(chi, in._odata[SE->_offset]);
+    } else {
+      chi = buf[SE->_offset];
    }
-    Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
-    spReconYm(result,Uchi);
+    Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
+    spReconYm(result, Uchi);
  }

  // Zm
-  if ( gamma == Zm ){
-    if (  SE->_is_local && SE->_permute ) {
-      spProjZm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
-    } else if ( SE->_is_local ) {
-      spProjZm(chi,in._odata[SE->_offset]);
-    } else { 
-      chi=buf[SE->_offset];
+  if (gamma == Zm) {
+    if (SE->_is_local && SE->_permute) {
+      spProjZm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
+    } else if (SE->_is_local) {
+      spProjZm(chi, in._odata[SE->_offset]);
+    } else {
+      chi = buf[SE->_offset];
    }
-    Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
-    spReconZm(result,Uchi);
-  }
-  
-  // Tm
-  if ( gamma==Tm ) {
-    if (  SE->_is_local && SE->_permute ) {
-      spProjTm(tmp,in._odata[SE->_offset]);
-      permute(chi,tmp,ptype);
-    } else if ( SE->_is_local ) {
-      spProjTm(chi,in._odata[SE->_offset]);
-    } else { 
-      chi=buf[SE->_offset];
-    }
-    Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
-    spReconTm(result,Uchi);
+    Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
+    spReconZm(result, Uchi);
  }

-  vstream(out._odata[sF],result);
+  // Tm
+  if (gamma == Tm) {
+    if (SE->_is_local && SE->_permute) {
+      spProjTm(tmp, in._odata[SE->_offset]);
+      permute(chi, tmp, ptype);
+    } else if (SE->_is_local) {
+      spProjTm(chi, in._odata[SE->_offset]);
+    } else {
+      chi = buf[SE->_offset];
+    }
+    Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
+    spReconTm(result, Uchi);
+  }
+
+  vstream(out._odata[sF], result);
 }

+FermOpTemplateInstantiate(WilsonKernels);
+AdjointFermOpTemplateInstantiate(WilsonKernels);

-  FermOpTemplateInstantiate(WilsonKernels);
-
-template class WilsonKernels<DomainWallRedBlack5dImplF>;		
+template class WilsonKernels<DomainWallRedBlack5dImplF>;
 template class WilsonKernels<DomainWallRedBlack5dImplD>;
-
-}}
+}
+}
--- a/lib/qcd/action/fermion/WilsonKernels.h
+++ b/lib/qcd/action/fermion/WilsonKernels.h
@@ -1,98 +1,183 @@
-    /*************************************************************************************
+/*************************************************************************************

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

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

-    Copyright (C) 2015
+Copyright (C) 2015

 Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: paboyle <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 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.
+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.
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
-#ifndef  GRID_QCD_DHOP_H
-#define  GRID_QCD_DHOP_H
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef GRID_QCD_DHOP_H
+#define GRID_QCD_DHOP_H

 namespace Grid {

-  namespace QCD {
+namespace QCD {

-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    // Helper routines that implement Wilson stencil for a single site.
-    // Common to both the WilsonFermion and WilsonFermion5D
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    class WilsonKernelsStatic { 
-    public:
-      // S-direction is INNERMOST and takes no part in the parity.
-      static int AsmOpt;  // these are a temporary hack
-      static int HandOpt; // these are a temporary hack
-    };
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Helper routines that implement Wilson stencil for a single site.
+// Common to both the WilsonFermion and WilsonFermion5D
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+class WilsonKernelsStatic {
+ public:
+  // S-direction is INNERMOST and takes no part in the parity.
+  static int AsmOpt;   // these are a temporary hack
+  static int HandOpt;  // these are a temporary hack
+};

-    template<class Impl> class WilsonKernels : public FermionOperator<Impl> , public WilsonKernelsStatic { 
-    public:
+template <class Impl>
+class WilsonKernels : public FermionOperator<Impl>, public WilsonKernelsStatic {
+ public:
+  INHERIT_IMPL_TYPES(Impl);
+  typedef FermionOperator<Impl> Base;

-     INHERIT_IMPL_TYPES(Impl);
-     typedef FermionOperator<Impl> Base;
-     
-    public:
-
-     void DiracOptDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
-			   std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-			   int sF, int sU,int Ls, int Ns, const FermionField &in, FermionField &out);
-      
-     void DiracOptDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
-			      std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-			      int sF,int sU,int Ls, int Ns, const FermionField &in,FermionField &out);
-
-     void DiracOptDhopDir(StencilImpl &st,DoubledGaugeField &U,
-			  std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-			  int sF,int sU,const FermionField &in, FermionField &out,int dirdisp,int gamma);
-
-    private:
-     // Specialised variants
-     void DiracOptGenericDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
-			   std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-			   int sF,int sU, const FermionField &in, FermionField &out);
-      
-     void DiracOptGenericDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
-			      std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-			      int sF,int sU,const FermionField &in,FermionField &out);
-
-     void DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
-			      std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-			      int sF,int sU,int Ls, int Ns, const FermionField &in, FermionField &out);
-
-
-     void DiracOptHandDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
-			      std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-			      int sF,int sU,const FermionField &in, FermionField &out);
-     
-     void DiracOptHandDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
-				 std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
-				 int sF,int sU,const FermionField &in, FermionField &out);
-    public:
-
-     WilsonKernels(const ImplParams &p= ImplParams());
-     
-    };
+ public:
+  template <bool EnableBool = true>
+  typename std::enable_if<Impl::Dimension == 3 && EnableBool, void>::type
+  DiracOptDhopSite(
+      StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+      std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
+      int sF, int sU, int Ls, int Ns, const FermionField &in,
+      FermionField &out) {
+#ifdef AVX512
+    if (AsmOpt) {
+      WilsonKernels<Impl>::DiracOptAsmDhopSite(st, lo, U, buf, sF, sU, Ls, Ns,
+                                               in, out);

+    } else {
+#else
+    {
+#endif
+      for (int site = 0; site < Ns; site++) {
+        for (int s = 0; s < Ls; s++) {
+          if (HandOpt)
+            WilsonKernels<Impl>::DiracOptHandDhopSite(st, lo, U, buf, sF, sU,
+                                                      in, out);
+          else
+            WilsonKernels<Impl>::DiracOptGenericDhopSite(st, lo, U, buf, sF, sU,
+                                                         in, out);
+          sF++;
+        }
+        sU++;
+      }
+    }
  }
+
+  template <bool EnableBool = true>
+  typename std::enable_if<Impl::Dimension != 3 && EnableBool, void>::type
+  DiracOptDhopSite(
+      StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+      std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
+      int sF, int sU, int Ls, int Ns, const FermionField &in,
+      FermionField &out) {
+    for (int site = 0; site < Ns; site++) {
+      for (int s = 0; s < Ls; s++) {
+        WilsonKernels<Impl>::DiracOptGenericDhopSite(st, lo, U, buf, sF, sU, in,
+                                                     out);
+        sF++;
+      }
+      sU++;
+    }
+  }
+
+  template <bool EnableBool = true>
+  typename std::enable_if<Impl::Dimension == 3 && EnableBool, void>::type
+  DiracOptDhopSiteDag(
+      StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+      std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
+      int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out) {
+    // No asm implementation yet.
+    //  if ( AsmOpt )
+    //  WilsonKernels<Impl>::DiracOptAsmDhopSiteDag(st,lo,U,buf,sF,sU,in,out);
+    //  else
+    for (int site = 0; site < Ns; site++) {
+      for (int s = 0; s < Ls; s++) {
+        if (HandOpt)
+          WilsonKernels<Impl>::DiracOptHandDhopSiteDag(st, lo, U, buf, sF, sU,
+                                                       in, out);
+        else
+          WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(st, lo, U, buf, sF,
+                                                          sU, in, out);
+        sF++;
+      }
+      sU++;
+    }
+  }
+
+  template <bool EnableBool = true>
+  typename std::enable_if<Impl::Dimension != 3 && EnableBool, void>::type
+  DiracOptDhopSiteDag(
+      StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+      std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
+      int sF, int sU, int Ls, int Ns, const FermionField &in, FermionField &out) {
+    for (int site = 0; site < Ns; site++) {
+      for (int s = 0; s < Ls; s++) {
+        WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(st, lo, U, buf, sF, sU,
+                                                        in, out);
+        sF++;
+      }
+      sU++;
+    }
+  }
+
+  void DiracOptDhopDir(
+      StencilImpl &st, DoubledGaugeField &U,
+      std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
+      int sF, int sU, const FermionField &in, FermionField &out, int dirdisp,
+      int gamma);
+
+ private:
+  // Specialised variants
+  void DiracOptGenericDhopSite(
+      StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+      std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
+      int sF, int sU, const FermionField &in, FermionField &out);
+
+  void DiracOptGenericDhopSiteDag(
+      StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+      std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
+      int sF, int sU, const FermionField &in, FermionField &out);
+
+  void DiracOptAsmDhopSite(
+      StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+      std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
+      int sF, int sU, int Ls, int Ns, const FermionField &in,
+      FermionField &out);
+
+  void DiracOptHandDhopSite(
+      StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+      std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
+      int sF, int sU, const FermionField &in, FermionField &out);
+
+  void DiracOptHandDhopSiteDag(
+      StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
+      std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
+      int sF, int sU, const FermionField &in, FermionField &out);
+
+ public:
+  WilsonKernels(const ImplParams &p = ImplParams());
+};
+}
 }
 #endif
--- a/lib/qcd/action/fermion/WilsonKernelsAsm.cc
+++ b/lib/qcd/action/fermion/WilsonKernelsAsm.cc
@@ -38,8 +38,8 @@ namespace QCD {
  ///////////////////////////////////////////////////////////
  // Default to no assembler implementation
  ///////////////////////////////////////////////////////////
-template<class Impl>
-void WilsonKernels<Impl >::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
+  template<class Impl>
+  void WilsonKernels<Impl>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
 					       std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
 					       int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
 {
--- a/lib/qcd/action/fermion/WilsonKernelsHand.cc
+++ b/lib/qcd/action/fermion/WilsonKernelsHand.cc
@@ -311,8 +311,8 @@ namespace Grid {
 namespace QCD {


-template<class Impl>
-void WilsonKernels<Impl >::DiracOptHandDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
+  template<class Impl>
+  void WilsonKernels<Impl>::DiracOptHandDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
 					       std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
 					       int ss,int sU,const FermionField &in, FermionField &out)
 {
@@ -554,8 +554,8 @@ void WilsonKernels<Impl >::DiracOptHandDhopSite(StencilImpl &st,LebesgueOrder &l
  }
 }

-template<class Impl>
-void WilsonKernels<Impl >::DiracOptHandDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
+  template<class Impl>
+  void WilsonKernels<Impl>::DiracOptHandDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
 					       std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  &buf,
 					       int ss,int sU,const FermionField &in, FermionField &out)
 {
--- a/lib/qcd/hmc/HmcRunner.h
+++ b/lib/qcd/hmc/HmcRunner.h
@@ -96,7 +96,7 @@ class NerscHmcRunnerTemplate {

    GridSerialRNG sRNG;
    GridParallelRNG pRNG(UGrid);
-    LatticeGaugeField U(UGrid);  // change this to an extended field (smearing class)
+    LatticeGaugeField U(UGrid);  // change this to an extended field (smearing class)?

    std::vector<int> SerSeed({1, 2, 3, 4, 5});
    std::vector<int> ParSeed({6, 7, 8, 9, 10});
--- a/lib/qcd/representations/hmc_types.h
+++ b/lib/qcd/representations/hmc_types.h
@@ -72,7 +72,7 @@ struct AccessTypes<A, TupleClass, 0, S...> {
  using elem = typename std::tuple_element<N, Rfields>::type;  // fields types

  typedef std::tuple<std::vector< A< elem<S> >* > ... > VectorCollection;
-  typedef std::tuple< A< elem<S> >* ... > ClassCollection;
+  typedef std::tuple< elem<S> ... > FieldTypeCollection;

  // Debug
  void return_size() {
--- a/tests/Test_hmc_WilsonAdjointFermionGauge.cc
+++ b/tests/Test_hmc_WilsonAdjointFermionGauge.cc
@@ -30,7 +30,6 @@ directory
 *************************************************************************************/
 /*  END LEGAL */
 #include "Grid.h"
-//#include "qcd/hmc/HmcRunner.h"

 using namespace std;
 using namespace Grid;
@@ -62,6 +61,7 @@ class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > {

    // temporarily need a gauge field
    LatticeGaugeField U(UGrid);
+    //AdjointRepresentation::LatticeField Ua(UGrid);

    // Gauge action
    WilsonGaugeActionR Waction(5.6);
@@ -69,7 +69,7 @@ class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > {
    Real mass = -0.77;
    FermionAction FermOp(U, *FGrid, *FrbGrid, mass);

-    ConjugateGradient<FermionField> CG(1.0e-8, 10000);
+    ConjugateGradient<FermionField> CG(1.0e-6, 10000);

    TwoFlavourPseudoFermionAction<ImplPolicy> Nf2(FermOp, CG, CG);