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

Mixed prec update

Browse Source
This commit is contained in:
Quadro 2021-06-09 16:33:02 -04:00
parent 096bb98f78
commit 4f5ad73501
1 changed files with 8 additions and 142 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

150
HMC/Mobius2p1fEOFA.cc
View File

@ -33,137 +33,8 @@ directory
#ifdef GRID_DEFAULT_PRECISION_DOUBLE #ifdef GRID_DEFAULT_PRECISION_DOUBLE
#define MIXED_PRECISION #define MIXED_PRECISION
#endif #endif
#include <Grid/qcd/utils/MixedPrecisionOperatorFunction.h>
NAMESPACE_BEGIN(Grid);
/*
* Need a plan for gauge field update for mixed precision in HMC (2x speed up)
* -- Store the single prec action operator.
* -- Clone the gauge field from the operator function argument.
* -- Build the mixed precision operator dynamically from the passed operator and single prec clone.
*/
template<class FermionOperatorD, class FermionOperatorF, class SchurOperatorD, class SchurOperatorF>
class MixedPrecisionConjugateGradientOperatorFunction : public OperatorFunction<typename FermionOperatorD::FermionField> {
public:
typedef typename FermionOperatorD::FermionField FieldD;
typedef typename FermionOperatorF::FermionField FieldF;
using OperatorFunction<FieldD>::operator();
RealD Tolerance;
RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed
Integer MaxInnerIterations;
Integer MaxOuterIterations;
GridBase* SinglePrecGrid4; //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
FermionOperatorF &FermOpF;
FermionOperatorD &FermOpD;;
SchurOperatorF &LinOpF;
SchurOperatorD &LinOpD;
Integer TotalInnerIterations; //Number of inner CG iterations
Integer TotalOuterIterations; //Number of restarts
Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step
MixedPrecisionConjugateGradientOperatorFunction(RealD tol,
Integer maxinnerit,
Integer maxouterit,
GridBase* _sp_grid4,
GridBase* _sp_grid5,
FermionOperatorF &_FermOpF,
FermionOperatorD &_FermOpD,
SchurOperatorF &_LinOpF,
SchurOperatorD &_LinOpD):
LinOpF(_LinOpF),
LinOpD(_LinOpD),
FermOpF(_FermOpF),
FermOpD(_FermOpD),
Tolerance(tol),
InnerTolerance(tol),
MaxInnerIterations(maxinnerit),
MaxOuterIterations(maxouterit),
SinglePrecGrid4(_sp_grid4),
SinglePrecGrid5(_sp_grid5),
OuterLoopNormMult(100.)
{
/* 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 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 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);
}
};
NAMESPACE_END(Grid);
int main(int argc, char **argv) { int main(int argc, char **argv) {
using namespace Grid; using namespace Grid;
@ -302,40 +173,37 @@ int main(int argc, char **argv) {
ConjugateGradient<FermionField> DerivativeCG(DerivativeStoppingCondition,MaxCGIterations); ConjugateGradient<FermionField> DerivativeCG(DerivativeStoppingCondition,MaxCGIterations);
#ifdef MIXED_PRECISION #ifdef MIXED_PRECISION
const int MX_inner = 1000; const int MX_inner = 1000;
const RealD MX_tol = 1.0e-6;
// Mixed precision EOFA // Mixed precision EOFA
LinearOperatorEOFAD Strange_LinOp_L (Strange_Op_L); LinearOperatorEOFAD Strange_LinOp_L (Strange_Op_L);
LinearOperatorEOFAD Strange_LinOp_R (Strange_Op_R); LinearOperatorEOFAD Strange_LinOp_R (Strange_Op_R);
LinearOperatorEOFAF Strange_LinOp_LF(Strange_Op_LF); LinearOperatorEOFAF Strange_LinOp_LF(Strange_Op_LF);
LinearOperatorEOFAF Strange_LinOp_RF(Strange_Op_RF); LinearOperatorEOFAF Strange_LinOp_RF(Strange_Op_RF);
MxPCG_EOFA ActionCGL(ActionStoppingCondition, MxPCG_EOFA ActionCGL(ActionStoppingCondition,MX_tol,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
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_tol,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
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_tol,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
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_tol,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
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);
@ -401,18 +269,16 @@ int main(int argc, char **argv) {
LinOpD.push_back(new LinearOperatorD(*Denominators[h])); LinOpD.push_back(new LinearOperatorD(*Denominators[h]));
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_tol,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
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_tol,
MX_inner, MX_inner,
MaxCGIterations, MaxCGIterations,
GridPtrF,
FrbGridF, FrbGridF,
*DenominatorsF[h],*Denominators[h], *DenominatorsF[h],*Denominators[h],
*LinOpF[h], *LinOpD[h]) ); *LinOpF[h], *LinOpD[h]) );