diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.h b/Grid/qcd/action/fermion/CayleyFermion5D.h index 77491def2..11019b4b2 100644 --- a/Grid/qcd/action/fermion/CayleyFermion5D.h +++ b/Grid/qcd/action/fermion/CayleyFermion5D.h @@ -81,6 +81,15 @@ public: mass_minus=_mass_minus; SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c); // Reset coeffs } ; + + // Everything SetMass needs, plus the boost, then rebuild the coefficient arrays in clone + void CloneCoefficientsInto(CayleyFermion5D &clone) + { + clone.mass_plus = mass_plus; + clone.mass_minus = mass_minus; + clone.qmu = qmu; + clone.SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c); + } void P(const FermionField &psi, FermionField &chi); void Pdag(const FermionField &psi, FermionField &chi); diff --git a/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h b/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h index 18287d8aa..aa130007f 100644 --- a/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h +++ b/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h @@ -189,6 +189,71 @@ public: void ImportGauge(const GaugeField& _Umu) override; + // Install already computed clover fields and inverses, and derive the checkerboarded parts + void ImportCloverTerms(const CloverDiagonalField& D, + const CloverTriangleField& T, + const CloverDiagonalField& Dinv, + const CloverTriangleField& Tinv) + { + Diagonal = D; + Triangle = T; + DiagonalInv = Dinv; + TriangleInv = Tinv; + pickCheckerboard(Even, DiagonalEven, Diagonal); + pickCheckerboard(Even, TriangleEven, Triangle); + pickCheckerboard(Odd, DiagonalOdd, Diagonal); + pickCheckerboard(Odd, TriangleOdd, Triangle); + pickCheckerboard(Even, DiagonalInvEven, DiagonalInv); + pickCheckerboard(Even, TriangleInvEven, TriangleInv); + pickCheckerboard(Odd, DiagonalInvOdd, DiagonalInv); + pickCheckerboard(Odd, TriangleInvOdd, TriangleInv); + } + + // Exact type only: derived operators must not inherit this + SplitOperator *SplitClone(const Coordinate& mpi_split) override + { + if ( typeid(*this) != typeid(CompactWilsonCloverFermion) ) { + return nullptr; + } + if ( Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field and clover fields are overwritten below. + // Boundary masks are rebuilt by the constructor from global coordinates. + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + CompactWilsonCloverFermion *clone = + new CompactWilsonCloverFermion(Uplaceholder, + *split->GaugeGrid, + *split->GaugeRBGrid, + this->mass, + csw_r, + csw_t, + cF, + this->anisotropyCoeff, + this->Params); + // The stored coefficients are post-scaled; the constructor would scale them again + clone->csw_r = csw_r; + clone->csw_t = csw_t; + + CloverDiagonalField sD(split->GaugeGrid); + CloverTriangleField sT(split->GaugeGrid); + CloverDiagonalField sDinv(split->GaugeGrid); + CloverTriangleField sTinv(split->GaugeGrid); + Grid_split(Diagonal,sD); + Grid_split(Triangle,sT); + Grid_split(DiagonalInv,sDinv); + Grid_split(TriangleInv,sTinv); + clone->ImportCloverTerms(sD,sT,sDinv,sTinv); + + this->SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + ///////////////////////////////////////////// // Helpers ///////////////////////////////////////////// diff --git a/Grid/qcd/action/fermion/DomainWallFermion.h b/Grid/qcd/action/fermion/DomainWallFermion.h index 680c203eb..38fbcdc36 100644 --- a/Grid/qcd/action/fermion/DomainWallFermion.h +++ b/Grid/qcd/action/fermion/DomainWallFermion.h @@ -132,6 +132,35 @@ public: Approx::zolotarev_free(zdata); } + // Exact type only: derived operators must not inherit this + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + if ( typeid(*this) != typeid(DomainWallFermion) ) { + return nullptr; + } + if ( this->Dirichlet || Impl::LsVectorised || Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field is overwritten below + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + DomainWallFermion *clone = new DomainWallFermion(Uplaceholder, + *split->FermionGrid, + *split->FermionRBGrid, + *split->GaugeGrid, + *split->GaugeRBGrid, + this->mass_plus, + this->M5, + this->Params); + this->CloneCoefficientsInto(*clone); + this->SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + }; NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/FermionOperator.h b/Grid/qcd/action/fermion/FermionOperator.h index a2abaa781..be12d46dc 100644 --- a/Grid/qcd/action/fermion/FermionOperator.h +++ b/Grid/qcd/action/fermion/FermionOperator.h @@ -59,6 +59,35 @@ public: virtual GridBase *GaugeGrid(void) =0; virtual GridBase *GaugeRedBlackGrid(void) =0; + //////////////////////////////////////////////////////////////// + // Split-communicator copies of this operator's gauge and fermion + // grids, returned in a new bundle with no Matrix or Linop yet. + // Partitions have MPI layout mpi_split. Collective. + //////////////////////////////////////////////////////////////// + SplitOperator *MakeSplitGrids(const Coordinate &mpi_split) + { + GridCartesian *U = dynamic_cast(GaugeGrid()); + GRID_ASSERT(U != nullptr); + GRID_ASSERT(mpi_split.size() == Nd); + + SplitOperator *split = new SplitOperator(); + + split->GaugeGrid = new GridCartesian(U->FullDimensions(),U->_simd_layout,mpi_split,*U); + split->GaugeRBGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(split->GaugeGrid); + split->Partition = GridSplitVectorIndex(U,split->GaugeGrid); + split->Partitions = U->ProcessorCount()/split->GaugeGrid->ProcessorCount(); + + if ( FermionGrid()->Nd() == Nd+1 ) { + int Ls = FermionGrid()->_fdimensions[0]; + split->FermionGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,split->GaugeGrid); + split->FermionRBGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,split->GaugeGrid); + } else { + split->FermionGrid = split->GaugeGrid; + split->FermionRBGrid = split->GaugeRBGrid; + } + return split; + } + // override multiply virtual void M (const FermionField &in, FermionField &out)=0; virtual void Mdag (const FermionField &in, FermionField &out)=0; @@ -88,6 +117,10 @@ public: virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; virtual void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; + // 1 if MoeDeriv/MeoDeriv return a force on a single 4D checkerboard (4D red-black), + // 0 if they return it on the full 4D gauge grid (a checkerboard that includes s) + virtual int CheckerboardedForce(void) { return 1; }; + virtual void Mdiag (const FermionField &in, FermionField &out) { Mooee(in,out);}; // Same as Mooee applied to both CB's virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac virtual void MdirAll(const FermionField &in, std::vector &out)=0; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac diff --git a/Grid/qcd/action/fermion/MobiusFermion.h b/Grid/qcd/action/fermion/MobiusFermion.h index cb3c85fd6..e1c3c5141 100644 --- a/Grid/qcd/action/fermion/MobiusFermion.h +++ b/Grid/qcd/action/fermion/MobiusFermion.h @@ -70,6 +70,37 @@ public: } + // Exact type only: derived operators (e.g. ScaledShamirFermion) must not inherit this + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + if ( typeid(*this) != typeid(MobiusFermion) ) { + return nullptr; + } + if ( this->Dirichlet || Impl::LsVectorised || Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field is overwritten below + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + MobiusFermion *clone = new MobiusFermion(Uplaceholder, + *split->FermionGrid, + *split->FermionRBGrid, + *split->GaugeGrid, + *split->GaugeRBGrid, + this->mass_plus, + this->M5, + this->_b, + this->_c, + this->Params); + this->CloneCoefficientsInto(*clone); + this->SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + }; NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/WilsonCloverFermion.h b/Grid/qcd/action/fermion/WilsonCloverFermion.h index 562f08f4a..ecad7f2c9 100644 --- a/Grid/qcd/action/fermion/WilsonCloverFermion.h +++ b/Grid/qcd/action/fermion/WilsonCloverFermion.h @@ -91,6 +91,60 @@ public: // Derivative parts unpreconditioned pseudofermions void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag); + // Install already computed clover term and inverse, and derive the checkerboarded parts + void ImportCloverTerms(const CloverField &C, const CloverField &Cinv) + { + CloverTerm = C; + CloverTermInv = Cinv; + pickCheckerboard(Even, CloverTermEven, CloverTerm); + pickCheckerboard(Odd, CloverTermOdd, CloverTerm); + pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm)); + pickCheckerboard(Odd, CloverTermDagOdd, adj(CloverTerm)); + pickCheckerboard(Even, CloverTermInvEven, CloverTermInv); + pickCheckerboard(Odd, CloverTermInvOdd, CloverTermInv); + pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv)); + pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv)); + } + + // Exact type only: derived operators must not inherit this + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + if ( typeid(*this) != typeid(WilsonCloverFermion) ) { + return nullptr; + } + if ( Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field and clover terms are overwritten below + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + WilsonCloverFermion *clone = + new WilsonCloverFermion(Uplaceholder, + *split->GaugeGrid, + *split->GaugeRBGrid, + this->mass, + csw_r, + csw_t, + this->anisotropyCoeff, + this->Params); + // The stored coefficients are post-scaled; the constructor would scale them again + clone->csw_r = csw_r; + clone->csw_t = csw_t; + + CloverField sC(split->GaugeGrid); + CloverField sCinv(split->GaugeGrid); + Grid_split(CloverTerm,sC); + Grid_split(CloverTermInv,sCinv); + clone->ImportCloverTerms(sC,sCinv); + + this->SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + public: // here fixing the 4 dimensions, make it more general? diff --git a/Grid/qcd/action/fermion/WilsonFermion.h b/Grid/qcd/action/fermion/WilsonFermion.h index 41a2471e2..0a2bb570e 100644 --- a/Grid/qcd/action/fermion/WilsonFermion.h +++ b/Grid/qcd/action/fermion/WilsonFermion.h @@ -147,6 +147,48 @@ public: // DoubleStore impl dependent void ImportGauge(const GaugeField &_Umu); + // Install an already doubled gauge field (phases, twists, anisotropy and -1/2 applied) + void ImportDoubledGauge(const DoubledGaugeField &U) + { + Umu = U; + pickCheckerboard(Even,UmuEven,Umu); + pickCheckerboard(Odd ,UmuOdd ,Umu); + } + + // Redistribute this operator's doubled gauge field onto clone's split grid. Collective. + void SplitDoubledGaugeInto(WilsonFermion &clone) + { + DoubledGaugeField sUmu(clone.GaugeGrid()); + Grid_split(Umu,sUmu); + clone.ImportDoubledGauge(sUmu); + } + + // Exact type only: derived operators (e.g. WilsonTMFermion) must not inherit this + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + if ( typeid(*this) != typeid(WilsonFermion) ) { + return nullptr; + } + if ( Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field is overwritten below + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + WilsonFermion *clone = new WilsonFermion(Uplaceholder, + *split->GaugeGrid, + *split->GaugeRBGrid, + mass, + this->Params, + anisotropyCoeff); + SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + /////////////////////////////////////////////////////////////// // Data members require to support the functionality /////////////////////////////////////////////////////////////// diff --git a/Grid/qcd/action/fermion/WilsonFermion5D.h b/Grid/qcd/action/fermion/WilsonFermion5D.h index ea0eaa32e..e65d086c8 100644 --- a/Grid/qcd/action/fermion/WilsonFermion5D.h +++ b/Grid/qcd/action/fermion/WilsonFermion5D.h @@ -181,6 +181,22 @@ public: // DoubleStore void ImportGauge(const GaugeField &_Umu); + // Install an already doubled gauge field (phases, twists and -1/2 applied) + void ImportDoubledGauge(const DoubledGaugeField &U) + { + Umu = U; + pickCheckerboard(Even,UmuEven,Umu); + pickCheckerboard(Odd ,UmuOdd ,Umu); + } + + // Redistribute this operator's doubled gauge field onto clone's split grid. Collective. + void SplitDoubledGaugeInto(WilsonFermion5D &clone) + { + DoubledGaugeField sUmu(clone.GaugeGrid()); + Grid_split(Umu,sUmu); + clone.ImportDoubledGauge(sUmu); + } + /////////////////////////////////////////////////////////////// // Data members require to support the functionality /////////////////////////////////////////////////////////////// diff --git a/Grid/qcd/action/fermion/ZMobiusFermion.h b/Grid/qcd/action/fermion/ZMobiusFermion.h index de832b7d7..2fa1b8d8f 100644 --- a/Grid/qcd/action/fermion/ZMobiusFermion.h +++ b/Grid/qcd/action/fermion/ZMobiusFermion.h @@ -67,6 +67,44 @@ public: this->SetCoefficientsInternal(1.0,zgamma,b,c); } + // Exact type only: derived operators must not inherit this + virtual SplitOperator *SplitClone(const Coordinate &mpi_split) + { + if ( typeid(*this) != typeid(ZMobiusFermion) ) { + return nullptr; + } + if ( this->Dirichlet || Impl::LsVectorised || Impl::isGparity ) { + return nullptr; + } + SplitOperator *split = this->MakeSplitGrids(mpi_split); + + // Placeholder links; the doubled field is overwritten below + GaugeField Uplaceholder(split->GaugeGrid); + Uplaceholder = Zero(); + + // Constructor argument only; CloneCoefficientsInto restores _gamma at full precision + std::vector gamma(this->Ls); + for ( int s=0; sLs; s++ ) { + gamma[s] = this->_gamma[s]; + } + + ZMobiusFermion *clone = new ZMobiusFermion(Uplaceholder, + *split->FermionGrid, + *split->FermionRBGrid, + *split->GaugeGrid, + *split->GaugeRBGrid, + this->mass_plus, + this->M5, + gamma, + this->_b, + this->_c, + this->Params); + this->CloneCoefficientsInto(*clone); + this->SplitDoubledGaugeInto(*clone); + split->Matrix = clone; + return split; + } + }; NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h b/Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h index 29408ea29..38e419bd8 100644 --- a/Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h +++ b/Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h @@ -64,6 +64,20 @@ public: GRID_ASSERT(U.Checkerboard()==Odd); GRID_ASSERT(V.Checkerboard()==U.Checkerboard()); + // Forces returned on the full 4D grid (checkerboard including s): same sequence, summed + if ( !this->_Mat.CheckerboardedForce() ) { + GaugeField ForceO(Force.Grid()); + GaugeField ForceE(Force.Grid()); + this->_Mat.Meooe (V,tmp1); + this->_Mat.MooeeInv(tmp1,tmp2); + this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo); + this->_Mat.MeooeDag (U,tmp1); + this->_Mat.MooeeInvDag(tmp1,tmp2); + this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerNo); + Force = -(ForceO+ForceE); + return; + } + // NOTE Guido: WE DO NOT WANT TO USE THE ucbgrid GRID FOR THE FORCE // it is not conformable with the HMC force field // Case: Ls vectorised fields @@ -109,6 +123,20 @@ public: GRID_ASSERT(V.Checkerboard()==Odd); GRID_ASSERT(V.Checkerboard()==V.Checkerboard()); + // Forces returned on the full 4D grid (checkerboard including s): same sequence, summed + if ( !this->_Mat.CheckerboardedForce() ) { + GaugeField ForceO(Force.Grid()); + GaugeField ForceE(Force.Grid()); + this->_Mat.MeooeDag (V,tmp1); + this->_Mat.MooeeInvDag(tmp1,tmp2); + this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerYes); + this->_Mat.Meooe (U,tmp1); + this->_Mat.MooeeInv(tmp1,tmp2); + this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerYes); + Force = -(ForceO+ForceE); + return; + } + // NOTE Guido: WE DO NOT WANT TO USE THE ucbgrid GRID FOR THE FORCE // it is not conformable with the HMC force field // INHERIT FROM THE Force field instead