expClover support via helpers template class

Grid/qcd/action/fermion/CloverHelpers.h

@@ -0,0 +1,371 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./lib/qcd/action/fermion/WilsonCloverFermionImplementation.h
+
+    Copyright (C) 2017 - 2022
+
+    Author: paboyle <paboyle@ph.ed.ac.uk>
+    Author: Daniel Richtmann <daniel.richtmann@gmail.com>
+    Author: Mattia Bruno <mattia.bruno@cern.ch>
+
+    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
+
+#include <Grid/Grid.h>
+#include <Grid/qcd/spin/Dirac.h>
+#include <Grid/qcd/action/fermion/WilsonCloverHelpers.h>
+
+////////////////////////////////////////////
+// Standard Clover 
+//   (4+m0) + csw * clover_term
+// Exp Clover
+//   (4+m0) * exp(csw/(4+m0) clover_term)
+//   = (4+m0) + csw * clover_term + ...
+////////////////////////////////////////////
+
+NAMESPACE_BEGIN(Grid);
+
+
+//////////////////////////////////
+// Generic Standard Clover 
+//////////////////////////////////
+
+template<class Impl> 
+class CloverHelpers: public WilsonCloverHelpers<Impl> {
+public:
+
+  INHERIT_IMPL_TYPES(Impl);
+  INHERIT_CLOVER_TYPES(Impl);
+
+  typedef WilsonCloverHelpers<Impl> Helpers;
+
+  static void Instantiate(CloverField& CloverTerm, CloverField& CloverTermInv, RealD csw_t, RealD diag_mass) {
+    GridBase *grid = CloverTerm.Grid();
+    CloverTerm += diag_mass;
+   
+    int lvol = grid->lSites();
+    int DimRep = Impl::Dimension;
+    {
+      autoView(CTv,CloverTerm,CpuRead);
+      autoView(CTIv,CloverTermInv,CpuWrite);
+      thread_for(site, lvol, {
+        Coordinate lcoor;
+        grid->LocalIndexToLocalCoor(site, lcoor);
+        Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
+        Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
+        typename SiteClover::scalar_object Qx = Zero(), Qxinv = Zero();
+        peekLocalSite(Qx, CTv, lcoor);
+        //if (csw!=0){
+        for (int j = 0; j < Ns; j++)
+          for (int k = 0; k < Ns; k++)
+            for (int a = 0; a < DimRep; a++)
+              for (int b = 0; b < DimRep; b++){
+                auto zz =  Qx()(j, k)(a, b);
+                EigenCloverOp(a + j * DimRep, b + k * DimRep) = std::complex<double>(zz);
+              }
+        //   if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl;
+        
+        EigenInvCloverOp = EigenCloverOp.inverse();
+        //std::cout << EigenInvCloverOp << std::endl;
+        for (int j = 0; j < Ns; j++)
+          for (int k = 0; k < Ns; k++)
+            for (int a = 0; a < DimRep; a++)
+              for (int b = 0; b < DimRep; b++)
+                Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep);
+               //    if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl;
+               //  }
+               pokeLocalSite(Qxinv, CTIv, lcoor);
+      });
+    }
+  }
+
+    
+  static void CloverTermDerivative(GaugeField& clover_force, 
+                                   const FermionField& X, 
+                                   const FermionField& Y,
+                                   const std::vector<GaugeLinkField>& U,
+                                   RealD csw_t, RealD csw_r) {
+    GaugeLinkField force_mu(clover_force.Grid()), lambda(clover_force.Grid());
+    PropagatorField Lambda(clover_force.Grid());
+    
+    ///////////////////////////////////////////////////////////
+    // Clover term derivative
+    ///////////////////////////////////////////////////////////
+    Impl::outerProductImpl(Lambda, X, Y);
+    //std::cout << "Lambda:" << Lambda << std::endl;
+  
+    Gamma::Algebra sigma[] = {
+      Gamma::Algebra::SigmaXY,
+      Gamma::Algebra::SigmaXZ,
+      Gamma::Algebra::SigmaXT,
+      Gamma::Algebra::MinusSigmaXY,
+      Gamma::Algebra::SigmaYZ,
+      Gamma::Algebra::SigmaYT,
+      Gamma::Algebra::MinusSigmaXZ,
+      Gamma::Algebra::MinusSigmaYZ,
+      Gamma::Algebra::SigmaZT,
+      Gamma::Algebra::MinusSigmaXT,
+      Gamma::Algebra::MinusSigmaYT,
+      Gamma::Algebra::MinusSigmaZT};
+  
+    /*
+    sigma_{\mu \nu}=
+    | 0         sigma[0]  sigma[1]  sigma[2] |
+    | sigma[3]    0       sigma[4]  sigma[5] |
+    | sigma[6]  sigma[7]     0      sigma[8] |
+    | sigma[9]  sigma[10] sigma[11]   0      |
+    */
+  
+    int count = 0;
+    clover_force = Zero();
+    for (int mu = 0; mu < 4; mu++)
+    {
+      force_mu = Zero();
+      for (int nu = 0; nu < 4; nu++)
+      {
+        if (mu == nu)
+          continue;
+  
+        RealD factor;
+        if (nu == 4 || mu == 4)
+        {
+          factor = 2.0 * csw_t;
+        }
+        else
+        {
+          factor = 2.0 * csw_r;
+        }
+        PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
+        Impl::TraceSpinImpl(lambda, Slambda);                   // traceSpin ok
+        force_mu -= factor*Helpers::Cmunu(U, lambda, mu, nu);   // checked
+        count++;
+      }
+  
+      pokeLorentz(clover_force, U[mu] * force_mu, mu);
+    }
+  }
+};
+
+
+//////////////////////////////////
+// Generic Exp Clover 
+//////////////////////////////////
+
+template<class Impl> 
+class ExpCloverHelpers: public WilsonCloverHelpers<Impl> {
+public:
+
+  INHERIT_IMPL_TYPES(Impl);
+  INHERIT_CLOVER_TYPES(Impl);
+
+  template <typename vtype> using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;
+  typedef WilsonCloverHelpers<Impl> Helpers;
+
+  static void plusIdentity(const CloverField& in) {
+    int DimRep = Impl::Dimension;
+  
+    autoView(in_v, in, AcceleratorWrite);
+    
+    accelerator_for(ss, in.Grid()->oSites(), 1, {
+      for (int sa=0; sa<Ns; sa++)
+        for (int ca=0; ca<DimRep; ca++)
+          in_v[ss]()(sa,sa)(ca,ca) += 1.0;
+    });
+  }
+
+  static int getNMAX(RealD prec, RealD R) {
+    /* compute stop condition for exponential */
+    int NMAX=1;
+    RealD cond=R*R/2.;
+
+    while (cond*std::exp(R)>prec) {
+      NMAX++;
+      cond*=R/(double)(NMAX+1);
+    }
+    return NMAX;
+  }
+
+  static int getNMAX(Lattice<iImplClover<vComplexD>> &t, RealD R) {return getNMAX(1e-12,R);}
+  static int getNMAX(Lattice<iImplClover<vComplexF>> &t, RealD R) {return getNMAX(1e-6,R);}
+
+  static void Instantiate(CloverField& Clover, CloverField& CloverInv, RealD csw_t, RealD diag_mass) {
+    GridBase* grid = Clover.Grid();
+    CloverField ExpClover(grid);
+  
+    int NMAX = getNMAX(Clover, 3.*csw_t/diag_mass);
+    
+    // csw/(diag_mass) * clover
+    Clover *= (1.0/diag_mass);
+    
+    ExpClover = Clover;
+    plusIdentity(ExpClover); // 1 + Clover
+    
+    // Taylor expansion, slow but generic
+    RealD pref = 1.0;
+    for (int i=2; i<=NMAX; i++) {
+      Clover = Clover * Clover;
+      pref /= (RealD)(i);
+      ExpClover += pref * Clover;
+    }
+  
+    // Convert the data layout of the clover terms
+    Clover = ExpClover * diag_mass;
+    CloverInv = adj(ExpClover * (1.0/diag_mass));
+  }
+
+  static void CloverTermDerivative(GaugeField& clover_force, const FermionField& X, const FermionField& Y,
+                                   const std::vector<GaugeLinkField>& U, RealD csw_t, RealD csw_r) {
+    assert(0);
+  }
+};
+
+
+//////////////////////////////////
+// Compact Standard Clover 
+//////////////////////////////////
+
+
+template<class Impl> 
+class CompactCloverHelpers: public CompactWilsonCloverHelpers<Impl>,
+                            public WilsonCloverHelpers<Impl> {
+public:
+
+  INHERIT_IMPL_TYPES(Impl);
+  INHERIT_CLOVER_TYPES(Impl);
+  INHERIT_COMPACT_CLOVER_TYPES(Impl);
+
+  typedef WilsonCloverHelpers<Impl> Helpers;
+  typedef CompactWilsonCloverHelpers<Impl> CompactHelpers;
+  
+  static void MassTerm(CloverField& Clover, RealD diag_mass) {
+    Clover += diag_mass;
+  }
+
+  static void Instantiate(CloverDiagonalField& Diagonal,
+                          CloverTriangleField& Triangle,
+                          CloverDiagonalField& DiagonalInv,
+                          CloverTriangleField& TriangleInv,
+                          RealD csw_t, RealD diag_mass) {
+    // Invert the clover term in the improved layout
+    CompactHelpers::Invert(Diagonal, Triangle, DiagonalInv, TriangleInv);
+  }
+
+  // TODO: implement Cmunu for better performances with compact layout, but don't do it
+  // here, but rather in WilsonCloverHelpers.h -> CompactWilsonCloverHelpers
+  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu) {
+    return Helpers::Cmunu(U, lambda, mu, nu);
+  }
+};
+
+//////////////////////////////////
+// Compact Exp Clover 
+//////////////////////////////////
+
+
+template<class Impl>
+class CompactExpCloverHelpers: public CompactWilsonCloverHelpers<Impl> {
+public:
+
+  INHERIT_IMPL_TYPES(Impl);
+  INHERIT_CLOVER_TYPES(Impl);
+  INHERIT_COMPACT_CLOVER_TYPES(Impl);
+  
+  template <typename vtype> using iImplClover = iScalar<iMatrix<iMatrix<vtype, Impl::Dimension>, Ns>>;  
+  typedef CompactWilsonCloverHelpers<Impl> CompactHelpers;
+  
+  static void plusIdentity(const CloverField& in) {
+    int DimRep = Impl::Dimension;
+    
+    autoView(in_v, in, AcceleratorWrite);
+    
+    accelerator_for(ss, in.Grid()->oSites(), 1, {
+      for (int sa=0; sa<Ns; sa++)
+        for (int ca=0; ca<DimRep; ca++)
+          in_v[ss]()(sa,sa)(ca,ca) += 1.0;
+    });
+  }
+  
+  static int getNMAX(RealD prec, RealD R) {
+    /* compute stop condition for exponential */
+    int NMAX=1;
+    RealD cond=R*R/2.;
+
+    while (cond*std::exp(R)>prec) {
+      NMAX++;
+      cond*=R/(double)(NMAX+1);
+    }
+    return NMAX;
+  }
+
+  static int getNMAX(Lattice<iImplCloverDiagonal<vComplexD>> &t, RealD R) {return getNMAX(1e-12,R);}
+  static int getNMAX(Lattice<iImplCloverDiagonal<vComplexF>> &t, RealD R) {return getNMAX(1e-6,R);}
+  
+  static void MassTerm(CloverField& Clover, RealD diag_mass) {
+    // csw/(diag_mass) * clover
+    Clover = Clover * (1.0/diag_mass);
+  }
+  
+  static void Instantiate(CloverDiagonalField& Diagonal, CloverTriangleField& Triangle, 
+                          CloverDiagonalField& DiagonalInv, CloverTriangleField& TriangleInv,
+                          RealD csw_t, RealD diag_mass) {
+    GridBase* grid = Diagonal.Grid();
+    int NMAX = getNMAX(Diagonal, 3.*csw_t/diag_mass);
+    // To be optimized: too much memory traffic; implement exp in improved layout
+    // Felix + Fabian: replace code below with
+    //
+    // ConvertLayout Clover -> Diagonal, Triangle
+    // ModifyBoundaries
+    // EvaluateExp
+    
+    // code to be replaced
+    CloverField Clover(grid), ExpClover(grid);
+    
+    CompactHelpers::ConvertLayout(Diagonal, Triangle, Clover);
+    
+    ExpClover = Clover;
+    plusIdentity(ExpClover); // 1 + Clover
+    
+    RealD pref = 1.0;
+    for (int i=2; i<=NMAX; i++) {
+      Clover = Clover * Clover;
+      pref /= (RealD)(i);
+      ExpClover += pref * Clover;
+    }
+  
+    // Convert the data layout of the clover terms
+    CompactHelpers::ConvertLayout(ExpClover, Diagonal, Triangle);
+    Diagonal = Diagonal * diag_mass;
+    Triangle = Triangle * diag_mass;
+    CompactHelpers::ConvertLayout(adj(ExpClover), DiagonalInv, TriangleInv);
+    DiagonalInv = DiagonalInv*(1.0/diag_mass);
+    TriangleInv = TriangleInv*(1.0/diag_mass);
+  }
+
+  
+  static GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu) {
+    assert(0);
+  }
+  
+};
+
+
+NAMESPACE_END(Grid);

Grid/qcd/action/fermion/CompactWilsonCloverFermion.h

@@ -31,6 +31,7 @@
 
 #include <Grid/qcd/action/fermion/WilsonCloverTypes.h>
 #include <Grid/qcd/action/fermion/WilsonCloverHelpers.h>
+#include <Grid/qcd/action/fermion/CloverHelpers.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -85,7 +86,7 @@ NAMESPACE_BEGIN(Grid);
 //           + (2 * 1 + 4 * 1/2) triangle parts = 4 triangle parts =  60 complex words per site
 //                                                                 =  84 complex words per site
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
 class CompactWilsonCloverFermion : public WilsonFermion<Impl>,
                                    public WilsonCloverHelpers<Impl>,
                                    public CompactWilsonCloverHelpers<Impl> {

Grid/qcd/action/fermion/Fermion.h

@@ -138,38 +138,62 @@ typedef WilsonTMFermion<WilsonImplF> WilsonTMFermionF;
 typedef WilsonTMFermion<WilsonImplD> WilsonTMFermionD;
 
 // Clover fermions
-typedef WilsonCloverFermion<WilsonImplR> WilsonCloverFermionR;
+typedef CloverHelpers<WilsonImplR> CloverR;
-typedef WilsonCloverFermion<WilsonImplF> WilsonCloverFermionF;
+typedef CloverHelpers<WilsonImplF> CloverF;
-typedef WilsonCloverFermion<WilsonImplD> WilsonCloverFermionD;
+typedef CloverHelpers<WilsonImplD> CloverD;
 
-typedef WilsonCloverFermion<WilsonAdjImplR> WilsonCloverAdjFermionR;
+typedef ExpCloverHelpers<WilsonImplR> ExpCloverR;
-typedef WilsonCloverFermion<WilsonAdjImplF> WilsonCloverAdjFermionF;
+typedef ExpCloverHelpers<WilsonImplF> ExpCloverF;
-typedef WilsonCloverFermion<WilsonAdjImplD> WilsonCloverAdjFermionD;
+typedef ExpCloverHelpers<WilsonImplD> ExpCloverD;
 
-typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplR> WilsonCloverTwoIndexSymmetricFermionR;
+typedef WilsonCloverFermion<WilsonImplR, CloverR> WilsonCloverFermionR;
-typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplF> WilsonCloverTwoIndexSymmetricFermionF;
+typedef WilsonCloverFermion<WilsonImplF, CloverF> WilsonCloverFermionF;
-typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplD> WilsonCloverTwoIndexSymmetricFermionD;
+typedef WilsonCloverFermion<WilsonImplD, CloverD> WilsonCloverFermionD;
 
-typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplR> WilsonCloverTwoIndexAntiSymmetricFermionR;
+typedef WilsonCloverFermion<WilsonImplR, ExpCloverR> WilsonExpCloverFermionR;
-typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplF> WilsonCloverTwoIndexAntiSymmetricFermionF;
+typedef WilsonCloverFermion<WilsonImplF, ExpCloverF> WilsonExpCloverFermionF;
-typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplD> WilsonCloverTwoIndexAntiSymmetricFermionD;
+typedef WilsonCloverFermion<WilsonImplD, ExpCloverD> WilsonExpCloverFermionD;
+
+typedef WilsonCloverFermion<WilsonAdjImplR, CloverR> WilsonCloverAdjFermionR;
+typedef WilsonCloverFermion<WilsonAdjImplF, CloverF> WilsonCloverAdjFermionF;
+typedef WilsonCloverFermion<WilsonAdjImplD, CloverD> WilsonCloverAdjFermionD;
+
+typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplR, CloverR> WilsonCloverTwoIndexSymmetricFermionR;
+typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplF, CloverF> WilsonCloverTwoIndexSymmetricFermionF;
+typedef WilsonCloverFermion<WilsonTwoIndexSymmetricImplD, CloverD> WilsonCloverTwoIndexSymmetricFermionD;
+
+typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplR, CloverR> WilsonCloverTwoIndexAntiSymmetricFermionR;
+typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplF, CloverF> WilsonCloverTwoIndexAntiSymmetricFermionF;
+typedef WilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplD, CloverD> WilsonCloverTwoIndexAntiSymmetricFermionD;
 
 // Compact Clover fermions
-typedef CompactWilsonCloverFermion<WilsonImplR> CompactWilsonCloverFermionR;
+typedef CompactCloverHelpers<WilsonImplR> CompactCloverR;
-typedef CompactWilsonCloverFermion<WilsonImplF> CompactWilsonCloverFermionF;
+typedef CompactCloverHelpers<WilsonImplF> CompactCloverF;
-typedef CompactWilsonCloverFermion<WilsonImplD> CompactWilsonCloverFermionD;
+typedef CompactCloverHelpers<WilsonImplD> CompactCloverD;
 
-typedef CompactWilsonCloverFermion<WilsonAdjImplR> CompactWilsonCloverAdjFermionR;
+typedef CompactExpCloverHelpers<WilsonImplR> CompactExpCloverR;
-typedef CompactWilsonCloverFermion<WilsonAdjImplF> CompactWilsonCloverAdjFermionF;
+typedef CompactExpCloverHelpers<WilsonImplF> CompactExpCloverF;
-typedef CompactWilsonCloverFermion<WilsonAdjImplD> CompactWilsonCloverAdjFermionD;
+typedef CompactExpCloverHelpers<WilsonImplD> CompactExpCloverD;
 
-typedef CompactWilsonCloverFermion<WilsonTwoIndexSymmetricImplR> CompactWilsonCloverTwoIndexSymmetricFermionR;
+typedef CompactWilsonCloverFermion<WilsonImplR, CompactCloverR> CompactWilsonCloverFermionR;
-typedef CompactWilsonCloverFermion<WilsonTwoIndexSymmetricImplF> CompactWilsonCloverTwoIndexSymmetricFermionF;
+typedef CompactWilsonCloverFermion<WilsonImplF, CompactCloverF> CompactWilsonCloverFermionF;
-typedef CompactWilsonCloverFermion<WilsonTwoIndexSymmetricImplD> CompactWilsonCloverTwoIndexSymmetricFermionD;
+typedef CompactWilsonCloverFermion<WilsonImplD, CompactCloverD> CompactWilsonCloverFermionD;
 
-typedef CompactWilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplR> CompactWilsonCloverTwoIndexAntiSymmetricFermionR;
+typedef CompactWilsonCloverFermion<WilsonImplR, CompactExpCloverR> CompactWilsonExpCloverFermionR;
-typedef CompactWilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplF> CompactWilsonCloverTwoIndexAntiSymmetricFermionF;
+typedef CompactWilsonCloverFermion<WilsonImplF, CompactExpCloverF> CompactWilsonExpCloverFermionF;
-typedef CompactWilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplD> CompactWilsonCloverTwoIndexAntiSymmetricFermionD;
+typedef CompactWilsonCloverFermion<WilsonImplD, CompactExpCloverD> CompactWilsonExpCloverFermionD;
+
+typedef CompactWilsonCloverFermion<WilsonAdjImplR, CompactCloverR> CompactWilsonCloverAdjFermionR;
+typedef CompactWilsonCloverFermion<WilsonAdjImplF, CompactCloverF> CompactWilsonCloverAdjFermionF;
+typedef CompactWilsonCloverFermion<WilsonAdjImplD, CompactCloverD> CompactWilsonCloverAdjFermionD;
+
+typedef CompactWilsonCloverFermion<WilsonTwoIndexSymmetricImplR, CompactCloverR> CompactWilsonCloverTwoIndexSymmetricFermionR;
+typedef CompactWilsonCloverFermion<WilsonTwoIndexSymmetricImplF, CompactCloverF> CompactWilsonCloverTwoIndexSymmetricFermionF;
+typedef CompactWilsonCloverFermion<WilsonTwoIndexSymmetricImplD, CompactCloverD> CompactWilsonCloverTwoIndexSymmetricFermionD;
+
+typedef CompactWilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplR, CompactCloverR> CompactWilsonCloverTwoIndexAntiSymmetricFermionR;
+typedef CompactWilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplF, CompactCloverF> CompactWilsonCloverTwoIndexAntiSymmetricFermionF;
+typedef CompactWilsonCloverFermion<WilsonTwoIndexAntiSymmetricImplD, CompactCloverD> CompactWilsonCloverTwoIndexAntiSymmetricFermionD;
 
 // Domain Wall fermions
 typedef DomainWallFermion<WilsonImplR> DomainWallFermionR;

Grid/qcd/action/fermion/WilsonCloverFermion.h

@@ -32,6 +32,7 @@
 
 #include <Grid/qcd/action/fermion/WilsonCloverTypes.h>
 #include <Grid/qcd/action/fermion/WilsonCloverHelpers.h>
+#include <Grid/qcd/action/fermion/CloverHelpers.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -51,7 +52,7 @@ NAMESPACE_BEGIN(Grid);
 // csw_r = csw_t to recover the isotropic version
 //////////////////////////////////////////////////////////////////
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
 class WilsonCloverFermion : public WilsonFermion<Impl>,
                             public WilsonCloverHelpers<Impl>
 {

Grid/qcd/action/fermion/WilsonCloverHelpers.h

@@ -209,6 +209,8 @@ public:
 };
 
 
+////////////////////////////////////////////////////////
+
 template<class Impl> class CompactWilsonCloverHelpers {
 public:
 

Grid/qcd/action/fermion/implementation/CompactWilsonCloverFermionImplementation.h

@@ -32,17 +32,18 @@
 #include <Grid/qcd/spin/Dirac.h>
 #include <Grid/qcd/action/fermion/CompactWilsonCloverFermion.h>
 
+
 NAMESPACE_BEGIN(Grid);
-template<class Impl>
+template<class Impl, class CloverHelpers>
-CompactWilsonCloverFermion<Impl>::CompactWilsonCloverFermion(GaugeField& _Umu,
+CompactWilsonCloverFermion<Impl, CloverHelpers>::CompactWilsonCloverFermion(GaugeField& _Umu,
-                                                             GridCartesian& Fgrid,
+                                                                            GridCartesian& Fgrid,
-                                                             GridRedBlackCartesian& Hgrid,
+                                                                            GridRedBlackCartesian& Hgrid,
-                                                             const RealD _mass,
+                                                                            const RealD _mass,
-                                                             const RealD _csw_r,
+                                                                            const RealD _csw_r,
-                                                             const RealD _csw_t,
+                                                                            const RealD _csw_t,
-                                                             const RealD _cF,
+                                                                            const RealD _cF,
-                                                             const WilsonAnisotropyCoefficients& clover_anisotropy,
+                                                                            const WilsonAnisotropyCoefficients& clover_anisotropy,
-                                                             const ImplParams& impl_p)
+                                                                            const ImplParams& impl_p)
   : WilsonBase(_Umu, Fgrid, Hgrid, _mass, impl_p, clover_anisotropy)
   , csw_r(_csw_r)
   , csw_t(_csw_t)
@@ -68,40 +69,40 @@ CompactWilsonCloverFermion<Impl>::CompactWilsonCloverFermion(GaugeField& _Umu,
     CompactHelpers::SetupMasks(this->BoundaryMask, this->BoundaryMaskEven, this->BoundaryMaskOdd);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::Dhop(const FermionField& in, FermionField& out, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::Dhop(const FermionField& in, FermionField& out, int dag) {
   WilsonBase::Dhop(in, out, dag);
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::DhopOE(const FermionField& in, FermionField& out, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopOE(const FermionField& in, FermionField& out, int dag) {
   WilsonBase::DhopOE(in, out, dag);
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::DhopEO(const FermionField& in, FermionField& out, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopEO(const FermionField& in, FermionField& out, int dag) {
   WilsonBase::DhopEO(in, out, dag);
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::DhopDir(const FermionField& in, FermionField& out, int dir, int disp) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopDir(const FermionField& in, FermionField& out, int dir, int disp) {
   WilsonBase::DhopDir(in, out, dir, disp);
   if(this->open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::DhopDirAll(const FermionField& in, std::vector<FermionField>& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::DhopDirAll(const FermionField& in, std::vector<FermionField>& out) {
   WilsonBase::DhopDirAll(in, out);
   if(this->open_boundaries) {
     for(auto& o : out) ApplyBoundaryMask(o);
   }
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::M(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::M(const FermionField& in, FermionField& out) {
   out.Checkerboard() = in.Checkerboard();
   WilsonBase::Dhop(in, out, DaggerNo); // call base to save applying bc
   Mooee(in, Tmp);
@@ -109,8 +110,8 @@ void CompactWilsonCloverFermion<Impl>::M(const FermionField& in, FermionField& o
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::Mdag(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::Mdag(const FermionField& in, FermionField& out) {
   out.Checkerboard() = in.Checkerboard();
   WilsonBase::Dhop(in, out, DaggerYes);  // call base to save applying bc
   MooeeDag(in, Tmp);
@@ -118,20 +119,20 @@ void CompactWilsonCloverFermion<Impl>::Mdag(const FermionField& in, FermionField
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::Meooe(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::Meooe(const FermionField& in, FermionField& out) {
   WilsonBase::Meooe(in, out);
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MeooeDag(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MeooeDag(const FermionField& in, FermionField& out) {
   WilsonBase::MeooeDag(in, out);
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::Mooee(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::Mooee(const FermionField& in, FermionField& out) {
   if(in.Grid()->_isCheckerBoarded) {
     if(in.Checkerboard() == Odd) {
       MooeeInternal(in, out, DiagonalOdd, TriangleOdd);
@@ -144,13 +145,13 @@ void CompactWilsonCloverFermion<Impl>::Mooee(const FermionField& in, FermionFiel
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MooeeDag(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeDag(const FermionField& in, FermionField& out) {
   Mooee(in, out); // blocks are hermitian
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MooeeInv(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeInv(const FermionField& in, FermionField& out) {
   if(in.Grid()->_isCheckerBoarded) {
     if(in.Checkerboard() == Odd) {
       MooeeInternal(in, out, DiagonalInvOdd, TriangleInvOdd);
@@ -163,23 +164,23 @@ void CompactWilsonCloverFermion<Impl>::MooeeInv(const FermionField& in, FermionF
   if(open_boundaries) ApplyBoundaryMask(out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MooeeInvDag(const FermionField& in, FermionField& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeInvDag(const FermionField& in, FermionField& out) {
   MooeeInv(in, out); // blocks are hermitian
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::Mdir(const FermionField& in, FermionField& out, int dir, int disp) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::Mdir(const FermionField& in, FermionField& out, int dir, int disp) {
   DhopDir(in, out, dir, disp);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MdirAll(const FermionField& in, std::vector<FermionField>& out) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MdirAll(const FermionField& in, std::vector<FermionField>& out) {
   DhopDirAll(in, out);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MDeriv(GaugeField& force, const FermionField& X, const FermionField& Y, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MDeriv(GaugeField& force, const FermionField& X, const FermionField& Y, int dag) {
   assert(!open_boundaries); // TODO check for changes required for open bc
 
   // NOTE: code copied from original clover term
@@ -251,7 +252,7 @@ void CompactWilsonCloverFermion<Impl>::MDeriv(GaugeField& force, const FermionFi
       }
       PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
       Impl::TraceSpinImpl(lambda, Slambda);                   // traceSpin ok
-      force_mu -= factor*Helpers::Cmunu(U, lambda, mu, nu);   // checked
+      force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu);   // checked
       count++;
     }
 
@@ -261,18 +262,18 @@ void CompactWilsonCloverFermion<Impl>::MDeriv(GaugeField& force, const FermionFi
   force += clover_force;
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MooDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
   assert(0);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MeeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MeeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
   assert(0);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::MooeeInternal(const FermionField&        in,
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::MooeeInternal(const FermionField&        in,
                     FermionField&              out,
                     const CloverDiagonalField& diagonal,
                     const CloverTriangleField& triangle) {
@@ -285,8 +286,8 @@ void CompactWilsonCloverFermion<Impl>::MooeeInternal(const FermionField&
   CompactHelpers::MooeeKernel(diagonal.oSites(), 1, in, out, diagonal, triangle);
 }
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-void CompactWilsonCloverFermion<Impl>::ImportGauge(const GaugeField& _Umu) {
+void CompactWilsonCloverFermion<Impl, CloverHelpers>::ImportGauge(const GaugeField& _Umu) {
   // NOTE: parts copied from original implementation
 
   // Import gauge into base class
@@ -318,7 +319,8 @@ void CompactWilsonCloverFermion<Impl>::ImportGauge(const GaugeField& _Umu) {
   TmpOriginal += Helpers::fillCloverXT(Ex) * csw_t;
   TmpOriginal += Helpers::fillCloverYT(Ey) * csw_t;
   TmpOriginal += Helpers::fillCloverZT(Ez) * csw_t;
-  TmpOriginal += this->diag_mass;
+  // Handle mass term based on clover policy
+  CloverHelpers::MassTerm(TmpOriginal, this->diag_mass);
   
   // Convert the data layout of the clover term
   double t4 = usecond();
@@ -328,9 +330,9 @@ void CompactWilsonCloverFermion<Impl>::ImportGauge(const GaugeField& _Umu) {
   double t5 = usecond();
   if(open_boundaries) CompactHelpers::ModifyBoundaries(Diagonal, Triangle, csw_t, cF, this->diag_mass);
 
-  // Invert the clover term in the improved layout
+  // Instantiate based on clover policy
   double t6 = usecond();
-  CompactHelpers::Invert(Diagonal, Triangle, DiagonalInv, TriangleInv);
+  CloverHelpers::Instantiate(Diagonal, Triangle, DiagonalInv, TriangleInv, csw_t, this->diag_mass);
 
   // Fill the remaining clover fields
   double t7 = usecond();
@@ -353,7 +355,7 @@ void CompactWilsonCloverFermion<Impl>::ImportGauge(const GaugeField& _Umu) {
             << ", fill clover = "               << (t4 - t3) / 1e6
             << ", convert = "                   << (t5 - t4) / 1e6
             << ", boundaries = "                << (t6 - t5) / 1e6
-            << ", inversions = "                << (t7 - t6) / 1e6
+            << ", instantiate = "               << (t7 - t6) / 1e6
             << ", pick cbs = "                  << (t8 - t7) / 1e6
             << ", total = "                     << (t8 - t0) / 1e6
             << std::endl;

Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h

@@ -34,8 +34,8 @@
 
 NAMESPACE_BEGIN(Grid);
 
-template<class Impl>
+template<class Impl, class CloverHelpers>
-WilsonCloverFermion<Impl>::WilsonCloverFermion(GaugeField&                         _Umu,
+WilsonCloverFermion<Impl, CloverHelpers>::WilsonCloverFermion(GaugeField&                         _Umu,
                                                GridCartesian&                      Fgrid,
                                                GridRedBlackCartesian&              Hgrid,
                                                const RealD                         _mass,
@@ -74,8 +74,8 @@ WilsonCloverFermion<Impl>::WilsonCloverFermion(GaugeField&
 }
 
 // *NOT* EO
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::M(const FermionField &in, FermionField &out)
 {
   FermionField temp(out.Grid());
 
@@ -89,8 +89,8 @@ void WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
   out += temp;
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::Mdag(const FermionField &in, FermionField &out)
 {
   FermionField temp(out.Grid());
 
@@ -104,8 +104,8 @@ void WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
   out += temp;
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
+void WilsonCloverFermion<Impl, CloverHelpers>::ImportGauge(const GaugeField &_Umu)
 {
   double t0 = usecond();
   WilsonFermion<Impl>::ImportGauge(_Umu);
@@ -131,47 +131,50 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
   CloverTerm += Helpers::fillCloverXT(Ex) * csw_t;
   CloverTerm += Helpers::fillCloverYT(Ey) * csw_t;
   CloverTerm += Helpers::fillCloverZT(Ez) * csw_t;
-  CloverTerm += diag_mass;
    
   double t4 = usecond();
-  int lvol = _Umu.Grid()->lSites();
+  CloverHelpers::Instantiate(CloverTerm, CloverTermInv, csw_t, this->diag_mass);
-  int DimRep = Impl::Dimension;
+//   CloverTerm += diag_mass;
+// 
+//   double t4 = usecond();
+//   int lvol = _Umu.Grid()->lSites();
+//   int DimRep = Impl::Dimension;
+// 
+//   double t5 = usecond();
+//   {
+//     autoView(CTv,CloverTerm,CpuRead);
+//     autoView(CTIv,CloverTermInv,CpuWrite);
+//     thread_for(site, lvol, {
+//       Coordinate lcoor;
+//       grid->LocalIndexToLocalCoor(site, lcoor);
+//       Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
+//       Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
+//       typename SiteClover::scalar_object Qx = Zero(), Qxinv = Zero();
+//       peekLocalSite(Qx, CTv, lcoor);
+//       //if (csw!=0){
+//       for (int j = 0; j < Ns; j++)
+// 	for (int k = 0; k < Ns; k++)
+// 	  for (int a = 0; a < DimRep; a++)
+// 	    for (int b = 0; b < DimRep; b++){
+// 	      auto zz =  Qx()(j, k)(a, b);
+// 	      EigenCloverOp(a + j * DimRep, b + k * DimRep) = std::complex<double>(zz);
+// 	    }
+//       //   if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl;
+//       
+//       EigenInvCloverOp = EigenCloverOp.inverse();
+//       //std::cout << EigenInvCloverOp << std::endl;
+//       for (int j = 0; j < Ns; j++)
+// 	for (int k = 0; k < Ns; k++)
+// 	  for (int a = 0; a < DimRep; a++)
+// 	    for (int b = 0; b < DimRep; b++)
+// 	      Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep);
+//       //    if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl;
+//       //  }
+//       pokeLocalSite(Qxinv, CTIv, lcoor);
+//     });
+//   }
 
   double t5 = usecond();
-  {
-    autoView(CTv,CloverTerm,CpuRead);
-    autoView(CTIv,CloverTermInv,CpuWrite);
-    thread_for(site, lvol, {
-      Coordinate lcoor;
-      grid->LocalIndexToLocalCoor(site, lcoor);
-      Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
-      Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep);
-      typename SiteClover::scalar_object Qx = Zero(), Qxinv = Zero();
-      peekLocalSite(Qx, CTv, lcoor);
-      //if (csw!=0){
-      for (int j = 0; j < Ns; j++)
-	for (int k = 0; k < Ns; k++)
-	  for (int a = 0; a < DimRep; a++)
-	    for (int b = 0; b < DimRep; b++){
-	      auto zz =  Qx()(j, k)(a, b);
-	      EigenCloverOp(a + j * DimRep, b + k * DimRep) = std::complex<double>(zz);
-	    }
-      //   if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl;
-      
-      EigenInvCloverOp = EigenCloverOp.inverse();
-      //std::cout << EigenInvCloverOp << std::endl;
-      for (int j = 0; j < Ns; j++)
-	for (int k = 0; k < Ns; k++)
-	  for (int a = 0; a < DimRep; a++)
-	    for (int b = 0; b < DimRep; b++)
-	      Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep);
-      //    if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl;
-      //  }
-      pokeLocalSite(Qxinv, CTIv, lcoor);
-    });
-  }
-
-  double t6 = usecond();
   // Separate the even and odd parts
   pickCheckerboard(Even, CloverTermEven, CloverTerm);
   pickCheckerboard(Odd, CloverTermOdd, CloverTerm);
@@ -184,7 +187,7 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
 
   pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv));
   pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv));
-  double t7 = usecond();
+  double t6 = usecond();
 
 #if 0
   std::cout << GridLogMessage << "WilsonCloverFermion::ImportGauge timings:"
@@ -192,40 +195,39 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
             << ", allocations = "               << (t2 - t1) / 1e6
             << ", field strength = "            << (t3 - t2) / 1e6
             << ", fill clover = "               << (t4 - t3) / 1e6
-            << ", misc = "                      << (t5 - t4) / 1e6
+            << ", instantiation = "             << (t5 - t4) / 1e6
-            << ", inversions = "                << (t6 - t5) / 1e6
+            << ", pick cbs = "                  << (t6 - t5) / 1e6
-            << ", pick cbs = "                  << (t7 - t6) / 1e6
+            << ", total = "                     << (t6 - t0) / 1e6
-            << ", total = "                     << (t7 - t0) / 1e6
             << std::endl;
 #endif
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::Mooee(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::Mooee(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerNo, InverseNo);
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::MooeeDag(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerYes, InverseNo);
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::MooeeInv(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerNo, InverseYes);
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out)
+void WilsonCloverFermion<Impl, CloverHelpers>::MooeeInvDag(const FermionField &in, FermionField &out)
 {
   this->MooeeInternal(in, out, DaggerYes, InverseYes);
 }
 
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv)
+void WilsonCloverFermion<Impl, CloverHelpers>::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv)
 {
   out.Checkerboard() = in.Checkerboard();
   CloverField *Clover;
@@ -278,8 +280,8 @@ void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionFie
 } // MooeeInternal
 
 // Derivative parts unpreconditioned pseudofermions
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag)
+void WilsonCloverFermion<Impl, CloverHelpers>::MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag)
 {
   conformable(X.Grid(), Y.Grid());
   conformable(X.Grid(), force.Grid());
@@ -349,7 +351,7 @@ void WilsonCloverFermion<Impl>::MDeriv(GaugeField &force, const FermionField &X,
       }
       PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
       Impl::TraceSpinImpl(lambda, Slambda);                   // traceSpin ok
-      force_mu -= factor*Helpers::Cmunu(U, lambda, mu, nu);                   // checked
+      force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu);                   // checked
       count++;
     }
 
@@ -360,15 +362,15 @@ void WilsonCloverFermion<Impl>::MDeriv(GaugeField &force, const FermionField &X,
 }
 
 // Derivative parts
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag)
+void WilsonCloverFermion<Impl, CloverHelpers>::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag)
 {
   assert(0);
 }
 
 // Derivative parts
-template <class Impl>
+template<class Impl, class CloverHelpers>
-void WilsonCloverFermion<Impl>::MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
+void WilsonCloverFermion<Impl, CloverHelpers>::MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
 {
   assert(0); // not implemented yet
 }

Grid/qcd/action/fermion/instantiation/CompactWilsonCloverFermionInstantiation.cc.master

@@ -9,6 +9,7 @@
     Author: paboyle <paboyle@ph.ed.ac.uk>
     Author: Guido Cossu <guido.cossu@ed.ac.uk>
     Author: Daniel Richtmann <daniel.richtmann@gmail.com>
+    Author: Mattia Bruno <mattia.bruno@cern.ch>
 
     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
@@ -32,10 +33,12 @@
 #include <Grid/qcd/spin/Dirac.h>
 #include <Grid/qcd/action/fermion/CompactWilsonCloverFermion.h>
 #include <Grid/qcd/action/fermion/implementation/CompactWilsonCloverFermionImplementation.h>
+#include <Grid/qcd/action/fermion/CloverHelpers.h>
 
 NAMESPACE_BEGIN(Grid);
 
 #include "impl.h"
-template class CompactWilsonCloverFermion<IMPLEMENTATION>; 
+template class CompactWilsonCloverFermion<IMPLEMENTATION, CompactCloverHelpers<IMPLEMENTATION>>; 
+template class CompactWilsonCloverFermion<IMPLEMENTATION, CompactExpCloverHelpers<IMPLEMENTATION>>; 
 
 NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/WilsonCloverFermionInstantiation.cc.master

@@ -8,6 +8,7 @@
 
     Author: paboyle <paboyle@ph.ed.ac.uk>
     Author: Guido Cossu <guido.cossu@ed.ac.uk>
+    Author: Mattia Bruno <mattia.bruno@cern.ch>
     
     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
@@ -31,10 +32,12 @@
 #include <Grid/qcd/spin/Dirac.h>
 #include <Grid/qcd/action/fermion/WilsonCloverFermion.h>
 #include <Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h>
+#include <Grid/qcd/action/fermion/CloverHelpers.h>
 
 NAMESPACE_BEGIN(Grid);
 
 #include "impl.h"
-template class WilsonCloverFermion<IMPLEMENTATION>; 
+template class WilsonCloverFermion<IMPLEMENTATION, CloverHelpers<IMPLEMENTATION>>; 
+template class WilsonCloverFermion<IMPLEMENTATION, ExpCloverHelpers<IMPLEMENTATION>>; 
 
 NAMESPACE_END(Grid);