portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Recovering mixed precision CG for Laplace

Browse Source 
Checking in to move to aurora
This commit is contained in:
Chulwoo Jung 2023-12-12 15:32:00 -05:00
parent 7af6022a2a
commit 076580c232
4 changed files with 147 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
Grid/qcd/hmc/integrators/Integrator.h
View File

@ -223,9 +223,11 @@ public:
as[level].apply(update_P_hireps, Representations, Mom, U, ep);
}
void implicit_update_P(Field& U, int level, double ep, bool intermediate = false) {
void implicit_update_P(Field& U, int level, double ep, double ep1, bool intermediate = false) {
t_P[level] += ep;
double ep2= ep-ep1;
std::cout << GridLogIntegrator << "[" << level << "] P "
<< " dt " << ep << " : t_P " << t_P[level] << std::endl;
std::cout << GridLogIntegrator << "U before implicit_update_P: " << std::sqrt(norm2(U)) << std::endl;
@ -267,7 +269,7 @@ public:
// std::cout << GridLogIntegrator << "MomDer1 implicit_update_P: " << std::sqrt(norm2(MomDer1)) << std::endl;
// Auxiliary fields
P.update_auxiliary_momenta(ep*0.5 );
P.update_auxiliary_momenta(ep1);
P.AuxiliaryFieldsDerivative(AuxDer);
Msum += AuxDer;
@ -285,7 +287,8 @@ public:
std::cout << GridLogIntegrator << "|Force| laplacian site average: " << force_abs
<< std::endl;
NewMom = P.Mom - ep* 0.5 * HMC_MOMENTUM_DENOMINATOR * (2.0*Msum + factor*MomDer + MomDer1);// simplify
// NewMom = P.Mom - ep* 0.5 * HMC_MOMENTUM_DENOMINATOR * (2.0*Msum + factor*MomDer + MomDer1);// simplify
NewMom = P.Mom - HMC_MOMENTUM_DENOMINATOR * (ep*Msum + ep1* factor*MomDer + ep2* MomDer1);// simplify
diff = NewMom - OldMom;
counter++;
RelativeError = std::sqrt(norm2(diff))/std::sqrt(norm2(NewMom));
@ -297,7 +300,11 @@ public:
std::cout << GridLogIntegrator << "NewMom implicit_update_P: " << std::sqrt(norm2(NewMom)) << std::endl;
// update the auxiliary fields momenta
P.update_auxiliary_momenta(ep*0.5 );
P.update_auxiliary_momenta(ep2);
}
void implicit_update_P(Field& U, int level, double ep, bool intermediate = false) {
implicit_update_P( U, level, ep, ep*0.5, intermediate );
}
void update_U(Field& U, double ep)

113
Grid/qcd/hmc/integrators/Integrator_algorithm.h
View File

@ -394,19 +394,11 @@ class ImplicitMinimumNorm2 : public Integrator<FieldImplementation, SmearingPoli
this->update_P(U, level, lambda * eps);
}
if (level == fl) { // lowest level
this->update_U(U, 0.5 * eps);
} else { // recursive function call
this->step(U, level + 1, first_step, 0);
}
this->update_P(U, level, (1.0 - 2.0 * lambda) * eps);
if (level == fl) { // lowest level
this->update_U(U, 0.5 * eps);
} else { // recursive function call
this->step(U, level + 1, 0, last_step);
}
int mm = (last_step) ? 1 : 2;
this->update_P(U, level, lambda * eps * mm);
@ -430,22 +422,11 @@ class ImplicitMinimumNorm2 : public Integrator<FieldImplementation, SmearingPoli
this->implicit_update_P(U, level, lambda * eps);
}
if (level == fl) { // lowest level
this->implicit_update_U(U, 0.5 * eps,lambda*eps);
} else { // recursive function call
this->step(U, level + 1, first_step, 0);
}
this->implicit_update_U(U, 0.5 * eps,lambda*eps);
this->implicit_update_P(U, level, (1.0 - 2.0 * lambda) * eps, true);
if (level == fl) { // lowest level
this->implicit_update_U(U, 0.5 * eps, (0.5-lambda)*eps);
} else { // recursive function call
this->step(U, level + 1, 0, last_step);
}
//int mm = (last_step) ? 1 : 2;
//this->update_P(U, level, lambda * eps * mm);
this->implicit_update_U(U, 0.5 * eps, (0.5-lambda)*eps);
if (last_step) {
this->update_P2(U, level, eps * lambda);
@ -458,6 +439,96 @@ class ImplicitMinimumNorm2 : public Integrator<FieldImplementation, SmearingPoli
}
};
template <class FieldImplementation, class SmearingPolicy,
class RepresentationPolicy =
Representations<FundamentalRepresentation> >
class ImplicitCampostrini : public Integrator<FieldImplementation, SmearingPolicy,
RepresentationPolicy> {
private:
// const RealD lambda = 0.1931833275037836;
public:
INHERIT_FIELD_TYPES(FieldImplementation);
ImplicitCampostrini(GridBase* grid, IntegratorParameters Par,
ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm, Metric<Field>& M)
: Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
grid, Par, Aset, Sm, M){};
std::string integrator_name(){return "ImplicitCampostrini";}
void step(Field& U, int level, int _first, int _last) {
// level : current level
// fl : final level
// eps : current step size
int fl = this->as.size() - 1;
// assert(Params.lambda.size()>level);
// RealD lambda= Params.lambda[level];
assert(level<3);
RealD lambda= this->Params.lambda0;
if (level>0) lambda= this->Params.lambda1;
if (level>1) lambda= this->Params.lambda2;
std::cout << GridLogMessage << "level: "<<level<< "lambda: "<<lambda<<std::endl;
RealD sigma=pow(2.0,1./3.);
if(level<fl){
//Still Omelyan. Needs to change step() to accept variable stepsize
RealD eps = this->Params.trajL/this->Params.MDsteps * 2.0;
for (int l = 0; l <= level; ++l) eps /= 2.0 * this->as[l].multiplier;
// Nesting: 2xupdate_U of size eps/2
// Next level is eps/2/multiplier
int multiplier = this->as[level].multiplier;
for (int e = 0; e < multiplier; ++e) { // steps per step
int first_step = _first && (e == 0);
int last_step = _last && (e == multiplier - 1);
if (first_step) { // initial half step
this->update_P(U, level, lambda * eps);
}
this->step(U, level + 1, first_step, 0);
this->update_P(U, level, (1.0 - 2.0 * lambda) * eps);
this->step(U, level + 1, 0, last_step);
int mm = (last_step) ? 1 : 2;
this->update_P(U, level, lambda * eps * mm);
}
}
else
{ // last level
RealD dt = this->Params.trajL/this->Params.MDsteps * 2.0;
for (int l = 0; l <= level; ++l) dt /= 2.0 * this->as[l].multiplier;
RealD epsilon = dt/(2.0 - sigma);
int multiplier = this->as[level].multiplier;
for (int e = 0; e < multiplier; ++e) { // steps per step
int first_step = _first && (e == 0);
int last_step = _last && (e == multiplier - 1);
// initial half step
if (first_step) { this->implicit_update_P(U, level, epsilon*0.5); }
this->implicit_update_U(U, epsilon,epsilon*0.5);
this->implicit_update_P(U, level, (1.0 - sigma) * epsilon *0.5, epsilon*0.5, true);
this->implicit_update_U(U, -epsilon*sigma, -epsilon*sigma*0.5);
this->implicit_update_P(U, level, (1.0 - sigma) * epsilon *0.5, -epsilon*sigma*0.5, true);
this->implicit_update_U(U, epsilon,epsilon*0.5);
if (last_step) { this->update_P2(U, level, epsilon*0.5 ); }
else
this->implicit_update_P(U, level, epsilon,epsilon*0.5);
}
}
}
};
NAMESPACE_END(Grid);
#endif // INTEGRATOR_INCLUDED

30
Grid/qcd/utils/CovariantLaplacianRat.h
View File

@ -27,7 +27,7 @@ directory
*************************************************************************************/
/* END LEGAL */
#pragma once
#undef MIXED_CG
#define MIXED_CG
//enable/disable push_back
#undef USE_CHRONO
@ -88,9 +88,13 @@ class LaplacianAdjointRat: public Metric<typename Impl::Field> {
GridBase *grid;
GridBase *grid_f;
CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField> LapStencil;
CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField> LapStencilF;
public:
INHERIT_GIMPL_TYPES(Impl);
// typedef typename GImpl::LinkField GaugeLinkField; \
// typedef typename GImpl::Field GaugeField;
typedef typename ImplF::Field GaugeFieldF;
typedef typename ImplF::LinkField GaugeLinkFieldF; \
GaugeField Usav;
GaugeFieldF UsavF;
std::vector< std::vector<GaugeLinkField> > prev_solnsM;
@ -99,7 +103,7 @@ public:
std::vector< std::vector<GaugeLinkField> > prev_solnsMinvDeriv;
LaplacianAdjointRat(GridBase* _grid, GridBase* _grid_f, OperatorFunction<GaugeField>& S, LaplacianRatParams& gpar, LaplacianRatParams& mpar)
: grid(_grid),grid_f(_grid_f), LapStencil(grid), U(Nd, _grid), Solver(S), Gparam(gpar), Mparam(mpar),Usav(_grid), UsavF(_grid_f),
: grid(_grid),grid_f(_grid_f), LapStencil(_grid), LapStencilF(_grid_f), U(Nd, _grid), Solver(S), Gparam(gpar), Mparam(mpar),Usav(_grid), UsavF(_grid_f),
prev_solnsM(4),prev_solnsMinv(4),prev_solnsMDeriv(4),prev_solnsMinvDeriv(4) {
// std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
this->triv=0;
@ -147,6 +151,7 @@ public:
RealD fac = - 1. / (double(4 * Nd)) ;
RealD coef=0.5;
LapStencil.GaugeImport(Usav);
LapStencilF.GaugeImport(UsavF);
for (int nu=0;nu<Nd;nu++){
@ -189,8 +194,9 @@ for (int nu=0;nu<Nd;nu++){
#ifndef MIXED_CG
CG(QuadOp,right_nu,MinvMom[i]);
#else
QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeLinkFieldF> QuadOpF(LaplacianF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeField,GaugeFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
// QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
MixedCG.InnerTolerance=par.tolerance;
MixedCG(right_nu,MinvMom[i]);
#endif
@ -211,11 +217,13 @@ for (int nu=0;nu<Nd;nu++){
LapStencil.M(MinvGMom, Gtemp2); LMinvGMom=fac*Gtemp2;
CG(QuadOp,right_nu,MinvMom[i]);
#else
QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeLinkFieldF> QuadOpF(LaplacianF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeField,GaugeFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
// QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
MixedCG.InnerTolerance=par.tolerance;
MixedCG(GMom,MinvGMom);
Laplacian.M(MinvGMom, LMinvGMom);
LapStencil.M(MinvGMom, Gtemp2); LMinvGMom=fac*Gtemp2;
// Laplacian.M(MinvGMom, LMinvGMom);
MixedCG(right_nu,MinvMom[i]);
#endif
#if USE_CHRONO
@ -277,6 +285,7 @@ for (int nu=0;nu<Nd;nu++){
RealD fac = -1. / (double(4 * Nd));
LapStencil.GaugeImport(Usav);
LapStencilF.GaugeImport(UsavF);
for(int nu=0; nu<Nd;nu++){
GaugeLinkField P_nu = PeekIndex<LorentzIndex>(P, nu);
GaugeLinkField Gp(P.Grid());
@ -297,10 +306,11 @@ for(int nu=0; nu<Nd;nu++){
#ifndef MIXED_CG
CG(QuadOp,P_nu,Gtemp);
#else
QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeFieldF> QuadOpF(LaplacianF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeField,GaugeFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
// QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
MixedCG.InnerTolerance=par.tolerance;
MixedCG(P,Gtemp[i]);
MixedCG(P_nu,Gtemp);
#endif
#if USE_CHRONO
prev_solns[nu].push_back(Gtemp);

48
HMC/Mobius2p1p1fEOFA_4Gev.cc
View File

@ -37,7 +37,7 @@ directory
// second level EOFA
#undef EOFA_H
#define USE_OBC
#undef DO_IMPLICIT
#define DO_IMPLICIT
NAMESPACE_BEGIN(Grid);
@ -269,11 +269,11 @@ int main(int argc, char **argv) {
Coordinate mpi = GridDefaultMpi();
Coordinate simdF = GridDefaultSimd(Nd,vComplexF::Nsimd());
Coordinate simdD = GridDefaultSimd(Nd,vComplexD::Nsimd());
auto GridPtrF = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi);
//auto UGrid_f = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi);
auto GridRBPtrF = SpaceTimeGrid::makeFourDimRedBlackGrid(GridPtrF);
auto FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtrF);
auto FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtrF);
// auto GridPtrF = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi);
auto UGrid_f = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi);
auto GridRBPtrF = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_f);
auto FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_f);
auto FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_f);
#ifndef USE_OBC
@ -287,7 +287,7 @@ int main(int argc, char **argv) {
// temporarily need a gauge field
LatticeGaugeField U(GridPtr);
LatticeGaugeFieldF UF(GridPtrF);
LatticeGaugeFieldF UF(UGrid_f);
// These lines are unecessary if BC are all periodic
#ifndef USE_OBC
@ -330,15 +330,15 @@ int main(int argc, char **argv) {
MobiusEOFAFermionD Strange_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , strange_mass, strange_mass, charm_mass, 0.0, -1, M5, b, c);
MobiusEOFAFermionF Strange_Op_LF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, strange_mass, strange_mass, charm_mass, 0.0, -1, M5, b, c);
MobiusEOFAFermionF Strange_Op_LF(UF, *FGridF, *FrbGridF, *UGrid_f, *GridRBPtrF, strange_mass, strange_mass, charm_mass, 0.0, -1, M5, b, c);
MobiusEOFAFermionD Strange_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , charm_mass, strange_mass, charm_mass, -1.0, 1, M5, b, c);
MobiusEOFAFermionF Strange_Op_RF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, charm_mass, strange_mass, charm_mass, -1.0, 1, M5, b, c);
MobiusEOFAFermionF Strange_Op_RF(UF, *FGridF, *FrbGridF, *UGrid_f, *GridRBPtrF, charm_mass, strange_mass, charm_mass, -1.0, 1, M5, b, c);
#ifdef EOFA_H
MobiusEOFAFermionD Strange2_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , eofa_mass, eofa_mass, charm_mass , 0.0, -1, M5, b, c);
MobiusEOFAFermionF Strange2_Op_LF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, eofa_mass, eofa_mass, charm_mass , 0.0, -1, M5, b, c);
MobiusEOFAFermionF Strange2_Op_LF(UF, *FGridF, *FrbGridF, *UGrid_f, *GridRBPtrF, eofa_mass, eofa_mass, charm_mass , 0.0, -1, M5, b, c);
MobiusEOFAFermionD Strange2_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , charm_mass , eofa_mass, charm_mass , -1.0, 1, M5, b, c);
MobiusEOFAFermionF Strange2_Op_RF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, charm_mass , eofa_mass, charm_mass , -1.0, 1, M5, b, c);
MobiusEOFAFermionF Strange2_Op_RF(UF, *FGridF, *FrbGridF, *UGrid_f, *GridRBPtrF, charm_mass , eofa_mass, charm_mass , -1.0, 1, M5, b, c);
#endif
ConjugateGradient<FermionField> ActionCG(ActionStoppingCondition,MaxCGIterations);
@ -363,7 +363,7 @@ int main(int argc, char **argv) {
MxPCG_EOFA ActionCGL(ActionStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
UGrid_f,
FrbGridF,
Strange_Op_LF,Strange_Op_L,
Strange_LinOp_LF,Strange_LinOp_L);
@ -372,7 +372,7 @@ int main(int argc, char **argv) {
MxPCG_EOFA ActionCGL2(ActionStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
UGrid_f,
FrbGridF,
Strange2_Op_LF,Strange2_Op_L,
Strange2_LinOp_LF,Strange2_LinOp_L);
@ -381,7 +381,7 @@ int main(int argc, char **argv) {
MxPCG_EOFA DerivativeCGL(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
UGrid_f,
FrbGridF,
Strange_Op_LF,Strange_Op_L,
Strange_LinOp_LF,Strange_LinOp_L);
@ -390,7 +390,7 @@ int main(int argc, char **argv) {
MxPCG_EOFA DerivativeCGL2(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
UGrid_f,
FrbGridF,
Strange2_Op_LF,Strange2_Op_L,
Strange2_LinOp_LF,Strange2_LinOp_L);
@ -399,7 +399,7 @@ int main(int argc, char **argv) {
MxPCG_EOFA ActionCGR(ActionStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
UGrid_f,
FrbGridF,
Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
@ -408,7 +408,7 @@ int main(int argc, char **argv) {
MxPCG_EOFA ActionCGR2(ActionStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
UGrid_f,
FrbGridF,
Strange2_Op_RF,Strange2_Op_R,
Strange2_LinOp_RF,Strange2_LinOp_R);
@ -417,7 +417,7 @@ int main(int argc, char **argv) {
MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
UGrid_f,
FrbGridF,
Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
@ -426,7 +426,7 @@ int main(int argc, char **argv) {
MxPCG_EOFA DerivativeCGR2(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
UGrid_f,
FrbGridF,
Strange2_Op_RF,Strange2_Op_R,
Strange2_LinOp_RF,Strange2_LinOp_R);
@ -515,7 +515,7 @@ int main(int argc, char **argv) {
////////////////////////////////////////////////////////////////////////////
double DerivativeStoppingConditionLoose = 1e-8;
DenominatorsF.push_back(new FermionActionF(UF,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,light_den[h],M5,b,c, ParamsF));
DenominatorsF.push_back(new FermionActionF(UF,*FGridF,*FrbGridF,*UGrid_f,*GridRBPtrF,light_den[h],M5,b,c, ParamsF));
LinOpD.push_back(new LinearOperatorD(*Denominators[h]));
LinOpF.push_back(new LinearOperatorF(*DenominatorsF[h]));
@ -524,7 +524,7 @@ int main(int argc, char **argv) {
MPCG.push_back(new MxPCG(conv,
MX_inner,
MaxCGIterations,
GridPtrF,
UGrid_f,
FrbGridF,
*DenominatorsF[h],*Denominators[h],
*LinOpF[h], *LinOpD[h]) );
@ -532,7 +532,7 @@ int main(int argc, char **argv) {
ActionMPCG.push_back(new MxPCG(ActionStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
UGrid_f,
FrbGridF,
*DenominatorsF[h],*Denominators[h],
*LinOpF[h], *LinOpD[h]) );
@ -619,12 +619,12 @@ int main(int argc, char **argv) {
std::cout << GridLogMessage << " b2= " << mpar.b2 <<std::endl;
// Assumes PeriodicGimplR or D at the moment
auto UGrid = TheHMC.Resources.GetCartesian("gauge");
// auto UGrid_f = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi);
// auto UGrid_f = GridPtrF;
// auto GridPtrF = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi);
// std::cout << GridLogMessage << " UGrid= " << UGrid <<std::endl;
// std::cout << GridLogMessage << " UGrid_f= " << UGrid_f <<std::endl;
LaplacianAdjointRat<HMCWrapper::ImplPolicy, PeriodicGimplF> Mtr(UGrid, GridPtrF ,CG, gpar, mpar);
LaplacianAdjointRat<HMCWrapper::ImplPolicy, PeriodicGimplF> Mtr(UGrid, UGrid_f ,CG, gpar, mpar);
#endif
std::cout << GridLogMessage << " Running the HMC "<< std::endl;