diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h index 5a8887f2..17bee86a 100644 --- a/Grid/qcd/action/gauge/GaugeImplTypes.h +++ b/Grid/qcd/action/gauge/GaugeImplTypes.h @@ -62,7 +62,7 @@ NAMESPACE_BEGIN(Grid); // hardcodes the exponential approximation in the template //template class GaugeImplTypes { -template class GaugeImplTypes { +template > class GaugeImplTypes { public: typedef S Simd; typedef typename Simd::scalar_type scalar_type; @@ -80,7 +80,7 @@ public: typedef Lattice LinkField; typedef Lattice Field; - typedef SU Group; + //typedef SU Group; // Guido: we can probably separate the types from the HMC functions // this will create 2 kind of implementations @@ -124,53 +124,31 @@ public: for (int mu = 0; mu < Nd; mu++) { - if (isSp2n == true) - { - //const int nSp = Nrepresentation/2; - Sp::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); - } else - { - - - Group::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); - } - + Group::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR) ; Pmu = Pmu*scale; PokeIndex(P, Pmu, mu); } } - - static inline Field projectForce(Field &P) { - if (isSp2n == true) - { - P = Ta(P); - //const int nsp = Nc / 2; - - Sp::iGroupMatrix gen; - - - auto Psum = P; - - Psum = Zero(); - - for (int a = 0; a < Sp::AlgebraDimension; a++) + /* // this works + iScalar > > gen; + auto Psum = P; + Psum = Zero(); + for (int a = 0; a < Group::AlgebraDimension; a++) { - Sp::generator(a, gen); - + Group::generator(a, gen); auto coeff = 2. * trace(P * gen); Psum += coeff * gen; - } - return Psum; - - } else - { - return Ta(P); - } + + return Psum;}*/ + //this doesnt + Field ret(P.Grid()); + Group::taProj(P, ret); + return ret; } static inline void update_field(Field& P, Field& U, double ep){ @@ -181,14 +159,8 @@ public: autoView(P_v,P,AcceleratorRead); accelerator_for(ss, P.Grid()->oSites(),1,{ for (int mu = 0; mu < Nd; mu++) { - if (isSp2n == true) - { - U_v[ss](mu) = ProjectOnSpGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu)); - } else - { - U_v[ss](mu) = ProjectOnGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu)); - } - + U_v[ss](mu) = Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu); // + Group::ProjectOnGaugeGroup(U_v[ss](mu)); } }); //auto end = std::chrono::high_resolution_clock::now(); @@ -209,52 +181,22 @@ public: static inline void Project(Field &U) { - if (isSp2n == true) - { - ProjectSp2n(U); - } else - { - ProjectSUn(U); - } + Group::ProjectGn(U); } - static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) { - if (isSp2n == true) - { - //const int nSp = Nrepresentation/2; - Sp::HotConfiguration(pRNG, U); - } else - { - Group::HotConfiguration(pRNG, U); - } + Group::HotConfiguration(pRNG, U); } static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) { - if (isSp2n == true) - { - //const int nSp = Nrepresentation/2; - Sp::TepidConfiguration(pRNG, U); - } else - { - Group::TepidConfiguration(pRNG, U); - } + Group::TepidConfiguration(pRNG, U); } - static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) { - if (isSp2n == true) - { - //const int nSp = Nrepresentation/2; - Sp::ColdConfiguration(pRNG, U); - } else - { - Group::ColdConfiguration(pRNG, U); - } - + Group::ColdConfiguration(pRNG, U); } }; @@ -264,9 +206,9 @@ typedef GaugeImplTypes GimplTypesR; typedef GaugeImplTypes GimplTypesF; typedef GaugeImplTypes GimplTypesD; -typedef GaugeImplTypes SymplGimplTypesR; -typedef GaugeImplTypes SymplGimplTypesF; -typedef GaugeImplTypes SymplGimplTypesD; +typedef GaugeImplTypes> SymplGimplTypesR; +typedef GaugeImplTypes> SymplGimplTypesF; +typedef GaugeImplTypes> SymplGimplTypesD; typedef GaugeImplTypes::AdjointDimension> GimplAdjointTypesR; typedef GaugeImplTypes::AdjointDimension> GimplAdjointTypesF;