diff --git a/lib/qcd/smearing/APEsmearing.h b/lib/qcd/smearing/APEsmearing.h index d3fe94f6..36e14423 100644 --- a/lib/qcd/smearing/APEsmearing.h +++ b/lib/qcd/smearing/APEsmearing.h @@ -25,134 +25,130 @@ with this program; if not, write to the Free Software Foundation, Inc., See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ -/* END LEGAL */ -/*! - @brief Declaration of Smear_APE class for APE smearing -*/ + /* END LEGAL */ + /*! + @brief Declaration of Smear_APE class for APE smearing + */ #ifndef APE_SMEAR_ #define APE_SMEAR_ - namespace Grid { - namespace QCD { +NAMESPACE_BEGIN(Grid); + +/*! @brief APE type smearing of link variables. */ +template +class Smear_APE: public Smear{ +private: + const std::vector rho;/*!< Array of weights */ + + //This member must be private - we do not want to control from outside + std::vector set_rho(const double common_rho) const { + std::vector res; + + for(int mn=0; mn - class Smear_APE: public Smear{ - private: - const std::vector rho;/*!< Array of weights */ + // Constructors and destructors + Smear_APE(const std::vector& 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(){} -//This member must be private - we do not want to control from outside - std::vector set_rho(const double common_rho) const { - std::vector res; + /////////////////////////////////////////////////////////////////////////////// + void smear(GaugeField& u_smr, const GaugeField& U)const{ + GridBase *grid = U._grid; + GaugeLinkField Cup(grid), tmp_stpl(grid); + WilsonLoops WL; + u_smr = zero; - for(int mn=0; mn& 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(){} - - /////////////////////////////////////////////////////////////////////////////// - void smear(GaugeField& u_smr, const GaugeField& U)const{ - GridBase *grid = U._grid; - GaugeLinkField Cup(grid), tmp_stpl(grid); - WilsonLoops WL; - u_smr = zero; - - for(int mu=0; mu 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; + WilsonLoops 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 mu = 0; mu < Nd; ++mu){ - U_mu = peekLorentz( U, mu); - iLambda_mu = peekLorentz(iLambda, mu); + for(int mu = 0; mu < Nd; ++mu){ + U_mu = peekLorentz( U, mu); + iLambda_mu = peekLorentz(iLambda, mu); - for(int nu = 0; nu < Nd; ++nu){ - if(nu==mu) continue; - U_nu = peekLorentz( U, nu); - iLambda_nu = peekLorentz(iLambda, nu); + for(int nu = 0; nu < Nd; ++nu){ + if(nu==mu) continue; + U_nu = peekLorentz( U, nu); + iLambda_nu = peekLorentz(iLambda, nu); - rho_munu = rho[mu + Nd * nu]; - rho_numu = rho[nu + Nd * mu]; + rho_munu = rho[mu + Nd * nu]; + rho_numu = rho[nu + Nd * mu]; - WL.StapleUpper(staple, U, mu, nu); + WL.StapleUpper(staple, U, mu, nu); - 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::AddLink(SigmaTerm, temp_Sigma, mu); + 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::AddLink(SigmaTerm, temp_Sigma, mu); - sh_field = Cshift(iLambda_nu, mu, 1);// general also for Gparity? + sh_field = Cshift(iLambda_nu, mu, 1);// general also for Gparity? - 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::AddLink(SigmaTerm, temp_Sigma, mu); + 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::AddLink(SigmaTerm, temp_Sigma, mu); - sh_field = Cshift(iLambda_mu, nu, 1); + sh_field = Cshift(iLambda_mu, nu, 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::AddLink(SigmaTerm, temp_Sigma, mu); + 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::AddLink(SigmaTerm, temp_Sigma, mu); - staple = zero; - sh_field = Cshift(U_nu, mu, 1); + staple = zero; + sh_field = Cshift(U_nu, mu, 1); - 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; + 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::AddLink(SigmaTerm, sh_field, mu); + 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::AddLink(SigmaTerm, sh_field, mu); - } - } - } - }; + } + } + } +}; +NAMESPACE_END(Grid); - - }// namespace QCD -}//namespace Grid #endif