From 0fa66e8f3ca4d7321224d2246d29acf3919b0c4e Mon Sep 17 00:00:00 2001
From: Guido Cossu <guido.cossu@ed.ac.uk>
Date: Mon, 4 Jul 2016 15:35:37 +0100
Subject: [PATCH] Debugged smearing for EOWilson, accepts

---
 lib/qcd/hmc/integrators/Integrator.h  | 238 +++++++++++++-------------
 lib/qcd/smearing/APEsmearing.h        | 121 +++++++------
 lib/qcd/smearing/GaugeConfiguration.h | 163 +++++++++---------
 lib/qcd/smearing/StoutSmearing.h      | 183 +++++---------------
 4 files changed, 305 insertions(+), 400 deletions(-)

diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h
index 58e1f585..db61b114 100644
--- a/lib/qcd/hmc/integrators/Integrator.h
+++ b/lib/qcd/hmc/integrators/Integrator.h
@@ -44,40 +44,40 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 #include <memory>
 
-namespace Grid{
-  namespace QCD{
+ namespace Grid{
+ 	namespace QCD{
 
-    struct IntegratorParameters{
+ 		struct IntegratorParameters{
 
-      int Nexp;
+ 			int Nexp;
       int MDsteps;  // number of outer steps
       RealD trajL;  // trajectory length 
       RealD stepsize;
 
       IntegratorParameters(int MDsteps_, 
-			   RealD trajL_=1.0,
-			   int Nexp_=12):
-        Nexp(Nexp_),
-	MDsteps(MDsteps_),
-	trajL(trajL_),
-	stepsize(trajL/MDsteps)
-        {
+      	RealD trajL_=1.0,
+      	int Nexp_=12):
+      Nexp(Nexp_),
+      MDsteps(MDsteps_),
+      trajL(trajL_),
+      stepsize(trajL/MDsteps)
+      {
 	  // empty body constructor
-	};
+      };
 
-    };
+  };
 
     /*! @brief Class for Molecular Dynamics management */   
     template<class GaugeField, class SmearingPolicy>
-    class Integrator {
+  class Integrator {
 
-    protected:
+  protected:
 
-      typedef IntegratorParameters ParameterType;
+  	typedef IntegratorParameters ParameterType;
 
-      IntegratorParameters Params;
+  	IntegratorParameters Params;
 
-      const ActionSet<GaugeField> as;
+  	const ActionSet<GaugeField> as;
 
       int levels;              //
       double t_U;              // Track time passing on each level and for U and for P
@@ -90,14 +90,14 @@ namespace Grid{
       // Should match any legal (SU(n)) gauge field
       // Need to use this template to match Ncol to pass to SU<N> class
       template<int Ncol,class vec> void generate_momenta(Lattice< iVector< iScalar< iMatrix<vec,Ncol> >, Nd> > & P,GridParallelRNG& pRNG){
-	typedef Lattice< iScalar< iScalar< iMatrix<vec,Ncol> > > > GaugeLinkField;
-	GaugeLinkField Pmu(P._grid);
-	Pmu = zero;
-	for(int mu=0;mu<Nd;mu++){
-	  SU<Ncol>::GaussianLieAlgebraMatrix(pRNG, Pmu);
-	  PokeIndex<LorentzIndex>(P, Pmu, mu);
-	}
+      typedef Lattice< iScalar< iScalar< iMatrix<vec,Ncol> > > > GaugeLinkField;
+      GaugeLinkField Pmu(P._grid);
+      Pmu = zero;
+      for(int mu=0;mu<Nd;mu++){
+      	SU<Ncol>::GaussianLieAlgebraMatrix(pRNG, Pmu);
+      	PokeIndex<LorentzIndex>(P, Pmu, mu);
       }
+  }
 
 
       //ObserverList observers; // not yet
@@ -105,126 +105,128 @@ namespace Grid{
       //      void register_observers();
       //      void notify_observers();
 
-      void update_P(GaugeField&U, int level,double ep){
-	t_P[level]+=ep;
-	update_P(P,U,level,ep);
+  void update_P(GaugeField&U, int level, double ep){
+  	t_P[level]+=ep;
+  	update_P(P,U,level,ep);
 
-	std::cout<<GridLogIntegrator<<"["<<level<<"] P " << " dt "<< ep <<" : t_P "<< t_P[level] <<std::endl;
-      }
+  	std::cout<<GridLogIntegrator<<"["<<level<<"] P " << " dt "<< ep <<" : t_P "<< t_P[level] <<std::endl;
+  }
 
-      void update_P(GaugeField &Mom,GaugeField&U, int level,double ep){
-	for(int a=0; a<as[level].actions.size(); ++a){
-	  GaugeField force(U._grid);
-	  GaugeField& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
-	  as[level].actions.at(a)->deriv(Us,force); // deriv should not include Ta
-	  std::cout<< GridLogIntegrator << "Smearing (on/off): "<<as[level].actions.at(a)->is_smeared <<std::endl;
-	  if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
-	  force = Ta(force);
-	  std::cout<< GridLogIntegrator << "Force average: "<< norm2(force)/(U._grid->gSites()) <<std::endl;
-	  Mom = Mom - force*ep;
+  void update_P(GaugeField &Mom,GaugeField&U, int level,double ep){
+  	// input U actually not used... 
+  	for(int a=0; a<as[level].actions.size(); ++a){
+  		GaugeField force(U._grid);
+  		GaugeField& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
+  		as[level].actions.at(a)->deriv(Us,force); // deriv should NOT include Ta
+
+	  	std::cout<< GridLogIntegrator << "Smearing (on/off): "<<as[level].actions.at(a)->is_smeared <<std::endl;
+	  	if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
+	  	force = Ta(force);
+	  	std::cout<< GridLogIntegrator << "Force average: "<< norm2(force)/(U._grid->gSites()) <<std::endl;
+	  	Mom -= force*ep;
+	  }
 	}
-      }
 
-      void update_U(GaugeField&U, double ep){
-	update_U(P,U,ep);
+	void update_U(GaugeField&U, double ep){
+		update_U(P,U,ep);
 
-	t_U+=ep;
-	int fl = levels-1;
-	std::cout<< GridLogIntegrator <<"   "<<"["<<fl<<"] U " << " dt "<< ep <<" : t_U "<< t_U <<std::endl;
+		t_U+=ep;
+		int fl = levels-1;
+		std::cout<< GridLogIntegrator <<"   "<<"["<<fl<<"] U " << " dt "<< ep <<" : t_U "<< t_U <<std::endl;
 
-      }
-      void update_U(GaugeField &Mom, GaugeField&U, double ep){
+	}
+	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);
-	  Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
-	  ProjectOnGroup(Umu);
-	  PokeIndex<LorentzIndex>(U, Umu, mu);
-	}
+		for (int mu = 0; mu < Nd; mu++){
+			auto Umu=PeekIndex<LorentzIndex>(U, mu);
+			auto Pmu=PeekIndex<LorentzIndex>(Mom, mu);
+			Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
+			ProjectOnGroup(Umu);
+			PokeIndex<LorentzIndex>(U, Umu, mu);
+		}
 	// Update the smeared fields, can be implemented as observer
-	Smearer.set_GaugeField(U);
-      }
-      
-      virtual void step (GaugeField& U,int level, int first,int last)=0;
+		Smearer.set_GaugeField(U);
+	}
 
-    public:
+	virtual void step (GaugeField& U,int level, int first,int last)=0;
 
-      Integrator(GridBase* grid, 
-		 IntegratorParameters Par,
-		 ActionSet<GaugeField> & Aset,
-		 SmearingPolicy &Sm):
-          Params(Par),
-    	  as(Aset),
-	  P(grid),
-	  levels(Aset.size()),
-	  Smearer(Sm)
-      {
-	t_P.resize(levels,0.0);
-	t_U=0.0;
+public:
+
+	Integrator(GridBase* grid, 
+		IntegratorParameters Par,
+		ActionSet<GaugeField> & Aset,
+		SmearingPolicy &Sm):
+	Params(Par),
+	as(Aset),
+	P(grid),
+	levels(Aset.size()),
+	Smearer(Sm)
+	{
+		t_P.resize(levels,0.0);
+		t_U=0.0;
 	// initialization of smearer delegated outside of Integrator
-      };
-      
-      virtual ~Integrator(){}
+	};
+
+	virtual ~Integrator(){}
 
       //Initialization of momenta and actions
-      void refresh(GaugeField& U,GridParallelRNG &pRNG){
-	std::cout<<GridLogIntegrator<< "Integrator refresh\n";
-	generate_momenta(P,pRNG);
-	for(int level=0; level< as.size(); ++level){
-	  for(int actionID=0; actionID<as[level].actions.size(); ++actionID){
+	void refresh(GaugeField& U,GridParallelRNG &pRNG){
+		std::cout<<GridLogIntegrator<< "Integrator refresh\n";
+		generate_momenta(P,pRNG);
+		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
 	    // based on the boolean is_smeared in actionID
-	    GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
-	    as[level].actions.at(actionID)->refresh(Us, pRNG);
-	  }
+				GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
+				as[level].actions.at(actionID)->refresh(Us, pRNG);
+			}
+		}
 	}
-      }
 
       // Calculate action
-      RealD S(GaugeField& U){
+	RealD S(GaugeField& U){// here also U not used
 
-	LatticeComplex Hloc(U._grid);	Hloc = zero;
+		LatticeComplex Hloc(U._grid);	Hloc = zero;
 	// Momenta
-	for (int mu=0; mu <Nd; mu++){
-	  auto Pmu = PeekIndex<LorentzIndex>(P, mu);
-	  Hloc -= trace(Pmu*Pmu);
-	}
-	Complex Hsum = sum(Hloc);
-	
-	RealD H = Hsum.real();
-	RealD Hterm;
-	std::cout<<GridLogMessage << "Momentum action H_p = "<< H << "\n";
+		for (int mu=0; mu <Nd; mu++){
+			auto Pmu = PeekIndex<LorentzIndex>(P, mu);
+			Hloc -= trace(Pmu*Pmu);
+		}
+		Complex Hsum = sum(Hloc);
+
+		RealD H = Hsum.real();
+		RealD Hterm;
+		std::cout<<GridLogMessage << "Momentum action H_p = "<< H << "\n";
 
 	// Actions
-	for(int level=0; level<as.size(); ++level){
-	  for(int actionID=0; actionID<as[level].actions.size(); ++actionID){
+		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
 	    // based on the boolean is_smeared in actionID
-	    GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
-	    Hterm = as[level].actions.at(actionID)->S(Us);
-	    std::cout<<GridLogMessage << "S Level "<<level<<" term "<<actionID<<" H = "<<Hterm<<std::endl;
-	    H += Hterm;
-	  }
-	}
-	
-	return H;
-      }
+				GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
+				Hterm = as[level].actions.at(actionID)->S(Us);
+				std::cout<<GridLogMessage << "S Level "<<level<<" term "<<actionID<<" H = "<<Hterm<<std::endl;
+				H += Hterm;
+			}
+		}
 
-      void integrate(GaugeField& U){
+		return H;
+	}
+
+	void integrate(GaugeField& U){
 
 	// reset the clocks
-	t_U=0;
-	for(int level=0; level<as.size(); ++level){
-	  t_P[level]=0;
-	}	
+		t_U=0;
+		for(int level=0; level<as.size(); ++level){
+			t_P[level]=0;
+		}	
 
 	for(int step=0; step< Params.MDsteps; ++step){   // MD step
-	  int first_step = (step==0);
-	  int  last_step = (step==Params.MDsteps-1);
-	  this->step(U,0,first_step,last_step);
+		int first_step = (step==0);
+		int  last_step = (step==Params.MDsteps-1);
+		this->step(U,0,first_step,last_step);
 	}
 
 	// Check the clocks all match on all levels
@@ -237,9 +239,9 @@ namespace Grid{
 	assert(fabs(t_U-Params.trajL) < 1.0e-6);
 
 
-      }
-    };
-    
-  }
+}
+};
+
+}
 }
 #endif//INTEGRATOR_INCLUDED
diff --git a/lib/qcd/smearing/APEsmearing.h b/lib/qcd/smearing/APEsmearing.h
index 894350af..4db11338 100644
--- a/lib/qcd/smearing/APEsmearing.h
+++ b/lib/qcd/smearing/APEsmearing.h
@@ -13,10 +13,10 @@
     template <class Gimpl> 
   		class Smear_APE: public Smear<Gimpl>{
   		private:
-      const std::vector<double> rho;/*!< Array of weights */
+      	const std::vector<double> rho;/*!< Array of weights */
 
-      //This member must be private - we do not want to control from outside 
-  			std::vector<double> set_rho(const double common_rho)const {
+//This member must be private - we do not want to control from outside 
+  			std::vector<double> set_rho(const double common_rho) const {
   				std::vector<double> res;
 
   				for(int mn=0; mn<Nd*Nd; ++mn) res.push_back(common_rho);
@@ -30,7 +30,7 @@
 
 
       // Constructors and destructors
-  				Smear_APE(const std::vector<double>& rho_):rho(rho_){}
+  				Smear_APE(const std::vector<double>& rho_):rho(rho_){} // check vector size
   				Smear_APE(double rho_val):rho(set_rho(rho_val)){}
   				Smear_APE():rho(set_rho(1.0)){}
   				~Smear_APE(){}
@@ -38,7 +38,6 @@
       ///////////////////////////////////////////////////////////////////////////////
   				void smear(GaugeField& u_smr, const GaugeField& U)const{
   					GridBase *grid = U._grid;
-  					double d_rho;
   					GaugeLinkField Cup(grid), tmp_stpl(grid);
   					WilsonLoops<Gimpl> WL;
   					u_smr = zero; 
@@ -47,21 +46,20 @@
   						Cup = zero;
   						for(int nu=0; nu<Nd; ++nu){
   							if (nu != mu) {
-  								d_rho = rho[mu + Nd * nu];
-	      // get the staple in direction mu, nu
-	      WL.Staple(tmp_stpl, U, mu, nu);  //nb staple conventions of IroIro and Grid differ by a dagger
-	      Cup += tmp_stpl*d_rho;
-	  }
-	}
-	  // save the Cup link-field on the u_smr gauge-field
-	  pokeLorentz(u_smr, adj(Cup), mu); // u_smr[mu] = Cup^dag
-	}
-}
+  								// get the staple in direction mu, nu
+	      						WL.Staple(tmp_stpl, U, mu, nu);  //nb staple conventions of IroIro and Grid differ by a dagger
+	      						Cup += tmp_stpl*rho[mu + Nd * nu];
+	      					}
+	      				}
+	  					// save the Cup link-field on the u_smr gauge-field
+	  					pokeLorentz(u_smr, adj(Cup), mu); // u_smr[mu] = Cup^dag   see conventions for Staple
+	  				}
+	  			}
 
 ////////////////////////////////////////////////////////////////////////////////
-void derivative(GaugeField& SigmaTerm,
-	const GaugeField& iLambda,
-	const GaugeField& U)const{
+	  			void derivative(GaugeField& SigmaTerm,
+	  				const GaugeField& iLambda,
+	  				const GaugeField& U)const{
 
 	// Reference 
 	// Morningstar, Peardon, Phys.Rev.D69,054501(2004)
@@ -69,62 +67,61 @@ void derivative(GaugeField& SigmaTerm,
     // Computing Sigma_mu, derivative of S[fat links] with respect to the thin links
     // Output SigmaTerm
 
-	GridBase *grid = U._grid;
-	int vol = U._grid->gSites();
-	
-	WilsonLoops<Gimpl> WL;
-	GaugeLinkField staple(grid), u_tmp(grid);
-	GaugeLinkField iLambda_mu(grid), iLambda_nu(grid);
-	GaugeLinkField U_mu(grid), U_nu(grid);
-	GaugeLinkField sh_field(grid), temp_Sigma(grid);
-	Real rho_munu, rho_numu;
+	  				GridBase *grid = U._grid;
 
-	for(int mu = 0; mu < Nd; ++mu){
-		U_mu       = PeekIndex<LorentzIndex>(      U, mu);
-		iLambda_mu = PeekIndex<LorentzIndex>(iLambda, mu);
+	  				WilsonLoops<Gimpl> WL;
+	  				GaugeLinkField staple(grid), u_tmp(grid);
+	  				GaugeLinkField iLambda_mu(grid), iLambda_nu(grid);
+	  				GaugeLinkField U_mu(grid), U_nu(grid);
+	  				GaugeLinkField sh_field(grid), temp_Sigma(grid);
+	  				Real rho_munu, rho_numu;
 
-		for(int nu = 0; nu < Nd; ++nu){
-			if(nu==mu) continue;
-			U_nu       = PeekIndex<LorentzIndex>(      U, nu);
-			iLambda_nu = PeekIndex<LorentzIndex>(iLambda, nu);
+	  				for(int mu = 0; mu < Nd; ++mu){
+	  					U_mu       = peekLorentz(      U, mu);
+	  					iLambda_mu = peekLorentz(iLambda, mu);
 
-			rho_munu = rho[mu + Nd * nu];
-			rho_numu = rho[nu + Nd * mu];
+	  					for(int nu = 0; nu < Nd; ++nu){
+	  						if(nu==mu) continue;
+	  						U_nu       = peekLorentz(      U, nu);
+	  						iLambda_nu = peekLorentz(iLambda, nu);
 
-			WL.StapleUpper(staple, U, mu, nu);
+	  						rho_munu = rho[mu + Nd * nu];
+	  						rho_numu = rho[nu + Nd * mu];
 
-			temp_Sigma = -rho_numu*staple*iLambda_nu;
-	    //-r_numu*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)*Lambda_nu(x)
-			Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
+	  						WL.StapleUpper(staple, U, mu, nu);
 
-	    sh_field = Cshift(iLambda_nu, mu, 1);// general also for Gparity?
+	  						temp_Sigma = -rho_numu*staple*iLambda_nu;  //ok
+	        				//-r_numu*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)*Lambda_nu(x)
+	  						Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
 
-	    temp_Sigma = rho_numu*sh_field*staple;
-	    //r_numu*Lambda_nu(mu)*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)
-	    Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
+	    					sh_field = Cshift(iLambda_nu, mu, 1);// general also for Gparity?
 
-	    sh_field = Cshift(iLambda_mu, nu, 1);
+	    					temp_Sigma = rho_numu*sh_field*staple; //ok
+	    					//r_numu*Lambda_nu(mu)*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)
+	    					Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
 
-	    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)
-	    Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
+	    					sh_field = Cshift(iLambda_mu, nu, 1);
 
-	    staple = zero;
-	    sh_field = Cshift(U_nu, mu, 1);
+	    					temp_Sigma = -rho_munu*staple*U_nu*sh_field*adj(U_nu); //ok
+	    					//-r_munu*U_nu(x+mu)*Udag_mu(x+nu)*Lambda_mu(x+nu)*Udag_nu(x)
+	    					Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
 
-	    temp_Sigma = -rho_munu*adj(sh_field)*adj(U_mu)*iLambda_mu*U_nu;
-	    temp_Sigma += rho_numu*adj(sh_field)*adj(U_mu)*iLambda_nu*U_nu;
+	    					staple = zero;
+	    					sh_field = Cshift(U_nu, mu, 1);
 
-	    u_tmp = adj(U_nu)*iLambda_nu;
-	    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);
-	    Gimpl::AddGaugeLink(SigmaTerm, sh_field, mu);
-	    
-	}
-}
-}
-};
+	    					temp_Sigma = -rho_munu*adj(sh_field)*adj(U_mu)*iLambda_mu*U_nu;
+	    					temp_Sigma += rho_numu*adj(sh_field)*adj(U_mu)*iLambda_nu*U_nu;
+
+	    					u_tmp = adj(U_nu)*iLambda_nu;
+	    					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);
+	    					Gimpl::AddGaugeLink(SigmaTerm, sh_field, mu);
+
+	    				}
+	    			}
+	    		}
+	    	};
 
 
 
diff --git a/lib/qcd/smearing/GaugeConfiguration.h b/lib/qcd/smearing/GaugeConfiguration.h
index 717763c2..369eed72 100644
--- a/lib/qcd/smearing/GaugeConfiguration.h
+++ b/lib/qcd/smearing/GaugeConfiguration.h
@@ -17,7 +17,7 @@
       the HMC update and integrators.
       An "advanced configuration" object that can provide not only the 
       data to store the gauge configuration but also operations to manipulate
-      it like smearing.
+      it, like smearing.
       
       It stores a list of smeared configurations.
     */
@@ -57,7 +57,7 @@
 
 	}
 }
-      //====================================================================
+//====================================================================
 GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime, 
 	const GaugeField& GaugeK) const{
 	GridBase *grid = GaugeK._grid;
@@ -68,12 +68,14 @@ GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime,
 	GaugeLinkField GaugeKmu(grid), Cmu(grid);
 
 	StoutSmearing.BaseSmear(C, GaugeK);
+	SigmaK = zero;
+	iLambda = zero;
 
 	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));
+		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);
@@ -81,14 +83,14 @@ GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime,
 	StoutSmearing.derivative(SigmaK, iLambda, GaugeK);// derivative of SmearBase
 	return SigmaK;
 }
-     
+
 
 
 /*! @brief Returns smeared configuration at level 'Level' */
 const GaugeField& get_smeared_conf(int Level) const{
 	return SmearedSet[Level];
 }
-      
+
 
 //====================================================================
 void set_iLambda(GaugeLinkField& iLambda, 
@@ -101,17 +103,17 @@ void set_iLambda(GaugeLinkField& iLambda,
 	GaugeLinkField unity(grid);
 	unity=1.0;
 	
-	LatticeReal u(grid), w(grid);
+	LatticeComplex u(grid), w(grid);
 	LatticeComplex f0(grid), f1(grid), f2(grid);
-	LatticeReal xi0(grid), xi1(grid), tmp(grid);
-	LatticeReal u2(grid), w2(grid), cosw(grid);
+	LatticeComplex xi0(grid), xi1(grid), tmp(grid);
+	LatticeComplex u2(grid), w2(grid), cosw(grid);
 	LatticeComplex emiu(grid), e2iu(grid), qt(grid), fden(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 LatticeUnitReal(grid);
+	LatticeComplex LatticeUnitComplex(grid);
 
-	LatticeUnitReal = 1.0;
+	LatticeUnitComplex = 1.0;
 	
 	// Exponential
 	iQ2 = iQ * iQ;
@@ -121,44 +123,44 @@ void set_iLambda(GaugeLinkField& iLambda,
 	e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2;
 
 	// Getting B1, B2, Gamma and Lambda
+	// simplify this part, reduntant calculations in set_fj
 	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(sin(u)));
-	e2iu = toComplex(cos(2.0*u)) + timesI(toComplex(sin(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)));
+	emiu = cos(u) - timesI(sin(u));
+	e2iu = cos(2.0*u) + timesI(sin(2.0*u));
+
+	r01 = (2.0*u + timesI(2.0*(u2-w2))) * e2iu
+	+ emiu * ((16.0*u*cosw + 2.0*u*(3.0*u2+w2)*xi0) +
+		timesI(-8.0*u2*cosw + 2.0*(9.0*u2+w2)*xi0));
 	
-	r11 = (toComplex(2.0*LatticeUnitReal) + 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)));
+	r11 = (2.0*LatticeUnitComplex + timesI(4.0*u))* e2iu
+	+ emiu * ((-2.0*cosw + (3.0*u2-w2)*xi0) +
+		timesI((2.0*u*cosw + 6.0*u*xi0)));
 
 	r21 = 2.0*timesI(e2iu)
-	+ emiu * (toComplex(-3.0*u*xi0) + timesI(toComplex(cosw - 3.0*xi0)));
+	+ emiu * (-3.0*u*xi0 + timesI(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))));
+	r02 = -2.0 * e2iu + emiu * (-8.0*u2*xi0 +
+		timesI(2.0*u*(cosw + xi0 + 3.0*u2*xi1)));
 
-	r12 = emiu * (toComplex(2.0*u*xi0) + timesI(toComplex(-cosw - xi0 + 3.0*u2*xi1)));
+	r12 = emiu * (2.0*u*xi0 + timesI(-cosw - xi0 + 3.0*u2*xi1));
 
-	r22 = emiu * (toComplex(xi0) - timesI(toComplex(3.0*u*xi1)));
+	r22 = emiu * (xi0 - timesI(3.0*u*xi1));
 
-	tmp = (2.0*(9.0*u2-w2)*(9.0*u2-w2));
-	fden = toComplex(pow(tmp, -1.0));  // 1/tmp
+	fden = LatticeUnitComplex/(2.0*(9.0*u2-w2)*(9.0*u2-w2));
+	
+	b10 = 2.0 * u * r01 + (3.0* u2 - w2)*r02 - (30.0 * u2 + 2.0 * w2)*f0;
+	b11 = 2.0 * u * r11 + (3.0* u2 - w2)*r12 - (30.0 * u2 + 2.0 * w2)*f1;
+	b12 = 2.0 * u * r21 + (3.0* u2 - w2)*r22 - (30.0 * u2 + 2.0 * w2)*f2;
 
-	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;
+	b20 = r01 - (3.0*u)*r02 - (24.0*u)*f0;
+	b21 = r11 - (3.0*u)*r12 - (24.0*u)*f1;
+	b22 = r21 - (3.0*u)*r22 - (24.0*u)*f2;
 
 	b10 *= fden;
 	b11 *= fden;
@@ -166,7 +168,8 @@ void set_iLambda(GaugeLinkField& iLambda,
 	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;
@@ -180,58 +183,60 @@ void set_iLambda(GaugeLinkField& iLambda,
 	GaugeLinkField iGamma = tr1 * timesMinusI(iQ) - tr2 * iQ2 +
 	f1 * USigmap + f2 * QUS + f2 * USQ;
 
-	iLambda = Ta(iGamma);
-	
+	iLambda = Ta(timesI(iGamma));
 
 }
-     
+
 //==================================================================== 
 public:
-      GaugeField* ThinLinks;      /*!< @brief Pointer to the thin 
+    GaugeField* ThinLinks;      /*!< @brief Pointer to the thin 
 				    					links configuration */
-      
+
       /*! @brief Standard constructor */
-      SmearedConfiguration(GridCartesian * UGrid,
-      	unsigned int Nsmear, 
-      	Smear_Stout<Gimpl>& Stout):
-      smearingLevels(Nsmear),
-      StoutSmearing(Stout),
-      ThinLinks(NULL){
-      	for (unsigned int i=0; i< smearingLevels; ++i)
-      		SmearedSet.push_back(*(new GaugeField(UGrid)));
-      }
-      
+    SmearedConfiguration(GridCartesian * UGrid,
+    	unsigned int Nsmear, 
+    	Smear_Stout<Gimpl>& Stout):
+    smearingLevels(Nsmear),
+    StoutSmearing(Stout),
+    ThinLinks(NULL){
+    	for (unsigned int i=0; i< smearingLevels; ++i)
+    		SmearedSet.push_back(*(new GaugeField(UGrid)));
+    }
+
       /*! For just thin links */
-      SmearedConfiguration():
-      smearingLevels(0),
-      StoutSmearing(),
-      SmearedSet(),
-      ThinLinks(NULL){}
-      
-      
+    SmearedConfiguration():
+    smearingLevels(0),
+    StoutSmearing(),
+    SmearedSet(),
+    ThinLinks(NULL){}
+
+
       // attach the smeared routines to the thin links U and fill the smeared set
-      void set_GaugeField(GaugeField& U){ fill_smearedSet(U);}
-      
-      //====================================================================
-      void smeared_force(GaugeField& SigmaTilde) const{
-      	if (smearingLevels > 0){
-	  GaugeField force = SigmaTilde;//actually = U*SigmaTilde
-	  GaugeLinkField tmp_mu(SigmaTilde._grid);
-	  
-	for (int mu = 0; mu < Nd; mu++){
-	    // to get SigmaTilde
-	  	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);
-	  }
+    void set_GaugeField(GaugeField& U){ fill_smearedSet(U);}
+
+//====================================================================
+    void smeared_force(GaugeField& SigmaTilde) const{
+
+    	if (smearingLevels > 0){
+    		GaugeField     force(SigmaTilde._grid); 
+    		GaugeLinkField tmp_mu(SigmaTilde._grid);
+	  		force = SigmaTilde;//actually = U*SigmaTilde
+
+	  		for (int mu = 0; mu < Nd; mu++){
+		    	// to get just SigmaTilde
+	  			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);
+	  		}
 	}// if smearingLevels = 0 do nothing
 }
       //====================================================================
@@ -246,7 +251,7 @@ GaugeField& get_U(bool smeared=false) {
 	if (smeared){
 		if (smearingLevels){ 
 			RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(SmearedSet[smearingLevels-1]);
-			std::cout<< GridLogDebug << "getting U Plaq: " << impl_plaq<< std::endl;
+			std::cout<< GridLogDebug << "getting Usmr Plaq: " << impl_plaq<< std::endl;
 			return get_SmearedU();
 
 		}
diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h
index 3c323c47..511a5c29 100644
--- a/lib/qcd/smearing/StoutSmearing.h
+++ b/lib/qcd/smearing/StoutSmearing.h
@@ -12,7 +12,6 @@
     template <class Gimpl> 
   		class Smear_Stout: public Smear<Gimpl> {
   		private:
-  			const std::vector<double> d_rho;
   			const Smear < Gimpl > * SmearBase;
 
   		public:
@@ -22,58 +21,48 @@
   				static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3");
   			}
 
-      /*! Default constructor */
+			/*! Default constructor */
   			Smear_Stout(double rho = 1.0):SmearBase(new Smear_APE < Gimpl > (rho)){
   				static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3");
   			}
 
-  			~Smear_Stout(){}
+  			~Smear_Stout(){} //delete SmearBase...
 
   			void smear(GaugeField& u_smr,const GaugeField& U) const{
-
-
   				GaugeField C(U._grid);
   				GaugeLinkField tmp(U._grid), iq_mu(U._grid), Umu(U._grid);
 
   				std::cout<< GridLogDebug << "Stout smearing started\n";
 
-	//Smear the configurations
+				//Smear the configurations
   				SmearBase->smear(C, U);
+
   				for (int mu = 0; mu<Nd; mu++)
   				{
   					tmp = peekLorentz(C,mu);
   					Umu = peekLorentz(U,mu);
-		  			std::cout << "source matrix " << Umu << std::endl;	
-
 		  			iq_mu = Ta(tmp * adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper
-
 		  			exponentiate_iQ(tmp, iq_mu);  
-		  			//Debug check
-		  			GaugeLinkField check = adj(tmp) * tmp - 1.0;
-		  			std::cout << "check " << check << std::endl;
 					pokeLorentz(u_smr, tmp*Umu, mu);// u_smr = exp(iQ_mu)*U_mu
-		  		}
-
-		  		std::cout<< GridLogDebug << "Stout smearing completed\n";
+				}
+				std::cout<< GridLogDebug << "Stout smearing completed\n";
+			};
 
 
-		  	};
+			void derivative(GaugeField& SigmaTerm,
+				const GaugeField& iLambda,
+				const GaugeField& Gauge) const{
+				SmearBase->derivative(SigmaTerm, iLambda, Gauge);
+			};
 
 
-		  	void derivative(GaugeField& SigmaTerm,
-		  		const GaugeField& iLambda,
-		  		const GaugeField& Gauge) const{
-		  		SmearBase->derivative(SigmaTerm, iLambda, Gauge);
-		  	};
+			void BaseSmear(GaugeField& C,
+				const GaugeField& U) const{
+				SmearBase->smear(C, U);
+			};
 
-
-		  	void BaseSmear(GaugeField& C,
-		  		const GaugeField& U) const{
-		  		SmearBase->smear(C, U);
-		  	};
-
-		  	void exponentiate_iQ(GaugeLinkField& e_iQ,
-		  		const GaugeLinkField& iQ) const{
+			void exponentiate_iQ(GaugeLinkField& e_iQ,
+				const GaugeLinkField& iQ) const{
 		// Put this outside 
 		// only valid for SU(3) matrices
 
@@ -84,124 +73,48 @@
 		// the i sign is coming from outside
 		// input matrix is anti-hermitian NOT hermitian
 
-		  		GridBase *grid = iQ._grid;
-		  		GaugeLinkField unity(grid);
-		  		unity=1.0;
+				GridBase *grid = iQ._grid;
+				GaugeLinkField unity(grid);
+				unity=1.0;
 
-		  		GaugeLinkField iQ2(grid), iQ3(grid);
-		  		LatticeComplex u(grid), w(grid);
-		  		LatticeComplex f0(grid), f1(grid), f2(grid);
+				GaugeLinkField iQ2(grid), iQ3(grid);
+				LatticeComplex u(grid), w(grid);
+				LatticeComplex f0(grid), f1(grid), f2(grid);
 
-		  		iQ2 = iQ * iQ;
-		  		iQ3 = iQ * iQ2;
+				iQ2 = iQ * iQ;
+				iQ3 = iQ * iQ2;
 
-		  		set_uw_complex(u, w, iQ2, iQ3);
-		  		set_fj_complex(f0, f1, f2, u, w);
+				set_uw(u, w, iQ2, iQ3);
+				set_fj(f0, f1, f2, u, w);
 
-		  		std::cout << "f0 " << f0 << std::endl;
-		  		std::cout << "f1 " << f1 << std::endl;
-		  		std::cout << "f2 " << f2 << std::endl;
-	  			std::cout << "iQ " << iQ << std::endl;	
-	  			std::cout << "iQ2 " << iQ2 << std::endl;	
-
-		  		e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2;
+				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;
+			void set_uw(LatticeComplex& u, LatticeComplex& w,
+				GaugeLinkField& iQ2, GaugeLinkField& iQ3) const{
+				Complex one_over_three = 1.0/3.0;
+				Complex one_over_two = 1.0/2.0;
 
-		  		GridBase *grid = u._grid;
-		  		LatticeReal c0(grid), c1(grid), tmp(grid), c0max(grid), theta(grid);
-
-				// sign in c0 from the conventions on the Ta
-				//	c0    = - toReal(imag(trace(iQ3))) * one_over_three;
-				c0    = - toReal(real(timesMinusI(trace(iQ3)))) * one_over_three; //slow and temporary, FIX the bug in imag
-				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);
-
-			}
-
-			void set_uw_complex(LatticeComplex& u, LatticeComplex& w,
-		  		GaugeLinkField& iQ2, GaugeLinkField& iQ3) const{
-		  		Complex one_over_three = 1.0/3.0;
-		  		Complex one_over_two = 1.0/2.0;
-
-		  		GridBase *grid = u._grid;
-		  		LatticeComplex c0(grid), c1(grid), tmp(grid), c0max(grid), theta(grid);
+				GridBase *grid = u._grid;
+				LatticeComplex c0(grid), c1(grid), tmp(grid), c0max(grid), theta(grid);
 
 				// sign in c0 from the conventions on the Ta
 				c0    = - real(timesMinusI(trace(iQ3))) * one_over_three; //temporary hack 
 				c1    = - real(trace(iQ2)) * one_over_two;
 
 				//Cayley Hamilton checks to machine precision, tested
-
-				std::cout << "c0 " << c0 << std::endl;
-				std::cout << "c1 " << c1 << std::endl;
-
 				tmp   = c1 * one_over_three;
 				c0max = 2.0 * pow(tmp, 1.5);
 
-				std::cout << "c0max " << c0max << std::endl;
-				LatticeComplex tempratio = c0/c0max;
-				std::cout << "ratio c0/c0max " << tempratio << std::endl;
-				theta = acos(c0/c0max); // divide by three here, now leave as it is
-				std::cout << "theta " << theta << std::endl;
-
-				u = sqrt(tmp) * cos( theta * one_over_three);
-				w = sqrt(c1)  * sin( theta * one_over_three);
-
-				std::cout << "u " << u << std::endl;
-				std::cout << "w " << w << std::endl;
-
+				theta = acos(c0/c0max)*one_over_three; // divide by three here, now leave as it is
+				u = sqrt(tmp) * cos( theta );
+				w = sqrt(c1)  * sin( theta );
 			}
 
-
 			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 * u;
-				w2    = w * w;
-				cosw  = cos(w);
-
-				ixi0  = timesI(toComplex(xi0));
-				emiu  = toComplex(cos(u)) - timesI(toComplex(sin(u)));
-				e2iu  = toComplex(cos(2.0*u)) + timesI(toComplex(sin(2.0*u)));
-
-				h0    = e2iu * toComplex(u2 - w2) + emiu *( toComplex(8.0*u2*cosw) +
-					toComplex(2.0*u*(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);
-
-				tmp   = 9.0*u2 - w2;
-				fden  = toComplex(pow(tmp, -1.0));
-				f0    = h0 * fden;
-				f1    = h1 * fden;
-				f2    = h2 * fden;	
-
-
-			}
-
-			void set_fj_complex(LatticeComplex& f0, LatticeComplex& f1, LatticeComplex& f2,
 				const LatticeComplex& u, const LatticeComplex& w) const{
 
 				GridBase *grid = u._grid;
@@ -212,47 +125,35 @@
 				LatticeComplex unity(grid);
 				unity = 1.0;
 
-				xi0   = sin(w)/w;//func_xi0(w);
-				std::cout << "xi0 " << xi0 << std::endl;
+				xi0   = func_xi0(w);
 				u2    = u * u;
-				std::cout << "u2 " << u2 << std::endl;
 				w2    = w * w;
-				std::cout << "w2 " << w2 << std::endl;
 				cosw  = cos(w);
-				std::cout << "cosw " << cosw << std::endl;
 
 				ixi0  = timesI(xi0);
 				emiu  = cos(u)     - timesI(sin(u));
 				e2iu  = cos(2.0*u) + timesI(sin(2.0*u));
-				std::cout << "emiu " << emiu << std::endl;
-				std::cout << "e2iu " << e2iu << std::endl;
 
 				h0    = e2iu * (u2 - w2) + emiu * ( (8.0*u2*cosw) + (2.0*u*(3.0*u2 + w2)*ixi0));
 				h1    = e2iu * (2.0 * u) - emiu * ( (2.0*u*cosw) - (3.0*u2-w2)*ixi0);
 				h2    = e2iu             - emiu * ( cosw + (3.0*u)*ixi0);
 
-				std::cout << "h0 " << h0 << std::endl;
-				std::cout << "h1 " << h1 << std::endl;
-				std::cout << "h2 " << h2 << std::endl;
-
 				fden   = unity/(9.0*u2 - w2);// reals
-				std::cout << "fden " << fden << std::endl;
 				f0    = h0 * fden;
 				f1    = h1 * fden;
 				f2    = h2 * fden;	
-
 			}
 
 
 
 
-			LatticeReal func_xi0(const LatticeReal& w) const{
+			LatticeComplex func_xi0(const LatticeComplex& 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(const LatticeReal& w) const{
+			LatticeComplex func_xi1(const LatticeComplex& 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);