From 09b233b82e5c302aaf543d86579b3f37a987d651 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 14 May 2021 22:21:23 -0400 Subject: [PATCH] RdagR inverse --- .../DomainDecomposedTwoFlavourBoundary.h | 131 ++++++------------ 1 file changed, 41 insertions(+), 90 deletions(-) diff --git a/Grid/qcd/action/pseudofermion/DomainDecomposedTwoFlavourBoundary.h b/Grid/qcd/action/pseudofermion/DomainDecomposedTwoFlavourBoundary.h index 7eee1229..24a7ffd8 100644 --- a/Grid/qcd/action/pseudofermion/DomainDecomposedTwoFlavourBoundary.h +++ b/Grid/qcd/action/pseudofermion/DomainDecomposedTwoFlavourBoundary.h @@ -38,7 +38,6 @@ public: INHERIT_IMPL_TYPES(Impl); private: - SchurFactoredFermionOperator & NumOp;// the basic operator SchurFactoredFermionOperator & DenOp;// the basic operator OperatorFunction &DerivativeSolver; @@ -46,16 +45,16 @@ private: FermionField Phi; // the pseudo fermion field for this trajectory + RealD refresh_action; public: - DomainBoundaryPseudoFermionAction(SchurFactoredFermionOperator &_NumOp, - SchurFactoredFermionOperator &_DenOp, + DomainBoundaryPseudoFermionAction(SchurFactoredFermionOperator &_DenOp, OperatorFunction & DS, OperatorFunction & AS - ) : NumOp(_NumOp), DenOp(_DenOp), + ) : DenOp(_DenOp), DerivativeSolver(DS), ActionSolver(AS), - Phi(_NumOp.FermionGrid()) {}; + Phi(_DenOp.FermOp.FermionGrid()) {}; - virtual std::string action_name(){return "DomainBoundaryPseudoFermionRatioAction";} + virtual std::string action_name(){return "DomainBoundaryPseudoFermionAction";} virtual std::string LogParameters(){ @@ -63,14 +62,13 @@ public: return sstream.str(); } - virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) + virtual void refresh(const GaugeField &U, GridSerialRNG& sRNG, GridParallelRNG& pRNG) { - // P(phi) = e^{- phi^dag P^dag Rdag^-1 R^-1 P phi} + // P(phi) = e^{- phi^dag Rdag^-1 R^-1 phi} // - // NumOp == P // DenOp == R // - // Take phi = P^{-1} R eta ; eta = R^-1 P Phi + // Take phi = R eta ; eta = R^-1 Phi // // P(eta) = e^{- eta^dag eta} // @@ -80,55 +78,40 @@ public: // RealD scale = std::sqrt(0.5); - NumOp.ImportGauge(U); DenOp.ImportGauge(U); - FermionField eta(NumOp.FermOp.FermionGrid()); - FermionField tmp(NumOp.FermOp.FermionGrid()); + FermionField eta(DenOp.FermOp.FermionGrid()); - gaussian(pRNG,eta); + gaussian(pRNG,eta); eta=eta*scale; + + DenOp.ProjectBoundaryBar(eta); + DenOp.R(eta,Phi); - NumOp.ProjectBoundaryBar(eta); - - NumOp.R(eta,tmp); - DenOp.RInv(tmp,Phi); - - Phi=Phi*scale; - - NumOp.ProjectBoundaryBar(Phi); + refresh_action = norm2(eta); }; ////////////////////////////////////////////////////// - // S = phi^dag V (Mdag M)^-1 Vdag phi + // S = phi^dag Rdag^-1 R^-1 phi ////////////////////////////////////////////////////// virtual RealD S(const GaugeField &U) { - NumOp.ImportGauge(U); DenOp.ImportGauge(U); - FermionField X(NumOp.FermOp.FermionGrid()); - FermionField Y(NumOp.FermOp.FermionGrid()); + FermionField X(DenOp.FermOp.FermionGrid()); - NumOp.R(Phi,Y); - DenOp.RInv(Y,X); + DenOp.RInv(Phi,X); RealD action = norm2(X); return action; }; - ////////////////////////////////////////////////////// - // dS/du = phi^dag dV (Mdag M)^-1 V^dag phi - // - phi^dag V (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 V^dag phi - // + phi^dag V (Mdag M)^-1 dV^dag phi - ////////////////////////////////////////////////////// - virtual void deriv(const GaugeField &U,GaugeField & dSdU) { - - NumOp.ImportGauge(U); + virtual void deriv(const GaugeField &U,GaugeField & dSdU) + { DenOp.ImportGauge(U); - GridBase *fgrid = NumOp.FermOp.FermionGrid(); - GridBase *ugrid = NumOp.FermOp.GaugeGrid(); + GridBase *fgrid = DenOp.FermOp.FermionGrid(); + GridBase *ugrid = DenOp.FermOp.GaugeGrid(); FermionField X(fgrid); FermionField Y(fgrid); @@ -136,66 +119,34 @@ public: GaugeField force(ugrid); - FermionField DobiDdbPhi(fgrid); // Vector A in my notes - FermionField DoiDdDobiDdbPhi(fgrid); // Vector B in my notes - FermionField DiDdbP_Phi(fgrid); // Vector C in my notes - FermionField DidRinvP_Phi(fgrid); // Vector D in my notes - FermionField DoidRinvDagRinvP_Phi(fgrid); // Vector E in my notes - FermionField DobidDddDoidRinvDagRinvP_Phi(fgrid); // Vector F in my notes + FermionField DiDdb_Phi(fgrid); // Vector C in my notes + FermionField DidRinv_Phi(fgrid); // Vector D in my notes + FermionField DdbdDidRinv_Phi(fgrid); - FermionField P_Phi(fgrid); - FermionField RinvP_Phi(fgrid); - FermionField RinvDagRinvP_Phi(fgrid); + FermionField Rinv_Phi(fgrid); + FermionField RinvDagRinv_Phi(fgrid); - // P term - NumOp.dBoundaryBar(Phi,tmp); - NumOp.dOmegaBarInv(tmp,DobiDdbPhi); // Vector A - NumOp.dBoundary(DobiDdbPhi,tmp); - NumOp.dOmegaInv(tmp,DoiDdDobiDdbPhi); // Vector B - P_Phi = Phi - DoiDdDobiDdbPhi; - NumOp.ProjectBoundaryBar(P_Phi); + // R^-1 term + DenOp.dBoundaryBar(Phi,tmp); + DenOp.Dinverse(tmp,DiDdb_Phi); // Vector C + Rinv_Phi = Phi - DiDdb_Phi; + DenOp.ProjectBoundaryBar(Rinv_Phi); + + // R^-dagger R^-1 term + DenOp.DinverseDag(Rinv_Phi,DidRinv_Phi); // Vector D + DenOp.dBoundaryBarDag(DidRinv_Phi,DdbdDidRinv_Phi); + RinvDagRinv_Phi = Rinv_Phi - DdbdDidRinv_Phi; + DenOp.ProjectBoundaryBar(RinvDagRinv_Phi); - // R^-1 P term - DenOp.dBoundaryBar(P_Phi,tmp); - DenOp.Dinverse(tmp,DiDdbP_Phi); // Vector C - RinvP_Phi = P_Phi - DiDdbP_Phi; - DenOp.ProjectBoundaryBar(RinvP_Phi); - - // R^-dagger R^-1 P term - DenOp.DinverseDag(RinvP_Phi,DidRinvP_Phi); // Vector D - RinvDagRinvP_Phi = RinvP_Phi - DidRinvP_Phi; - DenOp.ProjectBoundaryBar(RinvDagRinvP_Phi); - - // P^dag R^-dagger R^-1 P term - NumOp.dOmegaDagInv(RinvDagRinvP_Phi,DoidRinvDagRinvP_Phi); // Vector E - NumOp.dBoundaryDag(DoidRinvDagRinvP_Phi,tmp); - NumOp.dOmegaBarDagInv(tmp,DobidDddDoidRinvDagRinvP_Phi); // Vector F - - - dSdU=Zero(); - - // phi^dag V (Mdag M)^-1 dV^dag phi - X = DobiDdbPhi; - Y = DobidDddDoidRinvDagRinvP_Phi; - NumOp.DirichletFermOp.MDeriv(force,X,Y,DaggerYes); dSdU=dsdU+force; - NumOp.DirichletFermOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU+force; - - X = DoiDdDobiDdbPhi; - Y = DoidRinvDagRinvP_Phi; - NumOp.DirichletFermOp.MDeriv(force,X,Y,DaggerYes); dSdU=dsdU+force; - NumOp.DirichletFermOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU+force; - - X = DiDdbP_Phi; - Y = DidRinvP_Phi; - NumOp.DirichletFermOp.MDeriv(force,X,Y,DaggerYes); dSdU=dsdU+force; - NumOp.DirichletFermOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU+force; + X = DiDdb_Phi; + Y = DidRinv_Phi; + DenOp.FermOp.MDeriv(force,Y,X,DaggerNo); dSdU=force; + DenOp.FermOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU+force; dSdU *= -1.0; - //dSdU = - Ta(dSdU); }; }; NAMESPACE_END(Grid); -#endif