From 987801c86d291ebb8c4115cee87c8f2abf3f265e Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 21 Jul 2015 13:56:22 +0900 Subject: [PATCH] Merge --- lib/algorithms/approx/Remez.h | 2 +- lib/qcd/action/fermion/WilsonKernels.cc | 2 - lib/qcd/hmc/integrators/Integrator.cc | 2 +- lib/qcd/hmc/integrators/Integrator_base.h | 212 ++++++++++++++++++++++ tests/Test_multishift_sqrt.cc | 35 +++- 5 files changed, 247 insertions(+), 6 deletions(-) create mode 100644 lib/qcd/hmc/integrators/Integrator_base.h diff --git a/lib/algorithms/approx/Remez.h b/lib/algorithms/approx/Remez.h index e6a973df..336382d2 100644 --- a/lib/algorithms/approx/Remez.h +++ b/lib/algorithms/approx/Remez.h @@ -15,7 +15,7 @@ #ifndef INCLUDED_ALG_REMEZ_H #define INCLUDED_ALG_REMEZ_H -#include +#include #define JMAX 10000 //Maximum number of iterations of Newton's approximation #define SUM_MAX 10 // Maximum number of terms in exponential diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 4054280a..bdc18eee 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -1,5 +1,4 @@ #include - namespace Grid { namespace QCD { @@ -13,7 +12,6 @@ void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, vHalfSpinColourVector Uchi; int offset,local,perm, ptype; - // Xp int ss = sF; offset = st._offsets [Xp][ss]; diff --git a/lib/qcd/hmc/integrators/Integrator.cc b/lib/qcd/hmc/integrators/Integrator.cc index 8947405e..16a6513e 100644 --- a/lib/qcd/hmc/integrators/Integrator.cc +++ b/lib/qcd/hmc/integrators/Integrator.cc @@ -18,7 +18,7 @@ namespace Grid{ Pmu = zero; for(int mu=0;mu(P, Pmu, mu); } } diff --git a/lib/qcd/hmc/integrators/Integrator_base.h b/lib/qcd/hmc/integrators/Integrator_base.h new file mode 100644 index 00000000..38798d67 --- /dev/null +++ b/lib/qcd/hmc/integrators/Integrator_base.h @@ -0,0 +1,212 @@ +//-------------------------------------------------------------------- +/*! @file Integrator_base.h + * @brief Declaration of classes for the abstract Molecular Dynamics integrator + * + * @author Guido Cossu + */ +//-------------------------------------------------------------------- + +#ifndef INTEGRATOR_INCLUDED +#define INTEGRATOR_INCLUDED + +#include + +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; + + class LeapFrog; + + + 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 + //ObserverList observers; // not yet + std::unique_ptr P; + + IntegratorPolicy TheIntegrator;// contains parameters too + + + 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 with the lorentz index? + LatticeColourMatrix Umu(U._grid); + LatticeColourMatrix Pmu(U._grid); + for (int mu = 0; mu < Nd; mu++){ + Umu=PeekIndex(U, mu); + Pmu=PeekIndex(*P, mu); + Umu = expMat(Pmu, Complex(ep, 0.0))*Umu; + } + + } + + void register_observers(); + void notify_observers(); + + friend void IntegratorPolicy::step (LatticeLorentzColourMatrix& U, + int level, std::vector& clock, + Integrator* Integ); + public: + Integrator(IntegratorParameters Par, + ActionSet& Aset, std::vector Nrel_): + Params(Par),as(Aset),Nrel(Nrel_){ + assert(as.size() == Nrel.size()); + }; + + ~Integrator(){} + + + //Initialization of momenta and actions + void init(LatticeLorentzColourMatrix& U, + GridParallelRNG& pRNG){ + std::cout<< "Integrator init\n"; + if (!P) + P = new LatticeLorentzColourMatrix(U._grid); + MDutils::generate_momenta(*P,pRNG); + + for(int level=0; level< as.size(); ++level){ + for(int actionID=0; actionIDinit(U, pRNG); + } + } + } + + + + RealD S(LatticeLorentzColourMatrix& U){ + // Momenta + LatticeComplex Hloc = - trace((*P)*adj(*P)); + Complex Hsum = sum(Hloc); + + RealD H = Hsum.real(); + + // Actions + for(int level=0; levelS(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); + + }; + + class LeapFrog{ + public: + void step (LatticeLorentzColourMatrix& U, + int level, std::vector& clock, + Integrator* Integ){ + // cl : 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); + ++clock[level]; + for(int l=0; lupdate_U(U, eps); + for(int l=0; lupdate_P(U, level,eps/2); + + ++clock[level]; + for(int l=0; lupdate_P(U, level,eps); + + clock[level]+=2; + for(int l=0; l Cheby(0.1,40.0,200,InverseApproximation); + + std::cout<<"Chebuy approx vector "<