portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Namesapce

Browse Source
This commit is contained in:
paboyle 2018-01-14 22:49:48 +00:00
parent a0380fad72
commit ff7e773d5e
1 changed files with 100 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

201
lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h
View File

@ -1,4 +1,4 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
@ -24,122 +24,121 @@ Author: Guido Cossu <guido.cossu@ed.ac.uk>
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory 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 #ifndef QCD_EVEN_ODD_SCHUR_DIFFERENTIABLE_H
#define QCD_EVEN_ODD_SCHUR_DIFFERENTIABLE_H #define QCD_EVEN_ODD_SCHUR_DIFFERENTIABLE_H
namespace Grid{ NAMESPACE_BEGIN(Grid);
namespace QCD{
// Base even odd HMC on the normal Mee based schur decomposition. // Base even odd HMC on the normal Mee based schur decomposition.
// //
// M = (Mee Meo) = (1 0 ) (Mee 0 ) (1 Mee^{-1} Meo) // 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 ) // (Moe Moo) (Moe Mee^-1 1 ) (0 Moo-Moe Mee^-1 Meo) (0 1 )
// //
// Determinant is det of middle factor // Determinant is det of middle factor
// This assumes Mee is indept of U. // This assumes Mee is indept of U.
// //
template<class Impl> template<class Impl>
class SchurDifferentiableOperator : public SchurDiagMooeeOperator<FermionOperator<Impl>,typename Impl::FermionField> class SchurDifferentiableOperator : public SchurDiagMooeeOperator<FermionOperator<Impl>,typename Impl::FermionField>
{ {
public: public:
INHERIT_IMPL_TYPES(Impl); INHERIT_IMPL_TYPES(Impl);
typedef FermionOperator<Impl> Matrix; typedef FermionOperator<Impl> Matrix;
SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator<Matrix,FermionField>(Mat) {}; SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator<Matrix,FermionField>(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 *fgrid = this->_Mat.FermionGrid();
GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid();
FermionField tmp1(fcbgrid); FermionField tmp1(fcbgrid);
FermionField tmp2(fcbgrid); FermionField tmp2(fcbgrid);
conformable(fcbgrid,U._grid); conformable(fcbgrid,U._grid);
conformable(fcbgrid,V._grid); conformable(fcbgrid,V._grid);
// Assert the checkerboard?? or code for either // Assert the checkerboard?? or code for either
assert(U.checkerboard==Odd); assert(U.checkerboard==Odd);
assert(V.checkerboard==U.checkerboard); assert(V.checkerboard==U.checkerboard);
// NOTE Guido: WE DO NOT WANT TO USE THE ucbgrid GRID FOR THE FORCE // NOTE Guido: WE DO NOT WANT TO USE THE ucbgrid GRID FOR THE FORCE
// it is not conformable with the HMC force field // it is not conformable with the HMC force field
// Case: Ls vectorised fields // Case: Ls vectorised fields
// INHERIT FROM THE Force field instead // INHERIT FROM THE Force field instead
GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid); GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid);
GaugeField ForceO(forcecb); GaugeField ForceO(forcecb);
GaugeField ForceE(forcecb); GaugeField ForceE(forcecb);
// X^dag Der_oe MeeInv Meo Y // X^dag Der_oe MeeInv Meo Y
// Use Mooee as nontrivial but gauge field indept // Use Mooee as nontrivial but gauge field indept
this->_Mat.Meooe (V,tmp1); // odd->even -- implicit -0.5 factor to be applied this->_Mat.Meooe (V,tmp1); // odd->even -- implicit -0.5 factor to be applied
this->_Mat.MooeeInv(tmp1,tmp2); // even->even this->_Mat.MooeeInv(tmp1,tmp2); // even->even
this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo); this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo);
// Accumulate X^dag M_oe MeeInv Der_eo Y // 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.MeooeDag (U,tmp1); // even->odd -- implicit -0.5 factor to be applied
this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even
this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerNo); this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerNo);
assert(ForceE.checkerboard==Even); assert(ForceE.checkerboard==Even);
assert(ForceO.checkerboard==Odd); assert(ForceO.checkerboard==Odd);
setCheckerboard(Force,ForceE); setCheckerboard(Force,ForceE);
setCheckerboard(Force,ForceO); setCheckerboard(Force,ForceO);
Force=-Force; 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;
}
};
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 #endif