portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-04 19:25:56 +01:00
Code Releases Activity

Merge branch 'fix/removeQCDremnants' into fix/remQCDns_ignore_ws

Browse Source
This commit is contained in:
Nils Asmussen 2019-10-02 16:11:27 +01:00
commit d5a180d914
17 changed files with 534 additions and 547 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Grid/parallelIO/NerscIO.h
View File

@ -354,6 +354,6 @@ public:
} }
}; };
NAMESPACE_END(QCD); NAMESPACE_END(Grid);
#endif #endif

1
Grid/qcd/hmc/HMC_GridModules.h
View File

@ -97,7 +97,6 @@ protected:
//////////////////////////////////// ////////////////////////////////////
// Classes for the user // Classes for the user
//////////////////////////////////// ////////////////////////////////////
// Note: the space time grid should be out of the QCD namespace
template <class vector_type> template <class vector_type>
class GridFourDimModule : public GridModule class GridFourDimModule : public GridModule
{ {

6
Grid/qcd/utils/CovariantSmearing.h
View File

@ -27,8 +27,7 @@ directory
*************************************************************************************/ *************************************************************************************/
#pragma once #pragma once
namespace Grid { NAMESPACE_BEGIN(Grid);
namespace QCD {
template <class Gimpl> class CovariantSmearing : public Gimpl template <class Gimpl> class CovariantSmearing : public Gimpl
{ {
@ -84,4 +83,5 @@ public:
} }
} }
}; };
}}
NAMESPACE_END(Grid);

1
Grid/qcd/utils/LinalgUtils.h
View File

@ -201,7 +201,6 @@ void G5R5(Lattice<vobj> &z,const Lattice<vobj> &x)
}); });
} }
// I explicitly need these outside the QCD namespace
template<typename vobj> template<typename vobj>
void G5C(Lattice<vobj> &z, const Lattice<vobj> &x) void G5C(Lattice<vobj> &z, const Lattice<vobj> &x)
{ {

17
HMC/Mobius2p1f.cc
View File

@ -31,7 +31,6 @@ directory
int main(int argc, char **argv) { int main(int argc, char **argv) {
using namespace Grid; using namespace Grid;
using namespace Grid::QCD;
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
int threads = GridThread::GetThreads(); int threads = GridThread::GetThreads();
@ -44,18 +43,18 @@ int main(int argc, char **argv) {
typedef typename FermionAction::FermionField FermionField; typedef typename FermionAction::FermionField FermionField;
typedef Grid::XmlReader Serialiser; typedef Grid::XmlReader Serialiser;
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
IntegratorParameters MD; IntegratorParameters MD;
// typedef GenericHMCRunner<LeapFrog> HMCWrapper; // typedef GenericHMCRunner<LeapFrog> HMCWrapper;
// MD.name = std::string("Leap Frog"); // MD.name = std::string("Leap Frog");
// typedef GenericHMCRunner<ForceGradient> HMCWrapper; // typedef GenericHMCRunner<ForceGradient> HMCWrapper;
// MD.name = std::string("Force Gradient"); // MD.name = std::string("Force Gradient");
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
MD.name = std::string("MinimumNorm2"); MD.name = std::string("MinimumNorm2");
MD.MDsteps = 20; MD.MDsteps = 20;
MD.trajL = 1.0; MD.trajL = 1.0;
HMCparameters HMCparams; HMCparameters HMCparams;
HMCparams.StartTrajectory = 0; HMCparams.StartTrajectory = 0;
HMCparams.Trajectories = 200; HMCparams.Trajectories = 200;
@ -67,7 +66,7 @@ int main(int argc, char **argv) {
// Grid from the command line arguments --grid and --mpi // Grid from the command line arguments --grid and --mpi
TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
CheckpointerParameters CPparams; CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_EODWF_lat"; CPparams.config_prefix = "ckpoint_EODWF_lat";
CPparams.rng_prefix = "ckpoint_EODWF_rng"; CPparams.rng_prefix = "ckpoint_EODWF_rng";
@ -81,7 +80,7 @@ int main(int argc, char **argv) {
TheHMC.Resources.SetRNGSeeds(RNGpar); TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables // Construct observables
// here there is too much indirection // here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs; typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>(); TheHMC.Resources.AddObservable<PlaqObs>();
////////////////////////////////////////////// //////////////////////////////////////////////
@ -118,7 +117,7 @@ int main(int argc, char **argv) {
// These lines are unecessary if BC are all periodic // These lines are unecessary if BC are all periodic
std::vector<Complex> boundary = {1,1,1,-1}; std::vector<Complex> boundary = {1,1,1,-1};
FermionAction::ImplParams Params(boundary); FermionAction::ImplParams Params(boundary);
double StoppingCondition = 1e-10; double StoppingCondition = 1e-10;
double MaxCGIterations = 30000; double MaxCGIterations = 30000;
ConjugateGradient<FermionField> CG(StoppingCondition,MaxCGIterations); ConjugateGradient<FermionField> CG(StoppingCondition,MaxCGIterations);

373
HMC/Mobius2p1fEOFA.cc
View File

@ -2,7 +2,7 @@
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: Source file:
Copyright (C) 2015-2016 Copyright (C) 2015-2016
@ -34,140 +34,139 @@ directory
#define MIXED_PRECISION #define MIXED_PRECISION
#endif #endif
namespace Grid{ NAMESPACE_BEGIN(Grid);
namespace QCD{
/* /*
* Need a plan for gauge field update for mixed precision in HMC (2x speed up) * Need a plan for gauge field update for mixed precision in HMC (2x speed up)
* -- Store the single prec action operator. * -- Store the single prec action operator.
* -- Clone the gauge field from the operator function argument. * -- Clone the gauge field from the operator function argument.
* -- Build the mixed precision operator dynamically from the passed operator and single prec clone. * -- Build the mixed precision operator dynamically from the passed operator and single prec clone.
*/ */
template<class FermionOperatorD, class FermionOperatorF, class SchurOperatorD, class SchurOperatorF> template<class FermionOperatorD, class FermionOperatorF, class SchurOperatorD, class SchurOperatorF>
class MixedPrecisionConjugateGradientOperatorFunction : public OperatorFunction<typename FermionOperatorD::FermionField> { class MixedPrecisionConjugateGradientOperatorFunction : public OperatorFunction<typename FermionOperatorD::FermionField> {
public: public:
typedef typename FermionOperatorD::FermionField FieldD; typedef typename FermionOperatorD::FermionField FieldD;
typedef typename FermionOperatorF::FermionField FieldF; typedef typename FermionOperatorF::FermionField FieldF;
using OperatorFunction<FieldD>::operator(); using OperatorFunction<FieldD>::operator();
RealD Tolerance; RealD Tolerance;
RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed
Integer MaxInnerIterations; Integer MaxInnerIterations;
Integer MaxOuterIterations; Integer MaxOuterIterations;
GridBase* SinglePrecGrid4; //Grid for single-precision fields GridBase* SinglePrecGrid4; //Grid for single-precision fields
GridBase* SinglePrecGrid5; //Grid for single-precision fields GridBase* SinglePrecGrid5; //Grid for single-precision fields
RealD OuterLoopNormMult; //Stop the outer loop and move to a final double prec solve when the residual is OuterLoopNormMult * Tolerance RealD OuterLoopNormMult; //Stop the outer loop and move to a final double prec solve when the residual is OuterLoopNormMult * Tolerance
FermionOperatorF &FermOpF; FermionOperatorF &FermOpF;
FermionOperatorD &FermOpD;; FermionOperatorD &FermOpD;;
SchurOperatorF &LinOpF; SchurOperatorF &LinOpF;
SchurOperatorD &LinOpD; SchurOperatorD &LinOpD;
Integer TotalInnerIterations; //Number of inner CG iterations Integer TotalInnerIterations; //Number of inner CG iterations
Integer TotalOuterIterations; //Number of restarts Integer TotalOuterIterations; //Number of restarts
Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step
MixedPrecisionConjugateGradientOperatorFunction(RealD tol, MixedPrecisionConjugateGradientOperatorFunction(RealD tol,
Integer maxinnerit, Integer maxinnerit,
Integer maxouterit, Integer maxouterit,
GridBase* _sp_grid4, GridBase* _sp_grid4,
GridBase* _sp_grid5, GridBase* _sp_grid5,
FermionOperatorF &_FermOpF, FermionOperatorF &_FermOpF,
FermionOperatorD &_FermOpD, FermionOperatorD &_FermOpD,
SchurOperatorF &_LinOpF, SchurOperatorF &_LinOpF,
SchurOperatorD &_LinOpD): SchurOperatorD &_LinOpD):
LinOpF(_LinOpF), LinOpF(_LinOpF),
LinOpD(_LinOpD), LinOpD(_LinOpD),
FermOpF(_FermOpF), FermOpF(_FermOpF),
FermOpD(_FermOpD), FermOpD(_FermOpD),
Tolerance(tol), Tolerance(tol),
InnerTolerance(tol), InnerTolerance(tol),
MaxInnerIterations(maxinnerit), MaxInnerIterations(maxinnerit),
MaxOuterIterations(maxouterit), MaxOuterIterations(maxouterit),
SinglePrecGrid4(_sp_grid4), SinglePrecGrid4(_sp_grid4),
SinglePrecGrid5(_sp_grid5), SinglePrecGrid5(_sp_grid5),
OuterLoopNormMult(100.) OuterLoopNormMult(100.)
{ {
/* Debugging instances of objects; references are stored /* Debugging instances of objects; references are stored
std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpF " <<std::hex<< &LinOpF<<std::dec <<std::endl; std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpF " <<std::hex<< &LinOpF<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpD " <<std::hex<< &LinOpD<<std::dec <<std::endl; std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpD " <<std::hex<< &LinOpD<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpF " <<std::hex<< &FermOpF<<std::dec <<std::endl; std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpF " <<std::hex<< &FermOpF<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpD " <<std::hex<< &FermOpD<<std::dec <<std::endl; std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpD " <<std::hex<< &FermOpD<<std::dec <<std::endl;
*/ */
};
void operator()(LinearOperatorBase<FieldD> &LinOpU, const FieldD &src, FieldD &psi) {
std::cout << GridLogMessage << " Mixed precision CG wrapper operator() "<<std::endl;
SchurOperatorD * SchurOpU = static_cast<SchurOperatorD *>(&LinOpU);
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpU " <<std::hex<< &(SchurOpU->_Mat)<<std::dec <<std::endl;
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpD " <<std::hex<< &(LinOpD._Mat) <<std::dec <<std::endl;
// Assumption made in code to extract gauge field
// We could avoid storing LinopD reference alltogether ?
assert(&(SchurOpU->_Mat)==&(LinOpD._Mat));
////////////////////////////////////////////////////////////////////////////////////
// Must snarf a single precision copy of the gauge field in Linop_d argument
////////////////////////////////////////////////////////////////////////////////////
typedef typename FermionOperatorF::GaugeField GaugeFieldF;
typedef typename FermionOperatorF::GaugeLinkField GaugeLinkFieldF;
typedef typename FermionOperatorD::GaugeField GaugeFieldD;
typedef typename FermionOperatorD::GaugeLinkField GaugeLinkFieldD;
GridBase * GridPtrF = SinglePrecGrid4;
GridBase * GridPtrD = FermOpD.Umu.Grid();
GaugeFieldF U_f (GridPtrF);
GaugeLinkFieldF Umu_f(GridPtrF);
// std::cout << " Dim gauge field "<<GridPtrF->Nd()<<std::endl; // 4d
// std::cout << " Dim gauge field "<<GridPtrD->Nd()<<std::endl; // 4d
////////////////////////////////////////////////////////////////////////////////////
// Moving this to a Clone method of fermion operator would allow to duplicate the
// physics parameters and decrease gauge field copies
////////////////////////////////////////////////////////////////////////////////////
GaugeLinkFieldD Umu_d(GridPtrD);
for(int mu=0;mu<Nd*2;mu++){
Umu_d = PeekIndex<LorentzIndex>(FermOpD.Umu, mu);
precisionChange(Umu_f,Umu_d);
PokeIndex<LorentzIndex>(FermOpF.Umu, Umu_f, mu);
}
pickCheckerboard(Even,FermOpF.UmuEven,FermOpF.Umu);
pickCheckerboard(Odd ,FermOpF.UmuOdd ,FermOpF.Umu);
////////////////////////////////////////////////////////////////////////////////////
// Could test to make sure that LinOpF and LinOpD agree to single prec?
////////////////////////////////////////////////////////////////////////////////////
/*
GridBase *Fgrid = psi._grid;
FieldD tmp2(Fgrid);
FieldD tmp1(Fgrid);
LinOpU.Op(src,tmp1);
LinOpD.Op(src,tmp2);
std::cout << " Double gauge field "<< norm2(FermOpD.Umu)<<std::endl;
std::cout << " Single gauge field "<< norm2(FermOpF.Umu)<<std::endl;
std::cout << " Test of operators "<<norm2(tmp1)<<std::endl;
std::cout << " Test of operators "<<norm2(tmp2)<<std::endl;
tmp1=tmp1-tmp2;
std::cout << " Test of operators diff "<<norm2(tmp1)<<std::endl;
*/
////////////////////////////////////////////////////////////////////////////////////
// Make a mixed precision conjugate gradient
////////////////////////////////////////////////////////////////////////////////////
MixedPrecisionConjugateGradient<FieldD,FieldF> MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD);
std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" <<std::endl;
MPCG(src,psi);
}
}; };
}};
void operator()(LinearOperatorBase<FieldD> &LinOpU, const FieldD &src, FieldD &psi) {
std::cout << GridLogMessage << " Mixed precision CG wrapper operator() "<<std::endl;
SchurOperatorD * SchurOpU = static_cast<SchurOperatorD *>(&LinOpU);
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpU " <<std::hex<< &(SchurOpU->_Mat)<<std::dec <<std::endl;
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpD " <<std::hex<< &(LinOpD._Mat) <<std::dec <<std::endl;
// Assumption made in code to extract gauge field
// We could avoid storing LinopD reference alltogether ?
assert(&(SchurOpU->_Mat)==&(LinOpD._Mat));
////////////////////////////////////////////////////////////////////////////////////
// Must snarf a single precision copy of the gauge field in Linop_d argument
////////////////////////////////////////////////////////////////////////////////////
typedef typename FermionOperatorF::GaugeField GaugeFieldF;
typedef typename FermionOperatorF::GaugeLinkField GaugeLinkFieldF;
typedef typename FermionOperatorD::GaugeField GaugeFieldD;
typedef typename FermionOperatorD::GaugeLinkField GaugeLinkFieldD;
GridBase * GridPtrF = SinglePrecGrid4;
GridBase * GridPtrD = FermOpD.Umu.Grid();
GaugeFieldF U_f (GridPtrF);
GaugeLinkFieldF Umu_f(GridPtrF);
// std::cout << " Dim gauge field "<<GridPtrF->Nd()<<std::endl; // 4d
// std::cout << " Dim gauge field "<<GridPtrD->Nd()<<std::endl; // 4d
////////////////////////////////////////////////////////////////////////////////////
// Moving this to a Clone method of fermion operator would allow to duplicate the
// physics parameters and decrease gauge field copies
////////////////////////////////////////////////////////////////////////////////////
GaugeLinkFieldD Umu_d(GridPtrD);
for(int mu=0;mu<Nd*2;mu++){
Umu_d = PeekIndex<LorentzIndex>(FermOpD.Umu, mu);
precisionChange(Umu_f,Umu_d);
PokeIndex<LorentzIndex>(FermOpF.Umu, Umu_f, mu);
}
pickCheckerboard(Even,FermOpF.UmuEven,FermOpF.Umu);
pickCheckerboard(Odd ,FermOpF.UmuOdd ,FermOpF.Umu);
////////////////////////////////////////////////////////////////////////////////////
// Could test to make sure that LinOpF and LinOpD agree to single prec?
////////////////////////////////////////////////////////////////////////////////////
/*
GridBase *Fgrid = psi._grid;
FieldD tmp2(Fgrid);
FieldD tmp1(Fgrid);
LinOpU.Op(src,tmp1);
LinOpD.Op(src,tmp2);
std::cout << " Double gauge field "<< norm2(FermOpD.Umu)<<std::endl;
std::cout << " Single gauge field "<< norm2(FermOpF.Umu)<<std::endl;
std::cout << " Test of operators "<<norm2(tmp1)<<std::endl;
std::cout << " Test of operators "<<norm2(tmp2)<<std::endl;
tmp1=tmp1-tmp2;
std::cout << " Test of operators diff "<<norm2(tmp1)<<std::endl;
*/
////////////////////////////////////////////////////////////////////////////////////
// Make a mixed precision conjugate gradient
////////////////////////////////////////////////////////////////////////////////////
MixedPrecisionConjugateGradient<FieldD,FieldF> MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD);
std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" <<std::endl;
MPCG(src,psi);
}
};
NAMESPACE_END(Grid);
int main(int argc, char **argv) { int main(int argc, char **argv) {
using namespace Grid; using namespace Grid;
using namespace Grid::QCD;
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
int threads = GridThread::GetThreads(); int threads = GridThread::GetThreads();
@ -184,18 +183,18 @@ int main(int argc, char **argv) {
typedef typename FermionActionF::FermionField FermionFieldF; typedef typename FermionActionF::FermionField FermionFieldF;
typedef Grid::XmlReader Serialiser; typedef Grid::XmlReader Serialiser;
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
IntegratorParameters MD; IntegratorParameters MD;
// typedef GenericHMCRunner<LeapFrog> HMCWrapper; // typedef GenericHMCRunner<LeapFrog> HMCWrapper;
// MD.name = std::string("Leap Frog"); // MD.name = std::string("Leap Frog");
typedef GenericHMCRunner<ForceGradient> HMCWrapper; typedef GenericHMCRunner<ForceGradient> HMCWrapper;
MD.name = std::string("Force Gradient"); MD.name = std::string("Force Gradient");
// typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; // typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
// MD.name = std::string("MinimumNorm2"); // MD.name = std::string("MinimumNorm2");
MD.MDsteps = 6; MD.MDsteps = 6;
MD.trajL = 1.0; MD.trajL = 1.0;
HMCparameters HMCparams; HMCparameters HMCparams;
HMCparams.StartTrajectory = 590; HMCparams.StartTrajectory = 590;
HMCparams.Trajectories = 1000; HMCparams.Trajectories = 1000;
@ -208,7 +207,7 @@ int main(int argc, char **argv) {
// Grid from the command line arguments --grid and --mpi // Grid from the command line arguments --grid and --mpi
TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
CheckpointerParameters CPparams; CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_EODWF_lat"; CPparams.config_prefix = "ckpoint_EODWF_lat";
CPparams.rng_prefix = "ckpoint_EODWF_rng"; CPparams.rng_prefix = "ckpoint_EODWF_rng";
@ -222,7 +221,7 @@ int main(int argc, char **argv) {
TheHMC.Resources.SetRNGSeeds(RNGpar); TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables // Construct observables
// here there is too much indirection // here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs; typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>(); TheHMC.Resources.AddObservable<PlaqObs>();
////////////////////////////////////////////// //////////////////////////////////////////////
@ -233,7 +232,7 @@ int main(int argc, char **argv) {
Real strange_mass = 0.04; Real strange_mass = 0.04;
Real pv_mass = 1.0; Real pv_mass = 1.0;
RealD M5 = 1.8; RealD M5 = 1.8;
RealD b = 1.0; RealD b = 1.0;
RealD c = 0.0; RealD c = 0.0;
std::vector<Real> hasenbusch({ 0.1, 0.3, 0.6 }); std::vector<Real> hasenbusch({ 0.1, 0.3, 0.6 });
@ -262,7 +261,7 @@ int main(int argc, char **argv) {
std::vector<Complex> boundary = {1,1,1,-1}; std::vector<Complex> boundary = {1,1,1,-1};
FermionAction::ImplParams Params(boundary); FermionAction::ImplParams Params(boundary);
FermionActionF::ImplParams ParamsF(boundary); FermionActionF::ImplParams ParamsF(boundary);
double ActionStoppingCondition = 1e-10; double ActionStoppingCondition = 1e-10;
double DerivativeStoppingCondition = 1e-6; double DerivativeStoppingCondition = 1e-6;
double MaxCGIterations = 30000; double MaxCGIterations = 30000;
@ -293,7 +292,7 @@ int main(int argc, char **argv) {
OFRp.degree = 14; OFRp.degree = 14;
OFRp.precision= 50; OFRp.precision= 50;
MobiusEOFAFermionR Strange_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c); MobiusEOFAFermionR Strange_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c);
MobiusEOFAFermionF Strange_Op_LF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c); MobiusEOFAFermionF Strange_Op_LF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c);
MobiusEOFAFermionR Strange_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c); MobiusEOFAFermionR Strange_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c);
@ -310,50 +309,50 @@ int main(int argc, char **argv) {
LinearOperatorEOFAF Strange_LinOp_RF(Strange_Op_RF); LinearOperatorEOFAF Strange_LinOp_RF(Strange_Op_RF);
MxPCG_EOFA ActionCGL(ActionStoppingCondition, MxPCG_EOFA ActionCGL(ActionStoppingCondition,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
Strange_Op_LF,Strange_Op_L, Strange_Op_LF,Strange_Op_L,
Strange_LinOp_LF,Strange_LinOp_L); Strange_LinOp_LF,Strange_LinOp_L);
MxPCG_EOFA DerivativeCGL(DerivativeStoppingCondition, MxPCG_EOFA DerivativeCGL(DerivativeStoppingCondition,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
Strange_Op_LF,Strange_Op_L, Strange_Op_LF,Strange_Op_L,
Strange_LinOp_LF,Strange_LinOp_L); Strange_LinOp_LF,Strange_LinOp_L);
MxPCG_EOFA ActionCGR(ActionStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> MxPCG_EOFA ActionCGR(ActionStoppingCondition,
EOFA(Strange_Op_L, Strange_Op_R, MX_inner,
ActionCG, MaxCGIterations,
ActionCGL, ActionCGR, GridPtrF,
DerivativeCGL, DerivativeCGR, FrbGridF,
OFRp, true); Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy>
EOFA(Strange_Op_L, Strange_Op_R,
ActionCG,
ActionCGL, ActionCGR,
DerivativeCGL, DerivativeCGR,
OFRp, true);
#else #else
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy>
EOFA(Strange_Op_L, Strange_Op_R, EOFA(Strange_Op_L, Strange_Op_R,
ActionCG, ActionCG,
ActionCG, ActionCG, ActionCG, ActionCG,
DerivativeCG, DerivativeCG, DerivativeCG, DerivativeCG,
OFRp, true); OFRp, true);
#endif #endif
Level1.push_back(&EOFA); Level1.push_back(&EOFA);
@ -384,7 +383,7 @@ int main(int argc, char **argv) {
std::vector<MxPCG *> MPCG; std::vector<MxPCG *> MPCG;
std::vector<FermionActionF *> DenominatorsF; std::vector<FermionActionF *> DenominatorsF;
std::vector<LinearOperatorD *> LinOpD; std::vector<LinearOperatorD *> LinOpD;
std::vector<LinearOperatorF *> LinOpF; std::vector<LinearOperatorF *> LinOpF;
for(int h=0;h<n_hasenbusch+1;h++){ for(int h=0;h<n_hasenbusch+1;h++){
@ -403,20 +402,20 @@ int main(int argc, char **argv) {
LinOpF.push_back(new LinearOperatorF(*DenominatorsF[h])); LinOpF.push_back(new LinearOperatorF(*DenominatorsF[h]));
MPCG.push_back(new MxPCG(DerivativeStoppingCondition, MPCG.push_back(new MxPCG(DerivativeStoppingCondition,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
*DenominatorsF[h],*Denominators[h], *DenominatorsF[h],*Denominators[h],
*LinOpF[h], *LinOpD[h]) ); *LinOpF[h], *LinOpD[h]) );
ActionMPCG.push_back(new MxPCG(ActionStoppingCondition, ActionMPCG.push_back(new MxPCG(ActionStoppingCondition,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
*DenominatorsF[h],*Denominators[h], *DenominatorsF[h],*Denominators[h],
*LinOpF[h], *LinOpD[h]) ); *LinOpF[h], *LinOpD[h]) );
// Heatbath not mixed yet. As inverts numerators not so important as raised mass. // Heatbath not mixed yet. As inverts numerators not so important as raised mass.
Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],ActionCG)); Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],ActionCG));

342
HMC/Mobius2p1fEOFA_F1.cc
View File

@ -2,7 +2,7 @@
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: Source file:
Copyright (C) 2015-2016 Copyright (C) 2015-2016
@ -34,123 +34,123 @@ directory
#define MIXED_PRECISION #define MIXED_PRECISION
#endif #endif
namespace Grid{ NAMESPACE_BEGIN(Grid);
namespace QCD{
/* /*
* Need a plan for gauge field update for mixed precision in HMC (2x speed up) * Need a plan for gauge field update for mixed precision in HMC (2x speed up)
* -- Store the single prec action operator. * -- Store the single prec action operator.
* -- Clone the gauge field from the operator function argument. * -- Clone the gauge field from the operator function argument.
* -- Build the mixed precision operator dynamically from the passed operator and single prec clone. * -- Build the mixed precision operator dynamically from the passed operator and single prec clone.
*/ */
template<class FermionOperatorD, class FermionOperatorF, class SchurOperatorD, class SchurOperatorF> template<class FermionOperatorD, class FermionOperatorF, class SchurOperatorD, class SchurOperatorF>
class MixedPrecisionConjugateGradientOperatorFunction : public OperatorFunction<typename FermionOperatorD::FermionField> { class MixedPrecisionConjugateGradientOperatorFunction : public OperatorFunction<typename FermionOperatorD::FermionField> {
public: public:
typedef typename FermionOperatorD::FermionField FieldD; typedef typename FermionOperatorD::FermionField FieldD;
typedef typename FermionOperatorF::FermionField FieldF; typedef typename FermionOperatorF::FermionField FieldF;
using OperatorFunction<FieldD>::operator(); using OperatorFunction<FieldD>::operator();
RealD Tolerance; RealD Tolerance;
RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed
Integer MaxInnerIterations; Integer MaxInnerIterations;
Integer MaxOuterIterations; Integer MaxOuterIterations;
GridBase* SinglePrecGrid4; //Grid for single-precision fields GridBase* SinglePrecGrid4; //Grid for single-precision fields
GridBase* SinglePrecGrid5; //Grid for single-precision fields GridBase* SinglePrecGrid5; //Grid for single-precision fields
RealD OuterLoopNormMult; //Stop the outer loop and move to a final double prec solve when the residual is OuterLoopNormMult * Tolerance RealD OuterLoopNormMult; //Stop the outer loop and move to a final double prec solve when the residual is OuterLoopNormMult * Tolerance
FermionOperatorF &FermOpF; FermionOperatorF &FermOpF;
FermionOperatorD &FermOpD;; FermionOperatorD &FermOpD;;
SchurOperatorF &LinOpF; SchurOperatorF &LinOpF;
SchurOperatorD &LinOpD; SchurOperatorD &LinOpD;
Integer TotalInnerIterations; //Number of inner CG iterations Integer TotalInnerIterations; //Number of inner CG iterations
Integer TotalOuterIterations; //Number of restarts Integer TotalOuterIterations; //Number of restarts
Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step
MixedPrecisionConjugateGradientOperatorFunction(RealD tol, MixedPrecisionConjugateGradientOperatorFunction(RealD tol,
Integer maxinnerit, Integer maxinnerit,
Integer maxouterit, Integer maxouterit,
GridBase* _sp_grid4, GridBase* _sp_grid4,
GridBase* _sp_grid5, GridBase* _sp_grid5,
FermionOperatorF &_FermOpF, FermionOperatorF &_FermOpF,
FermionOperatorD &_FermOpD, FermionOperatorD &_FermOpD,
SchurOperatorF &_LinOpF, SchurOperatorF &_LinOpF,
SchurOperatorD &_LinOpD): SchurOperatorD &_LinOpD):
LinOpF(_LinOpF), LinOpF(_LinOpF),
LinOpD(_LinOpD), LinOpD(_LinOpD),
FermOpF(_FermOpF), FermOpF(_FermOpF),
FermOpD(_FermOpD), FermOpD(_FermOpD),
Tolerance(tol), Tolerance(tol),
InnerTolerance(tol), InnerTolerance(tol),
MaxInnerIterations(maxinnerit), MaxInnerIterations(maxinnerit),
MaxOuterIterations(maxouterit), MaxOuterIterations(maxouterit),
SinglePrecGrid4(_sp_grid4), SinglePrecGrid4(_sp_grid4),
SinglePrecGrid5(_sp_grid5), SinglePrecGrid5(_sp_grid5),
OuterLoopNormMult(100.) OuterLoopNormMult(100.)
{ {
/* Debugging instances of objects; references are stored /* Debugging instances of objects; references are stored
std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpF " <<std::hex<< &LinOpF<<std::dec <<std::endl; std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpF " <<std::hex<< &LinOpF<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpD " <<std::hex<< &LinOpD<<std::dec <<std::endl; std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpD " <<std::hex<< &LinOpD<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpF " <<std::hex<< &FermOpF<<std::dec <<std::endl; std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpF " <<std::hex<< &FermOpF<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpD " <<std::hex<< &FermOpD<<std::dec <<std::endl; std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpD " <<std::hex<< &FermOpD<<std::dec <<std::endl;
*/ */
};
void operator()(LinearOperatorBase<FieldD> &LinOpU, const FieldD &src, FieldD &psi) {
std::cout << GridLogMessage << " Mixed precision CG wrapper operator() "<<std::endl;
SchurOperatorD * SchurOpU = static_cast<SchurOperatorD *>(&LinOpU);
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpU " <<std::hex<< &(SchurOpU->_Mat)<<std::dec <<std::endl;
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpD " <<std::hex<< &(LinOpD._Mat) <<std::dec <<std::endl;
// Assumption made in code to extract gauge field
// We could avoid storing LinopD reference alltogether ?
assert(&(SchurOpU->_Mat)==&(LinOpD._Mat));
////////////////////////////////////////////////////////////////////////////////////
// Must snarf a single precision copy of the gauge field in Linop_d argument
////////////////////////////////////////////////////////////////////////////////////
typedef typename FermionOperatorF::GaugeField GaugeFieldF;
typedef typename FermionOperatorF::GaugeLinkField GaugeLinkFieldF;
typedef typename FermionOperatorD::GaugeField GaugeFieldD;
typedef typename FermionOperatorD::GaugeLinkField GaugeLinkFieldD;
GridBase * GridPtrF = SinglePrecGrid4;
GridBase * GridPtrD = FermOpD.Umu.Grid();
GaugeFieldF U_f (GridPtrF);
GaugeLinkFieldF Umu_f(GridPtrF);
// std::cout << " Dim gauge field "<<GridPtrF->Nd()<<std::endl; // 4d
// std::cout << " Dim gauge field "<<GridPtrD->Nd()<<std::endl; // 4d
////////////////////////////////////////////////////////////////////////////////////
// Moving this to a Clone method of fermion operator would allow to duplicate the
// physics parameters and decrease gauge field copies
////////////////////////////////////////////////////////////////////////////////////
GaugeLinkFieldD Umu_d(GridPtrD);
for(int mu=0;mu<Nd*2;mu++){
Umu_d = PeekIndex<LorentzIndex>(FermOpD.Umu, mu);
precisionChange(Umu_f,Umu_d);
PokeIndex<LorentzIndex>(FermOpF.Umu, Umu_f, mu);
}
pickCheckerboard(Even,FermOpF.UmuEven,FermOpF.Umu);
pickCheckerboard(Odd ,FermOpF.UmuOdd ,FermOpF.Umu);
////////////////////////////////////////////////////////////////////////////////////
// Make a mixed precision conjugate gradient
////////////////////////////////////////////////////////////////////////////////////
MixedPrecisionConjugateGradient<FieldD,FieldF> MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD);
std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" <<std::endl;
MPCG(src,psi);
}
}; };
}};
void operator()(LinearOperatorBase<FieldD> &LinOpU, const FieldD &src, FieldD &psi) {
std::cout << GridLogMessage << " Mixed precision CG wrapper operator() "<<std::endl;
SchurOperatorD * SchurOpU = static_cast<SchurOperatorD *>(&LinOpU);
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpU " <<std::hex<< &(SchurOpU->_Mat)<<std::dec <<std::endl;
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpD " <<std::hex<< &(LinOpD._Mat) <<std::dec <<std::endl;
// Assumption made in code to extract gauge field
// We could avoid storing LinopD reference alltogether ?
assert(&(SchurOpU->_Mat)==&(LinOpD._Mat));
////////////////////////////////////////////////////////////////////////////////////
// Must snarf a single precision copy of the gauge field in Linop_d argument
////////////////////////////////////////////////////////////////////////////////////
typedef typename FermionOperatorF::GaugeField GaugeFieldF;
typedef typename FermionOperatorF::GaugeLinkField GaugeLinkFieldF;
typedef typename FermionOperatorD::GaugeField GaugeFieldD;
typedef typename FermionOperatorD::GaugeLinkField GaugeLinkFieldD;
GridBase * GridPtrF = SinglePrecGrid4;
GridBase * GridPtrD = FermOpD.Umu.Grid();
GaugeFieldF U_f (GridPtrF);
GaugeLinkFieldF Umu_f(GridPtrF);
// std::cout << " Dim gauge field "<<GridPtrF->Nd()<<std::endl; // 4d
// std::cout << " Dim gauge field "<<GridPtrD->Nd()<<std::endl; // 4d
////////////////////////////////////////////////////////////////////////////////////
// Moving this to a Clone method of fermion operator would allow to duplicate the
// physics parameters and decrease gauge field copies
////////////////////////////////////////////////////////////////////////////////////
GaugeLinkFieldD Umu_d(GridPtrD);
for(int mu=0;mu<Nd*2;mu++){
Umu_d = PeekIndex<LorentzIndex>(FermOpD.Umu, mu);
precisionChange(Umu_f,Umu_d);
PokeIndex<LorentzIndex>(FermOpF.Umu, Umu_f, mu);
}
pickCheckerboard(Even,FermOpF.UmuEven,FermOpF.Umu);
pickCheckerboard(Odd ,FermOpF.UmuOdd ,FermOpF.Umu);
////////////////////////////////////////////////////////////////////////////////////
// Make a mixed precision conjugate gradient
////////////////////////////////////////////////////////////////////////////////////
MixedPrecisionConjugateGradient<FieldD,FieldF> MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD);
std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" <<std::endl;
MPCG(src,psi);
}
};
NAMESPACE_END(Grid);
int main(int argc, char **argv) { int main(int argc, char **argv) {
using namespace Grid; using namespace Grid;
using namespace Grid::QCD;
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
int threads = GridThread::GetThreads(); int threads = GridThread::GetThreads();
@ -167,12 +167,12 @@ int main(int argc, char **argv) {
typedef typename FermionActionF::FermionField FermionFieldF; typedef typename FermionActionF::FermionField FermionFieldF;
typedef Grid::XmlReader Serialiser; typedef Grid::XmlReader Serialiser;
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// typedef GenericHMCRunner<LeapFrog> HMCWrapper; // typedef GenericHMCRunner<LeapFrog> HMCWrapper;
// typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; // typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
typedef GenericHMCRunner<ForceGradient> HMCWrapper; typedef GenericHMCRunner<ForceGradient> HMCWrapper;
HMCparameters HMCparams; HMCparameters HMCparams;
{ {
@ -184,7 +184,7 @@ int main(int argc, char **argv) {
// Grid from the command line arguments --grid and --mpi // Grid from the command line arguments --grid and --mpi
TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
CheckpointerParameters CPparams; CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_EODWF_lat"; CPparams.config_prefix = "ckpoint_EODWF_lat";
CPparams.rng_prefix = "ckpoint_EODWF_rng"; CPparams.rng_prefix = "ckpoint_EODWF_rng";
@ -198,7 +198,7 @@ int main(int argc, char **argv) {
TheHMC.Resources.SetRNGSeeds(RNGpar); TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables // Construct observables
// here there is too much indirection // here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs; typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>(); TheHMC.Resources.AddObservable<PlaqObs>();
////////////////////////////////////////////// //////////////////////////////////////////////
@ -209,7 +209,7 @@ int main(int argc, char **argv) {
Real strange_mass = 0.02144; Real strange_mass = 0.02144;
Real pv_mass = 1.0; Real pv_mass = 1.0;
RealD M5 = 1.8; RealD M5 = 1.8;
RealD b = 1.5; RealD b = 1.5;
RealD c = 0.5; RealD c = 0.5;
// Copied from paper // Copied from paper
@ -222,7 +222,7 @@ int main(int argc, char **argv) {
/////////////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////////////////
//Bad choices with large dH. Equalising force L2 norm was not wise. //Bad choices with large dH. Equalising force L2 norm was not wise.
/////////////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////////////////
//std::vector<Real> hasenbusch({ 0.03, 0.2, 0.3, 0.5, 0.8 }); //std::vector<Real> hasenbusch({ 0.03, 0.2, 0.3, 0.5, 0.8 });
//std::vector<Real> hasenbusch({ 0.05, 0.2, 0.4, 0.6, 0.8 }); //std::vector<Real> hasenbusch({ 0.05, 0.2, 0.4, 0.6, 0.8 });
auto GridPtr = TheHMC.Resources.GetCartesian(); auto GridPtr = TheHMC.Resources.GetCartesian();
@ -249,7 +249,7 @@ int main(int argc, char **argv) {
std::vector<Complex> boundary = {1,1,1,-1}; std::vector<Complex> boundary = {1,1,1,-1};
FermionAction::ImplParams Params(boundary); FermionAction::ImplParams Params(boundary);
FermionActionF::ImplParams ParamsF(boundary); FermionActionF::ImplParams ParamsF(boundary);
double ActionStoppingCondition = 1e-10; double ActionStoppingCondition = 1e-10;
double DerivativeStoppingCondition = 1e-7; double DerivativeStoppingCondition = 1e-7;
double MaxCGIterations = 30000; double MaxCGIterations = 30000;
@ -280,7 +280,7 @@ int main(int argc, char **argv) {
OFRp.degree = 12; OFRp.degree = 12;
OFRp.precision= 50; OFRp.precision= 50;
MobiusEOFAFermionR Strange_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c); MobiusEOFAFermionR Strange_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c);
MobiusEOFAFermionF Strange_Op_LF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c); MobiusEOFAFermionF Strange_Op_LF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c);
MobiusEOFAFermionR Strange_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c); MobiusEOFAFermionR Strange_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c);
@ -298,51 +298,51 @@ int main(int argc, char **argv) {
LinearOperatorEOFAF Strange_LinOp_RF(Strange_Op_RF); LinearOperatorEOFAF Strange_LinOp_RF(Strange_Op_RF);
MxPCG_EOFA ActionCGL(ActionStoppingCondition, MxPCG_EOFA ActionCGL(ActionStoppingCondition,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
Strange_Op_LF,Strange_Op_L, Strange_Op_LF,Strange_Op_L,
Strange_LinOp_LF,Strange_LinOp_L); Strange_LinOp_LF,Strange_LinOp_L);
MxPCG_EOFA DerivativeCGL(DerivativeStoppingCondition, MxPCG_EOFA DerivativeCGL(DerivativeStoppingCondition,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
Strange_Op_LF,Strange_Op_L, Strange_Op_LF,Strange_Op_L,
Strange_LinOp_LF,Strange_LinOp_L); Strange_LinOp_LF,Strange_LinOp_L);
MxPCG_EOFA ActionCGR(ActionStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> MxPCG_EOFA ActionCGR(ActionStoppingCondition,
EOFA(Strange_Op_L, Strange_Op_R, MX_inner,
ActionCG, MaxCGIterations,
ActionCGL, ActionCGR, GridPtrF,
DerivativeCGL, DerivativeCGR, FrbGridF,
OFRp, true); Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy>
EOFA(Strange_Op_L, Strange_Op_R,
ActionCG,
ActionCGL, ActionCGR,
DerivativeCGL, DerivativeCGR,
OFRp, true);
#else #else
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy>
EOFA(Strange_Op_L, Strange_Op_R, EOFA(Strange_Op_L, Strange_Op_R,
ActionCG, ActionCG,
ActionCG, ActionCG, ActionCG, ActionCG,
ActionCG, ActionCG, ActionCG, ActionCG,
// DerivativeCG, DerivativeCG, // DerivativeCG, DerivativeCG,
OFRp, true); OFRp, true);
#endif #endif
Level1.push_back(&EOFA); Level1.push_back(&EOFA);
@ -373,7 +373,7 @@ int main(int argc, char **argv) {
std::vector<MxPCG *> MPCG; std::vector<MxPCG *> MPCG;
std::vector<FermionActionF *> DenominatorsF; std::vector<FermionActionF *> DenominatorsF;
std::vector<LinearOperatorD *> LinOpD; std::vector<LinearOperatorD *> LinOpD;
std::vector<LinearOperatorF *> LinOpF; std::vector<LinearOperatorF *> LinOpF;
for(int h=0;h<n_hasenbusch+1;h++){ for(int h=0;h<n_hasenbusch+1;h++){
@ -395,20 +395,20 @@ int main(int argc, char **argv) {
double conv = DerivativeStoppingCondition; double conv = DerivativeStoppingCondition;
if (h<3) conv= DerivativeStoppingConditionLoose; // Relax on first two hasenbusch factors if (h<3) conv= DerivativeStoppingConditionLoose; // Relax on first two hasenbusch factors
MPCG.push_back(new MxPCG(conv, MPCG.push_back(new MxPCG(conv,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
*DenominatorsF[h],*Denominators[h], *DenominatorsF[h],*Denominators[h],
*LinOpF[h], *LinOpD[h]) ); *LinOpF[h], *LinOpD[h]) );
ActionMPCG.push_back(new MxPCG(ActionStoppingCondition, ActionMPCG.push_back(new MxPCG(ActionStoppingCondition,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
*DenominatorsF[h],*Denominators[h], *DenominatorsF[h],*Denominators[h],
*LinOpF[h], *LinOpD[h]) ); *LinOpF[h], *LinOpD[h]) );
// Heatbath not mixed yet. As inverts numerators not so important as raised mass. // Heatbath not mixed yet. As inverts numerators not so important as raised mass.
Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],ActionCG)); Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],ActionCG));

23
HMC/Mobius2p1fRHMC.cc
View File

@ -31,7 +31,6 @@ directory
int main(int argc, char **argv) { int main(int argc, char **argv) {
using namespace Grid; using namespace Grid;
using namespace Grid::QCD;
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
int threads = GridThread::GetThreads(); int threads = GridThread::GetThreads();
@ -44,18 +43,18 @@ int main(int argc, char **argv) {
typedef typename FermionAction::FermionField FermionField; typedef typename FermionAction::FermionField FermionField;
typedef Grid::XmlReader Serialiser; typedef Grid::XmlReader Serialiser;
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
IntegratorParameters MD; IntegratorParameters MD;
// typedef GenericHMCRunner<LeapFrog> HMCWrapper; // typedef GenericHMCRunner<LeapFrog> HMCWrapper;
// MD.name = std::string("Leap Frog"); // MD.name = std::string("Leap Frog");
// typedef GenericHMCRunner<ForceGradient> HMCWrapper; // typedef GenericHMCRunner<ForceGradient> HMCWrapper;
// MD.name = std::string("Force Gradient"); // MD.name = std::string("Force Gradient");
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
MD.name = std::string("MinimumNorm2"); MD.name = std::string("MinimumNorm2");
MD.MDsteps = 20; MD.MDsteps = 20;
MD.trajL = 1.0; MD.trajL = 1.0;
HMCparameters HMCparams; HMCparameters HMCparams;
HMCparams.StartTrajectory = 30; HMCparams.StartTrajectory = 30;
HMCparams.Trajectories = 200; HMCparams.Trajectories = 200;
@ -68,7 +67,7 @@ int main(int argc, char **argv) {
// Grid from the command line arguments --grid and --mpi // Grid from the command line arguments --grid and --mpi
TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
CheckpointerParameters CPparams; CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_EODWF_lat"; CPparams.config_prefix = "ckpoint_EODWF_lat";
CPparams.rng_prefix = "ckpoint_EODWF_rng"; CPparams.rng_prefix = "ckpoint_EODWF_rng";
@ -82,7 +81,7 @@ int main(int argc, char **argv) {
TheHMC.Resources.SetRNGSeeds(RNGpar); TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables // Construct observables
// here there is too much indirection // here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs; typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>(); TheHMC.Resources.AddObservable<PlaqObs>();
////////////////////////////////////////////// //////////////////////////////////////////////
@ -93,11 +92,11 @@ int main(int argc, char **argv) {
Real strange_mass = 0.04; Real strange_mass = 0.04;
Real pv_mass = 1.0; Real pv_mass = 1.0;
RealD M5 = 1.8; RealD M5 = 1.8;
RealD b = 1.0; RealD b = 1.0;
RealD c = 0.0; RealD c = 0.0;
// FIXME: // FIXME:
// Same in MC and MD // Same in MC and MD
// Need to mix precision too // Need to mix precision too
OneFlavourRationalParams OFRp; OneFlavourRationalParams OFRp;
OFRp.lo = 4.0e-3; OFRp.lo = 4.0e-3;
@ -122,7 +121,7 @@ int main(int argc, char **argv) {
// These lines are unecessary if BC are all periodic // These lines are unecessary if BC are all periodic
std::vector<Complex> boundary = {1,1,1,-1}; std::vector<Complex> boundary = {1,1,1,-1};
FermionAction::ImplParams Params(boundary); FermionAction::ImplParams Params(boundary);
double StoppingCondition = 1e-10; double StoppingCondition = 1e-10;
double MaxCGIterations = 30000; double MaxCGIterations = 30000;
ConjugateGradient<FermionField> CG(StoppingCondition,MaxCGIterations); ConjugateGradient<FermionField> CG(StoppingCondition,MaxCGIterations);

1
Hadrons/Utilities/HadronsXmlValidate.cc
View File

@ -29,7 +29,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Hadrons/Application.hpp> #include <Hadrons/Application.hpp>
using namespace Grid; using namespace Grid;
using namespace QCD;
using namespace Hadrons; using namespace Hadrons;
int main(int argc, char *argv[]) int main(int argc, char *argv[])

1
tests/hadrons/Test_diskvector.cc
View File

@ -28,7 +28,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Hadrons/DiskVector.hpp> #include <Hadrons/DiskVector.hpp>
using namespace Grid; using namespace Grid;
using namespace Grid::QCD;
using namespace Grid::Hadrons; using namespace Grid::Hadrons;
GRID_SERIALIZABLE_ENUM(Enum, undef, red, 1, blue, 2, green, 3); GRID_SERIALIZABLE_ENUM(Enum, undef, red, 1, blue, 2, green, 3);

1
tests/hmc/Test_hmc_Mobius2p1f.cc
View File

@ -31,7 +31,6 @@ directory
int main(int argc, char **argv) { int main(int argc, char **argv) {
using namespace Grid; using namespace Grid;
using namespace Grid::QCD;
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
int threads = GridThread::GetThreads(); int threads = GridThread::GetThreads();

190
tests/qdpxx/Test_qdpxx_loops_staples.cc
View File

@ -27,15 +27,15 @@ Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
/* END LEGAL */ /* END LEGAL */
#include <Grid/Grid.h> #include <Grid/Grid.h>
double calc_grid_p (Grid::QCD::LatticeGaugeField & lat); double calc_grid_p (Grid::LatticeGaugeField & lat);
double calc_chroma_p (Grid::QCD::LatticeGaugeField & lat); double calc_chroma_p (Grid::LatticeGaugeField & lat);
double calc_grid_r (Grid::QCD::LatticeGaugeField & lat); double calc_grid_r (Grid::LatticeGaugeField & lat);
double calc_grid_IW (Grid::QCD::LatticeGaugeField & lat); double calc_grid_IW (Grid::LatticeGaugeField & lat);
double calc_grid_r_dir (Grid::QCD::LatticeGaugeField & lat); double calc_grid_r_dir (Grid::LatticeGaugeField & lat);
double calc_chroma_r (Grid::QCD::LatticeGaugeField & lat); double calc_chroma_r (Grid::LatticeGaugeField & lat);
double calc_chroma_IW (Grid::QCD::LatticeGaugeField & lat); double calc_chroma_IW (Grid::LatticeGaugeField & lat);
void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu); void check_grid_r_staple(Grid::LatticeGaugeField & Umu);
void check_grid_p_staple(Grid::QCD::LatticeGaugeField & Umu); void check_grid_p_staple(Grid::LatticeGaugeField & Umu);
const double beta=2.6; const double beta=2.6;
const double c1=-0.331; const double c1=-0.331;
@ -53,10 +53,10 @@ public:
typedef multi1d<LatticeColorMatrix> U; typedef multi1d<LatticeColorMatrix> U;
static void ImportGauge(Grid::QCD::LatticeGaugeField & gr, static void ImportGauge(Grid::LatticeGaugeField & gr,
QDP::multi1d<QDP::LatticeColorMatrix> & ch) QDP::multi1d<QDP::LatticeColorMatrix> & ch)
{ {
Grid::QCD::LorentzColourMatrix LCM; Grid::LorentzColourMatrix LCM;
Grid::Complex cc; Grid::Complex cc;
QDP::ColorMatrix cm; QDP::ColorMatrix cm;
QDP::Complex c; QDP::Complex c;
@ -127,9 +127,9 @@ int main (int argc,char **argv )
* Setup Grid * Setup Grid
*********************************************************/ *********************************************************/
Grid::Grid_init(&argc,&argv); Grid::Grid_init(&argc,&argv);
Grid::GridCartesian * UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(), Grid::GridCartesian * UGrid = Grid::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
Grid::GridDefaultSimd(Grid::QCD::Nd,Grid::vComplex::Nsimd()), Grid::GridDefaultSimd(Grid::Nd,Grid::vComplex::Nsimd()),
Grid::GridDefaultMpi()); Grid::GridDefaultMpi());
std::vector<int> gd = UGrid->GlobalDimensions(); std::vector<int> gd = UGrid->GlobalDimensions();
QDP::multi1d<int> nrow(QDP::Nd); QDP::multi1d<int> nrow(QDP::Nd);
@ -138,7 +138,7 @@ int main (int argc,char **argv )
QDP::Layout::setLattSize(nrow); QDP::Layout::setLattSize(nrow);
QDP::Layout::create(); QDP::Layout::create();
Grid::QCD::LatticeGaugeField lat(UGrid); Grid::LatticeGaugeField lat(UGrid);
double s_grid = calc_grid_p (lat); double s_grid = calc_grid_p (lat);
@ -181,7 +181,7 @@ int main (int argc,char **argv )
Chroma::finalize(); Chroma::finalize();
} }
double calc_chroma_IW(Grid::QCD::LatticeGaugeField & lat) double calc_chroma_IW(Grid::LatticeGaugeField & lat)
{ {
typedef QDP::multi1d<QDP::LatticeColorMatrix> U; typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
@ -203,7 +203,7 @@ double calc_chroma_IW(Grid::QCD::LatticeGaugeField & lat)
return s; return s;
} }
double calc_chroma_r(Grid::QCD::LatticeGaugeField & lat) double calc_chroma_r(Grid::LatticeGaugeField & lat)
{ {
typedef QDP::multi1d<QDP::LatticeColorMatrix> U; typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
@ -245,7 +245,7 @@ double calc_chroma_r(Grid::QCD::LatticeGaugeField & lat)
// action = beta * Nd*Nd-1*vol*0.5 - beta * Nd*Nd-1*vol*0.5*plaq // action = beta * Nd*Nd-1*vol*0.5 - beta * Nd*Nd-1*vol*0.5*plaq
// //
// plaq == sumplaq * 2/(V*Nd*(Nd-1)*Nc) // plaq == sumplaq * 2/(V*Nd*(Nd-1)*Nc)
double calc_chroma_p(Grid::QCD::LatticeGaugeField & lat) double calc_chroma_p(Grid::LatticeGaugeField & lat)
{ {
typedef QDP::multi1d<QDP::LatticeColorMatrix> U; typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
@ -270,60 +270,60 @@ double calc_chroma_p(Grid::QCD::LatticeGaugeField & lat)
double calc_grid_p(Grid::QCD::LatticeGaugeField & Umu) double calc_grid_p(Grid::LatticeGaugeField & Umu)
{ {
std::vector<int> seeds4({1,2,3,4}); std::vector<int> seeds4({1,2,3,4});
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid(); Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
Grid::QCD::SU3::HotConfiguration(RNG4,Umu); Grid::SU3::HotConfiguration(RNG4,Umu);
Grid::QCD::LatticeColourMatrix tmp(UGrid); Grid::LatticeColourMatrix tmp(UGrid);
tmp = Grid::zero; tmp = Grid::zero;
Grid::QCD::PokeIndex<LorentzIndex>(Umu,tmp,2); Grid::PokeIndex<LorentzIndex>(Umu,tmp,2);
Grid::QCD::PokeIndex<LorentzIndex>(Umu,tmp,3); Grid::PokeIndex<LorentzIndex>(Umu,tmp,3);
Grid::QCD::WilsonGaugeActionR Wilson(beta); // Just take beta = 1.0 Grid::WilsonGaugeActionR Wilson(beta); // Just take beta = 1.0
return Wilson.S(Umu); return Wilson.S(Umu);
} }
double calc_grid_r(Grid::QCD::LatticeGaugeField & Umu) double calc_grid_r(Grid::LatticeGaugeField & Umu)
{ {
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid(); Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
Grid::QCD::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0 Grid::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
return Wilson.S(Umu); return Wilson.S(Umu);
} }
double calc_grid_IW(Grid::QCD::LatticeGaugeField & Umu) double calc_grid_IW(Grid::LatticeGaugeField & Umu)
{ {
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid(); Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
Grid::QCD::IwasakiGaugeActionR Iwasaki(beta); Grid::IwasakiGaugeActionR Iwasaki(beta);
return Iwasaki.S(Umu); return Iwasaki.S(Umu);
} }
double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu) double calc_grid_r_dir(Grid::LatticeGaugeField & Umu)
{ {
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid(); Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
std::vector<Grid::QCD::LatticeColourMatrix> U(4,UGrid); std::vector<Grid::LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){ for(int mu=0;mu<Nd;mu++){
U[mu] = Grid::PeekIndex<LorentzIndex>(Umu,mu); U[mu] = Grid::PeekIndex<LorentzIndex>(Umu,mu);
} }
Grid::QCD::LatticeComplex rect(UGrid); Grid::LatticeComplex rect(UGrid);
Grid::QCD::TComplex trect; Grid::TComplex trect;
Grid::QCD::Complex crect; Grid::Complex crect;
Grid::RealD rrect; Grid::RealD rrect;
Grid::RealD vol = UGrid->gSites(); Grid::RealD vol = UGrid->gSites();
for(int mu=0;mu<Grid::QCD::Nd;mu++){ for(int mu=0;mu<Grid::Nd;mu++){
for(int nu=0;nu<Grid::QCD::Nd;nu++){ for(int nu=0;nu<Grid::Nd;nu++){
if ( mu!=nu ) { if ( mu!=nu ) {
Grid::QCD::ColourWilsonLoops::traceDirRectangle(rect,U,mu,nu); Grid::ColourWilsonLoops::traceDirRectangle(rect,U,mu,nu);
trect = Grid::sum(rect); trect = Grid::sum(rect);
crect = Grid::TensorRemove(trect); crect = Grid::TensorRemove(trect);
rrect = real(crect); rrect = real(crect);
@ -335,9 +335,9 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// Staple test // Staple test
Peter.Start(); Peter.Start();
{ {
Grid::QCD::LatticeColourMatrix Stap(UGrid); Grid::LatticeColourMatrix Stap(UGrid);
Grid::QCD::LatticeComplex SumTrStap(UGrid); Grid::LatticeComplex SumTrStap(UGrid);
Grid::QCD::LatticeComplex TrStap(UGrid); Grid::LatticeComplex TrStap(UGrid);
/* /*
* Make staple for loops centered at coor of link ; this one is ok. // | * Make staple for loops centered at coor of link ; this one is ok. // |
@ -346,10 +346,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// __ ___ // __ ___
// | __ | // | __ |
Stap = Stap =
Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[mu],mu, Grid::Cshift(Grid::PeriodicBC::CovShiftForward (U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu, Grid::PeriodicBC::CovShiftForward (U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1); Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap); TrStap = Grid::trace (U[mu]*Stap);
@ -364,10 +364,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// __ // __
// |__ __ | // |__ __ |
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[mu],mu, Stap = Grid::Cshift(Grid::PeriodicBC::CovShiftForward (U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu, Grid::PeriodicBC::CovShiftBackward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1); Grid::PeriodicBC::CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap); TrStap = Grid::trace (U[mu]*Stap);
@ -379,10 +379,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// __ // __
// |__ __ | // |__ __ |
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu, Stap = Grid::Cshift(Grid::PeriodicBC::CovShiftBackward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[mu])))) , mu, 1); Grid::PeriodicBC::CovShiftForward(U[nu],nu,U[mu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap); TrStap = Grid::trace (U[mu]*Stap);
@ -395,10 +395,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// __ ___ // __ ___
// |__ | // |__ |
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu, Stap = Grid::Cshift(Grid::PeriodicBC::CovShiftForward (U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1); Grid::PeriodicBC::CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap); TrStap = Grid::trace (U[mu]*Stap);
@ -418,12 +418,12 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
* Make staple for loops centered at coor of link ; this one is ok. // | * Make staple for loops centered at coor of link ; this one is ok. // |
*/ */
// Stap = // Stap =
// Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]),mu,1)* // ->|| // Grid::Cshift(Grid::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]),mu,1)* // ->||
// Grid::adj(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[mu]))) ; // Grid::adj(Grid::PeriodicBC::CovShiftForward(U[nu],nu,Grid::PeriodicBC::CovShiftForward(U[nu],nu,U[mu]))) ;
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu, Stap = Grid::Cshift(Grid::PeriodicBC::CovShiftForward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu, Grid::PeriodicBC::CovShiftForward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu, Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1); Grid::PeriodicBC::CovShiftBackward(U[nu],nu, Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap); TrStap = Grid::trace (U[mu]*Stap);
SumTrStap += TrStap; SumTrStap += TrStap;
@ -440,10 +440,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// | | // | |
// -- // --
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu, Stap = Grid::Cshift(Grid::PeriodicBC::CovShiftBackward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu, Grid::PeriodicBC::CovShiftBackward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu,U[nu])))) , mu, 1); Grid::PeriodicBC::CovShiftForward (U[nu],nu,U[nu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap); TrStap = Grid::trace (U[mu]*Stap);
trect = Grid::sum(TrStap); trect = Grid::sum(TrStap);
@ -459,20 +459,20 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
Peter.Stop(); Peter.Stop();
Azusa.Start(); Azusa.Start();
{ {
Grid::QCD::LatticeComplex RectPlaq_d(UGrid); Grid::LatticeComplex RectPlaq_d(UGrid);
Grid::QCD::LatticeColourMatrix ds_U(UGrid); Grid::LatticeColourMatrix ds_U(UGrid);
Grid::QCD::LatticeColourMatrix left_2(UGrid); Grid::LatticeColourMatrix left_2(UGrid);
Grid::QCD::LatticeColourMatrix upper_l(UGrid); Grid::LatticeColourMatrix upper_l(UGrid);
Grid::QCD::LatticeColourMatrix upper_staple(UGrid); Grid::LatticeColourMatrix upper_staple(UGrid);
Grid::QCD::LatticeColourMatrix down_l(UGrid); Grid::LatticeColourMatrix down_l(UGrid);
Grid::QCD::LatticeColourMatrix down_staple(UGrid); Grid::LatticeColourMatrix down_staple(UGrid);
Grid::QCD::LatticeColourMatrix tmp(UGrid); Grid::LatticeColourMatrix tmp(UGrid);
// 2 (mu)x1(nu) // 2 (mu)x1(nu)
left_2= Grid::QCD::PeriodicBC::CovShiftForward(U[mu],mu,U[mu]); // Umu(x) Umu(x+mu) left_2= Grid::PeriodicBC::CovShiftForward(U[mu],mu,U[mu]); // Umu(x) Umu(x+mu)
tmp=Grid::Cshift(U[nu],mu,2); // Unu(x+2mu) tmp=Grid::Cshift(U[nu],mu,2); // Unu(x+2mu)
upper_l= Grid::QCD::PeriodicBC::CovShiftForward(tmp,nu,Grid::adj(left_2)); // Unu(x+2mu) Umu^dag(x+mu+nu) Umu^dag(x+nu) upper_l= Grid::PeriodicBC::CovShiftForward(tmp,nu,Grid::adj(left_2)); // Unu(x+2mu) Umu^dag(x+mu+nu) Umu^dag(x+nu)
// __ __ // __ __
// = | // = |
@ -546,9 +546,9 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// _ // _
// | | // | |
// | | // | |
Grid::QCD::LatticeColourMatrix up2= Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]); Grid::LatticeColourMatrix up2= Grid::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]);
upper_l= Grid::QCD::PeriodicBC::CovShiftForward(Grid::Cshift(up2,mu,1),nu,Grid::Cshift(adj(U[mu]),nu,1)); upper_l= Grid::PeriodicBC::CovShiftForward(Grid::Cshift(up2,mu,1),nu,Grid::Cshift(adj(U[mu]),nu,1));
ds_U= upper_l*Grid::adj(up2); ds_U= upper_l*Grid::adj(up2);
RectPlaq_d = Grid::trace(U[mu]*ds_U); RectPlaq_d = Grid::trace(U[mu]*ds_U);
@ -569,7 +569,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
downer_l= | downer_l= |
(x)<----V (x)<----V
*/ */
down_l= Grid::adj(Grid::QCD::PeriodicBC::CovShiftForward(U[mu],mu,up2)); //downer_l down_l= Grid::adj(Grid::PeriodicBC::CovShiftForward(U[mu],mu,up2)); //downer_l
/* /*
^ | ^ |
down_staple = | V down_staple = | V
@ -601,23 +601,23 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
} }
}} }}
Grid::QCD::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0 Grid::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
return Wilson.S(Umu); return Wilson.S(Umu);
}; };
void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu) void check_grid_r_staple(Grid::LatticeGaugeField & Umu)
{ {
std::vector<int> seeds4({1,2,3,4}); std::vector<int> seeds4({1,2,3,4});
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid(); Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
Grid::QCD::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0 Grid::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
Grid::QCD::LatticeColourMatrix staple(UGrid); Grid::LatticeColourMatrix staple(UGrid);
Grid::QCD::LatticeColourMatrix link(UGrid); Grid::LatticeColourMatrix link(UGrid);
Grid::QCD::LatticeComplex Traced(UGrid); Grid::LatticeComplex Traced(UGrid);
Grid::Complex Rplaq(0.0); Grid::Complex Rplaq(0.0);
for(int mu=0;mu<Nd;mu++){ for(int mu=0;mu<Nd;mu++){
@ -630,12 +630,12 @@ void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu)
// Vol as for each site // Vol as for each site
Grid::RealD RectScale(1.0/vol/12.0/6.0/3.0); Grid::RealD RectScale(1.0/vol/12.0/6.0/3.0);
Grid::QCD::ColourWilsonLoops::RectStaple(staple,Umu,mu); Grid::ColourWilsonLoops::RectStaple(staple,Umu,mu);
link = Grid::QCD::PeekIndex<LorentzIndex>(Umu,mu); link = Grid::PeekIndex<LorentzIndex>(Umu,mu);
Traced = Grid::trace( link*staple) * RectScale; Traced = Grid::trace( link*staple) * RectScale;
Grid::QCD::TComplex Tp = Grid::sum(Traced); Grid::TComplex Tp = Grid::sum(Traced);
Grid::Complex p = Grid::TensorRemove(Tp); Grid::Complex p = Grid::TensorRemove(Tp);
std::cout<< "Rect from RectStaple "<<p<<std::endl; std::cout<< "Rect from RectStaple "<<p<<std::endl;
@ -645,18 +645,18 @@ void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu)
std::cout<< "Rect from RectStaple "<<Rplaq<<std::endl; std::cout<< "Rect from RectStaple "<<Rplaq<<std::endl;
} }
void check_grid_p_staple(Grid::QCD::LatticeGaugeField & Umu) void check_grid_p_staple(Grid::LatticeGaugeField & Umu)
{ {
std::vector<int> seeds4({1,2,3,4}); std::vector<int> seeds4({1,2,3,4});
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid(); Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
Grid::QCD::PlaqPlusRectangleActionR Wilson(1.0,0.0); // Just take c1 = 0.0 Grid::PlaqPlusRectangleActionR Wilson(1.0,0.0); // Just take c1 = 0.0
Grid::QCD::LatticeColourMatrix staple(UGrid); Grid::LatticeColourMatrix staple(UGrid);
Grid::QCD::LatticeColourMatrix link(UGrid); Grid::LatticeColourMatrix link(UGrid);
Grid::QCD::LatticeComplex Traced(UGrid); Grid::LatticeComplex Traced(UGrid);
Grid::Complex plaq(0.0); Grid::Complex plaq(0.0);
for(int mu=0;mu<Nd;mu++){ for(int mu=0;mu<Nd;mu++){
@ -669,12 +669,12 @@ void check_grid_p_staple(Grid::QCD::LatticeGaugeField & Umu)
// Vol as for each site // Vol as for each site
Grid::RealD Scale(1.0/vol/12.0/2.0/3.0); Grid::RealD Scale(1.0/vol/12.0/2.0/3.0);
Grid::QCD::ColourWilsonLoops::Staple(staple,Umu,mu); Grid::ColourWilsonLoops::Staple(staple,Umu,mu);
link = Grid::QCD::PeekIndex<LorentzIndex>(Umu,mu); link = Grid::PeekIndex<LorentzIndex>(Umu,mu);
Traced = Grid::trace( link*staple) * Scale; Traced = Grid::trace( link*staple) * Scale;
Grid::QCD::TComplex Tp = Grid::sum(Traced); Grid::TComplex Tp = Grid::sum(Traced);
Grid::Complex p = Grid::TensorRemove(Tp); Grid::Complex p = Grid::TensorRemove(Tp);
std::cout<< "Plaq from PlaqStaple "<<p<<std::endl; std::cout<< "Plaq from PlaqStaple "<<p<<std::endl;

60
tests/qdpxx/Test_qdpxx_munprec.cc
View File

@ -52,8 +52,8 @@ enum ChromaAction {
HtContFracZolo HtContFracZolo
}; };
void calc_grid (ChromaAction action,Grid::QCD::LatticeGaugeField & lat, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res,int dag); void calc_grid (ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag);
void calc_chroma (ChromaAction action,Grid::QCD::LatticeGaugeField & lat, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res,int dag); void calc_chroma (ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag);
#include <chroma.h> #include <chroma.h>
#include <actions/ferm/invert/syssolver_linop_cg_array.h> #include <actions/ferm/invert/syssolver_linop_cg_array.h>
@ -71,10 +71,10 @@ public:
typedef LatticeFermion T4; typedef LatticeFermion T4;
typedef multi1d<LatticeFermion> T5; typedef multi1d<LatticeFermion> T5;
static void ImportGauge(Grid::QCD::LatticeGaugeField & gr, static void ImportGauge(Grid::LatticeGaugeField & gr,
QDP::multi1d<QDP::LatticeColorMatrix> & ch) QDP::multi1d<QDP::LatticeColorMatrix> & ch)
{ {
Grid::QCD::LorentzColourMatrix LCM; Grid::LorentzColourMatrix LCM;
Grid::Complex cc; Grid::Complex cc;
QDP::ColorMatrix cm; QDP::ColorMatrix cm;
QDP::Complex c; QDP::Complex c;
@ -112,10 +112,10 @@ public:
}}}} }}}}
} }
static void ImportFermion(Grid::QCD::LatticeFermion & gr, static void ImportFermion(Grid::LatticeFermion & gr,
QDP::multi1d<QDP::LatticeFermion> & ch ) QDP::multi1d<QDP::LatticeFermion> & ch )
{ {
Grid::QCD::SpinColourVector F; Grid::SpinColourVector F;
Grid::Complex c; Grid::Complex c;
QDP::Fermion cF; QDP::Fermion cF;
@ -154,10 +154,10 @@ public:
QDP::pokeSite(ch[s],cF,cx); QDP::pokeSite(ch[s],cF,cx);
}}}}} }}}}}
} }
static void ExportFermion(Grid::QCD::LatticeFermion & gr, static void ExportFermion(Grid::LatticeFermion & gr,
QDP::multi1d<QDP::LatticeFermion> & ch ) QDP::multi1d<QDP::LatticeFermion> & ch )
{ {
Grid::QCD::SpinColourVector F; Grid::SpinColourVector F;
Grid::Complex c; Grid::Complex c;
QDP::Fermion cF; QDP::Fermion cF;
@ -384,9 +384,9 @@ int main (int argc,char **argv )
* Setup Grid * Setup Grid
*********************************************************/ *********************************************************/
Grid::Grid_init(&argc,&argv); Grid::Grid_init(&argc,&argv);
Grid::GridCartesian * UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(), Grid::GridCartesian * UGrid = Grid::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
Grid::GridDefaultSimd(Grid::QCD::Nd,Grid::vComplex::Nsimd()), Grid::GridDefaultSimd(Grid::Nd,Grid::vComplex::Nsimd()),
Grid::GridDefaultMpi()); Grid::GridDefaultMpi());
std::vector<int> gd = UGrid->GlobalDimensions(); std::vector<int> gd = UGrid->GlobalDimensions();
QDP::multi1d<int> nrow(QDP::Nd); QDP::multi1d<int> nrow(QDP::Nd);
@ -395,11 +395,11 @@ int main (int argc,char **argv )
QDP::Layout::setLattSize(nrow); QDP::Layout::setLattSize(nrow);
QDP::Layout::create(); QDP::Layout::create();
Grid::GridCartesian * FGrid = Grid::QCD::SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); Grid::GridCartesian * FGrid = Grid::SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
Grid::QCD::LatticeGaugeField lat(UGrid); Grid::LatticeGaugeField lat(UGrid);
Grid::QCD::LatticeFermion src(FGrid); Grid::LatticeFermion src(FGrid);
Grid::QCD::LatticeFermion res_chroma(FGrid); Grid::LatticeFermion res_chroma(FGrid);
Grid::QCD::LatticeFermion res_grid (FGrid); Grid::LatticeFermion res_grid (FGrid);
std::vector<ChromaAction> ActionList({ std::vector<ChromaAction> ActionList({
HtCayleyTanh, // Plain old DWF. HtCayleyTanh, // Plain old DWF.
@ -446,7 +446,7 @@ int main (int argc,char **argv )
Chroma::finalize(); Chroma::finalize();
} }
void calc_chroma(ChromaAction action,Grid::QCD::LatticeGaugeField & lat, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res,int dag) void calc_chroma(ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag)
{ {
QDP::multi1d<QDP::LatticeColorMatrix> u(4); QDP::multi1d<QDP::LatticeColorMatrix> u(4);
@ -483,7 +483,7 @@ void calc_chroma(ChromaAction action,Grid::QCD::LatticeGaugeField & lat, Grid::Q
void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res,int dag) void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag)
{ {
using namespace Grid; using namespace Grid;
; ;
@ -493,8 +493,8 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid(); Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
Grid::GridCartesian * FGrid = (Grid::GridCartesian *) src.Grid(); Grid::GridCartesian * FGrid = (Grid::GridCartesian *) src.Grid();
Grid::GridRedBlackCartesian * UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); Grid::GridRedBlackCartesian * UrbGrid = Grid::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
Grid::GridRedBlackCartesian * FrbGrid = Grid::QCD::SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); Grid::GridRedBlackCartesian * FrbGrid = Grid::SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
Grid::GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); Grid::GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
@ -502,12 +502,12 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
Grid::gaussian(RNG5,src); Grid::gaussian(RNG5,src);
Grid::gaussian(RNG5,res); Grid::gaussian(RNG5,res);
Grid::QCD::SU3::HotConfiguration(RNG4,Umu); Grid::SU3::HotConfiguration(RNG4,Umu);
/* /*
Grid::QCD::LatticeColourMatrix U(UGrid); Grid::LatticeColourMatrix U(UGrid);
U=Grid::zero; U=Grid::zero;
for(int nn=0;nn<Grid::QCD::Nd;nn++){ for(int nn=0;nn<Grid::Nd;nn++){
if ( nn>=4 ) { if ( nn>=4 ) {
Grid::PokeIndex<LorentzIndex>(Umu,U,nn); Grid::PokeIndex<LorentzIndex>(Umu,U,nn);
} }
@ -519,7 +519,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
if ( action == HtCayleyTanh ) { if ( action == HtCayleyTanh ) {
Grid::QCD::DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5); Grid::DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5);
std::cout << Grid::GridLogMessage <<" Calling domain wall multiply "<<std::endl; std::cout << Grid::GridLogMessage <<" Calling domain wall multiply "<<std::endl;
@ -535,7 +535,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
Grid::Real _b = 0.5*(mobius_scale +1.0); Grid::Real _b = 0.5*(mobius_scale +1.0);
Grid::Real _c = 0.5*(mobius_scale -1.0); Grid::Real _c = 0.5*(mobius_scale -1.0);
Grid::QCD::MobiusZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c,zolo_lo,zolo_hi); Grid::MobiusZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c,zolo_lo,zolo_hi);
std::cout << Grid::GridLogMessage <<" Calling mobius zolo multiply "<<std::endl; std::cout << Grid::GridLogMessage <<" Calling mobius zolo multiply "<<std::endl;
@ -549,7 +549,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
if ( action == HtCayleyZolo ) { if ( action == HtCayleyZolo ) {
Grid::QCD::ShamirZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi); Grid::ShamirZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
std::cout << Grid::GridLogMessage <<" Calling shamir zolo multiply "<<std::endl; std::cout << Grid::GridLogMessage <<" Calling shamir zolo multiply "<<std::endl;
@ -565,7 +565,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
if ( action == HmCayleyTanh ) { if ( action == HmCayleyTanh ) {
Grid::Real _b = 0.5*(mobius_scale +1.0); Grid::Real _b = 0.5*(mobius_scale +1.0);
Grid::Real _c = 0.5*(mobius_scale -1.0); Grid::Real _c = 0.5*(mobius_scale -1.0);
Grid::QCD::MobiusFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c); Grid::MobiusFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c);
std::cout << Grid::GridLogMessage <<" Calling mobius tanh multiply "<<std::endl; std::cout << Grid::GridLogMessage <<" Calling mobius tanh multiply "<<std::endl;
@ -581,7 +581,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
if ( action == HmCayleyTanh ) { if ( action == HmCayleyTanh ) {
Grid::QCD::ScaledShamirFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,mobius_scale); Grid::ScaledShamirFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,mobius_scale);
std::cout << Grid::GridLogMessage <<" Calling scaled shamir multiply "<<std::endl; std::cout << Grid::GridLogMessage <<" Calling scaled shamir multiply "<<std::endl;
@ -595,7 +595,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
if ( action == HwCayleyTanh ) { if ( action == HwCayleyTanh ) {
Grid::QCD::OverlapWilsonCayleyTanhFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0); Grid::OverlapWilsonCayleyTanhFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0);
if ( dag ) if ( dag )
D.Mdag(src,res); D.Mdag(src,res);
@ -607,7 +607,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
if ( action == HwCayleyZolo ) { if ( action == HwCayleyZolo ) {
Grid::QCD::OverlapWilsonCayleyZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi); Grid::OverlapWilsonCayleyZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
if ( dag ) if ( dag )
D.Mdag(src,res); D.Mdag(src,res);

24
tests/qdpxx/Test_qdpxx_stag.cc
View File

@ -30,8 +30,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
double mq=0.1; double mq=0.1;
typedef Grid::QCD::StaggeredImplR::FermionField FermionField; typedef Grid::StaggeredImplR::FermionField FermionField;
typedef Grid::QCD::LatticeGaugeField GaugeField; typedef Grid::LatticeGaugeField GaugeField;
void make_gauge (GaugeField & lat, FermionField &src); void make_gauge (GaugeField & lat, FermionField &src);
void calc_grid (GaugeField & lat, GaugeField & uthin,GaugeField & ufat, FermionField &src, FermionField &res,int dag); void calc_grid (GaugeField & lat, GaugeField & uthin,GaugeField & ufat, FermionField &src, FermionField &res,int dag);
@ -53,7 +53,7 @@ public:
static void ImportGauge(GaugeField & gr, static void ImportGauge(GaugeField & gr,
QDP::multi1d<QDP::LatticeColorMatrix> & ch) QDP::multi1d<QDP::LatticeColorMatrix> & ch)
{ {
Grid::QCD::LorentzColourMatrix LCM; Grid::LorentzColourMatrix LCM;
Grid::Complex cc; Grid::Complex cc;
QDP::ColorMatrix cm; QDP::ColorMatrix cm;
QDP::Complex c; QDP::Complex c;
@ -88,7 +88,7 @@ public:
static void ExportGauge(GaugeField & gr, static void ExportGauge(GaugeField & gr,
QDP::multi1d<QDP::LatticeColorMatrix> & ch) QDP::multi1d<QDP::LatticeColorMatrix> & ch)
{ {
Grid::QCD::LorentzColourMatrix LCM; Grid::LorentzColourMatrix LCM;
Grid::Complex cc; Grid::Complex cc;
QDP::ColorMatrix cm; QDP::ColorMatrix cm;
QDP::Complex c; QDP::Complex c;
@ -124,7 +124,7 @@ public:
static void ImportFermion(FermionField & gr, static void ImportFermion(FermionField & gr,
QDP::LatticeStaggeredFermion & ch ) QDP::LatticeStaggeredFermion & ch )
{ {
Grid::QCD::ColourVector F; Grid::ColourVector F;
Grid::Complex c; Grid::Complex c;
@ -157,7 +157,7 @@ public:
static void ExportFermion(FermionField & gr, static void ExportFermion(FermionField & gr,
QDP::LatticeStaggeredFermion & ch ) QDP::LatticeStaggeredFermion & ch )
{ {
Grid::QCD::ColourVector F; Grid::ColourVector F;
Grid::Complex c; Grid::Complex c;
std::vector<int> x(5); std::vector<int> x(5);
@ -222,9 +222,9 @@ int main (int argc,char **argv )
* Setup Grid * Setup Grid
*********************************************************/ *********************************************************/
Grid::Grid_init(&argc,&argv); Grid::Grid_init(&argc,&argv);
Grid::GridCartesian * UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(), Grid::GridCartesian * UGrid = Grid::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
Grid::GridDefaultSimd(Grid::QCD::Nd,Grid::vComplex::Nsimd()), Grid::GridDefaultSimd(Grid::Nd,Grid::vComplex::Nsimd()),
Grid::GridDefaultMpi()); Grid::GridDefaultMpi());
std::vector<int> gd = UGrid->GlobalDimensions(); std::vector<int> gd = UGrid->GlobalDimensions();
QDP::multi1d<int> nrow(QDP::Nd); QDP::multi1d<int> nrow(QDP::Nd);
@ -333,7 +333,7 @@ void make_gauge(GaugeField & Umu,FermionField &src)
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid(); Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
Grid::QCD::SU3::HotConfiguration(RNG4,Umu); Grid::SU3::HotConfiguration(RNG4,Umu);
Grid::gaussian(RNG4,src); Grid::gaussian(RNG4,src);
} }
@ -343,9 +343,9 @@ void calc_grid(GaugeField & Uthin, GaugeField & Utriple, GaugeField & Ufat, Ferm
; ;
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Uthin.Grid(); Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Uthin.Grid();
Grid::GridRedBlackCartesian * UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); Grid::GridRedBlackCartesian * UrbGrid = Grid::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
Grid::QCD::ImprovedStaggeredFermionR Dstag(Uthin,Utriple,Ufat,*UGrid,*UrbGrid,mq*2.0); Grid::ImprovedStaggeredFermionR Dstag(Uthin,Utriple,Ufat,*UGrid,*UrbGrid,mq*2.0);
std::cout << Grid::GridLogMessage <<" Calling Grid staggered multiply "<<std::endl; std::cout << Grid::GridLogMessage <<" Calling Grid staggered multiply "<<std::endl;

36
tests/qdpxx/Test_qdpxx_wilson.cc
View File

@ -35,8 +35,8 @@
double mq = 0.1; double mq = 0.1;
// Define Wilson Types // Define Wilson Types
typedef Grid::QCD::WilsonImplR::FermionField FermionField; typedef Grid::WilsonImplR::FermionField FermionField;
typedef Grid::QCD::LatticeGaugeField GaugeField; typedef Grid::LatticeGaugeField GaugeField;
enum ChromaAction enum ChromaAction
{ {
@ -56,7 +56,7 @@ public:
static void ImportGauge(GaugeField &gr, static void ImportGauge(GaugeField &gr,
QDP::multi1d<QDP::LatticeColorMatrix> &ch) QDP::multi1d<QDP::LatticeColorMatrix> &ch)
{ {
Grid::QCD::LorentzColourMatrix LCM; Grid::LorentzColourMatrix LCM;
Grid::Complex cc; Grid::Complex cc;
QDP::ColorMatrix cm; QDP::ColorMatrix cm;
QDP::Complex c; QDP::Complex c;
@ -101,7 +101,7 @@ public:
static void ExportGauge(GaugeField &gr, static void ExportGauge(GaugeField &gr,
QDP::multi1d<QDP::LatticeColorMatrix> &ch) QDP::multi1d<QDP::LatticeColorMatrix> &ch)
{ {
Grid::QCD::LorentzColourMatrix LCM; Grid::LorentzColourMatrix LCM;
Grid::Complex cc; Grid::Complex cc;
QDP::ColorMatrix cm; QDP::ColorMatrix cm;
QDP::Complex c; QDP::Complex c;
@ -145,10 +145,10 @@ public:
} }
// Specific for Wilson Fermions // Specific for Wilson Fermions
static void ImportFermion(Grid::QCD::LatticeFermion &gr, static void ImportFermion(Grid::LatticeFermion &gr,
QDP::LatticeFermion &ch) QDP::LatticeFermion &ch)
{ {
Grid::QCD::SpinColourVector F; Grid::SpinColourVector F;
Grid::Complex c; Grid::Complex c;
QDP::Fermion cF; QDP::Fermion cF;
@ -195,10 +195,10 @@ public:
} }
// Specific for 4d Wilson fermions // Specific for 4d Wilson fermions
static void ExportFermion(Grid::QCD::LatticeFermion &gr, static void ExportFermion(Grid::LatticeFermion &gr,
QDP::LatticeFermion &ch) QDP::LatticeFermion &ch)
{ {
Grid::QCD::SpinColourVector F; Grid::SpinColourVector F;
Grid::Complex c; Grid::Complex c;
QDP::Fermion cF; QDP::Fermion cF;
@ -342,19 +342,18 @@ void calc_chroma(ChromaAction action, GaugeField &lat, FermionField &src, Fermio
void make_gauge(GaugeField &Umu, FermionField &src) void make_gauge(GaugeField &Umu, FermionField &src)
{ {
using namespace Grid; using namespace Grid;
using namespace Grid::QCD;
std::vector<int> seeds4({1, 2, 3, 4}); std::vector<int> seeds4({1, 2, 3, 4});
Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid; Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid;
Grid::GridParallelRNG RNG4(UGrid); Grid::GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4); RNG4.SeedFixedIntegers(seeds4);
Grid::QCD::SU3::HotConfiguration(RNG4, Umu); Grid::SU3::HotConfiguration(RNG4, Umu);
// Fermion field // Fermion field
Grid::gaussian(RNG4, src); Grid::gaussian(RNG4, src);
/* /*
Grid::QCD::SpinColourVector F; Grid::SpinColourVector F;
Grid::Complex c; Grid::Complex c;
@ -391,13 +390,12 @@ void make_gauge(GaugeField &Umu, FermionField &src)
*/ */
} }
void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag) void calc_grid(ChromaAction action, Grid::LatticeGaugeField &Umu, Grid::LatticeFermion &src, Grid::LatticeFermion &res, int dag)
{ {
using namespace Grid; using namespace Grid;
using namespace Grid::QCD;
Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid; Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid;
Grid::GridRedBlackCartesian *UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); Grid::GridRedBlackCartesian *UrbGrid = Grid::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
Grid::RealD _mass = mq; Grid::RealD _mass = mq;
@ -409,7 +407,7 @@ void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD
anis.xi_0 = 2.0; anis.xi_0 = 2.0;
anis.nu = 1.0; anis.nu = 1.0;
WilsonImplParams iParam; WilsonImplParams iParam;
Grid::QCD::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, iParam, anis); Grid::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, iParam, anis);
std::cout << Grid::GridLogMessage << " Calling Grid Wilson Fermion multiply " << std::endl; std::cout << Grid::GridLogMessage << " Calling Grid Wilson Fermion multiply " << std::endl;
@ -430,7 +428,7 @@ void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD
anis.xi_0 = 2.0; anis.xi_0 = 2.0;
anis.nu = 1.0; anis.nu = 1.0;
WilsonImplParams CloverImplParam; WilsonImplParams CloverImplParam;
Grid::QCD::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw_r, _csw_t, anis, CloverImplParam); Grid::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw_r, _csw_t, anis, CloverImplParam);
Wf.ImportGauge(Umu); Wf.ImportGauge(Umu);
std::cout << Grid::GridLogMessage << " Calling Grid Wilson Clover Fermion multiply " << std::endl; std::cout << Grid::GridLogMessage << " Calling Grid Wilson Clover Fermion multiply " << std::endl;
@ -458,9 +456,9 @@ int main(int argc, char **argv)
* Setup Grid * Setup Grid
*********************************************************/ *********************************************************/
Grid::Grid_init(&argc, &argv); Grid::Grid_init(&argc, &argv);
Grid::GridCartesian *UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(), Grid::GridCartesian *UGrid = Grid::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
Grid::GridDefaultSimd(Grid::QCD::Nd, Grid::vComplex::Nsimd()), Grid::GridDefaultSimd(Grid::Nd, Grid::vComplex::Nsimd()),
Grid::GridDefaultMpi()); Grid::GridDefaultMpi());
std::vector<int> gd = UGrid->GlobalDimensions(); std::vector<int> gd = UGrid->GlobalDimensions();
QDP::multi1d<int> nrow(QDP::Nd); QDP::multi1d<int> nrow(QDP::Nd);

2
tests/smearing/Test_smearing.cc
View File

@ -30,8 +30,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
using namespace std; using namespace std;
using namespace Grid; using namespace Grid;
using namespace Grid::QCD;
int main (int argc, char ** argv) int main (int argc, char ** argv)
{ {

1
tests/solver/Test_wilsonclover_mg_lime.cc
View File

@ -32,7 +32,6 @@
using namespace std; using namespace std;
using namespace Grid; using namespace Grid;
using namespace Grid::QCD;
int main(int argc, char **argv) { int main(int argc, char **argv) {