diff --git a/RbcUkqcd.hpp b/RbcUkqcd.hpp index 1ed3c50..95c05e9 100644 --- a/RbcUkqcd.hpp +++ b/RbcUkqcd.hpp @@ -39,25 +39,104 @@ struct RbcUkqcd unsigned int L, T, Ls; }; - static constexpr EnsembleParameters m0UnitaryPar{0.000678, 0.02661, 1.8, 2., 64, 128, 12}; - static constexpr EnsembleParameters m0LCDPar{0.0006203, 0.02661, 1.8, 2., 64, 128, 12}; + // M0 + inline static constexpr EnsembleParameters m0UnitaryPar{0.000678, 0.02661, 1.8, 2., 64, 128, 12}; + inline static constexpr EnsembleParameters m0LCDPar{0.0006203, 0.02661, 1.8, 2., 64, 128, 12}; + + // C1M + inline static constexpr EnsembleParameters c1mIRLPar{0.005, 0.0362, 1.8, 2., 24, 64, 24}; + inline static constexpr EnsembleParameters c1m16IRLPar{0.005, 0.0362, 1.8, 2., 16, 64, 24}; + inline static constexpr EnsembleParameters c1m20IRLPar{0.005, 0.0362, 1.8, 2., 20, 64, 24}; + inline static constexpr EnsembleParameters c1m32IRLPar{0.005, 0.0362, 1.8, 2., 32, 64, 24}; + + // Runtime deflation parameters + struct DeflationParameters + { + double alpha, beta; + unsigned int nPoly, nStop, nK, nM; + }; + + inline static constexpr DeflationParameters c1mDeflPar{5.0e-04, 5.5, 101, 100, 110, 120}; + inline static constexpr DeflationParameters c1m16DeflPar{5.0e-03, 5.5, 101, 200, 220, 230}; + inline static constexpr DeflationParameters c1m20DeflPar{3.0e-03, 5.5, 101, 200, 220, 230}; + inline static constexpr DeflationParameters c1m32DeflPar{3.0e-06, 5.5, 101, 100, 110, 120}; + + // Light solver: load deflation from disk + static inline void + addM0LightLCDSolver(Application &app, const std::string solverName, const std::string gaugeName, + const std::string gaugeTransform, const std::string eigenpackPath, + const std::string boundary = "1 1 1 1", const double residual = 1.0e-8); + + // Light solver: deflation at runtime + static inline void addLightRuntimeIRLSolver(Application &app, + const RbcUkqcd::EnsembleParameters &par, + const RbcUkqcd::DeflationParameters &deflPar, + const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary, const double residual); + + static inline void addC1MLightRuntimeIRLSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary = "1 1 1 -1", + const double residual = 1.0e-8); + static inline void addC1M16LightRuntimeIRLSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary = "1 1 1 -1", + const double residual = 1.0e-8); + static inline void addC1M20LightRuntimeIRLSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary = "1 1 1 -1", + const double residual = 1.0e-8); + static inline void addC1M32LightRuntimeIRLSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary = "1 1 1 -1", + const double residual = 1.0e-8); + + // Strange solver (undeflated) + static inline void addStrangeSolver(Application &app, const RbcUkqcd::EnsembleParameters &par, + const std::string solverName, const std::string gaugeName, + const std::string gaugeTransform, const std::string boundary, + const double residual); - // Solvers (LCD: Local Coherence Deflation) - static inline void addM0LightLCDSolver(Application &app, const std::string solverName, - const std::string gaugeName, - const std::string gaugeTransform, - const std::string eigenpackPath, - const double residual = 1.0e-8); static inline void addM0StrangeSolver(Application &app, const std::string solverName, const std::string gaugeName, const std::string gaugeTransform, + const std::string boundary = "1 1 1 1", const double residual = 1.0e-8); + static inline void addC1MStrangeSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary = "1 1 1 -1", + const double residual = 1.0e-8); + static inline void addC1M16StrangeSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary = "1 1 1 -1", + const double residual = 1.0e-8); + static inline void addC1M20StrangeSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary = "1 1 1 -1", + const double residual = 1.0e-8); + static inline void addC1M32StrangeSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary = "1 1 1 -1", + const double residual = 1.0e-8); }; // Implementations ///////////////////////////////////////////////////////////////////////////////// + +// Light M0 (load deflation from disk) void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverName, const std::string gaugeName, const std::string gaugeTransform, - const std::string eigenpackPath, const double residual) + const std::string eigenpackPath, const std::string boundary, + const double residual) { const std::string prefix = solverName; @@ -75,7 +154,7 @@ void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverNam actionPar.M5 = RbcUkqcd::m0LCDPar.M5; actionPar.mass = RbcUkqcd::m0LCDPar.ml; actionPar.scale = RbcUkqcd::m0LCDPar.scale; - actionPar.boundary = "1 1 1 1"; + actionPar.boundary = boundary; actionPar.twist = "0. 0. 0. 0."; app.createModule(prefix + "_dwf", actionPar); actionPar.gauge = prefix + "_gauge_fp32"; @@ -117,9 +196,12 @@ void RbcUkqcd::addM0LightLCDSolver(Application &app, const std::string solverNam app.createModule(solverName, solverPar); } -void RbcUkqcd::addM0StrangeSolver(Application &app, const std::string solverName, - const std::string gaugeName, const std::string gaugeTransform, - const double residual) +// Deflation at runtime +void RbcUkqcd::addLightRuntimeIRLSolver(Application &app, const RbcUkqcd::EnsembleParameters &par, + const RbcUkqcd::DeflationParameters &deflPar, + const std::string solverName, const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary, const double residual) { const std::string prefix = solverName; @@ -133,17 +215,125 @@ void RbcUkqcd::addM0StrangeSolver(Application &app, const std::string solverName MAction::ScaledDWF::Par actionPar; actionPar.gauge = gaugeName; - actionPar.Ls = RbcUkqcd::m0LCDPar.Ls; - actionPar.M5 = RbcUkqcd::m0LCDPar.M5; - actionPar.mass = RbcUkqcd::m0LCDPar.ms; - actionPar.scale = RbcUkqcd::m0LCDPar.scale; - actionPar.boundary = "1 1 1 1"; + actionPar.Ls = par.Ls; + actionPar.M5 = par.M5; + actionPar.mass = par.ml; + actionPar.scale = par.scale; + actionPar.boundary = boundary; actionPar.twist = "0. 0. 0. 0."; app.createModule(prefix + "_dwf", actionPar); actionPar.gauge = prefix + "_gauge_fp32"; app.createModule(prefix + "_dwf_fp32", actionPar); - // Batched mixed-precision red-black preconditionned CG + // Guesser + MFermion::Operators::Par opPar; + opPar.action = prefix + "_dwf"; + app.createModule(prefix + "_dwf_op", opPar); + + MSolver::FermionImplicitlyRestartedLanczosIo32::Par lanPar; + lanPar.op = prefix + "_dwf_op_schur"; + lanPar.multiFile = false; + lanPar.redBlack = true; + lanPar.lanczosParams.Cheby.alpha = deflPar.alpha; + lanPar.lanczosParams.Cheby.beta = deflPar.beta; + lanPar.lanczosParams.Cheby.Npoly = deflPar.nPoly; + lanPar.lanczosParams.Nstop = deflPar.nStop; + lanPar.lanczosParams.Nk = deflPar.nK; + lanPar.lanczosParams.Nm = deflPar.nM; + lanPar.lanczosParams.resid = 3e-6; + lanPar.lanczosParams.MaxIt = 10000; + lanPar.lanczosParams.betastp = 0; + lanPar.lanczosParams.MinRes = 0; + lanPar.output = ""; + app.createModule(prefix + "_epack", lanPar); + + MGuesser::ExactDeflation::Par guessPar; + guessPar.eigenPack = prefix + "_epack"; + guessPar.size = deflPar.nStop; + app.createModule(prefix + "_defl", guessPar); + + // Mixed-precision red-black preconditionned CG + MSolver::MixedPrecisionRBPrecCG::Par solverPar; + + solverPar.innerAction = prefix + "_dwf_fp32"; + solverPar.outerAction = prefix + "_dwf"; + solverPar.maxInnerIteration = 30000; + solverPar.maxOuterIteration = 100; + solverPar.residual = residual; + solverPar.innerGuesser = ""; + solverPar.outerGuesser = prefix + "_defl"; + app.createModule(solverName, solverPar); +} + +// Light C1M +void RbcUkqcd::addC1MLightRuntimeIRLSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary, const double residual) +{ + RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1mIRLPar, RbcUkqcd::c1mDeflPar, solverName, + gaugeName, gaugeTransform, boundary, residual); +} + +// Light C1M16 +void RbcUkqcd::addC1M16LightRuntimeIRLSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary, const double residual) +{ + RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m16IRLPar, RbcUkqcd::c1m16DeflPar, solverName, + gaugeName, gaugeTransform, boundary, residual); +} + +// Light C1M20 +void RbcUkqcd::addC1M20LightRuntimeIRLSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary, const double residual) +{ + RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m20IRLPar, RbcUkqcd::c1m20DeflPar, solverName, + gaugeName, gaugeTransform, boundary, residual); +} + +// Light C1M32 +void RbcUkqcd::addC1M32LightRuntimeIRLSolver(Application &app, const std::string solverName, + const std::string gaugeName, + const std::string gaugeTransform, + const std::string boundary, const double residual) +{ + RbcUkqcd::addLightRuntimeIRLSolver(app, RbcUkqcd::c1m32IRLPar, RbcUkqcd::c1m32DeflPar, solverName, + gaugeName, gaugeTransform, boundary, residual); +} + +// Strange +void RbcUkqcd::addStrangeSolver(Application &app, const RbcUkqcd::EnsembleParameters &par, + const std::string solverName, const std::string gaugeName, + const std::string gaugeTransform, const std::string boundary, + const double residual) +{ + const std::string prefix = solverName; + + // Gauge field FP32 cast + MUtilities::GaugeSinglePrecisionCast::Par gaugeCastPar; + + gaugeCastPar.field = gaugeName; + app.createModule(prefix + "_gauge_fp32", gaugeCastPar); + + // Scaled DWF action + FP32 version + MAction::ScaledDWF::Par actionPar; + + actionPar.gauge = gaugeName; + actionPar.Ls = par.Ls; + actionPar.M5 = par.M5; + actionPar.mass = par.ms; + actionPar.scale = par.scale; + actionPar.boundary = boundary; + actionPar.twist = "0. 0. 0. 0."; + app.createModule(prefix + "_dwf", actionPar); + actionPar.gauge = prefix + "_gauge_fp32"; + app.createModule(prefix + "_dwf_fp32", actionPar); + + // Mixed-precision red-black preconditionned CG MSolver::MixedPrecisionRBPrecCG::Par solverPar; solverPar.innerAction = prefix + "_dwf_fp32"; @@ -156,4 +346,44 @@ void RbcUkqcd::addM0StrangeSolver(Application &app, const std::string solverName app.createModule(solverName, solverPar); } +void RbcUkqcd::addM0StrangeSolver(Application &app, const std::string solverName, + const std::string gaugeName, const std::string gaugeTransform, + const std::string boundary, const double residual) +{ + RbcUkqcd::addStrangeSolver(app, RbcUkqcd::m0LCDPar, solverName, gaugeName, gaugeTransform, + boundary, residual); +} + +void RbcUkqcd::addC1MStrangeSolver(Application &app, const std::string solverName, + const std::string gaugeName, const std::string gaugeTransform, + const std::string boundary, const double residual) +{ + RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1mIRLPar, solverName, gaugeName, gaugeTransform, + boundary, residual); +} + +void RbcUkqcd::addC1M16StrangeSolver(Application &app, const std::string solverName, + const std::string gaugeName, const std::string gaugeTransform, + const std::string boundary, const double residual) +{ + RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m16IRLPar, solverName, gaugeName, gaugeTransform, + boundary, residual); +} + +void RbcUkqcd::addC1M20StrangeSolver(Application &app, const std::string solverName, + const std::string gaugeName, const std::string gaugeTransform, + const std::string boundary, const double residual) +{ + RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m20IRLPar, solverName, gaugeName, gaugeTransform, + boundary, residual); +} + +void RbcUkqcd::addC1M32StrangeSolver(Application &app, const std::string solverName, + const std::string gaugeName, const std::string gaugeTransform, + const std::string boundary, const double residual) +{ + RbcUkqcd::addStrangeSolver(app, RbcUkqcd::c1m32IRLPar, solverName, gaugeName, gaugeTransform, + boundary, residual); +} + } // namespace hadpresets