diff --git a/benchmarks/Benchmark_su3.cc b/benchmarks/Benchmark_su3.cc index 656f816a..7e5436b1 100644 --- a/benchmarks/Benchmark_su3.cc +++ b/benchmarks/Benchmark_su3.cc @@ -35,18 +35,18 @@ using namespace Grid::QCD; int main (int argc, char ** argv) { Grid_init(&argc,&argv); -#define LMIN (16) -#define LMAX (40) -#define LINC (8) +#define LMAX (32) +#define LMIN (4) +#define LINC (4) - int64_t Nloop=200; + int64_t Nloop=2000; std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); int64_t threads = GridThread::GetThreads(); std::cout< latt_size ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]}); int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3]; @@ -190,20 +191,19 @@ int main (int argc, char ** argv) LatticeColourMatrix y(&Grid); random(pRNG,y); for(int mu=0;mu<4;mu++){ - double start=usecond(); - for(int64_t i=0;i & operator = (Lattice && r) + { + _grid = r._grid; + checkerboard = r.checkerboard; + _odata =std::move(r._odata); + return *this; + } + inline Lattice & operator = (const Lattice & r){ + _grid = r._grid; + checkerboard = r.checkerboard; + _odata.resize(_grid->oSites());// essential + + parallel_for(int ss=0;ss<_grid->oSites();ss++){ + _odata[ss]=r._odata[ss]; + } + return *this; + } + + template strong_inline Lattice & operator = (const Lattice & r){ + this->checkerboard = r.checkerboard; + conformable(*this,r); + + parallel_for(int ss=0;ss<_grid->oSites();ss++){ + this->_odata[ss]=r._odata[ss]; + } + return *this; + } + virtual ~Lattice(void) = default; void reset(GridBase* grid) { @@ -277,15 +309,6 @@ public: return *this; } - template strong_inline Lattice & operator = (const Lattice & r){ - this->checkerboard = r.checkerboard; - conformable(*this,r); - - parallel_for(int ss=0;ss<_grid->oSites();ss++){ - this->_odata[ss]=r._odata[ss]; - } - return *this; - } // *=,+=,-= operators inherit behvour from correspond */+/- operation template strong_inline Lattice &operator *=(const T &r) { diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h index 95b19904..232eba42 100644 --- a/lib/qcd/hmc/integrators/Integrator.h +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -114,18 +114,26 @@ class Integrator { // input U actually not used in the fundamental case // Fundamental updates, include smearing - for (int a = 0; a < as[level].actions.size(); ++a) { + for (int a = 0; a < as[level].actions.size(); ++a) { + double start_full = usecond(); Field force(U._grid); conformable(U._grid, Mom._grid); + Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared); + double start_force = usecond(); 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 = FieldImplementation::projectForce(force); // Ta for gauge fields + double end_force = usecond(); Real force_abs = std::sqrt(norm2(force)/U._grid->gSites()); - std::cout << GridLogIntegrator << "Force average: " << force_abs << std::endl; + std::cout << GridLogIntegrator << "["< -class NoSmearing { +//trivial class for no smearing +template +class NoSmearing +{ public: INHERIT_FIELD_TYPES(Impl); - Field* ThinField; + Field *ThinField; - NoSmearing(): ThinField(NULL) {} + NoSmearing() : ThinField(NULL) {} - void set_Field(Field& U) { ThinField = &U; } + void set_Field(Field &U) { ThinField = &U; } - void smeared_force(Field&) const {} + void smeared_force(Field &) const {} - Field& get_SmearedU() { return *ThinField; } + Field &get_SmearedU() { return *ThinField; } - Field& get_U(bool smeared = false) { + Field &get_U(bool smeared = false) + { return *ThinField; } - }; /*! @@ -44,32 +47,36 @@ public: It stores a list of smeared configurations. */ template -class SmearedConfiguration { - public: +class SmearedConfiguration +{ +public: INHERIT_GIMPL_TYPES(Gimpl); - private: +private: const unsigned int smearingLevels; Smear_Stout StoutSmearing; std::vector SmearedSet; // Member functions //==================================================================== - void fill_smearedSet(GaugeField& U) { - ThinLinks = &U; // attach the smearing routine to the field U + void fill_smearedSet(GaugeField &U) + { + ThinLinks = &U; // attach the smearing routine to the field U // check the pointer is not null if (ThinLinks == NULL) std::cout << GridLogError << "[SmearedConfiguration] Error in ThinLinks pointer\n"; - if (smearingLevels > 0) { + if (smearingLevels > 0) + { std::cout << GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n"; GaugeField previous_u(ThinLinks->_grid); previous_u = *ThinLinks; - for (int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl) { + for (int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl) + { StoutSmearing.smear(SmearedSet[smearLvl], previous_u); previous_u = SmearedSet[smearLvl]; @@ -81,9 +88,10 @@ class SmearedConfiguration { } } //==================================================================== - GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime, - const GaugeField& GaugeK) const { - GridBase* grid = GaugeK._grid; + 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); @@ -94,7 +102,8 @@ class SmearedConfiguration { SigmaK = zero; iLambda = zero; - for (int mu = 0; mu < Nd; mu++) { + for (int mu = 0; mu < Nd; mu++) + { Cmu = peekLorentz(C, mu); GaugeKmu = peekLorentz(GaugeK, mu); SigmaKPrime_mu = peekLorentz(SigmaKPrime, mu); @@ -104,20 +113,22 @@ class SmearedConfiguration { pokeLorentz(iLambda, iLambda_mu, mu); } StoutSmearing.derivative(SigmaK, iLambda, - GaugeK); // derivative of SmearBase + GaugeK); // derivative of SmearBase return SigmaK; } /*! @brief Returns smeared configuration at level 'Level' */ - const GaugeField& get_smeared_conf(int Level) const { + const GaugeField &get_smeared_conf(int Level) const + { return SmearedSet[Level]; } //==================================================================== - void set_iLambda(GaugeLinkField& iLambda, GaugeLinkField& e_iQ, - const GaugeLinkField& iQ, const GaugeLinkField& Sigmap, - const GaugeLinkField& GaugeK) const { - GridBase* grid = iQ._grid; + 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; @@ -206,15 +217,15 @@ class SmearedConfiguration { } //==================================================================== - public: - GaugeField* - ThinLinks; /*!< @brief Pointer to the thin - links configuration */ +public: + GaugeField * + ThinLinks; /* Pointer to the thin links configuration */ - /*! @brief Standard constructor */ - SmearedConfiguration(GridCartesian* UGrid, unsigned int Nsmear, - Smear_Stout& Stout) - : smearingLevels(Nsmear), StoutSmearing(Stout), ThinLinks(NULL) { + /* Standard constructor */ + SmearedConfiguration(GridCartesian *UGrid, unsigned int Nsmear, + Smear_Stout &Stout) + : smearingLevels(Nsmear), StoutSmearing(Stout), ThinLinks(NULL) + { for (unsigned int i = 0; i < smearingLevels; ++i) SmearedSet.push_back(*(new GaugeField(UGrid))); } @@ -223,21 +234,29 @@ class SmearedConfiguration { SmearedConfiguration() : smearingLevels(0), StoutSmearing(), SmearedSet(), ThinLinks(NULL) {} - - // attach the smeared routines to the thin links U and fill the smeared set - void set_Field(GaugeField& U) { fill_smearedSet(U); } + void set_Field(GaugeField &U) + { + double start = usecond(); + fill_smearedSet(U); + double end = usecond(); + double time = (end - start)/ 1e3; + std::cout << GridLogMessage << "Smearing in " << time << " ms" << std::endl; + } //==================================================================== - void smeared_force(GaugeField& SigmaTilde) const { - if (smearingLevels > 0) { + void smeared_force(GaugeField &SigmaTilde) const + { + if (smearingLevels > 0) + { + double start = usecond(); GaugeField force = SigmaTilde; // actually = U*SigmaTilde GaugeLinkField tmp_mu(SigmaTilde._grid); - for (int mu = 0; mu < Nd; mu++) { + for (int mu = 0; mu < Nd; mu++) + { // to get just SigmaTilde - tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels - 1], mu)) * - peekLorentz(force, mu); + tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels - 1], mu)) * peekLorentz(force, mu); pokeLorentz(force, tmp_mu, mu); } @@ -246,33 +265,43 @@ class SmearedConfiguration { force = AnalyticSmearedForce(force, *ThinLinks); - for (int mu = 0; mu < Nd; mu++) { + 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 + double end = usecond(); + double time = (end - start)/ 1e3; + std::cout << GridLogMessage << "Smearing force in " << time << " ms" << std::endl; + } // if smearingLevels = 0 do nothing } //==================================================================== - GaugeField& get_SmearedU() { return SmearedSet[smearingLevels - 1]; } + GaugeField &get_SmearedU() { return SmearedSet[smearingLevels - 1]; } - GaugeField& get_U(bool smeared = false) { + GaugeField &get_U(bool smeared = false) + { // get the config, thin links by default - if (smeared) { - if (smearingLevels) { + if (smeared) + { + if (smearingLevels) + { RealD impl_plaq = WilsonLoops::avgPlaquette(SmearedSet[smearingLevels - 1]); std::cout << GridLogDebug << "getting Usmr Plaq: " << impl_plaq << std::endl; return get_SmearedU(); - - } else { + } + else + { RealD impl_plaq = WilsonLoops::avgPlaquette(*ThinLinks); std::cout << GridLogDebug << "getting Thin Plaq: " << impl_plaq << std::endl; return *ThinLinks; } - } else { + } + else + { RealD impl_plaq = WilsonLoops::avgPlaquette(*ThinLinks); std::cout << GridLogDebug << "getting Thin Plaq: " << impl_plaq << std::endl; diff --git a/tests/Test_dwf_mixedcg_prec.cc b/tests/Test_dwf_mixedcg_prec.cc index a53d8921..92567b6f 100644 --- a/tests/Test_dwf_mixedcg_prec.cc +++ b/tests/Test_dwf_mixedcg_prec.cc @@ -49,6 +49,8 @@ int main (int argc, char ** argv) const int Ls=8; + std::cout << GridLogMessage << "::::: NB: to enable a quick bit reproducibility check use the --checksums flag. " << std::endl; + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexD::Nsimd()),GridDefaultMpi()); GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); @@ -90,24 +92,23 @@ int main (int argc, char ** argv) SchurDiagMooeeOperator HermOpEO(Ddwf); SchurDiagMooeeOperator HermOpEO_f(Ddwf_f); - std::cout << "Starting mixed CG" << std::endl; + std::cout << GridLogMessage << "::::::::::::: Starting mixed CG" << std::endl; MixedPrecisionConjugateGradient mCG(1.0e-8, 10000, 50, FrbGrid_f, HermOpEO_f, HermOpEO); mCG(src_o,result_o); - std::cout << "Starting regular CG" << std::endl; + std::cout << GridLogMessage << "::::::::::::: Starting regular CG" << std::endl; ConjugateGradient CG(1.0e-8,10000); CG(HermOpEO,src_o,result_o_2); LatticeFermionD diff_o(FrbGrid); RealD diff = axpy_norm(diff_o, -1.0, result_o, result_o_2); - std::cout << "Diff between mixed and regular CG: " << diff << std::endl; + std::cout << GridLogMessage << "::::::::::::: Diff between mixed and regular CG: " << diff << std::endl; #ifdef HAVE_LIME if( GridCmdOptionExists(argv,argv+argc,"--checksums") ){ std::string file1("./Propagator1"); - std::string file2("./Propagator2"); emptyUserRecord record; uint32_t nersc_csum; uint32_t scidac_csuma; @@ -121,12 +122,12 @@ int main (int argc, char ** argv) BinaryIO::writeLatticeObject(result_o,file1,munge, 0, format, nersc_csum,scidac_csuma,scidac_csumb); - std::cout << " Mixed checksums "<(result_o_2,file1,munge, 0, format, nersc_csum,scidac_csuma,scidac_csumb); - std::cout << " CG checksums "<