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

Merge branch 'develop' into bugfix/LatTransfer

Browse Source 
* develop:
  Pass serial RNG around
  Sycl happier
This commit is contained in:
Michael Marshall 2021-03-04 20:41:30 +00:00
commit f252d69eef
31 changed files with 105 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Grid/Makefile.am
View File

@ -54,9 +54,11 @@ Version.h: version-cache
include Make.inc include Make.inc
include Eigen.inc include Eigen.inc
extra_sources+=$(ZWILS_FERMION_FILES)
extra_sources+=$(WILS_FERMION_FILES) extra_sources+=$(WILS_FERMION_FILES)
extra_sources+=$(STAG_FERMION_FILES) extra_sources+=$(STAG_FERMION_FILES)
if BUILD_ZMOBIUS
extra_sources+=$(ZWILS_FERMION_FILES)
endif
if BUILD_GPARITY if BUILD_GPARITY
extra_sources+=$(GP_FERMION_FILES) extra_sources+=$(GP_FERMION_FILES)
endif endif

2
Grid/cartesian/Cartesian_red_black.h
View File

@ -36,7 +36,7 @@ static const int CbBlack=1;
static const int Even =CbRed; static const int Even =CbRed;
static const int Odd =CbBlack; static const int Odd =CbBlack;
accelerator_inline int RedBlackCheckerBoardFromOindex (int oindex, Coordinate &rdim, Coordinate &chk_dim_msk) accelerator_inline int RedBlackCheckerBoardFromOindex (int oindex,const Coordinate &rdim,const Coordinate &chk_dim_msk)
{ {
int nd=rdim.size(); int nd=rdim.size();
Coordinate coor(nd); Coordinate coor(nd);

2
Grid/qcd/action/ActionBase.h
View File

@ -41,7 +41,7 @@ class Action
public: public:
bool is_smeared = false; bool is_smeared = false;
// Heatbath? // Heatbath?
virtual void refresh(const GaugeField& U, GridParallelRNG& pRNG) = 0; // refresh pseudofermions virtual void refresh(const GaugeField& U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) = 0; // refresh pseudofermions
virtual RealD S(const GaugeField& U) = 0; // evaluate the action virtual RealD S(const GaugeField& U) = 0; // evaluate the action
virtual void deriv(const GaugeField& U, GaugeField& dSdU) = 0; // evaluate the action derivative virtual void deriv(const GaugeField& U, GaugeField& dSdU) = 0; // evaluate the action derivative
virtual std::string action_name() = 0; // return the action name virtual std::string action_name() = 0; // return the action name

4
Grid/qcd/action/fermion/FermionOperatorImpl.h
View File

@ -153,8 +153,8 @@ public:
typedef typename Impl::StencilImpl StencilImpl; \ typedef typename Impl::StencilImpl StencilImpl; \
typedef typename Impl::ImplParams ImplParams; \ typedef typename Impl::ImplParams ImplParams; \
typedef typename Impl::StencilImpl::View_type StencilView; \ typedef typename Impl::StencilImpl::View_type StencilView; \
typedef typename ViewMap<FermionField>::Type FermionFieldView; \ typedef const typename ViewMap<FermionField>::Type FermionFieldView; \
typedef typename ViewMap<DoubledGaugeField>::Type DoubledGaugeFieldView; typedef const typename ViewMap<DoubledGaugeField>::Type DoubledGaugeFieldView;
#define INHERIT_IMPL_TYPES(Base) \ #define INHERIT_IMPL_TYPES(Base) \
INHERIT_GIMPL_TYPES(Base) \ INHERIT_GIMPL_TYPES(Base) \

24
Grid/qcd/action/fermion/WilsonCompressor.h
View File

@ -61,7 +61,7 @@ public:
typedef typename SiteHalfSpinor::vector_type vComplexHigh; typedef typename SiteHalfSpinor::vector_type vComplexHigh;
constexpr static int Nw=sizeof(SiteHalfSpinor)/sizeof(vComplexHigh); constexpr static int Nw=sizeof(SiteHalfSpinor)/sizeof(vComplexHigh);
accelerator_inline int CommDatumSize(void) { accelerator_inline int CommDatumSize(void) const {
return sizeof(SiteHalfCommSpinor); return sizeof(SiteHalfCommSpinor);
} }
@ -69,7 +69,7 @@ public:
/* Compress includes precision change if mpi data is not same */ /* Compress includes precision change if mpi data is not same */
/*****************************************************/ /*****************************************************/
template<class _SiteHalfSpinor, class _SiteSpinor> template<class _SiteHalfSpinor, class _SiteSpinor>
accelerator_inline void Compress(_SiteHalfSpinor *buf,Integer o,const _SiteSpinor &in) { accelerator_inline void Compress(_SiteHalfSpinor *buf,Integer o,const _SiteSpinor &in) const {
_SiteHalfSpinor tmp; _SiteHalfSpinor tmp;
projector::Proj(tmp,in,mu,dag); projector::Proj(tmp,in,mu,dag);
vstream(buf[o],tmp); vstream(buf[o],tmp);
@ -81,7 +81,7 @@ public:
accelerator_inline void Exchange(SiteHalfSpinor *mp, accelerator_inline void Exchange(SiteHalfSpinor *mp,
const SiteHalfSpinor * __restrict__ vp0, const SiteHalfSpinor * __restrict__ vp0,
const SiteHalfSpinor * __restrict__ vp1, const SiteHalfSpinor * __restrict__ vp1,
Integer type,Integer o){ Integer type,Integer o) const {
SiteHalfSpinor tmp1; SiteHalfSpinor tmp1;
SiteHalfSpinor tmp2; SiteHalfSpinor tmp2;
exchange(tmp1,tmp2,vp0[o],vp1[o],type); exchange(tmp1,tmp2,vp0[o],vp1[o],type);
@ -93,7 +93,7 @@ public:
/* Have a decompression step if mpi data is not same */ /* Have a decompression step if mpi data is not same */
/*****************************************************/ /*****************************************************/
accelerator_inline void Decompress(SiteHalfSpinor * __restrict__ out, accelerator_inline void Decompress(SiteHalfSpinor * __restrict__ out,
SiteHalfSpinor * __restrict__ in, Integer o) { SiteHalfSpinor * __restrict__ in, Integer o) const {
assert(0); assert(0);
} }
@ -103,7 +103,7 @@ public:
accelerator_inline void CompressExchange(SiteHalfSpinor * __restrict__ out0, accelerator_inline void CompressExchange(SiteHalfSpinor * __restrict__ out0,
SiteHalfSpinor * __restrict__ out1, SiteHalfSpinor * __restrict__ out1,
const SiteSpinor * __restrict__ in, const SiteSpinor * __restrict__ in,
Integer j,Integer k, Integer m,Integer type) Integer j,Integer k, Integer m,Integer type) const
{ {
SiteHalfSpinor temp1, temp2; SiteHalfSpinor temp1, temp2;
SiteHalfSpinor temp3, temp4; SiteHalfSpinor temp3, temp4;
@ -117,7 +117,7 @@ public:
/*****************************************************/ /*****************************************************/
/* Pass the info to the stencil */ /* Pass the info to the stencil */
/*****************************************************/ /*****************************************************/
accelerator_inline bool DecompressionStep(void) { return false; } accelerator_inline bool DecompressionStep(void) const { return false; }
}; };
@ -142,7 +142,7 @@ public:
typedef typename SiteHalfSpinor::vector_type vComplexHigh; typedef typename SiteHalfSpinor::vector_type vComplexHigh;
constexpr static int Nw=sizeof(SiteHalfSpinor)/sizeof(vComplexHigh); constexpr static int Nw=sizeof(SiteHalfSpinor)/sizeof(vComplexHigh);
accelerator_inline int CommDatumSize(void) { accelerator_inline int CommDatumSize(void) const {
return sizeof(SiteHalfCommSpinor); return sizeof(SiteHalfCommSpinor);
} }
@ -150,7 +150,7 @@ public:
/* Compress includes precision change if mpi data is not same */ /* Compress includes precision change if mpi data is not same */
/*****************************************************/ /*****************************************************/
template<class _SiteHalfSpinor, class _SiteSpinor> template<class _SiteHalfSpinor, class _SiteSpinor>
accelerator_inline void Compress(_SiteHalfSpinor *buf,Integer o,const _SiteSpinor &in) { accelerator_inline void Compress(_SiteHalfSpinor *buf,Integer o,const _SiteSpinor &in) const {
_SiteHalfSpinor hsp; _SiteHalfSpinor hsp;
SiteHalfCommSpinor *hbuf = (SiteHalfCommSpinor *)buf; SiteHalfCommSpinor *hbuf = (SiteHalfCommSpinor *)buf;
projector::Proj(hsp,in,mu,dag); projector::Proj(hsp,in,mu,dag);
@ -163,7 +163,7 @@ public:
accelerator_inline void Exchange(SiteHalfSpinor *mp, accelerator_inline void Exchange(SiteHalfSpinor *mp,
SiteHalfSpinor *vp0, SiteHalfSpinor *vp0,
SiteHalfSpinor *vp1, SiteHalfSpinor *vp1,
Integer type,Integer o){ Integer type,Integer o) const {
SiteHalfSpinor vt0,vt1; SiteHalfSpinor vt0,vt1;
SiteHalfCommSpinor *vpp0 = (SiteHalfCommSpinor *)vp0; SiteHalfCommSpinor *vpp0 = (SiteHalfCommSpinor *)vp0;
SiteHalfCommSpinor *vpp1 = (SiteHalfCommSpinor *)vp1; SiteHalfCommSpinor *vpp1 = (SiteHalfCommSpinor *)vp1;
@ -175,7 +175,7 @@ public:
/*****************************************************/ /*****************************************************/
/* Have a decompression step if mpi data is not same */ /* Have a decompression step if mpi data is not same */
/*****************************************************/ /*****************************************************/
accelerator_inline void Decompress(SiteHalfSpinor *out, SiteHalfSpinor *in, Integer o){ accelerator_inline void Decompress(SiteHalfSpinor *out, SiteHalfSpinor *in, Integer o) const {
SiteHalfCommSpinor *hin=(SiteHalfCommSpinor *)in; SiteHalfCommSpinor *hin=(SiteHalfCommSpinor *)in;
precisionChange((vComplexHigh *)&out[o],(vComplexLow *)&hin[o],Nw); precisionChange((vComplexHigh *)&out[o],(vComplexLow *)&hin[o],Nw);
} }
@ -186,7 +186,7 @@ public:
accelerator_inline void CompressExchange(SiteHalfSpinor *out0, accelerator_inline void CompressExchange(SiteHalfSpinor *out0,
SiteHalfSpinor *out1, SiteHalfSpinor *out1,
const SiteSpinor *in, const SiteSpinor *in,
Integer j,Integer k, Integer m,Integer type){ Integer j,Integer k, Integer m,Integer type) const {
SiteHalfSpinor temp1, temp2,temp3,temp4; SiteHalfSpinor temp1, temp2,temp3,temp4;
SiteHalfCommSpinor *hout0 = (SiteHalfCommSpinor *)out0; SiteHalfCommSpinor *hout0 = (SiteHalfCommSpinor *)out0;
SiteHalfCommSpinor *hout1 = (SiteHalfCommSpinor *)out1; SiteHalfCommSpinor *hout1 = (SiteHalfCommSpinor *)out1;
@ -200,7 +200,7 @@ public:
/*****************************************************/ /*****************************************************/
/* Pass the info to the stencil */ /* Pass the info to the stencil */
/*****************************************************/ /*****************************************************/
accelerator_inline bool DecompressionStep(void) { return true; } accelerator_inline bool DecompressionStep(void) const { return true; }
}; };

2
Grid/qcd/action/fermion/WilsonImpl.h
View File

@ -72,7 +72,7 @@ public:
typedef WilsonCompressor<SiteHalfCommSpinor,SiteHalfSpinor, SiteSpinor> Compressor; typedef WilsonCompressor<SiteHalfCommSpinor,SiteHalfSpinor, SiteSpinor> Compressor;
typedef WilsonImplParams ImplParams; typedef WilsonImplParams ImplParams;
typedef WilsonStencil<SiteSpinor, SiteHalfSpinor,ImplParams> StencilImpl; typedef WilsonStencil<SiteSpinor, SiteHalfSpinor,ImplParams> StencilImpl;
typedef typename StencilImpl::View_type StencilView; typedef const typename StencilImpl::View_type StencilView;
ImplParams Params; ImplParams Params;

28
Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h
View File

@ -93,7 +93,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
Chimu_32=coalescedReadPermute<ptype>(ref()(3)(2),perm); } Chimu_32=coalescedReadPermute<ptype>(ref()(3)(2),perm); }
#define PERMUTE_DIR(dir) ; #define PERMUTE_DIR(dir) ;
#else #else
#define LOAD_CHIMU \ #define LOAD_CHIMU(ptype) \
{const SiteSpinor & ref (in[offset]); \ {const SiteSpinor & ref (in[offset]); \
Chimu_00=ref()(0)(0);\ Chimu_00=ref()(0)(0);\
Chimu_01=ref()(0)(1);\ Chimu_01=ref()(0)(1);\
@ -482,19 +482,19 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
Simd U_11; \ Simd U_11; \
Simd U_21; Simd U_21;
#define ZERO_RESULT \ #define ZERO_RESULT \
result_00=S(0.0,0.0); \ zeroit(result_00); \
result_01=S(0.0,0.0); \ zeroit(result_01); \
result_02=S(0.0,0.0); \ zeroit(result_02); \
result_10=S(0.0,0.0); \ zeroit(result_10); \
result_11=S(0.0,0.0); \ zeroit(result_11); \
result_12=S(0.0,0.0); \ zeroit(result_12); \
result_20=S(0.0,0.0); \ zeroit(result_20); \
result_21=S(0.0,0.0); \ zeroit(result_21); \
result_22=S(0.0,0.0); \ zeroit(result_22); \
result_30=S(0.0,0.0); \ zeroit(result_30); \
result_31=S(0.0,0.0); \ zeroit(result_31); \
result_32=S(0.0,0.0); zeroit(result_32);
#define Chimu_00 Chi_00 #define Chimu_00 Chi_00
#define Chimu_01 Chi_01 #define Chimu_01 Chi_01

2
Grid/qcd/action/gauge/GaugeImplTypes.h
View File

@ -96,7 +96,7 @@ public:
/////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////
// Move these to another class // Move these to another class
// HMC auxiliary functions // HMC auxiliary functions
static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) static inline void generate_momenta(Field &P, GridSerialRNG & sRNG, GridParallelRNG &pRNG)
{ {
// Zbigniew Srocinsky thesis: // Zbigniew Srocinsky thesis:
// //

2
Grid/qcd/action/gauge/PlaqPlusRectangleAction.h
View File

@ -49,7 +49,7 @@ public:
virtual std::string action_name(){return "PlaqPlusRectangleAction";} virtual std::string action_name(){return "PlaqPlusRectangleAction";}
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {}; // noop as no pseudoferms virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {}; // noop as no pseudoferms
virtual std::string LogParameters(){ virtual std::string LogParameters(){
std::stringstream sstream; std::stringstream sstream;

3
Grid/qcd/action/gauge/WilsonGaugeAction.h
View File

@ -54,8 +54,7 @@ public:
return sstream.str(); return sstream.str();
} }
virtual void refresh(const GaugeField &U, virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG){}; // noop as no pseudoferms
GridParallelRNG &pRNG){}; // noop as no pseudoferms
virtual RealD S(const GaugeField &U) { virtual RealD S(const GaugeField &U) {
RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U); RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);

2
Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h
View File

@ -124,7 +124,7 @@ NAMESPACE_BEGIN(Grid);
// //
// As a check of rational require \Phi^dag M_{EOFA} \Phi == eta^dag M^-1/2^dag M M^-1/2 eta = eta^dag eta // As a check of rational require \Phi^dag M_{EOFA} \Phi == eta^dag M^-1/2^dag M M^-1/2 eta = eta^dag eta
// //
virtual void refresh(const GaugeField& U, GridParallelRNG& pRNG) virtual void refresh(const GaugeField& U, GridSerialRNG &sRNG, GridParallelRNG& pRNG)
{ {
Lop.ImportGauge(U); Lop.ImportGauge(U);
Rop.ImportGauge(U); Rop.ImportGauge(U);

11
Grid/qcd/action/pseudofermion/OneFlavourEvenOddRational.h
View File

@ -1,4 +1,3 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
@ -43,8 +42,7 @@ NAMESPACE_BEGIN(Grid);
// //
template <class Impl> template <class Impl>
class OneFlavourEvenOddRationalPseudoFermionAction class OneFlavourEvenOddRationalPseudoFermionAction : public Action<typename Impl::GaugeField> {
: public Action<typename Impl::GaugeField> {
public: public:
INHERIT_IMPL_TYPES(Impl); INHERIT_IMPL_TYPES(Impl);
@ -103,7 +101,7 @@ public:
return sstream.str(); return sstream.str();
} }
virtual void refresh(const GaugeField &U, GridParallelRNG &pRNG) { virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG) {
// P(phi) = e^{- phi^dag (MpcdagMpc)^-1/2 phi} // P(phi) = e^{- phi^dag (MpcdagMpc)^-1/2 phi}
// = e^{- phi^dag (MpcdagMpc)^-1/4 (MpcdagMpc)^-1/4 phi} // = e^{- phi^dag (MpcdagMpc)^-1/4 (MpcdagMpc)^-1/4 phi}
// Phi = MpcdagMpc^{1/4} eta // Phi = MpcdagMpc^{1/4} eta
@ -156,7 +154,10 @@ public:
msCG(Mpc, PhiOdd, Y); msCG(Mpc, PhiOdd, Y);
if ( (rand()%param.BoundsCheckFreq)==0 ) { auto grid = FermOp.FermionGrid();
auto r=rand();
grid->Broadcast(0,r);
if ( (r%param.BoundsCheckFreq)==0 ) {
FermionField gauss(FermOp.FermionRedBlackGrid()); FermionField gauss(FermOp.FermionRedBlackGrid());
gauss = PhiOdd; gauss = PhiOdd;
HighBoundCheck(Mpc,gauss,param.hi); HighBoundCheck(Mpc,gauss,param.hi);

7
Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h
View File

@ -101,7 +101,7 @@ NAMESPACE_BEGIN(Grid);
} }
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) { virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
// //
@ -170,7 +170,10 @@ NAMESPACE_BEGIN(Grid);
msCG_M(MdagM,X,Y); msCG_M(MdagM,X,Y);
// Randomly apply rational bounds checks. // Randomly apply rational bounds checks.
if ( (rand()%param.BoundsCheckFreq)==0 ) { auto grid = NumOp.FermionGrid();
auto r=rand();
grid->Broadcast(0,r);
if ( (r%param.BoundsCheckFreq)==0 ) {
FermionField gauss(NumOp.FermionRedBlackGrid()); FermionField gauss(NumOp.FermionRedBlackGrid());
gauss = PhiOdd; gauss = PhiOdd;
HighBoundCheck(MdagM,gauss,param.hi); HighBoundCheck(MdagM,gauss,param.hi);

7
Grid/qcd/action/pseudofermion/OneFlavourRational.h
View File

@ -98,7 +98,7 @@ NAMESPACE_BEGIN(Grid);
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) { virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag (MdagM)^-1/2 phi} // P(phi) = e^{- phi^dag (MdagM)^-1/2 phi}
@ -142,7 +142,10 @@ NAMESPACE_BEGIN(Grid);
msCG(MdagMOp,Phi,Y); msCG(MdagMOp,Phi,Y);
if ( (rand()%param.BoundsCheckFreq)==0 ) { auto grid = FermOp.FermionGrid();
auto r=rand();
grid->Broadcast(0,r);
if ( (r%param.BoundsCheckFreq)==0 ) {
FermionField gauss(FermOp.FermionGrid()); FermionField gauss(FermOp.FermionGrid());
gauss = Phi; gauss = Phi;
HighBoundCheck(MdagMOp,gauss,param.hi); HighBoundCheck(MdagMOp,gauss,param.hi);

7
Grid/qcd/action/pseudofermion/OneFlavourRationalRatio.h
View File

@ -95,7 +95,7 @@ NAMESPACE_BEGIN(Grid);
} }
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) { virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
// S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
// //
@ -156,7 +156,10 @@ NAMESPACE_BEGIN(Grid);
msCG_M(MdagM,X,Y); msCG_M(MdagM,X,Y);
// Randomly apply rational bounds checks. // Randomly apply rational bounds checks.
if ( (rand()%param.BoundsCheckFreq)==0 ) { auto grid = NumOp.FermionGrid();
auto r=rand();
grid->Broadcast(0,r);
if ( (r%param.BoundsCheckFreq)==0 ) {
FermionField gauss(NumOp.FermionGrid()); FermionField gauss(NumOp.FermionGrid());
gauss = Phi; gauss = Phi;
HighBoundCheck(MdagM,gauss,param.hi); HighBoundCheck(MdagM,gauss,param.hi);

2
Grid/qcd/action/pseudofermion/TwoFlavour.h
View File

@ -73,7 +73,7 @@ public:
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
// Push the gauge field in to the dops. Assume any BC's and smearing already applied // Push the gauge field in to the dops. Assume any BC's and smearing already applied
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
virtual void refresh(const GaugeField &U, GridParallelRNG &pRNG) { virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG) {
// P(phi) = e^{- phi^dag (MdagM)^-1 phi} // P(phi) = e^{- phi^dag (MdagM)^-1 phi}
// Phi = Mdag eta // Phi = Mdag eta
// P(eta) = e^{- eta^dag eta} // P(eta) = e^{- eta^dag eta}

2
Grid/qcd/action/pseudofermion/TwoFlavourEvenOdd.h
View File

@ -77,7 +77,7 @@ public:
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
// Push the gauge field in to the dops. Assume any BC's and smearing already applied // Push the gauge field in to the dops. Assume any BC's and smearing already applied
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) { virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag (MpcdagMpc)^-1 phi} // P(phi) = e^{- phi^dag (MpcdagMpc)^-1 phi}
// Phi = McpDag eta // Phi = McpDag eta

2
Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h
View File

@ -84,7 +84,7 @@ NAMESPACE_BEGIN(Grid);
} }
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) { virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag Vpc (MpcdagMpc)^-1 Vpcdag phi} // P(phi) = e^{- phi^dag Vpc (MpcdagMpc)^-1 Vpcdag phi}
// //

2
Grid/qcd/action/pseudofermion/TwoFlavourRatio.h
View File

@ -64,7 +64,7 @@ public:
return sstream.str(); return sstream.str();
} }
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) { virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi} // P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi}
// //

2
Grid/qcd/action/scalar/ScalarAction.h
View File

@ -55,7 +55,7 @@ public:
} }
virtual std::string action_name() {return "ScalarAction";} virtual std::string action_name() {return "ScalarAction";}
virtual void refresh(const Field &U, GridParallelRNG &pRNG) {} // noop as no pseudoferms virtual void refresh(const Field &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG) {} // noop as no pseudoferms
virtual RealD S(const Field &p) { virtual RealD S(const Field &p) {
return (mass_square * 0.5 + Nd) * ScalarObs<Impl>::sumphisquared(p) + return (mass_square * 0.5 + Nd) * ScalarObs<Impl>::sumphisquared(p) +

4
Grid/qcd/action/scalar/ScalarImpl.h
View File

@ -27,7 +27,7 @@ public:
typedef Field FermionField; typedef Field FermionField;
typedef Field PropagatorField; typedef Field PropagatorField;
static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){ static inline void generate_momenta(Field& P, GridSerialRNG &sRNG, GridParallelRNG& pRNG){
RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR); // CPS/UKQCD momentum rescaling RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR); // CPS/UKQCD momentum rescaling
gaussian(pRNG, P); gaussian(pRNG, P);
P *= scale; P *= scale;
@ -151,7 +151,7 @@ public:
out = one / out; out = one / out;
} }
static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) static inline void generate_momenta(Field &P, GridSerialRNG & sRNG, GridParallelRNG &pRNG)
{ {
RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR); // CPS/UKQCD momentum rescaling RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR); // CPS/UKQCD momentum rescaling
#ifndef USE_FFT_ACCELERATION #ifndef USE_FFT_ACCELERATION

2
Grid/qcd/action/scalar/ScalarInteractionAction.h
View File

@ -77,7 +77,7 @@ public:
virtual std::string action_name() { return "ScalarAction"; } virtual std::string action_name() { return "ScalarAction"; }
virtual void refresh(const Field &U, GridParallelRNG &pRNG) {} virtual void refresh(const Field &U, GridSerialRNG & sRNG, GridParallelRNG &pRNG) {}
virtual RealD S(const Field &p) virtual RealD S(const Field &p)
{ {

2
Grid/qcd/hmc/HMC.h
View File

@ -139,7 +139,7 @@ private:
// Evolution // Evolution
///////////////////////////////////////////////////////// /////////////////////////////////////////////////////////
RealD evolve_hmc_step(Field &U) { RealD evolve_hmc_step(Field &U) {
TheIntegrator.refresh(U, pRNG); // set U and initialize P and phi's TheIntegrator.refresh(U, sRNG, pRNG); // set U and initialize P and phi's
RealD H0 = TheIntegrator.S(U); // initial state action RealD H0 = TheIntegrator.S(U); // initial state action

13
Grid/qcd/hmc/integrators/Integrator.h
View File

@ -236,10 +236,9 @@ public:
// over the representations // over the representations
struct _refresh { struct _refresh {
template <class FieldType, class Repr> template <class FieldType, class Repr>
void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep, void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep, GridSerialRNG & sRNG, GridParallelRNG& pRNG) {
GridParallelRNG& pRNG) {
for (int a = 0; a < repr_set.size(); ++a){ for (int a = 0; a < repr_set.size(); ++a){
repr_set.at(a)->refresh(Rep.U, pRNG); repr_set.at(a)->refresh(Rep.U, sRNG, pRNG);
std::cout << GridLogDebug << "Hirep refreshing pseudofermions" << std::endl; std::cout << GridLogDebug << "Hirep refreshing pseudofermions" << std::endl;
} }
@ -247,12 +246,12 @@ public:
} refresh_hireps{}; } refresh_hireps{};
// Initialization of momenta and actions // Initialization of momenta and actions
void refresh(Field& U, GridParallelRNG& pRNG) void refresh(Field& U, GridSerialRNG & sRNG, GridParallelRNG& pRNG)
{ {
assert(P.Grid() == U.Grid()); assert(P.Grid() == U.Grid());
std::cout << GridLogIntegrator << "Integrator refresh\n"; std::cout << GridLogIntegrator << "Integrator refresh\n";
FieldImplementation::generate_momenta(P, pRNG); FieldImplementation::generate_momenta(P, sRNG, pRNG);
// Update the smeared fields, can be implemented as observer // Update the smeared fields, can be implemented as observer
// necessary to keep the fields updated even after a reject // necessary to keep the fields updated even after a reject
@ -269,11 +268,11 @@ public:
// get gauge field from the SmearingPolicy and // get gauge field from the SmearingPolicy and
// based on the boolean is_smeared in actionID // based on the boolean is_smeared in actionID
Field& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared); Field& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
as[level].actions.at(actionID)->refresh(Us, pRNG); as[level].actions.at(actionID)->refresh(Us, sRNG, pRNG);
} }
// Refresh the higher representation actions // Refresh the higher representation actions
as[level].apply(refresh_hireps, Representations, pRNG); as[level].apply(refresh_hireps, Representations, sRNG, pRNG);
} }
MomFilter->applyFilter(P); MomFilter->applyFilter(P);

1
Grid/simd/Grid_gpu_vec.h
View File

@ -67,6 +67,7 @@ public:
accelerator_inline GpuComplex(const GpuComplex &zz) { z = zz.z;}; accelerator_inline GpuComplex(const GpuComplex &zz) { z = zz.z;};
accelerator_inline Real real(void) const { return z.x; }; accelerator_inline Real real(void) const { return z.x; };
accelerator_inline Real imag(void) const { return z.y; }; accelerator_inline Real imag(void) const { return z.y; };
accelerator_inline GpuComplex &operator=(const Zero &zz) { z.x = 0; z.y=0; return *this; };
accelerator_inline GpuComplex &operator*=(const GpuComplex &r) { accelerator_inline GpuComplex &operator*=(const GpuComplex &r) {
*this = (*this) * r; *this = (*this) * r;
return *this; return *this;

14
Grid/stencil/SimpleCompressor.h
View File

@ -7,20 +7,20 @@ template<class vobj>
class SimpleCompressor { class SimpleCompressor {
public: public:
void Point(int) {}; void Point(int) {};
accelerator_inline int CommDatumSize(void) { return sizeof(vobj); } accelerator_inline int CommDatumSize(void) const { return sizeof(vobj); }
accelerator_inline bool DecompressionStep(void) { return false; } accelerator_inline bool DecompressionStep(void) const { return false; }
template<class cobj> accelerator_inline void Compress(cobj *buf,int o,const cobj &in) { buf[o]=in; } template<class cobj> accelerator_inline void Compress(cobj *buf,int o,const cobj &in) const { buf[o]=in; }
accelerator_inline void Exchange(vobj *mp,vobj *vp0,vobj *vp1,Integer type,Integer o){ accelerator_inline void Exchange(vobj *mp,vobj *vp0,vobj *vp1,Integer type,Integer o) const {
exchange(mp[2*o],mp[2*o+1],vp0[o],vp1[o],type); exchange(mp[2*o],mp[2*o+1],vp0[o],vp1[o],type);
} }
accelerator_inline void Decompress(vobj *out,vobj *in, int o){ assert(0); } accelerator_inline void Decompress(vobj *out,vobj *in, int o) const { assert(0); }
accelerator_inline void CompressExchange(vobj *out0,vobj *out1,const vobj *in, accelerator_inline void CompressExchange(vobj *out0,vobj *out1,const vobj *in,
int j,int k, int m,int type){ int j,int k, int m,int type) const {
exchange(out0[j],out1[j],in[k],in[m],type); exchange(out0[j],out1[j],in[k],in[m],type);
} }
// For cshift. Cshift should drop compressor coupling altogether // For cshift. Cshift should drop compressor coupling altogether
// because I had to decouple the code from the Stencil anyway // because I had to decouple the code from the Stencil anyway
accelerator_inline vobj operator() (const vobj &arg) { accelerator_inline vobj operator() (const vobj &arg) const {
return arg; return arg;
} }
}; };

14
Grid/stencil/Stencil.h
View File

@ -147,16 +147,16 @@ class CartesianStencilAccelerator {
cobj* u_recv_buf_p; cobj* u_recv_buf_p;
cobj* u_send_buf_p; cobj* u_send_buf_p;
accelerator_inline cobj *CommBuf(void) { return u_recv_buf_p; } accelerator_inline cobj *CommBuf(void) const { return u_recv_buf_p; }
accelerator_inline int GetNodeLocal(int osite,int point) { accelerator_inline int GetNodeLocal(int osite,int point) const {
return this->_entries_p[point+this->_npoints*osite]._is_local; return this->_entries_p[point+this->_npoints*osite]._is_local;
} }
accelerator_inline StencilEntry * GetEntry(int &ptype,int point,int osite) { accelerator_inline StencilEntry * GetEntry(int &ptype,int point,int osite) const {
ptype = this->_permute_type[point]; return & this->_entries_p[point+this->_npoints*osite]; ptype = this->_permute_type[point]; return & this->_entries_p[point+this->_npoints*osite];
} }
accelerator_inline uint64_t GetInfo(int &ptype,int &local,int &perm,int point,int ent,uint64_t base) { accelerator_inline uint64_t GetInfo(int &ptype,int &local,int &perm,int point,int ent,uint64_t base) const {
uint64_t cbase = (uint64_t)&u_recv_buf_p[0]; uint64_t cbase = (uint64_t)&u_recv_buf_p[0];
local = this->_entries_p[ent]._is_local; local = this->_entries_p[ent]._is_local;
perm = this->_entries_p[ent]._permute; perm = this->_entries_p[ent]._permute;
@ -168,14 +168,14 @@ class CartesianStencilAccelerator {
} }
} }
accelerator_inline uint64_t GetPFInfo(int ent,uint64_t base) { accelerator_inline uint64_t GetPFInfo(int ent,uint64_t base) const {
uint64_t cbase = (uint64_t)&u_recv_buf_p[0]; uint64_t cbase = (uint64_t)&u_recv_buf_p[0];
int local = this->_entries_p[ent]._is_local; int local = this->_entries_p[ent]._is_local;
if (local) return base + this->_entries_p[ent]._byte_offset; if (local) return base + this->_entries_p[ent]._byte_offset;
else return cbase + this->_entries_p[ent]._byte_offset; else return cbase + this->_entries_p[ent]._byte_offset;
} }
accelerator_inline void iCoorFromIindex(Coordinate &coor,int lane) accelerator_inline void iCoorFromIindex(Coordinate &coor,int lane) const
{ {
Lexicographic::CoorFromIndex(coor,lane,this->_simd_layout); Lexicographic::CoorFromIndex(coor,lane,this->_simd_layout);
} }
@ -221,7 +221,7 @@ public:
typedef typename cobj::vector_type vector_type; typedef typename cobj::vector_type vector_type;
typedef typename cobj::scalar_type scalar_type; typedef typename cobj::scalar_type scalar_type;
typedef typename cobj::scalar_object scalar_object; typedef typename cobj::scalar_object scalar_object;
typedef CartesianStencilView<vobj,cobj,Parameters> View_type; typedef const CartesianStencilView<vobj,cobj,Parameters> View_type;
typedef typename View_type::StencilVector StencilVector; typedef typename View_type::StencilVector StencilVector;
/////////////////////////////////////////// ///////////////////////////////////////////
// Helper structs // Helper structs

4
Grid/tensors/Tensor_SIMT.h
View File

@ -66,7 +66,7 @@ void coalescedWriteNonTemporal(vobj & __restrict__ vec,const vobj & __restrict__
#ifndef GRID_SYCL #ifndef GRID_SYCL
// Use the scalar as our own complex on GPU // Use the scalar as our own complex on GPU ... thrust::complex or std::complex
template<class vsimd,IfSimd<vsimd> = 0> accelerator_inline template<class vsimd,IfSimd<vsimd> = 0> accelerator_inline
typename vsimd::scalar_type typename vsimd::scalar_type
coalescedRead(const vsimd & __restrict__ vec,int lane=acceleratorSIMTlane(vsimd::Nsimd())) coalescedRead(const vsimd & __restrict__ vec,int lane=acceleratorSIMTlane(vsimd::Nsimd()))
@ -96,6 +96,8 @@ void coalescedWrite(vsimd & __restrict__ vec,
p[lane]=extracted; p[lane]=extracted;
} }
#else #else
// For SyCL have option to use GpuComplex from inside the vector type in SIMT loops
// Faster for some reason
template<class vsimd,IfSimd<vsimd> = 0> accelerator_inline template<class vsimd,IfSimd<vsimd> = 0> accelerator_inline
typename vsimd::vector_type::datum typename vsimd::vector_type::datum
coalescedRead(const vsimd & __restrict__ vec,int lane=acceleratorSIMTlane(vsimd::Nsimd())) coalescedRead(const vsimd & __restrict__ vec,int lane=acceleratorSIMTlane(vsimd::Nsimd()))

2
Grid/threads/Accelerator.h
View File

@ -456,7 +456,7 @@ accelerator_inline void acceleratorSynchronise(void)
__syncwarp(); __syncwarp();
#endif #endif
#ifdef GRID_SYCL #ifdef GRID_SYCL
// No barrier call on SYCL?? // Option get __spir:: stuff to do warp barrier cl::sycl::detail::workGroupBarrier();
#endif #endif
#ifdef GRID_HIP #ifdef GRID_HIP
__syncthreads(); __syncthreads();

2
benchmarks/Benchmark_dwf_fp32.cc
View File

@ -53,7 +53,7 @@ int main (int argc, char ** argv)
int threads = GridThread::GetThreads(); int threads = GridThread::GetThreads();
Coordinate latt4 = GridDefaultLatt(); Coordinate latt4 = GridDefaultLatt();
int Ls=8; int Ls=16;
for(int i=0;i<argc;i++) for(int i=0;i<argc;i++)
if(std::string(argv[i]) == "-Ls"){ if(std::string(argv[i]) == "-Ls"){
std::stringstream ss(argv[i+1]); ss >> Ls; std::stringstream ss(argv[i+1]); ss >> Ls;

11
configure.ac
View File

@ -140,12 +140,23 @@ AC_ARG_ENABLE([gparity],
[ac_GPARITY=${enable_gparity}], [ac_GPARITY=yes]) [ac_GPARITY=${enable_gparity}], [ac_GPARITY=yes])
AM_CONDITIONAL(BUILD_GPARITY, [ test "${ac_GPARITY}X" == "yesX" ]) AM_CONDITIONAL(BUILD_GPARITY, [ test "${ac_GPARITY}X" == "yesX" ])
AC_ARG_ENABLE([zmobius],
[AC_HELP_STRING([--enable-zmobius=yes|no], [enable Zmobius support])],
[ac_ZMOBIUS=${enable_zmobius}], [ac_ZMOBIUS=yes])
AM_CONDITIONAL(BUILD_ZMOBIUS, [ test "${ac_ZMOBIUS}X" == "yesX" ])
case ${ac_FERMION_REPS} in case ${ac_FERMION_REPS} in
yes) AC_DEFINE([ENABLE_FERMION_REPS],[1],[non QCD fermion reps]);; yes) AC_DEFINE([ENABLE_FERMION_REPS],[1],[non QCD fermion reps]);;
esac esac
case ${ac_GPARITY} in case ${ac_GPARITY} in
yes) AC_DEFINE([ENABLE_GPARITY],[1],[fermion actions with GPARITY BCs]);; yes) AC_DEFINE([ENABLE_GPARITY],[1],[fermion actions with GPARITY BCs]);;
esac esac
case ${ac_ZMOBIUS} in
yes) AC_DEFINE([ENABLE_ZMOBIUS],[1],[Zmobius fermion actions]);;
esac
############### Nc ############### Nc
AC_ARG_ENABLE([Nc], AC_ARG_ENABLE([Nc],
[AC_HELP_STRING([--enable-Nc=2|3|4], [enable number of colours])], [AC_HELP_STRING([--enable-Nc=2|3|4], [enable number of colours])],