From aa52fdadccfbfdb92cd7af5b1e3ac3c107cef03a Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 30 Aug 2015 12:18:34 +0100 Subject: [PATCH] Global edit on HMC sector -- making GaugeField a template parameter and preparing to pass integrator, smearing, bc's as policy classes to hmc. Propose to unify "integrator" and integrator algorithm in a base/derived way to override step. Want to read through ForceGradient to ensure that abstraction covers the force gradient case. --- lib/Make.inc | 4 +- lib/parallelIO/NerscIO.h | 13 +- lib/qcd/QCD.h | 7 +- lib/qcd/action/Actions.h | 15 +- lib/qcd/action/gauge/WilsonGaugeAction.h | 19 +- lib/qcd/hmc/HMC.h | 46 +++-- lib/qcd/hmc/integrators/Integrator.cc | 28 --- lib/qcd/hmc/integrators/Integrator.h | 148 ++++++++------ .../hmc/integrators/Integrator_algorithm.h | 193 +++++++----------- lib/qcd/utils/WilsonLoops.h | 13 +- lib/tensors/Tensor_class.h | 3 + tests/Test_hmc_EODWFRatio.cc | 14 +- tests/Test_hmc_EOWilsonFermionGauge.cc | 22 +- tests/Test_hmc_EOWilsonRatio.cc | 21 +- tests/Test_hmc_WilsonFermionGauge.cc | 23 ++- tests/Test_hmc_WilsonGauge.cc | 22 +- tests/Test_hmc_WilsonRatio.cc | 21 +- tests/Test_rhmc_EOWilson1p1.cc | 22 +- tests/Test_rhmc_EOWilsonRatio.cc | 22 +- tests/Test_rhmc_Wilson1p1.cc | 22 +- tests/Test_rhmc_WilsonRatio.cc | 22 +- 21 files changed, 379 insertions(+), 321 deletions(-) delete mode 100644 lib/qcd/hmc/integrators/Integrator.cc diff --git a/lib/Make.inc b/lib/Make.inc index eba274ca..7bce455d 100644 --- a/lib/Make.inc +++ b/lib/Make.inc @@ -1,4 +1,4 @@ -HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Config.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/NerscIO.h ./pugixml/pugixml.h ./qcd/action/ActionBase.h ./qcd/action/ActionParams.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/EvenOddSchurDifferentiable.h ./qcd/action/pseudofermion/OneFlavourEvenOddRational.h ./qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h ./qcd/action/pseudofermion/OneFlavourRational.h ./qcd/action/pseudofermion/OneFlavourRationalRatio.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./serialisation/BinaryIO.h ./serialisation/MacroMagic.h ./serialisation/Serialisation.h ./serialisation/XmlIO.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h +HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Config.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/BinaryIO.h ./parallelIO/NerscIO.h ./pugixml/pugixml.h ./qcd/action/ActionBase.h ./qcd/action/ActionParams.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/EvenOddSchurDifferentiable.h ./qcd/action/pseudofermion/OneFlavourEvenOddRational.h ./qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h ./qcd/action/pseudofermion/OneFlavourRational.h ./qcd/action/pseudofermion/OneFlavourRationalRatio.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./serialisation/BaseIO.h ./serialisation/BinaryIO.h ./serialisation/MacroMagic.h ./serialisation/Serialisation.h ./serialisation/TextIO.h ./serialisation/XmlIO.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h -CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./Init.cc ./Log.cc ./pugixml/pugixml.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/hmc/HMC.cc ./qcd/hmc/integrators/Integrator.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc +CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./Init.cc ./Log.cc ./pugixml/pugixml.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/hmc/HMC.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc diff --git a/lib/parallelIO/NerscIO.h b/lib/parallelIO/NerscIO.h index ac17ede2..595588dd 100644 --- a/lib/parallelIO/NerscIO.h +++ b/lib/parallelIO/NerscIO.h @@ -66,11 +66,11 @@ inline void NerscGrid(GridBase *grid,NerscField &header) header.boundary[d] = std::string("PERIODIC"); } } -template +template inline void NerscStatistics(GaugeField & data,NerscField &header) { - header.link_trace=Grid::QCD::WilsonLoops::linkTrace(data); - header.plaquette =Grid::QCD::WilsonLoops::avgPlaquette(data); + header.link_trace=Grid::QCD::WilsonLoops::linkTrace(data); + header.plaquette =Grid::QCD::WilsonLoops::avgPlaquette(data); } inline void NerscMachineCharacteristics(NerscField &header) @@ -298,7 +298,6 @@ template static inline void readConfiguration(Lattice > &Umu,NerscField& header,std::string file) { typedef Lattice > GaugeField; - typedef Lattice > GaugeMatrix; GridBase *grid = Umu._grid; int offset = readHeader(file,Umu._grid,header); @@ -341,7 +340,7 @@ static inline void readConfiguration(Lattice > &Umu, assert(0); } - NerscStatistics(Umu,clone); + NerscStatistics(Umu,clone); assert(fabs(clone.plaquette -header.plaquette ) < 1.0e-5 ); assert(fabs(clone.link_trace-header.link_trace) < 1.0e-6 ); @@ -354,7 +353,7 @@ template static inline void writeConfiguration(Lattice > &Umu,std::string file, int two_row,int bits32) { typedef Lattice > GaugeField; - typedef Lattice > GaugeMatrix; + typedef iLorentzColourMatrix vobj; typedef typename vobj::scalar_object sobj; @@ -372,7 +371,7 @@ static inline void writeConfiguration(Lattice > &Umu GridBase *grid = Umu._grid; NerscGrid(grid,header); - NerscStatistics(Umu,header); + NerscStatistics(Umu,header); NerscMachineCharacteristics(header); uint32_t csum; diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index 1286c6bd..60a00243 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -62,7 +62,6 @@ namespace QCD { template using iGparitySpinColourVector = iVector, Nhs>, Ngp >; template using iGparityHalfSpinColourVector = iVector, Nhs>, Ngp >; - // Spin matrix typedef iSpinMatrix SpinMatrix; typedef iSpinMatrix SpinMatrixF; @@ -154,7 +153,7 @@ namespace QCD { typedef iHalfSpinColourVector vHalfSpinColourVectorD; // singlets - typedef iSinglet TComplex; // FIXME This is painful. Tensor singlet complex type. + typedef iSinglet TComplex; // FIXME This is painful. Tensor singlet complex type. typedef iSinglet TComplexF; // FIXME This is painful. Tensor singlet complex type. typedef iSinglet TComplexD; // FIXME This is painful. Tensor singlet complex type. @@ -162,7 +161,7 @@ namespace QCD { typedef iSinglet vTComplexF; // what if we don't know the tensor structure typedef iSinglet vTComplexD; // what if we don't know the tensor structure - typedef iSinglet TReal; // Shouldn't need these; can I make it work without? + typedef iSinglet TReal; // Shouldn't need these; can I make it work without? typedef iSinglet TRealF; // Shouldn't need these; can I make it work without? typedef iSinglet TRealD; // Shouldn't need these; can I make it work without? @@ -251,6 +250,8 @@ namespace QCD { typedef LatticeDoubleStoredColourMatrixF LatticeDoubledGaugeFieldF; typedef LatticeDoubleStoredColourMatrixD LatticeDoubledGaugeFieldD; + template using LorentzScalar = Lattice >; + // Uhgg... typing this hurt ;) // (my keyboard got burning hot when I typed this, must be the anti-Fermion) typedef Lattice LatticeStaggeredFermion; diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index d732277d..9fe0f7a7 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -16,10 +16,6 @@ #include #include -//////////////////////////////////////////// -// Gauge Actions -//////////////////////////////////////////// -#include //////////////////////////////////////////// @@ -31,6 +27,17 @@ #include //used by all wilson type fermions +//////////////////////////////////////////// +// Gauge Actions +//////////////////////////////////////////// +#include +namespace Grid { +namespace QCD { +typedef WilsonGaugeAction WilsonGaugeActionR; +typedef WilsonGaugeAction WilsonGaugeActionF; +typedef WilsonGaugeAction WilsonGaugeActionD; +}} + //////////////////////////////////////////////////////////////////////////////////////////////////// // Explicit explicit template instantiation is still required in the .cc files // diff --git a/lib/qcd/action/gauge/WilsonGaugeAction.h b/lib/qcd/action/gauge/WilsonGaugeAction.h index 9e0b79d9..7c492927 100644 --- a/lib/qcd/action/gauge/WilsonGaugeAction.h +++ b/lib/qcd/action/gauge/WilsonGaugeAction.h @@ -7,8 +7,12 @@ namespace Grid{ //////////////////////////////////////////////////////////////////////// // Wilson Gauge Action .. should I template the Nc etc.. //////////////////////////////////////////////////////////////////////// - template - class WilsonGaugeAction : public Action { + template + class WilsonGaugeAction : public Action { + public: + + typedef LorentzScalar GaugeLinkField; + private: RealD beta; public: @@ -17,23 +21,26 @@ namespace Grid{ virtual void init(const GaugeField &U, GridParallelRNG& pRNG) {}; virtual RealD S(const GaugeField &U) { - RealD plaq = WilsonLoops::avgPlaquette(U); + RealD plaq = WilsonLoops::avgPlaquette(U); std::cout<gSites(); RealD action=beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5; std::cout << GridLogMessage << "WilsonGauge action "<(U,mu); + // Staple in direction mu - WilsonLoops::Staple(dSdU_mu,U,mu); + WilsonLoops::Staple(dSdU_mu,U,mu); dSdU_mu = Ta(Umu*adj(dSdU_mu))*factor; pokeLorentz(dSdU, dSdU_mu, mu); } diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h index 5e4a22e1..e2e4f376 100644 --- a/lib/qcd/hmc/HMC.h +++ b/lib/qcd/hmc/HMC.h @@ -15,6 +15,7 @@ namespace Grid{ namespace QCD{ + struct HMCparameters{ Integer Nsweeps; /* @brief Number of sweeps in this run */ Integer TotalSweeps; /* @brief If provided, the total number of sweeps */ @@ -26,14 +27,17 @@ namespace Grid{ HMCparameters(); }; - template - class HybridMonteCarlo{ + // template + template + class HybridMonteCarlo { + private: const HMCparameters Params; - - GridSerialRNG sRNG; // Fixme: need a RNG management strategy. - - Integrator& MD; + + GridSerialRNG &sRNG; // Fixme: need a RNG management strategy. + GridParallelRNG &pRNG; // Fixme: need a RNG management strategy. + typedef Integrator IntegratorType; + IntegratorType &TheIntegrator; ///////////////////////////////////////////////////////// // Metropolis step @@ -63,16 +67,20 @@ namespace Grid{ ///////////////////////////////////////////////////////// // Evolution ///////////////////////////////////////////////////////// - RealD evolve_step(LatticeGaugeField& U){ - MD.init(U); // set U and initialize P and phi's + RealD evolve_step(GaugeField& U){ + + TheIntegrator.init(U); // set U and initialize P and phi's + + RealD H0 = TheIntegrator.S(U); // initial state action - RealD H0 = MD.S(U); // initial state action std::cout<& MolDyn): Params(Pms),MD(MolDyn) { - - //FIXME... initialize RNGs also with seed ; RNG management strategy - sRNG.SeedRandomDevice(); - + HybridMonteCarlo(HMCparameters Pms, IntegratorType &_Int, GridSerialRNG &_sRNG, GridParallelRNG &_pRNG ) : + Params(Pms), + TheIntegrator(_Int), + sRNG(_sRNG), + pRNG(_pRNG) + { } ~HybridMonteCarlo(){}; - void evolve(LatticeGaugeField& Uin){ + void evolve(GaugeField& Uin){ + Real DeltaH; // Thermalizations @@ -102,7 +112,7 @@ namespace Grid{ } // Actual updates (evolve a copy Ucopy then copy back eventually) - LatticeGaugeField Ucopy(Uin._grid); + GaugeField Ucopy(Uin._grid); for(int iter=Params.StartingConfig; iter < Params.Nsweeps+Params.StartingConfig; ++iter){ std::cout< - -namespace Grid{ - namespace QCD{ - - void MDutils::generate_momenta(LatticeLorentzColourMatrix& P,GridParallelRNG& pRNG){ - // for future support of different groups - MDutils::generate_momenta_su3(P, pRNG); - } - - void MDutils::generate_momenta_su3(LatticeLorentzColourMatrix& P,GridParallelRNG& pRNG){ - LatticeColourMatrix Pmu(P._grid); - Pmu = zero; - for(int mu=0;mu(P, Pmu, mu); - } - - } - - - } -} diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h index aa9b44bc..56d83727 100644 --- a/lib/qcd/hmc/integrators/Integrator.h +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -10,27 +10,15 @@ #ifndef INTEGRATOR_INCLUDED #define INTEGRATOR_INCLUDED -class Observer; +//class Observer; #include namespace Grid{ namespace QCD{ - typedef Action* ActPtr; // now force the same colours as the rest of the code - struct ActionLevel{ - int multiplier; - public: - std::vector actions; - explicit ActionLevel(int mul = 1):multiplier(mul){assert (mul > 0);}; - void push_back(ActPtr ptr){ - actions.push_back(ptr); - } - }; - typedef std::vector ActionSet; - typedef std::vector ObserverList; - struct IntegratorParameters{ + int Nexp; int MDsteps; // number of outer steps RealD trajL; // trajectory length @@ -39,94 +27,133 @@ namespace Grid{ IntegratorParameters(int Nexp_, int MDsteps_, RealD trajL_): - Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){}; + Nexp(Nexp_), + MDsteps(MDsteps_), + trajL(trajL_), + stepsize(trajL/MDsteps) + { + // empty body constructor + }; + }; - namespace MDutils{ - void generate_momenta(LatticeGaugeField&,GridParallelRNG&); - void generate_momenta_su3(LatticeGaugeField&,GridParallelRNG&); + // Should match any legal (SU(n)) gauge field + // Need to use this template to match Ncol to pass to SU class + template void generate_momenta(Lattice< iVector< iScalar< iMatrix >, Nd> > & P,GridParallelRNG& pRNG){ + typedef Lattice< iScalar< iScalar< iMatrix > > > GaugeLinkField; + GaugeLinkField Pmu(P._grid); + Pmu = zero; + for(int mu=0;mu::GaussianLieAlgebraMatrix(pRNG, Pmu); + PokeIndex(P, Pmu, mu); + } } + template struct ActionLevel{ + public: + + typedef Action* ActPtr; // now force the same colours as the rest of the code + + int multiplier; + + std::vector actions; + + ActionLevel(int mul = 1) : multiplier(mul) { + assert (mul > 0); + }; + + void push_back(ActPtr ptr){ + actions.push_back(ptr); + } + }; + + template using ActionSet = std::vector >; + /*! @brief Class for Molecular Dynamics management */ - template< class IntegratorAlgorithm > - class Integrator{ + template + class Integrator : public Algorithm { + private: - int levels; - std::vector t_P; - double t_U; - IntegratorParameters Params; - const ActionSet as; - std::unique_ptr P; - GridParallelRNG pRNG; + + GridParallelRNG pRNG; // Store this somewhere more sensible and pass as reference + const ActionSet as; + + int levels; // + double t_U; // Track time passing on each level and for U and for P + std::vector t_P; // + + GaugeField P; // is a pointer really necessary? //ObserverList observers; // not yet - - IntegratorAlgorithm TheIntegrator; + // typedef std::vector ObserverList; + // void register_observers(); + // void notify_observers(); - void register_observers(); - void notify_observers(); - - void update_P(LatticeGaugeField&U, int level,double ep){ + void update_P(GaugeField&U, int level,double ep){ t_P[level]+=ep; for(int a=0; aderiv(U,force); - *P -= force*ep; + P = P - force*ep; } + std::cout<(U, mu); - Pmu=PeekIndex(*P, mu); + auto Umu=PeekIndex(U, mu); + auto Pmu=PeekIndex(P, mu); Umu = expMat(Pmu, ep, Params.Nexp)*Umu; PokeIndex(U, Umu, mu); } + t_U+=ep; int fl = levels-1; std::cout<& clock, - Integrator* Integ); + friend void Algorithm::step (GaugeField& U, + int level, + std::vector& clock, + Integrator* Integ); + public: Integrator(GridBase* grid, IntegratorParameters Par, - ActionSet& Aset): + ActionSet & Aset): Params(Par), as(Aset), - P(new LatticeGaugeField(grid)), + P(grid), pRNG(grid), levels(Aset.size()) { - - std::vector seeds({1,2,3,4,5}); + std::vector seeds({1,2,3,4,5}); // Fixme; Pass it the RNG as a ref pRNG.SeedFixedIntegers(seeds); t_P.resize(levels,0.0); t_U=0.0; - }; ~Integrator(){} //Initialization of momenta and actions - void init(LatticeGaugeField& U){ + void init(GaugeField& U){ std::cout<init(U, pRNG); @@ -135,12 +162,12 @@ namespace Grid{ } // Calculate action - RealD S(LatticeGaugeField& U){ - LatticeComplex Hloc(U._grid); - Hloc = zero; + RealD S(GaugeField& U){ + + LatticeComplex Hloc(U._grid); Hloc = zero; // Momenta for (int mu=0; mu (P, mu); Hloc -= trace(Pmu*Pmu); } Complex Hsum = sum(Hloc); @@ -162,16 +189,23 @@ namespace Grid{ return H; } - void integrate(LatticeGaugeField& U){ + void integrate(GaugeField& U){ + std::vector clock; + clock.resize(as.size(),0); + + // All the clock stuff is removed if we pass first, last to the step down the way for(int step=0; step< Params.MDsteps; ++step){ // MD step - TheIntegrator.step(U,0,clock, (this)); + Algorithm::step(U,0,clock, (this)); } + + // Check the clocks all match for(int level=0; level& clock, - Integrator* Integ){ - - // level : current level - // fl : final level - // eps : current step size - - int fl = Integ->as.size() -1; - - double eps = Integ->Params.stepsize; - - for(int l=0; l<=level; ++l) eps/= 2.0*Integ->as[l].multiplier; - - // which is final half step - int fin = Integ->as[0].multiplier; - for(int l=1; l<=level; ++l) fin*= 2.0*Integ->as[l].multiplier; - fin = 3*Integ->Params.MDsteps*fin -1; - - int multiplier = Integ->as[level].multiplier; - for(int e=0; eupdate_P(U,level,lambda*eps); ++clock[level]; - } - - if(level == fl){ // lowest level - Integ->update_U(U,0.5*eps); - }else{ // recursive function call - step(U,level+1,clock, Integ); - } - - Integ->update_P(U,level,(1.0-2.0*lambda)*eps); ++clock[level]; - - if(level == fl){ // lowest level - Integ->update_U(U,0.5*eps); - }else{ // recursive function call - step(U,level+1,clock, Integ); - } - - // Handle the final half step - std::cout << GridLogMessage << " P[" <update_P(U,level,lambda*eps*mm); clock[level]+=mm; - - } - - } - - }; - - - - /* - -Chroma: Recursive leapfrog - - -00124 Real dtau = traj_length / n_steps; -00125 Real dtauby2 = dtau/2; -00126 -00127 // Its sts so: -00128 expSdt(s, dtauby2); // First half step -00129 for(int i=0; i < n_steps-1; i++) { // N-1 full steps -00130 subIntegrator(s, dtau); <--- either leapU or next integrator -00131 expSdt(s, dtau); -00132 } -00133 subIntegrator(s, dtau); // Last Full Step -00134 expSdt(s, dtauby2); // Last Half Step - + /* PAB: + * + * Recursive leapfrog; explanation of nested stepping + * * Nested 1:4; units in dt for top level integrator - * CHROMA GUIDO + * + * CHROMA IroIro * 0 1 0 * P 1/2 P 1/2 * P 1/16 P1/16 @@ -129,7 +31,7 @@ Chroma: Recursive leapfrog * U 1/8 U1/8 * P 1/16 P1/8 * P 1 P 1 - * P 1/16 * skipped + * P 1/16 * skipped --- avoids revaluating force * U 1/8 U1/8 * P 1/8 P1/8 * U 1/8 U1/8 @@ -159,13 +61,16 @@ Chroma: Recursive leapfrog * U 1/8 U1/8 * P 1/16 P1/16 * P 1/2 P 1/2 - * Total */ - class LeapFrog{ + + template class LeapFrog { public: - void step (LatticeLorentzColourMatrix& U, + + typedef LeapFrog Algorithm; + + void step (GaugeField& U, int level, std::vector& clock, - Integrator* Integ){ + Integrator * Integ){ // level : current level // fl : final level @@ -186,13 +91,7 @@ Chroma: Recursive leapfrog int first_step,last_step; - if ( level==0 ) { - first_step = (clock[level]==0); - last_step = (clock[level]==fin); - } else { - first_step = (e==0); - last_step = (e==multiplier-1); - } + first_step = (clock[level]==0); if(first_step){ // initial half step Integ->update_P(U, level,eps/2.0); ++clock[level]; @@ -204,8 +103,7 @@ Chroma: Recursive leapfrog step(U, level+1,clock, Integ); } - // Handle the final half step - std::cout << GridLogMessage << " P[" <update_P(U, level,mm*eps/2.0); clock[level]+=mm; @@ -214,6 +112,65 @@ Chroma: Recursive leapfrog } }; + template class MinimumNorm2 { + public: + typedef MinimumNorm2 Algorithm; + + private: + const double lambda = 0.1931833275037836; + public: + + void step (GaugeField& U, + int level, std::vector& clock, + Integrator* Integ){ + + // level : current level + // fl : final level + // eps : current step size + + int fl = Integ->as.size() -1; + + double eps = Integ->Params.stepsize; + + for(int l=0; l<=level; ++l) eps/= 2.0*Integ->as[l].multiplier; + + // which is final half step + int fin = Integ->as[0].multiplier; + for(int l=1; l<=level; ++l) fin*= 2.0*Integ->as[l].multiplier; + fin = 3*Integ->Params.MDsteps*fin -1; + + int multiplier = Integ->as[level].multiplier; + for(int e=0; eupdate_P(U,level,lambda*eps); ++clock[level]; + } + + if(level == fl){ // lowest level + Integ->update_U(U,0.5*eps); + }else{ // recursive function call + step(U,level+1,clock, Integ); + } + + Integ->update_P(U,level,(1.0-2.0*lambda)*eps); ++clock[level]; + + if(level == fl){ // lowest level + Integ->update_U(U,0.5*eps); + }else{ // recursive function call + step(U,level+1,clock, Integ); + } + + last_step = (clock[level]==fin); + int mm = (last_step) ? 1 : 2; + Integ->update_P(U,level,lambda*eps*mm); clock[level]+=mm; + + } + } + }; } } diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index 26975121..b166cdc8 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -4,9 +4,12 @@ namespace Grid { namespace QCD { // Common wilson loop observables -template +template class WilsonLoops { public: + + typedef LorentzScalar GaugeMat; + ////////////////////////////////////////////////// // directed plaquette oriented in mu,nu plane ////////////////////////////////////////////////// @@ -141,10 +144,10 @@ void siteRectangle(GaugeMat &plaq,const std::vector &U, const int mu, }; - typedef WilsonLoops ColourWilsonLoops; - typedef WilsonLoops U1WilsonLoops; - typedef WilsonLoops SU2WilsonLoops; - typedef WilsonLoops SU3WilsonLoops; + typedef WilsonLoops ColourWilsonLoops; + typedef WilsonLoops U1WilsonLoops; + typedef WilsonLoops SU2WilsonLoops; + typedef WilsonLoops SU3WilsonLoops; }} diff --git a/lib/tensors/Tensor_class.h b/lib/tensors/Tensor_class.h index 80247164..4f741f1d 100644 --- a/lib/tensors/Tensor_class.h +++ b/lib/tensors/Tensor_class.h @@ -23,6 +23,7 @@ template class iScalar public: vtype _internal; + typedef vtype element; typedef typename GridTypeMapper::scalar_type scalar_type; typedef typename GridTypeMapper::vector_type vector_type; typedef typename GridTypeMapper::tensor_reduced tensor_reduced_v; @@ -124,6 +125,7 @@ template class iVector public: vtype _internal[N]; + typedef vtype element; typedef typename GridTypeMapper::scalar_type scalar_type; typedef typename GridTypeMapper::vector_type vector_type; typedef typename GridTypeMapper::tensor_reduced tensor_reduced_v; @@ -219,6 +221,7 @@ template class iMatrix public: vtype _internal[N][N]; + typedef vtype element; typedef typename GridTypeMapper::scalar_type scalar_type; typedef typename GridTypeMapper::vector_type vector_type; typedef typename GridTypeMapper::tensor_reduced tensor_reduced_v; diff --git a/tests/Test_hmc_EODWFRatio.cc b/tests/Test_hmc_EODWFRatio.cc index b4b67dd9..e6a5cbb2 100644 --- a/tests/Test_hmc_EODWFRatio.cc +++ b/tests/Test_hmc_EODWFRatio.cc @@ -17,7 +17,9 @@ int main (int argc, char ** argv) GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + GridSerialRNG sRNG; GridParallelRNG pRNG(UGrid); + sRNG.SeedRandomDevice(); pRNG.SeedRandomDevice(); LatticeLorentzColourMatrix U(UGrid); @@ -25,7 +27,7 @@ int main (int argc, char ** argv) SU3::HotConfiguration(pRNG, U); // simplify template declaration? Strip the lorentz from the second template - WilsonGaugeAction Waction(5.6); + WilsonGaugeActionR Waction(5.6); Real mass=0.04; Real pv =1.0; @@ -38,21 +40,21 @@ int main (int argc, char ** argv) TwoFlavourEvenOddRatioPseudoFermionAction Nf2(NumOp, DenOp,CG,CG); //Collect actions - ActionLevel Level1; + ActionLevel Level1; Level1.push_back(&Nf2); Level1.push_back(&Waction); - ActionSet FullSet; + ActionSet FullSet; FullSet.push_back(Level1); // Create integrator // typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm - typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm + typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm IntegratorParameters MDpar(12,20,1.0); - Integrator MDynamics(UGrid,MDpar, FullSet); + Integrator MDynamics(UGrid,MDpar, FullSet); // Create HMC HMCparameters HMCpar; - HybridMonteCarlo HMC(HMCpar, MDynamics); + HybridMonteCarlo HMC(HMCpar, MDynamics,sRNG,pRNG); // Run it HMC.evolve(U); diff --git a/tests/Test_hmc_EOWilsonFermionGauge.cc b/tests/Test_hmc_EOWilsonFermionGauge.cc index ebe6d63e..0d20e207 100644 --- a/tests/Test_hmc_EOWilsonFermionGauge.cc +++ b/tests/Test_hmc_EOWilsonFermionGauge.cc @@ -15,15 +15,21 @@ int main (int argc, char ** argv) GridCartesian Fine(latt_size,simd_layout,mpi_layout); GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout); - GridParallelRNG pRNG(&Fine); + std::vector seeds({6,7,8,80}); + GridParallelRNG pRNG(&Fine); pRNG.SeedFixedIntegers(seeds); + + std::vector seedsS({1,2,3,4}); + GridSerialRNG sRNG; + sRNG.SeedFixedIntegers(seedsS); + LatticeLorentzColourMatrix U(&Fine); SU3::HotConfiguration(pRNG, U); // simplify template declaration? Strip the lorentz from the second template - WilsonGaugeAction Waction(5.6); + WilsonGaugeActionR Waction(5.6); Real mass=-0.77; WilsonFermionR FermOp(U,Fine,RBFine,mass); @@ -33,25 +39,25 @@ int main (int argc, char ** argv) TwoFlavourEvenOddPseudoFermionAction WilsonNf2(FermOp,CG,CG); //Collect actions - ActionLevel Level1(1); - ActionLevel Level2(4); + ActionLevel Level1(1); + ActionLevel Level2(4); Level1.push_back(&WilsonNf2); Level2.push_back(&Waction); - ActionSet FullSet; + ActionSet FullSet; FullSet.push_back(Level1); FullSet.push_back(Level2); // Create integrator // typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm // IntegratorParameters MDpar(12,40,1.0); - typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm + typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm IntegratorParameters MDpar(12,10,1.0); - Integrator MDynamics(&Fine,MDpar, FullSet); + Integrator MDynamics(&Fine,MDpar, FullSet); // Create HMC HMCparameters HMCpar; - HybridMonteCarlo HMC(HMCpar, MDynamics); + HybridMonteCarlo HMC(HMCpar, MDynamics,sRNG,pRNG); HMC.evolve(U); diff --git a/tests/Test_hmc_EOWilsonRatio.cc b/tests/Test_hmc_EOWilsonRatio.cc index 6a4e2bbb..9f47770d 100644 --- a/tests/Test_hmc_EOWilsonRatio.cc +++ b/tests/Test_hmc_EOWilsonRatio.cc @@ -15,14 +15,21 @@ int main (int argc, char ** argv) GridCartesian Fine(latt_size,simd_layout,mpi_layout); GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout); + + std::vector seeds({6,7,8,80}); GridParallelRNG pRNG(&Fine); - pRNG.SeedRandomDevice(); + pRNG.SeedFixedIntegers(seeds); + + std::vector seedsS({1,2,3,4}); + GridSerialRNG sRNG; + sRNG.SeedFixedIntegers(seedsS); + LatticeLorentzColourMatrix U(&Fine); SU3::HotConfiguration(pRNG, U); // simplify template declaration? Strip the lorentz from the second template - WilsonGaugeAction Waction(5.6); + WilsonGaugeActionR Waction(5.6); Real mass=-0.77; Real pv =0.0; @@ -34,22 +41,22 @@ int main (int argc, char ** argv) TwoFlavourEvenOddRatioPseudoFermionAction WilsonNf2(NumOp, DenOp,CG,CG); //Collect actions - ActionLevel Level1; + ActionLevel Level1; Level1.push_back(&WilsonNf2); Level1.push_back(&Waction); - ActionSet FullSet; + ActionSet FullSet; FullSet.push_back(Level1); // Create integrator - typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm + typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm // typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm IntegratorParameters MDpar(12,20,1.0); - Integrator MDynamics(&Fine,MDpar, FullSet); + Integrator MDynamics(&Fine,MDpar, FullSet); // Create HMC HMCparameters HMCpar; - HybridMonteCarlo HMC(HMCpar, MDynamics); + HybridMonteCarlo HMC(HMCpar, MDynamics,sRNG,pRNG); HMC.evolve(U); diff --git a/tests/Test_hmc_WilsonFermionGauge.cc b/tests/Test_hmc_WilsonFermionGauge.cc index 75c78e99..93412536 100644 --- a/tests/Test_hmc_WilsonFermionGauge.cc +++ b/tests/Test_hmc_WilsonFermionGauge.cc @@ -15,16 +15,21 @@ int main (int argc, char ** argv) GridCartesian Fine(latt_size,simd_layout,mpi_layout); GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout); - GridParallelRNG pRNG(&Fine); - std::vector seeds({5,6,7,8}); + std::vector seeds({6,7,8,80}); + GridParallelRNG pRNG(&Fine); pRNG.SeedFixedIntegers(seeds); + + std::vector seedsS({1,2,3,4}); + GridSerialRNG sRNG; + sRNG.SeedFixedIntegers(seedsS); + LatticeLorentzColourMatrix U(&Fine); SU3::HotConfiguration(pRNG, U); // simplify template declaration? Strip the lorentz from the second template - WilsonGaugeAction Waction(5.6); + WilsonGaugeActionR Waction(5.6); Real mass=-0.77; WilsonFermionR FermOp(U,Fine,RBFine,mass); @@ -34,25 +39,25 @@ int main (int argc, char ** argv) TwoFlavourPseudoFermionAction WilsonNf2(FermOp,CG,CG); //Collect actions - ActionLevel Level1(1); - ActionLevel Level2(4); + ActionLevel Level1(1); + ActionLevel Level2(4); Level1.push_back(&WilsonNf2); Level1.push_back(&Waction); // Level1.push_back(&Waction); - ActionSet FullSet; + ActionSet FullSet; FullSet.push_back(Level1); // FullSet.push_back(Level2); // Create integrator - typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm + typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm // typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm IntegratorParameters MDpar(12,20,1.0); - Integrator MDynamics(&Fine,MDpar, FullSet); + Integrator MDynamics(&Fine,MDpar, FullSet); // Create HMC HMCparameters HMCpar; - HybridMonteCarlo HMC(HMCpar, MDynamics); + HybridMonteCarlo HMC(HMCpar, MDynamics,sRNG,pRNG); HMC.evolve(U); diff --git a/tests/Test_hmc_WilsonGauge.cc b/tests/Test_hmc_WilsonGauge.cc index 725a114c..448b11cd 100644 --- a/tests/Test_hmc_WilsonGauge.cc +++ b/tests/Test_hmc_WilsonGauge.cc @@ -22,31 +22,37 @@ int main (int argc, char ** argv) double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3]; GridCartesian Fine(latt_size,simd_layout,mpi_layout); + + + std::vector seeds({6,7,8,80}); GridParallelRNG pRNG(&Fine); - std::vector seeds({1,2,3,4,5,6,7,8}); pRNG.SeedFixedIntegers(seeds); + + std::vector seedsS({1,2,3,4}); + GridSerialRNG sRNG; + sRNG.SeedFixedIntegers(seedsS); + LatticeGaugeField U(&Fine); SU3::HotConfiguration(pRNG, U); - // simplify template declaration? Strip the lorentz from the second template - WilsonGaugeAction Waction(6.0); + WilsonGaugeActionR Waction(6.0); //Collect actions - ActionLevel Level1; + ActionLevel Level1; Level1.push_back(&Waction); - ActionSet FullSet; + ActionSet FullSet; FullSet.push_back(Level1); // Create integrator - typedef MinimumNorm2 IntegratorAlgorithm;// change here to modify the algorithm + typedef MinimumNorm2 IntegratorAlgorithm;// change here to modify the algorithm IntegratorParameters MDpar(12,20,1.0); - Integrator MDynamics(&Fine,MDpar, FullSet); + Integrator MDynamics(&Fine,MDpar, FullSet); // Create HMC HMCparameters HMCpar; - HybridMonteCarlo HMC(HMCpar, MDynamics); + HybridMonteCarlo HMC(HMCpar, MDynamics, sRNG, pRNG); HMC.evolve(U); diff --git a/tests/Test_hmc_WilsonRatio.cc b/tests/Test_hmc_WilsonRatio.cc index c7cdec39..16bee903 100644 --- a/tests/Test_hmc_WilsonRatio.cc +++ b/tests/Test_hmc_WilsonRatio.cc @@ -15,14 +15,21 @@ int main (int argc, char ** argv) GridCartesian Fine(latt_size,simd_layout,mpi_layout); GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout); + + std::vector seeds({6,7,8,80}); GridParallelRNG pRNG(&Fine); - pRNG.SeedRandomDevice(); + pRNG.SeedFixedIntegers(seeds); + + std::vector seedsS({1,2,3,4}); + GridSerialRNG sRNG; + sRNG.SeedFixedIntegers(seedsS); + LatticeLorentzColourMatrix U(&Fine); SU3::HotConfiguration(pRNG, U); // simplify template declaration? Strip the lorentz from the second template - WilsonGaugeAction Waction(5.6); + WilsonGaugeActionR Waction(5.6); Real mass=-0.77; Real pv =0.0; @@ -34,22 +41,22 @@ int main (int argc, char ** argv) TwoFlavourRatioPseudoFermionAction WilsonNf2(NumOp,FermOp,CG,CG); //Collect actions - ActionLevel Level1; + ActionLevel Level1; Level1.push_back(&WilsonNf2); Level1.push_back(&Waction); - ActionSet FullSet; + ActionSet FullSet; FullSet.push_back(Level1); // Create integrator - typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm + typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm // typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm IntegratorParameters MDpar(12,20,1.0); - Integrator MDynamics(&Fine,MDpar, FullSet); + Integrator MDynamics(&Fine,MDpar, FullSet); // Create HMC HMCparameters HMCpar; - HybridMonteCarlo HMC(HMCpar, MDynamics); + HybridMonteCarlo HMC(HMCpar, MDynamics, sRNG, pRNG); HMC.evolve(U); diff --git a/tests/Test_rhmc_EOWilson1p1.cc b/tests/Test_rhmc_EOWilson1p1.cc index 96be8b20..86f2bfe2 100644 --- a/tests/Test_rhmc_EOWilson1p1.cc +++ b/tests/Test_rhmc_EOWilson1p1.cc @@ -15,14 +15,22 @@ int main (int argc, char ** argv) GridCartesian Fine(latt_size,simd_layout,mpi_layout); GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout); + + + std::vector seeds({6,7,8,80}); GridParallelRNG pRNG(&Fine); - pRNG.SeedRandomDevice(); + pRNG.SeedFixedIntegers(seeds); + + std::vector seedsS({1,2,3,4}); + GridSerialRNG sRNG; + sRNG.SeedFixedIntegers(seedsS); + LatticeLorentzColourMatrix U(&Fine); SU3::HotConfiguration(pRNG, U); // simplify template declaration? Strip the lorentz from the second template - WilsonGaugeAction Waction(5.6); + WilsonGaugeActionR Waction(5.6); Real mass=-0.77; WilsonFermionR FermOp(U,Fine,RBFine,mass); @@ -33,23 +41,23 @@ int main (int argc, char ** argv) OneFlavourEvenOddRationalPseudoFermionAction WilsonNf1b(FermOp,Params); //Collect actions - ActionLevel Level1; + ActionLevel Level1; Level1.push_back(&WilsonNf1a); Level1.push_back(&WilsonNf1b); Level1.push_back(&Waction); - ActionSet FullSet; + ActionSet FullSet; FullSet.push_back(Level1); // Create integrator - typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm + typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm IntegratorParameters MDpar(12,20,1.0); - Integrator MDynamics(&Fine,MDpar, FullSet); + Integrator MDynamics(&Fine,MDpar, FullSet); // Create HMC HMCparameters HMCpar; - HybridMonteCarlo HMC(HMCpar, MDynamics); + HybridMonteCarlo HMC(HMCpar, MDynamics, sRNG, pRNG); HMC.evolve(U); } diff --git a/tests/Test_rhmc_EOWilsonRatio.cc b/tests/Test_rhmc_EOWilsonRatio.cc index fcb66498..c28d30a7 100644 --- a/tests/Test_rhmc_EOWilsonRatio.cc +++ b/tests/Test_rhmc_EOWilsonRatio.cc @@ -15,14 +15,22 @@ int main (int argc, char ** argv) GridCartesian Fine(latt_size,simd_layout,mpi_layout); GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout); + + + std::vector seeds({6,7,8,80}); GridParallelRNG pRNG(&Fine); - pRNG.SeedRandomDevice(); + pRNG.SeedFixedIntegers(seeds); + + std::vector seedsS({1,2,3,4}); + GridSerialRNG sRNG; + sRNG.SeedFixedIntegers(seedsS); + LatticeLorentzColourMatrix U(&Fine); SU3::HotConfiguration(pRNG, U); // simplify template declaration? Strip the lorentz from the second template - WilsonGaugeAction Waction(5.6); + WilsonGaugeActionR Waction(5.6); Real mass=-0.77; Real pv =0.0; @@ -35,23 +43,23 @@ int main (int argc, char ** argv) OneFlavourEvenOddRatioRationalPseudoFermionAction WilsonNf1b(NumOp,DenOp,Params); //Collect actions - ActionLevel Level1; + ActionLevel Level1; Level1.push_back(&WilsonNf1a); Level1.push_back(&WilsonNf1b); Level1.push_back(&Waction); - ActionSet FullSet; + ActionSet FullSet; FullSet.push_back(Level1); // Create integrator - typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm + typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm // typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm IntegratorParameters MDpar(12,20,1.0); - Integrator MDynamics(&Fine,MDpar, FullSet); + Integrator MDynamics(&Fine,MDpar, FullSet); // Create HMC HMCparameters HMCpar; - HybridMonteCarlo HMC(HMCpar, MDynamics); + HybridMonteCarlo HMC(HMCpar, MDynamics, sRNG, pRNG); HMC.evolve(U); diff --git a/tests/Test_rhmc_Wilson1p1.cc b/tests/Test_rhmc_Wilson1p1.cc index f44e121a..bee457f6 100644 --- a/tests/Test_rhmc_Wilson1p1.cc +++ b/tests/Test_rhmc_Wilson1p1.cc @@ -15,14 +15,22 @@ int main (int argc, char ** argv) GridCartesian Fine(latt_size,simd_layout,mpi_layout); GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout); + + + std::vector seeds({6,7,8,80}); GridParallelRNG pRNG(&Fine); - pRNG.SeedRandomDevice(); + pRNG.SeedFixedIntegers(seeds); + + std::vector seedsS({1,2,3,4}); + GridSerialRNG sRNG; + sRNG.SeedFixedIntegers(seedsS); + LatticeLorentzColourMatrix U(&Fine); SU3::HotConfiguration(pRNG, U); // simplify template declaration? Strip the lorentz from the second template - WilsonGaugeAction Waction(5.6); + WilsonGaugeActionR Waction(5.6); Real mass=-0.77; WilsonFermionR FermOp(U,Fine,RBFine,mass); @@ -33,23 +41,23 @@ int main (int argc, char ** argv) OneFlavourRationalPseudoFermionAction WilsonNf1b(FermOp,Params); //Collect actions - ActionLevel Level1; + ActionLevel Level1; Level1.push_back(&WilsonNf1a); Level1.push_back(&WilsonNf1b); Level1.push_back(&Waction); - ActionSet FullSet; + ActionSet FullSet; FullSet.push_back(Level1); // Create integrator - typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm + typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm IntegratorParameters MDpar(12,20,1.0); - Integrator MDynamics(&Fine,MDpar, FullSet); + Integrator MDynamics(&Fine,MDpar, FullSet); // Create HMC HMCparameters HMCpar; - HybridMonteCarlo HMC(HMCpar, MDynamics); + HybridMonteCarlo HMC(HMCpar, MDynamics, sRNG, pRNG); HMC.evolve(U); diff --git a/tests/Test_rhmc_WilsonRatio.cc b/tests/Test_rhmc_WilsonRatio.cc index 9b27fa70..a47d1eaa 100644 --- a/tests/Test_rhmc_WilsonRatio.cc +++ b/tests/Test_rhmc_WilsonRatio.cc @@ -15,14 +15,22 @@ int main (int argc, char ** argv) GridCartesian Fine(latt_size,simd_layout,mpi_layout); GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout); + + + std::vector seeds({6,7,8,80}); GridParallelRNG pRNG(&Fine); - pRNG.SeedRandomDevice(); + pRNG.SeedFixedIntegers(seeds); + + std::vector seedsS({1,2,3,4}); + GridSerialRNG sRNG; + sRNG.SeedFixedIntegers(seedsS); + LatticeLorentzColourMatrix U(&Fine); SU3::HotConfiguration(pRNG, U); // simplify template declaration? Strip the lorentz from the second template - WilsonGaugeAction Waction(5.6); + WilsonGaugeActionR Waction(5.6); Real mass=-0.77; Real pv =0.0; @@ -35,23 +43,23 @@ int main (int argc, char ** argv) OneFlavourRatioRationalPseudoFermionAction WilsonNf1b(NumOp,DenOp,Params); //Collect actions - ActionLevel Level1; + ActionLevel Level1; Level1.push_back(&WilsonNf1a); Level1.push_back(&WilsonNf1b); Level1.push_back(&Waction); - ActionSet FullSet; + ActionSet FullSet; FullSet.push_back(Level1); // Create integrator - typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm + typedef MinimumNorm2 IntegratorAlgorithm;// change here to change the algorithm // typedef LeapFrog IntegratorAlgorithm;// change here to change the algorithm IntegratorParameters MDpar(12,20,1.0); - Integrator MDynamics(&Fine,MDpar, FullSet); + Integrator MDynamics(&Fine,MDpar, FullSet); // Create HMC HMCparameters HMCpar; - HybridMonteCarlo HMC(HMCpar, MDynamics); + HybridMonteCarlo HMC(HMCpar, MDynamics, sRNG, pRNG); HMC.evolve(U);