From fd8b6a23a6b7d1d7a9c12ef8ca5129cccafd980c Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 12 Aug 2026 12:49:21 -0400 Subject: [PATCH] PvDagM or other left prec precon two flavour ratio --- .../pseudofermion/TwoFlavourRatioLeftPrec.h | 206 ++++++++++++++++++ 1 file changed, 206 insertions(+) create mode 100644 Grid/qcd/action/pseudofermion/TwoFlavourRatioLeftPrec.h diff --git a/Grid/qcd/action/pseudofermion/TwoFlavourRatioLeftPrec.h b/Grid/qcd/action/pseudofermion/TwoFlavourRatioLeftPrec.h new file mode 100644 index 000000000..687d784d9 --- /dev/null +++ b/Grid/qcd/action/pseudofermion/TwoFlavourRatioLeftPrec.h @@ -0,0 +1,206 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/pseudofermion/TwoFlavourRatioLeftPrec.h + + Copyright (C) 2026 + +Author: Peter Boyle + + 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 */ +#pragma once + +NAMESPACE_BEGIN(Grid); + +/////////////////////////////////////////////////////////////////////////////// +// Two flavour ratio with LEFT-PRECONDITIONED solves. +// +// Same action content as TwoFlavourRatio.h: +// +// S = phi^dag V (Mdag M)^-1 Vdag phi ==> det[ Mdag M / Vdag V ] +// +// (V = NumOp the heavier / Pauli-Villars operator, M = DenOp the lighter), +// but organised around the composite +// +// F = Vdag M +// +// which is the 2-hop-coarsenable operator the non-Hermitian multigrid +// serves. Solving M X = b as F X = Vdag b is LEFT PRECONDITIONING by +// Vdag; the determinant/action layer is the standard quotient, and all +// novelty is confined to the solver contract. +// +// TwoFlavourRatio.h is tied to a normal-equations solver: one (MdagM)^-1 +// solve, then Y = M X gives Mdag^-1 Vdag phi almost free. The left- +// preconditioned idiom is DIFFERENT: the chain +// +// b = Vdag phi +// z : Fdag z = b (adjoint F solve) +// Y = V z (= Mdag^-1 Vdag phi -- harvested from solve 1) +// s = Vdag Y (= Vdag V z) +// X : F X = s (forward F solve; X = (MdagM)^-1 Vdag phi) +// +// yields Y BEFORE X (so S(U) needs only the adjoint solve), with Y's +// accuracy independent of the second solve. Force terms are then the +// standard four MDeriv insertions of TwoFlavourRatio. +// +// Solver slots are LinearFunctions with the F-SOLVE contract (solution +// overwritten, zero guess imposed internally): +// ForwardSolver(b,x) : F x = b +// AdjointSolver(b,z) : Fdag z = b +// implemented in production by the multigrid-GCR stack (forward cycle and +// adjoint cycle); in tests by CG on the composite normal equations. +// HeatbathSolver(b,x) : x = (Vdag V)^-1 b -- heavy operator, plain CG. +// +// Heatbath is exact by operator algebra: phi = V (VdagV)^-1 Mdag eta +// ==> S = | Mdag^-1 Vdag phi |^2 = |eta|^2 (to solver tolerance); the +// deterministic refresh(U,eta) hook below is the test point. +// +// Hasenbusch: nothing requires V to have mass one; any (heavier,lighter) +// pair works, F(V,M) = Vdag M coarsenable by the same machinery, rungs' +// solves are F-family (mrhs-batchable, mass-shared coarse space). +/////////////////////////////////////////////////////////////////////////////// +template +class TwoFlavourRatioLeftPrecPseudoFermionAction : public Action { +public: + INHERIT_IMPL_TYPES(Impl); + +private: + FermionOperator & NumOp;// V + FermionOperator & DenOp;// M + + LinearFunction &DerivForwardSolver; // F x = b, MD tolerance + LinearFunction &DerivAdjointSolver; // Fdag z = b, MD tolerance + LinearFunction &ActionAdjointSolver; // Fdag z = b, accept/reject tolerance + LinearFunction &HeatbathSolver; // (VdagV)^-1 b, heavy op + + FermionField Phi; // the pseudo fermion field for this trajectory + +public: + TwoFlavourRatioLeftPrecPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + LinearFunction & DFS, + LinearFunction & DAS, + LinearFunction & AAS, + LinearFunction & HS + ) : NumOp(_NumOp), + DenOp(_DenOp), + DerivForwardSolver(DFS), + DerivAdjointSolver(DAS), + ActionAdjointSolver(AAS), + HeatbathSolver(HS), + Phi(_NumOp.FermionGrid()) + {}; + + virtual std::string action_name(){return "TwoFlavourRatioLeftPrecPseudoFermionAction";} + + virtual std::string LogParameters(){ + std::stringstream sstream; + sstream << GridLogMessage << "["< sig^2 = 0.5 ; eta enters with width 1/sqrt(2). + RealD scale = std::sqrt(0.5); + FermionField eta(NumOp.FermionGrid()); + gaussian(pRNG,eta); + eta = eta * scale; + refresh(U,eta); + } + + // Deterministic-noise variant (test hook): + // after this, S(U) == norm2(eta) exactly (to solver tolerance). + void refresh(const GaugeField &U, const FermionField &eta) { + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + FermionField tmp(NumOp.FermionGrid()); + FermionField w (NumOp.FermionGrid()); + + DenOp.Mdag(eta,tmp); // tmp = Mdag eta + w = Zero(); + HeatbathSolver(tmp,w); // w = (VdagV)^-1 Mdag eta + NumOp.M(w,Phi); // Phi = V (VdagV)^-1 Mdag eta = Vdag^-1 Mdag eta + std::cout << GridLogMessage << action_name() << " refresh |Phi|^2 = "<< norm2(Phi)< X = (MdagM)^-1 Vdag phi + + // phi^dag V (MdagM)^-1 dVdag phi + NumOp.MDeriv(force , X, Phi, DaggerYes); dSdU = force; + // phi^dag dV (MdagM)^-1 Vdag phi + NumOp.MDeriv(force , Phi, X, DaggerNo ); dSdU = dSdU+force; + // - phi^dag V (MdagM)^-1 Mdag dM (MdagM)^-1 Vdag phi + // - phi^dag V (MdagM)^-1 dMdag M (MdagM)^-1 Vdag phi + DenOp.MDeriv(force, Y, X, DaggerNo ); dSdU = dSdU-force; + DenOp.MDeriv(force, X, Y, DaggerYes); dSdU = dSdU-force; + + dSdU *= -1.0; + }; +}; + +NAMESPACE_END(Grid);