diff --git a/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h b/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h index 90b6dbaa..b19f608f 100644 --- a/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h +++ b/lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h @@ -1,4 +1,4 @@ - /************************************************************************************* +/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -24,122 +24,121 @@ Author: Guido Cossu 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 */ +*************************************************************************************/ +/* END LEGAL */ #ifndef QCD_EVEN_ODD_SCHUR_DIFFERENTIABLE_H #define QCD_EVEN_ODD_SCHUR_DIFFERENTIABLE_H -namespace Grid{ - namespace QCD{ +NAMESPACE_BEGIN(Grid); - // Base even odd HMC on the normal Mee based schur decomposition. - // - // M = (Mee Meo) = (1 0 ) (Mee 0 ) (1 Mee^{-1} Meo) - // (Moe Moo) (Moe Mee^-1 1 ) (0 Moo-Moe Mee^-1 Meo) (0 1 ) - // - // Determinant is det of middle factor - // This assumes Mee is indept of U. - // - template - class SchurDifferentiableOperator : public SchurDiagMooeeOperator,typename Impl::FermionField> - { - public: - INHERIT_IMPL_TYPES(Impl); +// Base even odd HMC on the normal Mee based schur decomposition. +// +// M = (Mee Meo) = (1 0 ) (Mee 0 ) (1 Mee^{-1} Meo) +// (Moe Moo) (Moe Mee^-1 1 ) (0 Moo-Moe Mee^-1 Meo) (0 1 ) +// +// Determinant is det of middle factor +// This assumes Mee is indept of U. +// +template +class SchurDifferentiableOperator : public SchurDiagMooeeOperator,typename Impl::FermionField> +{ +public: + INHERIT_IMPL_TYPES(Impl); - typedef FermionOperator Matrix; + typedef FermionOperator Matrix; - SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator(Mat) {}; + SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator(Mat) {}; - void MpcDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { + void MpcDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { - GridBase *fgrid = this->_Mat.FermionGrid(); - GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); + GridBase *fgrid = this->_Mat.FermionGrid(); + GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); - FermionField tmp1(fcbgrid); - FermionField tmp2(fcbgrid); + FermionField tmp1(fcbgrid); + FermionField tmp2(fcbgrid); - conformable(fcbgrid,U._grid); - conformable(fcbgrid,V._grid); + conformable(fcbgrid,U._grid); + conformable(fcbgrid,V._grid); - // Assert the checkerboard?? or code for either - assert(U.checkerboard==Odd); - assert(V.checkerboard==U.checkerboard); + // Assert the checkerboard?? or code for either + assert(U.checkerboard==Odd); + assert(V.checkerboard==U.checkerboard); - // 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 - // INHERIT FROM THE Force field instead - GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid); - GaugeField ForceO(forcecb); - GaugeField ForceE(forcecb); + // 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 + // INHERIT FROM THE Force field instead + GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid); + GaugeField ForceO(forcecb); + GaugeField ForceE(forcecb); - // X^dag Der_oe MeeInv Meo Y - // Use Mooee as nontrivial but gauge field indept - this->_Mat.Meooe (V,tmp1); // odd->even -- implicit -0.5 factor to be applied - this->_Mat.MooeeInv(tmp1,tmp2); // even->even - this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo); - // Accumulate X^dag M_oe MeeInv Der_eo Y - this->_Mat.MeooeDag (U,tmp1); // even->odd -- implicit -0.5 factor to be applied - this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even - this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerNo); + // X^dag Der_oe MeeInv Meo Y + // Use Mooee as nontrivial but gauge field indept + this->_Mat.Meooe (V,tmp1); // odd->even -- implicit -0.5 factor to be applied + this->_Mat.MooeeInv(tmp1,tmp2); // even->even + this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo); + // Accumulate X^dag M_oe MeeInv Der_eo Y + this->_Mat.MeooeDag (U,tmp1); // even->odd -- implicit -0.5 factor to be applied + this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even + this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerNo); - assert(ForceE.checkerboard==Even); - assert(ForceO.checkerboard==Odd); + assert(ForceE.checkerboard==Even); + assert(ForceO.checkerboard==Odd); - setCheckerboard(Force,ForceE); - setCheckerboard(Force,ForceO); - Force=-Force; - - delete forcecb; - } - - - void MpcDagDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { - - GridBase *fgrid = this->_Mat.FermionGrid(); - GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); - - FermionField tmp1(fcbgrid); - FermionField tmp2(fcbgrid); - - conformable(fcbgrid,U._grid); - conformable(fcbgrid,V._grid); - - // Assert the checkerboard?? or code for either - assert(V.checkerboard==Odd); - assert(V.checkerboard==V.checkerboard); - - // 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 - GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid); - GaugeField ForceO(forcecb); - GaugeField ForceE(forcecb); - - // X^dag Der_oe MeeInv Meo Y - // Use Mooee as nontrivial but gauge field indept - this->_Mat.MeooeDag (V,tmp1); // odd->even -- implicit -0.5 factor to be applied - this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even - this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerYes); - - // Accumulate X^dag M_oe MeeInv Der_eo Y - this->_Mat.Meooe (U,tmp1); // even->odd -- implicit -0.5 factor to be applied - this->_Mat.MooeeInv(tmp1,tmp2); // even->even - this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerYes); - - assert(ForceE.checkerboard==Even); - assert(ForceO.checkerboard==Odd); - - setCheckerboard(Force,ForceE); - setCheckerboard(Force,ForceO); - Force=-Force; - - delete forcecb; - } - - }; + setCheckerboard(Force,ForceE); + setCheckerboard(Force,ForceO); + Force=-Force; + delete forcecb; } -} + + + void MpcDagDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { + + GridBase *fgrid = this->_Mat.FermionGrid(); + GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); + + FermionField tmp1(fcbgrid); + FermionField tmp2(fcbgrid); + + conformable(fcbgrid,U._grid); + conformable(fcbgrid,V._grid); + + // Assert the checkerboard?? or code for either + assert(V.checkerboard==Odd); + assert(V.checkerboard==V.checkerboard); + + // 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 + GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid); + GaugeField ForceO(forcecb); + GaugeField ForceE(forcecb); + + // X^dag Der_oe MeeInv Meo Y + // Use Mooee as nontrivial but gauge field indept + this->_Mat.MeooeDag (V,tmp1); // odd->even -- implicit -0.5 factor to be applied + this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even + this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerYes); + + // Accumulate X^dag M_oe MeeInv Der_eo Y + this->_Mat.Meooe (U,tmp1); // even->odd -- implicit -0.5 factor to be applied + this->_Mat.MooeeInv(tmp1,tmp2); // even->even + this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerYes); + + assert(ForceE.checkerboard==Even); + assert(ForceO.checkerboard==Odd); + + setCheckerboard(Force,ForceE); + setCheckerboard(Force,ForceO); + Force=-Force; + + delete forcecb; + } + +}; + +NAMESPACE_END(Grid); + #endif