diff --git a/lib/lattice/Lattice_ET.h b/lib/lattice/Lattice_ET.h
index 742c4537..77e0ce5d 100644
--- a/lib/lattice/Lattice_ET.h
+++ b/lib/lattice/Lattice_ET.h
@@ -213,6 +213,7 @@ inline void CBFromExpression( int &cb,const LatticeTrinaryExpression<Op,T1,T2,T3
  GridUnopClass(UnaryToReal,toReal(a));
  GridUnopClass(UnaryToComplex,toComplex(a));
  GridUnopClass(UnaryTimesI,timesI(a));
+ GridUnopClass(UnaryTimesMinusI,timesMinusI(a));
  GridUnopClass(UnaryAbs,abs(a));
  GridUnopClass(UnarySqrt,sqrt(a));
  GridUnopClass(UnaryRsqrt,rsqrt(a));
@@ -329,6 +330,7 @@ GRID_DEF_UNOP(imag,UnaryImag);
 GRID_DEF_UNOP(toReal,UnaryToReal);
 GRID_DEF_UNOP(toComplex,UnaryToComplex);
 GRID_DEF_UNOP(timesI,UnaryTimesI);
+GRID_DEF_UNOP(timesMinusI,UnaryTimesMinusI);
 GRID_DEF_UNOP(abs  ,UnaryAbs); //abs overloaded in cmath C++98; DON'T do the abs-fabs-dabs-labs thing
 GRID_DEF_UNOP(sqrt ,UnarySqrt);
 GRID_DEF_UNOP(rsqrt,UnaryRsqrt);
diff --git a/lib/qcd/action/gauge/GaugeImpl.h b/lib/qcd/action/gauge/GaugeImpl.h
index 2181ec6a..808feca2 100644
--- a/lib/qcd/action/gauge/GaugeImpl.h
+++ b/lib/qcd/action/gauge/GaugeImpl.h
@@ -61,7 +61,7 @@ template<class Gimpl> class WilsonLoops;
       typedef Lattice<SiteGaugeField>                   GaugeField;
 
     // Move this elsewhere?
-    void AddGaugeLink(GaugeField& U, GaugeLinkField& W, int mu){  // U[mu] += W 
+    static inline void AddGaugeLink(GaugeField& U, GaugeLinkField& W, int mu){  // U[mu] += W 
 PARALLEL_FOR_LOOP
       for(auto ss=0;ss<U._grid->oSites();ss++){
 	U._odata[ss]._internal[mu] = U._odata[ss]._internal[mu] + W._odata[ss]._internal;
diff --git a/lib/qcd/smearing/APEsmearing.h b/lib/qcd/smearing/APEsmearing.h
index 24ac4b96..2580d296 100644
--- a/lib/qcd/smearing/APEsmearing.h
+++ b/lib/qcd/smearing/APEsmearing.h
@@ -86,19 +86,19 @@ namespace Grid {
 
 	    temp_Sigma = -rho_numu*staple*iLambda_nu;
 	    //-r_numu*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)*Lambda_nu(x)
-	    AddGaugeLink(SigmaTerm, temp_Sigma, mu);
+	    Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
 
 	    sh_field = Cshift(iLambda_nu, mu, 1);// general also for Gparity?
 	
 	    temp_Sigma = rho_numu*sh_field*staple;
 	    //r_numu*Lambda_nu(mu)*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)
-	    AddGaugeLink(SigmaTerm, temp_Sigma, mu);
+	    Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
 
 	    sh_field = Cshift(iLambda_mu, nu, 1);
 
 	    temp_Sigma = -rho_munu*staple*U_nu*sh_field*adj(U_nu);
 	    //-r_munu*U_nu(x+mu)*Udag_mu(x+nu)*Lambda_mu(x+nu)*Udag_nu(x)
-	    AddGaugeLink(SigmaTerm, temp_Sigma, mu);
+	    Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
 
 	    staple = zero;
 	    sh_field = Cshift(U_nu, mu, 1);
@@ -110,7 +110,7 @@ namespace Grid {
 	    sh_field = Cshift(u_tmp, mu, 1);
 	    temp_Sigma += -rho_numu*sh_field*adj(U_mu)*U_nu;
 	    sh_field = Cshift(temp_Sigma, nu, -1);
-	    AddGaugeLink(SigmaTerm, sh_field, mu);
+	    Gimpl::AddGaugeLink(SigmaTerm, sh_field, mu);
 	    
 	  }
 	}
diff --git a/lib/qcd/smearing/GaugeConfiguration.h b/lib/qcd/smearing/GaugeConfiguration.h
index f38af279..42baced8 100644
--- a/lib/qcd/smearing/GaugeConfiguration.h
+++ b/lib/qcd/smearing/GaugeConfiguration.h
@@ -24,34 +24,149 @@ namespace Grid {
     template <class Gimpl>
     class SmearedConfiguration {
     public:
-      INHERIT_GIMPL_TYPES(Gimpl)
-      private:
+      INHERIT_GIMPL_TYPES(Gimpl) ;
+      
+    private:
       const unsigned int smearingLevels;
-      Smear_Stout StoutSmearing;
+      Smear_Stout<Gimpl> StoutSmearing;
       std::vector<GaugeField> SmearedSet;
       
       // Member functions
-      void fill_smearedSet();
-      GaugeField AnalyticSmearedForce(const GaugeField&, 
-				      const GaugeField&) const;
-      const GaugeField& get_smeared_conf(int) const;
+      //====================================================================
+      void fill_smearedSet(){
+	GaugeField previous_u;
+	
+	std::cout<< GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n";
+	
+	previous_u = *ThinLinks;
+	for(int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl){
+	  StoutSmearing.smear(SmearedSet[smearLvl],previous_u);
+	  previous_u = SmearedSet[smearLvl];
+	}
+	
+      }
+      //====================================================================
+      GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime, 
+				      const GaugeField& GaugeK) const{
+	GridBase *grid = GaugeK._grid;
+	GaugeField C(grid), SigmaK(grid), iLambda(grid);
+	GaugeLinkField iLambda_mu(grid);
+	GaugeLinkField iQ(grid), e_iQ(grid);
+	GaugeLinkField SigmaKPrime_mu(grid);
+	GaugeLinkField GaugeKmu(grid), Cmu(grid);
+
+	StoutSmearing.BaseSmear(C, GaugeK);
+
+	for (int mu = 0; mu < Nd; mu++){
+	  Cmu            = peekLorentz(     C,mu);
+	  GaugeKmu       = peekLorentz(GaugeK,mu);
+	  SigmaKPrime_mu = peekLorentz(SigmaKPrime,mu);
+	  iQ = Ta(Cmu*adj(GaugeKmu));
+	  set_iLambda(iLambda_mu, e_iQ, iQ, SigmaKPrime_mu, GaugeKmu);
+	  pokeLorentz(SigmaK, SigmaKPrime_mu*e_iQ + adj(Cmu)*iLambda_mu, mu);
+	  pokeLorentz(iLambda, iLambda_mu, mu);
+	}
+	StoutSmearing.derivative(SigmaK, iLambda, GaugeK);
+	return SigmaK;
+      }
+      /*! @brief Returns smeared configuration at level 'Level' */
+      const GaugeField& get_smeared_conf(int Level) const{
+	return SmearedSet[Level];
+      }
       
-      void set_iLambda(GaugeField& iLambda, 
-		       GaugeField& e_iQ,
-		       const GaugeField& iQ, 
-		       const GaugeField& Sigmap,
-		       const GaugeField& U)const;
-      
-      /* Check these types (do I need to pass iQ1,2 ? )
-      void set_uw(RealD& u, RealD& w,
-		  const SUNmat& iQ1, const SUNmat& iQ2)const ;
-      void set_fj(ComplexD& f0, ComplexD& f1,
-		  CompledD& f2, const RealD& u,
-		  const RealD& w)const;
-      */
-      
-      RealD func_xi0(RealD w)const;
-      RealD func_xi1(RealD w)const; 
+      void set_iLambda(GaugeLinkField& iLambda, 
+		       GaugeLinkField& e_iQ,
+		       const GaugeLinkField& iQ, 
+		       const GaugeLinkField& Sigmap,
+		       const GaugeLinkField& GaugeK)const{
+	GridBase *grid = iQ._grid;
+	GaugeLinkField iQ2(grid), iQ3(grid), B1(grid), B2(grid), USigmap(grid);
+	GaugeLinkField unity(grid);
+	unity=1.0;
+	
+	LatticeReal u(grid), w(grid);
+	LatticeComplex f0(grid), f1(grid), f2(grid);
+	LatticeReal xi0(grid), xi1(grid), fden(grid), tmp(grid);
+	LatticeReal u2(grid), w2(grid), cosw(grid);
+	LatticeComplex emiu(grid), e2iu(grid), qt(grid);
+	LatticeComplex r01(grid), r11(grid), r21(grid), r02(grid), r12(grid);
+	LatticeComplex r22(grid), tr1(grid), tr2(grid);
+	LatticeComplex b10(grid), b11(grid), b12(grid), b20(grid), b21(grid), b22(grid);
+	LatticeReal unitReal(grid);
+
+	unitReal = 1.0;
+	
+	// Exponential
+	iQ2 = iQ * iQ;
+	iQ3 = iQ * iQ2;
+	StoutSmearing.set_uw(u,w,iQ2,iQ3);
+	StoutSmearing.set_fj(f0,f1,f2,u,w);
+	e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2;
+
+	// Getting B1, B2, Gamma and Lambda
+	xi0 = StoutSmearing.func_xi0(w);
+	xi1 = StoutSmearing.func_xi1(w);
+	u2 = u * u;
+	w2 = w * w;
+	cosw = cos(w);
+	
+	emiu = toComplex(cos(u)) - timesI(toComplex(u));
+	e2iu = toComplex(cos(2.0*u)) + timesI(toComplex(2.0*u));
+
+	r01 = (toComplex(2.0*u) + timesI(toComplex(2.0*(u2-w2)))) * e2iu
+	  + emiu * (toComplex(16.0*u*cosw + 2.0*u*(3.0*u2+w2)*xi0) +
+		    timesI(toComplex(-8.0*u2*cosw + 2.0*(9.0*u2+w2)*xi0)));
+	
+	r11 = (toComplex(2.0*unitReal) + timesI(toComplex(4.0*u)))* e2iu
+	  + emiu * (toComplex(-2.0*cosw + (3.0*u2-w2)*xi0) +
+		    timesI(toComplex(2.0*u*cosw + 6.0*u*xi0)));
+
+	r21 = timesI(toComplex(2.0*unitReal)) * e2iu
+	  + emiu * (toComplex(-3.0*u*xi0) + timesI(toComplex(cosw - 3.0*xi0)));
+
+	
+	r02 = -2.0 * e2iu + emiu * (toComplex(-8.0*u2*xi0) +
+				    timesI(toComplex(2.0*u*(cosw + xi0 + 3.0*u2*xi1))));
+
+	r12 = emiu * (toComplex(2.0*u*xi0) + timesI(toComplex(-cosw - xi0 + 3.0*u2*xi1)));
+
+	r22 = emiu * (toComplex(xi0) - timesI(toComplex(3.0*u*xi1)));
+
+	tmp = (2.0*(9.0*u2-w2)*(9.0*u2-w2));
+	fden = toComplex(pow(tmp, -1.0));  // 1/tmp
+
+	b10 = toComplex(2.0*u) * r01 + toComplex(3.0*u2 - w2)*r02 - toComplex(30.0*u2 + 2.0*w2)*f0;
+	b11 = toComplex(2.0*u) * r11 + toComplex(3.0*u2 - w2)*r12 - toComplex(30.0*u2 + 2.0*w2)*f1;
+	b12 = toComplex(2.0*u) * r21 + toComplex(3.0*u2 - w2)*r22 - toComplex(30.0*u2 + 2.0*w2)*f2;
+
+	b20 = r01 - toComplex(3.0*u)*r02 - toComplex(24.0*u)*f0;
+	b21 = r11 - toComplex(3.0*u)*r12 - toComplex(24.0*u)*f1;
+	b22 = r21 - toComplex(3.0*u)*r22 - toComplex(24.0*u)*f2;
+
+	b10 *= fden;
+	b11 *= fden;
+	b12 *= fden;
+	b20 *= fden;
+	b21 *= fden;
+	b22 *= fden;
+	
+	B1 = b10*unity + timesMinusI(b11) * iQ - b12 * iQ2;
+	B2 = b20*unity + timesMinusI(b21) * iQ - b22 * iQ2;
+	USigmap = GaugeK * Sigmap;
+
+	tr1 = trace(USigmap*B1);
+	tr2 = trace(USigmap*B2);
+
+	GaugeLinkField QUS = timesMinusI(iQ) * USigmap;
+	GaugeLinkField USQ = USigmap * timesMinusI(iQ);
+
+	GaugeLinkField iGamma = tr1 * timesMinusI(iQ) - tr2 * iQ2 +
+	  f1 * USigmap + f2 * QUS + f2 * USQ;
+
+	iLambda = Ta(iGamma);
+	
+		
+      }
       
     public:
       GaugeField* ThinLinks;      /*!< @brief Pointer to the thin 
@@ -59,11 +174,11 @@ namespace Grid {
       
       /*! @brief Standard constructor */
       SmearedConfiguration(GridCartesian * UGrid,
-			 unsigned int Nsmear, 
-			 Smear_Stout& Stout):
+			   unsigned int Nsmear, 
+			   Smear_Stout<Gimpl>& Stout):
 	smearingLevels(Nsmear),
 	StoutSmearing(Stout),
-	ThinLinks(new GaugeField){
+	ThinLinks(new GaugeField(UGrid)){
 	for (unsigned int i=0; i< smearingLevels; ++i)
 	  SmearedSet.push_back(*(new GaugeField(UGrid)));
       }
@@ -76,7 +191,26 @@ namespace Grid {
 	ThinLinks(new GaugeField(UGrid)){}
       
       void set_GaugeField(){ fill_smearedSet(); }
-      void smeared_force(GaugeField&) const;
+      void smeared_force(GaugeField& SigmaTilde) const{
+	GaugeField force = SigmaTilde;//actually = U*SigmaTilde, check this for Grid
+	GaugeLinkField tmp_mu(SigmaTilde._grid);
+	
+	for (int mu = 0; mu < Nd; mu++){
+	  tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels-1], mu)) * peekLorentz(force,mu);
+	  pokeLorentz(force, tmp_mu, mu);
+	}
+	for(int ismr = smearingLevels - 1; ismr > 0; --ismr)
+	  force = AnalyticSmearedForce(force,get_smeared_conf(ismr-1));
+
+	force = AnalyticSmearedForce(force,*ThinLinks);
+
+	for (int mu = 0; mu < Nd; mu++){
+	  tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu);
+	  pokeLorentz(SigmaTilde, tmp_mu, mu);
+	}
+      }
+
+
       GaugeField* get_SmearedU() const{ 
 	return const_cast<GaugeField*>(&(SmearedSet[smearingLevels-1]));
       }
diff --git a/lib/qcd/smearing/Smearing.h b/lib/qcd/smearing/Smearing.h
index 76d7d935..89842c80 100644
--- a/lib/qcd/smearing/Smearing.h
+++ b/lib/qcd/smearing/Smearing.h
@@ -4,6 +4,6 @@
 #include <qcd/smearing/BaseSmearing.h>
 #include <qcd/smearing/APEsmearing.h>
 #include <qcd/smearing/StoutSmearing.h>
-
+#include <qcd/smearing/GaugeConfiguration.h>
 
 #endif
diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h
index a9faed23..3538417a 100644
--- a/lib/qcd/smearing/StoutSmearing.h
+++ b/lib/qcd/smearing/StoutSmearing.h
@@ -15,11 +15,7 @@ namespace Grid {
       const std::vector<double> d_rho;
       const Smear < Gimpl > * SmearBase;
       
-      LatticeReal func_xi0(LatticeReal w) const{
-	// Define a function to do the check
-	//if( w < 1e-4 ) std::cout << GridLogWarning << "[Smear_stout] w too small: "<< w <<"\n";
-	return  sin(w)/w;
-      }
+
       
     public:
       INHERIT_GIMPL_TYPES(Gimpl)
@@ -32,20 +28,20 @@ namespace Grid {
       ~Smear_Stout(){}
       
       void smear(GaugeField& u_smr,const GaugeField& U) const{
-	long double timing;
-	
-	GaugeField u_tmp1, q_mu;
+	GaugeField C(U._grid);
+	GaugeLinkField tmp(U._grid), q_mu(U._grid), Umu(U._grid);
 	
 	std::cout<< GridLogDebug << "Stout smearing started\n";
 	
 	//Smear the configurations
-	SmearBase->smear(u_tmp1, U);
-	
-	q_mu = Ta(u_tmp1*adj(u_tmp1)); // q_mu = Ta(Omega_mu)
-	
-	exponentiate_iQ(u_tmp1, q_mu);
-
-	u_smr = u_tmp1*U;
+	SmearBase->smear(C, U);
+	for (int mu = 0; mu<Nd; mu++){
+	  tmp = peekLorentz(C,mu);
+	  Umu = peekLorentz(U,mu);
+	  q_mu = Ta(tmp * adj(Umu)); // q_mu = Ta(Omega_mu)
+	  exponentiate_iQ(tmp, q_mu);  
+	  pokeLorentz(u_smr, tmp*Umu, mu);// u_smr = exp(iQ_mu)*U_mu
+	}
 	
 	std::cout<< GridLogDebug << "Stout smearing completed\n";
       }
@@ -61,72 +57,95 @@ namespace Grid {
 	SmearBase->smear(C, U);
       }
       
-      void exponentiate_iQ(GaugeField& e_iQ,
-			   const GaugeField& iQ) const{
+      void exponentiate_iQ(GaugeLinkField& e_iQ,
+			   const GaugeLinkField& iQ) const{
 	// Put this outside 
 	// only valid for SU(3) matrices
 
+	// only one Lorentz direction at a time 
+
 	GridBase *grid = iQ._grid;
-	Real one_over_three = 1.0/3.0;
-	Real one_over_two = 1.0/2.0;
+	GaugeLinkField unity(grid);
+	unity=1.0;
 
-	GaugeField unity;
-	GaugeLinkField Umu(iQ._grid);
-	Umu=1.0;
-	for(int mu=0;mu<Nd;mu++){
-	  pokeLorentz(unity,Umu,mu);
-	}
-
-	
-	GaugeField iQ2, iQ3;
-	LatticeReal c0(grid), c1(grid), c0max(grid), u_val(grid), tmp(grid);
-	LatticeReal w(grid), theta(grid), xi0(grid), u2(grid), w2(grid), cosw(grid);
-	LatticeComplex fden(grid);
-	LatticeComplex f0(grid), f1(grid), f2(grid), h0(grid), h1(grid), h2(grid);
-	LatticeComplex e2iu(grid), emiu(grid), ixi0(grid), qt(grid);
+	GaugeLinkField iQ2(grid), iQ3(grid);
+	LatticeReal u(grid), w(grid);
+	LatticeComplex f0(grid), f1(grid), f2(grid);
 
 	iQ2 = iQ * iQ;
 	iQ3 = iQ * iQ2;
 
-	c0    = - imag(trace(iQ3)) * one_over_three;
-	c1    = - real(trace(iQ2)) * one_over_two;
+	set_uw(u, w, iQ2, iQ3);
+	set_fj(f0, f1, f2, u, w);
+	e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2;
+      };
+
+
+      void set_uw(LatticeReal& u, LatticeReal& w,
+		  GaugeLinkField& iQ2, GaugeLinkField& iQ3) const{
+	Real one_over_three = 1.0/3.0;
+	Real one_over_two = 1.0/2.0;
+
+	GridBase *grid = u._grid;
+	LatticeReal c0(grid), c1(grid), tmp(grid), c0max(grid), theta(grid);
+	
+	
+	c0    = - toReal(imag(trace(iQ3))) * one_over_three;
+	c1    = - toReal(real(trace(iQ2))) * one_over_two;
 	tmp   = c1 * one_over_three;
     	c0max = 2.0 * pow(tmp, 1.5);
-
+	
 	theta = acos(c0/c0max);
+	u = sqrt(tmp) * cos( theta * one_over_three);
+	w = sqrt(c1) * sin ( theta * one_over_three);
+      }
 
-	u_val = sqrt(tmp) * cos( theta * one_over_three);
-	w     = sqrt(c1) * sin ( theta * one_over_three);
+      void set_fj(LatticeComplex& f0, LatticeComplex& f1, LatticeComplex& f2,
+		  const LatticeReal& u, const LatticeReal& w) const{
+
+	GridBase *grid = u._grid;
+	LatticeReal xi0(grid), u2(grid), w2(grid), cosw(grid), tmp(grid);
+	LatticeComplex fden(grid);
+	LatticeComplex h0(grid), h1(grid), h2(grid);
+	LatticeComplex e2iu(grid), emiu(grid), ixi0(grid), qt(grid);
+	
 	xi0   = func_xi0(w);
-	u2    = u_val * u_val;
+	u2    = u * u;
 	w2    = w * w;
 	cosw  = cos(w);
-
+	
 	ixi0  = timesI(toComplex(xi0));
-	emiu  = toComplex(cos(u_val)) - timesI(toComplex(u_val));
-	e2iu  = toComplex(cos(2.0*u_val)) + timesI(toComplex(2.0*u_val));
+	emiu  = toComplex(cos(u)) - timesI(toComplex(u));
+	e2iu  = toComplex(cos(2.0*u)) + timesI(toComplex(2.0*u));
 	
 	h0    = e2iu * toComplex(u2 - w2) + emiu *( toComplex(8.0*u2*cosw) +
-						    toComplex(2.0*u_val*(3.0*u2 + w2))*ixi0);
+						    toComplex(2.0*u*(3.0*u2 + w2))*ixi0);
 	
-	h1    = toComplex(2.0*u_val) * e2iu - emiu*( toComplex(2.0*u_val*cosw) -
-						     toComplex(3.0*u2-w2)*ixi0);
+	h1    = toComplex(2.0*u) * e2iu - emiu*( toComplex(2.0*u*cosw) -
+						 toComplex(3.0*u2-w2)*ixi0);
+	
+	h2    = e2iu - emiu * (toComplex(cosw) + toComplex(3.0*u)*ixi0);
 	
-	h2    = e2iu - emiu * (toComplex(cosw) + toComplex(3.0*u_val)*ixi0);
-
 	
 	tmp   = 9.0*u2 - w2;
 	fden  = toComplex(pow(tmp, -1.0));
 	f0    = h0 * fden;
 	f1    = h1 * fden;
-	f2    = h2 * fden;
+	f2    = h2 * fden;	
+      }
 
-	
-	e_iQ = f0*unity + f1 * timesMinusI(iQ) - f2 * iQ2;
-	
-	
-	
-      };
+
+      LatticeReal func_xi0(LatticeReal w) const{
+	// Define a function to do the check
+	//if( w < 1e-4 ) std::cout << GridLogWarning<< "[Smear_stout] w too small: "<< w <<"\n";
+	return  sin(w)/w;
+      }
+
+      LatticeReal func_xi1(LatticeReal w) const{
+	// Define a function to do the check
+	//if( w < 1e-4 ) std::cout << GridLogWarning << "[Smear_stout] w too small: "<< w <<"\n";
+	return  cos(w)/(w*w) - sin(w)/(w*w*w);
+      }
       
     };
 
diff --git a/lib/simd/Grid_vector_unops.h b/lib/simd/Grid_vector_unops.h
index bb3c44b0..2ce8590d 100644
--- a/lib/simd/Grid_vector_unops.h
+++ b/lib/simd/Grid_vector_unops.h
@@ -118,21 +118,21 @@ namespace Grid {
   };
 
   template<class scalar> struct RealFunctor {
-    scalar operator()(const scalar &a)  const {
+    scalar operator()(const std::complex<scalar> &a)  const {
       return real(a);
     }
   };
   template<class scalar> struct ImagFunctor {
-    scalar operator()(const scalar &a)  const {
+    scalar operator()(const std::complex<scalar> &a)  const {
       return imag(a);
     }
   };
   template < class S, class V > 
-  inline Grid_simd<S,V> real(const Grid_simd<S,V> &r) {
+    inline Grid_simd<S,V> real(const Grid_simd<S,V> &r) {
     return SimdApply(RealFunctor<S>(),r);
   }
   template < class S, class V > 
-  inline Grid_simd<S,V> imag(const Grid_simd<S,V> &r) {
+    inline Grid_simd<S,V> imag(const Grid_simd<S,V> &r) {
     return SimdApply(ImagFunctor<S>(),r);
   }
 
diff --git a/tests/Test_main.cc b/tests/Test_main.cc
index 7904ff21..467c8a5e 100644
--- a/tests/Test_main.cc
+++ b/tests/Test_main.cc
@@ -609,8 +609,10 @@ int main (int argc, char ** argv)
 
 
  // Testing Smearing routine compilation
+ GridCartesian           Fine(latt_size,simd_layout,mpi_layout);
  Smear_APE< PeriodicGimplR > APEsmearing; // periodic gauge implemetation
- 
+ Smear_Stout< PeriodicGimplR > StoutSmearing(&APEsmearing);
+ SmearedConfiguration< PeriodicGimplR > SmartConf(&Fine, 3, StoutSmearing);
  
  std::cout<<GridLogMessage << sizeof(vComplexF) << std::endl;