/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid Source file: Copyright (C) 2015-2016 Author: Peter Boyle Author: Guido Cossu Author: David Murphy This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ #include #ifdef GRID_DEFAULT_PRECISION_DOUBLE #define MIXED_PRECISION #endif /* * 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. */ #include int main(int argc, char **argv) { using namespace Grid; Grid_init(&argc, &argv); int threads = GridThread::GetThreads(); // here make a routine to print all the relevant information on the run std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; // Typedefs to simplify notation typedef WilsonImplR FermionImplPolicy; typedef MobiusFermionR FermionAction; typedef MobiusFermionF FermionActionF; typedef MobiusEOFAFermionR FermionEOFAAction; typedef MobiusEOFAFermionF FermionEOFAActionF; typedef typename FermionAction::FermionField FermionField; typedef typename FermionActionF::FermionField FermionFieldF; typedef Grid::XmlReader Serialiser; //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: // typedef GenericHMCRunner HMCWrapper; // typedef GenericHMCRunner HMCWrapper; typedef GenericHMCRunner HMCWrapper; HMCparameters HMCparams; { XmlReader HMCrd("HMCparameters.xml"); read(HMCrd,"HMCparameters",HMCparams); std::cout << GridLogMessage<< HMCparams < PlaqObs; TheHMC.Resources.AddObservable(); ////////////////////////////////////////////// const int Ls = 12; Real beta = 2.31; Real light_mass = 0.002144; Real strange_mass = 0.02144; Real pv_mass = 1.0; RealD M5 = 1.8; RealD b = 1.5; RealD c = 0.5; // Copied from paper std::vector hasenbusch({ 0.005, 0.017, 0.07, 0.18, 0.45 }); // Paper values from F1 incorrect run // std::vector hasenbusch({ 0.004, 0.016, 0.07, 0.18, 0.45 }); // Paper values from F1 incorrect run 120-130 // std::vector hasenbusch({ 0.004, 0.015, 0.07, 0.18, 0.45 }); // Paper values from F1 incorrect run 110-120 // std::vector hasenbusch({ 0.005, 0.017, 0.07, 0.18, 0.45 }); // Paper values from F1 incorrect run // std::vector hasenbusch({ 0.01, 0.07, 0.18, 0.45 }); // Experiment? /////////////////////////////////////////////////////////////////////////////////////////////// //Bad choices with large dH. Equalising force L2 norm was not wise. /////////////////////////////////////////////////////////////////////////////////////////////// //std::vector hasenbusch({ 0.03, 0.2, 0.3, 0.5, 0.8 }); //std::vector hasenbusch({ 0.05, 0.2, 0.4, 0.6, 0.8 }); auto GridPtr = TheHMC.Resources.GetCartesian(); auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr); auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr); Coordinate latt = GridDefaultLatt(); Coordinate mpi = GridDefaultMpi(); Coordinate simdF = GridDefaultSimd(Nd,vComplexF::Nsimd()); Coordinate simdD = GridDefaultSimd(Nd,vComplexD::Nsimd()); auto GridPtrF = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi); auto GridRBPtrF = SpaceTimeGrid::makeFourDimRedBlackGrid(GridPtrF); auto FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtrF); auto FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtrF); IwasakiGaugeActionR GaugeAction(beta); // temporarily need a gauge field LatticeGaugeField U(GridPtr); LatticeGaugeFieldF UF(GridPtrF); // These lines are unecessary if BC are all periodic std::vector boundary = {1,1,1,-1}; FermionAction::ImplParams Params(boundary); FermionActionF::ImplParams ParamsF(boundary); double ActionStoppingCondition = 1e-10; double DerivativeStoppingCondition = 1e-7; double MaxCGIterations = 30000; //////////////////////////////////// // Collect actions //////////////////////////////////// ActionLevel Level1(1); ActionLevel Level2(16); //////////////////////////////////// // Strange action //////////////////////////////////// typedef SchurDiagMooeeOperator LinearOperatorF; typedef SchurDiagMooeeOperator LinearOperatorD; typedef SchurDiagMooeeOperator LinearOperatorEOFAF; typedef SchurDiagMooeeOperator LinearOperatorEOFAD; typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG; typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG_EOFA; // DJM: setup for EOFA ratio (Mobius) OneFlavourRationalParams OFRp; OFRp.lo = 0.1; // How do I know this on F1? OFRp.hi = 25.0; OFRp.MaxIter = 10000; OFRp.tolerance= 1.0e-9; OFRp.degree = 12; OFRp.precision= 50; 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); MobiusEOFAFermionR Strange_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c); MobiusEOFAFermionF Strange_Op_RF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c); ConjugateGradient ActionCG(ActionStoppingCondition,MaxCGIterations); ConjugateGradient DerivativeCG(DerivativeStoppingCondition,MaxCGIterations); #ifdef MIXED_PRECISION const int MX_inner = 5000; const RealD MX_tol = 1.0e-6; // Mixed precision EOFA LinearOperatorEOFAD Strange_LinOp_L (Strange_Op_L); LinearOperatorEOFAD Strange_LinOp_R (Strange_Op_R); LinearOperatorEOFAF Strange_LinOp_LF(Strange_Op_LF); LinearOperatorEOFAF Strange_LinOp_RF(Strange_Op_RF); MxPCG_EOFA ActionCGL(ActionStoppingCondition, MX_tol, MX_inner, MaxCGIterations, FrbGridF, Strange_Op_LF,Strange_Op_L, Strange_LinOp_LF,Strange_LinOp_L); MxPCG_EOFA DerivativeCGL(DerivativeStoppingCondition, MX_tol, MX_inner, MaxCGIterations, FrbGridF, Strange_Op_LF,Strange_Op_L, Strange_LinOp_LF,Strange_LinOp_L); MxPCG_EOFA ActionCGR(ActionStoppingCondition, MX_tol, MX_inner, MaxCGIterations, FrbGridF, Strange_Op_RF,Strange_Op_R, Strange_LinOp_RF,Strange_LinOp_R); MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition, MX_tol, MX_inner, MaxCGIterations, FrbGridF, Strange_Op_RF,Strange_Op_R, Strange_LinOp_RF,Strange_LinOp_R); ExactOneFlavourRatioPseudoFermionAction EOFA(Strange_Op_L, Strange_Op_R, ActionCG, ActionCGL, ActionCGR, DerivativeCGL, DerivativeCGR, OFRp, true); #else ExactOneFlavourRatioPseudoFermionAction EOFA(Strange_Op_L, Strange_Op_R, ActionCG, ActionCG, ActionCG, ActionCG, ActionCG, // DerivativeCG, DerivativeCG, OFRp, true); #endif Level1.push_back(&EOFA); //////////////////////////////////// // up down action //////////////////////////////////// std::vector light_den; std::vector light_num; int n_hasenbusch = hasenbusch.size(); light_den.push_back(light_mass); for(int h=0;h Numerators; std::vector Denominators; std::vector *> Quotients; std::vector ActionMPCG; std::vector MPCG; std::vector DenominatorsF; std::vector LinOpD; std::vector LinOpF; for(int h=0;h(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],ActionCG)); #else //////////////////////////////////////////////////////////////////////////// // Standard CG for 2f force //////////////////////////////////////////////////////////////////////////// Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction(*Numerators[h],*Denominators[h],DerivativeCG,ActionCG)); #endif } for(int h=0;h