diff --git a/lib/Make.inc b/lib/Make.inc index 85ed84c2..dfa8401b 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/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 +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.h ./qcd/hmc/integrators/Integrator_algorithm.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 +CCFILES=./qcd/hmc/integrators/Integrator.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 a61102a3..1227ad74 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -413,7 +413,8 @@ namespace QCD { #include #include #include -#include +#include +#include #include diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h index f2bb4bcb..c0f06b64 100644 --- a/lib/qcd/hmc/HMC.h +++ b/lib/qcd/hmc/HMC.h @@ -25,11 +25,11 @@ namespace Grid{ HMCparameters(); }; - template + template class HybridMonteCarlo{ const HMCparameters Params; GridSerialRNG sRNG; - Integrator& MD; + Integrator& MD; bool metropolis_test(const RealD DeltaH){ RealD rn_test; @@ -51,23 +51,21 @@ namespace Grid{ RealD evolve_step(LatticeLorentzColourMatrix& U){ - MD.init(U); // set U and initialize P and phi's - RealD H0 = MD.S(U); // current state + MD.init(U); // set U and initialize P and phi's + RealD H0 = MD.S(U); // initial state action std::cout<<"Total H before = "<< H0 << "\n"; - MD.integrate(U,0); + MD.integrate(U); - RealD H1 = MD.S(U); // updated state + RealD H1 = MD.S(U); // updated state action std::cout<<"Total H after = "<< H1 << "\n"; return (H1-H0); } - - public: HybridMonteCarlo(HMCparameters Pms, - Integrator& MolDyn): + Integrator& MolDyn): Params(Pms),MD(MolDyn){ //FIXME @@ -87,7 +85,6 @@ namespace Grid{ DeltaH = evolve_step(Uin); std::cout<< " dH = "<< DeltaH << "\n"; - } // Actual updates (evolve a copy Ucopy then copy back eventually) diff --git a/lib/qcd/hmc/integrators/Integrator_base.cc b/lib/qcd/hmc/integrators/Integrator.cc similarity index 100% rename from lib/qcd/hmc/integrators/Integrator_base.cc rename to lib/qcd/hmc/integrators/Integrator.cc diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h new file mode 100644 index 00000000..5eb319b2 --- /dev/null +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -0,0 +1,148 @@ +//-------------------------------------------------------------------- +/*! @file Integrator.h + * @brief Declaration of classes for the Molecular Dynamics integrator + * + * @author Guido Cossu + */ +//-------------------------------------------------------------------- + +#ifndef INTEGRATOR_INCLUDED +#define INTEGRATOR_INCLUDED + +class Observer; + +#include + +namespace Grid{ + namespace QCD{ + + typedef Action* ActPtr; // now force the same colours as the rest of the code + typedef std::vector ActionLevel; + typedef std::vector ActionSet; + typedef std::vector ObserverList; + + struct IntegratorParameters{ + int Nexp; + int MDsteps; // number of outer steps + RealD trajL; // trajectory length + RealD stepsize; + + IntegratorParameters(int Nexp_, + int MDsteps_, + RealD trajL_): + Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){}; + }; + + + namespace MDutils{ + void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&); + void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&); + } + + /*! @brief Class for Molecular Dynamics management */ + template< class IntegratorAlgorithm > + class Integrator{ + private: + IntegratorParameters Params; + const ActionSet as; + const std::vector Nrel; //relative step size per level + std::unique_ptr P; + GridParallelRNG pRNG; + //ObserverList observers; // not yet + + IntegratorAlgorithm TheIntegrator; + + void register_observers(); + void notify_observers(); + + void update_P(LatticeLorentzColourMatrix&U, int level,double ep){ + for(int a=0; aderiv(U,force); + *P -= force*ep; + } + } + + + void update_U(LatticeLorentzColourMatrix&U, double ep){ + //rewrite exponential to deal automatically with the lorentz index? + LatticeColourMatrix Umu(U._grid); + LatticeColourMatrix Pmu(U._grid); + for (int mu = 0; mu < Nd; mu++){ + Umu=peekLorentz(U, mu); + Pmu=peekLorentz(*P, mu); + Umu = expMat(Pmu, ep, Params.Nexp)*Umu; + pokeLorentz(U, Umu, mu); + } + + } + + + + friend void IntegratorAlgorithm::step (LatticeLorentzColourMatrix& U, + int level, std::vector& clock, + Integrator* Integ); + public: + Integrator(GridBase* grid, IntegratorParameters Par, + ActionSet& Aset, std::vector Nrel_): + Params(Par),as(Aset),Nrel(Nrel_),P(new LatticeLorentzColourMatrix(grid)),pRNG(grid){ + assert(as.size() == Nrel.size()); + pRNG.SeedRandomDevice(); + }; + + ~Integrator(){} + + + //Initialization of momenta and actions + void init(LatticeLorentzColourMatrix& U){ + std::cout<< "Integrator init\n"; + + MDutils::generate_momenta(*P,pRNG); + for(int level=0; level< as.size(); ++level){ + for(int actionID=0; actionIDinit(U, pRNG); + } + } + } + + + // Calculate action + RealD S(LatticeLorentzColourMatrix& U){ + LatticeComplex Hloc(U._grid); + Hloc = zero; + // Momenta + for (int mu=0; mu S(U); + + return H; + } + + void integrate(LatticeLorentzColourMatrix& U){ + std::vector clock; + clock.resize(as.size(),0); + for(int step=0; step< Params.MDsteps; ++step) // MD step + TheIntegrator.step(U,0,clock, (this)); + } + }; + + + + + + + + } +} +#endif//INTEGRATOR_INCLUDED diff --git a/lib/qcd/hmc/integrators/Integrator_algorithm.h b/lib/qcd/hmc/integrators/Integrator_algorithm.h new file mode 100644 index 00000000..0dce2978 --- /dev/null +++ b/lib/qcd/hmc/integrators/Integrator_algorithm.h @@ -0,0 +1,152 @@ +//-------------------------------------------------------------------- +/*! @file Integrator_algorithm.h + * @brief Declaration of classes for the Molecular Dynamics algorithms + * + * @author Guido Cossu + */ +//-------------------------------------------------------------------- + +#ifndef INTEGRATOR_ALG_INCLUDED +#define INTEGRATOR_ALG_INCLUDED + + +namespace Grid{ + namespace QCD{ + + + class MinimumNorm2{ + const double lambda = 0.1931833275037836; + public: + void step (LatticeLorentzColourMatrix& 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->Nrel[l]; + + int fin = Integ->Nrel[0]; + for(int l=1; l<=level; ++l) fin*= 2.0*Integ->Nrel[l]; + fin = 3*Integ->Params.MDsteps*fin -1; + + + for(int e=0; eNrel[level]; ++e){ + + if(clock[level] == 0){ // initial half step + Integ->update_P(U,level,lambda*eps); + ++clock[level]; + for(int l=0; lupdate_U(U,0.5*eps); + + for(int l=0; lupdate_P(U,level,(1.0-2.0*lambda)*eps); + ++clock[level]; + for(int l=0; lupdate_U(U,0.5*eps); + + for(int l=0; lupdate_P(U,level,lambda*eps); + + ++clock[level]; + for(int l=0; lupdate_P(U,level,lambda*2.0*eps); + + clock[level]+=2; + for(int l=0; l& clock, + Integrator* Integ){ + // level : current level + // fl : final level + // eps : current step size + + int fl = Integ->as.size() -1; + double eps = Integ->Params.stepsize; + + // Get current level step size + for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l]; + + int fin = 1; + for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l]; + fin = 2*Integ->Params.MDsteps*fin - 1; + + for(int e=0; eNrel[level]; ++e){ + + if(clock[level] == 0){ // initial half step + Integ->update_P(U, level,eps/2.0); + ++clock[level]; + for(int l=0; lupdate_U(U, eps); + for(int l=0; lupdate_P(U, level,eps/2.0); + + ++clock[level]; + for(int l=0; lupdate_P(U, level,eps); + + clock[level]+=2; + for(int l=0; l - -class Observer; - - -/*! @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; - - struct IntegratorParameters{ - int Nexp; - int MDsteps; // number of outer steps - RealD trajL; // trajectory length - RealD stepsize; - - IntegratorParameters(int Nexp_, - int MDsteps_, - RealD trajL_): - Nexp(Nexp_),MDsteps(MDsteps_),trajL(trajL_),stepsize(trajL/MDsteps){}; - }; - - - namespace MDutils{ - void generate_momenta(LatticeLorentzColourMatrix&,GridParallelRNG&); - void generate_momenta_su3(LatticeLorentzColourMatrix&,GridParallelRNG&); - } - - - template< class IntegratorPolicy > - class Integrator{ - private: - IntegratorParameters Params; - const ActionSet as; - const std::vector Nrel; //relative step size per level - std::unique_ptr P; - GridParallelRNG pRNG; - //ObserverList observers; // not yet - - IntegratorPolicy TheIntegrator; - - void register_observers(); - void notify_observers(); - - void update_P(LatticeLorentzColourMatrix&U, int level,double ep){ - for(int a=0; aderiv(U,force); - *P -= force*ep; - } - } - - - void update_U(LatticeLorentzColourMatrix&U, double ep){ - //rewrite exponential to deal automatically with the lorentz index? - LatticeColourMatrix Umu(U._grid); - LatticeColourMatrix Pmu(U._grid); - for (int mu = 0; mu < Nd; mu++){ - Umu=peekLorentz(U, mu); - Pmu=peekLorentz(*P, mu); - Umu = expMat(Pmu, ep, Params.Nexp)*Umu; - pokeLorentz(U, Umu, mu); - } - - } - - - - friend void IntegratorPolicy::step (LatticeLorentzColourMatrix& U, - int level, std::vector& clock, - Integrator* Integ); - public: - Integrator(GridBase* grid, IntegratorParameters Par, - ActionSet& Aset, std::vector Nrel_): - Params(Par),as(Aset),Nrel(Nrel_),P(new LatticeLorentzColourMatrix(grid)),pRNG(grid){ - assert(as.size() == Nrel.size()); - pRNG.SeedRandomDevice(); - }; - - ~Integrator(){} - - - //Initialization of momenta and actions - void init(LatticeLorentzColourMatrix& U){ - std::cout<< "Integrator init\n"; - - MDutils::generate_momenta(*P,pRNG); - for(int level=0; level< as.size(); ++level){ - for(int actionID=0; actionIDinit(U, pRNG); - } - } - } - - - // Calculate action - RealD S(LatticeLorentzColourMatrix& U){ - LatticeComplex Hloc(U._grid); - Hloc = zero; - // Momenta - for (int mu=0; mu S(U); - - return H; - } - - - - void integrate(LatticeLorentzColourMatrix& U, int level){ - std::vector clock; - clock.resize(as.size(),0); - for(int step=0; step< Params.MDsteps; ++step) // MD step - TheIntegrator.step(U,0,clock, (this)); - } - }; - - - class MinimumNorm2{ - const double lambda = 0.1931833275037836; - public: - void step (LatticeLorentzColourMatrix& 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->Nrel[l]; - - int fin = Integ->Nrel[0]; - for(int l=1; l<=level; ++l) fin*= 2.0*Integ->Nrel[l]; - fin = 3*Integ->Params.MDsteps*fin -1; - - - for(int e=0; eNrel[level]; ++e){ - - if(clock[level] == 0){ // initial half step - Integ->update_P(U,level,lambda*eps); - ++clock[level]; - for(int l=0; lupdate_U(U,0.5*eps); - - for(int l=0; lupdate_P(U,level,(1.0-2.0*lambda)*eps); - ++clock[level]; - for(int l=0; lupdate_U(U,0.5*eps); - - for(int l=0; lupdate_P(U,level,lambda*eps); - - ++clock[level]; - for(int l=0; lupdate_P(U,level,lambda*2.0*eps); - - clock[level]+=2; - for(int l=0; l& clock, - Integrator* Integ){ - // fl : final level - // eps : current step size - - int fl = Integ->as.size() -1; - double eps = Integ->Params.stepsize; - - // Get current level step size - for(int l=0; l<=level; ++l) eps/= Integ->Nrel[l]; - - int fin = 1; - for(int l=0; l<=level; ++l) fin*= Integ->Nrel[l]; - fin = 2*Integ->Params.MDsteps*fin - 1; - - for(int e=0; eNrel[level]; ++e){ - - if(clock[level] == 0){ // initial half step - Integ->update_P(U, level,eps/2.0); - ++clock[level]; - for(int l=0; lupdate_U(U, eps); - for(int l=0; lupdate_P(U, level,eps/2.0); - - ++clock[level]; - for(int l=0; lupdate_P(U, level,eps); - - clock[level]+=2; - for(int l=0; l rel ={1}; Integrator MDynamics(&Fine,MDpar, FullSet,rel);