mirror of
https://github.com/paboyle/Grid.git
synced 2025-06-13 04:37:05 +01:00
Merge branch 'develop' into feature/feynman-rules
# Conflicts: # lib/Threads.h # lib/qcd/action/fermion/WilsonFermion.cc # lib/qcd/action/fermion/WilsonFermion.h # lib/qcd/utils/SUn.h # lib/simd/Grid_avx.h # lib/simd/Intel512common.h
This commit is contained in:
@ -1,87 +1,153 @@
|
||||
/*************************************************************************************
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/qcd/action/ActionBase.h
|
||||
Source file: ./lib/qcd/action/ActionBase.h
|
||||
|
||||
Copyright (C) 2015
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: neo <cossu@post.kek.jp>
|
||||
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
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.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
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.,
|
||||
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 */
|
||||
See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#ifndef QCD_ACTION_BASE
|
||||
#define QCD_ACTION_BASE
|
||||
namespace Grid {
|
||||
namespace QCD{
|
||||
|
||||
template<class GaugeField>
|
||||
class Action {
|
||||
namespace QCD {
|
||||
|
||||
template <class GaugeField>
|
||||
class Action {
|
||||
public:
|
||||
bool is_smeared = false;
|
||||
// Boundary conditions? // Heatbath?
|
||||
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) = 0;// refresh pseudofermions
|
||||
virtual RealD S (const GaugeField &U) = 0; // evaluate the action
|
||||
virtual void deriv(const GaugeField &U,GaugeField & dSdU ) = 0; // evaluate the action derivative
|
||||
virtual ~Action() {};
|
||||
virtual void refresh(const GaugeField& U,
|
||||
GridParallelRNG& pRNG) = 0; // refresh pseudofermions
|
||||
virtual RealD S(const GaugeField& U) = 0; // evaluate the action
|
||||
virtual void deriv(const GaugeField& U,
|
||||
GaugeField& dSdU) = 0; // evaluate the action derivative
|
||||
virtual ~Action(){};
|
||||
};
|
||||
|
||||
// Indexing of tuple types
|
||||
template <class T, class Tuple>
|
||||
struct Index;
|
||||
|
||||
template <class T, class... Types>
|
||||
struct Index<T, std::tuple<T, Types...>> {
|
||||
static const std::size_t value = 0;
|
||||
};
|
||||
|
||||
template <class T, class U, class... Types>
|
||||
struct Index<T, std::tuple<U, Types...>> {
|
||||
static const std::size_t value = 1 + Index<T, std::tuple<Types...>>::value;
|
||||
};
|
||||
|
||||
// Could derive PseudoFermion action with a PF field, FermionField, and a Grid; implement refresh
|
||||
/*
|
||||
template<class GaugeField, class FermionField>
|
||||
class PseudoFermionAction : public Action<GaugeField> {
|
||||
template <class GaugeField>
|
||||
struct ActionLevel {
|
||||
public:
|
||||
FermionField Phi;
|
||||
GridParallelRNG &pRNG;
|
||||
GridBase &Grid;
|
||||
typedef Action<GaugeField>*
|
||||
ActPtr; // now force the same colours as the rest of the code
|
||||
|
||||
PseudoFermionAction(GridBase &_Grid,GridParallelRNG &_pRNG) : Grid(_Grid), Phi(&_Grid), pRNG(_pRNG) {
|
||||
};
|
||||
//Add supported representations here
|
||||
|
||||
virtual void refresh(const GaugeField &gauge) {
|
||||
gaussian(Phi,pRNG);
|
||||
};
|
||||
|
||||
};
|
||||
*/
|
||||
|
||||
template<class GaugeField> struct ActionLevel{
|
||||
public:
|
||||
|
||||
typedef Action<GaugeField>* ActPtr; // now force the same colours as the rest of the code
|
||||
|
||||
int multiplier;
|
||||
unsigned int multiplier;
|
||||
|
||||
std::vector<ActPtr> actions;
|
||||
|
||||
ActionLevel(int mul = 1) : multiplier(mul) {
|
||||
assert (mul > 0);
|
||||
ActionLevel(unsigned int mul = 1) : actions(0), multiplier(mul) {
|
||||
assert(mul >= 1);
|
||||
};
|
||||
|
||||
void push_back(ActPtr ptr){
|
||||
actions.push_back(ptr);
|
||||
|
||||
void push_back(ActPtr ptr) { actions.push_back(ptr); }
|
||||
};
|
||||
*/
|
||||
|
||||
template <class GaugeField, class Repr = NoHirep >
|
||||
struct ActionLevel {
|
||||
public:
|
||||
unsigned int multiplier;
|
||||
|
||||
// Fundamental repr actions separated because of the smearing
|
||||
typedef Action<GaugeField>* ActPtr;
|
||||
|
||||
// construct a tuple of vectors of the actions for the corresponding higher
|
||||
// representation fields
|
||||
typedef typename AccessTypes<Action, Repr>::VectorCollection action_collection;
|
||||
action_collection actions_hirep;
|
||||
typedef typename AccessTypes<Action, Repr>::FieldTypeCollection action_hirep_types;
|
||||
|
||||
std::vector<ActPtr>& actions;
|
||||
|
||||
// Temporary conversion between ActionLevel and ActionLevelHirep
|
||||
//ActionLevelHirep(ActionLevel<GaugeField>& AL ):actions(AL.actions), multiplier(AL.multiplier){}
|
||||
|
||||
ActionLevel(unsigned int mul = 1) : actions(std::get<0>(actions_hirep)), multiplier(mul) {
|
||||
// initialize the hirep vectors to zero.
|
||||
//apply(this->resize, actions_hirep, 0); //need a working resize
|
||||
assert(mul >= 1);
|
||||
};
|
||||
|
||||
//void push_back(ActPtr ptr) { actions.push_back(ptr); }
|
||||
|
||||
|
||||
|
||||
template < class Field >
|
||||
void push_back(Action<Field>* ptr) {
|
||||
// insert only in the correct vector
|
||||
std::get< Index < Field, action_hirep_types>::value >(actions_hirep).push_back(ptr);
|
||||
};
|
||||
|
||||
|
||||
|
||||
template < class ActPtr>
|
||||
static void resize(ActPtr ap, unsigned int n){
|
||||
ap->resize(n);
|
||||
|
||||
}
|
||||
|
||||
//template <std::size_t I>
|
||||
//auto getRepresentation(Repr& R)->decltype(std::get<I>(R).U) {return std::get<I>(R).U;}
|
||||
|
||||
// Loop on tuple for a callable function
|
||||
template <std::size_t I = 1, typename Callable, typename ...Args>
|
||||
inline typename std::enable_if<I == std::tuple_size<action_collection>::value, void>::type apply(
|
||||
Callable, Repr& R,Args&...) const {}
|
||||
|
||||
template <std::size_t I = 1, typename Callable, typename ...Args>
|
||||
inline typename std::enable_if<I < std::tuple_size<action_collection>::value, void>::type apply(
|
||||
Callable fn, Repr& R, Args&... arguments) const {
|
||||
fn(std::get<I>(actions_hirep), std::get<I>(R.rep), arguments...);
|
||||
apply<I + 1>(fn, R, arguments...);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
template<class GaugeField> using ActionSet = std::vector<ActionLevel< GaugeField > >;
|
||||
|
||||
//template <class GaugeField>
|
||||
//using ActionSet = std::vector<ActionLevel<GaugeField> >;
|
||||
|
||||
}}
|
||||
template <class GaugeField, class R>
|
||||
using ActionSet = std::vector<ActionLevel<GaugeField, R> >;
|
||||
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
@ -116,6 +116,14 @@ typedef SymanzikGaugeAction<ConjugateGimplD> ConjugateSymanzikGaugeAction
|
||||
template class A<GparityWilsonImplF>; \
|
||||
template class A<GparityWilsonImplD>;
|
||||
|
||||
#define AdjointFermOpTemplateInstantiate(A) \
|
||||
template class A<WilsonAdjImplF>; \
|
||||
template class A<WilsonAdjImplD>;
|
||||
|
||||
#define TwoIndexFermOpTemplateInstantiate(A) \
|
||||
template class A<WilsonTwoIndexSymmetricImplF>; \
|
||||
template class A<WilsonTwoIndexSymmetricImplD>;
|
||||
|
||||
#define FermOp5dVecTemplateInstantiate(A) \
|
||||
template class A<DomainWallVec5dImplF>; \
|
||||
template class A<DomainWallVec5dImplD>; \
|
||||
@ -126,6 +134,7 @@ typedef SymanzikGaugeAction<ConjugateGimplD> ConjugateSymanzikGaugeAction
|
||||
FermOp4dVecTemplateInstantiate(A) \
|
||||
FermOp5dVecTemplateInstantiate(A)
|
||||
|
||||
|
||||
#define GparityFermOpTemplateInstantiate(A)
|
||||
|
||||
////////////////////////////////////////////
|
||||
@ -171,6 +180,14 @@ typedef WilsonFermion<WilsonImplR> WilsonFermionR;
|
||||
typedef WilsonFermion<WilsonImplF> WilsonFermionF;
|
||||
typedef WilsonFermion<WilsonImplD> WilsonFermionD;
|
||||
|
||||
typedef WilsonFermion<WilsonAdjImplR> WilsonAdjFermionR;
|
||||
typedef WilsonFermion<WilsonAdjImplF> WilsonAdjFermionF;
|
||||
typedef WilsonFermion<WilsonAdjImplD> WilsonAdjFermionD;
|
||||
|
||||
typedef WilsonFermion<WilsonTwoIndexSymmetricImplR> WilsonTwoIndexSymmetricFermionR;
|
||||
typedef WilsonFermion<WilsonTwoIndexSymmetricImplF> WilsonTwoIndexSymmetricFermionF;
|
||||
typedef WilsonFermion<WilsonTwoIndexSymmetricImplD> WilsonTwoIndexSymmetricFermionD;
|
||||
|
||||
typedef WilsonTMFermion<WilsonImplR> WilsonTMFermionR;
|
||||
typedef WilsonTMFermion<WilsonImplF> WilsonTMFermionF;
|
||||
typedef WilsonTMFermion<WilsonImplD> WilsonTMFermionD;
|
||||
|
||||
@ -1,35 +1,36 @@
|
||||
/*************************************************************************************
|
||||
/*************************************************************************************
|
||||
|
||||
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 <paboyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
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.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
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.,
|
||||
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 */
|
||||
#ifndef GRID_QCD_FERMION_OPERATOR_IMPL_H
|
||||
#define GRID_QCD_FERMION_OPERATOR_IMPL_H
|
||||
See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#ifndef GRID_QCD_FERMION_OPERATOR_IMPL_H
|
||||
#define GRID_QCD_FERMION_OPERATOR_IMPL_H
|
||||
|
||||
namespace Grid {
|
||||
|
||||
@ -99,258 +100,281 @@ namespace Grid {
|
||||
typedef typename Impl::SiteSpinor SiteSpinor; \
|
||||
typedef typename Impl::SiteHalfSpinor SiteHalfSpinor; \
|
||||
typedef typename Impl::Compressor Compressor; \
|
||||
typedef typename Impl::StencilImpl StencilImpl; \
|
||||
typedef typename Impl::ImplParams ImplParams; \
|
||||
typedef typename Impl::StencilImpl StencilImpl; \
|
||||
typedef typename Impl::ImplParams ImplParams; \
|
||||
typedef typename Impl::Coeff_t Coeff_t;
|
||||
|
||||
#define INHERIT_IMPL_TYPES(Base) \
|
||||
INHERIT_GIMPL_TYPES(Base)\
|
||||
INHERIT_GIMPL_TYPES(Base) \
|
||||
INHERIT_FIMPL_TYPES(Base)
|
||||
|
||||
|
||||
///////
|
||||
// Single flavour four spinors with colour index
|
||||
///////
|
||||
template<class S,int Nrepresentation=Nc,class _Coeff_t = RealD>
|
||||
class WilsonImpl : public PeriodicGaugeImpl< GaugeImplTypes< S, Nrepresentation> > {
|
||||
template <class S, class Representation = FundamentalRepresentation,class _Coeff_t = RealD >
|
||||
class WilsonImpl
|
||||
: public PeriodicGaugeImpl<GaugeImplTypes<S, Representation::Dimension > > {
|
||||
public:
|
||||
|
||||
static const int Dimension = Representation::Dimension;
|
||||
typedef PeriodicGaugeImpl<GaugeImplTypes<S, Dimension > > Gimpl;
|
||||
|
||||
//Necessary?
|
||||
constexpr bool is_fundamental() const{return Dimension == Nc ? 1 : 0;}
|
||||
|
||||
const bool LsVectorised=false;
|
||||
|
||||
typedef _Coeff_t Coeff_t;
|
||||
typedef PeriodicGaugeImpl< GaugeImplTypes< S,Nrepresentation> > Gimpl;
|
||||
|
||||
|
||||
INHERIT_GIMPL_TYPES(Gimpl);
|
||||
|
||||
template<typename vtype> using iImplSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
|
||||
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 iImplHalfSpinor<Simd> SiteHalfSpinor;
|
||||
typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
|
||||
|
||||
typedef Lattice<SiteSpinor> FermionField;
|
||||
|
||||
template <typename vtype> using iImplSpinor = iScalar<iVector<iVector<vtype, Dimension>, Ns> >;
|
||||
template <typename vtype> using iImplHalfSpinor = iScalar<iVector<iVector<vtype, Dimension>, Nhs> >;
|
||||
template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>;
|
||||
|
||||
typedef iImplSpinor<Simd> SiteSpinor;
|
||||
typedef iImplHalfSpinor<Simd> SiteHalfSpinor;
|
||||
typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
|
||||
|
||||
typedef Lattice<SiteSpinor> FermionField;
|
||||
typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
|
||||
|
||||
typedef WilsonCompressor<SiteHalfSpinor,SiteSpinor> Compressor;
|
||||
|
||||
typedef WilsonCompressor<SiteHalfSpinor, SiteSpinor> Compressor;
|
||||
typedef WilsonImplParams ImplParams;
|
||||
typedef WilsonStencil<SiteSpinor,SiteHalfSpinor> StencilImpl;
|
||||
|
||||
typedef WilsonStencil<SiteSpinor, SiteHalfSpinor> StencilImpl;
|
||||
|
||||
ImplParams Params;
|
||||
|
||||
WilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {};
|
||||
|
||||
|
||||
WilsonImpl(const ImplParams &p = ImplParams()) : Params(p){};
|
||||
|
||||
bool overlapCommsCompute(void) { return Params.overlapCommsCompute; };
|
||||
|
||||
inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,StencilImpl &St){
|
||||
mult(&phi(),&U(mu),&chi());
|
||||
}
|
||||
|
||||
template<class ref>
|
||||
inline void loadLinkElement(Simd & reg,ref &memory){
|
||||
reg = memory;
|
||||
}
|
||||
inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
|
||||
{
|
||||
conformable(Uds._grid,GaugeGrid);
|
||||
conformable(Umu._grid,GaugeGrid);
|
||||
GaugeLinkField U(GaugeGrid);
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U = PeekIndex<LorentzIndex>(Umu,mu);
|
||||
PokeIndex<LorentzIndex>(Uds,U,mu);
|
||||
U = adj(Cshift(U,mu,-1));
|
||||
PokeIndex<LorentzIndex>(Uds,U,mu+4);
|
||||
}
|
||||
|
||||
inline void multLink(SiteHalfSpinor &phi,
|
||||
const SiteDoubledGaugeField &U,
|
||||
const SiteHalfSpinor &chi,
|
||||
int mu,
|
||||
StencilEntry *SE,
|
||||
StencilImpl &St) {
|
||||
mult(&phi(), &U(mu), &chi());
|
||||
}
|
||||
|
||||
template <class ref>
|
||||
inline void loadLinkElement(Simd ®,
|
||||
ref &memory) {
|
||||
reg = memory;
|
||||
}
|
||||
|
||||
inline void DoubleStore(GridBase *GaugeGrid,
|
||||
DoubledGaugeField &Uds,
|
||||
const GaugeField &Umu) {
|
||||
conformable(Uds._grid, GaugeGrid);
|
||||
conformable(Umu._grid, GaugeGrid);
|
||||
GaugeLinkField U(GaugeGrid);
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
U = PeekIndex<LorentzIndex>(Umu, mu);
|
||||
PokeIndex<LorentzIndex>(Uds, U, mu);
|
||||
U = adj(Cshift(U, mu, -1));
|
||||
PokeIndex<LorentzIndex>(Uds, U, mu + 4);
|
||||
}
|
||||
}
|
||||
|
||||
inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){
|
||||
GaugeLinkField link(mat._grid);
|
||||
link = TraceIndex<SpinIndex>(outerProduct(Btilde,A));
|
||||
PokeIndex<LorentzIndex>(mat,link,mu);
|
||||
}
|
||||
|
||||
|
||||
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){
|
||||
|
||||
|
||||
int Ls=Btilde._grid->_fdimensions[0];
|
||||
|
||||
GaugeLinkField tmp(mat._grid);
|
||||
tmp = zero;
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int sss=0;sss<tmp._grid->oSites();sss++){
|
||||
int sU=sss;
|
||||
for(int s=0;s<Ls;s++){
|
||||
int sF = s+Ls*sU;
|
||||
tmp[sU] = tmp[sU]+ traceIndex<SpinIndex>(outerProduct(Btilde[sF],Atilde[sF])); // ordering here
|
||||
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int sss=0;sss<tmp._grid->oSites();sss++){
|
||||
int sU=sss;
|
||||
for(int s=0;s<Ls;s++){
|
||||
int sF = s+Ls*sU;
|
||||
tmp[sU] = tmp[sU]+ traceIndex<SpinIndex>(outerProduct(Btilde[sF],Atilde[sF])); // ordering here
|
||||
}
|
||||
}
|
||||
}
|
||||
PokeIndex<LorentzIndex>(mat,tmp,mu);
|
||||
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
///////
|
||||
// Single flavour four spinors with colour index, 5d redblack
|
||||
///////
|
||||
template<class S,int Nrepresentation=Nc,class _Coeff_t = RealD>
|
||||
class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepresentation> > {
|
||||
public:
|
||||
|
||||
|
||||
static const int Dimension = Nrepresentation;
|
||||
const bool LsVectorised=true;
|
||||
|
||||
typedef _Coeff_t Coeff_t;
|
||||
typedef PeriodicGaugeImpl< GaugeImplTypes< S,Nrepresentation> > Gimpl;
|
||||
typedef _Coeff_t Coeff_t;
|
||||
typedef PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > Gimpl;
|
||||
|
||||
INHERIT_GIMPL_TYPES(Gimpl);
|
||||
|
||||
template<typename vtype> using iImplSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
|
||||
template<typename vtype> using iImplHalfSpinor = 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 iImplHalfSpinor<Simd> SiteHalfSpinor;
|
||||
typedef Lattice<SiteSpinor> FermionField;
|
||||
|
||||
template <typename vtype> using iImplSpinor = iScalar<iVector<iVector<vtype, Nrepresentation>, Ns> >;
|
||||
template <typename vtype> using iImplHalfSpinor = 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 iImplHalfSpinor<Simd> SiteHalfSpinor;
|
||||
typedef Lattice<SiteSpinor> FermionField;
|
||||
|
||||
// Make the doubled gauge field a *scalar*
|
||||
typedef iImplDoubledGaugeField<typename Simd::scalar_type> SiteDoubledGaugeField; // This is a scalar
|
||||
typedef iImplGaugeField<typename Simd::scalar_type> SiteScalarGaugeField; // scalar
|
||||
typedef iImplGaugeLink <typename Simd::scalar_type> SiteScalarGaugeLink; // scalar
|
||||
|
||||
typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
|
||||
|
||||
typedef WilsonCompressor<SiteHalfSpinor,SiteSpinor> Compressor;
|
||||
typedef iImplDoubledGaugeField<typename Simd::scalar_type>
|
||||
SiteDoubledGaugeField; // This is a scalar
|
||||
typedef iImplGaugeField<typename Simd::scalar_type>
|
||||
SiteScalarGaugeField; // scalar
|
||||
typedef iImplGaugeLink<typename Simd::scalar_type>
|
||||
SiteScalarGaugeLink; // scalar
|
||||
|
||||
typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
|
||||
|
||||
typedef WilsonCompressor<SiteHalfSpinor, SiteSpinor> Compressor;
|
||||
typedef WilsonImplParams ImplParams;
|
||||
typedef WilsonStencil<SiteSpinor,SiteHalfSpinor> StencilImpl;
|
||||
|
||||
typedef WilsonStencil<SiteSpinor, SiteHalfSpinor> StencilImpl;
|
||||
|
||||
ImplParams Params;
|
||||
|
||||
DomainWallVec5dImpl(const ImplParams &p= ImplParams()) : Params(p) {};
|
||||
|
||||
|
||||
DomainWallVec5dImpl(const ImplParams &p = ImplParams()) : Params(p){};
|
||||
|
||||
bool overlapCommsCompute(void) { return false; };
|
||||
|
||||
template<class ref>
|
||||
inline void loadLinkElement(Simd & reg,ref &memory){
|
||||
vsplat(reg,memory);
|
||||
|
||||
template <class ref>
|
||||
inline void loadLinkElement(Simd ®, ref &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;
|
||||
for(int i=0;i<Nrepresentation;i++){
|
||||
for(int j=0;j<Nrepresentation;j++){
|
||||
vsplat(UU()()(i,j),U(mu)()(i,j));
|
||||
for (int i = 0; i < Nrepresentation; i++) {
|
||||
for (int j = 0; j < Nrepresentation; 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)
|
||||
{
|
||||
SiteScalarGaugeField ScalarUmu;
|
||||
|
||||
inline void DoubleStore(GridBase *GaugeGrid, DoubledGaugeField &Uds,
|
||||
const GaugeField &Umu) {
|
||||
SiteScalarGaugeField ScalarUmu;
|
||||
SiteDoubledGaugeField ScalarUds;
|
||||
|
||||
GaugeLinkField U (Umu._grid);
|
||||
GaugeField Uadj(Umu._grid);
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U = PeekIndex<LorentzIndex>(Umu,mu);
|
||||
U = adj(Cshift(U,mu,-1));
|
||||
PokeIndex<LorentzIndex>(Uadj,U,mu);
|
||||
}
|
||||
|
||||
for(int lidx=0;lidx<GaugeGrid->lSites();lidx++){
|
||||
std::vector<int> lcoor;
|
||||
GaugeGrid->LocalIndexToLocalCoor(lidx,lcoor);
|
||||
|
||||
peekLocalSite(ScalarUmu,Umu,lcoor);
|
||||
for(int mu=0;mu<4;mu++) ScalarUds(mu) = ScalarUmu(mu);
|
||||
|
||||
peekLocalSite(ScalarUmu,Uadj,lcoor);
|
||||
for(int mu=0;mu<4;mu++) ScalarUds(mu+4) = ScalarUmu(mu);
|
||||
|
||||
pokeLocalSite(ScalarUds,Uds,lcoor);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){
|
||||
assert(0);
|
||||
}
|
||||
|
||||
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){
|
||||
GaugeLinkField U(Umu._grid);
|
||||
GaugeField Uadj(Umu._grid);
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
U = PeekIndex<LorentzIndex>(Umu, mu);
|
||||
U = adj(Cshift(U, mu, -1));
|
||||
PokeIndex<LorentzIndex>(Uadj, U, mu);
|
||||
}
|
||||
|
||||
for (int lidx = 0; lidx < GaugeGrid->lSites(); lidx++) {
|
||||
std::vector<int> lcoor;
|
||||
GaugeGrid->LocalIndexToLocalCoor(lidx, lcoor);
|
||||
|
||||
peekLocalSite(ScalarUmu, Umu, lcoor);
|
||||
for (int mu = 0; mu < 4; mu++) ScalarUds(mu) = ScalarUmu(mu);
|
||||
|
||||
peekLocalSite(ScalarUmu, Uadj, lcoor);
|
||||
for (int mu = 0; mu < 4; mu++) ScalarUds(mu + 4) = ScalarUmu(mu);
|
||||
|
||||
pokeLocalSite(ScalarUds, Uds, lcoor);
|
||||
}
|
||||
}
|
||||
|
||||
inline void InsertForce4D(GaugeField &mat, FermionField &Btilde,
|
||||
FermionField &A, int mu) {
|
||||
assert(0);
|
||||
}
|
||||
|
||||
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde,
|
||||
FermionField Ã, int mu) {
|
||||
assert(0);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////////////
|
||||
// Flavour doubled spinors; is Gparity the only? what about C*?
|
||||
////////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
template<class S,int Nrepresentation,class _Coeff_t = RealD>
|
||||
class GparityWilsonImpl : public ConjugateGaugeImpl< GaugeImplTypes<S,Nrepresentation> >{
|
||||
|
||||
template <class S, int Nrepresentation,class _Coeff_t = RealD>
|
||||
class GparityWilsonImpl
|
||||
: public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresentation> > {
|
||||
public:
|
||||
static const int Dimension = Nrepresentation;
|
||||
|
||||
const bool LsVectorised=false;
|
||||
|
||||
typedef _Coeff_t Coeff_t;
|
||||
typedef ConjugateGaugeImpl< GaugeImplTypes<S,Nrepresentation> > Gimpl;
|
||||
|
||||
|
||||
INHERIT_GIMPL_TYPES(Gimpl);
|
||||
|
||||
template<typename vtype> using iImplSpinor = 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 iImplHalfSpinor<Simd> SiteHalfSpinor;
|
||||
typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
|
||||
|
||||
typedef Lattice<SiteSpinor> FermionField;
|
||||
|
||||
template <typename vtype>
|
||||
using iImplSpinor =
|
||||
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 iImplHalfSpinor<Simd> SiteHalfSpinor;
|
||||
typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
|
||||
|
||||
typedef Lattice<SiteSpinor> FermionField;
|
||||
typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
|
||||
|
||||
typedef WilsonCompressor<SiteHalfSpinor,SiteSpinor> Compressor;
|
||||
typedef WilsonStencil<SiteSpinor,SiteHalfSpinor> StencilImpl;
|
||||
|
||||
typedef WilsonCompressor<SiteHalfSpinor, SiteSpinor> Compressor;
|
||||
typedef WilsonStencil<SiteSpinor, SiteHalfSpinor> StencilImpl;
|
||||
|
||||
typedef GparityWilsonImplParams ImplParams;
|
||||
|
||||
|
||||
ImplParams Params;
|
||||
|
||||
GparityWilsonImpl(const ImplParams &p= ImplParams()) : Params(p) {};
|
||||
|
||||
|
||||
GparityWilsonImpl(const ImplParams &p = ImplParams()) : Params(p){};
|
||||
|
||||
bool overlapCommsCompute(void) { return Params.overlapCommsCompute; };
|
||||
|
||||
// provide the multiply by link that is differentiated between Gparity (with flavour index) and non-Gparity
|
||||
inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE,StencilImpl &St){
|
||||
|
||||
// provide the multiply by link that is differentiated between Gparity (with
|
||||
// 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 typename SiteHalfSpinor::scalar_object sobj;
|
||||
|
||||
|
||||
vobj vtmp;
|
||||
sobj stmp;
|
||||
|
||||
GridBase *grid = St._grid;
|
||||
|
||||
|
||||
const int Nsimd = grid->Nsimd();
|
||||
|
||||
int direction = St._directions[mu];
|
||||
int distance = St._distances[mu];
|
||||
int ptype = St._permute_type[mu];
|
||||
int sl = St._grid->_simd_layout[direction];
|
||||
|
||||
int direction = St._directions[mu];
|
||||
int distance = St._distances[mu];
|
||||
int ptype = St._permute_type[mu];
|
||||
int sl = St._grid->_simd_layout[direction];
|
||||
|
||||
// Fixme X.Y.Z.T hardcode in stencil
|
||||
int mmu = mu % Nd;
|
||||
|
||||
int mmu = mu % Nd;
|
||||
|
||||
// assert our assumptions
|
||||
assert((distance==1)||(distance==-1)); // nearest neighbour stencil hard code
|
||||
assert((sl==1)||(sl==2));
|
||||
assert((distance == 1) || (distance == -1)); // nearest neighbour stencil hard code
|
||||
assert((sl == 1) || (sl == 2));
|
||||
|
||||
std::vector<int> icoor;
|
||||
|
||||
|
||||
if ( SE->_around_the_world && Params.twists[mmu] ) {
|
||||
|
||||
if ( sl == 2 ) {
|
||||
@ -391,7 +415,7 @@ PARALLEL_FOR_LOOP
|
||||
mult(&phi(1),&U(1)(mu),&chi(1));
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu)
|
||||
{
|
||||
@ -404,7 +428,7 @@ PARALLEL_FOR_LOOP
|
||||
GaugeLinkField Uconj(GaugeGrid);
|
||||
|
||||
Lattice<iScalar<vInteger> > coor(GaugeGrid);
|
||||
|
||||
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
|
||||
@ -412,19 +436,19 @@ PARALLEL_FOR_LOOP
|
||||
|
||||
U = PeekIndex<LorentzIndex>(Umu,mu);
|
||||
Uconj = conjugate(U);
|
||||
|
||||
|
||||
// This phase could come from a simple bc 1,1,-1,1 ..
|
||||
int neglink = GaugeGrid->GlobalDimensions()[mu]-1;
|
||||
if ( Params.twists[mu] ) {
|
||||
Uconj = where(coor==neglink,-Uconj,Uconj);
|
||||
}
|
||||
|
||||
|
||||
PARALLEL_FOR_LOOP
|
||||
for(auto ss=U.begin();ss<U.end();ss++){
|
||||
Uds[ss](0)(mu) = U[ss]();
|
||||
Uds[ss](1)(mu) = Uconj[ss]();
|
||||
}
|
||||
|
||||
PARALLEL_FOR_LOOP
|
||||
for(auto ss=U.begin();ss<U.end();ss++){
|
||||
Uds[ss](0)(mu) = U[ss]();
|
||||
Uds[ss](1)(mu) = Uconj[ss]();
|
||||
}
|
||||
|
||||
U = adj(Cshift(U ,mu,-1)); // correct except for spanning the boundary
|
||||
Uconj = adj(Cshift(Uconj,mu,-1));
|
||||
@ -434,80 +458,86 @@ PARALLEL_FOR_LOOP
|
||||
Utmp = where(coor==0,Uconj,Utmp);
|
||||
}
|
||||
|
||||
PARALLEL_FOR_LOOP
|
||||
for(auto ss=U.begin();ss<U.end();ss++){
|
||||
Uds[ss](0)(mu+4) = Utmp[ss]();
|
||||
}
|
||||
PARALLEL_FOR_LOOP
|
||||
for(auto ss=U.begin();ss<U.end();ss++){
|
||||
Uds[ss](0)(mu+4) = Utmp[ss]();
|
||||
}
|
||||
|
||||
Utmp = Uconj;
|
||||
if ( Params.twists[mu] ) {
|
||||
Utmp = where(coor==0,U,Utmp);
|
||||
}
|
||||
|
||||
PARALLEL_FOR_LOOP
|
||||
for(auto ss=U.begin();ss<U.end();ss++){
|
||||
Uds[ss](1)(mu+4) = Utmp[ss]();
|
||||
}
|
||||
PARALLEL_FOR_LOOP
|
||||
for(auto ss=U.begin();ss<U.end();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.
|
||||
GaugeLinkField link(mat._grid);
|
||||
// use lorentz for flavour as hack.
|
||||
PARALLEL_FOR_LOOP
|
||||
for(auto ss=link.begin();ss<link.end();ss++){
|
||||
auto ttmp = traceIndex<SpinIndex>(outerProduct(Btilde[ss],A[ss]));
|
||||
link[ss]() = ttmp(0,0) + conjugate(ttmp(1,1)) ;
|
||||
}
|
||||
PokeIndex<LorentzIndex>(mat,link,mu);
|
||||
auto tmp = TraceIndex<SpinIndex>(outerProduct(Btilde, A));
|
||||
PARALLEL_FOR_LOOP
|
||||
for (auto ss = tmp.begin(); ss < tmp.end(); ss++) {
|
||||
link[ss]() = tmp[ss](0, 0) - conjugate(tmp[ss](1, 1));
|
||||
}
|
||||
PokeIndex<LorentzIndex>(mat, link, mu);
|
||||
return;
|
||||
}
|
||||
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){
|
||||
|
||||
int Ls=Btilde._grid->_fdimensions[0];
|
||||
|
||||
|
||||
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde,
|
||||
FermionField Ã, int mu) {
|
||||
int Ls = Btilde._grid->_fdimensions[0];
|
||||
|
||||
GaugeLinkField tmp(mat._grid);
|
||||
tmp = zero;
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int ss=0;ss<tmp._grid->oSites();ss++){
|
||||
for(int s=0;s<Ls;s++){
|
||||
int sF = s+Ls*ss;
|
||||
auto ttmp = traceIndex<SpinIndex>(outerProduct(Btilde[sF],Atilde[sF]));
|
||||
tmp[ss]() = tmp[ss]()+ ttmp(0,0) + conjugate(ttmp(1,1));
|
||||
PARALLEL_FOR_LOOP
|
||||
for (int ss = 0; ss < tmp._grid->oSites(); ss++) {
|
||||
for (int s = 0; s < Ls; s++) {
|
||||
int sF = s + Ls * ss;
|
||||
auto ttmp = traceIndex<SpinIndex>(outerProduct(Btilde[sF], Atilde[sF]));
|
||||
tmp[ss]() = tmp[ss]() + ttmp(0, 0) + conjugate(ttmp(1, 1));
|
||||
}
|
||||
}
|
||||
}
|
||||
PokeIndex<LorentzIndex>(mat,tmp,mu);
|
||||
PokeIndex<LorentzIndex>(mat, tmp, mu);
|
||||
return;
|
||||
}
|
||||
};
|
||||
|
||||
typedef WilsonImpl<vComplex ,Nc> WilsonImplR; // Real.. whichever prec
|
||||
typedef WilsonImpl<vComplexF,Nc> WilsonImplF; // Float
|
||||
typedef WilsonImpl<vComplexD,Nc> WilsonImplD; // Double
|
||||
typedef WilsonImpl<vComplex, FundamentalRepresentation > WilsonImplR; // Real.. whichever prec
|
||||
typedef WilsonImpl<vComplexF, FundamentalRepresentation > WilsonImplF; // Float
|
||||
typedef WilsonImpl<vComplexD, FundamentalRepresentation > WilsonImplD; // Double
|
||||
|
||||
typedef WilsonImpl<vComplex ,Nc,ComplexD> ZWilsonImplR; // Real.. whichever prec
|
||||
typedef WilsonImpl<vComplexF,Nc,ComplexD> ZWilsonImplF; // Float
|
||||
typedef WilsonImpl<vComplexD,Nc,ComplexD> ZWilsonImplD; // Double
|
||||
|
||||
typedef WilsonImpl<vComplex, FundamentalRepresentation, ComplexD > ZWilsonImplR; // Real.. whichever prec
|
||||
typedef WilsonImpl<vComplexF, FundamentalRepresentation, ComplexD > ZWilsonImplF; // Float
|
||||
typedef WilsonImpl<vComplexD, FundamentalRepresentation, ComplexD > ZWilsonImplD; // Double
|
||||
|
||||
typedef WilsonImpl<vComplex, AdjointRepresentation > WilsonAdjImplR; // Real.. whichever prec
|
||||
typedef WilsonImpl<vComplexF, AdjointRepresentation > WilsonAdjImplF; // Float
|
||||
typedef WilsonImpl<vComplexD, AdjointRepresentation > WilsonAdjImplD; // Double
|
||||
|
||||
typedef WilsonImpl<vComplex, TwoIndexSymmetricRepresentation > WilsonTwoIndexSymmetricImplR; // Real.. whichever prec
|
||||
typedef WilsonImpl<vComplexF, TwoIndexSymmetricRepresentation > WilsonTwoIndexSymmetricImplF; // Float
|
||||
typedef WilsonImpl<vComplexD, TwoIndexSymmetricRepresentation > WilsonTwoIndexSymmetricImplD; // Double
|
||||
|
||||
typedef DomainWallVec5dImpl<vComplex ,Nc> DomainWallVec5dImplR; // Real.. whichever prec
|
||||
typedef DomainWallVec5dImpl<vComplexF,Nc> DomainWallVec5dImplF; // Float
|
||||
typedef DomainWallVec5dImpl<vComplexD,Nc> DomainWallVec5dImplD; // Double
|
||||
|
||||
|
||||
typedef DomainWallVec5dImpl<vComplex ,Nc,ComplexD> ZDomainWallVec5dImplR; // Real.. whichever prec
|
||||
typedef DomainWallVec5dImpl<vComplexF,Nc,ComplexD> ZDomainWallVec5dImplF; // Float
|
||||
typedef DomainWallVec5dImpl<vComplexD,Nc,ComplexD> ZDomainWallVec5dImplD; // Double
|
||||
|
||||
typedef DomainWallVec5dImpl<vComplex ,Nc> DomainWallVec5dImplR; // Real.. whichever prec
|
||||
typedef DomainWallVec5dImpl<vComplexF,Nc> DomainWallVec5dImplF; // Float
|
||||
typedef DomainWallVec5dImpl<vComplexD,Nc> DomainWallVec5dImplD; // Double
|
||||
|
||||
typedef GparityWilsonImpl<vComplex ,Nc> GparityWilsonImplR; // Real.. whichever prec
|
||||
typedef GparityWilsonImpl<vComplexF,Nc> GparityWilsonImplF; // Float
|
||||
typedef GparityWilsonImpl<vComplexD,Nc> GparityWilsonImplD; // Double
|
||||
|
||||
}
|
||||
typedef GparityWilsonImpl<vComplex, Nc> GparityWilsonImplR; // Real.. whichever prec
|
||||
typedef GparityWilsonImpl<vComplexF, Nc> GparityWilsonImplF; // Float
|
||||
typedef GparityWilsonImpl<vComplexD, Nc> GparityWilsonImplD; // Double
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
@ -1,130 +1,129 @@
|
||||
/*************************************************************************************
|
||||
/*************************************************************************************
|
||||
|
||||
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 <paboyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
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.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
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.,
|
||||
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 */
|
||||
See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid.h>
|
||||
|
||||
namespace Grid {
|
||||
namespace QCD {
|
||||
|
||||
const std::vector<int> WilsonFermionStatic::directions ({0,1,2,3, 0, 1, 2, 3});
|
||||
const std::vector<int> WilsonFermionStatic::displacements({1,1,1,1,-1,-1,-1,-1});
|
||||
int WilsonFermionStatic::HandOptDslash;
|
||||
const std::vector<int> WilsonFermionStatic::directions({0, 1, 2, 3, 0, 1, 2,
|
||||
3});
|
||||
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>
|
||||
WilsonFermion<Impl>::WilsonFermion(GaugeField &_Umu,
|
||||
GridCartesian &Fgrid,
|
||||
GridRedBlackCartesian &Hgrid,
|
||||
RealD _mass,const ImplParams &p) :
|
||||
Kernels(p),
|
||||
_grid(&Fgrid),
|
||||
_cbgrid(&Hgrid),
|
||||
Stencil (&Fgrid,npoint,Even,directions,displacements),
|
||||
StencilEven(&Hgrid,npoint,Even,directions,displacements), // source is Even
|
||||
StencilOdd (&Hgrid,npoint,Odd ,directions,displacements), // source is Odd
|
||||
mass(_mass),
|
||||
Lebesgue(_grid),
|
||||
LebesgueEvenOdd(_cbgrid),
|
||||
Umu(&Fgrid),
|
||||
UmuEven(&Hgrid),
|
||||
UmuOdd (&Hgrid)
|
||||
{
|
||||
// Allocate the required comms buffer
|
||||
ImportGauge(_Umu);
|
||||
template <class Impl>
|
||||
WilsonFermion<Impl>::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
|
||||
GridRedBlackCartesian &Hgrid, RealD _mass,
|
||||
const ImplParams &p)
|
||||
: Kernels(p),
|
||||
_grid(&Fgrid),
|
||||
_cbgrid(&Hgrid),
|
||||
Stencil(&Fgrid, npoint, Even, directions, displacements),
|
||||
StencilEven(&Hgrid, npoint, Even, directions,
|
||||
displacements), // source is Even
|
||||
StencilOdd(&Hgrid, npoint, Odd, directions,
|
||||
displacements), // source is Odd
|
||||
mass(_mass),
|
||||
Lebesgue(_grid),
|
||||
LebesgueEvenOdd(_cbgrid),
|
||||
Umu(&Fgrid),
|
||||
UmuEven(&Hgrid),
|
||||
UmuOdd(&Hgrid) {
|
||||
// Allocate the required comms buffer
|
||||
ImportGauge(_Umu);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::ImportGauge(const GaugeField &_Umu) {
|
||||
GaugeField HUmu(_Umu._grid);
|
||||
HUmu = _Umu * (-0.5);
|
||||
Impl::DoubleStore(GaugeGrid(), Umu, HUmu);
|
||||
pickCheckerboard(Even, UmuEven, Umu);
|
||||
pickCheckerboard(Odd, UmuOdd, Umu);
|
||||
}
|
||||
|
||||
/////////////////////////////
|
||||
// Implement the interface
|
||||
/////////////////////////////
|
||||
|
||||
template <class Impl>
|
||||
RealD WilsonFermion<Impl>::M(const FermionField &in, FermionField &out) {
|
||||
out.checkerboard = in.checkerboard;
|
||||
Dhop(in, out, DaggerNo);
|
||||
return axpy_norm(out, 4 + mass, in, out);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
RealD WilsonFermion<Impl>::Mdag(const FermionField &in, FermionField &out) {
|
||||
out.checkerboard = in.checkerboard;
|
||||
Dhop(in, out, DaggerYes);
|
||||
return axpy_norm(out, 4 + mass, in, out);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::Meooe(const FermionField &in, FermionField &out) {
|
||||
if (in.checkerboard == Odd) {
|
||||
DhopEO(in, out, DaggerNo);
|
||||
} else {
|
||||
DhopOE(in, out, DaggerNo);
|
||||
}
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::ImportGauge(const GaugeField &_Umu)
|
||||
{
|
||||
GaugeField HUmu(_Umu._grid);
|
||||
HUmu = _Umu*(-0.5);
|
||||
Impl::DoubleStore(GaugeGrid(),Umu,HUmu);
|
||||
pickCheckerboard(Even,UmuEven,Umu);
|
||||
pickCheckerboard(Odd ,UmuOdd,Umu);
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
|
||||
if (in.checkerboard == Odd) {
|
||||
DhopEO(in, out, DaggerYes);
|
||||
} else {
|
||||
DhopOE(in, out, DaggerYes);
|
||||
}
|
||||
}
|
||||
|
||||
/////////////////////////////
|
||||
// Implement the interface
|
||||
/////////////////////////////
|
||||
|
||||
template<class Impl>
|
||||
RealD WilsonFermion<Impl>::M(const FermionField &in, FermionField &out)
|
||||
{
|
||||
out.checkerboard=in.checkerboard;
|
||||
Dhop(in,out,DaggerNo);
|
||||
return axpy_norm(out,4+mass,in,out);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
RealD WilsonFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
|
||||
{
|
||||
out.checkerboard=in.checkerboard;
|
||||
Dhop(in,out,DaggerYes);
|
||||
return axpy_norm(out,4+mass,in,out);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::Meooe(const FermionField &in, FermionField &out)
|
||||
{
|
||||
if ( in.checkerboard == Odd ) {
|
||||
DhopEO(in,out,DaggerNo);
|
||||
} else {
|
||||
DhopOE(in,out,DaggerNo);
|
||||
}
|
||||
}
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out)
|
||||
{
|
||||
if ( in.checkerboard == Odd ) {
|
||||
DhopEO(in,out,DaggerYes);
|
||||
} else {
|
||||
DhopOE(in,out,DaggerYes);
|
||||
}
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::Mooee(const FermionField &in, FermionField &out) {
|
||||
out.checkerboard = in.checkerboard;
|
||||
typename FermionField::scalar_type scal(4.0+mass);
|
||||
out = scal*in;
|
||||
typename FermionField::scalar_type scal(4.0 + mass);
|
||||
out = scal * in;
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) {
|
||||
out.checkerboard = in.checkerboard;
|
||||
Mooee(in,out);
|
||||
Mooee(in, out);
|
||||
}
|
||||
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) {
|
||||
out.checkerboard = in.checkerboard;
|
||||
@ -136,8 +135,9 @@ namespace QCD {
|
||||
out.checkerboard = in.checkerboard;
|
||||
MooeeInv(in,out);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>:: MomentumSpacePropagator(FermionField &out, const FermionField &in,RealD _m) {
|
||||
void WilsonFermion<Impl>::MomentumSpacePropagator(FermionField &out, const FermionField &in,RealD _m) {
|
||||
|
||||
// what type LatticeComplex
|
||||
conformable(_grid,out._grid);
|
||||
@ -190,184 +190,182 @@ namespace QCD {
|
||||
out = num*denom; // [ -i gmu sin k + 2 sin^2 k/2 + m] / [ sin^2 k + (2 sin^2 k/2 + m)^2 ]
|
||||
|
||||
}
|
||||
|
||||
///////////////////////////////////
|
||||
// Internal
|
||||
///////////////////////////////////
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::DerivInternal(StencilImpl & st,
|
||||
DoubledGaugeField & U,
|
||||
GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,int dag) {
|
||||
|
||||
assert((dag==DaggerNo) ||(dag==DaggerYes));
|
||||
|
||||
Compressor compressor(dag);
|
||||
|
||||
FermionField Btilde(B._grid);
|
||||
FermionField Atilde(B._grid);
|
||||
Atilde = A;
|
||||
|
||||
st.HaloExchange(B,compressor);
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// Flip gamma (1+g)<->(1-g) if dag
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
int gamma = mu;
|
||||
if ( !dag ) gamma+= Nd;
|
||||
|
||||
////////////////////////
|
||||
// Call the single hop
|
||||
////////////////////////
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int sss=0;sss<B._grid->oSites();sss++){
|
||||
Kernels::DiracOptDhopDir(st,U,st.comm_buf,sss,sss,B,Btilde,mu,gamma);
|
||||
}
|
||||
|
||||
//////////////////////////////////////////////////
|
||||
// spin trace outer product
|
||||
//////////////////////////////////////////////////
|
||||
Impl::InsertForce4D(mat,Btilde,Atilde,mu);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)
|
||||
{
|
||||
conformable(U._grid,_grid);
|
||||
conformable(U._grid,V._grid);
|
||||
conformable(U._grid,mat._grid);
|
||||
|
||||
mat.checkerboard = U.checkerboard;
|
||||
|
||||
DerivInternal(Stencil,Umu,mat,U,V,dag);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)
|
||||
{
|
||||
conformable(U._grid,_cbgrid);
|
||||
conformable(U._grid,V._grid);
|
||||
conformable(U._grid,mat._grid);
|
||||
|
||||
assert(V.checkerboard==Even);
|
||||
assert(U.checkerboard==Odd);
|
||||
mat.checkerboard = Odd;
|
||||
|
||||
DerivInternal(StencilEven,UmuOdd,mat,U,V,dag);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)
|
||||
{
|
||||
conformable(U._grid,_cbgrid);
|
||||
conformable(U._grid,V._grid);
|
||||
conformable(U._grid,mat._grid);
|
||||
|
||||
assert(V.checkerboard==Odd);
|
||||
assert(U.checkerboard==Even);
|
||||
mat.checkerboard = Even;
|
||||
|
||||
DerivInternal(StencilOdd,UmuEven,mat,U,V,dag);
|
||||
}
|
||||
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out,int dag) {
|
||||
conformable(in._grid,_grid); // verifies full grid
|
||||
conformable(in._grid,out._grid);
|
||||
|
||||
out.checkerboard = in.checkerboard;
|
||||
|
||||
DhopInternal(Stencil,Lebesgue,Umu,in,out,dag);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out,int dag) {
|
||||
conformable(in._grid,_cbgrid); // verifies half grid
|
||||
conformable(in._grid,out._grid); // drops the cb check
|
||||
|
||||
assert(in.checkerboard==Even);
|
||||
out.checkerboard = Odd;
|
||||
|
||||
DhopInternal(StencilEven,LebesgueEvenOdd,UmuOdd,in,out,dag);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::DhopEO(const FermionField &in, FermionField &out,int dag) {
|
||||
conformable(in._grid,_cbgrid); // verifies half grid
|
||||
conformable(in._grid,out._grid); // drops the cb check
|
||||
|
||||
assert(in.checkerboard==Odd);
|
||||
out.checkerboard = Even;
|
||||
|
||||
DhopInternal(StencilOdd,LebesgueEvenOdd,UmuEven,in,out,dag);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::Mdir (const FermionField &in, FermionField &out,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;
|
||||
|
||||
DhopDirDisp(in,out,dirdisp,gamma,DaggerNo);
|
||||
|
||||
};
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::DhopDirDisp(const FermionField &in, FermionField &out,int dirdisp,int gamma,int dag) {
|
||||
|
||||
Compressor compressor(dag);
|
||||
|
||||
Stencil.HaloExchange(in,compressor);
|
||||
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int sss=0;sss<in._grid->oSites();sss++){
|
||||
Kernels::DiracOptDhopDir(Stencil,Umu,Stencil.comm_buf,sss,sss,in,out,dirdisp,gamma);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion<Impl>::DhopInternal(StencilImpl & st,LebesgueOrder& lo,DoubledGaugeField & U,
|
||||
const FermionField &in, FermionField &out,int dag)
|
||||
{
|
||||
assert((dag==DaggerNo) ||(dag==DaggerYes));
|
||||
|
||||
Compressor compressor(dag);
|
||||
st.HaloExchange(in,compressor);
|
||||
|
||||
if ( dag == DaggerYes ) {
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int sss=0;sss<in._grid->oSites();sss++){
|
||||
Kernels::DiracOptDhopSiteDag(st,lo,U,st.comm_buf,sss,sss,1,1,in,out);
|
||||
}
|
||||
} else {
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int sss=0;sss<in._grid->oSites();sss++){
|
||||
Kernels::DiracOptDhopSite(st,lo,U,st.comm_buf,sss,sss,1,1,in,out);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
FermOpTemplateInstantiate(WilsonFermion);
|
||||
GparityFermOpTemplateInstantiate(WilsonFermion);
|
||||
|
||||
///////////////////////////////////
|
||||
// Internal
|
||||
///////////////////////////////////
|
||||
|
||||
}}
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U,
|
||||
GaugeField &mat, const FermionField &A,
|
||||
const FermionField &B, int dag) {
|
||||
assert((dag == DaggerNo) || (dag == DaggerYes));
|
||||
|
||||
Compressor compressor(dag);
|
||||
|
||||
FermionField Btilde(B._grid);
|
||||
FermionField Atilde(B._grid);
|
||||
Atilde = A;
|
||||
|
||||
st.HaloExchange(B, compressor);
|
||||
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// Flip gamma (1+g)<->(1-g) if dag
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
int gamma = mu;
|
||||
if (!dag) gamma += Nd;
|
||||
|
||||
////////////////////////
|
||||
// Call the single hop
|
||||
////////////////////////
|
||||
PARALLEL_FOR_LOOP
|
||||
for (int sss = 0; sss < B._grid->oSites(); sss++) {
|
||||
Kernels::DiracOptDhopDir(st, U, st.comm_buf, sss, sss, B, Btilde, mu,
|
||||
gamma);
|
||||
}
|
||||
|
||||
//////////////////////////////////////////////////
|
||||
// spin trace outer product
|
||||
//////////////////////////////////////////////////
|
||||
Impl::InsertForce4D(mat, Btilde, Atilde, mu);
|
||||
}
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U,
|
||||
const FermionField &V, int dag) {
|
||||
conformable(U._grid, _grid);
|
||||
conformable(U._grid, V._grid);
|
||||
conformable(U._grid, mat._grid);
|
||||
|
||||
mat.checkerboard = U.checkerboard;
|
||||
|
||||
DerivInternal(Stencil, Umu, mat, U, V, dag);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U,
|
||||
const FermionField &V, int dag) {
|
||||
conformable(U._grid, _cbgrid);
|
||||
conformable(U._grid, V._grid);
|
||||
conformable(U._grid, mat._grid);
|
||||
|
||||
assert(V.checkerboard == Even);
|
||||
assert(U.checkerboard == Odd);
|
||||
mat.checkerboard = Odd;
|
||||
|
||||
DerivInternal(StencilEven, UmuOdd, mat, U, V, dag);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U,
|
||||
const FermionField &V, int dag) {
|
||||
conformable(U._grid, _cbgrid);
|
||||
conformable(U._grid, V._grid);
|
||||
conformable(U._grid, mat._grid);
|
||||
|
||||
assert(V.checkerboard == Odd);
|
||||
assert(U.checkerboard == Even);
|
||||
mat.checkerboard = Even;
|
||||
|
||||
DerivInternal(StencilOdd, UmuEven, mat, U, V, dag);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::Dhop(const FermionField &in, FermionField &out,
|
||||
int dag) {
|
||||
conformable(in._grid, _grid); // verifies full grid
|
||||
conformable(in._grid, out._grid);
|
||||
|
||||
out.checkerboard = in.checkerboard;
|
||||
|
||||
DhopInternal(Stencil, Lebesgue, Umu, in, out, dag);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::DhopOE(const FermionField &in, FermionField &out,
|
||||
int dag) {
|
||||
conformable(in._grid, _cbgrid); // verifies half grid
|
||||
conformable(in._grid, out._grid); // drops the cb check
|
||||
|
||||
assert(in.checkerboard == Even);
|
||||
out.checkerboard = Odd;
|
||||
|
||||
DhopInternal(StencilEven, LebesgueEvenOdd, UmuOdd, in, out, dag);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::DhopEO(const FermionField &in, FermionField &out,
|
||||
int dag) {
|
||||
conformable(in._grid, _cbgrid); // verifies half grid
|
||||
conformable(in._grid, out._grid); // drops the cb check
|
||||
|
||||
assert(in.checkerboard == Odd);
|
||||
out.checkerboard = Even;
|
||||
|
||||
DhopInternal(StencilOdd, LebesgueEvenOdd, UmuEven, in, out, dag);
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::Mdir(const FermionField &in, FermionField &out,
|
||||
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;
|
||||
|
||||
DhopDirDisp(in, out, dirdisp, gamma, DaggerNo);
|
||||
};
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::DhopDirDisp(const FermionField &in, FermionField &out,
|
||||
int dirdisp, int gamma, int dag) {
|
||||
Compressor compressor(dag);
|
||||
|
||||
Stencil.HaloExchange(in, compressor);
|
||||
|
||||
PARALLEL_FOR_LOOP
|
||||
for (int sss = 0; sss < in._grid->oSites(); sss++) {
|
||||
Kernels::DiracOptDhopDir(Stencil, Umu, Stencil.comm_buf, sss, sss, in, out,
|
||||
dirdisp, gamma);
|
||||
}
|
||||
};
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder &lo,
|
||||
DoubledGaugeField &U,
|
||||
const FermionField &in,
|
||||
FermionField &out, int dag) {
|
||||
assert((dag == DaggerNo) || (dag == DaggerYes));
|
||||
|
||||
Compressor compressor(dag);
|
||||
st.HaloExchange(in, compressor);
|
||||
|
||||
if (dag == DaggerYes) {
|
||||
PARALLEL_FOR_LOOP
|
||||
for (int sss = 0; sss < in._grid->oSites(); sss++) {
|
||||
Kernels::DiracOptDhopSiteDag(st, lo, U, st.comm_buf, sss, sss, 1, 1, in,
|
||||
out);
|
||||
}
|
||||
} else {
|
||||
PARALLEL_FOR_LOOP
|
||||
for (int sss = 0; sss < in._grid->oSites(); sss++) {
|
||||
Kernels::DiracOptDhopSite(st, lo, U, st.comm_buf, sss, sss, 1, 1, in,
|
||||
out);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
FermOpTemplateInstantiate(WilsonFermion);
|
||||
AdjointFermOpTemplateInstantiate(WilsonFermion);
|
||||
TwoIndexFermOpTemplateInstantiate(WilsonFermion);
|
||||
GparityFermOpTemplateInstantiate(WilsonFermion);
|
||||
}
|
||||
}
|
||||
|
||||
@ -1,163 +1,154 @@
|
||||
/*************************************************************************************
|
||||
/*************************************************************************************
|
||||
|
||||
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 <paboyle@ph.ed.ac.uk>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
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.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
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.,
|
||||
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 */
|
||||
#ifndef GRID_QCD_WILSON_FERMION_H
|
||||
#define GRID_QCD_WILSON_FERMION_H
|
||||
See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#ifndef GRID_QCD_WILSON_FERMION_H
|
||||
#define GRID_QCD_WILSON_FERMION_H
|
||||
|
||||
namespace Grid {
|
||||
|
||||
namespace QCD {
|
||||
namespace QCD {
|
||||
|
||||
class WilsonFermionStatic {
|
||||
public:
|
||||
static int HandOptDslash; // these are a temporary hack
|
||||
static int MortonOrder;
|
||||
static const std::vector<int> directions ;
|
||||
static const std::vector<int> displacements;
|
||||
static const int npoint=8;
|
||||
};
|
||||
class WilsonFermionStatic {
|
||||
public:
|
||||
static int HandOptDslash; // these are a temporary hack
|
||||
static int MortonOrder;
|
||||
static const std::vector<int> directions;
|
||||
static const std::vector<int> displacements;
|
||||
static const int npoint = 8;
|
||||
};
|
||||
|
||||
template<class Impl>
|
||||
class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic
|
||||
{
|
||||
public:
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
typedef WilsonKernels<Impl> Kernels;
|
||||
template <class Impl>
|
||||
class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
|
||||
public:
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
typedef WilsonKernels<Impl> Kernels;
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Implement the abstract base
|
||||
///////////////////////////////////////////////////////////////
|
||||
GridBase *GaugeGrid(void) { return _grid ;}
|
||||
GridBase *GaugeRedBlackGrid(void) { return _cbgrid ;}
|
||||
GridBase *FermionGrid(void) { return _grid;}
|
||||
GridBase *FermionRedBlackGrid(void) { return _cbgrid;}
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Implement the abstract base
|
||||
///////////////////////////////////////////////////////////////
|
||||
GridBase *GaugeGrid(void) { return _grid; }
|
||||
GridBase *GaugeRedBlackGrid(void) { return _cbgrid; }
|
||||
GridBase *FermionGrid(void) { return _grid; }
|
||||
GridBase *FermionRedBlackGrid(void) { return _cbgrid; }
|
||||
|
||||
//////////////////////////////////////////////////////////////////
|
||||
// override multiply; cut number routines if pass dagger argument
|
||||
// and also make interface more uniformly consistent
|
||||
//////////////////////////////////////////////////////////////////
|
||||
RealD M(const FermionField &in, FermionField &out);
|
||||
RealD Mdag(const FermionField &in, FermionField &out);
|
||||
//////////////////////////////////////////////////////////////////
|
||||
// override multiply; cut number routines if pass dagger argument
|
||||
// and also make interface more uniformly consistent
|
||||
//////////////////////////////////////////////////////////////////
|
||||
RealD M(const FermionField &in, FermionField &out);
|
||||
RealD Mdag(const FermionField &in, FermionField &out);
|
||||
|
||||
/////////////////////////////////////////////////////////
|
||||
// half checkerboard operations
|
||||
// could remain virtual so we can derive Clover from Wilson base
|
||||
/////////////////////////////////////////////////////////
|
||||
void Meooe(const FermionField &in, FermionField &out) ;
|
||||
void MeooeDag(const FermionField &in, FermionField &out) ;
|
||||
/////////////////////////////////////////////////////////
|
||||
// half checkerboard operations
|
||||
// could remain virtual so we can derive Clover from Wilson base
|
||||
/////////////////////////////////////////////////////////
|
||||
void Meooe(const FermionField &in, FermionField &out);
|
||||
void MeooeDag(const FermionField &in, FermionField &out);
|
||||
|
||||
// allow override for twisted mass and clover
|
||||
virtual void Mooee(const FermionField &in, FermionField &out) ;
|
||||
virtual void MooeeDag(const FermionField &in, FermionField &out) ;
|
||||
virtual void MooeeInv(const FermionField &in, FermionField &out) ;
|
||||
virtual void MooeeInvDag(const FermionField &in, FermionField &out) ;
|
||||
// allow override for twisted mass and clover
|
||||
virtual void Mooee(const FermionField &in, FermionField &out);
|
||||
virtual void MooeeDag(const FermionField &in, FermionField &out);
|
||||
virtual void MooeeInv(const FermionField &in, FermionField &out);
|
||||
virtual void MooeeInvDag(const FermionField &in, FermionField &out);
|
||||
|
||||
virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _mass) ;
|
||||
virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _mass) ;
|
||||
|
||||
////////////////////////
|
||||
// Derivative interface
|
||||
////////////////////////
|
||||
// Interface calls an internal routine
|
||||
void DhopDeriv(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);
|
||||
////////////////////////
|
||||
// Derivative interface
|
||||
////////////////////////
|
||||
// Interface calls an internal routine
|
||||
void DhopDeriv(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
|
||||
///////////////////////////////////////////////////////////////
|
||||
void Dhop(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);
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Multigrid assistance; force term uses too
|
||||
///////////////////////////////////////////////////////////////
|
||||
void Mdir(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);
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Extra methods added by derived
|
||||
///////////////////////////////////////////////////////////////
|
||||
void DerivInternal(StencilImpl &st, DoubledGaugeField &U, GaugeField &mat,
|
||||
const FermionField &A, const FermionField &B, int dag);
|
||||
|
||||
void DhopInternal(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
const FermionField &in, FermionField &out, int dag);
|
||||
|
||||
// Constructor
|
||||
WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid,
|
||||
GridRedBlackCartesian &Hgrid, RealD _mass,
|
||||
const ImplParams &p = ImplParams());
|
||||
|
||||
// DoubleStore impl dependent
|
||||
void ImportGauge(const GaugeField &_Umu);
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Data members require to support the functionality
|
||||
///////////////////////////////////////////////////////////////
|
||||
|
||||
// protected:
|
||||
public:
|
||||
RealD mass;
|
||||
|
||||
GridBase *_grid;
|
||||
GridBase *_cbgrid;
|
||||
|
||||
// Defines the stencils for even and odd
|
||||
StencilImpl Stencil;
|
||||
StencilImpl StencilEven;
|
||||
StencilImpl StencilOdd;
|
||||
|
||||
// Copy of the gauge field , with even and odd subsets
|
||||
DoubledGaugeField Umu;
|
||||
DoubledGaugeField UmuEven;
|
||||
DoubledGaugeField UmuOdd;
|
||||
|
||||
LebesgueOrder Lebesgue;
|
||||
LebesgueOrder LebesgueEvenOdd;
|
||||
};
|
||||
|
||||
typedef WilsonFermion<WilsonImplF> WilsonFermionF;
|
||||
typedef WilsonFermion<WilsonImplD> WilsonFermionD;
|
||||
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// non-hermitian hopping term; half cb or both
|
||||
///////////////////////////////////////////////////////////////
|
||||
void Dhop(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) ;
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Multigrid assistance; force term uses too
|
||||
///////////////////////////////////////////////////////////////
|
||||
void Mdir (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) ;
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Extra methods added by derived
|
||||
///////////////////////////////////////////////////////////////
|
||||
void DerivInternal(StencilImpl & st,
|
||||
DoubledGaugeField & U,
|
||||
GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag);
|
||||
|
||||
void DhopInternal(StencilImpl & st,LebesgueOrder & lo,DoubledGaugeField & U,
|
||||
const FermionField &in, FermionField &out,int dag) ;
|
||||
|
||||
// Constructor
|
||||
WilsonFermion(GaugeField &_Umu,
|
||||
GridCartesian &Fgrid,
|
||||
GridRedBlackCartesian &Hgrid,
|
||||
RealD _mass,
|
||||
const ImplParams &p= ImplParams()
|
||||
) ;
|
||||
|
||||
// DoubleStore impl dependent
|
||||
void ImportGauge(const GaugeField &_Umu);
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Data members require to support the functionality
|
||||
///////////////////////////////////////////////////////////////
|
||||
|
||||
// protected:
|
||||
public:
|
||||
|
||||
RealD mass;
|
||||
|
||||
GridBase * _grid;
|
||||
GridBase * _cbgrid;
|
||||
|
||||
//Defines the stencils for even and odd
|
||||
StencilImpl Stencil;
|
||||
StencilImpl StencilEven;
|
||||
StencilImpl StencilOdd;
|
||||
|
||||
// Copy of the gauge field , with even and odd subsets
|
||||
DoubledGaugeField Umu;
|
||||
DoubledGaugeField UmuEven;
|
||||
DoubledGaugeField UmuOdd;
|
||||
|
||||
LebesgueOrder Lebesgue;
|
||||
LebesgueOrder LebesgueEvenOdd;
|
||||
|
||||
|
||||
};
|
||||
|
||||
typedef WilsonFermion<WilsonImplF> WilsonFermionF;
|
||||
typedef WilsonFermion<WilsonImplD> WilsonFermionD;
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
@ -42,11 +42,11 @@ const std::vector<int> WilsonFermion5DStatic::displacements({1,1,1,1,-1,-1,-1,-1
|
||||
// 5d lattice for DWF.
|
||||
template<class Impl>
|
||||
WilsonFermion5D<Impl>::WilsonFermion5D(GaugeField &_Umu,
|
||||
GridCartesian &FiveDimGrid,
|
||||
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
||||
GridCartesian &FourDimGrid,
|
||||
GridRedBlackCartesian &FourDimRedBlackGrid,
|
||||
RealD _M5,const ImplParams &p) :
|
||||
GridCartesian &FiveDimGrid,
|
||||
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
||||
GridCartesian &FourDimGrid,
|
||||
GridRedBlackCartesian &FourDimRedBlackGrid,
|
||||
RealD _M5,const ImplParams &p) :
|
||||
Kernels(p),
|
||||
_FiveDimGrid (&FiveDimGrid),
|
||||
_FiveDimRedBlackGrid(&FiveDimRedBlackGrid),
|
||||
@ -135,10 +135,10 @@ WilsonFermion5D<Impl>::WilsonFermion5D(GaugeField &_Umu,
|
||||
/*
|
||||
template<class Impl>
|
||||
WilsonFermion5D<Impl>::WilsonFermion5D(int simd,GaugeField &_Umu,
|
||||
GridCartesian &FiveDimGrid,
|
||||
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
||||
GridCartesian &FourDimGrid,
|
||||
RealD _M5,const ImplParams &p) :
|
||||
GridCartesian &FiveDimGrid,
|
||||
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
||||
GridCartesian &FourDimGrid,
|
||||
RealD _M5,const ImplParams &p) :
|
||||
{
|
||||
int nsimd = Simd::Nsimd();
|
||||
|
||||
@ -175,6 +175,73 @@ WilsonFermion5D<Impl>::WilsonFermion5D(int simd,GaugeField &_Umu,
|
||||
}
|
||||
*/
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion5D<Impl>::Report(void)
|
||||
{
|
||||
std::vector<int> latt = GridDefaultLatt();
|
||||
RealD volume = Ls; for(int mu=0;mu<Nd;mu++) volume=volume*latt[mu];
|
||||
RealD NP = _FourDimGrid->_Nprocessors;
|
||||
|
||||
if ( DhopCalls > 0 ) {
|
||||
std::cout << GridLogMessage << "#### Dhop calls report " << std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D Number of Dhop Calls : " << DhopCalls << std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D Total Communication time : " << DhopCommTime
|
||||
<< " us" << std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D CommTime/Calls : "
|
||||
<< DhopCommTime / DhopCalls << " us" << std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D Total Compute time : "
|
||||
<< DhopComputeTime << " us" << std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D ComputeTime/Calls : "
|
||||
<< DhopComputeTime / DhopCalls << " us" << std::endl;
|
||||
|
||||
RealD mflops = 1344*volume*DhopCalls/DhopComputeTime;
|
||||
std::cout << GridLogMessage << "Average mflops/s per call : " << mflops << std::endl;
|
||||
std::cout << GridLogMessage << "Average mflops/s per call per node : " << mflops/NP << std::endl;
|
||||
|
||||
}
|
||||
|
||||
if ( DerivCalls > 0 ) {
|
||||
std::cout << GridLogMessage << "#### Deriv calls report "<< std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D Number of Deriv Calls : " <<DerivCalls <<std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D Total Communication time : " <<DerivCommTime <<" us"<<std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D CommTime/Calls : " <<DerivCommTime/DerivCalls<<" us" <<std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D Total Compute time : " <<DerivComputeTime <<" us"<<std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D ComputeTime/Calls : " <<DerivComputeTime/DerivCalls<<" us" <<std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D Total Dhop Compute time : " <<DerivDhopComputeTime <<" us"<<std::endl;
|
||||
std::cout << GridLogMessage << "WilsonFermion5D Dhop ComputeTime/Calls : " <<DerivDhopComputeTime/DerivCalls<<" us" <<std::endl;
|
||||
|
||||
|
||||
|
||||
RealD mflops = 144*volume*DerivCalls/DerivDhopComputeTime;
|
||||
std::cout << GridLogMessage << "Average mflops/s per call : " << mflops << std::endl;
|
||||
std::cout << GridLogMessage << "Average mflops/s per call per node : " << mflops/NP << std::endl;
|
||||
|
||||
}
|
||||
|
||||
if (DerivCalls > 0 || DhopCalls > 0){
|
||||
std::cout << GridLogMessage << "WilsonFermion5D Stencil"<<std::endl; Stencil.Report();
|
||||
std::cout << GridLogMessage << "WilsonFermion5D StencilEven"<<std::endl; StencilEven.Report();
|
||||
std::cout << GridLogMessage << "WilsonFermion5D StencilOdd"<<std::endl; StencilOdd.Report();
|
||||
}
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion5D<Impl>::ZeroCounters(void) {
|
||||
DhopCalls = 0;
|
||||
DhopCommTime = 0;
|
||||
DhopComputeTime = 0;
|
||||
|
||||
DerivCalls = 0;
|
||||
DerivCommTime = 0;
|
||||
DerivComputeTime = 0;
|
||||
DerivDhopComputeTime = 0;
|
||||
|
||||
Stencil.ZeroCounters();
|
||||
StencilEven.ZeroCounters();
|
||||
StencilOdd.ZeroCounters();
|
||||
}
|
||||
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion5D<Impl>::ImportGauge(const GaugeField &_Umu)
|
||||
{
|
||||
@ -215,12 +282,13 @@ PARALLEL_FOR_LOOP
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion5D<Impl>::DerivInternal(StencilImpl & st,
|
||||
DoubledGaugeField & U,
|
||||
GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
DoubledGaugeField & U,
|
||||
GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
{
|
||||
DerivCalls++;
|
||||
assert((dag==DaggerNo) ||(dag==DaggerYes));
|
||||
|
||||
conformable(st._grid,A._grid);
|
||||
@ -231,51 +299,53 @@ void WilsonFermion5D<Impl>::DerivInternal(StencilImpl & st,
|
||||
FermionField Btilde(B._grid);
|
||||
FermionField Atilde(B._grid);
|
||||
|
||||
DerivCommTime-=usecond();
|
||||
st.HaloExchange(B,compressor);
|
||||
DerivCommTime+=usecond();
|
||||
|
||||
Atilde=A;
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
|
||||
DerivComputeTime-=usecond();
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// Flip gamma if dag
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
int gamma = mu;
|
||||
if ( !dag ) gamma+= Nd;
|
||||
if (!dag) gamma += Nd;
|
||||
|
||||
////////////////////////
|
||||
// Call the single hop
|
||||
////////////////////////
|
||||
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int sss=0;sss<U._grid->oSites();sss++){
|
||||
for(int s=0;s<Ls;s++){
|
||||
int sU=sss;
|
||||
int sF = s+Ls*sU;
|
||||
DerivDhopComputeTime -= usecond();
|
||||
PARALLEL_FOR_LOOP
|
||||
for (int sss = 0; sss < U._grid->oSites(); sss++) {
|
||||
for (int s = 0; s < Ls; s++) {
|
||||
int sU = sss;
|
||||
int sF = s + Ls * sU;
|
||||
|
||||
assert ( sF< B._grid->oSites());
|
||||
assert ( sU< U._grid->oSites());
|
||||
assert(sF < B._grid->oSites());
|
||||
assert(sU < U._grid->oSites());
|
||||
|
||||
Kernels::DiracOptDhopDir(st,U,st.comm_buf,sF,sU,B,Btilde,mu,gamma);
|
||||
|
||||
////////////////////////////
|
||||
// spin trace outer product
|
||||
////////////////////////////
|
||||
Kernels::DiracOptDhopDir(st, U, st.comm_buf, sF, sU, B, Btilde, mu,
|
||||
gamma);
|
||||
|
||||
////////////////////////////
|
||||
// spin trace outer product
|
||||
////////////////////////////
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
Impl::InsertForce5D(mat,Btilde,Atilde,mu);
|
||||
|
||||
DerivDhopComputeTime += usecond();
|
||||
Impl::InsertForce5D(mat, Btilde, Atilde, mu);
|
||||
}
|
||||
DerivComputeTime += usecond();
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion5D<Impl>::DhopDeriv( GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
{
|
||||
conformable(A._grid,FermionGrid());
|
||||
conformable(A._grid,B._grid);
|
||||
@ -288,9 +358,9 @@ void WilsonFermion5D<Impl>::DhopDeriv( GaugeField &mat,
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion5D<Impl>::DhopDerivEO(GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
{
|
||||
conformable(A._grid,FermionRedBlackGrid());
|
||||
conformable(GaugeRedBlackGrid(),mat._grid);
|
||||
@ -306,9 +376,9 @@ void WilsonFermion5D<Impl>::DhopDerivEO(GaugeField &mat,
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion5D<Impl>::DhopDerivOE(GaugeField &mat,
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
const FermionField &A,
|
||||
const FermionField &B,
|
||||
int dag)
|
||||
{
|
||||
conformable(A._grid,FermionRedBlackGrid());
|
||||
conformable(GaugeRedBlackGrid(),mat._grid);
|
||||
@ -323,32 +393,61 @@ void WilsonFermion5D<Impl>::DhopDerivOE(GaugeField &mat,
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion5D<Impl>::DhopInternal(StencilImpl & st, LebesgueOrder &lo,
|
||||
DoubledGaugeField & U,
|
||||
const FermionField &in, FermionField &out,int dag)
|
||||
DoubledGaugeField & U,
|
||||
const FermionField &in, FermionField &out,int dag)
|
||||
{
|
||||
DhopCalls++;
|
||||
// assert((dag==DaggerNo) ||(dag==DaggerYes));
|
||||
Compressor compressor(dag);
|
||||
|
||||
int LLs = in._grid->_rdimensions[0];
|
||||
|
||||
DhopCommTime-=usecond();
|
||||
st.HaloExchange(in,compressor);
|
||||
DhopCommTime+=usecond();
|
||||
|
||||
DhopComputeTime-=usecond();
|
||||
// Dhop takes the 4d grid from U, and makes a 5d index for fermion
|
||||
if ( dag == DaggerYes ) {
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int ss=0;ss<U._grid->oSites();ss++){
|
||||
int sU=ss;
|
||||
int sF=LLs*sU;
|
||||
Kernels::DiracOptDhopSiteDag(st,lo,U,st.comm_buf,sF,sU,LLs,1,in,out);
|
||||
if (dag == DaggerYes) {
|
||||
PARALLEL_FOR_LOOP
|
||||
for (int ss = 0; ss < U._grid->oSites(); ss++) {
|
||||
int sU = ss;
|
||||
int sF = LLs * sU;
|
||||
Kernels::DiracOptDhopSiteDag(st, lo, U, st.comm_buf, sF, sU, LLs, 1, in,
|
||||
out);
|
||||
}
|
||||
#ifdef AVX512
|
||||
} else if (stat.is_init() ) {
|
||||
|
||||
int nthreads;
|
||||
stat.start();
|
||||
#pragma omp parallel
|
||||
{
|
||||
#pragma omp master
|
||||
nthreads = omp_get_num_threads();
|
||||
int mythread = omp_get_thread_num();
|
||||
stat.enter(mythread);
|
||||
#pragma omp for nowait
|
||||
for(int ss=0;ss<U._grid->oSites();ss++)
|
||||
{
|
||||
int sU=ss;
|
||||
int sF=LLs*sU;
|
||||
Kernels::DiracOptDhopSite(st,lo,U,st.comm_buf,sF,sU,LLs,1,in,out);
|
||||
}
|
||||
stat.exit(mythread);
|
||||
}
|
||||
stat.accum(nthreads);
|
||||
#endif
|
||||
} else {
|
||||
PARALLEL_FOR_LOOP
|
||||
for(int ss=0;ss<U._grid->oSites();ss++){
|
||||
int sU=ss;
|
||||
int sF=LLs*sU;
|
||||
Kernels::DiracOptDhopSite(st,lo,U,st.comm_buf,sF,sU,LLs,1,in,out);
|
||||
PARALLEL_FOR_LOOP
|
||||
for (int ss = 0; ss < U._grid->oSites(); ss++) {
|
||||
int sU = ss;
|
||||
int sF = LLs * sU;
|
||||
Kernels::DiracOptDhopSite(st, lo, U, st.comm_buf, sF, sU, LLs, 1, in,
|
||||
out);
|
||||
}
|
||||
}
|
||||
DhopComputeTime+=usecond();
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -31,6 +31,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
#ifndef GRID_QCD_WILSON_FERMION_5D_H
|
||||
#define GRID_QCD_WILSON_FERMION_5D_H
|
||||
|
||||
#include <Grid/Stat.h>
|
||||
|
||||
namespace Grid {
|
||||
|
||||
namespace QCD {
|
||||
@ -60,6 +62,18 @@ namespace Grid {
|
||||
public:
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
typedef WilsonKernels<Impl> Kernels;
|
||||
PmuStat stat;
|
||||
|
||||
void Report(void);
|
||||
void ZeroCounters(void);
|
||||
double DhopCalls;
|
||||
double DhopCommTime;
|
||||
double DhopComputeTime;
|
||||
|
||||
double DerivCalls;
|
||||
double DerivCommTime;
|
||||
double DerivComputeTime;
|
||||
double DerivDhopComputeTime;
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
// Implement the abstract base
|
||||
|
||||
@ -1,103 +1,54 @@
|
||||
/*************************************************************************************
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
|
||||
Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
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.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
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.,
|
||||
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 */
|
||||
See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid.h>
|
||||
namespace Grid {
|
||||
namespace QCD {
|
||||
|
||||
int WilsonKernelsStatic::HandOpt;
|
||||
int WilsonKernelsStatic::AsmOpt;
|
||||
int WilsonKernelsStatic::HandOpt;
|
||||
int WilsonKernelsStatic::AsmOpt;
|
||||
|
||||
template<class Impl>
|
||||
WilsonKernels<Impl>::WilsonKernels(const ImplParams &p): Base(p) {};
|
||||
template <class Impl>
|
||||
WilsonKernels<Impl>::WilsonKernels(const ImplParams &p) : Base(p){};
|
||||
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl>::DiracOptDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,int Ls, int Ns, const FermionField &in, FermionField &out)
|
||||
{
|
||||
#ifdef AVX512
|
||||
if ( AsmOpt ) {
|
||||
////////////////////////////////////////////
|
||||
// Generic implementation; move to different file?
|
||||
////////////////////////////////////////////
|
||||
|
||||
WilsonKernels<Impl>::DiracOptAsmDhopSite(st,lo,U,buf,sF,sU,Ls,Ns,in,out);
|
||||
|
||||
} else {
|
||||
#else
|
||||
{
|
||||
#endif
|
||||
for(int site=0;site<Ns;site++) {
|
||||
for(int s=0;s<Ls;s++) {
|
||||
if (HandOpt) WilsonKernels<Impl>::DiracOptHandDhopSite(st,lo,U,buf,sF,sU,in,out);
|
||||
else WilsonKernels<Impl>::DiracOptGenericDhopSite(st,lo,U,buf,sF,sU,in,out);
|
||||
sF++;
|
||||
}
|
||||
sU++;
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl>::DiracOptDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,int Ls, int Ns, const FermionField &in, FermionField &out)
|
||||
{
|
||||
#ifdef AVX512
|
||||
if ( AsmOpt ) {
|
||||
WilsonKernels<Impl>::DiracOptAsmDhopSiteDag(st,lo,U,buf,sF,sU,Ls,Ns,in,out);
|
||||
} else {
|
||||
#else
|
||||
{
|
||||
#endif
|
||||
for(int site=0;site<Ns;site++) {
|
||||
for(int s=0;s<Ls;s++) {
|
||||
if (HandOpt) WilsonKernels<Impl>::DiracOptHandDhopSiteDag(st,lo,U,buf,sF,sU,in,out);
|
||||
else WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(st,lo,U,buf,sF,sU,in,out);
|
||||
sF++;
|
||||
}
|
||||
sU++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
////////////////////////////////////////////
|
||||
// Generic implementation; move to different file?
|
||||
////////////////////////////////////////////
|
||||
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,const FermionField &in, FermionField &out)
|
||||
{
|
||||
SiteHalfSpinor tmp;
|
||||
SiteHalfSpinor chi;
|
||||
template <class Impl>
|
||||
void WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf, int sF,
|
||||
int sU, const FermionField &in, FermionField &out) {
|
||||
SiteHalfSpinor tmp;
|
||||
SiteHalfSpinor chi;
|
||||
SiteHalfSpinor *chi_p;
|
||||
SiteHalfSpinor Uchi;
|
||||
SiteSpinor result;
|
||||
@ -107,176 +58,175 @@ void WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(StencilImpl &st,LebesgueOrd
|
||||
///////////////////////////
|
||||
// Xp
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Xp,sF);
|
||||
SE = st.GetEntry(ptype, Xp, sF);
|
||||
|
||||
if (SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjXp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjXp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjXp(chi,in._odata[SE->_offset]);
|
||||
spProjXp(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
} else {
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Xp,SE,st);
|
||||
spReconXp(result,Uchi);
|
||||
|
||||
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Xp, SE, st);
|
||||
spReconXp(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Yp
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Yp,sF);
|
||||
SE = st.GetEntry(ptype, Yp, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjYp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjYp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjYp(chi,in._odata[SE->_offset]);
|
||||
spProjYp(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
} else {
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Yp,SE,st);
|
||||
accumReconYp(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Yp, SE, st);
|
||||
accumReconYp(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Zp
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Zp,sF);
|
||||
SE = st.GetEntry(ptype, Zp, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjZp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjZp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjZp(chi,in._odata[SE->_offset]);
|
||||
spProjZp(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
} else {
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Zp,SE,st);
|
||||
accumReconZp(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Zp, SE, st);
|
||||
accumReconZp(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Tp
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Tp,sF);
|
||||
SE = st.GetEntry(ptype, Tp, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjTp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjTp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjTp(chi,in._odata[SE->_offset]);
|
||||
spProjTp(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Tp,SE,st);
|
||||
accumReconTp(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Tp, SE, st);
|
||||
accumReconTp(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Xm
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Xm,sF);
|
||||
SE = st.GetEntry(ptype, Xm, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjXm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjXm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjXm(chi,in._odata[SE->_offset]);
|
||||
spProjXm(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Xm,SE,st);
|
||||
accumReconXm(result,Uchi);
|
||||
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Xm, SE, st);
|
||||
accumReconXm(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Ym
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Ym,sF);
|
||||
SE = st.GetEntry(ptype, Ym, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjYm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjYm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjYm(chi,in._odata[SE->_offset]);
|
||||
spProjYm(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Ym,SE,st);
|
||||
accumReconYm(result,Uchi);
|
||||
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Ym, SE, st);
|
||||
accumReconYm(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Zm
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Zm,sF);
|
||||
SE = st.GetEntry(ptype, Zm, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjZm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjZm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjZm(chi,in._odata[SE->_offset]);
|
||||
spProjZm(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Zm,SE,st);
|
||||
accumReconZm(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Zm, SE, st);
|
||||
accumReconZm(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Tm
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Tm,sF);
|
||||
SE = st.GetEntry(ptype, Tm, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjTm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
} else {
|
||||
spProjTm(chi,in._odata[SE->_offset]);
|
||||
if (SE->_permute) {
|
||||
spProjTm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjTm(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Tm,SE,st);
|
||||
accumReconTm(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Tm, SE, st);
|
||||
accumReconTm(result, Uchi);
|
||||
|
||||
vstream(out._odata[sF],result);
|
||||
vstream(out._odata[sF], result);
|
||||
};
|
||||
|
||||
|
||||
// Need controls to do interior, exterior, or both
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl>::DiracOptGenericDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,const FermionField &in, FermionField &out)
|
||||
{
|
||||
SiteHalfSpinor tmp;
|
||||
SiteHalfSpinor chi;
|
||||
SiteHalfSpinor *chi_p;
|
||||
// Need controls to do interior, exterior, or both
|
||||
template <class Impl>
|
||||
void WilsonKernels<Impl>::DiracOptGenericDhopSite(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf, int sF,
|
||||
int sU, const FermionField &in, FermionField &out) {
|
||||
SiteHalfSpinor tmp;
|
||||
SiteHalfSpinor chi;
|
||||
SiteHalfSpinor *chi_p;
|
||||
SiteHalfSpinor Uchi;
|
||||
SiteSpinor result;
|
||||
StencilEntry *SE;
|
||||
@ -285,296 +235,298 @@ void WilsonKernels<Impl>::DiracOptGenericDhopSite(StencilImpl &st,LebesgueOrder
|
||||
///////////////////////////
|
||||
// Xp
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Xm,sF);
|
||||
SE = st.GetEntry(ptype, Xm, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjXp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjXp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjXp(chi,in._odata[SE->_offset]);
|
||||
spProjXp(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
} else {
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Xm,SE,st);
|
||||
spReconXp(result,Uchi);
|
||||
|
||||
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Xm, SE, st);
|
||||
spReconXp(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Yp
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Ym,sF);
|
||||
SE = st.GetEntry(ptype, Ym, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjYp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjYp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjYp(chi,in._odata[SE->_offset]);
|
||||
spProjYp(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
} else {
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Ym,SE,st);
|
||||
accumReconYp(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Ym, SE, st);
|
||||
accumReconYp(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Zp
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Zm,sF);
|
||||
SE = st.GetEntry(ptype, Zm, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjZp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjZp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjZp(chi,in._odata[SE->_offset]);
|
||||
spProjZp(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
} else {
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Zm,SE,st);
|
||||
accumReconZp(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Zm, SE, st);
|
||||
accumReconZp(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Tp
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Tm,sF);
|
||||
SE = st.GetEntry(ptype, Tm, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjTp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjTp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjTp(chi,in._odata[SE->_offset]);
|
||||
spProjTp(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Tm,SE,st);
|
||||
accumReconTp(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Tm, SE, st);
|
||||
accumReconTp(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Xm
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Xp,sF);
|
||||
SE = st.GetEntry(ptype, Xp, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjXm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjXm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjXm(chi,in._odata[SE->_offset]);
|
||||
spProjXm(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Xp,SE,st);
|
||||
accumReconXm(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Xp, SE, st);
|
||||
accumReconXm(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Ym
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Yp,sF);
|
||||
SE = st.GetEntry(ptype, Yp, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjYm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjYm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjYm(chi,in._odata[SE->_offset]);
|
||||
spProjYm(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Yp,SE,st);
|
||||
accumReconYm(result,Uchi);
|
||||
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Yp, SE, st);
|
||||
accumReconYm(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Zm
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Zp,sF);
|
||||
SE = st.GetEntry(ptype, Zp, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjZm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
if (SE->_permute) {
|
||||
spProjZm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjZm(chi,in._odata[SE->_offset]);
|
||||
spProjZm(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Zp,SE,st);
|
||||
accumReconZm(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Zp, SE, st);
|
||||
accumReconZm(result, Uchi);
|
||||
|
||||
///////////////////////////
|
||||
// Tm
|
||||
///////////////////////////
|
||||
SE=st.GetEntry(ptype,Tp,sF);
|
||||
SE = st.GetEntry(ptype, Tp, sF);
|
||||
|
||||
if ( SE->_is_local ) {
|
||||
if (SE->_is_local) {
|
||||
chi_p = χ
|
||||
if ( SE->_permute ) {
|
||||
spProjTm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
} else {
|
||||
spProjTm(chi,in._odata[SE->_offset]);
|
||||
if (SE->_permute) {
|
||||
spProjTm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else {
|
||||
spProjTm(chi, in._odata[SE->_offset]);
|
||||
}
|
||||
} else {
|
||||
chi_p=&buf[SE->_offset];
|
||||
chi_p = &buf[SE->_offset];
|
||||
}
|
||||
|
||||
Impl::multLink(Uchi,U._odata[sU],*chi_p,Tp,SE,st);
|
||||
accumReconTm(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], *chi_p, Tp, SE, st);
|
||||
accumReconTm(result, Uchi);
|
||||
|
||||
vstream(out._odata[sF],result);
|
||||
vstream(out._odata[sF], result);
|
||||
};
|
||||
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl>::DiracOptDhopDir(StencilImpl &st,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,const FermionField &in, FermionField &out,int dir,int gamma)
|
||||
{
|
||||
SiteHalfSpinor tmp;
|
||||
SiteHalfSpinor chi;
|
||||
SiteSpinor result;
|
||||
template <class Impl>
|
||||
void WilsonKernels<Impl>::DiracOptDhopDir(
|
||||
StencilImpl &st, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf, int sF,
|
||||
int sU, const FermionField &in, FermionField &out, int dir, int gamma) {
|
||||
SiteHalfSpinor tmp;
|
||||
SiteHalfSpinor chi;
|
||||
SiteSpinor result;
|
||||
SiteHalfSpinor Uchi;
|
||||
StencilEntry *SE;
|
||||
int ptype;
|
||||
|
||||
SE=st.GetEntry(ptype,dir,sF);
|
||||
SE = st.GetEntry(ptype, dir, sF);
|
||||
|
||||
// Xp
|
||||
if(gamma==Xp){
|
||||
if ( SE->_is_local && SE->_permute ) {
|
||||
spProjXp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
} else if ( SE->_is_local ) {
|
||||
spProjXp(chi,in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi=buf[SE->_offset];
|
||||
if (gamma == Xp) {
|
||||
if (SE->_is_local && SE->_permute) {
|
||||
spProjXp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else if (SE->_is_local) {
|
||||
spProjXp(chi, in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi = buf[SE->_offset];
|
||||
}
|
||||
Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
|
||||
spReconXp(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
|
||||
spReconXp(result, Uchi);
|
||||
}
|
||||
|
||||
// Yp
|
||||
if ( gamma==Yp ){
|
||||
if ( SE->_is_local && SE->_permute ) {
|
||||
spProjYp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
} else if ( SE->_is_local ) {
|
||||
spProjYp(chi,in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi=buf[SE->_offset];
|
||||
if (gamma == Yp) {
|
||||
if (SE->_is_local && SE->_permute) {
|
||||
spProjYp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else if (SE->_is_local) {
|
||||
spProjYp(chi, in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi = buf[SE->_offset];
|
||||
}
|
||||
Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
|
||||
spReconYp(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
|
||||
spReconYp(result, Uchi);
|
||||
}
|
||||
|
||||
|
||||
// Zp
|
||||
if ( gamma ==Zp ){
|
||||
if ( SE->_is_local && SE->_permute ) {
|
||||
spProjZp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
} else if ( SE->_is_local ) {
|
||||
spProjZp(chi,in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi=buf[SE->_offset];
|
||||
if (gamma == Zp) {
|
||||
if (SE->_is_local && SE->_permute) {
|
||||
spProjZp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else if (SE->_is_local) {
|
||||
spProjZp(chi, in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi = buf[SE->_offset];
|
||||
}
|
||||
Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
|
||||
spReconZp(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
|
||||
spReconZp(result, Uchi);
|
||||
}
|
||||
|
||||
|
||||
// Tp
|
||||
if ( gamma ==Tp ){
|
||||
if ( SE->_is_local && SE->_permute ) {
|
||||
spProjTp(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
} else if ( SE->_is_local ) {
|
||||
spProjTp(chi,in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi=buf[SE->_offset];
|
||||
if (gamma == Tp) {
|
||||
if (SE->_is_local && SE->_permute) {
|
||||
spProjTp(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else if (SE->_is_local) {
|
||||
spProjTp(chi, in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi = buf[SE->_offset];
|
||||
}
|
||||
Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
|
||||
spReconTp(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
|
||||
spReconTp(result, Uchi);
|
||||
}
|
||||
|
||||
// Xm
|
||||
if ( gamma==Xm ){
|
||||
if ( SE->_is_local && SE->_permute ) {
|
||||
spProjXm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
} else if ( SE->_is_local ) {
|
||||
spProjXm(chi,in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi=buf[SE->_offset];
|
||||
if (gamma == Xm) {
|
||||
if (SE->_is_local && SE->_permute) {
|
||||
spProjXm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else if (SE->_is_local) {
|
||||
spProjXm(chi, in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi = buf[SE->_offset];
|
||||
}
|
||||
Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
|
||||
spReconXm(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
|
||||
spReconXm(result, Uchi);
|
||||
}
|
||||
|
||||
// Ym
|
||||
if ( gamma == Ym ){
|
||||
if ( SE->_is_local && SE->_permute ) {
|
||||
spProjYm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
} else if ( SE->_is_local ) {
|
||||
spProjYm(chi,in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi=buf[SE->_offset];
|
||||
if (gamma == Ym) {
|
||||
if (SE->_is_local && SE->_permute) {
|
||||
spProjYm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else if (SE->_is_local) {
|
||||
spProjYm(chi, in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi = buf[SE->_offset];
|
||||
}
|
||||
Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
|
||||
spReconYm(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
|
||||
spReconYm(result, Uchi);
|
||||
}
|
||||
|
||||
// Zm
|
||||
if ( gamma == Zm ){
|
||||
if ( SE->_is_local && SE->_permute ) {
|
||||
spProjZm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
} else if ( SE->_is_local ) {
|
||||
spProjZm(chi,in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi=buf[SE->_offset];
|
||||
if (gamma == Zm) {
|
||||
if (SE->_is_local && SE->_permute) {
|
||||
spProjZm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else if (SE->_is_local) {
|
||||
spProjZm(chi, in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi = buf[SE->_offset];
|
||||
}
|
||||
Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
|
||||
spReconZm(result,Uchi);
|
||||
}
|
||||
|
||||
// Tm
|
||||
if ( gamma==Tm ) {
|
||||
if ( SE->_is_local && SE->_permute ) {
|
||||
spProjTm(tmp,in._odata[SE->_offset]);
|
||||
permute(chi,tmp,ptype);
|
||||
} else if ( SE->_is_local ) {
|
||||
spProjTm(chi,in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi=buf[SE->_offset];
|
||||
}
|
||||
Impl::multLink(Uchi,U._odata[sU],chi,dir,SE,st);
|
||||
spReconTm(result,Uchi);
|
||||
Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
|
||||
spReconZm(result, Uchi);
|
||||
}
|
||||
|
||||
vstream(out._odata[sF],result);
|
||||
// Tm
|
||||
if (gamma == Tm) {
|
||||
if (SE->_is_local && SE->_permute) {
|
||||
spProjTm(tmp, in._odata[SE->_offset]);
|
||||
permute(chi, tmp, ptype);
|
||||
} else if (SE->_is_local) {
|
||||
spProjTm(chi, in._odata[SE->_offset]);
|
||||
} else {
|
||||
chi = buf[SE->_offset];
|
||||
}
|
||||
Impl::multLink(Uchi, U._odata[sU], chi, dir, SE, st);
|
||||
spReconTm(result, Uchi);
|
||||
}
|
||||
|
||||
vstream(out._odata[sF], result);
|
||||
}
|
||||
|
||||
|
||||
FermOpTemplateInstantiate(WilsonKernels);
|
||||
FermOpTemplateInstantiate(WilsonKernels);
|
||||
AdjointFermOpTemplateInstantiate(WilsonKernels);
|
||||
TwoIndexFermOpTemplateInstantiate(WilsonKernels);
|
||||
|
||||
}}
|
||||
|
||||
|
||||
@ -1,34 +1,35 @@
|
||||
/*************************************************************************************
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/qcd/action/fermion/WilsonKernels.h
|
||||
Source file: ./lib/qcd/action/fermion/WilsonKernels.h
|
||||
|
||||
Copyright (C) 2015
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
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.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
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.,
|
||||
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 */
|
||||
#ifndef GRID_QCD_DHOP_H
|
||||
#define GRID_QCD_DHOP_H
|
||||
See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#ifndef GRID_QCD_DHOP_H
|
||||
#define GRID_QCD_DHOP_H
|
||||
|
||||
namespace Grid {
|
||||
|
||||
@ -48,55 +49,158 @@ namespace Grid {
|
||||
template<class Impl> class WilsonKernels : public FermionOperator<Impl> , public WilsonKernelsStatic {
|
||||
public:
|
||||
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
typedef FermionOperator<Impl> Base;
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
typedef FermionOperator<Impl> Base;
|
||||
|
||||
public:
|
||||
|
||||
void DiracOptDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU,int Ls, int Ns, const FermionField &in, FermionField &out);
|
||||
|
||||
void DiracOptDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,int Ls, int Ns, const FermionField &in,FermionField &out);
|
||||
template <bool EnableBool = true>
|
||||
typename std::enable_if<Impl::Dimension == 3 && Nc == 3 &&EnableBool, void>::type
|
||||
DiracOptDhopSite(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU, int Ls, int Ns, const FermionField &in,
|
||||
FermionField &out) {
|
||||
#ifdef AVX512
|
||||
if (AsmOpt) {
|
||||
WilsonKernels<Impl>::DiracOptAsmDhopSite(st, lo, U, buf, sF, sU, Ls, Ns,
|
||||
in, out);
|
||||
|
||||
} else {
|
||||
#else
|
||||
{
|
||||
#endif
|
||||
for (int site = 0; site < Ns; site++) {
|
||||
for (int s = 0; s < Ls; s++) {
|
||||
if (HandOpt)
|
||||
WilsonKernels<Impl>::DiracOptHandDhopSite(st, lo, U, buf, sF, sU,
|
||||
in, out);
|
||||
else
|
||||
WilsonKernels<Impl>::DiracOptGenericDhopSite(st, lo, U, buf, sF, sU,
|
||||
in, out);
|
||||
sF++;
|
||||
}
|
||||
sU++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <bool EnableBool = true>
|
||||
typename std::enable_if<(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool, void>::type
|
||||
DiracOptDhopSite(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU, int Ls, int Ns, const FermionField &in,
|
||||
FermionField &out) {
|
||||
for (int site = 0; site < Ns; site++) {
|
||||
for (int s = 0; s < Ls; s++) {
|
||||
WilsonKernels<Impl>::DiracOptGenericDhopSite(st, lo, U, buf, sF, sU, in,
|
||||
out);
|
||||
sF++;
|
||||
}
|
||||
sU++;
|
||||
}
|
||||
}
|
||||
|
||||
template <bool EnableBool = true>
|
||||
typename std::enable_if<Impl::Dimension == 3 && Nc == 3 && EnableBool,
|
||||
void>::type
|
||||
DiracOptDhopSiteDag(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU, int Ls, int Ns, const FermionField &in,
|
||||
FermionField &out) {
|
||||
#ifdef AVX512
|
||||
if (AsmOpt) {
|
||||
WilsonKernels<Impl>::DiracOptAsmDhopSiteDag(st, lo, U, buf, sF, sU, Ls,
|
||||
Ns, in, out);
|
||||
} else {
|
||||
#else
|
||||
{
|
||||
#endif
|
||||
for (int site = 0; site < Ns; site++) {
|
||||
for (int s = 0; s < Ls; s++) {
|
||||
if (HandOpt)
|
||||
WilsonKernels<Impl>::DiracOptHandDhopSiteDag(st, lo, U, buf, sF, sU,
|
||||
in, out);
|
||||
else
|
||||
WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(st, lo, U, buf, sF,
|
||||
sU, in, out);
|
||||
sF++;
|
||||
}
|
||||
sU++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <bool EnableBool = true>
|
||||
typename std::enable_if<
|
||||
(Impl::Dimension != 3 || (Impl::Dimension == 3 && Nc != 3)) && EnableBool,
|
||||
void>::type
|
||||
DiracOptDhopSiteDag(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU, int Ls, int Ns, const FermionField &in,
|
||||
FermionField &out) {
|
||||
for (int site = 0; site < Ns; site++) {
|
||||
for (int s = 0; s < Ls; s++) {
|
||||
WilsonKernels<Impl>::DiracOptGenericDhopSiteDag(st, lo, U, buf, sF, sU,
|
||||
in, out);
|
||||
sF++;
|
||||
}
|
||||
sU++;
|
||||
}
|
||||
}
|
||||
|
||||
void DiracOptDhopDir(
|
||||
StencilImpl &st, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU, const FermionField &in, FermionField &out, int dirdisp,
|
||||
int gamma);
|
||||
|
||||
private:
|
||||
// Specialised variants
|
||||
void DiracOptGenericDhopSite(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU, const FermionField &in, FermionField &out);
|
||||
|
||||
void DiracOptGenericDhopSiteDag(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU, const FermionField &in, FermionField &out);
|
||||
|
||||
void DiracOptAsmDhopSite(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU, int Ls, int Ns, const FermionField &in,
|
||||
FermionField &out);
|
||||
|
||||
void DiracOptAsmDhopSiteDag(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU, int Ls, int Ns, const FermionField &in,
|
||||
FermionField &out);
|
||||
|
||||
void DiracOptHandDhopSite(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU, const FermionField &in, FermionField &out);
|
||||
|
||||
void DiracOptHandDhopSiteDag(
|
||||
StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor, alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF, int sU, const FermionField &in, FermionField &out);
|
||||
|
||||
public:
|
||||
WilsonKernels(const ImplParams &p = ImplParams());
|
||||
};
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void DiracOptDhopDir(StencilImpl &st,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,const FermionField &in, FermionField &out,int dirdisp,int gamma);
|
||||
|
||||
private:
|
||||
// Specialised variants
|
||||
void DiracOptGenericDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU, const FermionField &in, FermionField &out);
|
||||
|
||||
void DiracOptGenericDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,const FermionField &in,FermionField &out);
|
||||
|
||||
void DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,int Ls, int Ns, const FermionField &in, FermionField &out);
|
||||
|
||||
void DiracOptAsmDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,int Ls, int Ns, const FermionField &in, FermionField &out);
|
||||
|
||||
|
||||
void DiracOptHandDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,const FermionField &in, FermionField &out);
|
||||
|
||||
void DiracOptHandDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int sF,int sU,const FermionField &in, FermionField &out);
|
||||
public:
|
||||
|
||||
WilsonKernels(const ImplParams &p= ImplParams());
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
/*************************************************************************************
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
@ -26,68 +26,77 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
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 */
|
||||
|
||||
#include <Grid.h>
|
||||
|
||||
|
||||
namespace Grid {
|
||||
namespace QCD {
|
||||
namespace QCD {
|
||||
|
||||
///////////////////////////////////////////////////////////
|
||||
// Default to no assembler implementation
|
||||
///////////////////////////////////////////////////////////
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl >::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
|
||||
{
|
||||
assert(0);
|
||||
}
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl >::DiracOptAsmDhopSiteDag(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
|
||||
{
|
||||
assert(0);
|
||||
}
|
||||
|
||||
|
||||
///////////////////////////////////////////////////////////
|
||||
// Default to no assembler implementation
|
||||
///////////////////////////////////////////////////////////
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl >::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
|
||||
{
|
||||
assert(0);
|
||||
}
|
||||
|
||||
#if defined(AVX512)
|
||||
|
||||
|
||||
///////////////////////////////////////////////////////////
|
||||
// If we are AVX512 specialise the single precision routine
|
||||
///////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
|
||||
///////////////////////////////////////////////////////////
|
||||
// If we are AVX512 specialise the single precision routine
|
||||
///////////////////////////////////////////////////////////
|
||||
|
||||
#include <simd/Intel512wilson.h>
|
||||
#include <simd/Intel512single.h>
|
||||
|
||||
static Vector<vComplexF> signs;
|
||||
|
||||
int setupSigns(void ){
|
||||
Vector<vComplexF> bother(2);
|
||||
signs = bother;
|
||||
vrsign(signs[0]);
|
||||
visign(signs[1]);
|
||||
return 1;
|
||||
}
|
||||
static int signInit = setupSigns();
|
||||
|
||||
|
||||
static Vector<vComplexF> signs;
|
||||
|
||||
int setupSigns(void ){
|
||||
Vector<vComplexF> bother(2);
|
||||
signs = bother;
|
||||
vrsign(signs[0]);
|
||||
visign(signs[1]);
|
||||
return 1;
|
||||
}
|
||||
static int signInit = setupSigns();
|
||||
|
||||
#define label(A) ilabel(A)
|
||||
#define ilabel(A) ".globl\n" #A ":\n"
|
||||
|
||||
|
||||
#define MAYBEPERM(A,perm) if (perm) { A ; }
|
||||
#define MULT_2SPIN(ptr,pf) MULT_ADDSUB_2SPIN(ptr,pf)
|
||||
#define FX(A) WILSONASM_ ##A
|
||||
|
||||
|
||||
#undef KERNEL_DAG
|
||||
template<>
|
||||
void WilsonKernels<WilsonImplF>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
|
||||
template<>
|
||||
void WilsonKernels<WilsonImplF>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
|
||||
#include <qcd/action/fermion/WilsonKernelsAsmBody.h>
|
||||
|
||||
|
||||
#define KERNEL_DAG
|
||||
template<>
|
||||
void WilsonKernels<WilsonImplF>::DiracOptAsmDhopSiteDag(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
|
||||
template<>
|
||||
void WilsonKernels<WilsonImplF>::DiracOptAsmDhopSiteDag(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
|
||||
#include <qcd/action/fermion/WilsonKernelsAsmBody.h>
|
||||
|
||||
|
||||
#undef VMOVIDUP
|
||||
#undef VMOVRDUP
|
||||
#undef MAYBEPERM
|
||||
@ -98,43 +107,43 @@ void WilsonKernels<WilsonImplF>::DiracOptAsmDhopSiteDag(StencilImpl &st,Lebesgue
|
||||
#define VMOVIDUP(A,B,C) VBCASTIDUPf(A,B,C)
|
||||
#define VMOVRDUP(A,B,C) VBCASTRDUPf(A,B,C)
|
||||
#define MULT_2SPIN(ptr,pf) MULT_ADDSUB_2SPIN_LS(ptr,pf)
|
||||
|
||||
|
||||
#undef KERNEL_DAG
|
||||
template<>
|
||||
void WilsonKernels<DomainWallVec5dImplF>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
|
||||
template<>
|
||||
void WilsonKernels<DomainWallVec5dImplF>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
|
||||
#include <qcd/action/fermion/WilsonKernelsAsmBody.h>
|
||||
|
||||
|
||||
#define KERNEL_DAG
|
||||
template<>
|
||||
void WilsonKernels<DomainWallVec5dImplF>::DiracOptAsmDhopSiteDag(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
|
||||
template<>
|
||||
void WilsonKernels<DomainWallVec5dImplF>::DiracOptAsmDhopSiteDag(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out)
|
||||
#include <qcd/action/fermion/WilsonKernelsAsmBody.h>
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
#define INSTANTIATE_ASM(A)\
|
||||
template void WilsonKernels<A>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,\
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,\
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out);\
|
||||
template void WilsonKernels<A>::DiracOptAsmDhopSiteDag(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,\
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,\
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out);\
|
||||
|
||||
template void WilsonKernels<WilsonImplF>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out);
|
||||
|
||||
template void WilsonKernels<WilsonImplD>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out);
|
||||
template void WilsonKernels<GparityWilsonImplF>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out);
|
||||
template void WilsonKernels<GparityWilsonImplD>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out);
|
||||
template void WilsonKernels<DomainWallVec5dImplF>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out);
|
||||
template void WilsonKernels<DomainWallVec5dImplD>::DiracOptAsmDhopSite(StencilImpl &st,LebesgueOrder & lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int ssU,int Ls,int Ns,const FermionField &in, FermionField &out);
|
||||
}}
|
||||
INSTANTIATE_ASM(WilsonImplF);
|
||||
INSTANTIATE_ASM(WilsonImplD);
|
||||
INSTANTIATE_ASM(ZWilsonImplF);
|
||||
INSTANTIATE_ASM(ZWilsonImplD);
|
||||
INSTANTIATE_ASM(GparityWilsonImplF);
|
||||
INSTANTIATE_ASM(GparityWilsonImplD);
|
||||
INSTANTIATE_ASM(DomainWallVec5dImplF);
|
||||
INSTANTIATE_ASM(DomainWallVec5dImplD);
|
||||
INSTANTIATE_ASM(ZDomainWallVec5dImplF);
|
||||
INSTANTIATE_ASM(ZDomainWallVec5dImplD);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -134,7 +134,9 @@
|
||||
////////////////////////////////
|
||||
// Xm
|
||||
////////////////////////////////
|
||||
#ifndef STREAM_STORE
|
||||
basep= (uint64_t) &out._odata[ss];
|
||||
#endif
|
||||
// basep= st.GetPFInfo(nent,plocal); nent++;
|
||||
if ( local ) {
|
||||
LOAD64(%r10,isigns); // times i => shuffle and xor the real part sign bit
|
||||
@ -229,7 +231,9 @@
|
||||
LOAD_CHI(base);
|
||||
}
|
||||
base= (uint64_t) &out._odata[ss];
|
||||
#ifndef STREAM_STORE
|
||||
PREFETCH_CHIMU(base);
|
||||
#endif
|
||||
{
|
||||
MULT_2SPIN_DIR_PFTM(Tm,basep);
|
||||
}
|
||||
|
||||
@ -311,8 +311,8 @@ namespace Grid {
|
||||
namespace QCD {
|
||||
|
||||
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl >::DiracOptHandDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl>::DiracOptHandDhopSite(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int sU,const FermionField &in, FermionField &out)
|
||||
{
|
||||
@ -554,8 +554,8 @@ void WilsonKernels<Impl >::DiracOptHandDhopSite(StencilImpl &st,LebesgueOrder &l
|
||||
}
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl >::DiracOptHandDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
template<class Impl>
|
||||
void WilsonKernels<Impl>::DiracOptHandDhopSiteDag(StencilImpl &st,LebesgueOrder &lo,DoubledGaugeField &U,
|
||||
std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> > &buf,
|
||||
int ss,int sU,const FermionField &in, FermionField &out)
|
||||
{
|
||||
|
||||
@ -30,7 +30,6 @@ directory
|
||||
#define GRID_QCD_GAUGE_IMPL_H
|
||||
|
||||
namespace Grid {
|
||||
|
||||
namespace QCD {
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
@ -52,7 +51,7 @@ public:
|
||||
typedef S Simd;
|
||||
|
||||
template <typename vtype>
|
||||
using iImplGaugeLink = iScalar<iScalar<iMatrix<vtype, Nrepresentation>>>;
|
||||
using iImplGaugeLink = iScalar<iScalar<iMatrix<vtype, Nrepresentation>>>;
|
||||
template <typename vtype>
|
||||
using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation>>, Nd>;
|
||||
|
||||
@ -64,7 +63,7 @@ public:
|
||||
// ugly
|
||||
typedef Lattice<SiteGaugeField> GaugeField;
|
||||
|
||||
// Move this elsewhere?
|
||||
// Move this elsewhere? FIXME
|
||||
static inline void AddGaugeLink(GaugeField &U, GaugeLinkField &W,
|
||||
int mu) { // U[mu] += W
|
||||
PARALLEL_FOR_LOOP
|
||||
@ -174,12 +173,19 @@ typedef GaugeImplTypes<vComplex, Nc> GimplTypesR;
|
||||
typedef GaugeImplTypes<vComplexF, Nc> GimplTypesF;
|
||||
typedef GaugeImplTypes<vComplexD, Nc> GimplTypesD;
|
||||
|
||||
typedef GaugeImplTypes<vComplex, SU<Nc>::AdjointDimension> GimplAdjointTypesR;
|
||||
typedef GaugeImplTypes<vComplexF, SU<Nc>::AdjointDimension> GimplAdjointTypesF;
|
||||
typedef GaugeImplTypes<vComplexD, SU<Nc>::AdjointDimension> GimplAdjointTypesD;
|
||||
|
||||
typedef PeriodicGaugeImpl<GimplTypesR> PeriodicGimplR; // Real.. whichever prec
|
||||
typedef PeriodicGaugeImpl<GimplTypesF> PeriodicGimplF; // Float
|
||||
typedef PeriodicGaugeImpl<GimplTypesD> PeriodicGimplD; // Double
|
||||
|
||||
typedef ConjugateGaugeImpl<GimplTypesR>
|
||||
ConjugateGimplR; // Real.. whichever prec
|
||||
typedef PeriodicGaugeImpl<GimplAdjointTypesR> PeriodicGimplAdjR; // Real.. whichever prec
|
||||
typedef PeriodicGaugeImpl<GimplAdjointTypesF> PeriodicGimplAdjF; // Float
|
||||
typedef PeriodicGaugeImpl<GimplAdjointTypesD> PeriodicGimplAdjD; // Double
|
||||
|
||||
typedef ConjugateGaugeImpl<GimplTypesR> ConjugateGimplR; // Real.. whichever prec
|
||||
typedef ConjugateGaugeImpl<GimplTypesF> ConjugateGimplF; // Float
|
||||
typedef ConjugateGaugeImpl<GimplTypesD> ConjugateGimplD; // Double
|
||||
}
|
||||
|
||||
@ -1,149 +1,151 @@
|
||||
/*************************************************************************************
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/qcd/action/pseudofermion/TwoFlavour.h
|
||||
Source file: ./lib/qcd/action/pseudofermion/TwoFlavour.h
|
||||
|
||||
Copyright (C) 2015
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
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.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
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.,
|
||||
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 */
|
||||
See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_H
|
||||
#define QCD_PSEUDOFERMION_TWO_FLAVOUR_H
|
||||
|
||||
namespace Grid{
|
||||
namespace QCD{
|
||||
namespace Grid {
|
||||
namespace QCD {
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// Two flavour pseudofermion action for any dop
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
template<class Impl>
|
||||
class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
|
||||
public:
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// Two flavour pseudofermion action for any dop
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
template <class Impl>
|
||||
class TwoFlavourPseudoFermionAction : public Action<typename Impl::GaugeField> {
|
||||
public:
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
|
||||
private:
|
||||
|
||||
FermionOperator<Impl> & FermOp;// the basic operator
|
||||
private:
|
||||
FermionOperator<Impl> &FermOp; // the basic operator
|
||||
|
||||
OperatorFunction<FermionField> &DerivativeSolver;
|
||||
OperatorFunction<FermionField> &DerivativeSolver;
|
||||
|
||||
OperatorFunction<FermionField> &ActionSolver;
|
||||
OperatorFunction<FermionField> &ActionSolver;
|
||||
|
||||
FermionField Phi; // the pseudo fermion field for this trajectory
|
||||
FermionField Phi; // the pseudo fermion field for this trajectory
|
||||
|
||||
public:
|
||||
/////////////////////////////////////////////////
|
||||
// Pass in required objects.
|
||||
/////////////////////////////////////////////////
|
||||
TwoFlavourPseudoFermionAction(FermionOperator<Impl> &Op,
|
||||
OperatorFunction<FermionField> & DS,
|
||||
OperatorFunction<FermionField> & AS
|
||||
) : FermOp(Op), DerivativeSolver(DS), ActionSolver(AS), Phi(Op.FermionGrid()) {
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
// Push the gauge field in to the dops. Assume any BC's and smearing already applied
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
|
||||
public:
|
||||
/////////////////////////////////////////////////
|
||||
// Pass in required objects.
|
||||
/////////////////////////////////////////////////
|
||||
TwoFlavourPseudoFermionAction(FermionOperator<Impl> &Op,
|
||||
OperatorFunction<FermionField> &DS,
|
||||
OperatorFunction<FermionField> &AS)
|
||||
: FermOp(Op),
|
||||
DerivativeSolver(DS),
|
||||
ActionSolver(AS),
|
||||
Phi(Op.FermionGrid()){};
|
||||
|
||||
// P(phi) = e^{- phi^dag (MdagM)^-1 phi}
|
||||
// Phi = Mdag eta
|
||||
// P(eta) = e^{- eta^dag eta}
|
||||
//
|
||||
// e^{x^2/2 sig^2} => sig^2 = 0.5.
|
||||
//
|
||||
// So eta should be of width sig = 1/sqrt(2).
|
||||
// and must multiply by 0.707....
|
||||
//
|
||||
// Chroma has this scale factor: two_flavor_monomial_w.h
|
||||
// IroIro: does not use this scale. It is absorbed by a change of vars
|
||||
// in the Phi integral, and thus is only an irrelevant prefactor for the partition function.
|
||||
//
|
||||
RealD scale = std::sqrt(0.5);
|
||||
FermionField eta(FermOp.FermionGrid());
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
// Push the gauge field in to the dops. Assume any BC's and smearing already
|
||||
// applied
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
virtual void refresh(const GaugeField &U, GridParallelRNG &pRNG) {
|
||||
// P(phi) = e^{- phi^dag (MdagM)^-1 phi}
|
||||
// Phi = Mdag eta
|
||||
// P(eta) = e^{- eta^dag eta}
|
||||
//
|
||||
// e^{x^2/2 sig^2} => sig^2 = 0.5.
|
||||
//
|
||||
// So eta should be of width sig = 1/sqrt(2).
|
||||
// and must multiply by 0.707....
|
||||
//
|
||||
// Chroma has this scale factor: two_flavor_monomial_w.h
|
||||
// IroIro: does not use this scale. It is absorbed by a change of vars
|
||||
// in the Phi integral, and thus is only an irrelevant prefactor for
|
||||
// the partition function.
|
||||
//
|
||||
RealD scale = std::sqrt(0.5);
|
||||
FermionField eta(FermOp.FermionGrid());
|
||||
|
||||
gaussian(pRNG,eta);
|
||||
gaussian(pRNG, eta);
|
||||
|
||||
FermOp.ImportGauge(U);
|
||||
FermOp.Mdag(eta,Phi);
|
||||
FermOp.ImportGauge(U);
|
||||
FermOp.Mdag(eta, Phi);
|
||||
|
||||
Phi=Phi*scale;
|
||||
|
||||
};
|
||||
Phi = Phi * scale;
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////
|
||||
// S = phi^dag (Mdag M)^-1 phi
|
||||
//////////////////////////////////////////////////////
|
||||
virtual RealD S(const GaugeField &U) {
|
||||
//////////////////////////////////////////////////////
|
||||
// S = phi^dag (Mdag M)^-1 phi
|
||||
//////////////////////////////////////////////////////
|
||||
virtual RealD S(const GaugeField &U) {
|
||||
FermOp.ImportGauge(U);
|
||||
|
||||
FermOp.ImportGauge(U);
|
||||
FermionField X(FermOp.FermionGrid());
|
||||
FermionField Y(FermOp.FermionGrid());
|
||||
|
||||
FermionField X(FermOp.FermionGrid());
|
||||
FermionField Y(FermOp.FermionGrid());
|
||||
|
||||
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
|
||||
X=zero;
|
||||
ActionSolver(MdagMOp,Phi,X);
|
||||
MdagMOp.Op(X,Y);
|
||||
MdagMLinearOperator<FermionOperator<Impl>, FermionField> MdagMOp(FermOp);
|
||||
X = zero;
|
||||
ActionSolver(MdagMOp, Phi, X);
|
||||
MdagMOp.Op(X, Y);
|
||||
|
||||
RealD action = norm2(Y);
|
||||
std::cout << GridLogMessage << "Pseudofermion action "<<action<<std::endl;
|
||||
return action;
|
||||
};
|
||||
RealD action = norm2(Y);
|
||||
std::cout << GridLogMessage << "Pseudofermion action " << action
|
||||
<< std::endl;
|
||||
return action;
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////
|
||||
// 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
|
||||
//
|
||||
// = - Ydag dM X - Xdag dMdag Y
|
||||
//
|
||||
//////////////////////////////////////////////////////
|
||||
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
|
||||
//////////////////////////////////////////////////////
|
||||
// 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
|
||||
//
|
||||
// = - Ydag dM X - Xdag dMdag Y
|
||||
//
|
||||
//////////////////////////////////////////////////////
|
||||
virtual void deriv(const GaugeField &U, GaugeField &dSdU) {
|
||||
FermOp.ImportGauge(U);
|
||||
|
||||
FermOp.ImportGauge(U);
|
||||
FermionField X(FermOp.FermionGrid());
|
||||
FermionField Y(FermOp.FermionGrid());
|
||||
GaugeField tmp(FermOp.GaugeGrid());
|
||||
|
||||
FermionField X(FermOp.FermionGrid());
|
||||
FermionField Y(FermOp.FermionGrid());
|
||||
GaugeField tmp(FermOp.GaugeGrid());
|
||||
MdagMLinearOperator<FermionOperator<Impl>, FermionField> MdagMOp(FermOp);
|
||||
|
||||
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(FermOp);
|
||||
X = zero;
|
||||
DerivativeSolver(MdagMOp, Phi, X); // X = (MdagM)^-1 phi
|
||||
MdagMOp.Op(X, Y); // Y = M X = (Mdag)^-1 phi
|
||||
|
||||
X=zero;
|
||||
DerivativeSolver(MdagMOp,Phi,X);
|
||||
MdagMOp.Op(X,Y);
|
||||
// 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.
|
||||
|
||||
// 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.
|
||||
FermOp.MDeriv(tmp, Y, X, DaggerNo);
|
||||
dSdU = tmp;
|
||||
FermOp.MDeriv(tmp, X, Y, DaggerYes);
|
||||
dSdU = dSdU + tmp;
|
||||
|
||||
FermOp.MDeriv(tmp , Y, X,DaggerNo ); dSdU=tmp;
|
||||
FermOp.MDeriv(tmp , X, Y,DaggerYes); dSdU=dSdU+tmp;
|
||||
|
||||
//dSdU = Ta(dSdU);
|
||||
|
||||
};
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
// not taking here the traceless antihermitian component
|
||||
};
|
||||
};
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
@ -1,70 +1,66 @@
|
||||
/*************************************************************************************
|
||||
/*************************************************************************************
|
||||
|
||||
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 <paboyle@ph.ed.ac.uk>
|
||||
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
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
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
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.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
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.,
|
||||
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 */
|
||||
See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#ifndef QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H
|
||||
#define QCD_PSEUDOFERMION_TWO_FLAVOUR_EVEN_ODD_H
|
||||
|
||||
namespace Grid{
|
||||
namespace QCD{
|
||||
namespace Grid {
|
||||
namespace QCD {
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// Two flavour pseudofermion action for any EO prec dop
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
template <class Impl>
|
||||
class TwoFlavourEvenOddPseudoFermionAction
|
||||
: public Action<typename Impl::GaugeField> {
|
||||
public:
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
|
||||
private:
|
||||
FermionOperator<Impl> &FermOp; // the basic operator
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// Two flavour pseudofermion action for any EO prec dop
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
template<class Impl>
|
||||
class TwoFlavourEvenOddPseudoFermionAction : public Action<typename Impl::GaugeField> {
|
||||
OperatorFunction<FermionField> &DerivativeSolver;
|
||||
OperatorFunction<FermionField> &ActionSolver;
|
||||
|
||||
public:
|
||||
FermionField PhiOdd; // the pseudo fermion field for this trajectory
|
||||
FermionField PhiEven; // the pseudo fermion field for this trajectory
|
||||
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
|
||||
private:
|
||||
|
||||
FermionOperator<Impl> & FermOp;// the basic operator
|
||||
|
||||
OperatorFunction<FermionField> &DerivativeSolver;
|
||||
OperatorFunction<FermionField> &ActionSolver;
|
||||
|
||||
FermionField PhiOdd; // the pseudo fermion field for this trajectory
|
||||
FermionField PhiEven; // the pseudo fermion field for this trajectory
|
||||
|
||||
public:
|
||||
/////////////////////////////////////////////////
|
||||
// Pass in required objects.
|
||||
/////////////////////////////////////////////////
|
||||
TwoFlavourEvenOddPseudoFermionAction(FermionOperator<Impl> &Op,
|
||||
OperatorFunction<FermionField> & DS,
|
||||
OperatorFunction<FermionField> & AS
|
||||
) :
|
||||
FermOp(Op),
|
||||
DerivativeSolver(DS),
|
||||
ActionSolver(AS),
|
||||
public:
|
||||
/////////////////////////////////////////////////
|
||||
// Pass in required objects.
|
||||
/////////////////////////////////////////////////
|
||||
TwoFlavourEvenOddPseudoFermionAction(FermionOperator<Impl> &Op,
|
||||
OperatorFunction<FermionField> &DS,
|
||||
OperatorFunction<FermionField> &AS)
|
||||
: FermOp(Op),
|
||||
DerivativeSolver(DS),
|
||||
ActionSolver(AS),
|
||||
PhiEven(Op.FermionRedBlackGrid()),
|
||||
PhiOdd(Op.FermionRedBlackGrid())
|
||||
{};
|
||||
|
||||
@ -131,9 +131,11 @@ namespace Grid{
|
||||
Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi
|
||||
X=zero;
|
||||
ActionSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi
|
||||
Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi
|
||||
//Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi
|
||||
// Multiply by Ydag
|
||||
RealD action = real(innerProduct(Y,X));
|
||||
|
||||
RealD action = norm2(Y);
|
||||
//RealD action = norm2(Y);
|
||||
|
||||
// The EE factorised block; normally can replace with zero if det is constant (gauge field indept)
|
||||
// Only really clover term that creates this. Leave the EE portion as a future to do to make most
|
||||
|
||||
Reference in New Issue
Block a user