portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00
Code Releases Activity

remove namespace QCD from directory HMC

Browse Source
This commit is contained in:
Nils Asmussen 2019-08-20 15:21:23 +01:00
parent ad01290545
commit 9210b0aa6e
4 changed files with 376 additions and 379 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
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();

333
HMC/Mobius2p1fEOFA.cc
View File

@ -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();
@ -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, MxPCG_EOFA ActionCGR(ActionStoppingCondition,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
Strange_Op_RF,Strange_Op_R, Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R); Strange_LinOp_RF,Strange_LinOp_R);
MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition, MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
Strange_Op_RF,Strange_Op_R, Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R); Strange_LinOp_RF,Strange_LinOp_R);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy>
EOFA(Strange_Op_L, Strange_Op_R, EOFA(Strange_Op_L, Strange_Op_R,
ActionCG, ActionCG,
ActionCGL, ActionCGR, ActionCGL, ActionCGR,
DerivativeCGL, DerivativeCGR, DerivativeCGL, DerivativeCGR,
OFRp, true); 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);
@ -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));

302
HMC/Mobius2p1fEOFA_F1.cc
View File

@ -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();
@ -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, MxPCG_EOFA ActionCGR(ActionStoppingCondition,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
Strange_Op_RF,Strange_Op_R, Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R); Strange_LinOp_RF,Strange_LinOp_R);
MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition, MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF, GridPtrF,
FrbGridF, FrbGridF,
Strange_Op_RF,Strange_Op_R, Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R); Strange_LinOp_RF,Strange_LinOp_R);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy>
EOFA(Strange_Op_L, Strange_Op_R, EOFA(Strange_Op_L, Strange_Op_R,
ActionCG, ActionCG,
ActionCGL, ActionCGR, ActionCGL, ActionCGR,
DerivativeCGL, DerivativeCGR, DerivativeCGL, DerivativeCGR,
OFRp, true); 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);
@ -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));

1
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();