Added representations definitions for the HMC

lib/qcd/QCD.h

@@ -484,17 +484,23 @@ namespace QCD {
 }   //namespace QCD
 } // Grid
 
-#include <qcd/utils/SpaceTimeGrid.h>
 #include <qcd/spin/Dirac.h>
+
 #include <qcd/spin/TwoSpinor.h>
+
+#include <qcd/utils/SpaceTimeGrid.h>
 #include <qcd/utils/LinalgUtils.h>
 #include <qcd/utils/CovariantCshift.h>
 #include <qcd/utils/SUn.h>
+#include <qcd/utils/SUnAdjoint.h>
+
 #include <qcd/action/Actions.h>
+
 #include <qcd/hmc/integrators/Integrator.h>
 #include <qcd/hmc/integrators/Integrator_algorithm.h>
 #include <qcd/hmc/HMC.h>
 
 #include <qcd/smearing/Smearing.h>
 
+#include <qcd/representations/hmc_types.h>
 #endif

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -116,6 +116,7 @@ class WilsonImpl
     : public PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > {
  public:
   typedef PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > Gimpl;
+  constexpr bool is_fundamental() const{return Nrepresentation == Nc ? 1 : 0;}
 
   INHERIT_GIMPL_TYPES(Gimpl);
 
@@ -504,6 +505,10 @@ typedef WilsonImpl<vComplex, Nc> WilsonImplR;   // Real.. whichever prec
 typedef WilsonImpl<vComplexF, Nc> WilsonImplF;  // Float
 typedef WilsonImpl<vComplexD, Nc> WilsonImplD;  // Double
 
+typedef WilsonImpl<vComplex,  SU_Adjoint<Nc>::Dimension > WilsonAdjImplR;   // Real.. whichever prec
+typedef WilsonImpl<vComplexF, SU_Adjoint<Nc>::Dimension > WilsonAdjImplF;  // Float
+typedef WilsonImpl<vComplexD, SU_Adjoint<Nc>::Dimension > WilsonAdjImplD;  // Double
+
 typedef DomainWallRedBlack5dImpl<vComplex, Nc>
     DomainWallRedBlack5dImplR;  // Real.. whichever prec
 typedef DomainWallRedBlack5dImpl<vComplexF, Nc>

lib/qcd/action/gauge/GaugeImpl.h

@@ -30,7 +30,6 @@ directory
 #define GRID_QCD_GAUGE_IMPL_H
 
 namespace Grid {
-
 namespace QCD {
 
 ////////////////////////////////////////////////////////////////////////
@@ -64,7 +63,7 @@ public:
                                                  // ugly
   typedef Lattice<SiteGaugeField> GaugeField;
 
-  // Move this elsewhere?
+  // Move this elsewhere? FIXME
   static inline void AddGaugeLink(GaugeField &U, GaugeLinkField &W,
                                   int mu) { // U[mu] += W
     PARALLEL_FOR_LOOP
@@ -174,12 +173,19 @@ typedef GaugeImplTypes<vComplex, Nc> GimplTypesR;
 typedef GaugeImplTypes<vComplexF, Nc> GimplTypesF;
 typedef GaugeImplTypes<vComplexD, Nc> GimplTypesD;
 
+typedef GaugeImplTypes<vComplex, SU<Nc>::AdjointDimension> GimplAdjointTypesR;
+typedef GaugeImplTypes<vComplexF, SU<Nc>::AdjointDimension> GimplAdjointTypesF;
+typedef GaugeImplTypes<vComplexD, SU<Nc>::AdjointDimension> GimplAdjointTypesD;
+
 typedef PeriodicGaugeImpl<GimplTypesR> PeriodicGimplR; // Real.. whichever prec
 typedef PeriodicGaugeImpl<GimplTypesF> PeriodicGimplF; // Float
 typedef PeriodicGaugeImpl<GimplTypesD> PeriodicGimplD; // Double
 
-typedef ConjugateGaugeImpl<GimplTypesR>
+typedef PeriodicGaugeImpl<GimplAdjointTypesR> PeriodicGimplAdjR; // Real.. whichever prec
-    ConjugateGimplR; // Real.. whichever prec
+typedef PeriodicGaugeImpl<GimplAdjointTypesF> PeriodicGimplAdjF; // Float
+typedef PeriodicGaugeImpl<GimplAdjointTypesD> PeriodicGimplAdjD; // Double
+
+typedef ConjugateGaugeImpl<GimplTypesR> ConjugateGimplR; // Real.. whichever prec
 typedef ConjugateGaugeImpl<GimplTypesF> ConjugateGimplF; // Float
 typedef ConjugateGaugeImpl<GimplTypesD> ConjugateGimplD; // Double
 }

lib/qcd/action/pseudofermion/TwoFlavour.h

@@ -23,7 +23,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 with this program; if not, write to the Free Software Foundation, Inc.,
 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
+See the full license in the file "LICENSE" in the top level distribution
+directory
 *************************************************************************************/
 /*  END LEGAL */
 #ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_H
@@ -41,7 +42,6 @@ namespace Grid{
   INHERIT_IMPL_TYPES(Impl);
 
  private:
-      
   FermionOperator<Impl> &FermOp;  // the basic operator
 
   OperatorFunction<FermionField> &DerivativeSolver;
@@ -56,15 +56,17 @@ namespace Grid{
   /////////////////////////////////////////////////
   TwoFlavourPseudoFermionAction(FermionOperator<Impl> &Op,
                                 OperatorFunction<FermionField> &DS,
-				  OperatorFunction<FermionField> & AS
+                                OperatorFunction<FermionField> &AS)
-				  ) : FermOp(Op), DerivativeSolver(DS), ActionSolver(AS), Phi(Op.FermionGrid()) {
+      : FermOp(Op),
-      };
+        DerivativeSolver(DS),
+        ActionSolver(AS),
+        Phi(Op.FermionGrid()){};
 
   //////////////////////////////////////////////////////////////////////////////////////
-      // Push the gauge field in to the dops. Assume any BC's and smearing already applied
+  // Push the gauge field in to the dops. Assume any BC's and smearing already
+  // applied
   //////////////////////////////////////////////////////////////////////////////////////
   virtual void refresh(const GaugeField &U, GridParallelRNG &pRNG) {
-
     // P(phi) = e^{- phi^dag (MdagM)^-1 phi}
     // Phi = Mdag eta
     // P(eta) = e^{- eta^dag eta}
@@ -76,7 +78,8 @@ namespace Grid{
     //
     // Chroma has this scale factor: two_flavor_monomial_w.h
     // IroIro: does not use this scale. It is absorbed by a change of vars
-	//         in the Phi integral, and thus is only an irrelevant prefactor for the partition function.
+    //         in the Phi integral, and thus is only an irrelevant prefactor for
+    //         the partition function.
     //
     RealD scale = std::sqrt(0.5);
     FermionField eta(FermOp.FermionGrid());
@@ -87,14 +90,12 @@ namespace Grid{
     FermOp.Mdag(eta, Phi);
 
     Phi = Phi * scale;
-	
   };
 
   //////////////////////////////////////////////////////
   // S = phi^dag (Mdag M)^-1 phi
   //////////////////////////////////////////////////////
   virtual RealD S(const GaugeField &U) {
-
     FermOp.ImportGauge(U);
 
     FermionField X(FermOp.FermionGrid());
@@ -106,19 +107,20 @@ namespace Grid{
     MdagMOp.Op(X, Y);
 
     RealD action = norm2(Y);
-	std::cout << GridLogMessage << "Pseudofermion action "<<action<<std::endl;
+    std::cout << GridLogMessage << "Pseudofermion action " << action
+              << std::endl;
     return action;
   };
 
   //////////////////////////////////////////////////////
   // dS/du = - phi^dag  (Mdag M)^-1 [ Mdag dM + dMdag M ]  (Mdag M)^-1 phi
-      //       = - phi^dag M^-1 dM (MdagM)^-1 phi -  phi^dag (MdagM)^-1 dMdag dM (Mdag)^-1 phi 
+  //       = - phi^dag M^-1 dM (MdagM)^-1 phi -  phi^dag (MdagM)^-1 dMdag dM
+  //       (Mdag)^-1 phi
   //
   //       = - Ydag dM X  - Xdag dMdag Y
   //
   //////////////////////////////////////////////////////
   virtual void deriv(const GaugeField &U, GaugeField &dSdU) {
-
     FermOp.ImportGauge(U);
 
     FermionField X(FermOp.FermionGrid());
@@ -131,18 +133,18 @@ namespace Grid{
     DerivativeSolver(MdagMOp, Phi, X);
     MdagMOp.Op(X, Y);
 
-	// Our conventions really make this UdSdU; We do not differentiate wrt Udag here.
+    // Our conventions really make this UdSdU; We do not differentiate wrt Udag
+    // here.
     // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt.
 
-	FermOp.MDeriv(tmp , Y, X,DaggerNo );  dSdU=tmp;
+    FermOp.MDeriv(tmp, Y, X, DaggerNo);
-	FermOp.MDeriv(tmp , X, Y,DaggerYes);  dSdU=dSdU+tmp;
+    dSdU = tmp;
-	
+    FermOp.MDeriv(tmp, X, Y, DaggerYes);
-	//dSdU = Ta(dSdU);
+    dSdU = dSdU + tmp;
 
+    // not taking here the traceless antihermitian component
   };
-
 };
-    
 }
 }
 

lib/qcd/hmc/HmcRunner.h

@@ -22,7 +22,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 with this program; if not, write to the Free Software Foundation, Inc.,
 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
+See the full license in the file "LICENSE" in the top level distribution
+directory
 *************************************************************************************/
 /*  END LEGAL */
 #ifndef HMC_RUNNER
@@ -31,11 +32,9 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 namespace Grid {
 namespace QCD {
 
-
 template <class Gimpl>
 class NerscHmcRunnerTemplate {
  public:
-
   INHERIT_GIMPL_TYPES(Gimpl);
 
   enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart };
@@ -49,20 +48,26 @@ public:
 
   virtual void BuildTheAction(int argc, char **argv) = 0;  // necessary?
 
-  
   void Run(int argc, char **argv) {
-
     StartType_t StartType = HotStart;
 
     std::string arg;
 
     if (GridCmdOptionExists(argv, argv + argc, "--StartType")) {
       arg = GridCmdOptionPayload(argv, argv + argc, "--StartType");
-      if ( arg == "HotStart" ) { StartType = HotStart; }
+      if (arg == "HotStart") {
-      else if ( arg == "ColdStart" ) { StartType = ColdStart; }
+        StartType = HotStart;
-      else if ( arg == "TepidStart" ) { StartType = TepidStart; }
+      } else if (arg == "ColdStart") {
-      else if ( arg == "CheckpointStart" ) { StartType = CheckpointStart; }
+        StartType = ColdStart;
-      else assert(0);
+      } else if (arg == "TepidStart") {
+        StartType = TepidStart;
+      } else if (arg == "CheckpointStart") {
+        StartType = CheckpointStart;
+      } else {
+        std::cout << GridLogError << "Unrecognized option in --StartType\n";
+        std::cout << GridLogError << "Valid [HotStart, ColdStart, TepidStart, CheckpointStart]\n";
+        assert(0);
+      }
     }
 
     int StartTraj = 0;
@@ -89,7 +94,6 @@ public:
       NumThermalizations = ivec[0];
     }
 
-
     GridSerialRNG sRNG;
     GridParallelRNG pRNG(UGrid);
     LatticeGaugeField U(UGrid);  // change this to an extended field (smearing class)
@@ -97,7 +101,6 @@ public:
     std::vector<int> SerSeed({1, 2, 3, 4, 5});
     std::vector<int> ParSeed({6, 7, 8, 9, 10});
 
-    
     // Create integrator, including the smearing policy
     // Smearing policy
     std::cout << GridLogDebug << " Creating the Stout class\n";
@@ -108,13 +111,14 @@ public:
     SmearedConfiguration<Gimpl> SmearingPolicy(UGrid, Nsmear, Stout);
     std::cout << GridLogDebug << " done\n";
     //////////////
-    typedef MinimumNorm2<GaugeField, SmearedConfiguration<Gimpl> >  IntegratorType;// change here to change the algorithm
+    typedef MinimumNorm2<GaugeField, SmearedConfiguration<Gimpl> >
+        IntegratorType;  // change here to change the algorithm
     IntegratorParameters MDpar(20);
     IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);
 
-    
     // Checkpoint strategy
-    NerscHmcCheckpointer<Gimpl> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
+    NerscHmcCheckpointer<Gimpl> Checkpoint(std::string("ckpoint_lat"),
+                                           std::string("ckpoint_rng"), 1);
     PlaquetteLogger<Gimpl> PlaqLog(std::string("plaq"));
 
     HMCparameters HMCpar;
@@ -122,7 +126,6 @@ public:
     HMCpar.Trajectories = NumTraj;
     HMCpar.NoMetropolisUntil = NumThermalizations;
 
-
     if (StartType == HotStart) {
       // Hot start
       HMCpar.MetropolisTest = true;
@@ -147,20 +150,20 @@ public:
       Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
     }
 
-    // Attach the gauge field to the smearing Policy and create the fill the smeared set
+    // Attach the gauge field to the smearing Policy and create the fill the
+    // smeared set
     // notice that the unit configuration is singular in this procedure
     std::cout << GridLogMessage << "Filling the smeared set\n";
     SmearingPolicy.set_GaugeField(U);
 
-    HybridMonteCarlo<GaugeField,IntegratorType>  HMC(HMCpar, MDynamics,sRNG,pRNG,U); 
+    HybridMonteCarlo<GaugeField, IntegratorType> HMC(HMCpar, MDynamics, sRNG,
+                                                     pRNG, U);
     HMC.AddObservable(&Checkpoint);
     HMC.AddObservable(&PlaqLog);
 
     // Run it
     HMC.evolve();
-    
   }
-  
 };
 
 typedef NerscHmcRunnerTemplate<PeriodicGimplR> NerscHmcRunner;
@@ -174,6 +177,6 @@ public:
 typedef NerscHmcRunnerTemplate<ConjugateGimplR> ConjugateNerscHmcRunner;
 typedef NerscHmcRunnerTemplate<ConjugateGimplF> ConjugateNerscHmcRunnerF;
 typedef NerscHmcRunnerTemplate<ConjugateGimplD> ConjugateNerscHmcRunnerD;
-
+}
-}}
+}
 #endif

lib/qcd/hmc/integrators/Integrator.h

@@ -138,8 +138,6 @@ class Integrator {
   }
   void update_U(GaugeField& Mom, GaugeField& U, double ep) {
     // rewrite exponential to deal automatically  with the lorentz index?
-    //	GaugeLinkField Umu(U._grid);
-    //	GaugeLinkField Pmu(U._grid);
     for (int mu = 0; mu < Nd; mu++) {
       auto Umu = PeekIndex<LorentzIndex>(U, mu);
       auto Pmu = PeekIndex<LorentzIndex>(Mom, mu);
@@ -168,6 +166,9 @@ class Integrator {
   void refresh(GaugeField& U, GridParallelRNG& pRNG) {
     std::cout << GridLogIntegrator << "Integrator refresh\n";
     generate_momenta(P, pRNG);
+    // The Smearer is attached to a pointer of the gauge field
+    // automatically gets the updated field
+    // whether or not has been accepted in the previous sweep
     for (int level = 0; level < as.size(); ++level) {
       for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
         // get gauge field from the SmearingPolicy and

lib/qcd/representations/adjoint.h

@@ -0,0 +1,61 @@
+/*
+ *	Policy classes for the HMC
+ *	Author: Guido Cossu
+*/	
+
+#ifndef ADJOINT_H
+#define ADJOINT_H
+
+
+namespace Grid {
+namespace QCD {
+
+/*
+* This is an helper class for the HMC
+* Should contain only the data for the adjoint representation
+* and the facility to convert from the fundamental -> adjoint
+*/
+
+template <int ncolour>
+class AdjointRep {
+ public:
+  typename SU_Adjoint<ncolour>::LatticeAdjMatrix U;
+  const int Dimension = ncolour * ncolour - 1;
+
+  explicit AdjointRep(GridBase* grid):U(grid) {}
+  void update_representation(const LatticeGaugeField& Uin) {
+    // Uin is in the fundamental representation
+    // get the U in AdjointRep
+    // (U_adj)_B = tr[e^a U e^b U^dag]
+    // e^a = t^a/sqrt(T_F)
+    // where t^a is the generator in the fundamental
+    // T_F is 1/2 for the fundamental representation
+    conformable(U, Uin);
+    U = zero;
+    LatticeGaugeField tmp(Uin._grid);
+
+    Vector<typename SU<ncolour>::Matrix > ta(ncolour * ncolour - 1);
+
+    // FIXME probably not very efficient to get all the generators everytime
+    for (int a = 0; a < Dimension; a++) SU<ncolour>::generator(a, ta[a]);
+
+    for (int a = 0; a < Dimension; a++) {
+    	tmp = 2.0 * adj(Uin) * ta[a] * Uin;
+      for (int b = 0; b < (ncolour * ncolour - 1); b++) {
+        auto Tr = TensorRemove(trace(tmp * ta[b]));
+        pokeColour(U, Tr, a,b);
+      }
+    }  	 
+
+  }
+};
+
+typedef	 AdjointRep<Nc> AdjointRepresentation;
+
+}
+}
+
+
+
+
+#endif

lib/qcd/representations/hmc_types.h

@@ -0,0 +1,46 @@
+#ifndef HMC_TYPES_H
+#define HMC_TYPES_H
+
+#include <tuple>
+#include <utility>
+#include <qcd/representations/adjoint.h>
+
+namespace Grid {
+namespace QCD {
+
+// Utility to add support for representations other than the fundamental
+
+template<class... Reptypes>
+class Representations{
+public:
+  typedef std::tuple<Reptypes...> Representation_type;
+  Representation_type rep;
+
+  // Multiple types constructor
+  explicit Representations(GridBase *grid):rep(Reptypes(grid)...){};
+
+  int size(){
+    return std::tuple_size< Representation_type >::value;
+  }
+
+  // update the fields
+  template <std::size_t I = 0>
+  inline typename std::enable_if< I == sizeof...(Reptypes), void >::type update(LatticeGaugeField& U) {}
+
+  template <std::size_t I = 0>
+      inline typename std::enable_if <
+      I<sizeof...(Reptypes), void >::type update(LatticeGaugeField& U) {
+    std::get<I>(rep).update_representation(U);
+    update<I + 1>(U);
+  }  
+};
+
+}
+}
+
+
+
+#endif
+
+
+

lib/qcd/utils/SUn.h

@@ -38,7 +38,8 @@ namespace QCD {
 template <int ncolour>
 class SU {
  public:
-  static int generators(void) { return ncolour * ncolour - 1; }
+  static const int Dimension = ncolour;
+  static const int AdjointDimension = ncolour * ncolour - 1;
   static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; }
 
   template <typename vtype>
@@ -46,11 +47,8 @@ class SU {
   template <typename vtype>
   using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >;
   template <typename vtype>
-  using iSUnAdjointMatrix = iScalar<iScalar<iMatrix<vtype, (ncolour*ncolour - 1)> > >;
+  using iSUnAlgebraVector =
-  template <typename vtype>
+      iScalar<iScalar<iVector<vtype, (ncolour * ncolour - 1)> > >;
-  using iSUnAlgebraVector = iScalar<iScalar<iVector<vtype , (ncolour*ncolour -1)> > >;
-
-
 
   //////////////////////////////////////////////////////////////////////////////////////////////////
   // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix,
@@ -64,16 +62,6 @@ class SU {
   typedef iSUnMatrix<vComplexF> vMatrixF;
   typedef iSUnMatrix<vComplexD> vMatrixD;
 
-  // Actually the adjoint matrices are real...
-  // Consider this overhead... FIXME
-  typedef iSUnAdjointMatrix<Complex> AMatrix;
-  typedef iSUnAdjointMatrix<ComplexF> AMatrixF;
-  typedef iSUnAdjointMatrix<ComplexD> AMatrixD;
-
-  typedef iSUnAdjointMatrix<vComplex> vAMatrix;
-  typedef iSUnAdjointMatrix<vComplexF> vAMatrixF;
-  typedef iSUnAdjointMatrix<vComplexD> vAMatrixD;
-
   // For the projectors to the algebra
   // these should be real...
   // keeping complex for consistency with the SIMD vector types
@@ -152,19 +140,6 @@ class SU {
   //   (    1   ) / sqrt(3) /2  = 1/2 lambda_8
   //   (      -2)
   //
-  //
-  // * Adjoint representation generators
-  //  
-  //   base for NxN hermitian traceless matrices
-  //   normalized to 1:
-  //
-  //   (e_Adj)^a = t^a / sqrt(T_F)
-  //   
-  //   then the real, antisymmetric generators for the adjoint representations
-  //   are computed ( shortcut: e^a == (e_Adj)^a )
-  //
-  //   (iT_adj)^d_ab = i tr[e^a t^d e^b - t^d e^a e^b]
-  // 
   ////////////////////////////////////////////////////////////////////////
   template <class cplx>
   static void generator(int lieIndex, iSUnMatrix<cplx> &ta) {
@@ -217,30 +192,6 @@ class SU {
     ta = ta * nrm;
   }
 
-  template <class cplx>
-  static void generatorAdjoint(int Index, iSUnAdjointMatrix<cplx> &iAdjTa){
-    // returns i(T_Adj)^index necessary for the projectors
-    // see definitions above
-    iAdjTa = zero;
-    Vector< iSUnMatrix<cplx> > ta(ncolour*ncolour -1);
-    iSUnMatrix<cplx> tmp;
-
-    // FIXME not very efficient to get all the generators everytime
-    for (int a = 0; a < (ncolour * ncolour - 1); a++) 
-      generator(a, ta[a]);
-
-
-    for (int a = 0; a < (ncolour*ncolour - 1); a++){
-      tmp = ta[a] * ta[Index] - ta[Index] * ta[a];
-      for (int b = 0; b < (ncolour*ncolour - 1); b++){
-        iSUnMatrix<cplx> tmp1 = 2.0 * tmp * ta[b];  // 2.0 from the normalization
-        Complex iTr = TensorRemove(timesI(trace(tmp1)));
-        iAdjTa()()(b,a) = iTr;
-      }
-    }
-
-
-  }
   ////////////////////////////////////////////////////////////////////////
   // Map a su2 subgroup number to the pair of rows that are non zero
   ////////////////////////////////////////////////////////////////////////
@@ -600,7 +551,7 @@ class SU {
   }
 
   static void printGenerators(void) {
-    for (int gen = 0; gen < generators(); gen++) {
+    for (int gen = 0; gen < AdjointDimension; gen++) {
       Matrix ta;
       generator(gen, ta);
       std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
@@ -609,71 +560,44 @@ class SU {
     }
   }
 
-  static void printAdjointGenerators(void) {
-    for (int gen = 0; gen < generators(); gen++) {
-      AMatrix ta;
-      generatorAdjoint(gen, ta);
-      std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
-                << std::endl;
-      std::cout << GridLogMessage << ta << std::endl;
-    }
-  }
-
 
 
   static void testGenerators(void) {
     Matrix ta;
     Matrix tb;
-    std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is 0.5 delta_ab"
+    std::cout << GridLogMessage
+              << "Fundamental - Checking trace ta tb is 0.5 delta_ab"
               << std::endl;
-    for (int a = 0; a < generators(); a++) {
+    for (int a = 0; a < AdjointDimension; a++) {
-      for (int b = 0; b < generators(); b++) {
+      for (int b = 0; b < AdjointDimension; b++) {
         generator(a, ta);
         generator(b, tb);
         Complex tr = TensorRemove(trace(ta * tb));
-        std::cout << GridLogMessage << "("<< a << "," << b << ") =  "<< tr << std::endl;
+        std::cout << GridLogMessage << "(" << a << "," << b << ") =  " << tr
+                  << std::endl;
         if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6);
         if (a != b) assert(abs(tr) < 1.0e-6);
       }
       std::cout << GridLogMessage << std::endl;
     }
-    std::cout << GridLogMessage << "Fundamental - Checking if hermitian" << std::endl;
+    std::cout << GridLogMessage << "Fundamental - Checking if hermitian"
-    for (int a = 0; a < generators(); a++) {
+              << std::endl;
+    for (int a = 0; a < AdjointDimension; a++) {
       generator(a, ta);
       std::cout << GridLogMessage << a << std::endl;
       assert(norm2(ta - adj(ta)) < 1.0e-6);
     }
     std::cout << GridLogMessage << std::endl;
 
-    std::cout << GridLogMessage << "Fundamental - Checking if traceless" << std::endl;
+    std::cout << GridLogMessage << "Fundamental - Checking if traceless"
-    for (int a = 0; a < generators(); a++) {
+              << std::endl;
+    for (int a = 0; a < AdjointDimension; a++) {
       generator(a, ta);
       Complex tr = TensorRemove(trace(ta));
       std::cout << GridLogMessage << a << " " << std::endl;
       assert(abs(tr) < 1.0e-6);
     }
     std::cout << GridLogMessage << std::endl;
-
-    AMatrix adjTa;
-    std::cout << GridLogMessage << "Adjoint - Checking if real" << std::endl;
-    for (int a = 0; a < generators(); a++) {
-      generatorAdjoint(a, adjTa);
-      std::cout << GridLogMessage << a << std::endl;
-      assert(norm2(adjTa - conjugate(adjTa)) < 1.0e-6);
-    }
-    std::cout << GridLogMessage << std::endl;
-
-    std::cout << GridLogMessage << "Adjoint - Checking if antisymmetric" << std::endl;
-    for (int a = 0; a < generators(); a++) {
-      generatorAdjoint(a, adjTa);
-      std::cout << GridLogMessage << a << std::endl;
-      assert(norm2(adjTa + transpose(adjTa)) < 1.0e-6);
-    }
-    std::cout << GridLogMessage << std::endl;
-
-
-
-
   }
 
   // reunitarise??
@@ -699,7 +623,7 @@ class SU {
     MatrixType ta;
 
     lie = zero;
-    for (int a = 0; a < generators(); a++) {
+    for (int a = 0; a < AdjointDimension; a++) {
       random(pRNG, ca);
 
       ca = (ca + conjugate(ca)) * 0.5;
@@ -724,7 +648,7 @@ class SU {
     Matrix ta;
 
     out = zero;
-    for (int a = 0; a < generators(); a++) {
+    for (int a = 0; a < AdjointDimension; a++) {
       gaussian(pRNG, ca);
       generator(a, ta);
       la = toComplex(ca) * ci * ta;
@@ -732,30 +656,23 @@ class SU {
     }
   }
 
-  static void FundamentalLieAlgebraMatrix(LatticeAlgebraVector &h, LatticeMatrix &out,
+  static void FundamentalLieAlgebraMatrix(LatticeAlgebraVector &h,
+                                          LatticeMatrix &out,
                                           Real scale = 1.0) {
+    conformable(h, out);
     GridBase *grid = out._grid;
     LatticeMatrix la(grid);
     Matrix ta;
 
     out = zero;
-    for (int a = 0; a < generators(); a++) {
+    for (int a = 0; a < AdjointDimension; a++) {
       generator(a, ta);
       la = peekColour(h, a) * scale * ta;
       out += la;
     }
   }
 
-  static void projectAdjointAlgebra(LatticeAlgebraVector &h_out, LatticeMatrix &in){
-    GridBase *grid = in._grid;
-    AMatrix iTa;
 
-    for (int a = 0; a< generators(); a++){
-      generatorAdjoint(a, iTa);
-      AlgebraVector tmp = real(trace(iTa * in));//*factor
-      pokeColour(h_out, tmp, a);
-    }
-  }
 
   template <typename GaugeField>
   static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) {

lib/qcd/utils/SUnAdjoint.h

@@ -0,0 +1,149 @@
+#ifndef QCD_UTIL_SUNADJOINT_H
+#define QCD_UTIL_SUNADJOINT_H
+
+////////////////////////////////////////////////////////////////////////
+//
+// * Adjoint representation generators
+//
+// * Normalisation for the fundamental generators: 
+//   trace ta tb = 1/2 delta_ab = T_F delta_ab
+//   T_F = 1/2  for SU(N) groups
+//
+//
+//   base for NxN hermitian traceless matrices
+//   normalized to 1:
+//
+//   (e_Adj)^a = t^a / sqrt(T_F)
+//
+//   then the real, antisymmetric generators for the adjoint representations
+//   are computed ( shortcut: e^a == (e_Adj)^a )
+//
+//   (iT_adj)^d_ba = i tr[e^a t^d e^b - t^d e^a e^b]
+//
+////////////////////////////////////////////////////////////////////////
+
+namespace Grid {
+namespace QCD {
+
+template <int ncolour>
+class SU_Adjoint : public SU<ncolour> {
+ public:
+  static const int Dimension = ncolour * ncolour - 1;
+
+  template <typename vtype>
+  using iSUnAdjointMatrix =
+      iScalar<iScalar<iMatrix<vtype, Dimension > > >;
+
+  // Actually the adjoint matrices are real...
+  // Consider this overhead... FIXME
+  typedef iSUnAdjointMatrix<Complex> AMatrix;
+  typedef iSUnAdjointMatrix<ComplexF> AMatrixF;
+  typedef iSUnAdjointMatrix<ComplexD> AMatrixD;
+
+  typedef iSUnAdjointMatrix<vComplex> vAMatrix;
+  typedef iSUnAdjointMatrix<vComplexF> vAMatrixF;
+  typedef iSUnAdjointMatrix<vComplexD> vAMatrixD;
+
+  typedef Lattice<vAMatrix>  LatticeAdjMatrix;
+  typedef Lattice<vAMatrixF> LatticeAdjMatrixF;
+  typedef Lattice<vAMatrixD> LatticeAdjMatrixD;
+
+
+  template <class cplx>
+  static void generator(int Index, iSUnAdjointMatrix<cplx> &iAdjTa) {
+    // returns i(T_Adj)^index necessary for the projectors
+    // see definitions above
+    iAdjTa = zero;
+    Vector<typename SU<ncolour>::template iSUnMatrix<cplx> > ta(ncolour * ncolour - 1);
+    typename SU<ncolour>::template iSUnMatrix<cplx> tmp;
+
+    // FIXME not very efficient to get all the generators everytime
+    for (int a = 0; a < Dimension; a++) SU<ncolour>::generator(a, ta[a]);
+
+    for (int a = 0; a < Dimension; a++) {
+      tmp = ta[a] * ta[Index] - ta[Index] * ta[a];
+      for (int b = 0; b < (ncolour * ncolour - 1); b++) {
+        typename SU<ncolour>::template iSUnMatrix<cplx> tmp1 =
+            2.0 * tmp * ta[b];  // 2.0 from the normalization
+        Complex iTr = TensorRemove(timesI(trace(tmp1)));
+        iAdjTa()()(b, a) = iTr;
+      }
+    }
+  }
+
+  static void printGenerators(void) {
+    for (int gen = 0; gen < Dimension; gen++) {
+      AMatrix ta;
+      generator(gen, ta);
+      std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
+                << std::endl;
+      std::cout << GridLogMessage << ta << std::endl;
+    }
+  }
+
+  static void testGenerators(void) {
+    AMatrix adjTa;
+    std::cout << GridLogMessage << "Adjoint - Checking if real" << std::endl;
+    for (int a = 0; a < Dimension; a++) {
+      generator(a, adjTa);
+      std::cout << GridLogMessage << a << std::endl;
+      assert(norm2(adjTa - conjugate(adjTa)) < 1.0e-6);
+    }
+    std::cout << GridLogMessage << std::endl;
+
+    std::cout << GridLogMessage << "Adjoint - Checking if antisymmetric"
+              << std::endl;
+    for (int a = 0; a < Dimension; a++) {
+      generator(a, adjTa);
+      std::cout << GridLogMessage << a << std::endl;
+      assert(norm2(adjTa + transpose(adjTa)) < 1.0e-6);
+    }
+    std::cout << GridLogMessage << std::endl;
+  }
+
+  // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 )
+  static void projectOnAlgebra(typename SU<ncolour>::LatticeAlgebraVector &h_out, LatticeAdjMatrix &in, Real scale = 1.0) {
+    conformable(h_out, in);
+    h_out = zero;
+    AMatrix iTa;
+
+    for (int a = 0; a < Dimension; a++) {
+      generator(a, iTa);
+      auto tmp = real(trace(iTa * in)) * scale; 
+      pokeColour(h_out, tmp, a);
+    }
+  }
+
+  // a projector that keeps the generators stored to avoid the overhead of recomputing. 
+  static void projector(typename SU<ncolour>::LatticeAlgebraVector &h_out, LatticeAdjMatrix &in, Real scale = 1.0) {
+    conformable(h_out, in);
+    static std::vector<AMatrix> iTa(Dimension);  // to store the generators
+    h_out = zero;
+    static bool precalculated = false; 
+    if (!precalculated){
+    	precalculated = true;
+        for (int a = 0; a < Dimension; a++) generator(a, iTa[a]);
+    }
+
+    for (int a = 0; a < Dimension; a++) {
+      auto tmp = real(trace(iTa[a] * in)) * scale; 
+      pokeColour(h_out, tmp, a);
+    }
+  }
+
+
+};
+
+
+
+
+
+
+typedef SU_Adjoint<2> SU2Adjoint;
+typedef SU_Adjoint<3> SU3Adjoint;
+typedef SU_Adjoint<4> SU4Adjoint;
+typedef SU_Adjoint<5> SU5Adjoint;
+}
+}
+
+#endif

tests/Make.inc

@@ -1,5 +1,5 @@
 
-bin_PROGRAMS += Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_cshift_red_black_rotate Test_cshift_rotate Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_gpforce Test_dwf_hdcr Test_dwf_lanczos Test_dwf_mixedcg_prec Test_dwf_rb5d Test_gamma Test_GaugeAction Test_gparity Test_gpdwf_force Test_gp_rect_force Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_rect_force Test_RectPlaq Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_wilson_tm_even_odd 
+bin_PROGRAMS += Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_cshift_red_black_rotate Test_cshift_rotate Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_gpforce Test_dwf_hdcr Test_dwf_lanczos Test_dwf_mixedcg_prec Test_dwf_rb5d Test_gamma Test_GaugeAction Test_gparity Test_gpdwf_force Test_gp_rect_force Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_rect_force Test_RectPlaq Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_wilson_tm_even_odd 
 
 
 Test_cayley_cg_SOURCES=Test_cayley_cg.cc
@@ -162,6 +162,10 @@ Test_hmc_RectGauge_SOURCES=Test_hmc_RectGauge.cc
 Test_hmc_RectGauge_LDADD=-lGrid
 
 
+Test_hmc_WilsonAdjointFermionGauge_SOURCES=Test_hmc_WilsonAdjointFermionGauge.cc
+Test_hmc_WilsonAdjointFermionGauge_LDADD=-lGrid
+
+
 Test_hmc_WilsonFermionGauge_SOURCES=Test_hmc_WilsonFermionGauge.cc
 Test_hmc_WilsonFermionGauge_LDADD=-lGrid
 

tests/Test_hmc_WilsonAdjointFermionGauge.cc

@@ -0,0 +1,100 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonAdjointFermionGauge.cc
+
+Copyright (C) 2015
+
+Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: neo <cossu@post.kek.jp>
+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 distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include "Grid.h"
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+namespace Grid {
+namespace QCD {
+
+class HmcRunner : public NerscHmcRunner {
+ public:
+  void BuildTheAction(int argc, char **argv)
+
+  {
+    typedef WilsonImplR ImplPolicy; // gauge field implemetation
+    typedef WilsonFermionR FermionAction; // type of lattice fermions
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid = SpaceTimeGrid::makeFourDimGrid(
+        GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
+        GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+    FGrid = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    Real mass = -0.77;
+    FermionAction FermOp(U, *FGrid, *FrbGrid, mass);
+
+    ConjugateGradient<FermionField> CG(1.0e-8, 10000);
+
+    TwoFlavourPseudoFermionAction<ImplPolicy> Nf2(FermOp, CG, CG);
+
+    // Set smearing (true/false), default: false
+    Nf2.is_smeared = true;
+
+    // Collect actions
+    ActionLevel<LatticeGaugeField> Level1(1);
+    Level1.push_back(&Nf2);
+
+    ActionLevel<LatticeGaugeField> Level2(4);
+    Level2.push_back(&Waction);
+
+    TheAction.push_back(Level1);
+    TheAction.push_back(Level2);
+
+    Run(argc, argv);
+  };
+};
+}
+}
+
+int main(int argc, char **argv) {
+  Grid_init(&argc, &argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout << GridLogMessage << "Grid is setup to use " << threads
+            << " threads" << std::endl;
+
+  HmcRunner TheHMC;
+
+  TheHMC.BuildTheAction(argc, argv);
+}

tests/Test_hmc_WilsonFermionGauge.cc

@@ -25,7 +25,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 with this program; if not, write to the Free Software Foundation, Inc.,
 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
+See the full license in the file "LICENSE" in the top level distribution
+directory
 *************************************************************************************/
 /*  END LEGAL */
 #include "Grid.h"
@@ -37,10 +38,8 @@ using namespace Grid::QCD;
 namespace Grid {
 namespace QCD {
 
-
 class HmcRunner : public NerscHmcRunner {
  public:
-
   void BuildTheAction(int argc, char **argv)
 
   {
@@ -48,7 +47,9 @@ public:
     typedef WilsonFermionR FermionAction;
     typedef typename FermionAction::FermionField FermionField;
 
-    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UGrid = SpaceTimeGrid::makeFourDimGrid(
+        GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
+        GridDefaultMpi());
     UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
 
     FGrid = UGrid;
@@ -70,7 +71,6 @@ public:
     // Set smearing (true/false), default: false
     Nf2.is_smeared = true;
 
-
     // Collect actions
     ActionLevel<LatticeGaugeField> Level1(1);
     Level1.push_back(&Nf2);
@@ -83,22 +83,18 @@ public:
 
     Run(argc, argv);
   };
-
 };
+}
+}
 
-}}
+int main(int argc, char **argv) {
-
-int main (int argc, char ** argv)
-{
   Grid_init(&argc, &argv);
 
   int threads = GridThread::GetThreads();
-  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+  std::cout << GridLogMessage << "Grid is setup to use " << threads
+            << " threads" << std::endl;
 
   HmcRunner TheHMC;
 
   TheHMC.BuildTheAction(argc, argv);
-
 }
-
-

tests/Test_lie_generators.cc

@@ -31,6 +31,8 @@ directory
 
 #include <qcd/utils/CovariantCshift.h>
 #include <qcd/utils/SUn.h>
+#include <qcd/utils/SUnAdjoint.h>
+#include <qcd/representations/adjoint.h>
 #include <qcd/utils/WilsonLoops.h>
 
 using namespace std;
@@ -52,8 +54,10 @@ int main(int argc, char** argv) {
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
   SU2::printGenerators();
-  SU2::printAdjointGenerators();
+  std::cout << "Dimension of adjoint representation: "<< SU2Adjoint::Dimension << std::endl;
+  SU2Adjoint::printGenerators();
   SU2::testGenerators();
+  SU2Adjoint::testGenerators();
 
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
@@ -61,14 +65,19 @@ int main(int argc, char** argv) {
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
   SU3::printGenerators();
-  SU3::printAdjointGenerators();
+  std::cout << "Dimension of adjoint representation: "<< SU3Adjoint::Dimension << std::endl;
+  SU3Adjoint::printGenerators();
   SU3::testGenerators();
+  SU3Adjoint::testGenerators();
 
   std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
   std::cout<<GridLogMessage<<"* Generators for SU(4)"<<std::endl;
   std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
   SU4::printGenerators();
+  std::cout << "Dimension of adjoint representation: "<< SU4Adjoint::Dimension << std::endl;
+  SU4Adjoint::printGenerators();
   SU4::testGenerators();
+  SU4Adjoint::testGenerators();
 
   //  std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
   //  std::cout<<GridLogMessage<<"* Generators for SU(5)"<<std::endl;
@@ -76,5 +85,34 @@ int main(int argc, char** argv) {
   //  SU5::printGenerators();
   //  SU5::testGenerators();
 
+
+  // Projectors 
+  GridParallelRNG gridRNG(grid);
+  gridRNG.SeedRandomDevice();
+  SU3Adjoint::LatticeAdjMatrix Gauss(grid);
+  SU3::LatticeAlgebraVector ha(grid);
+  SU3::LatticeAlgebraVector hb(grid);
+  random(gridRNG,Gauss);
+
+  std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl;
+  SU3Adjoint::projectOnAlgebra(ha, Gauss);
+  std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl;
+  std::cout << GridLogMessage << "Start projector" << std::endl;
+  SU3Adjoint::projector(hb, Gauss);
+  std::cout << GridLogMessage << "end projector" << std::endl;
+
+  std::cout << GridLogMessage << "ReStart projector" << std::endl;
+  SU3Adjoint::projector(hb, Gauss);
+  std::cout << GridLogMessage << "end projector" << std::endl;
+  SU3::LatticeAlgebraVector diff = ha -hb;
+  std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl;
+
+
+  // Testing HMC representation classes
+  AdjointRep<3> AdjRep(grid);
+
+  // AdjointRepresentation has the predefined number of colours Nc
+  Representations<AdjointRepresentation> RepresentationTypes(grid);  
+
   Grid_finalize();
 }

tests/Test_wilson_force.cc

@@ -76,6 +76,8 @@ int main (int argc, char ** argv)
 
   Dw.MDeriv(tmp , Mphi,  phi,DaggerNo );  UdSdU=tmp;
   Dw.MDeriv(tmp , phi,  Mphi,DaggerYes ); UdSdU=(UdSdU+tmp);
+  // Take the trace
+  UdSdU = Ta(UdSdU);
   
   LatticeFermion Ftmp      (&Grid);
 

tests/Test_zmm.cc

@@ -22,16 +22,14 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 with this program; if not, write to the Free Software Foundation, Inc.,
 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
+See the full license in the file "LICENSE" in the top level distribution
+directory
 *************************************************************************************/
 /*  END LEGAL */
 #include <Grid.h>
 #include <PerfCount.h>
 
-int main(int argc,char **argv)
+int main(int argc, char **argv) { return 0; }
-{
-  return 0;
-}
 #if 0
 #include <simd/Intel512wilson.h>
 using namespace Grid;