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

Minor modifications

Browse Source
This commit is contained in:
Guido Cossu 2016-07-06 11:41:27 +01:00
parent fc4a043663
commit ffedeb1c58
4 changed files with 1029 additions and 1025 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

564
lib/qcd/action/fermion/FermionOperatorImpl.h
View File

@ -1,99 +1,104 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/fermion/FermionOperatorImpl.h Source file: ./lib/qcd/action/fermion/FermionOperatorImpl.h
Copyright (C) 2015 Copyright (C) 2015
Author: Peter Boyle <pabobyle@ph.ed.ac.uk> Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local> Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
Author: paboyle <paboyle@ph.ed.ac.uk> Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along 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., with this program; if not, write to the Free Software Foundation, Inc.,
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 GRID_QCD_FERMION_OPERATOR_IMPL_H #ifndef GRID_QCD_FERMION_OPERATOR_IMPL_H
#define GRID_QCD_FERMION_OPERATOR_IMPL_H #define GRID_QCD_FERMION_OPERATOR_IMPL_H
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
//////////////////////////////////////////////
// Template parameter class constructs to package
// externally control Fermion implementations
// in orthogonal directions
//
// Ultimately need Impl to always define types where XXX is opaque
//
// typedef typename XXX Simd;
// typedef typename XXX GaugeLinkField;
// typedef typename XXX GaugeField;
// typedef typename XXX GaugeActField;
// typedef typename XXX FermionField;
// typedef typename XXX DoubledGaugeField;
// typedef typename XXX SiteSpinor;
// typedef typename XXX SiteHalfSpinor;
// typedef typename XXX Compressor;
//
// and Methods:
// void ImportGauge(GridBase *GaugeGrid,DoubledGaugeField &Uds,const
// GaugeField &Umu)
// void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const
// GaugeField &Umu)
// void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const
// SiteHalfSpinor &chi,int mu,StencilEntry *SE,StencilImpl &St)
// void InsertForce4D(GaugeField &mat,const FermionField &Btilde,const
// FermionField &A,int mu)
// void InsertForce5D(GaugeField &mat,const FermionField &Btilde,const
// FermionField &A,int mu)
//
//
// To acquire the typedefs from "Base" (either a base class or template param)
// use:
//
// INHERIT_GIMPL_TYPES(Base)
// INHERIT_FIMPL_TYPES(Base)
// INHERIT_IMPL_TYPES(Base)
//
// The Fermion operators will do the following:
//
// struct MyOpParams {
// RealD mass;
// };
//
//
// template<class Impl>
// class MyOp : pubic<Impl> {
// public:
//
// INHERIT_ALL_IMPL_TYPES(Impl);
//
// MyOp(MyOpParams Myparm, ImplParams &ImplParam) : Impl(ImplParam)
// {
//
// };
//
// }
//////////////////////////////////////////////
////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
// Template parameter class constructs to package // Implementation dependent fermion types
// externally control Fermion implementations ////////////////////////////////////////////////////////////////////////
// in orthogonal directions
//
// Ultimately need Impl to always define types where XXX is opaque
//
// typedef typename XXX Simd;
// typedef typename XXX GaugeLinkField;
// typedef typename XXX GaugeField;
// typedef typename XXX GaugeActField;
// typedef typename XXX FermionField;
// typedef typename XXX DoubledGaugeField;
// typedef typename XXX SiteSpinor;
// typedef typename XXX SiteHalfSpinor;
// typedef typename XXX Compressor;
//
// and Methods:
// void ImportGauge(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
// void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
// void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,StencilImpl &St)
// void InsertForce4D(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu)
// void InsertForce5D(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu)
//
//
// To acquire the typedefs from "Base" (either a base class or template param) use:
//
// INHERIT_GIMPL_TYPES(Base)
// INHERIT_FIMPL_TYPES(Base)
// INHERIT_IMPL_TYPES(Base)
//
// The Fermion operators will do the following:
//
// struct MyOpParams {
// RealD mass;
// };
//
//
// template<class Impl>
// class MyOp : pubic<Impl> {
// public:
//
// INHERIT_ALL_IMPL_TYPES(Impl);
//
// MyOp(MyOpParams Myparm, ImplParams &ImplParam) : Impl(ImplParam)
// {
//
// };
//
// }
//////////////////////////////////////////////
#define INHERIT_FIMPL_TYPES(Impl) \
////////////////////////////////////////////////////////////////////////
// Implementation dependent fermion types
////////////////////////////////////////////////////////////////////////
#define INHERIT_FIMPL_TYPES(Impl)\
typedef typename Impl::FermionField FermionField; \ typedef typename Impl::FermionField FermionField; \
typedef typename Impl::DoubledGaugeField DoubledGaugeField; \ typedef typename Impl::DoubledGaugeField DoubledGaugeField; \
typedef typename Impl::SiteSpinor SiteSpinor; \ typedef typename Impl::SiteSpinor SiteSpinor; \
@ -103,219 +108,243 @@ namespace Grid {
typedef typename Impl::ImplParams ImplParams; typedef typename Impl::ImplParams ImplParams;
#define INHERIT_IMPL_TYPES(Base) \ #define INHERIT_IMPL_TYPES(Base) \
INHERIT_GIMPL_TYPES(Base)\ INHERIT_GIMPL_TYPES(Base) \
INHERIT_FIMPL_TYPES(Base) INHERIT_FIMPL_TYPES(Base)
/////// ///////
// Single flavour four spinors with colour index // Single flavour four spinors with colour index
/////// ///////
template<class S,int Nrepresentation=Nc> template <class S, int Nrepresentation = Nc>
class WilsonImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepresentation> > { class WilsonImpl
: public PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > {
public: public:
typedef PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > Gimpl;
typedef PeriodicGaugeImpl< GaugeImplTypes< S,Nrepresentation> > Gimpl;
INHERIT_GIMPL_TYPES(Gimpl); INHERIT_GIMPL_TYPES(Gimpl);
template<typename vtype> using iImplSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >; template <typename vtype>
template<typename vtype> using iImplHalfSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Nhs> >; using iImplSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
template<typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds >; template <typename vtype>
using iImplHalfSpinor =
iScalar<iVector<iVector<vtype, Nrepresentation>, Nhs> >;
template <typename vtype>
using iImplDoubledGaugeField =
iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>;
typedef iImplSpinor <Simd> SiteSpinor; typedef iImplSpinor<Simd> SiteSpinor;
typedef iImplHalfSpinor<Simd> SiteHalfSpinor; typedef iImplHalfSpinor<Simd> SiteHalfSpinor;
typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField; typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
typedef Lattice<SiteSpinor> FermionField; typedef Lattice<SiteSpinor> FermionField;
typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField; typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
typedef WilsonCompressor<SiteHalfSpinor,SiteSpinor> Compressor; typedef WilsonCompressor<SiteHalfSpinor, SiteSpinor> Compressor;
typedef WilsonImplParams ImplParams; typedef WilsonImplParams ImplParams;
typedef WilsonStencil<SiteSpinor,SiteHalfSpinor> StencilImpl; typedef WilsonStencil<SiteSpinor, SiteHalfSpinor> StencilImpl;
ImplParams Params; ImplParams Params;
WilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {}; WilsonImpl(const ImplParams &p = ImplParams()) : Params(p){};
bool overlapCommsCompute(void) { return Params.overlapCommsCompute; }; bool overlapCommsCompute(void) { return Params.overlapCommsCompute; };
inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,StencilImpl &St){ inline void multLink(SiteHalfSpinor &phi, const SiteDoubledGaugeField &U,
mult(&phi(),&U(mu),&chi()); const SiteHalfSpinor &chi, int mu, StencilEntry *SE,
StencilImpl &St) {
mult(&phi(), &U(mu), &chi());
} }
template<class ref> template <class ref>
inline void loadLinkElement(Simd & reg,ref &memory){ inline void loadLinkElement(Simd &reg, ref &memory) {
reg = memory; reg = memory;
} }
inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) inline void DoubleStore(GridBase *GaugeGrid, DoubledGaugeField &Uds,
{ const GaugeField &Umu) {
conformable(Uds._grid,GaugeGrid); conformable(Uds._grid, GaugeGrid);
conformable(Umu._grid,GaugeGrid); conformable(Umu._grid, GaugeGrid);
GaugeLinkField U(GaugeGrid); GaugeLinkField U(GaugeGrid);
for(int mu=0;mu<Nd;mu++){ for (int mu = 0; mu < Nd; mu++) {
U = PeekIndex<LorentzIndex>(Umu,mu); U = PeekIndex<LorentzIndex>(Umu, mu);
PokeIndex<LorentzIndex>(Uds,U,mu); PokeIndex<LorentzIndex>(Uds, U, mu);
U = adj(Cshift(U,mu,-1)); U = adj(Cshift(U, mu, -1));
PokeIndex<LorentzIndex>(Uds,U,mu+4); PokeIndex<LorentzIndex>(Uds, U, mu + 4);
} }
} }
inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){ inline void InsertForce4D(GaugeField &mat, FermionField &Btilde,
FermionField &A, int mu) {
GaugeLinkField link(mat._grid); GaugeLinkField link(mat._grid);
link = TraceIndex<SpinIndex>(outerProduct(Btilde,A)); link = TraceIndex<SpinIndex>(outerProduct(Btilde, A));
PokeIndex<LorentzIndex>(mat,link,mu); PokeIndex<LorentzIndex>(mat, link, mu);
} }
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){ inline void InsertForce5D(GaugeField &mat, FermionField &Btilde,
FermionField &Atilde, int mu) {
int Ls=Btilde._grid->_fdimensions[0]; int Ls = Btilde._grid->_fdimensions[0];
GaugeLinkField tmp(mat._grid); GaugeLinkField tmp(mat._grid);
tmp = zero; tmp = zero;
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int sss=0;sss<tmp._grid->oSites();sss++){ for (int sss = 0; sss < tmp._grid->oSites(); sss++) {
int sU=sss; int sU = sss;
for(int s=0;s<Ls;s++){ for (int s = 0; s < Ls; s++) {
int sF = s+Ls*sU; int sF = s + Ls * sU;
tmp[sU] = tmp[sU]+ traceIndex<SpinIndex>(outerProduct(Btilde[sF],Atilde[sF])); // ordering here tmp[sU] = tmp[sU] + traceIndex<SpinIndex>(outerProduct(
Btilde[sF], Atilde[sF])); // ordering here
} }
} }
PokeIndex<LorentzIndex>(mat,tmp,mu); PokeIndex<LorentzIndex>(mat, tmp, mu);
} }
};
}; ///////
// Single flavour four spinors with colour index, 5d redblack
///////
template <class S, int Nrepresentation = Nc>
/////// class DomainWallRedBlack5dImpl
// Single flavour four spinors with colour index, 5d redblack : public PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > {
///////
template<class S,int Nrepresentation=Nc>
class DomainWallRedBlack5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepresentation> > {
public: public:
typedef PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > Gimpl;
typedef PeriodicGaugeImpl< GaugeImplTypes< S,Nrepresentation> > Gimpl;
INHERIT_GIMPL_TYPES(Gimpl); INHERIT_GIMPL_TYPES(Gimpl);
template<typename vtype> using iImplSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >; template <typename vtype>
template<typename vtype> using iImplHalfSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Nhs> >; using iImplSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
template<typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds >; template <typename vtype>
template<typename vtype> using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd >; using iImplHalfSpinor =
template<typename vtype> using iImplGaugeLink = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >; iScalar<iVector<iVector<vtype, Nrepresentation>, Nhs> >;
template <typename vtype>
using iImplDoubledGaugeField =
iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>;
template <typename vtype>
using iImplGaugeField =
iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd>;
template <typename vtype>
using iImplGaugeLink = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
typedef iImplSpinor <Simd> SiteSpinor; typedef iImplSpinor<Simd> SiteSpinor;
typedef iImplHalfSpinor<Simd> SiteHalfSpinor; typedef iImplHalfSpinor<Simd> SiteHalfSpinor;
typedef Lattice<SiteSpinor> FermionField; typedef Lattice<SiteSpinor> FermionField;
// Make the doubled gauge field a *scalar* // Make the doubled gauge field a *scalar*
typedef iImplDoubledGaugeField<typename Simd::scalar_type> SiteDoubledGaugeField; // This is a scalar typedef iImplDoubledGaugeField<typename Simd::scalar_type>
typedef iImplGaugeField<typename Simd::scalar_type> SiteScalarGaugeField; // scalar SiteDoubledGaugeField; // This is a scalar
typedef iImplGaugeLink <typename Simd::scalar_type> SiteScalarGaugeLink; // scalar typedef iImplGaugeField<typename Simd::scalar_type>
SiteScalarGaugeField; // scalar
typedef iImplGaugeLink<typename Simd::scalar_type>
SiteScalarGaugeLink; // scalar
typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField; typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
typedef WilsonCompressor<SiteHalfSpinor,SiteSpinor> Compressor; typedef WilsonCompressor<SiteHalfSpinor, SiteSpinor> Compressor;
typedef WilsonImplParams ImplParams; typedef WilsonImplParams ImplParams;
typedef WilsonStencil<SiteSpinor,SiteHalfSpinor> StencilImpl; typedef WilsonStencil<SiteSpinor, SiteHalfSpinor> StencilImpl;
ImplParams Params; ImplParams Params;
DomainWallRedBlack5dImpl(const ImplParams &p= ImplParams()) : Params(p) {}; DomainWallRedBlack5dImpl(const ImplParams &p = ImplParams()) : Params(p){};
bool overlapCommsCompute(void) { return false; }; bool overlapCommsCompute(void) { return false; };
template<class ref> template <class ref>
inline void loadLinkElement(Simd & reg,ref &memory){ inline void loadLinkElement(Simd &reg, ref &memory) {
vsplat(reg,memory); vsplat(reg, memory);
} }
inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,StencilImpl &St) inline void multLink(SiteHalfSpinor &phi, const SiteDoubledGaugeField &U,
{ const SiteHalfSpinor &chi, int mu, StencilEntry *SE,
StencilImpl &St) {
SiteGaugeLink UU; SiteGaugeLink UU;
for(int i=0;i<Nrepresentation;i++){ for (int i = 0; i < Nrepresentation; i++) {
for(int j=0;j<Nrepresentation;j++){ for (int j = 0; j < Nrepresentation; j++) {
vsplat(UU()()(i,j),U(mu)()(i,j)); vsplat(UU()()(i, j), U(mu)()(i, j));
} }
} }
mult(&phi(),&UU(),&chi()); mult(&phi(), &UU(), &chi());
} }
inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) inline void DoubleStore(GridBase *GaugeGrid, DoubledGaugeField &Uds,
{ const GaugeField &Umu) {
SiteScalarGaugeField ScalarUmu; SiteScalarGaugeField ScalarUmu;
SiteDoubledGaugeField ScalarUds; SiteDoubledGaugeField ScalarUds;
GaugeLinkField U (Umu._grid); GaugeLinkField U(Umu._grid);
GaugeField Uadj(Umu._grid); GaugeField Uadj(Umu._grid);
for(int mu=0;mu<Nd;mu++){ for (int mu = 0; mu < Nd; mu++) {
U = PeekIndex<LorentzIndex>(Umu,mu); U = PeekIndex<LorentzIndex>(Umu, mu);
U = adj(Cshift(U,mu,-1)); U = adj(Cshift(U, mu, -1));
PokeIndex<LorentzIndex>(Uadj,U,mu); PokeIndex<LorentzIndex>(Uadj, U, mu);
} }
for(int lidx=0;lidx<GaugeGrid->lSites();lidx++){ for (int lidx = 0; lidx < GaugeGrid->lSites(); lidx++) {
std::vector<int> lcoor; std::vector<int> lcoor;
GaugeGrid->LocalIndexToLocalCoor(lidx,lcoor); GaugeGrid->LocalIndexToLocalCoor(lidx, lcoor);
peekLocalSite(ScalarUmu,Umu,lcoor); peekLocalSite(ScalarUmu, Umu, lcoor);
for(int mu=0;mu<4;mu++) ScalarUds(mu) = ScalarUmu(mu); for (int mu = 0; mu < 4; mu++) ScalarUds(mu) = ScalarUmu(mu);
peekLocalSite(ScalarUmu,Uadj,lcoor); peekLocalSite(ScalarUmu, Uadj, lcoor);
for(int mu=0;mu<4;mu++) ScalarUds(mu+4) = ScalarUmu(mu); for (int mu = 0; mu < 4; mu++) ScalarUds(mu + 4) = ScalarUmu(mu);
pokeLocalSite(ScalarUds,Uds,lcoor); pokeLocalSite(ScalarUds, Uds, lcoor);
}
} }
} inline void InsertForce4D(GaugeField &mat, FermionField &Btilde,
FermionField &A, int mu) {
inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){
assert(0); assert(0);
} }
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){ inline void InsertForce5D(GaugeField &mat, FermionField &Btilde,
FermionField &Atilde, int mu) {
assert(0); assert(0);
} }
};
}; ////////////////////////////////////////////////////////////////////////////////////////
// Flavour doubled spinors; is Gparity the only? what about C*?
////////////////////////////////////////////////////////////////////////////////////////
template <class S, int Nrepresentation>
//////////////////////////////////////////////////////////////////////////////////////// class GparityWilsonImpl
// Flavour doubled spinors; is Gparity the only? what about C*? : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresentation> > {
////////////////////////////////////////////////////////////////////////////////////////
template<class S,int Nrepresentation>
class GparityWilsonImpl : public ConjugateGaugeImpl< GaugeImplTypes<S,Nrepresentation> >{
public: public:
typedef ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresentation> > Gimpl;
typedef ConjugateGaugeImpl< GaugeImplTypes<S,Nrepresentation> > Gimpl;
INHERIT_GIMPL_TYPES(Gimpl); INHERIT_GIMPL_TYPES(Gimpl);
template<typename vtype> using iImplSpinor = iVector<iVector<iVector<vtype, Nrepresentation>, Ns>, Ngp >; template <typename vtype>
template<typename vtype> using iImplHalfSpinor = iVector<iVector<iVector<vtype, Nrepresentation>, Nhs>, Ngp >; using iImplSpinor =
template<typename vtype> using iImplDoubledGaugeField = iVector<iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds >, Ngp >; iVector<iVector<iVector<vtype, Nrepresentation>, Ns>, Ngp>;
template <typename vtype>
using iImplHalfSpinor =
iVector<iVector<iVector<vtype, Nrepresentation>, Nhs>, Ngp>;
template <typename vtype>
using iImplDoubledGaugeField =
iVector<iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nds>, Ngp>;
typedef iImplSpinor <Simd> SiteSpinor; typedef iImplSpinor<Simd> SiteSpinor;
typedef iImplHalfSpinor<Simd> SiteHalfSpinor; typedef iImplHalfSpinor<Simd> SiteHalfSpinor;
typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField; typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
typedef Lattice<SiteSpinor> FermionField; typedef Lattice<SiteSpinor> FermionField;
typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField; typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
typedef WilsonCompressor<SiteHalfSpinor,SiteSpinor> Compressor; typedef WilsonCompressor<SiteHalfSpinor, SiteSpinor> Compressor;
typedef WilsonStencil<SiteSpinor,SiteHalfSpinor> StencilImpl; typedef WilsonStencil<SiteSpinor, SiteHalfSpinor> StencilImpl;
typedef GparityWilsonImplParams ImplParams; typedef GparityWilsonImplParams ImplParams;
ImplParams Params; ImplParams Params;
GparityWilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {}; GparityWilsonImpl(const ImplParams &p = ImplParams()) : Params(p){};
bool overlapCommsCompute(void) { return Params.overlapCommsCompute; }; bool overlapCommsCompute(void) { return Params.overlapCommsCompute; };
// provide the multiply by link that is differentiated between Gparity (with flavour index) and non-Gparity // provide the multiply by link that is differentiated between Gparity (with
inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,StencilImpl &St){ // flavour index) and non-Gparity
inline void multLink(SiteHalfSpinor &phi, const SiteDoubledGaugeField &U,
const SiteHalfSpinor &chi, int mu, StencilEntry *SE,
StencilImpl &St) {
typedef SiteHalfSpinor vobj; typedef SiteHalfSpinor vobj;
typedef typename SiteHalfSpinor::scalar_object sobj; typedef typename SiteHalfSpinor::scalar_object sobj;
@ -335,156 +364,151 @@ PARALLEL_FOR_LOOP
int mmu = mu % Nd; int mmu = mu % Nd;
// assert our assumptions // assert our assumptions
assert((distance==1)||(distance==-1)); // nearest neighbour stencil hard code assert((distance == 1) ||
assert((sl==1)||(sl==2)); (distance == -1)); // nearest neighbour stencil hard code
assert((sl == 1) || (sl == 2));
std::vector<int> icoor; std::vector<int> icoor;
if ( SE->_around_the_world && Params.twists[mmu] ) { if (SE->_around_the_world && Params.twists[mmu]) {
if (sl == 2) {
if ( sl == 2 ) {
std::vector<sobj> vals(Nsimd); std::vector<sobj> vals(Nsimd);
extract(chi,vals); extract(chi, vals);
for(int s=0;s<Nsimd;s++){ for (int s = 0; s < Nsimd; s++) {
grid->iCoorFromIindex(icoor, s);
grid->iCoorFromIindex(icoor,s); assert((icoor[direction] == 0) || (icoor[direction] == 1));
assert((icoor[direction]==0)||(icoor[direction]==1));
int permute_lane; int permute_lane;
if ( distance == 1) { if (distance == 1) {
permute_lane = icoor[direction]?1:0; permute_lane = icoor[direction] ? 1 : 0;
} else { } else {
permute_lane = icoor[direction]?0:1; permute_lane = icoor[direction] ? 0 : 1;
} }
if ( permute_lane ) { if (permute_lane) {
stmp(0) = vals[s](1); stmp(0) = vals[s](1);
stmp(1) = vals[s](0); stmp(1) = vals[s](0);
vals[s] = stmp; vals[s] = stmp;
} }
} }
merge(vtmp,vals); merge(vtmp, vals);
} else { } else {
vtmp(0) = chi(1); vtmp(0) = chi(1);
vtmp(1) = chi(0); vtmp(1) = chi(0);
} }
mult(&phi(0),&U(0)(mu),&vtmp(0)); mult(&phi(0), &U(0)(mu), &vtmp(0));
mult(&phi(1),&U(1)(mu),&vtmp(1)); mult(&phi(1), &U(1)(mu), &vtmp(1));
} else { } else {
mult(&phi(0),&U(0)(mu),&chi(0)); mult(&phi(0), &U(0)(mu), &chi(0));
mult(&phi(1),&U(1)(mu),&chi(1)); mult(&phi(1), &U(1)(mu), &chi(1));
}
} }
} inline void DoubleStore(GridBase *GaugeGrid, DoubledGaugeField &Uds,
const GaugeField &Umu) {
conformable(Uds._grid, GaugeGrid);
conformable(Umu._grid, GaugeGrid);
inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) GaugeLinkField Utmp(GaugeGrid);
{ GaugeLinkField U(GaugeGrid);
conformable(Uds._grid,GaugeGrid);
conformable(Umu._grid,GaugeGrid);
GaugeLinkField Utmp (GaugeGrid);
GaugeLinkField U (GaugeGrid);
GaugeLinkField Uconj(GaugeGrid); GaugeLinkField Uconj(GaugeGrid);
Lattice<iScalar<vInteger> > coor(GaugeGrid); Lattice<iScalar<vInteger> > coor(GaugeGrid);
for (int mu = 0; mu < Nd; mu++) {
LatticeCoordinate(coor, mu);
for(int mu=0;mu<Nd;mu++){ U = PeekIndex<LorentzIndex>(Umu, mu);
LatticeCoordinate(coor,mu);
U = PeekIndex<LorentzIndex>(Umu,mu);
Uconj = conjugate(U); Uconj = conjugate(U);
// This phase could come from a simple bc 1,1,-1,1 .. // This phase could come from a simple bc 1,1,-1,1 ..
int neglink = GaugeGrid->GlobalDimensions()[mu]-1; int neglink = GaugeGrid->GlobalDimensions()[mu] - 1;
if ( Params.twists[mu] ) { if (Params.twists[mu]) {
Uconj = where(coor==neglink,-Uconj,Uconj); Uconj = where(coor == neglink, -Uconj, Uconj);
} }
PARALLEL_FOR_LOOP
PARALLEL_FOR_LOOP for (auto ss = U.begin(); ss < U.end(); ss++) {
for(auto ss=U.begin();ss<U.end();ss++){
Uds[ss](0)(mu) = U[ss](); Uds[ss](0)(mu) = U[ss]();
Uds[ss](1)(mu) = Uconj[ss](); Uds[ss](1)(mu) = Uconj[ss]();
} }
U = adj(Cshift(U ,mu,-1)); // correct except for spanning the boundary U = adj(Cshift(U, mu, -1)); // correct except for spanning the boundary
Uconj = adj(Cshift(Uconj,mu,-1)); Uconj = adj(Cshift(Uconj, mu, -1));
Utmp = U; Utmp = U;
if ( Params.twists[mu] ) { if (Params.twists[mu]) {
Utmp = where(coor==0,Uconj,Utmp); Utmp = where(coor == 0, Uconj, Utmp);
} }
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(auto ss=U.begin();ss<U.end();ss++){ for (auto ss = U.begin(); ss < U.end(); ss++) {
Uds[ss](0)(mu+4) = Utmp[ss](); Uds[ss](0)(mu + 4) = Utmp[ss]();
} }
Utmp = Uconj; Utmp = Uconj;
if ( Params.twists[mu] ) { if (Params.twists[mu]) {
Utmp = where(coor==0,U,Utmp); Utmp = where(coor == 0, U, Utmp);
} }
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(auto ss=U.begin();ss<U.end();ss++){ for (auto ss = U.begin(); ss < U.end(); ss++) {
Uds[ss](1)(mu+4) = Utmp[ss](); Uds[ss](1)(mu + 4) = Utmp[ss]();
} }
} }
} }
inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){ inline void InsertForce4D(GaugeField &mat, FermionField &Btilde,
FermionField &A, int mu) {
// DhopDir provides U or Uconj depending on coor/flavour. // DhopDir provides U or Uconj depending on coor/flavour.
GaugeLinkField link(mat._grid); GaugeLinkField link(mat._grid);
// use lorentz for flavour as hack. // use lorentz for flavour as hack.
auto tmp = TraceIndex<SpinIndex>(outerProduct(Btilde,A)); auto tmp = TraceIndex<SpinIndex>(outerProduct(Btilde, A));
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(auto ss=tmp.begin();ss<tmp.end();ss++){ for (auto ss = tmp.begin(); ss < tmp.end(); ss++) {
link[ss]() = tmp[ss](0,0) - conjugate(tmp[ss](1,1)) ; link[ss]() = tmp[ss](0, 0) - conjugate(tmp[ss](1, 1));
} }
PokeIndex<LorentzIndex>(mat,link,mu); PokeIndex<LorentzIndex>(mat, link, mu);
return; return;
} }
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){ inline void InsertForce5D(GaugeField &mat, FermionField &Btilde,
FermionField &Atilde, int mu) {
int Ls=Btilde._grid->_fdimensions[0]; int Ls = Btilde._grid->_fdimensions[0];
GaugeLinkField tmp(mat._grid); GaugeLinkField tmp(mat._grid);
tmp = zero; tmp = zero;
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int ss=0;ss<tmp._grid->oSites();ss++){ for (int ss = 0; ss < tmp._grid->oSites(); ss++) {
for(int s=0;s<Ls;s++){ for (int s = 0; s < Ls; s++) {
int sF = s+Ls*ss; int sF = s + Ls * ss;
auto ttmp = traceIndex<SpinIndex>(outerProduct(Btilde[sF],Atilde[sF])); auto ttmp = traceIndex<SpinIndex>(outerProduct(Btilde[sF], Atilde[sF]));
tmp[ss]() = tmp[ss]()+ ttmp(0,0) + conjugate(ttmp(1,1)); tmp[ss]() = tmp[ss]() + ttmp(0, 0) + conjugate(ttmp(1, 1));
} }
} }
PokeIndex<LorentzIndex>(mat,tmp,mu); PokeIndex<LorentzIndex>(mat, tmp, mu);
return; return;
} }
}; };
typedef WilsonImpl<vComplex ,Nc> WilsonImplR; // Real.. whichever prec typedef WilsonImpl<vComplex, Nc> WilsonImplR; // Real.. whichever prec
typedef WilsonImpl<vComplexF,Nc> WilsonImplF; // Float typedef WilsonImpl<vComplexF, Nc> WilsonImplF; // Float
typedef WilsonImpl<vComplexD,Nc> WilsonImplD; // Double typedef WilsonImpl<vComplexD, Nc> WilsonImplD; // Double
typedef DomainWallRedBlack5dImpl<vComplex ,Nc> DomainWallRedBlack5dImplR; // Real.. whichever prec typedef DomainWallRedBlack5dImpl<vComplex, Nc>
typedef DomainWallRedBlack5dImpl<vComplexF,Nc> DomainWallRedBlack5dImplF; // Float DomainWallRedBlack5dImplR; // Real.. whichever prec
typedef DomainWallRedBlack5dImpl<vComplexD,Nc> DomainWallRedBlack5dImplD; // Double typedef DomainWallRedBlack5dImpl<vComplexF, Nc>
DomainWallRedBlack5dImplF; // Float
typedef DomainWallRedBlack5dImpl<vComplexD, Nc>
DomainWallRedBlack5dImplD; // Double
typedef GparityWilsonImpl<vComplex ,Nc> GparityWilsonImplR; // Real.. whichever prec typedef GparityWilsonImpl<vComplex, Nc>
typedef GparityWilsonImpl<vComplexF,Nc> GparityWilsonImplF; // Float GparityWilsonImplR; // Real.. whichever prec
typedef GparityWilsonImpl<vComplexD,Nc> GparityWilsonImplD; // Double typedef GparityWilsonImpl<vComplexF, Nc> GparityWilsonImplF; // Float
typedef GparityWilsonImpl<vComplexD, Nc> GparityWilsonImplD; // Double
} }
} }
#endif #endif

410
lib/qcd/action/fermion/WilsonFermion.cc
View File

@ -1,154 +1,150 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/fermion/WilsonFermion.cc Source file: ./lib/qcd/action/fermion/WilsonFermion.cc
Copyright (C) 2015 Copyright (C) 2015
Author: Peter Boyle <pabobyle@ph.ed.ac.uk> Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local> Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
Author: paboyle <paboyle@ph.ed.ac.uk> Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along 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., with this program; if not, write to the Free Software Foundation, Inc.,
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 */
#include <Grid.h> #include <Grid.h>
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
const std::vector<int> WilsonFermionStatic::directions ({0,1,2,3, 0, 1, 2, 3}); const std::vector<int> WilsonFermionStatic::directions({0, 1, 2, 3, 0, 1, 2,
const std::vector<int> WilsonFermionStatic::displacements({1,1,1,1,-1,-1,-1,-1}); 3});
int WilsonFermionStatic::HandOptDslash; const std::vector<int> WilsonFermionStatic::displacements({1, 1, 1, 1, -1, -1,
-1, -1});
int WilsonFermionStatic::HandOptDslash;
///////////////////////////////// /////////////////////////////////
// Constructor and gauge import // Constructor and gauge import
///////////////////////////////// /////////////////////////////////
template<class Impl> template <class Impl>
WilsonFermion<Impl>::WilsonFermion(GaugeField &_Umu, WilsonFermion<Impl>::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, RealD _mass,
GridRedBlackCartesian &Hgrid, const ImplParams &p)
RealD _mass,const ImplParams &p) : : Kernels(p),
Kernels(p),
_grid(&Fgrid), _grid(&Fgrid),
_cbgrid(&Hgrid), _cbgrid(&Hgrid),
Stencil (&Fgrid,npoint,Even,directions,displacements), Stencil(&Fgrid, npoint, Even, directions, displacements),
StencilEven(&Hgrid,npoint,Even,directions,displacements), // source is Even StencilEven(&Hgrid, npoint, Even, directions,
StencilOdd (&Hgrid,npoint,Odd ,directions,displacements), // source is Odd displacements), // source is Even
StencilOdd(&Hgrid, npoint, Odd, directions,
displacements), // source is Odd
mass(_mass), mass(_mass),
Lebesgue(_grid), Lebesgue(_grid),
LebesgueEvenOdd(_cbgrid), LebesgueEvenOdd(_cbgrid),
Umu(&Fgrid), Umu(&Fgrid),
UmuEven(&Hgrid), UmuEven(&Hgrid),
UmuOdd (&Hgrid) UmuOdd(&Hgrid) {
{
// Allocate the required comms buffer // Allocate the required comms buffer
ImportGauge(_Umu); ImportGauge(_Umu);
} }
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::ImportGauge(const GaugeField &_Umu) void WilsonFermion<Impl>::ImportGauge(const GaugeField &_Umu) {
{
GaugeField HUmu(_Umu._grid); GaugeField HUmu(_Umu._grid);
HUmu = _Umu*(-0.5); HUmu = _Umu * (-0.5);
Impl::DoubleStore(GaugeGrid(),Umu,HUmu); Impl::DoubleStore(GaugeGrid(), Umu, HUmu);
pickCheckerboard(Even,UmuEven,Umu); pickCheckerboard(Even, UmuEven, Umu);
pickCheckerboard(Odd ,UmuOdd,Umu); pickCheckerboard(Odd, UmuOdd, Umu);
} }
///////////////////////////// /////////////////////////////
// Implement the interface // Implement the interface
///////////////////////////// /////////////////////////////
template<class Impl> template <class Impl>
RealD WilsonFermion<Impl>::M(const FermionField &in, FermionField &out) RealD WilsonFermion<Impl>::M(const FermionField &in, FermionField &out) {
{ out.checkerboard = in.checkerboard;
out.checkerboard=in.checkerboard; Dhop(in, out, DaggerNo);
Dhop(in,out,DaggerNo); return axpy_norm(out, 4 + mass, in, out);
return axpy_norm(out,4+mass,in,out); }
}
template<class Impl> template <class Impl>
RealD WilsonFermion<Impl>::Mdag(const FermionField &in, FermionField &out) RealD WilsonFermion<Impl>::Mdag(const FermionField &in, FermionField &out) {
{ out.checkerboard = in.checkerboard;
out.checkerboard=in.checkerboard; Dhop(in, out, DaggerYes);
Dhop(in,out,DaggerYes); return axpy_norm(out, 4 + mass, in, out);
return axpy_norm(out,4+mass,in,out); }
}
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::Meooe(const FermionField &in, FermionField &out) void WilsonFermion<Impl>::Meooe(const FermionField &in, FermionField &out) {
{ if (in.checkerboard == Odd) {
if ( in.checkerboard == Odd ) { DhopEO(in, out, DaggerNo);
DhopEO(in,out,DaggerNo);
} else { } else {
DhopOE(in,out,DaggerNo); DhopOE(in, out, DaggerNo);
} }
} }
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
{ if (in.checkerboard == Odd) {
if ( in.checkerboard == Odd ) { DhopEO(in, out, DaggerYes);
DhopEO(in,out,DaggerYes);
} else { } else {
DhopOE(in,out,DaggerYes); DhopOE(in, out, DaggerYes);
}
} }
}
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::Mooee(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::Mooee(const FermionField &in, FermionField &out) {
out.checkerboard = in.checkerboard; out.checkerboard = in.checkerboard;
typename FermionField::scalar_type scal(4.0+mass); typename FermionField::scalar_type scal(4.0 + mass);
out = scal*in; out = scal * in;
} }
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) {
out.checkerboard = in.checkerboard; out.checkerboard = in.checkerboard;
Mooee(in,out); Mooee(in, out);
} }
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) {
out.checkerboard = in.checkerboard; out.checkerboard = in.checkerboard;
out = (1.0/(4.0+mass))*in; out = (1.0 / (4.0 + mass)) * in;
} }
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out) { void WilsonFermion<Impl>::MooeeInvDag(const FermionField &in,
FermionField &out) {
out.checkerboard = in.checkerboard; out.checkerboard = in.checkerboard;
MooeeInv(in,out); MooeeInv(in, out);
} }
/////////////////////////////////// ///////////////////////////////////
// Internal // Internal
/////////////////////////////////// ///////////////////////////////////
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::DerivInternal(StencilImpl & st, void WilsonFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U,
DoubledGaugeField & U, GaugeField &mat, const FermionField &A,
GaugeField &mat, const FermionField &B, int dag) {
const FermionField &A, assert((dag == DaggerNo) || (dag == DaggerYes));
const FermionField &B,int dag) {
assert((dag==DaggerNo) ||(dag==DaggerYes));
Compressor compressor(dag); Compressor compressor(dag);
@ -156,164 +152,162 @@ namespace QCD {
FermionField Atilde(B._grid); FermionField Atilde(B._grid);
Atilde = A; Atilde = A;
st.HaloExchange(B,compressor); st.HaloExchange(B, compressor);
for(int mu=0;mu<Nd;mu++){
for (int mu = 0; mu < Nd; mu++) {
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
// Flip gamma (1+g)<->(1-g) if dag // Flip gamma (1+g)<->(1-g) if dag
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
int gamma = mu; int gamma = mu;
if ( !dag ) gamma+= Nd; if (!dag) gamma += Nd;
//////////////////////// ////////////////////////
// Call the single hop // Call the single hop
//////////////////////// ////////////////////////
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int sss=0;sss<B._grid->oSites();sss++){ for (int sss = 0; sss < B._grid->oSites(); sss++) {
Kernels::DiracOptDhopDir(st,U,st.comm_buf,sss,sss,B,Btilde,mu,gamma); Kernels::DiracOptDhopDir(st, U, st.comm_buf, sss, sss, B, Btilde, mu,
gamma);
} }
////////////////////////////////////////////////// //////////////////////////////////////////////////
// spin trace outer product // spin trace outer product
////////////////////////////////////////////////// //////////////////////////////////////////////////
Impl::InsertForce4D(mat,Btilde,Atilde,mu); Impl::InsertForce4D(mat, Btilde, Atilde, mu);
}
} }
}
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag) void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U,
{ const FermionField &V, int dag) {
conformable(U._grid,_grid); conformable(U._grid, _grid);
conformable(U._grid,V._grid); conformable(U._grid, V._grid);
conformable(U._grid,mat._grid); conformable(U._grid, mat._grid);
mat.checkerboard = U.checkerboard; mat.checkerboard = U.checkerboard;
DerivInternal(Stencil,Umu,mat,U,V,dag); DerivInternal(Stencil, Umu, mat, U, V, dag);
} }
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag) void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U,
{ const FermionField &V, int dag) {
conformable(U._grid,_cbgrid); conformable(U._grid, _cbgrid);
conformable(U._grid,V._grid); conformable(U._grid, V._grid);
conformable(U._grid,mat._grid); conformable(U._grid, mat._grid);
assert(V.checkerboard==Even); assert(V.checkerboard == Even);
assert(U.checkerboard==Odd); assert(U.checkerboard == Odd);
mat.checkerboard = Odd; mat.checkerboard = Odd;
DerivInternal(StencilEven,UmuOdd,mat,U,V,dag); DerivInternal(StencilEven, UmuOdd, mat, U, V, dag);
} }
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag) void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U,
{ const FermionField &V, int dag) {
conformable(U._grid,_cbgrid); conformable(U._grid, _cbgrid);
conformable(U._grid,V._grid); conformable(U._grid, V._grid);
conformable(U._grid,mat._grid); conformable(U._grid, mat._grid);
assert(V.checkerboard==Odd); assert(V.checkerboard == Odd);
assert(U.checkerboard==Even); assert(U.checkerboard == Even);
mat.checkerboard = Even; mat.checkerboard = Even;
DerivInternal(StencilOdd,UmuEven,mat,U,V,dag); DerivInternal(StencilOdd, UmuEven, mat, U, V, dag);
} }
template <class Impl>
template<class Impl> void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out,
void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out,int dag) { int dag) {
conformable(in._grid,_grid); // verifies full grid conformable(in._grid, _grid); // verifies full grid
conformable(in._grid,out._grid); conformable(in._grid, out._grid);
out.checkerboard = in.checkerboard; out.checkerboard = in.checkerboard;
DhopInternal(Stencil,Lebesgue,Umu,in,out,dag); DhopInternal(Stencil, Lebesgue, Umu, in, out, dag);
} }
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out,int dag) { void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out,
conformable(in._grid,_cbgrid); // verifies half grid int dag) {
conformable(in._grid,out._grid); // drops the cb check conformable(in._grid, _cbgrid); // verifies half grid
conformable(in._grid, out._grid); // drops the cb check
assert(in.checkerboard==Even); assert(in.checkerboard == Even);
out.checkerboard = Odd; out.checkerboard = Odd;
DhopInternal(StencilEven,LebesgueEvenOdd,UmuOdd,in,out,dag); DhopInternal(StencilEven, LebesgueEvenOdd, UmuOdd, in, out, dag);
} }
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopEO(const FermionField &in, FermionField &out,int dag) { void WilsonFermion<Impl>::DhopEO(const FermionField &in, FermionField &out,
conformable(in._grid,_cbgrid); // verifies half grid int dag) {
conformable(in._grid,out._grid); // drops the cb check conformable(in._grid, _cbgrid); // verifies half grid
conformable(in._grid, out._grid); // drops the cb check
assert(in.checkerboard==Odd); assert(in.checkerboard == Odd);
out.checkerboard = Even; out.checkerboard = Even;
DhopInternal(StencilOdd,LebesgueEvenOdd,UmuEven,in,out,dag); DhopInternal(StencilOdd, LebesgueEvenOdd, UmuEven, in, out, dag);
} }
template<class Impl> template <class Impl>
void WilsonFermion<Impl>::Mdir (const FermionField &in, FermionField &out,int dir,int disp) { void WilsonFermion<Impl>::Mdir(const FermionField &in, FermionField &out,
DhopDir(in,out,dir,disp); int dir, int disp) {
} DhopDir(in, out, dir, disp);
}
template <class Impl>
void WilsonFermion<Impl>::DhopDir(const FermionField &in, FermionField &out,
int dir, int disp) {
int skip = (disp == 1) ? 0 : 1;
int dirdisp = dir + skip * 4;
int gamma = dir + (1 - skip) * 4;
template<class Impl> DhopDirDisp(in, out, dirdisp, gamma, DaggerNo);
void WilsonFermion<Impl>::DhopDir(const FermionField &in, FermionField &out,int dir,int disp){ };
int skip = (disp==1) ? 0 : 1;
int dirdisp = dir+skip*4;
int gamma = dir+(1-skip)*4;
DhopDirDisp(in,out,dirdisp,gamma,DaggerNo);
};
template<class Impl>
void WilsonFermion<Impl>::DhopDirDisp(const FermionField &in, FermionField &out,int dirdisp,int gamma,int dag) {
template <class Impl>
void WilsonFermion<Impl>::DhopDirDisp(const FermionField &in, FermionField &out,
int dirdisp, int gamma, int dag) {
Compressor compressor(dag); Compressor compressor(dag);
Stencil.HaloExchange(in,compressor); Stencil.HaloExchange(in, compressor);
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int sss=0;sss<in._grid->oSites();sss++){ for (int sss = 0; sss < in._grid->oSites(); sss++) {
Kernels::DiracOptDhopDir(Stencil,Umu,Stencil.comm_buf,sss,sss,in,out,dirdisp,gamma); Kernels::DiracOptDhopDir(Stencil, Umu, Stencil.comm_buf, sss, sss, in, out,
dirdisp, gamma);
} }
};
}; template <class Impl>
void WilsonFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder &lo,
template<class Impl> DoubledGaugeField &U,
void WilsonFermion<Impl>::DhopInternal(StencilImpl & st,LebesgueOrder& lo,DoubledGaugeField & U, const FermionField &in,
const FermionField &in, FermionField &out,int dag) FermionField &out, int dag) {
{ assert((dag == DaggerNo) || (dag == DaggerYes));
assert((dag==DaggerNo) ||(dag==DaggerYes));
Compressor compressor(dag); Compressor compressor(dag);
st.HaloExchange(in,compressor); st.HaloExchange(in, compressor);
if ( dag == DaggerYes ) { if (dag == DaggerYes) {
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int sss=0;sss<in._grid->oSites();sss++){ for (int sss = 0; sss < in._grid->oSites(); sss++) {
Kernels::DiracOptDhopSiteDag(st,lo,U,st.comm_buf,sss,sss,1,1,in,out); Kernels::DiracOptDhopSiteDag(st, lo, U, st.comm_buf, sss, sss, 1, 1, in,
out);
} }
} else { } else {
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int sss=0;sss<in._grid->oSites();sss++){ for (int sss = 0; sss < in._grid->oSites(); sss++) {
Kernels::DiracOptDhopSite(st,lo,U,st.comm_buf,sss,sss,1,1,in,out); Kernels::DiracOptDhopSite(st, lo, U, st.comm_buf, sss, sss, 1, 1, in,
out);
} }
} }
}; };
FermOpTemplateInstantiate(WilsonFermion);
GparityFermOpTemplateInstantiate(WilsonFermion);
}}
FermOpTemplateInstantiate(WilsonFermion);
GparityFermOpTemplateInstantiate(WilsonFermion);
}
}

134
lib/qcd/action/fermion/WilsonFermion.h
View File

@ -1,51 +1,51 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/fermion/WilsonFermion.h Source file: ./lib/qcd/action/fermion/WilsonFermion.h
Copyright (C) 2015 Copyright (C) 2015
Author: Peter Boyle <pabobyle@ph.ed.ac.uk> Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk> Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along 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., with this program; if not, write to the Free Software Foundation, Inc.,
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 GRID_QCD_WILSON_FERMION_H #ifndef GRID_QCD_WILSON_FERMION_H
#define GRID_QCD_WILSON_FERMION_H #define GRID_QCD_WILSON_FERMION_H
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
class WilsonFermionStatic { class WilsonFermionStatic {
public: public:
static int HandOptDslash; // these are a temporary hack static int HandOptDslash; // these are a temporary hack
static int MortonOrder; static int MortonOrder;
static const std::vector<int> directions ; static const std::vector<int> directions;
static const std::vector<int> displacements; static const std::vector<int> displacements;
static const int npoint=8; static const int npoint = 8;
}; };
template<class Impl> template <class Impl>
class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
{
public: public:
INHERIT_IMPL_TYPES(Impl); INHERIT_IMPL_TYPES(Impl);
typedef WilsonKernels<Impl> Kernels; typedef WilsonKernels<Impl> Kernels;
@ -53,10 +53,10 @@ namespace Grid {
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
// Implement the abstract base // Implement the abstract base
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
GridBase *GaugeGrid(void) { return _grid ;} GridBase *GaugeGrid(void) { return _grid; }
GridBase *GaugeRedBlackGrid(void) { return _cbgrid ;} GridBase *GaugeRedBlackGrid(void) { return _cbgrid; }
GridBase *FermionGrid(void) { return _grid;} GridBase *FermionGrid(void) { return _grid; }
GridBase *FermionRedBlackGrid(void) { return _cbgrid;} GridBase *FermionRedBlackGrid(void) { return _cbgrid; }
////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////
// override multiply; cut number routines if pass dagger argument // override multiply; cut number routines if pass dagger argument
@ -69,58 +69,54 @@ namespace Grid {
// half checkerboard operations // half checkerboard operations
// could remain virtual so we can derive Clover from Wilson base // could remain virtual so we can derive Clover from Wilson base
///////////////////////////////////////////////////////// /////////////////////////////////////////////////////////
void Meooe(const FermionField &in, FermionField &out) ; void Meooe(const FermionField &in, FermionField &out);
void MeooeDag(const FermionField &in, FermionField &out) ; void MeooeDag(const FermionField &in, FermionField &out);
// allow override for twisted mass and clover // allow override for twisted mass and clover
virtual void Mooee(const FermionField &in, FermionField &out) ; virtual void Mooee(const FermionField &in, FermionField &out);
virtual void MooeeDag(const FermionField &in, FermionField &out) ; virtual void MooeeDag(const FermionField &in, FermionField &out);
virtual void MooeeInv(const FermionField &in, FermionField &out) ; virtual void MooeeInv(const FermionField &in, FermionField &out);
virtual void MooeeInvDag(const FermionField &in, FermionField &out) ; virtual void MooeeInvDag(const FermionField &in, FermionField &out);
//////////////////////// ////////////////////////
// Derivative interface // Derivative interface
//////////////////////// ////////////////////////
// Interface calls an internal routine // Interface calls an internal routine
void DhopDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); void DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V,
void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); int dag);
void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); void DhopDerivOE(GaugeField &mat, const FermionField &U,
const FermionField &V, int dag);
void DhopDerivEO(GaugeField &mat, const FermionField &U,
const FermionField &V, int dag);
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
// non-hermitian hopping term; half cb or both // non-hermitian hopping term; half cb or both
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
void Dhop(const FermionField &in, FermionField &out,int dag) ; void Dhop(const FermionField &in, FermionField &out, int dag);
void DhopOE(const FermionField &in, FermionField &out,int dag) ; void DhopOE(const FermionField &in, FermionField &out, int dag);
void DhopEO(const FermionField &in, FermionField &out,int dag) ; void DhopEO(const FermionField &in, FermionField &out, int dag);
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
// Multigrid assistance; force term uses too // Multigrid assistance; force term uses too
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
void Mdir (const FermionField &in, FermionField &out,int dir,int disp) ; void Mdir(const FermionField &in, FermionField &out, int dir, int disp);
void DhopDir(const FermionField &in, FermionField &out,int dir,int disp); void DhopDir(const FermionField &in, FermionField &out, int dir, int disp);
void DhopDirDisp(const FermionField &in, FermionField &out,int dirdisp,int gamma,int dag) ; void DhopDirDisp(const FermionField &in, FermionField &out, int dirdisp,
int gamma, int dag);
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
// Extra methods added by derived // Extra methods added by derived
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
void DerivInternal(StencilImpl & st, void DerivInternal(StencilImpl &st, DoubledGaugeField &U, GaugeField &mat,
DoubledGaugeField & U, const FermionField &A, const FermionField &B, int dag);
GaugeField &mat,
const FermionField &A,
const FermionField &B,
int dag);
void DhopInternal(StencilImpl & st,LebesgueOrder & lo,DoubledGaugeField & U, void DhopInternal(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
const FermionField &in, FermionField &out,int dag) ; const FermionField &in, FermionField &out, int dag);
// Constructor // Constructor
WilsonFermion(GaugeField &_Umu, WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, RealD _mass,
GridRedBlackCartesian &Hgrid, const ImplParams &p = ImplParams());
RealD _mass,
const ImplParams &p= ImplParams()
) ;
// DoubleStore impl dependent // DoubleStore impl dependent
void ImportGauge(const GaugeField &_Umu); void ImportGauge(const GaugeField &_Umu);
@ -131,13 +127,12 @@ namespace Grid {
// protected: // protected:
public: public:
RealD mass; RealD mass;
GridBase * _grid; GridBase *_grid;
GridBase * _cbgrid; GridBase *_cbgrid;
//Defines the stencils for even and odd // Defines the stencils for even and odd
StencilImpl Stencil; StencilImpl Stencil;
StencilImpl StencilEven; StencilImpl StencilEven;
StencilImpl StencilOdd; StencilImpl StencilOdd;
@ -149,13 +144,10 @@ namespace Grid {
LebesgueOrder Lebesgue; LebesgueOrder Lebesgue;
LebesgueOrder LebesgueEvenOdd; LebesgueOrder LebesgueEvenOdd;
};
typedef WilsonFermion<WilsonImplF> WilsonFermionF;
}; typedef WilsonFermion<WilsonImplD> WilsonFermionD;
}
typedef WilsonFermion<WilsonImplF> WilsonFermionF;
typedef WilsonFermion<WilsonImplD> WilsonFermionD;
}
} }
#endif #endif

138
lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h
View File

@ -1,52 +1,49 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h Source file: ./lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h
Copyright (C) 2015 Copyright (C) 2015
Author: Peter Boyle <pabobyle@ph.ed.ac.uk> Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along 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., with this program; if not, write to the Free Software Foundation, Inc.,
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_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H #ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H
#define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H #define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H
namespace Grid{ namespace Grid {
namespace QCD{ namespace QCD {
////////////////////////////////////////////////////////////////////////
// Two flavour pseudofermion action for any EO prec dop
////////////////////////////////////////////////////////////////////////
template<class Impl>
class TwoFlavourEvenOddPseudoFermionAction : public Action<typename Impl::GaugeField> {
////////////////////////////////////////////////////////////////////////
// Two flavour pseudofermion action for any EO prec dop
////////////////////////////////////////////////////////////////////////
template <class Impl>
class TwoFlavourEvenOddPseudoFermionAction
: public Action<typename Impl::GaugeField> {
public: public:
INHERIT_IMPL_TYPES(Impl); INHERIT_IMPL_TYPES(Impl);
private: private:
FermionOperator<Impl> &FermOp; // the basic operator
FermionOperator<Impl> & FermOp;// the basic operator
OperatorFunction<FermionField> &DerivativeSolver; OperatorFunction<FermionField> &DerivativeSolver;
OperatorFunction<FermionField> &ActionSolver; OperatorFunction<FermionField> &ActionSolver;
@ -59,21 +56,19 @@ namespace Grid{
// Pass in required objects. // Pass in required objects.
///////////////////////////////////////////////// /////////////////////////////////////////////////
TwoFlavourEvenOddPseudoFermionAction(FermionOperator<Impl> &Op, TwoFlavourEvenOddPseudoFermionAction(FermionOperator<Impl> &Op,
OperatorFunction<FermionField> & DS, OperatorFunction<FermionField> &DS,
OperatorFunction<FermionField> & AS OperatorFunction<FermionField> &AS)
) : : FermOp(Op),
FermOp(Op),
DerivativeSolver(DS), DerivativeSolver(DS),
ActionSolver(AS), ActionSolver(AS),
PhiEven(Op.FermionRedBlackGrid()), PhiEven(Op.FermionRedBlackGrid()),
PhiOdd(Op.FermionRedBlackGrid()) PhiOdd(Op.FermionRedBlackGrid()){};
{};
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
// Push the gauge field in to the dops. Assume any BC's and smearing already applied // Push the gauge field in to the dops. Assume any BC's and smearing already
// applied
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) { virtual void refresh(const GaugeField &U, GridParallelRNG &pRNG) {
// P(phi) = e^{- phi^dag (MpcdagMpc)^-1 phi} // P(phi) = e^{- phi^dag (MpcdagMpc)^-1 phi}
// Phi = McpDag eta // Phi = McpDag eta
// P(eta) = e^{- eta^dag eta} // P(eta) = e^{- eta^dag eta}
@ -82,25 +77,23 @@ namespace Grid{
RealD scale = std::sqrt(0.5); RealD scale = std::sqrt(0.5);
FermionField eta (FermOp.FermionGrid()); FermionField eta(FermOp.FermionGrid());
FermionField etaOdd (FermOp.FermionRedBlackGrid()); FermionField etaOdd(FermOp.FermionRedBlackGrid());
FermionField etaEven(FermOp.FermionRedBlackGrid()); FermionField etaEven(FermOp.FermionRedBlackGrid());
gaussian(pRNG,eta); gaussian(pRNG, eta);
pickCheckerboard(Even,etaEven,eta); pickCheckerboard(Even, etaEven, eta);
pickCheckerboard(Odd,etaOdd,eta); pickCheckerboard(Odd, etaOdd, eta);
FermOp.ImportGauge(U); FermOp.ImportGauge(U);
SchurDifferentiableOperator<Impl> PCop(FermOp); SchurDifferentiableOperator<Impl> PCop(FermOp);
PCop.MpcDag(etaOdd, PhiOdd);
PCop.MpcDag(etaOdd,PhiOdd); FermOp.MooeeDag(etaEven, PhiEven);
FermOp.MooeeDag(etaEven,PhiEven);
PhiOdd =PhiOdd*scale;
PhiEven=PhiEven*scale;
PhiOdd = PhiOdd * scale;
PhiEven = PhiEven * scale;
}; };
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
@ -108,7 +101,6 @@ namespace Grid{
// + phi^dag (Mdag M)^-1 phi (even) // + phi^dag (Mdag M)^-1 phi (even)
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) { virtual RealD S(const GaugeField &U) {
FermOp.ImportGauge(U); FermOp.ImportGauge(U);
FermionField X(FermOp.FermionRedBlackGrid()); FermionField X(FermOp.FermionRedBlackGrid());
@ -116,30 +108,32 @@ namespace Grid{
SchurDifferentiableOperator<Impl> PCop(FermOp); SchurDifferentiableOperator<Impl> PCop(FermOp);
X=zero; X = zero;
ActionSolver(PCop,PhiOdd,X); ActionSolver(PCop, PhiOdd, X);
PCop.Op(X,Y); PCop.Op(X, Y);
RealD action = norm2(Y); RealD action = norm2(Y);
// The EE factorised block; normally can replace with zero if det is constant (gauge field indept) // The EE factorised block; normally can replace with zero if det is
// constant (gauge field indept)
// Only really clover term that creates this. // Only really clover term that creates this.
FermOp.MooeeInvDag(PhiEven,Y); FermOp.MooeeInvDag(PhiEven, Y);
action = action + norm2(Y); action = action + norm2(Y);
std::cout << GridLogMessage << "Pseudofermion EO action "<<action<<std::endl; std::cout << GridLogMessage << "Pseudofermion EO action " << action
<< std::endl;
return action; return action;
}; };
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
// //
// dS/du = - phi^dag (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 phi // dS/du = - phi^dag (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 phi
// = - phi^dag M^-1 dM (MdagM)^-1 phi - phi^dag (MdagM)^-1 dMdag dM (Mdag)^-1 phi // = - phi^dag M^-1 dM (MdagM)^-1 phi - phi^dag (MdagM)^-1 dMdag dM
// (Mdag)^-1 phi
// //
// = - Ydag dM X - Xdag dMdag Y // = - Ydag dM X - Xdag dMdag Y
// //
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U,GaugeField & dSdU) { virtual void deriv(const GaugeField &U, GaugeField &dSdU) {
FermOp.ImportGauge(U); FermOp.ImportGauge(U);
FermionField X(FermOp.FermionRedBlackGrid()); FermionField X(FermOp.FermionRedBlackGrid());
@ -148,14 +142,17 @@ namespace Grid{
SchurDifferentiableOperator<Impl> Mpc(FermOp); SchurDifferentiableOperator<Impl> Mpc(FermOp);
// Our conventions really make this UdSdU; We do not differentiate wrt Udag here. // Our conventions really make this UdSdU; We do not differentiate wrt Udag
// here.
// So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt. // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt.
X=zero; X = zero;
DerivativeSolver(Mpc,PhiOdd,X); DerivativeSolver(Mpc, PhiOdd, X);
Mpc.Mpc(X,Y); Mpc.Mpc(X, Y);
Mpc.MpcDeriv(tmp , Y, X ); dSdU=tmp; Mpc.MpcDeriv(tmp, Y, X);
Mpc.MpcDagDeriv(tmp , X, Y); dSdU=dSdU+tmp; dSdU = tmp;
Mpc.MpcDagDeriv(tmp, X, Y);
dSdU = dSdU + tmp;
// Treat the EE case. (MdagM)^-1 = Minv Minvdag // Treat the EE case. (MdagM)^-1 = Minv Minvdag
// Deriv defaults to zero. // Deriv defaults to zero.
@ -174,12 +171,9 @@ namespace Grid{
*/ */
dSdU = Ta(dSdU); dSdU = Ta(dSdU);
}; };
};
}; }
}
} }
#endif #endif