From 30c9dc473dcfa751364e6563055c2d42220b36c3 Mon Sep 17 00:00:00 2001 From: neo Date: Sat, 4 Jul 2015 02:43:14 +0900 Subject: [PATCH] More progress in the HMC construction --- lib/Make.inc | 2 +- lib/qcd/QCD.h | 10 +-- .../{gauge/GaugeActionBase.h => ActionBase.h} | 11 +-- lib/qcd/action/Actions.h | 9 ++ lib/qcd/action/gauge/WilsonGaugeAction.h | 52 +++++++----- lib/qcd/hmc/HMC.cc | 76 +---------------- lib/qcd/hmc/HMC.h | 84 +++++++++++++++++-- lib/qcd/hmc/integrators/Integrator_base.cc | 16 ++-- lib/qcd/hmc/integrators/Integrator_base.h | 67 +++++++++++---- lib/qcd/utils/SUn.h | 4 +- tests/Make.inc | 6 +- tests/Test_hmc.cc | 36 ++++++++ 12 files changed, 234 insertions(+), 139 deletions(-) rename lib/qcd/action/{gauge/GaugeActionBase.h => ActionBase.h} (65%) create mode 100644 tests/Test_hmc.cc diff --git a/lib/Make.inc b/lib/Make.inc index 18a29228..85ed84c2 100644 --- a/lib/Make.inc +++ b/lib/Make.inc @@ -1,4 +1,4 @@ -HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_exp.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_determinant.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/hmc/HMC.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Old/Tensor_poke.h ./Old/Tensor_peek.h ./Threads.h ./Grid.h ./algorithms/Preconditioner.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/AdefGeneric.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.h +HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_exp.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_determinant.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/hmc/integrators/Integrator_base.h ./qcd/hmc/HMC.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/ActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Old/Tensor_poke.h ./Old/Tensor_peek.h ./Threads.h ./Grid.h ./algorithms/Preconditioner.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/AdefGeneric.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.h CCFILES=./qcd/hmc/integrators/Integrator_base.cc ./qcd/hmc/HMC.cc ./qcd/utils/SpaceTimeGrid.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/spin/Dirac.cc ./GridInit.cc ./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index cf8ff6d8..a318141d 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -303,35 +303,35 @@ namespace QCD { const Lattice(lhs._odata[0],0))> & rhs, int i) { - pokeIndex(lhs,rhs,i); + PokeIndex(lhs,rhs,i); } template void pokeColour(Lattice &lhs, const Lattice(lhs._odata[0],0,0))> & rhs, int i,int j) { - pokeIndex(lhs,rhs,i,j); + PokeIndex(lhs,rhs,i,j); } template void pokeSpin(Lattice &lhs, const Lattice(lhs._odata[0],0))> & rhs, int i) { - pokeIndex(lhs,rhs,i); + PokeIndex(lhs,rhs,i); } template void pokeSpin(Lattice &lhs, const Lattice(lhs._odata[0],0,0))> & rhs, int i,int j) { - pokeIndex(lhs,rhs,i,j); + PokeIndex(lhs,rhs,i,j); } template void pokeLorentz(Lattice &lhs, const Lattice(lhs._odata[0],0))> & rhs, int i) { - pokeIndex(lhs,rhs,i); + PokeIndex(lhs,rhs,i); } ////////////////////////////////////////////// diff --git a/lib/qcd/action/gauge/GaugeActionBase.h b/lib/qcd/action/ActionBase.h similarity index 65% rename from lib/qcd/action/gauge/GaugeActionBase.h rename to lib/qcd/action/ActionBase.h index 55141965..e3b13f0e 100644 --- a/lib/qcd/action/gauge/GaugeActionBase.h +++ b/lib/qcd/action/ActionBase.h @@ -1,19 +1,20 @@ -#ifndef QCD_GAUGE_ACTION_BASE -#define QCD_GAUGE_ACTION_BASE +#ifndef QCD_ACTION_BASE +#define QCD_ACTION_BASE namespace Grid { namespace QCD{ template -class GaugeActionBase { // derive this from TermInAction? +class Action { public: + virtual void init(const GaugeField &U) = 0; virtual RealD S(const GaugeField &U) = 0; // evaluate the action virtual void deriv(const GaugeField &U,GaugeField & dSdU ) = 0; // evaluate the action derivative virtual void staple(const GaugeField &U,GaugeField & dSdU ) = 0; // evaluate the action derivative - virtual void refresh(const GaugeField & ) {}; + //virtual void refresh(const GaugeField & ) {} ; // Boundary conditions? // Heatbath? - virtual ~GaugeActionBase() {}; + virtual ~Action() {}; }; }} diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index 56b453c6..94809723 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -15,8 +15,17 @@ //////////////////////////////////////////// // Abstract base interface //////////////////////////////////////////// +#include + #include + +//////////////////////////////////////////// +// Gauge Actions +//////////////////////////////////////////// +#include + + //////////////////////////////////////////// // Utility functions //////////////////////////////////////////// diff --git a/lib/qcd/action/gauge/WilsonGaugeAction.h b/lib/qcd/action/gauge/WilsonGaugeAction.h index 8ecb1f81..7726db80 100644 --- a/lib/qcd/action/gauge/WilsonGaugeAction.h +++ b/lib/qcd/action/gauge/WilsonGaugeAction.h @@ -1,24 +1,38 @@ #ifndef QCD_WILSON_GAUGE_ACTION_H #define QCD_WILSON_GAUGE_ACTION_H -//////////////////////////////////////////////////////////////////////// -// Wilson Gauge Action .. should I template the Nc etc.. -//////////////////////////////////////////////////////////////////////// -template -class WilsonGaugeAction : public GaugeActionBase { - public: - virtual RealD S(const GaugeField &U) { - return WilsonLoops::sumPlaquette(U); - }; - virtual RealD deriv(GaugeField &U,GaugeField & dSdU ) { - WilsonLoops::Staple(dSdU,U,mu); - }; - virtual void staple(const MatrixField &stap,GaugeField & U,int mu ) { - WilsonLoops::Staple(stap,U,mu); - }; -}; - - -#endif +namespace Grid{ + namespace QCD{ + + //////////////////////////////////////////////////////////////////////// + // Wilson Gauge Action .. should I template the Nc etc.. + //////////////////////////////////////////////////////////////////////// + template + class WilsonGaugeAction : public Action { + private: + RealD beta; + public: + WilsonGaugeAction(RealD b):beta(b){}; + + virtual void init(const GaugeField &U) {//FIXME + }; + + virtual RealD S(const GaugeField &U) { + return WilsonLoops::sumPlaquette(U); + }; + virtual void deriv(const GaugeField &U,GaugeField & dSdU) { + //FIXME loop on directions + MatrixField dSdU_mu(U._grid); + WilsonLoops::Staple(dSdU_mu,U,0); + }; + virtual void staple(const GaugeField &stap,GaugeField & U) { + //FIXME loop on directions + MatrixField stap_mu(U._grid); + WilsonLoops::Staple(stap_mu,U,0); + }; + }; + + } +} #endif diff --git a/lib/qcd/hmc/HMC.cc b/lib/qcd/hmc/HMC.cc index 32ef8ca0..859d6646 100644 --- a/lib/qcd/hmc/HMC.cc +++ b/lib/qcd/hmc/HMC.cc @@ -18,81 +18,7 @@ namespace Grid{ } - ///////////////////////////////////////////////////////////////// - - HybridMonteCarlo::HybridMonteCarlo(GridSerialRNG& R):RNG(R){ - //FIXME - } - - - void HybridMonteCarlo::evolve(LatticeColourMatrix& Uin)const{ - Real DeltaH; - Real timer; - - // Thermalizations - for(int iter=1; iter <= Params.ThermalizationSteps; ++iter){ - std::cout << "-- # Thermalization step = "<< iter << "\n"; - - DeltaH = evolve_step(Uin); - - std::cout<< "[Timing] Trajectory time (s) : "<< timer/1000.0 << "\n"; - std::cout<< "dH = "<< DeltaH << "\n"; - - // Update matrix - //Uin = md_->get_U(); //accept every time - } - - // Actual updates - for(int iter=Params.StartingConfig; - iter < Params.Nsweeps+Params.StartingConfig; ++iter){ - std::cout << "-- # Sweep = "<< iter << "\n"; - - DeltaH = evolve_step(Uin); - - if(metropolis_test(DeltaH)) {};//Uin = md_->get_U(); - // need sync? - } - - - } - - - RealD HybridMonteCarlo::evolve_step(LatticeColourMatrix& Uin)const{ - /* - md_->init(Uin,RNG); // set U and initialize P and phi's - RealD H0 = md_->S(); // current state - std::cout<<"Total H_before = "<< H0 << "\n"; - - md_->integrator(); - - RealD H1 = md_->calc_H(); // updated state - std::cout<<"Total H_after = "<< H1 << "\n"; - - return (H1-H0); - */ - return 0; - } - - - - - bool HybridMonteCarlo::metropolis_test(const RealD DeltaH)const{ - RealD rn_test; - RealD prob = std::exp(-DeltaH); - random(RNG,rn_test); - - std::cout<< "--------------------------------------------\n"; - std::cout<< "dH = "<1.0) || (rn_test <= prob)){ // accepted - std::cout <<"-- ACCEPTED\n"; - return true; - } else { // rejected - std::cout <<"-- REJECTED\n"; - return false; - } - } + } } diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h index e25bb3ed..e4d8740c 100644 --- a/lib/qcd/hmc/HMC.h +++ b/lib/qcd/hmc/HMC.h @@ -25,21 +25,91 @@ namespace Grid{ HMCparameters(); }; - + template class HybridMonteCarlo{ const HMCparameters Params; - GridSerialRNG& RNG; - // FIXME need the integrator + GridSerialRNG sRNG; + GridParallelRNG pRNG; + std::unique_ptr< Integrator > MD; - bool metropolis_test(const RealD DeltaH)const; - RealD evolve_step(LatticeColourMatrix&)const; + + bool metropolis_test(const RealD DeltaH){ + RealD rn_test; + RealD prob = std::exp(-DeltaH); + random(sRNG,rn_test); + + std::cout<< "--------------------------------------------\n"; + std::cout<< "dH = "<1.0) || (rn_test <= prob)){ // accepted + std::cout <<"-- ACCEPTED\n"; + return true; + } else { // rejected + std::cout <<"-- REJECTED\n"; + return false; + } + } + + RealD evolve_step(LatticeColourMatrix& Uin){ + + MD->init(Uin,pRNG); // set U and initialize P and phi's + RealD H0 = MD->S(); // current state + std::cout<<"Total H_before = "<< H0 << "\n"; + + MD->integrate(0); + + RealD H1 = MD->S(); // updated state + std::cout<<"Total H_after = "<< H1 << "\n"; + + return (H1-H0); + } + public: - HybridMonteCarlo(GridSerialRNG&); + HybridMonteCarlo(Integrator& MolDyn, GridBase* grid):MD(&MolDyn),pRNG(grid){ + //FIXME + + // initialize RNGs + sRNG.SeedRandomDevice(); + pRNG.SeedRandomDevice(); + } ~HybridMonteCarlo(){}; - void evolve(LatticeColourMatrix& Uin)const; + + + void evolve(LatticeColourMatrix& Uin){ + Real DeltaH; + Real timer; + + // Thermalizations + for(int iter=1; iter <= Params.ThermalizationSteps; ++iter){ + std::cout << "-- # Thermalization step = "<< iter << "\n"; + + DeltaH = evolve_step(Uin); + + std::cout<< "[Timing] Trajectory time (s) : "<< timer/1000.0 << "\n"; + std::cout<< "dH = "<< DeltaH << "\n"; + + // Update matrix + Uin = MD->get_U(); //accept every time + } + + // Actual updates + for(int iter=Params.StartingConfig; + iter < Params.Nsweeps+Params.StartingConfig; ++iter){ + std::cout << "-- # Sweep = "<< iter << "\n"; + + DeltaH = evolve_step(Uin); + + if(metropolis_test(DeltaH)) Uin = MD->get_U(); + // need sync? + } + + + } + }; diff --git a/lib/qcd/hmc/integrators/Integrator_base.cc b/lib/qcd/hmc/integrators/Integrator_base.cc index 5fc89482..0a256444 100644 --- a/lib/qcd/hmc/integrators/Integrator_base.cc +++ b/lib/qcd/hmc/integrators/Integrator_base.cc @@ -3,18 +3,22 @@ @brief utilities for MD including funcs to generate initial HMC momentum */ - #include -static const double sq3i = 1.0/sqrt(3.0); - 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(LatticeColourMatrix& P,GridParallelRNG& pRNG){ - SU3::GaussianLieAlgebraMatrix(pRNG, P); + void MDutils::generate_momenta_su3(LatticeLorentzColourMatrix& P,GridParallelRNG& pRNG){ + LatticeColourMatrix Pmu(P._grid); + for(int mu=0;mu ActionLevel; -typedef std::vector ActionSet; -typedef std::vector ObserverList; /*! @brief Abstract base class for Molecular Dynamics management */ namespace Grid{ namespace QCD{ + + typedef Action* ActPtr; // now force the same size as the rest of the code + typedef std::vector ActionLevel; + typedef std::vector ActionSet; + typedef std::vector ObserverList; + + class Integrator2MN{ + const double lambda = 0.1931833275037836; + void step (LatticeColourMatrix&, LatticeColourMatrix&, + int, std::vector&); + + }; + + class IntegratorLeapFrog{ + void step (LatticeColourMatrix&, LatticeColourMatrix&, + int, std::vector&); + }; + + + + template< class IntegratorPolicy > class Integrator{ private: - virtual void update_P(int lv,double ep) = 0; - virtual void update_U(double ep) = 0; + int Nexp; + int MDsteps; // number of outer steps + RealD trajL; // trajectory length + RealD stepsize; + + const std::vector Nrel; // relative steps per level + const ActionSet as; + ObserverList observers; + // LatticeColourMatrix* const U; // is shared among all actions - or use a singleton... + LatticeColourMatrix P; - virtual void register_observers() = 0; - virtual void notify_observers() = 0; + IntegratorPolicy TheIntegrator;// contains parameters too + void update_P(int lv,double ep); + void update_U(double ep); + + void register_observers(); + void notify_observers(); + void integrator_step(int level ,std::vector& clock); + public: - virtual ~Integrator(){} - virtual void init(const LatticeColourMatrix&, - const GridParallelRNG& RNG)=0; - virtual double S()const =0; - virtual void integrate(int level) =0; - virtual const LatticeColourMatrix get_U() const =0; - - void generate_momenta(LatticeColourMatrix& P,const RandNum& rand); + Integrator(int Nexp_, int MDsteps_, RealD trajL_, + ActionSet& Aset, ObserverList obs):as(Aset), observers(obs){}; + ~Integrator(){} + void init(LatticeLorentzColourMatrix&, + GridParallelRNG&); + double S(); + void integrate(int level); + LatticeColourMatrix get_U(); }; namespace MDutils{ - void generate_momenta_su3(LatticeColourMatrix& P,GridParallelRNG& RNG); + void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&); + void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&); } } diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h index 9a34f3e3..70d0dcb2 100644 --- a/lib/qcd/utils/SUn.h +++ b/lib/qcd/utils/SUn.h @@ -574,14 +574,14 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g LatticeMatrix Umu(out._grid); for(int mu=0;mu latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + latt_size.resize(4); + + latt_size[0] = 8; + latt_size[1] = 8; + latt_size[2] = 8; + latt_size[3] = 8; + double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3]; + + GridCartesian Fine(latt_size,simd_layout,mpi_layout); + GridParallelRNG FineRNG(&Fine); + GridSerialRNG SerialRNG; + FineRNG.SeedRandomDevice(); + + WilsonGaugeAction Waction(6.0); + //Collect actions + ActionLevel Level; + Level.push_back(&Waction); + + + //Integrator