diff --git a/Grid/Makefile.am b/Grid/Makefile.am index ded6d146..7c3c151b 100644 --- a/Grid/Makefile.am +++ b/Grid/Makefile.am @@ -54,9 +54,11 @@ Version.h: version-cache include Make.inc include Eigen.inc -extra_sources+=$(ZWILS_FERMION_FILES) extra_sources+=$(WILS_FERMION_FILES) extra_sources+=$(STAG_FERMION_FILES) +if BUILD_ZMOBIUS + extra_sources+=$(ZWILS_FERMION_FILES) +endif if BUILD_GPARITY extra_sources+=$(GP_FERMION_FILES) endif diff --git a/Grid/cartesian/Cartesian_red_black.h b/Grid/cartesian/Cartesian_red_black.h index b71981f5..092d4910 100644 --- a/Grid/cartesian/Cartesian_red_black.h +++ b/Grid/cartesian/Cartesian_red_black.h @@ -36,7 +36,7 @@ static const int CbBlack=1; static const int Even =CbRed; 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(); Coordinate coor(nd); diff --git a/Grid/communicator/Communicator_base.h b/Grid/communicator/Communicator_base.h index bb06d43f..a15f9789 100644 --- a/Grid/communicator/Communicator_base.h +++ b/Grid/communicator/Communicator_base.h @@ -1,4 +1,3 @@ - /************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -108,6 +107,8 @@ public: //////////////////////////////////////////////////////////// // Reduction //////////////////////////////////////////////////////////// + void GlobalMax(RealD &); + void GlobalMax(RealF &); void GlobalSum(RealF &); void GlobalSumVector(RealF *,int N); void GlobalSum(RealD &); diff --git a/Grid/communicator/Communicator_mpi3.cc b/Grid/communicator/Communicator_mpi3.cc index c6543851..5713fe35 100644 --- a/Grid/communicator/Communicator_mpi3.cc +++ b/Grid/communicator/Communicator_mpi3.cc @@ -275,6 +275,16 @@ void CartesianCommunicator::GlobalXOR(uint64_t &u){ int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT64_T,MPI_BXOR,communicator); assert(ierr==0); } +void CartesianCommunicator::GlobalMax(float &f) +{ + int ierr=MPI_Allreduce(MPI_IN_PLACE,&f,1,MPI_FLOAT,MPI_MAX,communicator); + assert(ierr==0); +} +void CartesianCommunicator::GlobalMax(double &d) +{ + int ierr = MPI_Allreduce(MPI_IN_PLACE,&d,1,MPI_DOUBLE,MPI_MAX,communicator); + assert(ierr==0); +} void CartesianCommunicator::GlobalSum(float &f){ int ierr=MPI_Allreduce(MPI_IN_PLACE,&f,1,MPI_FLOAT,MPI_SUM,communicator); assert(ierr==0); diff --git a/Grid/communicator/Communicator_none.cc b/Grid/communicator/Communicator_none.cc index 6cb431a2..beb2cc97 100644 --- a/Grid/communicator/Communicator_none.cc +++ b/Grid/communicator/Communicator_none.cc @@ -67,6 +67,8 @@ CartesianCommunicator::CartesianCommunicator(const Coordinate &processors) CartesianCommunicator::~CartesianCommunicator(){} +void CartesianCommunicator::GlobalMax(float &){} +void CartesianCommunicator::GlobalMax(double &){} void CartesianCommunicator::GlobalSum(float &){} void CartesianCommunicator::GlobalSumVector(float *,int N){} void CartesianCommunicator::GlobalSum(double &){} diff --git a/Grid/communicator/SharedMemoryMPI.cc b/Grid/communicator/SharedMemoryMPI.cc index a12418e6..466f6a1e 100644 --- a/Grid/communicator/SharedMemoryMPI.cc +++ b/Grid/communicator/SharedMemoryMPI.cc @@ -7,6 +7,7 @@ Copyright (C) 2015 Author: Peter Boyle +Author: Christoph Lehner 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 @@ -169,6 +170,23 @@ static inline int divides(int a,int b) } void GlobalSharedMemory::GetShmDims(const Coordinate &WorldDims,Coordinate &ShmDims) { + //////////////////////////////////////////////////////////////// + // Allow user to configure through environment variable + //////////////////////////////////////////////////////////////// + char* str = getenv(("GRID_SHM_DIMS_" + std::to_string(ShmDims.size())).c_str()); + if ( str ) { + std::vector IntShmDims; + GridCmdOptionIntVector(std::string(str),IntShmDims); + assert(IntShmDims.size() == WorldDims.size()); + long ShmSize = 1; + for (int dim=0;dim &rhs, int n1=rhs.Grid()->_slice_stride[dimension]; if ( cbmask ==0x3){ -#ifdef ACCELERATOR_CSHIFT +#ifdef ACCELERATOR_CSHIFT autoView(rhs_v , rhs, AcceleratorRead); - accelerator_for2d(n,e1,b,e2,1,{ + accelerator_for(nn,e1*e2,1,{ + int n = nn%e1; + int b = nn/e1; int o = n*n1; int offset = b+n*e2; @@ -135,7 +137,9 @@ Gather_plane_extract(const Lattice &rhs, std::cout << " Dense packed buffer WARNING " < void Scatter_plane_merge(Lattice &rhs,ExtractPointerA int _slice_block = rhs.Grid()->_slice_block[dimension]; #ifdef ACCELERATOR_CSHIFT autoView( rhs_v , rhs, AcceleratorWrite); - accelerator_for2d(n,e1,b,e2,1,{ + accelerator_for(nn,e1*e2,1,{ + int n = nn%e1; + int b = nn/e1; int o = n*_slice_stride; int offset = b+n*_slice_block; merge(rhs_v[so+o+b],pointers,offset); @@ -274,7 +280,7 @@ template void Scatter_plane_merge(Lattice &rhs,ExtractPointerA // Case of SIMD split AND checker dim cannot currently be hit, except in // Test_cshift_red_black code. - // std::cout << "Scatter_plane merge assert(0); think this is buggy FIXME "<< std::endl;// think this is buggy FIXME + std::cout << "Scatter_plane merge assert(0); think this is buggy FIXME "<< std::endl;// think this is buggy FIXME std::cout<<" Unthreaded warning -- buffer is not densely packed ??"< void Cshift_comms(Lattice &ret,const Lattice &r assert(shift_slice_nblock[dimension]*rhs.Grid()->_slice_block[dimension]; - cshiftVector send_buf(buffer_size); - cshiftVector recv_buf(buffer_size); + static cshiftVector send_buf; send_buf.resize(buffer_size); + static cshiftVector recv_buf; recv_buf.resize(buffer_size); int cb= (cbmask==0x2)? Odd : Even; int sshift= rhs.Grid()->CheckerBoardShiftForCB(rhs.Checkerboard(),dimension,shift,cb); @@ -198,8 +198,8 @@ template void Cshift_comms_simd(Lattice &ret,const Lattice_slice_nblock[dimension]*grid->_slice_block[dimension]; // int words = sizeof(vobj)/sizeof(vector_type); - std::vector > send_buf_extract(Nsimd); - std::vector > recv_buf_extract(Nsimd); + static std::vector > send_buf_extract; send_buf_extract.resize(Nsimd); + static std::vector > recv_buf_extract; recv_buf_extract.resize(Nsimd); scalar_object * recv_buf_extract_mpi; scalar_object * send_buf_extract_mpi; @@ -294,8 +294,8 @@ template void Cshift_comms(Lattice &ret,const Lattice &r assert(shift_slice_nblock[dimension]*rhs.Grid()->_slice_block[dimension]; - cshiftVector send_buf_v(buffer_size); - cshiftVector recv_buf_v(buffer_size); + static cshiftVector send_buf_v; send_buf_v.resize(buffer_size); + static cshiftVector recv_buf_v; recv_buf_v.resize(buffer_size); vobj *send_buf; vobj *recv_buf; { @@ -381,8 +381,8 @@ template void Cshift_comms_simd(Lattice &ret,const Lattice_slice_nblock[dimension]*grid->_slice_block[dimension]; // int words = sizeof(vobj)/sizeof(vector_type); - std::vector > send_buf_extract(Nsimd); - std::vector > recv_buf_extract(Nsimd); + static std::vector > send_buf_extract; send_buf_extract.resize(Nsimd); + static std::vector > recv_buf_extract; recv_buf_extract.resize(Nsimd); scalar_object * recv_buf_extract_mpi; scalar_object * send_buf_extract_mpi; { diff --git a/Grid/lattice/Lattice_reduction.h b/Grid/lattice/Lattice_reduction.h index c2955485..0a5fbcb6 100644 --- a/Grid/lattice/Lattice_reduction.h +++ b/Grid/lattice/Lattice_reduction.h @@ -96,8 +96,34 @@ inline typename vobj::scalar_objectD sumD_cpu(const vobj *arg, Integer osites) ssobj ret = ssum; return ret; } +/* +Threaded max, don't use for now +template +inline Double max(const Double *arg, Integer osites) +{ + // const int Nsimd = vobj::Nsimd(); + const int nthread = GridThread::GetThreads(); - + std::vector maxarray(nthread); + + thread_for(thr,nthread, { + int nwork, mywork, myoff; + nwork = osites; + GridThread::GetWork(nwork,thr,mywork,myoff); + Double max=arg[0]; + for(int ss=myoff;ss max ) max = arg[ss]; + } + maxarray[thr]=max; + }); + + Double tmax=maxarray[0]; + for(int i=0;itmax) tmax = maxarray[i]; + } + return tmax; +} +*/ template inline typename vobj::scalar_object sum(const vobj *arg, Integer osites) { @@ -141,6 +167,32 @@ template inline RealD norm2(const Lattice &arg){ return real(nrm); } +//The global maximum of the site norm2 +template inline RealD maxLocalNorm2(const Lattice &arg) +{ + typedef typename vobj::tensor_reduced vscalar; //iScalar > > + typedef typename vscalar::scalar_object scalar; //iScalar > > + + Lattice inner = localNorm2(arg); + + auto grid = arg.Grid(); + + RealD max; + for(int l=0;llSites();l++){ + Coordinate coor; + scalar val; + RealD r; + grid->LocalIndexToLocalCoor(l,coor); + peekLocalSite(val,inner,coor); + r=real(TensorRemove(val)); + if( (l==0) || (r>max)){ + max=r; + } + } + grid->GlobalMax(max); + return max; +} + // Double inner product template inline ComplexD rankInnerProduct(const Lattice &left,const Lattice &right) diff --git a/Grid/lattice/Lattice_transfer.h b/Grid/lattice/Lattice_transfer.h index 202103d1..e7ab77a8 100644 --- a/Grid/lattice/Lattice_transfer.h +++ b/Grid/lattice/Lattice_transfer.h @@ -97,6 +97,20 @@ accelerator_inline void convertType(ComplexF & out, const std::complex & out = in; } +template +accelerator_inline EnableIf> convertType(T & out, const T & in) { + out = in; +} + +// This would allow for conversions between GridFundamental types, but is not strictly needed as yet +/*template +accelerator_inline typename std::enable_if::value && isGridFundamental::value>::type +// Or to make this very broad, conversions between anything that's not a GridTensor could be allowed +//accelerator_inline typename std::enable_if::value && !isGridTensor::value>::type +convertType(T1 & out, const T2 & in) { + out = in; +}*/ + #ifdef GRID_SIMT accelerator_inline void convertType(vComplexF & out, const ComplexF & in) { ((ComplexF*)&out)[acceleratorSIMTlane(vComplexF::Nsimd())] = in; @@ -117,23 +131,18 @@ accelerator_inline void convertType(vComplexD2 & out, const vComplexF & in) { Optimization::PrecisionChange::StoD(in.v,out._internal[0].v,out._internal[1].v); } -template - accelerator_inline void convertType(iMatrix & out, const iMatrix & in); -template - accelerator_inline void convertType(iVector & out, const iVector & in); - -template::value, T1>::type* = nullptr> -accelerator_inline void convertType(T1 & out, const iScalar & in) { - convertType(out,in._internal); +template +accelerator_inline void convertType(iScalar & out, const iScalar & in) { + convertType(out._internal,in._internal); } -template::value, T1>::type* = nullptr> -accelerator_inline void convertType(T1 & out, const iScalar & in) { +template +accelerator_inline NotEnableIf> convertType(T1 & out, const iScalar & in) { convertType(out,in._internal); } template -accelerator_inline void convertType(iScalar & out, const T2 & in) { +accelerator_inline NotEnableIf> convertType(iScalar & out, const T2 & in) { convertType(out._internal,in); } @@ -150,11 +159,6 @@ accelerator_inline void convertType(iVector & out, const iVector & i convertType(out._internal[i],in._internal[i]); } -template::value, T>::type* = nullptr> -accelerator_inline void convertType(T & out, const T & in) { - out = in; -} - template accelerator_inline void convertType(Lattice & out, const Lattice & in) { autoView( out_v , out,AcceleratorWrite); diff --git a/Grid/lattice/Lattice_view.h b/Grid/lattice/Lattice_view.h index 3b76b921..cb568abd 100644 --- a/Grid/lattice/Lattice_view.h +++ b/Grid/lattice/Lattice_view.h @@ -67,9 +67,14 @@ public: accelerator_inline const vobj & operator()(size_t i) const { return this->_odata[i]; } #endif +#if 1 + // accelerator_inline const vobj & operator[](size_t i) const { return this->_odata[i]; }; + accelerator_inline vobj & operator[](size_t i) const { return this->_odata[i]; }; +#else accelerator_inline const vobj & operator[](size_t i) const { return this->_odata[i]; }; accelerator_inline vobj & operator[](size_t i) { return this->_odata[i]; }; - +#endif + accelerator_inline uint64_t begin(void) const { return 0;}; accelerator_inline uint64_t end(void) const { return this->_odata_size; }; accelerator_inline uint64_t size(void) const { return this->_odata_size; }; diff --git a/Grid/lattice/Lattice_where.h b/Grid/lattice/Lattice_where.h index 6686d1b3..777f4015 100644 --- a/Grid/lattice/Lattice_where.h +++ b/Grid/lattice/Lattice_where.h @@ -43,7 +43,7 @@ inline void whereWolf(Lattice &ret,const Lattice &predicate,Lattice< conformable(iftrue,predicate); conformable(iftrue,ret); - GridBase *grid=iftrue._grid; + GridBase *grid=iftrue.Grid(); typedef typename vobj::scalar_object scalar_object; typedef typename vobj::scalar_type scalar_type; @@ -52,22 +52,23 @@ inline void whereWolf(Lattice &ret,const Lattice &predicate,Lattice< const int Nsimd = grid->Nsimd(); - std::vector mask(Nsimd); - std::vector truevals (Nsimd); - std::vector falsevals(Nsimd); - - parallel_for(int ss=0;ssoSites(); ss++){ - - extract(iftrue._odata[ss] ,truevals); - extract(iffalse._odata[ss] ,falsevals); - extract(TensorRemove(predicate._odata[ss]),mask); - - for(int s=0;soSites(); + thread_for(ss,NN,{ + Integer mask; + scalar_object trueval; + scalar_object falseval; + for(int l=0;l @@ -76,9 +77,9 @@ inline Lattice whereWolf(const Lattice &predicate,Lattice &ift conformable(iftrue,iffalse); conformable(iftrue,predicate); - Lattice ret(iftrue._grid); + Lattice ret(iftrue.Grid()); - where(ret,predicate,iftrue,iffalse); + whereWolf(ret,predicate,iftrue,iffalse); return ret; } diff --git a/Grid/qcd/action/ActionBase.h b/Grid/qcd/action/ActionBase.h index bff21d1d..17980ee0 100644 --- a/Grid/qcd/action/ActionBase.h +++ b/Grid/qcd/action/ActionBase.h @@ -41,7 +41,7 @@ class Action public: bool is_smeared = false; // 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 void deriv(const GaugeField& U, GaugeField& dSdU) = 0; // evaluate the action derivative virtual std::string action_name() = 0; // return the action name diff --git a/Grid/qcd/action/fermion/FermionOperatorImpl.h b/Grid/qcd/action/fermion/FermionOperatorImpl.h index b444f6dc..9345c0e6 100644 --- a/Grid/qcd/action/fermion/FermionOperatorImpl.h +++ b/Grid/qcd/action/fermion/FermionOperatorImpl.h @@ -153,8 +153,8 @@ public: typedef typename Impl::StencilImpl StencilImpl; \ typedef typename Impl::ImplParams ImplParams; \ typedef typename Impl::StencilImpl::View_type StencilView; \ - typedef typename ViewMap::Type FermionFieldView; \ - typedef typename ViewMap::Type DoubledGaugeFieldView; + typedef const typename ViewMap::Type FermionFieldView; \ + typedef const typename ViewMap::Type DoubledGaugeFieldView; #define INHERIT_IMPL_TYPES(Base) \ INHERIT_GIMPL_TYPES(Base) \ diff --git a/Grid/qcd/action/fermion/MADWF.h b/Grid/qcd/action/fermion/MADWF.h index 6b3c6e71..5d17e865 100644 --- a/Grid/qcd/action/fermion/MADWF.h +++ b/Grid/qcd/action/fermion/MADWF.h @@ -85,7 +85,7 @@ class MADWF maxiter =_maxiter; }; - void operator() (const FermionFieldo &src4,FermionFieldo &sol5) + void operator() (const FermionFieldo &src,FermionFieldo &sol5) { std::cout << GridLogMessage<< " ************************************************" << std::endl; std::cout << GridLogMessage<< " MADWF-like algorithm " << std::endl; @@ -114,8 +114,16 @@ class MADWF /////////////////////////////////////// //Import source, include Dminus factors /////////////////////////////////////// - Mato.ImportPhysicalFermionSource(src4,b); - std::cout << GridLogMessage << " src4 " < + static accelerator_inline void multLink(_Spinor &phi, const SiteDoubledGaugeField &U, - const SiteSpinor &chi, + const _Spinor &chi, int mu) { - mult(&phi(), &U(mu), &chi()); + auto UU = coalescedRead(U(mu)); + mult(&phi(), &UU, &chi()); } - static accelerator_inline void multLinkAdd(SiteSpinor &phi, + template + static accelerator_inline void multLinkAdd(_Spinor &phi, const SiteDoubledGaugeField &U, - const SiteSpinor &chi, + const _Spinor &chi, int mu) { - mac(&phi(), &U(mu), &chi()); + auto UU = coalescedRead(U(mu)); + mac(&phi(), &UU, &chi()); } template diff --git a/Grid/qcd/action/fermion/WilsonCloverFermion.h b/Grid/qcd/action/fermion/WilsonCloverFermion.h index 91ad6d6d..92af7111 100644 --- a/Grid/qcd/action/fermion/WilsonCloverFermion.h +++ b/Grid/qcd/action/fermion/WilsonCloverFermion.h @@ -245,7 +245,7 @@ public: return out; } -private: +protected: // here fixing the 4 dimensions, make it more general? RealD csw_r; // Clover coefficient - spatial diff --git a/Grid/qcd/action/fermion/WilsonCompressor.h b/Grid/qcd/action/fermion/WilsonCompressor.h index 10e98f33..0760bcba 100644 --- a/Grid/qcd/action/fermion/WilsonCompressor.h +++ b/Grid/qcd/action/fermion/WilsonCompressor.h @@ -61,7 +61,7 @@ public: typedef typename SiteHalfSpinor::vector_type vComplexHigh; constexpr static int Nw=sizeof(SiteHalfSpinor)/sizeof(vComplexHigh); - accelerator_inline int CommDatumSize(void) { + accelerator_inline int CommDatumSize(void) const { return sizeof(SiteHalfCommSpinor); } @@ -69,7 +69,7 @@ public: /* Compress includes precision change if mpi data is not same */ /*****************************************************/ template - 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; projector::Proj(tmp,in,mu,dag); vstream(buf[o],tmp); @@ -81,7 +81,7 @@ public: accelerator_inline void Exchange(SiteHalfSpinor *mp, const SiteHalfSpinor * __restrict__ vp0, const SiteHalfSpinor * __restrict__ vp1, - Integer type,Integer o){ + Integer type,Integer o) const { SiteHalfSpinor tmp1; SiteHalfSpinor tmp2; exchange(tmp1,tmp2,vp0[o],vp1[o],type); @@ -93,7 +93,7 @@ public: /* Have a decompression step if mpi data is not same */ /*****************************************************/ accelerator_inline void Decompress(SiteHalfSpinor * __restrict__ out, - SiteHalfSpinor * __restrict__ in, Integer o) { + SiteHalfSpinor * __restrict__ in, Integer o) const { assert(0); } @@ -103,7 +103,7 @@ public: accelerator_inline void CompressExchange(SiteHalfSpinor * __restrict__ out0, SiteHalfSpinor * __restrict__ out1, 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 temp3, temp4; @@ -117,7 +117,7 @@ public: /*****************************************************/ /* 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; constexpr static int Nw=sizeof(SiteHalfSpinor)/sizeof(vComplexHigh); - accelerator_inline int CommDatumSize(void) { + accelerator_inline int CommDatumSize(void) const { return sizeof(SiteHalfCommSpinor); } @@ -150,7 +150,7 @@ public: /* Compress includes precision change if mpi data is not same */ /*****************************************************/ template - 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; SiteHalfCommSpinor *hbuf = (SiteHalfCommSpinor *)buf; projector::Proj(hsp,in,mu,dag); @@ -163,7 +163,7 @@ public: accelerator_inline void Exchange(SiteHalfSpinor *mp, SiteHalfSpinor *vp0, SiteHalfSpinor *vp1, - Integer type,Integer o){ + Integer type,Integer o) const { SiteHalfSpinor vt0,vt1; SiteHalfCommSpinor *vpp0 = (SiteHalfCommSpinor *)vp0; SiteHalfCommSpinor *vpp1 = (SiteHalfCommSpinor *)vp1; @@ -175,7 +175,7 @@ public: /*****************************************************/ /* 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; precisionChange((vComplexHigh *)&out[o],(vComplexLow *)&hin[o],Nw); } @@ -186,7 +186,7 @@ public: accelerator_inline void CompressExchange(SiteHalfSpinor *out0, SiteHalfSpinor *out1, 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; SiteHalfCommSpinor *hout0 = (SiteHalfCommSpinor *)out0; SiteHalfCommSpinor *hout1 = (SiteHalfCommSpinor *)out1; @@ -200,7 +200,7 @@ public: /*****************************************************/ /* Pass the info to the stencil */ /*****************************************************/ - accelerator_inline bool DecompressionStep(void) { return true; } + accelerator_inline bool DecompressionStep(void) const { return true; } }; diff --git a/Grid/qcd/action/fermion/WilsonImpl.h b/Grid/qcd/action/fermion/WilsonImpl.h index d7941d1f..2ff6feba 100644 --- a/Grid/qcd/action/fermion/WilsonImpl.h +++ b/Grid/qcd/action/fermion/WilsonImpl.h @@ -72,7 +72,7 @@ public: typedef WilsonCompressor Compressor; typedef WilsonImplParams ImplParams; typedef WilsonStencil StencilImpl; - typedef typename StencilImpl::View_type StencilView; + typedef const typename StencilImpl::View_type StencilView; ImplParams Params; @@ -184,18 +184,22 @@ public: mat = TraceIndex(P); } - inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds){ + inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds) + { for (int mu = 0; mu < Nd; mu++) mat[mu] = PeekIndex(Uds, mu); } - - inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu){ - + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu) + { +#undef USE_OLD_INSERT_FORCE int Ls=Btilde.Grid()->_fdimensions[0]; + autoView( mat_v , mat, AcceleratorWrite); +#ifdef USE_OLD_INSERT_FORCE GaugeLinkField tmp(mat.Grid()); tmp = Zero(); { + const int Nsimd = SiteSpinor::Nsimd(); autoView( tmp_v , tmp, AcceleratorWrite); autoView( Btilde_v , Btilde, AcceleratorRead); autoView( Atilde_v , Atilde, AcceleratorRead); @@ -208,6 +212,29 @@ public: }); } PokeIndex(mat,tmp,mu); +#else + { + const int Nsimd = SiteSpinor::Nsimd(); + autoView( Btilde_v , Btilde, AcceleratorRead); + autoView( Atilde_v , Atilde, AcceleratorRead); + accelerator_for(sss,mat.Grid()->oSites(),Nsimd,{ + int sU=sss; + typedef decltype(coalescedRead(mat_v[sU](mu)() )) ColorMatrixType; + ColorMatrixType sum; + zeroit(sum); + for(int s=0;s Base; - + typedef AcceleratorVector StencilVector; public: +#ifdef GRID_SYCL +#define SYCL_HACK +#endif +#ifdef SYCL_HACK + static void HandDhopSiteSycl(StencilVector st_perm,StencilEntry *st_p, SiteDoubledGaugeField *U,SiteHalfSpinor *buf, + int ss,int sU,const SiteSpinor *in, SiteSpinor *out); +#endif + static void DhopKernel(int Opt,StencilImpl &st, DoubledGaugeField &U, SiteHalfSpinor * buf, int Ls, int Nsite, const FermionField &in, FermionField &out, int interior=1,int exterior=1) ; diff --git a/Grid/qcd/action/fermion/implementation/StaggeredKernelsImplementation.h b/Grid/qcd/action/fermion/implementation/StaggeredKernelsImplementation.h index 0b6f9fb0..dd62e109 100644 --- a/Grid/qcd/action/fermion/implementation/StaggeredKernelsImplementation.h +++ b/Grid/qcd/action/fermion/implementation/StaggeredKernelsImplementation.h @@ -35,39 +35,32 @@ NAMESPACE_BEGIN(Grid); #define GENERIC_STENCIL_LEG(U,Dir,skew,multLink) \ SE = st.GetEntry(ptype, Dir+skew, sF); \ if (SE->_is_local ) { \ - if (SE->_permute) { \ - chi_p = χ \ - permute(chi, in[SE->_offset], ptype); \ - } else { \ - chi_p = &in[SE->_offset]; \ - } \ + int perm= SE->_permute; \ + chi = coalescedReadPermute(in[SE->_offset],ptype,perm,lane);\ } else { \ - chi_p = &buf[SE->_offset]; \ + chi = coalescedRead(buf[SE->_offset],lane); \ } \ - multLink(Uchi, U[sU], *chi_p, Dir); + acceleratorSynchronise(); \ + multLink(Uchi, U[sU], chi, Dir); #define GENERIC_STENCIL_LEG_INT(U,Dir,skew,multLink) \ SE = st.GetEntry(ptype, Dir+skew, sF); \ if (SE->_is_local ) { \ - if (SE->_permute) { \ - chi_p = χ \ - permute(chi, in[SE->_offset], ptype); \ - } else { \ - chi_p = &in[SE->_offset]; \ - } \ + int perm= SE->_permute; \ + chi = coalescedReadPermute(in[SE->_offset],ptype,perm,lane);\ } else if ( st.same_node[Dir] ) { \ - chi_p = &buf[SE->_offset]; \ + chi = coalescedRead(buf[SE->_offset],lane); \ } \ if (SE->_is_local || st.same_node[Dir] ) { \ - multLink(Uchi, U[sU], *chi_p, Dir); \ + multLink(Uchi, U[sU], chi, Dir); \ } #define GENERIC_STENCIL_LEG_EXT(U,Dir,skew,multLink) \ SE = st.GetEntry(ptype, Dir+skew, sF); \ if ((!SE->_is_local) && (!st.same_node[Dir]) ) { \ nmu++; \ - chi_p = &buf[SE->_offset]; \ - multLink(Uchi, U[sU], *chi_p, Dir); \ + chi = coalescedRead(buf[SE->_offset],lane); \ + multLink(Uchi, U[sU], chi, Dir); \ } template @@ -84,12 +77,14 @@ void StaggeredKernels::DhopSiteGeneric(StencilView &st, SiteSpinor *buf, int sF, int sU, const FermionFieldView &in, FermionFieldView &out, int dag) { - const SiteSpinor *chi_p; - SiteSpinor chi; - SiteSpinor Uchi; + typedef decltype(coalescedRead(in[0])) calcSpinor; + calcSpinor chi; + calcSpinor Uchi; StencilEntry *SE; int ptype; int skew; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); // for(int s=0;s::DhopSiteGeneric(StencilView &st, if ( dag ) { Uchi = - Uchi; } - vstream(out[sF], Uchi); + coalescedWrite(out[sF], Uchi,lane); } }; @@ -130,13 +125,16 @@ template accelerator_inline void StaggeredKernels::DhopSiteGenericInt(StencilView &st, DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU, SiteSpinor *buf, int sF, int sU, - const FermionFieldView &in, FermionFieldView &out,int dag) { - const SiteSpinor *chi_p; - SiteSpinor chi; - SiteSpinor Uchi; + const FermionFieldView &in, FermionFieldView &out,int dag) +{ + typedef decltype(coalescedRead(in[0])) calcSpinor; + calcSpinor chi; + calcSpinor Uchi; StencilEntry *SE; int ptype; int skew ; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); // for(int s=0;s::DhopSiteGenericInt(StencilView &st, if ( dag ) { Uchi = - Uchi; } - vstream(out[sF], Uchi); + coalescedWrite(out[sF], Uchi,lane); } }; @@ -178,14 +176,17 @@ template accelerator_inline void StaggeredKernels::DhopSiteGenericExt(StencilView &st, DoubledGaugeFieldView &U, DoubledGaugeFieldView &UUU, SiteSpinor *buf, int sF, int sU, - const FermionFieldView &in, FermionFieldView &out,int dag) { - const SiteSpinor *chi_p; - // SiteSpinor chi; - SiteSpinor Uchi; + const FermionFieldView &in, FermionFieldView &out,int dag) +{ + typedef decltype(coalescedRead(in[0])) calcSpinor; + calcSpinor chi; + calcSpinor Uchi; StencilEntry *SE; int ptype; int nmu=0; int skew ; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); // for(int s=0;s::DhopSiteGenericExt(StencilView &st, GENERIC_STENCIL_LEG_EXT(UUU,Zm,skew,Impl::multLinkAdd); GENERIC_STENCIL_LEG_EXT(UUU,Tm,skew,Impl::multLinkAdd); } - if ( nmu ) { - if ( dag ) { - out[sF] = out[sF] - Uchi; + if ( nmu ) { + auto _out = coalescedRead(out[sF],lane); + if ( dag ) { + coalescedWrite(out[sF], _out-Uchi,lane); } else { - out[sF] = out[sF] + Uchi; + coalescedWrite(out[sF], _out+Uchi,lane); } } } @@ -261,6 +263,8 @@ void StaggeredKernels::DhopImproved(StencilImpl &st, LebesgueOrder &lo, GridBase *FGrid=in.Grid(); GridBase *UGrid=U.Grid(); typedef StaggeredKernels ThisKernel; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); autoView( UUU_v , UUU, AcceleratorRead); autoView( U_v , U, AcceleratorRead); autoView( in_v , in, AcceleratorRead); @@ -301,6 +305,8 @@ void StaggeredKernels::DhopNaive(StencilImpl &st, LebesgueOrder &lo, GridBase *FGrid=in.Grid(); GridBase *UGrid=U.Grid(); typedef StaggeredKernels ThisKernel; + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); autoView( UUU_v , U, AcceleratorRead); autoView( U_v , U, AcceleratorRead); autoView( in_v , in, AcceleratorRead); diff --git a/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h index 4977ea68..84ac25c1 100644 --- a/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h +++ b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h @@ -397,6 +397,7 @@ void WilsonFermion::DhopDerivEO(GaugeField &mat, const FermionField &U, co template void WilsonFermion::Dhop(const FermionField &in, FermionField &out, int dag) { + DhopCalls+=2; conformable(in.Grid(), _grid); // verifies full grid conformable(in.Grid(), out.Grid()); @@ -408,6 +409,7 @@ void WilsonFermion::Dhop(const FermionField &in, FermionField &out, int da template void WilsonFermion::DhopOE(const FermionField &in, FermionField &out, int dag) { + DhopCalls++; conformable(in.Grid(), _cbgrid); // verifies half grid conformable(in.Grid(), out.Grid()); // drops the cb check @@ -420,6 +422,7 @@ void WilsonFermion::DhopOE(const FermionField &in, FermionField &out, int template void WilsonFermion::DhopEO(const FermionField &in, FermionField &out,int dag) { + DhopCalls++; conformable(in.Grid(), _cbgrid); // verifies half grid conformable(in.Grid(), out.Grid()); // drops the cb check diff --git a/Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h index 89ae5668..0703b613 100644 --- a/Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h +++ b/Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h @@ -76,7 +76,24 @@ Author: paboyle #define REGISTER -#define LOAD_CHIMU \ +#ifdef GRID_SIMT +#define LOAD_CHIMU(ptype) \ + {const SiteSpinor & ref (in[offset]); \ + Chimu_00=coalescedReadPermute(ref()(0)(0),perm,lane); \ + Chimu_01=coalescedReadPermute(ref()(0)(1),perm,lane); \ + Chimu_02=coalescedReadPermute(ref()(0)(2),perm,lane); \ + Chimu_10=coalescedReadPermute(ref()(1)(0),perm,lane); \ + Chimu_11=coalescedReadPermute(ref()(1)(1),perm,lane); \ + Chimu_12=coalescedReadPermute(ref()(1)(2),perm,lane); \ + Chimu_20=coalescedReadPermute(ref()(2)(0),perm,lane); \ + Chimu_21=coalescedReadPermute(ref()(2)(1),perm,lane); \ + Chimu_22=coalescedReadPermute(ref()(2)(2),perm,lane); \ + Chimu_30=coalescedReadPermute(ref()(3)(0),perm,lane); \ + Chimu_31=coalescedReadPermute(ref()(3)(1),perm,lane); \ + Chimu_32=coalescedReadPermute(ref()(3)(2),perm,lane); } +#define PERMUTE_DIR(dir) ; +#else +#define LOAD_CHIMU(ptype) \ {const SiteSpinor & ref (in[offset]); \ Chimu_00=ref()(0)(0);\ Chimu_01=ref()(0)(1);\ @@ -91,55 +108,55 @@ Author: paboyle Chimu_31=ref()(3)(1);\ Chimu_32=ref()(3)(2);} -#define LOAD_CHI\ - {const SiteHalfSpinor &ref(buf[offset]); \ - Chi_00 = ref()(0)(0);\ - Chi_01 = ref()(0)(1);\ - Chi_02 = ref()(0)(2);\ - Chi_10 = ref()(1)(0);\ - Chi_11 = ref()(1)(1);\ - Chi_12 = ref()(1)(2);} - -// To splat or not to splat depends on the implementation -#define MULT_2SPIN(A)\ - {auto & ref(U[sU](A)); \ - Impl::loadLinkElement(U_00,ref()(0,0)); \ - Impl::loadLinkElement(U_10,ref()(1,0)); \ - Impl::loadLinkElement(U_20,ref()(2,0)); \ - Impl::loadLinkElement(U_01,ref()(0,1)); \ - Impl::loadLinkElement(U_11,ref()(1,1)); \ - Impl::loadLinkElement(U_21,ref()(2,1)); \ - UChi_00 = U_00*Chi_00;\ - UChi_10 = U_00*Chi_10;\ - UChi_01 = U_10*Chi_00;\ - UChi_11 = U_10*Chi_10;\ - UChi_02 = U_20*Chi_00;\ - UChi_12 = U_20*Chi_10;\ - UChi_00+= U_01*Chi_01;\ - UChi_10+= U_01*Chi_11;\ - UChi_01+= U_11*Chi_01;\ - UChi_11+= U_11*Chi_11;\ - UChi_02+= U_21*Chi_01;\ - UChi_12+= U_21*Chi_11;\ - Impl::loadLinkElement(U_00,ref()(0,2)); \ - Impl::loadLinkElement(U_10,ref()(1,2)); \ - Impl::loadLinkElement(U_20,ref()(2,2)); \ - UChi_00+= U_00*Chi_02;\ - UChi_10+= U_00*Chi_12;\ - UChi_01+= U_10*Chi_02;\ - UChi_11+= U_10*Chi_12;\ - UChi_02+= U_20*Chi_02;\ - UChi_12+= U_20*Chi_12;} - - #define PERMUTE_DIR(dir) \ - permute##dir(Chi_00,Chi_00);\ + permute##dir(Chi_00,Chi_00); \ permute##dir(Chi_01,Chi_01);\ permute##dir(Chi_02,Chi_02);\ - permute##dir(Chi_10,Chi_10);\ + permute##dir(Chi_10,Chi_10); \ permute##dir(Chi_11,Chi_11);\ permute##dir(Chi_12,Chi_12); +#endif + +#define MULT_2SPIN(A)\ + {auto & ref(U[sU](A)); \ + U_00=coalescedRead(ref()(0,0),lane); \ + U_10=coalescedRead(ref()(1,0),lane); \ + U_20=coalescedRead(ref()(2,0),lane); \ + U_01=coalescedRead(ref()(0,1),lane); \ + U_11=coalescedRead(ref()(1,1),lane); \ + U_21=coalescedRead(ref()(2,1),lane); \ + UChi_00 = U_00*Chi_00; \ + UChi_10 = U_00*Chi_10; \ + UChi_01 = U_10*Chi_00; \ + UChi_11 = U_10*Chi_10; \ + UChi_02 = U_20*Chi_00; \ + UChi_12 = U_20*Chi_10; \ + UChi_00+= U_01*Chi_01; \ + UChi_10+= U_01*Chi_11; \ + UChi_01+= U_11*Chi_01; \ + UChi_11+= U_11*Chi_11; \ + UChi_02+= U_21*Chi_01; \ + UChi_12+= U_21*Chi_11; \ + U_00=coalescedRead(ref()(0,2),lane); \ + U_10=coalescedRead(ref()(1,2),lane); \ + U_20=coalescedRead(ref()(2,2),lane); \ + UChi_00+= U_00*Chi_02; \ + UChi_10+= U_00*Chi_12; \ + UChi_01+= U_10*Chi_02; \ + UChi_11+= U_10*Chi_12; \ + UChi_02+= U_20*Chi_02; \ + UChi_12+= U_20*Chi_12;} + +#define LOAD_CHI \ + {const SiteHalfSpinor &ref(buf[offset]); \ + Chi_00 = coalescedRead(ref()(0)(0),lane); \ + Chi_01 = coalescedRead(ref()(0)(1),lane); \ + Chi_02 = coalescedRead(ref()(0)(2),lane); \ + Chi_10 = coalescedRead(ref()(1)(0),lane); \ + Chi_11 = coalescedRead(ref()(1)(1),lane); \ + Chi_12 = coalescedRead(ref()(1)(2),lane);} + // hspin(0)=fspin(0)+timesI(fspin(3)); // hspin(1)=fspin(1)+timesI(fspin(2)); #define XP_PROJ \ @@ -353,13 +370,13 @@ Author: paboyle result_31-= UChi_11; \ result_32-= UChi_12; -#define HAND_STENCIL_LEG(PROJ,PERM,DIR,RECON) \ +#define HAND_STENCIL_LEGB(PROJ,PERM,DIR,RECON) \ SE=st.GetEntry(ptype,DIR,ss); \ offset = SE->_offset; \ local = SE->_is_local; \ perm = SE->_permute; \ if ( local ) { \ - LOAD_CHIMU; \ + LOAD_CHIMU(PERM); \ PROJ; \ if ( perm) { \ PERMUTE_DIR(PERM); \ @@ -367,6 +384,37 @@ Author: paboyle } else { \ LOAD_CHI; \ } \ + acceleratorSynchronise(); \ + MULT_2SPIN(DIR); \ + RECON; + +#define HAND_STENCIL_LEG(PROJ,PERM,DIR,RECON) \ + SE=&st_p[DIR+8*ss]; \ + ptype=st_perm[DIR]; \ + offset = SE->_offset; \ + local = SE->_is_local; \ + perm = SE->_permute; \ + if ( local ) { \ + LOAD_CHIMU(PERM); \ + PROJ; \ + if ( perm) { \ + PERMUTE_DIR(PERM); \ + } \ + } else { \ + LOAD_CHI; \ + } \ + acceleratorSynchronise(); \ + MULT_2SPIN(DIR); \ + RECON; + +#define HAND_STENCIL_LEGA(PROJ,PERM,DIR,RECON) \ + SE=&st_p[DIR+8*ss]; \ + ptype=st_perm[DIR]; \ + /*SE=st.GetEntry(ptype,DIR,ss);*/ \ + offset = SE->_offset; \ + perm = SE->_permute; \ + LOAD_CHIMU(PERM); \ + PROJ; \ MULT_2SPIN(DIR); \ RECON; @@ -376,7 +424,7 @@ Author: paboyle local = SE->_is_local; \ perm = SE->_permute; \ if ( local ) { \ - LOAD_CHIMU; \ + LOAD_CHIMU(PERM); \ PROJ; \ if ( perm) { \ PERMUTE_DIR(PERM); \ @@ -384,10 +432,12 @@ Author: paboyle } else if ( st.same_node[DIR] ) { \ LOAD_CHI; \ } \ + acceleratorSynchronise(); \ if (local || st.same_node[DIR] ) { \ MULT_2SPIN(DIR); \ RECON; \ - } + } \ + acceleratorSynchronise(); #define HAND_STENCIL_LEG_EXT(PROJ,PERM,DIR,RECON) \ SE=st.GetEntry(ptype,DIR,ss); \ @@ -397,44 +447,44 @@ Author: paboyle MULT_2SPIN(DIR); \ RECON; \ nmu++; \ - } + } \ + acceleratorSynchronise(); #define HAND_RESULT(ss) \ { \ - SiteSpinor & ref (out[ss]); \ - vstream(ref()(0)(0),result_00); \ - vstream(ref()(0)(1),result_01); \ - vstream(ref()(0)(2),result_02); \ - vstream(ref()(1)(0),result_10); \ - vstream(ref()(1)(1),result_11); \ - vstream(ref()(1)(2),result_12); \ - vstream(ref()(2)(0),result_20); \ - vstream(ref()(2)(1),result_21); \ - vstream(ref()(2)(2),result_22); \ - vstream(ref()(3)(0),result_30); \ - vstream(ref()(3)(1),result_31); \ - vstream(ref()(3)(2),result_32); \ + SiteSpinor & ref (out[ss]); \ + coalescedWrite(ref()(0)(0),result_00,lane); \ + coalescedWrite(ref()(0)(1),result_01,lane); \ + coalescedWrite(ref()(0)(2),result_02,lane); \ + coalescedWrite(ref()(1)(0),result_10,lane); \ + coalescedWrite(ref()(1)(1),result_11,lane); \ + coalescedWrite(ref()(1)(2),result_12,lane); \ + coalescedWrite(ref()(2)(0),result_20,lane); \ + coalescedWrite(ref()(2)(1),result_21,lane); \ + coalescedWrite(ref()(2)(2),result_22,lane); \ + coalescedWrite(ref()(3)(0),result_30,lane); \ + coalescedWrite(ref()(3)(1),result_31,lane); \ + coalescedWrite(ref()(3)(2),result_32,lane); \ } -#define HAND_RESULT_EXT(ss) \ - if (nmu){ \ - SiteSpinor & ref (out[ss]); \ - ref()(0)(0)+=result_00; \ - ref()(0)(1)+=result_01; \ - ref()(0)(2)+=result_02; \ - ref()(1)(0)+=result_10; \ - ref()(1)(1)+=result_11; \ - ref()(1)(2)+=result_12; \ - ref()(2)(0)+=result_20; \ - ref()(2)(1)+=result_21; \ - ref()(2)(2)+=result_22; \ - ref()(3)(0)+=result_30; \ - ref()(3)(1)+=result_31; \ - ref()(3)(2)+=result_32; \ +#define HAND_RESULT_EXT(ss) \ + { \ + SiteSpinor & ref (out[ss]); \ + coalescedWrite(ref()(0)(0),coalescedRead(ref()(0)(0))+result_00,lane); \ + coalescedWrite(ref()(0)(1),coalescedRead(ref()(0)(1))+result_01,lane); \ + coalescedWrite(ref()(0)(2),coalescedRead(ref()(0)(2))+result_02,lane); \ + coalescedWrite(ref()(1)(0),coalescedRead(ref()(1)(0))+result_10,lane); \ + coalescedWrite(ref()(1)(1),coalescedRead(ref()(1)(1))+result_11,lane); \ + coalescedWrite(ref()(1)(2),coalescedRead(ref()(1)(2))+result_12,lane); \ + coalescedWrite(ref()(2)(0),coalescedRead(ref()(2)(0))+result_20,lane); \ + coalescedWrite(ref()(2)(1),coalescedRead(ref()(2)(1))+result_21,lane); \ + coalescedWrite(ref()(2)(2),coalescedRead(ref()(2)(2))+result_22,lane); \ + coalescedWrite(ref()(3)(0),coalescedRead(ref()(3)(0))+result_30,lane); \ + coalescedWrite(ref()(3)(1),coalescedRead(ref()(3)(1))+result_31,lane); \ + coalescedWrite(ref()(3)(2),coalescedRead(ref()(3)(2))+result_32,lane); \ } - -#define HAND_DECLARATIONS(a) \ +#define HAND_DECLARATIONS(Simd) \ Simd result_00; \ Simd result_01; \ Simd result_02; \ @@ -466,19 +516,19 @@ Author: paboyle Simd U_11; \ Simd U_21; -#define ZERO_RESULT \ - result_00=Zero(); \ - result_01=Zero(); \ - result_02=Zero(); \ - result_10=Zero(); \ - result_11=Zero(); \ - result_12=Zero(); \ - result_20=Zero(); \ - result_21=Zero(); \ - result_22=Zero(); \ - result_30=Zero(); \ - result_31=Zero(); \ - result_32=Zero(); +#define ZERO_RESULT \ + zeroit(result_00); \ + zeroit(result_01); \ + zeroit(result_02); \ + zeroit(result_10); \ + zeroit(result_11); \ + zeroit(result_12); \ + zeroit(result_20); \ + zeroit(result_21); \ + zeroit(result_22); \ + zeroit(result_30); \ + zeroit(result_31); \ + zeroit(result_32); #define Chimu_00 Chi_00 #define Chimu_01 Chi_01 @@ -495,15 +545,53 @@ Author: paboyle NAMESPACE_BEGIN(Grid); + +#ifdef SYCL_HACK template accelerator_inline void -WilsonKernels::HandDhopSite(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf, - int ss,int sU,const FermionFieldView &in, FermionFieldView &out) +WilsonKernels::HandDhopSiteSycl(StencilVector st_perm,StencilEntry *st_p, SiteDoubledGaugeField *U,SiteHalfSpinor *buf, + int ss,int sU,const SiteSpinor *in, SiteSpinor *out) { // T==0, Z==1, Y==2, Z==3 expect 1,2,2,2 simd layout etc... typedef typename Simd::scalar_type S; typedef typename Simd::vector_type V; + typedef iSinglet vCplx; + // typedef decltype( coalescedRead( vCplx()()() )) Simt; + typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt; - HAND_DECLARATIONS(ignore); + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); + + HAND_DECLARATIONS(Simt); + + int offset,local,perm, ptype; + StencilEntry *SE; + HAND_STENCIL_LEG(XM_PROJ,3,Xp,XM_RECON); + HAND_STENCIL_LEG(YM_PROJ,2,Yp,YM_RECON_ACCUM); + HAND_STENCIL_LEG(ZM_PROJ,1,Zp,ZM_RECON_ACCUM); + HAND_STENCIL_LEG(TM_PROJ,0,Tp,TM_RECON_ACCUM); + HAND_STENCIL_LEG(XP_PROJ,3,Xm,XP_RECON_ACCUM); + HAND_STENCIL_LEG(YP_PROJ,2,Ym,YP_RECON_ACCUM); + HAND_STENCIL_LEG(ZP_PROJ,1,Zm,ZP_RECON_ACCUM); + HAND_STENCIL_LEG(TP_PROJ,0,Tm,TP_RECON_ACCUM); + HAND_RESULT(ss); +} +#endif + +template accelerator_inline void +WilsonKernels::HandDhopSite(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor *buf, + int ss,int sU,const FermionFieldView &in, FermionFieldView &out) +{ + auto st_p = st._entries_p; + auto st_perm = st._permute_type; +// T==0, Z==1, Y==2, Z==3 expect 1,2,2,2 simd layout etc... + typedef typename Simd::scalar_type S; + typedef typename Simd::vector_type V; + typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt; + + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); + + HAND_DECLARATIONS(Simt); int offset,local,perm, ptype; StencilEntry *SE; @@ -523,10 +611,16 @@ template accelerator_inline void WilsonKernels::HandDhopSiteDag(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf, int ss,int sU,const FermionFieldView &in, FermionFieldView &out) { + auto st_p = st._entries_p; + auto st_perm = st._permute_type; typedef typename Simd::scalar_type S; typedef typename Simd::vector_type V; + typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt; - HAND_DECLARATIONS(ignore); + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); + + HAND_DECLARATIONS(Simt); StencilEntry *SE; int offset,local,perm, ptype; @@ -546,11 +640,17 @@ template accelerator_inline void WilsonKernels::HandDhopSiteInt(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf, int ss,int sU,const FermionFieldView &in, FermionFieldView &out) { + auto st_p = st._entries_p; + auto st_perm = st._permute_type; // T==0, Z==1, Y==2, Z==3 expect 1,2,2,2 simd layout etc... typedef typename Simd::scalar_type S; typedef typename Simd::vector_type V; + typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt; - HAND_DECLARATIONS(ignore); + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); + + HAND_DECLARATIONS(Simt); int offset,local,perm, ptype; StencilEntry *SE; @@ -570,10 +670,16 @@ template accelerator_inline void WilsonKernels::HandDhopSiteDagInt(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf, int ss,int sU,const FermionFieldView &in, FermionFieldView &out) { + auto st_p = st._entries_p; + auto st_perm = st._permute_type; typedef typename Simd::scalar_type S; typedef typename Simd::vector_type V; + typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt; - HAND_DECLARATIONS(ignore); + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); + + HAND_DECLARATIONS(Simt); StencilEntry *SE; int offset,local,perm, ptype; @@ -593,11 +699,17 @@ template accelerator_inline void WilsonKernels::HandDhopSiteExt(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf, int ss,int sU,const FermionFieldView &in, FermionFieldView &out) { + auto st_p = st._entries_p; + auto st_perm = st._permute_type; // T==0, Z==1, Y==2, Z==3 expect 1,2,2,2 simd layout etc... typedef typename Simd::scalar_type S; typedef typename Simd::vector_type V; + typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt; - HAND_DECLARATIONS(ignore); + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); + + HAND_DECLARATIONS(Simt); int offset, ptype; StencilEntry *SE; @@ -618,10 +730,16 @@ template accelerator_inline void WilsonKernels::HandDhopSiteDagExt(StencilView &st,DoubledGaugeFieldView &U,SiteHalfSpinor *buf, int ss,int sU,const FermionFieldView &in, FermionFieldView &out) { + auto st_p = st._entries_p; + auto st_perm = st._permute_type; typedef typename Simd::scalar_type S; typedef typename Simd::vector_type V; + typedef decltype( coalescedRead( in[0]()(0)(0) )) Simt; - HAND_DECLARATIONS(ignore); + const int Nsimd = SiteHalfSpinor::Nsimd(); + const int lane=acceleratorSIMTlane(Nsimd); + + HAND_DECLARATIONS(Simt); StencilEntry *SE; int offset, ptype; @@ -682,3 +800,4 @@ NAMESPACE_END(Grid); #undef HAND_RESULT #undef HAND_RESULT_INT #undef HAND_RESULT_EXT +#undef HAND_DECLARATIONS diff --git a/Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h index c5f50bbb..9228b84c 100644 --- a/Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h +++ b/Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h @@ -416,7 +416,21 @@ void WilsonKernels::DhopDirKernel( StencilImpl &st, DoubledGaugeField &U,S #undef LoopBody } -#define KERNEL_CALLNB(A) \ +#define KERNEL_CALL_TMP(A) \ + const uint64_t NN = Nsite*Ls; \ + auto U_p = & U_v[0]; \ + auto in_p = & in_v[0]; \ + auto out_p = & out_v[0]; \ + auto st_p = st_v._entries_p; \ + auto st_perm = st_v._permute_type; \ + accelerator_forNB( ss, NN, Simd::Nsimd(), { \ + int sF = ss; \ + int sU = ss/Ls; \ + WilsonKernels::A(st_perm,st_p,U_p,buf,sF,sU,in_p,out_p); \ + }); \ + accelerator_barrier(); + +#define KERNEL_CALLNB(A) \ const uint64_t NN = Nsite*Ls; \ accelerator_forNB( ss, NN, Simd::Nsimd(), { \ int sF = ss; \ @@ -445,20 +459,24 @@ void WilsonKernels::DhopKernel(int Opt,StencilImpl &st, DoubledGaugeField if( interior && exterior ) { if (Opt == WilsonKernelsStatic::OptGeneric ) { KERNEL_CALL(GenericDhopSite); return;} -#ifndef GRID_CUDA +#ifdef SYCL_HACK + if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL_TMP(HandDhopSiteSycl); return; } +#else if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSite); return;} +#endif +#ifndef GRID_CUDA if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSite); return;} #endif } else if( interior ) { if (Opt == WilsonKernelsStatic::OptGeneric ) { KERNEL_CALLNB(GenericDhopSiteInt); return;} -#ifndef GRID_CUDA if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALLNB(HandDhopSiteInt); return;} +#ifndef GRID_CUDA if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSiteInt); return;} #endif } else if( exterior ) { if (Opt == WilsonKernelsStatic::OptGeneric ) { KERNEL_CALL(GenericDhopSiteExt); return;} -#ifndef GRID_CUDA if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteExt); return;} +#ifndef GRID_CUDA if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSiteExt); return;} #endif } @@ -476,20 +494,20 @@ void WilsonKernels::DhopKernel(int Opt,StencilImpl &st, DoubledGaugeField if( interior && exterior ) { if (Opt == WilsonKernelsStatic::OptGeneric ) { KERNEL_CALL(GenericDhopSiteDag); return;} -#ifndef GRID_CUDA if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDag); return;} +#ifndef GRID_CUDA if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSiteDag); return;} #endif } else if( interior ) { if (Opt == WilsonKernelsStatic::OptGeneric ) { KERNEL_CALL(GenericDhopSiteDagInt); return;} -#ifndef GRID_CUDA if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDagInt); return;} +#ifndef GRID_CUDA if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSiteDagInt); return;} #endif } else if( exterior ) { if (Opt == WilsonKernelsStatic::OptGeneric ) { KERNEL_CALL(GenericDhopSiteDagExt); return;} -#ifndef GRID_CUDA if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDagExt); return;} +#ifndef GRID_CUDA if (Opt == WilsonKernelsStatic::OptInlineAsm ) { ASM_CALL(AsmDhopSiteDagExt); return;} #endif } diff --git a/Grid/qcd/action/gauge/GaugeImplTypes.h b/Grid/qcd/action/gauge/GaugeImplTypes.h index 55a20eca..2499e0e9 100644 --- a/Grid/qcd/action/gauge/GaugeImplTypes.h +++ b/Grid/qcd/action/gauge/GaugeImplTypes.h @@ -96,7 +96,7 @@ public: /////////////////////////////////////////////////////////// // Move these to another class // 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: // diff --git a/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h b/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h index 639aca19..7690092d 100644 --- a/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h +++ b/Grid/qcd/action/gauge/PlaqPlusRectangleAction.h @@ -49,7 +49,7 @@ public: 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(){ std::stringstream sstream; diff --git a/Grid/qcd/action/gauge/WilsonGaugeAction.h b/Grid/qcd/action/gauge/WilsonGaugeAction.h index 40d600d2..f535b54f 100644 --- a/Grid/qcd/action/gauge/WilsonGaugeAction.h +++ b/Grid/qcd/action/gauge/WilsonGaugeAction.h @@ -54,8 +54,7 @@ public: return sstream.str(); } - 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 RealD S(const GaugeField &U) { RealD plaq = WilsonLoops::avgPlaquette(U); diff --git a/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h b/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h index 9fc0a3b0..576a8cf6 100644 --- a/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h +++ b/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h @@ -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 // - virtual void refresh(const GaugeField& U, GridParallelRNG& pRNG) + virtual void refresh(const GaugeField& U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) { Lop.ImportGauge(U); Rop.ImportGauge(U); diff --git a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRational.h b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRational.h index 56dff94d..656e9b2f 100644 --- a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRational.h +++ b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRational.h @@ -1,4 +1,3 @@ - /************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -43,8 +42,7 @@ NAMESPACE_BEGIN(Grid); // template -class OneFlavourEvenOddRationalPseudoFermionAction - : public Action { +class OneFlavourEvenOddRationalPseudoFermionAction : public Action { public: INHERIT_IMPL_TYPES(Impl); @@ -103,7 +101,7 @@ public: 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} // = e^{- phi^dag (MpcdagMpc)^-1/4 (MpcdagMpc)^-1/4 phi} // Phi = MpcdagMpc^{1/4} eta @@ -156,7 +154,10 @@ public: 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()); gauss = PhiOdd; HighBoundCheck(Mpc,gauss,param.hi); diff --git a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h index e5f0b602..e968b8e4 100644 --- a/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h +++ b/Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h @@ -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 // @@ -170,7 +170,10 @@ NAMESPACE_BEGIN(Grid); msCG_M(MdagM,X,Y); // 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()); gauss = PhiOdd; HighBoundCheck(MdagM,gauss,param.hi); diff --git a/Grid/qcd/action/pseudofermion/OneFlavourRational.h b/Grid/qcd/action/pseudofermion/OneFlavourRational.h index f6c823c9..aa647445 100644 --- a/Grid/qcd/action/pseudofermion/OneFlavourRational.h +++ b/Grid/qcd/action/pseudofermion/OneFlavourRational.h @@ -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} @@ -142,7 +142,10 @@ NAMESPACE_BEGIN(Grid); 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()); gauss = Phi; HighBoundCheck(MdagMOp,gauss,param.hi); diff --git a/Grid/qcd/action/pseudofermion/OneFlavourRationalRatio.h b/Grid/qcd/action/pseudofermion/OneFlavourRationalRatio.h index 5fae2fe9..128c869a 100644 --- a/Grid/qcd/action/pseudofermion/OneFlavourRationalRatio.h +++ b/Grid/qcd/action/pseudofermion/OneFlavourRationalRatio.h @@ -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 // @@ -156,7 +156,10 @@ NAMESPACE_BEGIN(Grid); msCG_M(MdagM,X,Y); // 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()); gauss = Phi; HighBoundCheck(MdagM,gauss,param.hi); diff --git a/Grid/qcd/action/pseudofermion/TwoFlavour.h b/Grid/qcd/action/pseudofermion/TwoFlavour.h index f905a675..2ac97ddd 100644 --- a/Grid/qcd/action/pseudofermion/TwoFlavour.h +++ b/Grid/qcd/action/pseudofermion/TwoFlavour.h @@ -73,7 +73,7 @@ public: ////////////////////////////////////////////////////////////////////////////////////// // 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} // Phi = Mdag eta // P(eta) = e^{- eta^dag eta} diff --git a/Grid/qcd/action/pseudofermion/TwoFlavourEvenOdd.h b/Grid/qcd/action/pseudofermion/TwoFlavourEvenOdd.h index a3cf8f08..2e5208a8 100644 --- a/Grid/qcd/action/pseudofermion/TwoFlavourEvenOdd.h +++ b/Grid/qcd/action/pseudofermion/TwoFlavourEvenOdd.h @@ -77,7 +77,7 @@ public: ////////////////////////////////////////////////////////////////////////////////////// // 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} // Phi = McpDag eta diff --git a/Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h b/Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h index d1d6f336..da628c75 100644 --- a/Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h +++ b/Grid/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h @@ -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} // diff --git a/Grid/qcd/action/pseudofermion/TwoFlavourRatio.h b/Grid/qcd/action/pseudofermion/TwoFlavourRatio.h index 4d72faba..f584706d 100644 --- a/Grid/qcd/action/pseudofermion/TwoFlavourRatio.h +++ b/Grid/qcd/action/pseudofermion/TwoFlavourRatio.h @@ -64,7 +64,7 @@ public: 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} // diff --git a/Grid/qcd/action/scalar/ScalarAction.h b/Grid/qcd/action/scalar/ScalarAction.h index 34fc4fac..8b4f4f79 100644 --- a/Grid/qcd/action/scalar/ScalarAction.h +++ b/Grid/qcd/action/scalar/ScalarAction.h @@ -55,7 +55,7 @@ public: } 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) { return (mass_square * 0.5 + Nd) * ScalarObs::sumphisquared(p) + diff --git a/Grid/qcd/action/scalar/ScalarImpl.h b/Grid/qcd/action/scalar/ScalarImpl.h index 403ea573..13bd6c90 100644 --- a/Grid/qcd/action/scalar/ScalarImpl.h +++ b/Grid/qcd/action/scalar/ScalarImpl.h @@ -27,7 +27,7 @@ public: typedef Field FermionField; 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 gaussian(pRNG, P); P *= scale; @@ -151,7 +151,7 @@ public: 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 #ifndef USE_FFT_ACCELERATION diff --git a/Grid/qcd/action/scalar/ScalarInteractionAction.h b/Grid/qcd/action/scalar/ScalarInteractionAction.h index 5a5f9251..e04dd486 100644 --- a/Grid/qcd/action/scalar/ScalarInteractionAction.h +++ b/Grid/qcd/action/scalar/ScalarInteractionAction.h @@ -77,7 +77,7 @@ public: 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) { diff --git a/Grid/qcd/hmc/HMC.h b/Grid/qcd/hmc/HMC.h index f168b69a..44674ea5 100644 --- a/Grid/qcd/hmc/HMC.h +++ b/Grid/qcd/hmc/HMC.h @@ -139,7 +139,7 @@ private: // Evolution ///////////////////////////////////////////////////////// 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 diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index 70055754..aa28c6c8 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -33,6 +33,7 @@ directory #define INTEGRATOR_INCLUDED #include +#include "MomentumFilter.h" NAMESPACE_BEGIN(Grid); @@ -78,8 +79,19 @@ protected: RepresentationPolicy Representations; IntegratorParameters Params; + //Filters allow the user to manipulate the conjugate momentum, for example to freeze links in DDHMC + //It is applied whenever the momentum is updated / refreshed + //The default filter does nothing + MomentumFilterBase const* MomFilter; + const ActionSet as; + //Get a pointer to a shared static instance of the "do-nothing" momentum filter to serve as a default + static MomentumFilterBase const* getDefaultMomFilter(){ + static MomentumFilterNone filter; + return &filter; + } + void update_P(Field& U, int level, double ep) { t_P[level] += ep; @@ -135,6 +147,8 @@ protected: // Force from the other representations as[level].apply(update_P_hireps, Representations, Mom, U, ep); + + MomFilter->applyFilter(Mom); } void update_U(Field& U, double ep) @@ -174,11 +188,23 @@ public: t_P.resize(levels, 0.0); t_U = 0.0; // initialization of smearer delegated outside of Integrator + + //Default the momentum filter to "do-nothing" + MomFilter = getDefaultMomFilter(); }; virtual ~Integrator() {} virtual std::string integrator_name() = 0; + + //Set the momentum filter allowing for manipulation of the conjugate momentum + void setMomentumFilter(const MomentumFilterBase &filter){ + MomFilter = &filter; + } + + //Access the conjugate momentum + const MomentaField & getMomentum() const{ return P; } + void print_parameters() { @@ -210,10 +236,9 @@ public: // over the representations struct _refresh { template - void operator()(std::vector*> repr_set, Repr& Rep, - GridParallelRNG& pRNG) { + void operator()(std::vector*> repr_set, Repr& Rep, GridSerialRNG & sRNG, GridParallelRNG& pRNG) { 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; } @@ -221,12 +246,12 @@ public: } refresh_hireps{}; // Initialization of momenta and actions - void refresh(Field& U, GridParallelRNG& pRNG) + void refresh(Field& U, GridSerialRNG & sRNG, GridParallelRNG& pRNG) { assert(P.Grid() == U.Grid()); 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 // necessary to keep the fields updated even after a reject @@ -243,12 +268,14 @@ public: // get gauge field from the SmearingPolicy and // based on the boolean is_smeared in actionID 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 - as[level].apply(refresh_hireps, Representations, pRNG); + as[level].apply(refresh_hireps, Representations, sRNG, pRNG); } + + MomFilter->applyFilter(P); } // to be used by the actionlevel class to iterate diff --git a/Grid/qcd/hmc/integrators/MomentumFilter.h b/Grid/qcd/hmc/integrators/MomentumFilter.h new file mode 100644 index 00000000..2a15d80c --- /dev/null +++ b/Grid/qcd/hmc/integrators/MomentumFilter.h @@ -0,0 +1,94 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/hmc/integrators/MomentumFilter.h + +Copyright (C) 2015 + +Author: Christopher Kelly +Author: Peter Boyle + +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. + +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 MOMENTUM_FILTER +#define MOMENTUM_FILTER + +NAMESPACE_BEGIN(Grid); + +//These filter objects allow the user to manipulate the conjugate momentum as part of the update / refresh + +template +struct MomentumFilterBase{ + virtual void applyFilter(MomentaField &P) const; +}; + +//Do nothing +template +struct MomentumFilterNone: public MomentumFilterBase{ + void applyFilter(MomentaField &P) const override{} +}; + +//Multiply each site/direction by a Lorentz vector complex number field +//Can be used to implement a mask, zeroing out sites +template +struct MomentumFilterApplyPhase: public MomentumFilterBase{ + typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type + typedef typename MomentaField::scalar_type scalar_type; //scalar complex type + typedef iVector >, Nd > LorentzScalarType; //complex phase for each site/direction + typedef Lattice LatticeLorentzScalarType; + + LatticeLorentzScalarType phase; + + MomentumFilterApplyPhase(const LatticeLorentzScalarType _phase): phase(_phase){} + + //Default to uniform field of (1,0) + MomentumFilterApplyPhase(GridBase* _grid): phase(_grid){ + LorentzScalarType one; + for(int mu=0;mu #include NAMESPACE_BEGIN(Grid); @@ -40,28 +39,21 @@ public: typedef typename FImpl::FermionField FermionField; typedef typename FImpl::PropagatorField PropagatorField; - typedef typename FImpl::SitePropagator pobj; - typedef typename ComplexField::vector_object vobj; - typedef Lattice> SpinMatrixField; - typedef typename SpinMatrixField::vector_object sobj; - - static const int epsilon[6][3] ; - static const Real epsilon_sgn[6]; private: - template + template accelerator_inline static void BaryonSite(const mobj &D1, - const mobj &D2, - const mobj &D3, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const int parity, - const bool * wick_contractions, - robj &result); - template + const mobj &D2, + const mobj &D3, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const int parity, + const int wick_contractions, + robj &result); + template accelerator_inline static void BaryonSiteMatrix(const mobj &D1, const mobj &D2, const mobj &D3, @@ -69,22 +61,22 @@ public: const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool * wick_contractions, - robj &result); + const int wick_contractions, + robj &result); public: static void WickContractions(std::string qi, std::string qf, - bool* wick_contractions); + int &wick_contractions); static void ContractBaryons(const PropagatorField &q1_left, - const PropagatorField &q2_left, - const PropagatorField &q3_left, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const bool* wick_contractions, - const int parity, - ComplexField &baryon_corr); + const PropagatorField &q2_left, + const PropagatorField &q3_left, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const int wick_contractions, + const int parity, + ComplexField &baryon_corr); static void ContractBaryonsMatrix(const PropagatorField &q1_left, const PropagatorField &q2_left, const PropagatorField &q3_left, @@ -92,20 +84,20 @@ public: const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool* wick_contractions, + const int wick_contractions, SpinMatrixField &baryon_corr); template static void ContractBaryonsSliced(const mobj &D1, - const mobj &D2, - const mobj &D3, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const bool* wick_contractions, - const int parity, - const int nt, - robj &result); + const mobj &D2, + const mobj &D3, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const int wick_contractions, + const int parity, + const int nt, + robj &result); template static void ContractBaryonsSlicedMatrix(const mobj &D1, const mobj &D2, @@ -114,43 +106,43 @@ public: const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const bool* wick_contractions, + const int wick_contractions, const int nt, robj &result); private: - template + template accelerator_inline static void BaryonGamma3ptGroup1Site( const mobj &Dq1_ti, const mobj2 &Dq2_spec, const mobj2 &Dq3_spec, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result); - template + template accelerator_inline static void BaryonGamma3ptGroup2Site( const mobj2 &Dq1_spec, const mobj &Dq2_ti, const mobj2 &Dq3_spec, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result); - template + template accelerator_inline static void BaryonGamma3ptGroup3Site( const mobj2 &Dq1_spec, const mobj2 &Dq2_spec, const mobj &Dq3_ti, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result); public: @@ -162,321 +154,352 @@ public: const PropagatorField &q_tf, int group, int wick_contraction, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, SpinMatrixField &stn_corr); private: - template + template accelerator_inline static void SigmaToNucleonQ1EyeSite(const mobj &Dq_loop, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result); - template + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); + template accelerator_inline static void SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti, - const mobj &Du_tf, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result); + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); - template + template accelerator_inline static void SigmaToNucleonQ2EyeSite(const mobj &Dq_loop, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result); - template + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); + template accelerator_inline static void SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti, - const mobj &Du_tf, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result); + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); + template accelerator_inline + static void XiToSigmaQ1EyeSite(const mobj &Dq_loop, + const mobj2 &Dd_spec, + const mobj2 &Ds_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); + template accelerator_inline + static void XiToSigmaQ2EyeSite(const mobj &Dq_loop, + const mobj2 &Dd_spec, + const mobj2 &Ds_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); public: template static void SigmaToNucleonEye(const PropagatorField &qq_loop, - const mobj &Du_spec, - const PropagatorField &qd_tf, - const PropagatorField &qs_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - const std::string op, - SpinMatrixField &stn_corr); + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr); template static void SigmaToNucleonNonEye(const PropagatorField &qq_ti, - const PropagatorField &qq_tf, - const mobj &Du_spec, - const PropagatorField &qd_tf, - const PropagatorField &qs_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - const std::string op, - SpinMatrixField &stn_corr); + const PropagatorField &qq_tf, + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr); + template + static void XiToSigmaEye(const PropagatorField &qq_loop, + const mobj &Dd_spec, + const mobj &Ds_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &xts_corr); }; - -template -const int BaryonUtils::epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; -/*template -const Complex BaryonUtils::epsilon_sgn[6] = {Complex(1), - Complex(1), - Complex(1), - Complex(-1), - Complex(-1), - Complex(-1)}; -*/ -template -const Real BaryonUtils::epsilon_sgn[6] = {1.,1.,1.,-1.,-1.,-1.}; - -//This is the old version +//This computes a baryon contraction on a lattice site, including the spin-trace of the correlation matrix template -template +template accelerator_inline void BaryonUtils::BaryonSite(const mobj &D1, const mobj &D2, const mobj &D3, - const Gamma GammaA_i, - const Gamma GammaB_i, - const Gamma GammaA_f, - const Gamma GammaB_f, + const Gamma GammaA_i, + const Gamma GammaB_i, + const Gamma GammaA_f, + const Gamma GammaB_f, const int parity, - const bool * wick_contraction, + const int wick_contraction, robj &result) { - Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4) + Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4) - auto D1_GAi = D1 * GammaA_i; - auto D1_GAi_g4 = D1_GAi * g4; - auto D1_GAi_P = 0.5*(D1_GAi + (Real)parity * D1_GAi_g4); - auto GAf_D1_GAi_P = GammaA_f * D1_GAi_P; - auto GBf_D1_GAi_P = GammaB_f * D1_GAi_P; + auto D1_GAi = D1 * GammaA_i; + auto D1_GAi_g4 = D1_GAi * g4; + auto D1_GAi_P = 0.5*(D1_GAi + (Real)parity * D1_GAi_g4); + auto GAf_D1_GAi_P = GammaA_f * D1_GAi_P; + auto GBf_D1_GAi_P = GammaB_f * D1_GAi_P; - auto D2_GBi = D2 * GammaB_i; - auto GBf_D2_GBi = GammaB_f * D2_GBi; - auto GAf_D2_GBi = GammaA_f * D2_GBi; + auto D2_GBi = D2 * GammaB_i; + auto GBf_D2_GBi = GammaB_f * D2_GBi; + auto GAf_D2_GBi = GammaA_f * D2_GBi; - auto GBf_D3 = GammaB_f * D3; - auto GAf_D3 = GammaA_f * D3; + auto GBf_D3 = GammaB_f * D3; + auto GAf_D3 = GammaA_f * D3; - for (int ie_f=0; ie_f < 6 ; ie_f++){ - int a_f = epsilon[ie_f][0]; //a - int b_f = epsilon[ie_f][1]; //b - int c_f = epsilon[ie_f][2]; //c + Real ee; + + for (int ie_f=0; ie_f < 6 ; ie_f++){ + int a_f = (ie_f < 3 ? ie_f : (6-ie_f)%3 ); //epsilon[ie_n][0]; //a + int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b + int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c + int eSgn_f = (ie_f < 3 ? 1 : -1); for (int ie_i=0; ie_i < 6 ; ie_i++){ - int a_i = epsilon[ie_i][0]; //a' - int b_i = epsilon[ie_i][1]; //b' - int c_i = epsilon[ie_i][2]; //c' + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - Real ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; - //This is the \delta_{456}^{123} part - if (wick_contraction[0]){ - for (int rho=0; rho -template +template accelerator_inline void BaryonUtils::BaryonSiteMatrix(const mobj &D1, const mobj &D2, const mobj &D3, - const Gamma GammaA_i, - const Gamma GammaB_i, - const Gamma GammaA_f, - const Gamma GammaB_f, - const bool * wick_contraction, + const Gamma GammaA_i, + const Gamma GammaB_i, + const Gamma GammaA_f, + const Gamma GammaB_f, + const int wick_contraction, robj &result) { - auto D1_GAi = D1 * GammaA_i; - auto GAf_D1_GAi = GammaA_f * D1_GAi; - auto GBf_D1_GAi = GammaB_f * D1_GAi; + auto D1_GAi = D1 * GammaA_i; + auto GAf_D1_GAi = GammaA_f * D1_GAi; + auto GBf_D1_GAi = GammaB_f * D1_GAi; - auto D2_GBi = D2 * GammaB_i; - auto GBf_D2_GBi = GammaB_f * D2_GBi; - auto GAf_D2_GBi = GammaA_f * D2_GBi; + auto D2_GBi = D2 * GammaB_i; + auto GBf_D2_GBi = GammaB_f * D2_GBi; + auto GAf_D2_GBi = GammaA_f * D2_GBi; - auto GBf_D3 = GammaB_f * D3; - auto GAf_D3 = GammaA_f * D3; + auto GBf_D3 = GammaB_f * D3; + auto GAf_D3 = GammaA_f * D3; - for (int ie_f=0; ie_f < 6 ; ie_f++){ - int a_f = epsilon[ie_f][0]; //a - int b_f = epsilon[ie_f][1]; //b - int c_f = epsilon[ie_f][2]; //c + Real ee; + + for (int ie_f=0; ie_f < 6 ; ie_f++){ + int a_f = (ie_f < 3 ? ie_f : (6-ie_f)%3 ); //epsilon[ie_n][0]; //a + int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b + int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c + int eSgn_f = (ie_f < 3 ? 1 : -1); for (int ie_i=0; ie_i < 6 ; ie_i++){ - int a_i = epsilon[ie_i][0]; //a' - int b_i = epsilon[ie_i][1]; //b' - int c_i = epsilon[ie_i][2]; //c' + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - Real ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; - //This is the \delta_{456}^{123} part - if (wick_contraction[0]){ - for (int rho_i=0; rho_i::BaryonSiteMatrix(const mobj &D1, * flavours. * * The array wick_contractions must be of length 6 */ template -void BaryonUtils::WickContractions(std::string qi, std::string qf, bool* wick_contractions) { +void BaryonUtils::WickContractions(std::string qi, std::string qf, int &wick_contractions) { + assert(qi.size() == 3 && qf.size() == 3 && "Only sets of 3 quarks accepted."); const int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; + wick_contractions=0; for (int ie=0; ie < 6 ; ie++) { - wick_contractions[ie] = (qi.size() == 3 && qf.size() == 3 - && qi[0] == qf[epsilon[ie][0]] + wick_contractions += ( ( qi[0] == qf[epsilon[ie][0]] && qi[1] == qf[epsilon[ie][1]] - && qi[2] == qf[epsilon[ie][2]]); + && qi[2] == qf[epsilon[ie][2]]) ? 1 : 0) << ie; } } @@ -501,15 +525,15 @@ void BaryonUtils::WickContractions(std::string qi, std::string qf, bool* * Wick_Contractions function above */ template void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, - const PropagatorField &q2_left, - const PropagatorField &q3_left, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const bool* wick_contractions, - const int parity, - ComplexField &baryon_corr) + const PropagatorField &q2_left, + const PropagatorField &q3_left, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const int wick_contractions, + const int parity, + ComplexField &baryon_corr) { assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); @@ -519,31 +543,31 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, GridBase *grid = q1_left.Grid(); - autoView(vbaryon_corr, baryon_corr,CpuWrite); - autoView( v1 , q1_left, CpuRead); - autoView( v2 , q2_left, CpuRead); - autoView( v3 , q3_left, CpuRead); + autoView(vbaryon_corr , baryon_corr , AcceleratorWrite); + autoView( v1 , q1_left , AcceleratorRead); + autoView( v2 , q2_left , AcceleratorRead); + autoView( v3 , q3_left , AcceleratorRead); Real bytes =0.; bytes += grid->oSites() * (432.*sizeof(vComplex) + 126.*sizeof(int) + 36.*sizeof(Real)); for (int ie=0; ie < 6 ; ie++){ if(ie==0 or ie==3){ - bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contractions[ie]; - } - else{ - bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contractions[ie]; + bytes += ( wick_contractions & (1 << ie) ) ? grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) : 0.; + } else{ + bytes += ( wick_contractions & (1 << ie) ) ? grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) : 0.; } } Real t=0.; t =-usecond(); accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto D1 = v1[ss]; - auto D2 = v2[ss]; - auto D3 = v3[ss]; - vobj result=Zero(); + auto D1 = v1(ss); + auto D2 = v2(ss); + auto D3 = v3(ss); + typedef decltype(coalescedRead(vbaryon_corr[0])) cVec; + cVec result=Zero(); BaryonSite(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contractions,result); - vbaryon_corr[ss] = result; + coalescedWrite(vbaryon_corr[ss],result); } );//end loop over lattice sites t += usecond(); @@ -555,50 +579,33 @@ template void BaryonUtils::ContractBaryonsMatrix(const PropagatorField &q1_left, const PropagatorField &q2_left, const PropagatorField &q3_left, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const bool* wick_contractions, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const int wick_contractions, SpinMatrixField &baryon_corr) { assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); assert(Nc==3 && "Baryon code only implemented for N_colour = 3"); - - GridBase *grid = q1_left.Grid(); - - autoView(vbaryon_corr, baryon_corr,CpuWrite); - autoView( v1 , q1_left, CpuRead); - autoView( v2 , q2_left, CpuRead); - autoView( v3 , q3_left, CpuRead); - // Real bytes =0.; - // bytes += grid->oSites() * (432.*sizeof(vComplex) + 126.*sizeof(int) + 36.*sizeof(Real)); - // for (int ie=0; ie < 6 ; ie++){ - // if(ie==0 or ie==3){ - // bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contractions[ie]; - // } - // else{ - // bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contractions[ie]; - // } - // } - // Real t=0.; - // t =-usecond(); + GridBase *grid = q1_left.Grid(); + + autoView(vbaryon_corr , baryon_corr , AcceleratorWrite); + autoView( v1 , q1_left , AcceleratorRead); + autoView( v2 , q2_left , AcceleratorRead); + autoView( v3 , q3_left , AcceleratorRead); accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto D1 = v1[ss]; - auto D2 = v2[ss]; - auto D3 = v3[ss]; - sobj result=Zero(); + auto D1 = v1(ss); + auto D2 = v2(ss); + auto D3 = v3(ss); + typedef decltype(coalescedRead(vbaryon_corr[0])) spinor; + spinor result=Zero(); BaryonSiteMatrix(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,wick_contractions,result); - vbaryon_corr[ss] = result; + coalescedWrite(vbaryon_corr[ss],result); } );//end loop over lattice sites - - // t += usecond(); - - // std::cout << GridLogDebug << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s " << std::endl; - } /* The array wick_contractions must be of length 6. The order * @@ -609,16 +616,16 @@ void BaryonUtils::ContractBaryonsMatrix(const PropagatorField &q1_left, template template void BaryonUtils::ContractBaryonsSliced(const mobj &D1, - const mobj &D2, - const mobj &D3, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const bool* wick_contractions, - const int parity, - const int nt, - robj &result) + const mobj &D2, + const mobj &D3, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const int wick_contractions, + const int parity, + const int nt, + robj &result) { assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); @@ -636,11 +643,11 @@ template void BaryonUtils::ContractBaryonsSlicedMatrix(const mobj &D1, const mobj &D2, const mobj &D3, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, - const bool* wick_contractions, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const int wick_contractions, const int nt, robj &result) { @@ -664,91 +671,93 @@ void BaryonUtils::ContractBaryonsSlicedMatrix(const mobj &D1, * Dq3_spec is a quark line from t_i to t_f * Dq4_tf is a quark line from t_f to t_J */ template -template +template accelerator_inline void BaryonUtils::BaryonGamma3ptGroup1Site( const mobj &Dq1_ti, const mobj2 &Dq2_spec, const mobj2 &Dq3_spec, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result) { - Gamma g5(Gamma::Algebra::Gamma5); + Gamma g5(Gamma::Algebra::Gamma5); - auto adjD4_g_D1 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq1_ti; - auto Gf_adjD4_g_D1 = GammaBf * adjD4_g_D1; - auto D2_Gi = Dq2_spec * GammaBi; - auto Gf_D2_Gi = GammaBf * D2_Gi; - auto Gf_D3 = GammaBf * Dq3_spec; + auto adjD4 = g5 * adj(Dq4_tf) * g5 ; + auto adjD4_g_D1 = adjD4 * GammaJ * Dq1_ti; + auto Gf_adjD4_g_D1 = GammaBf * adjD4_g_D1; + auto D2_Gi = Dq2_spec * GammaBi; + auto Gf_D2_Gi = GammaBf * D2_Gi; + auto Gf_D3 = GammaBf * Dq3_spec; - int a_f, b_f, c_f; - int a_i, b_i, c_i; + Real ee; - Real ee; - - for (int ie_f=0; ie_f < 6 ; ie_f++){ - a_f = epsilon[ie_f][0]; //a - b_f = epsilon[ie_f][1]; //b - c_f = epsilon[ie_f][2]; //c + for (int ie_f=0; ie_f < 6 ; ie_f++){ + int a_f = (ie_f < 3 ? ie_f : (6-ie_f)%3 ); //epsilon[ie_n][0]; //a + int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b + int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c + int eSgn_f = (ie_f < 3 ? 1 : -1); for (int ie_i=0; ie_i < 6 ; ie_i++){ - a_i = epsilon[ie_i][0]; //a' - b_i = epsilon[ie_i][1]; //b' - c_i = epsilon[ie_i][2]; //c' + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + ee = Real(eSgn_f * eSgn_i); - for (int alpha_f=0; alpha_f::BaryonGamma3ptGroup1Site( * Dq3_spec is a quark line from t_i to t_f * Dq4_tf is a quark line from t_f to t_J */ template -template +template accelerator_inline void BaryonUtils::BaryonGamma3ptGroup2Site( const mobj2 &Dq1_spec, const mobj &Dq2_ti, const mobj2 &Dq3_spec, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result) { - Gamma g5(Gamma::Algebra::Gamma5); + Gamma g5(Gamma::Algebra::Gamma5); - auto adjD4_g_D2_Gi = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq2_ti * GammaBi; - auto Gf_adjD4_g_D2_Gi = GammaBf * adjD4_g_D2_Gi; - auto Gf_D1 = GammaBf * Dq1_spec; - auto Gf_D3 = GammaBf * Dq3_spec; + auto adjD4_g_D2_Gi = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq2_ti * GammaBi; + auto Gf_adjD4_g_D2_Gi = GammaBf * adjD4_g_D2_Gi; + auto Gf_D1 = GammaBf * Dq1_spec; + auto Gf_D3 = GammaBf * Dq3_spec; - int a_f, b_f, c_f; - int a_i, b_i, c_i; + Real ee; - Real ee; - - for (int ie_f=0; ie_f < 6 ; ie_f++){ - a_f = epsilon[ie_f][0]; //a - b_f = epsilon[ie_f][1]; //b - c_f = epsilon[ie_f][2]; //c + for (int ie_f=0; ie_f < 6 ; ie_f++){ + int a_f = (ie_f < 3 ? ie_f : (6-ie_f)%3 ); //epsilon[ie_n][0]; //a + int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b + int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c + int eSgn_f = (ie_f < 3 ? 1 : -1); for (int ie_i=0; ie_i < 6 ; ie_i++){ - a_i = epsilon[ie_i][0]; //a' - b_i = epsilon[ie_i][1]; //b' - c_i = epsilon[ie_i][2]; //c' + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; - for (int alpha_f=0; alpha_f::BaryonGamma3ptGroup2Site( * Dq3_ti is a quark line from t_i to t_J * Dq4_tf is a quark line from t_f to t_J */ template -template +template accelerator_inline void BaryonUtils::BaryonGamma3ptGroup3Site( const mobj2 &Dq1_spec, const mobj2 &Dq2_spec, const mobj &Dq3_ti, const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, int wick_contraction, robj &result) { - Gamma g5(Gamma::Algebra::Gamma5); + Gamma g5(Gamma::Algebra::Gamma5); - auto adjD4_g_D3 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq3_ti; - auto Gf_adjD4_g_D3 = GammaBf * adjD4_g_D3; - auto Gf_D1 = GammaBf * Dq1_spec; - auto D2_Gi = Dq2_spec * GammaBi; - auto Gf_D2_Gi = GammaBf * D2_Gi; + auto adjD4_g_D3 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq3_ti; + auto Gf_adjD4_g_D3 = GammaBf * adjD4_g_D3; + auto Gf_D1 = GammaBf * Dq1_spec; + auto D2_Gi = Dq2_spec * GammaBi; + auto Gf_D2_Gi = GammaBf * D2_Gi; - int a_f, b_f, c_f; - int a_i, b_i, c_i; + Real ee; - Real ee; - - for (int ie_f=0; ie_f < 6 ; ie_f++){ - a_f = epsilon[ie_f][0]; //a - b_f = epsilon[ie_f][1]; //b - c_f = epsilon[ie_f][2]; //c + for (int ie_f=0; ie_f < 6 ; ie_f++){ + int a_f = (ie_f < 3 ? ie_f : (6-ie_f)%3 ); //epsilon[ie_n][0]; //a + int b_f = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b + int c_f = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c + int eSgn_f = (ie_f < 3 ? 1 : -1); for (int ie_i=0; ie_i < 6 ; ie_i++){ - a_i = epsilon[ie_i][0]; //a' - b_i = epsilon[ie_i][1]; //b' - c_i = epsilon[ie_i][2]; //c' + int a_i = (ie_i < 3 ? ie_i : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a' + int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b' + int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_i = (ie_i < 3 ? 1 : -1); - ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; - for (int alpha_f=0; alpha_f::BaryonGamma3pt( const PropagatorField &q_tf, int group, int wick_contraction, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, SpinMatrixField &stn_corr) { - GridBase *grid = q_tf.Grid(); + assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); + assert(Nc==3 && "Baryon code only implemented for N_colour = 3"); - autoView( vcorr, stn_corr, CpuWrite); - autoView( vq_ti , q_ti, CpuRead); - autoView( vq_tf , q_tf, CpuRead); + GridBase *grid = q_tf.Grid(); - if (group == 1) { - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti[ss]; - auto Dq_tf = vq_tf[ss]; - sobj result=Zero(); - BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - vcorr[ss] += result; - });//end loop over lattice sites - } else if (group == 2) { - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti[ss]; - auto Dq_tf = vq_tf[ss]; - sobj result=Zero(); - BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - vcorr[ss] += result; - });//end loop over lattice sites - } else if (group == 3) { - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti[ss]; - auto Dq_tf = vq_tf[ss]; - sobj result=Zero(); - BaryonGamma3ptGroup3Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + autoView( vcorr , stn_corr , AcceleratorWrite); + autoView( vq_ti , q_ti , AcceleratorRead); + autoView( vq_tf , q_tf , AcceleratorRead); + + Vector my_Dq_spec{Dq_spec1,Dq_spec2}; + mobj * Dq_spec_p = &my_Dq_spec[0]; + + if (group == 1) { + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec_p[0],Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); + });//end loop over lattice sites + } else if (group == 2) { + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + BaryonGamma3ptGroup2Site(Dq_spec_p[0],Dq_ti,Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); + });//end loop over lattice sites + } else if (group == 3) { + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + BaryonGamma3ptGroup3Site(Dq_spec_p[0],Dq_spec_p[1],Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + coalescedWrite(vcorr[ss],coalescedRead(vcorr[ss])+result); + });//end loop over lattice sites + } - vcorr[ss] += result; - });//end loop over lattice sites - } } /*********************************************************************** * End of BaryonGamma3pt-function code. * - * * + * * * The following code is for Sigma -> N rare hypeon decays * **********************************************************************/ @@ -997,46 +1017,56 @@ void BaryonUtils::BaryonGamma3pt( * Dd_tf is a quark line from t_f to t_H * Ds_ti is a quark line from t_i to t_H */ template -template +template accelerator_inline void BaryonUtils::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result) + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) { Gamma g5(Gamma::Algebra::Gamma5); - auto DuG = Du_spec * GammaB_nucl; - // Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) - auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5; - // Dq_loop * \gamma_\mu^L - auto DqG = Dq_loop * Gamma_H; + auto adjDd_GH_Ds = g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti; + auto Gn_adjDd_GH_Ds = GammaB_nucl * adjDd_GH_Ds; + auto Du_Gs = Du_spec * GammaB_sigma; + auto Dq_GH = Dq_loop * Gamma_H; + auto Tr_Dq_GH = trace(Dq_GH)()()(); + + Real ee; for (int ie_n=0; ie_n < 6 ; ie_n++){ - int a_n = epsilon[ie_n][0]; //a - int b_n = epsilon[ie_n][1]; //b - int c_n = epsilon[ie_n][2]; //c + int a_n = (ie_n < 3 ? ie_n : (6-ie_n)%3 ); //epsilon[ie_n][0]; //a + int b_n = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b + int c_n = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c + int eSgn_n = (ie_n < 3 ? 1 : -1); for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = epsilon[ie_s][0]; //a' - int b_s = epsilon[ie_s][1]; //b' - int c_s = epsilon[ie_s][2]; //c' - for (int alpha_s=0; alpha_s::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop, * Dd_tf is a quark line from t_f to t_H * Ds_ti is a quark line from t_i to t_H */ template -template +template accelerator_inline void BaryonUtils::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti, - const mobj &Du_tf, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result) + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) { Gamma g5(Gamma::Algebra::Gamma5); - auto DuG = Du_spec * GammaB_nucl; - auto adjDu = g5 * adj(Du_tf) * g5; - auto adjDuG = adjDu * GammaB_nucl; - // Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) - auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5; - // Dq_loop * \gamma_\mu^L - auto DuGH = Du_ti * Gamma_H; + auto Du_Gs = Du_spec * GammaB_sigma; + auto adjDd_GH_Ds = g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti; + auto Gn_adjDd_GH_Ds = GammaB_nucl * adjDd_GH_Ds; + auto adjDu_GH_Du = g5 * adj(Du_tf) * g5 * Gamma_H * Du_ti; + auto adjDu_GH_Du_Gs = adjDu_GH_Du * GammaB_sigma; + + Real ee; for (int ie_n=0; ie_n < 6 ; ie_n++){ - int a_n = epsilon[ie_n][0]; //a - int b_n = epsilon[ie_n][1]; //b - int c_n = epsilon[ie_n][2]; //c + int a_n = (ie_n < 3 ? ie_n : (6-ie_n)%3 ); //epsilon[ie_n][0]; //a + int b_n = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b + int c_n = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c + int eSgn_n = (ie_n < 3 ? 1 : -1); for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = epsilon[ie_s][0]; //a' - int b_s = epsilon[ie_s][1]; //b' - int c_s = epsilon[ie_s][2]; //c' - for (int alpha_s=0; alpha_s::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti, * Dd_tf is a quark line from t_f to t_H * Ds_ti is a quark line from t_i to t_H */ template -template +template accelerator_inline void BaryonUtils::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop, - const mobj2 &Du_spec, - const mobj &Dd_tf, - const mobj &Ds_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - robj &result) + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) { Gamma g5(Gamma::Algebra::Gamma5); - auto DuG = Du_spec * GammaB_nucl; - // Gamma^B * Ds * \gamma_\mu^L - auto GDsG = GammaB_sigma * Ds_ti * Gamma_H; - // Dq_loop * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) - auto DqGDd = Dq_loop * Gamma_H * g5 * adj(Dd_tf) * g5; + auto adjDd_GH_Duloop_GH_Ds = g5 * adj(Dd_tf) * g5 * Gamma_H * Dq_loop * Gamma_H * Ds_ti; + auto Gn_adjDd_GH_Duloop_GH_Ds = GammaB_nucl * adjDd_GH_Duloop_GH_Ds; + auto Du_Gs = Du_spec * GammaB_sigma; + + Real ee; for (int ie_n=0; ie_n < 6 ; ie_n++){ - int a_n = epsilon[ie_n][0]; //a - int b_n = epsilon[ie_n][1]; //b - int c_n = epsilon[ie_n][2]; //c + int a_n = (ie_n < 3 ? ie_n : (6-ie_n)%3 ); //epsilon[ie_n][0]; //a + int b_n = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b + int c_n = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c + int eSgn_n = (ie_n < 3 ? 1 : -1); for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = epsilon[ie_s][0]; //a' - int b_s = epsilon[ie_s][1]; //b' - int c_s = epsilon[ie_s][2]; //c' - for (int alpha_s=0; alpha_s::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop, * Dd_tf is a quark line from t_f to t_H * Ds_ti is a quark line from t_i to t_H */ template -template +template accelerator_inline void BaryonUtils::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti, - const mobj &Du_tf, - const mobj2 &Du_spec, + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) +{ + + Gamma g5(Gamma::Algebra::Gamma5); + + auto Du_Gs = Du_spec * GammaB_sigma; + auto adjDu_GH_Ds = g5 * adj(Du_tf) * g5 * Gamma_H * Ds_ti; + auto adjDd_GH_Du = g5 * adj(Dd_tf) * g5 * Gamma_H * Du_ti; + auto Gn_adjDd_GH_Du = GammaB_nucl * adjDd_GH_Du; // for some reason I needed to split this into two lines to avoid the compilation error 'error: identifier "Grid::Gamma::mul" is undefined in device code' + + auto Gn_adjDd_GH_Du_Gs = Gn_adjDd_GH_Du * GammaB_sigma; + + Real ee; + + for (int ie_n=0; ie_n < 6 ; ie_n++){ + int a_n = (ie_n < 3 ? ie_n : (6-ie_n)%3 ); //epsilon[ie_n][0]; //a + int b_n = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b + int c_n = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c + int eSgn_n = (ie_n < 3 ? 1 : -1); + for (int ie_s=0; ie_s < 6 ; ie_s++){ + int a_s = (ie_s < 3 ? ie_s : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a' + int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b' + int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_s = (ie_s < 3 ? 1 : -1); + + ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s]; + + for (int alpha_n=0; alpha_n +template +void BaryonUtils::SigmaToNucleonEye(const PropagatorField &qq_loop, + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr) +{ + + assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); + assert(Nc==3 && "Baryon code only implemented for N_colour = 3"); + + GridBase *grid = qs_ti.Grid(); + + autoView( vcorr , stn_corr , AcceleratorWrite); + autoView( vq_loop , qq_loop , AcceleratorRead); + autoView( vd_tf , qd_tf , AcceleratorRead); + autoView( vs_ti , qs_ti , AcceleratorRead); + + Vector my_Dq_spec{Du_spec}; + mobj * Dq_spec_p = &my_Dq_spec[0]; + + if(op == "Q1"){ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_loop = vq_loop(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + SigmaToNucleonQ1EyeSite(Dq_loop,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites + } else if(op == "Q2"){ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_loop = vq_loop(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + SigmaToNucleonQ2EyeSite(Dq_loop,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites + } else { + assert(0 && "Weak Operator not correctly specified"); + } +} + +template +template +void BaryonUtils::SigmaToNucleonNonEye(const PropagatorField &qq_ti, + const PropagatorField &qq_tf, + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr) +{ + + assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); + assert(Nc==3 && "Baryon code only implemented for N_colour = 3"); + + GridBase *grid = qs_ti.Grid(); + + autoView( vcorr , stn_corr , AcceleratorWrite ); + autoView( vq_ti , qq_ti , AcceleratorRead ); + autoView( vq_tf , qq_tf , AcceleratorRead ); + autoView( vd_tf , qd_tf , AcceleratorRead ); + autoView( vs_ti , qs_ti , AcceleratorRead ); + + Vector my_Dq_spec{Du_spec}; + mobj * Dq_spec_p = &my_Dq_spec[0]; + + if(op == "Q1"){ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + SigmaToNucleonQ1NonEyeSite(Dq_ti,Dq_tf,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites + } else if(op == "Q2"){ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti(ss); + auto Dq_tf = vq_tf(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + SigmaToNucleonQ2NonEyeSite(Dq_ti,Dq_tf,Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + coalescedWrite(vcorr[ss],result); + });//end loop over lattice sites + } else { + assert(0 && "Weak Operator not correctly specified"); + } +} + + +/*********************************************************************** + * The following code is for Xi -> Sigma rare hypeon decays * + **********************************************************************/ + +/* Dq_loop is a quark line from t_H to t_H + * Dd_spec is a quark line from t_i to t_f + * Ds_spec is a quark line from t_i to t_f + * Dd_tf is a quark line from t_f to t_H + * Ds_ti is a quark line from t_i to t_H */ +template +template accelerator_inline +void BaryonUtils::XiToSigmaQ1EyeSite(const mobj &Dq_loop, + const mobj2 &Dd_spec, + const mobj2 &Ds_spec, const mobj &Dd_tf, const mobj &Ds_ti, const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma GammaB_xi, + const Gamma GammaB_sigma, robj &result) { Gamma g5(Gamma::Algebra::Gamma5); - auto DuG = Du_spec * GammaB_nucl; - auto adjDu = g5 * adj(Du_tf) * g5; - auto adjDuG = adjDu * GammaB_nucl; - // Gamma^B * Ds * \gamma_\mu^L - auto GDsG = GammaB_sigma * Ds_ti * Gamma_H; - // Du * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) - auto DuGDd = Du_ti * Gamma_H * g5 * adj(Dd_tf) * g5; + auto DdG = Dd_spec * GammaB_sigma; + auto GDs = GammaB_xi * Ds_spec; + // Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) + auto DsGDd = Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5; + // DsGDd * GammaB + auto DsGDdG = DsGDd * GammaB_sigma; + // GammaB * DsGDd + auto GDsGDd = GammaB_xi * DsGDd; + // GammaB * DsGDd * GammaB + auto GDsGDdG = GDsGDd * GammaB_sigma; + // \gamma_\mu^L * Dq_loop + auto trGDq = TensorRemove(trace(Gamma_H * Dq_loop)); - for (int ie_n=0; ie_n < 6 ; ie_n++){ - int a_n = epsilon[ie_n][0]; //a - int b_n = epsilon[ie_n][1]; //b - int c_n = epsilon[ie_n][2]; //c - for (int ie_s=0; ie_s < 6 ; ie_s++){ - int a_s = epsilon[ie_s][0]; //a' - int b_s = epsilon[ie_s][1]; //b' - int c_s = epsilon[ie_s][2]; //c' - for (int alpha_s=0; alpha_s +template accelerator_inline +void BaryonUtils::XiToSigmaQ2EyeSite(const mobj &Dq_loop, + const mobj2 &Dd_spec, + const mobj2 &Ds_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_xi, + const Gamma GammaB_sigma, + robj &result) +{ + + Gamma g5(Gamma::Algebra::Gamma5); + + auto DdG = Dd_spec * GammaB_sigma; + auto GDs = GammaB_xi * Ds_spec; + // Ds * \gamma_\mu^L * Dq_loop * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) + auto DsGDqGDd = Ds_ti * Gamma_H * Dq_loop * Gamma_H * g5 * adj(Dd_tf) * g5; + // DsGDd * GammaB + auto DsGDqGDdG = DsGDqGDd * GammaB_sigma; + // GammaB * DsGDd + auto GDsGDqGDd = GammaB_xi * DsGDqGDd; + // GammaB * DsGDd * GammaB + auto GDsGDqGDdG = GDsGDqGDd * GammaB_sigma; + + Real ee; + + for (int ie_s=0; ie_s < 6 ; ie_s++){ + int a_s = (ie_s < 3 ? ie_s : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a' + int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b' + int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c' + int eSgn_s = (ie_s < 3 ? 1 : -1); + for (int ie_x=0; ie_x < 6 ; ie_x++){ + int a_x = (ie_x < 3 ? ie_x : (6-ie_x)%3 ); //epsilon[ie_x][0]; //a' + int b_x = (ie_x < 3 ? (ie_x+1)%3 : (8-ie_x)%3 ); //epsilon[ie_x][1]; //b' + int c_x = (ie_x < 3 ? (ie_x+2)%3 : (7-ie_x)%3 ); //epsilon[ie_x][2]; //c' + int eSgn_x = (ie_x < 3 ? 1 : -1); + ee = Real(eSgn_s * eSgn_x); + for (int alpha_x=0; alpha_x template -void BaryonUtils::SigmaToNucleonEye(const PropagatorField &qq_loop, - const mobj &Du_spec, +void BaryonUtils::XiToSigmaEye(const PropagatorField &qq_loop, + const mobj &Dd_spec, + const mobj &Ds_spec, const PropagatorField &qd_tf, const PropagatorField &qs_ti, const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, + const Gamma GammaB_xi, + const Gamma GammaB_sigma, const std::string op, - SpinMatrixField &stn_corr) + SpinMatrixField &xts_corr) { assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); @@ -1229,67 +1539,38 @@ void BaryonUtils::SigmaToNucleonEye(const PropagatorField &qq_loop, GridBase *grid = qs_ti.Grid(); - autoView( vcorr, stn_corr, CpuWrite); - autoView( vq_loop , qq_loop, CpuRead); - autoView( vd_tf , qd_tf, CpuRead); - autoView( vs_ti , qs_ti, CpuRead); + autoView( vcorr , xts_corr , AcceleratorWrite); + autoView( vq_loop , qq_loop , AcceleratorRead); + autoView( vd_tf , qd_tf , AcceleratorRead); + autoView( vs_ti , qs_ti , AcceleratorRead); - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_loop = vq_loop[ss]; - auto Dd_tf = vd_tf[ss]; - auto Ds_ti = vs_ti[ss]; - sobj result=Zero(); - if(op == "Q1"){ - SigmaToNucleonQ1EyeSite(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); - } else if(op == "Q2"){ - SigmaToNucleonQ2EyeSite(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); - } else { - assert(0 && "Weak Operator not correctly specified"); - } - vcorr[ss] = result; - } );//end loop over lattice sites + Vector my_Dq_spec{Dd_spec,Ds_spec}; + mobj * Dq_spec_p = &my_Dq_spec[0]; + + if(op == "Q1"){ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_loop = vq_loop(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + XiToSigmaQ1EyeSite(Dq_loop,Dq_spec_p[0],Dq_spec_p[1],Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result); + coalescedWrite(vcorr[ss],result); + } );//end loop over lattice sites + } else if(op == "Q2"){ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_loop = vq_loop(ss); + auto Dd_tf = vd_tf(ss); + auto Ds_ti = vs_ti(ss); + typedef decltype(coalescedRead(vcorr[0])) spinor; + spinor result=Zero(); + XiToSigmaQ2EyeSite(Dq_loop,Dq_spec_p[0],Dq_spec_p[1],Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result); + coalescedWrite(vcorr[ss],result); + } );//end loop over lattice sites + } else { + assert(0 && "Weak Operator not correctly specified"); + } } -template -template -void BaryonUtils::SigmaToNucleonNonEye(const PropagatorField &qq_ti, - const PropagatorField &qq_tf, - const mobj &Du_spec, - const PropagatorField &qd_tf, - const PropagatorField &qs_ti, - const Gamma Gamma_H, - const Gamma GammaB_sigma, - const Gamma GammaB_nucl, - const std::string op, - SpinMatrixField &stn_corr) -{ - - assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); - assert(Nc==3 && "Baryon code only implemented for N_colour = 3"); - - GridBase *grid = qs_ti.Grid(); - - autoView( vcorr , stn_corr, CpuWrite); - autoView( vq_ti , qq_ti, CpuRead); - autoView( vq_tf , qq_tf, CpuRead); - autoView( vd_tf , qd_tf, CpuRead); - autoView( vs_ti , qs_ti, CpuRead); - // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - thread_for(ss,grid->oSites(),{ - auto Dq_ti = vq_ti[ss]; - auto Dq_tf = vq_tf[ss]; - auto Dd_tf = vd_tf[ss]; - auto Ds_ti = vs_ti[ss]; - sobj result=Zero(); - if(op == "Q1"){ - SigmaToNucleonQ1NonEyeSite(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); - } else if(op == "Q2"){ - SigmaToNucleonQ2NonEyeSite(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); - } else { - assert(0 && "Weak Operator not correctly specified"); - } - vcorr[ss] = result; - } );//end loop over lattice sites -} NAMESPACE_END(Grid); diff --git a/Grid/qcd/utils/Metric.h b/Grid/qcd/utils/Metric.h index 10d06de8..d8ae27dc 100644 --- a/Grid/qcd/utils/Metric.h +++ b/Grid/qcd/utils/Metric.h @@ -93,13 +93,13 @@ public: GeneralisedMomenta(GridBase* grid, Metric& M): M(M), Mom(grid), AuxMom(grid), AuxField(grid){} // Correct - void MomentaDistribution(GridParallelRNG& pRNG){ + void MomentaDistribution(GridSerialRNG & sRNG, GridParallelRNG& pRNG){ // Generate a distribution for // P^dag G P // where G = M^-1 // Generate gaussian momenta - Implementation::generate_momenta(Mom, pRNG); + Implementation::generate_momenta(Mom, sRNG, pRNG); // Modify the distribution with the metric M.MSquareRoot(Mom); @@ -107,8 +107,8 @@ public: // Auxiliary momenta // do nothing if trivial, so hide in the metric MomentaField AuxMomTemp(Mom.Grid()); - Implementation::generate_momenta(AuxMom, pRNG); - Implementation::generate_momenta(AuxField, pRNG); + Implementation::generate_momenta(AuxMom, sRNG, pRNG); + Implementation::generate_momenta(AuxField, sRNG, pRNG); // Modify the distribution with the metric // Aux^dag M Aux M.MInvSquareRoot(AuxMom); // AuxMom = M^{-1/2} AuxMomTemp diff --git a/Grid/simd/Grid_gpu_vec.h b/Grid/simd/Grid_gpu_vec.h index 8e55ce2f..b2c7588f 100644 --- a/Grid/simd/Grid_gpu_vec.h +++ b/Grid/simd/Grid_gpu_vec.h @@ -60,11 +60,26 @@ template class GpuComplex { public: pair z; - typedef decltype(z.x) real; + typedef decltype(z.x) Real; public: accelerator_inline GpuComplex() = default; - accelerator_inline GpuComplex(real re,real im) { z.x=re; z.y=im; }; + accelerator_inline GpuComplex(Real re,Real im) { z.x=re; z.y=im; }; accelerator_inline GpuComplex(const GpuComplex &zz) { z = zz.z;}; + accelerator_inline Real real(void) const { return z.x; }; + 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) { + *this = (*this) * r; + return *this; + } + accelerator_inline GpuComplex &operator+=(const GpuComplex &r) { + *this = (*this) + r; + return *this; + } + accelerator_inline GpuComplex &operator-=(const GpuComplex &r) { + *this = (*this) - r; + return *this; + } friend accelerator_inline GpuComplex operator+(const GpuComplex &lhs,const GpuComplex &rhs) { GpuComplex r ; r.z.x = lhs.z.x + rhs.z.x; @@ -157,6 +172,11 @@ typedef GpuVector GpuVectorRD; typedef GpuVector GpuVectorCD; typedef GpuVector GpuVectorI; +accelerator_inline GpuComplexF timesI(const GpuComplexF &r) { return(GpuComplexF(-r.imag(),r.real()));} +accelerator_inline GpuComplexD timesI(const GpuComplexD &r) { return(GpuComplexD(-r.imag(),r.real()));} +accelerator_inline GpuComplexF timesMinusI(const GpuComplexF &r){ return(GpuComplexF(r.imag(),-r.real()));} +accelerator_inline GpuComplexD timesMinusI(const GpuComplexD &r){ return(GpuComplexD(r.imag(),-r.real()));} + accelerator_inline float half2float(half h) { float f; diff --git a/Grid/simd/Grid_vector_types.h b/Grid/simd/Grid_vector_types.h index c07077a3..4f952bb2 100644 --- a/Grid/simd/Grid_vector_types.h +++ b/Grid/simd/Grid_vector_types.h @@ -208,8 +208,8 @@ struct RealPart > { ////////////////////////////////////// // type alias used to simplify the syntax of std::enable_if template using Invoke = typename T::type; -template using EnableIf = Invoke >; -template using NotEnableIf = Invoke >; +template using EnableIf = Invoke >; +template using NotEnableIf = Invoke >; //////////////////////////////////////////////////////// // Check for complexity with type traits diff --git a/Grid/simd/Simd.h b/Grid/simd/Simd.h index 1dc86c1b..76ca3bef 100644 --- a/Grid/simd/Simd.h +++ b/Grid/simd/Simd.h @@ -148,10 +148,14 @@ accelerator_inline void sub (ComplexF * __restrict__ y,const ComplexF * __restri accelerator_inline void add (ComplexF * __restrict__ y,const ComplexF * __restrict__ l,const ComplexF *__restrict__ r){ *y = (*l) + (*r); } //conjugate already supported for complex -accelerator_inline ComplexF timesI(const ComplexF &r) { return(r*ComplexF(0.0,1.0));} -accelerator_inline ComplexD timesI(const ComplexD &r) { return(r*ComplexD(0.0,1.0));} -accelerator_inline ComplexF timesMinusI(const ComplexF &r){ return(r*ComplexF(0.0,-1.0));} -accelerator_inline ComplexD timesMinusI(const ComplexD &r){ return(r*ComplexD(0.0,-1.0));} +accelerator_inline ComplexF timesI(const ComplexF &r) { return(ComplexF(-r.imag(),r.real()));} +accelerator_inline ComplexD timesI(const ComplexD &r) { return(ComplexD(-r.imag(),r.real()));} +accelerator_inline ComplexF timesMinusI(const ComplexF &r){ return(ComplexF(r.imag(),-r.real()));} +accelerator_inline ComplexD timesMinusI(const ComplexD &r){ return(ComplexD(r.imag(),-r.real()));} +//accelerator_inline ComplexF timesI(const ComplexF &r) { return(r*ComplexF(0.0,1.0));} +//accelerator_inline ComplexD timesI(const ComplexD &r) { return(r*ComplexD(0.0,1.0));} +//accelerator_inline ComplexF timesMinusI(const ComplexF &r){ return(r*ComplexF(0.0,-1.0));} +//accelerator_inline ComplexD timesMinusI(const ComplexD &r){ return(r*ComplexD(0.0,-1.0));} // define projections to real and imaginay parts accelerator_inline ComplexF projReal(const ComplexF &r){return( ComplexF(r.real(), 0.0));} diff --git a/Grid/stencil/SimpleCompressor.h b/Grid/stencil/SimpleCompressor.h index be7c89c0..2ce48369 100644 --- a/Grid/stencil/SimpleCompressor.h +++ b/Grid/stencil/SimpleCompressor.h @@ -7,20 +7,20 @@ template class SimpleCompressor { public: void Point(int) {}; - accelerator_inline int CommDatumSize(void) { return sizeof(vobj); } - accelerator_inline bool DecompressionStep(void) { return false; } - template accelerator_inline void Compress(cobj *buf,int o,const cobj &in) { buf[o]=in; } - accelerator_inline void Exchange(vobj *mp,vobj *vp0,vobj *vp1,Integer type,Integer o){ + accelerator_inline int CommDatumSize(void) const { return sizeof(vobj); } + accelerator_inline bool DecompressionStep(void) const { return false; } + template 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) const { 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, - 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); } // For cshift. Cshift should drop compressor coupling altogether // 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; } }; diff --git a/Grid/stencil/Stencil.h b/Grid/stencil/Stencil.h index 23fc8203..58cebed3 100644 --- a/Grid/stencil/Stencil.h +++ b/Grid/stencil/Stencil.h @@ -147,16 +147,16 @@ class CartesianStencilAccelerator { cobj* u_recv_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; } - 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]; } - 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]; local = this->_entries_p[ent]._is_local; 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]; int local = this->_entries_p[ent]._is_local; if (local) return base + 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); } @@ -221,7 +221,7 @@ public: typedef typename cobj::vector_type vector_type; typedef typename cobj::scalar_type scalar_type; typedef typename cobj::scalar_object scalar_object; - typedef CartesianStencilView View_type; + typedef const CartesianStencilView View_type; typedef typename View_type::StencilVector StencilVector; /////////////////////////////////////////// // Helper structs diff --git a/Grid/tensors/Tensor_SIMT.h b/Grid/tensors/Tensor_SIMT.h index ec57a679..672f385f 100644 --- a/Grid/tensors/Tensor_SIMT.h +++ b/Grid/tensors/Tensor_SIMT.h @@ -64,6 +64,70 @@ void coalescedWriteNonTemporal(vobj & __restrict__ vec,const vobj & __restrict__ } #else + +#ifndef GRID_SYCL +// Use the scalar as our own complex on GPU ... thrust::complex or std::complex +template = 0> accelerator_inline +typename vsimd::scalar_type +coalescedRead(const vsimd & __restrict__ vec,int lane=acceleratorSIMTlane(vsimd::Nsimd())) +{ + typedef typename vsimd::scalar_type S; + S * __restrict__ p=(S *)&vec; + return p[lane]; +} +template = 0> accelerator_inline +typename vsimd::scalar_type +coalescedReadPermute(const vsimd & __restrict__ vec,int doperm,int lane=acceleratorSIMTlane(vsimd::Nsimd())) +{ + typedef typename vsimd::scalar_type S; + + S * __restrict__ p=(S *)&vec; + int mask = vsimd::Nsimd() >> (ptype + 1); + int plane= doperm ? lane ^ mask : lane; + return p[plane]; +} +template = 0> accelerator_inline +void coalescedWrite(vsimd & __restrict__ vec, + const typename vsimd::scalar_type & __restrict__ extracted, + int lane=acceleratorSIMTlane(vsimd::Nsimd())) +{ + typedef typename vsimd::scalar_type S; + S * __restrict__ p=(S *)&vec; + p[lane]=extracted; +} +#else +// For SyCL have option to use GpuComplex from inside the vector type in SIMT loops +// Faster for some reason +template = 0> accelerator_inline +typename vsimd::vector_type::datum +coalescedRead(const vsimd & __restrict__ vec,int lane=acceleratorSIMTlane(vsimd::Nsimd())) +{ + typedef typename vsimd::vector_type::datum S; + S * __restrict__ p=(S *)&vec; + return p[lane]; +} +template = 0> accelerator_inline +typename vsimd::vector_type::datum +coalescedReadPermute(const vsimd & __restrict__ vec,int doperm,int lane=acceleratorSIMTlane(vsimd::Nsimd())) +{ + typedef typename vsimd::vector_type::datum S; + + S * __restrict__ p=(S *)&vec; + int mask = vsimd::Nsimd() >> (ptype + 1); + int plane= doperm ? lane ^ mask : lane; + return p[plane]; +} +template = 0> accelerator_inline +void coalescedWrite(vsimd & __restrict__ vec, + const typename vsimd::vector_type::datum & __restrict__ extracted, + int lane=acceleratorSIMTlane(vsimd::Nsimd())) +{ + typedef typename vsimd::vector_type::datum S; + S * __restrict__ p=(S *)&vec; + p[lane]=extracted; +} +#endif + ////////////////////////////////////////// // Extract and insert slices on the GPU ////////////////////////////////////////// diff --git a/Grid/tensors/Tensor_exp.h b/Grid/tensors/Tensor_exp.h index 0a1d6389..1f637d5f 100644 --- a/Grid/tensors/Tensor_exp.h +++ b/Grid/tensors/Tensor_exp.h @@ -28,7 +28,7 @@ Author: neo #ifndef GRID_MATH_EXP_H #define GRID_MATH_EXP_H -#define DEFAULT_MAT_EXP 12 +#define DEFAULT_MAT_EXP 20 NAMESPACE_BEGIN(Grid); diff --git a/Grid/tensors/Tensor_outer.h b/Grid/tensors/Tensor_outer.h index 4902c22f..a32a2a91 100644 --- a/Grid/tensors/Tensor_outer.h +++ b/Grid/tensors/Tensor_outer.h @@ -34,6 +34,16 @@ NAMESPACE_BEGIN(Grid); // outerProduct Scalar x Scalar -> Scalar // Vector x Vector -> Matrix /////////////////////////////////////////////////////////////////////////////////////// +template = 0> +accelerator_inline CC outerProduct(const CC &l, const CC& r) +{ + return l*conj(r); +} +template = 0> +accelerator_inline RR outerProduct(const RR &l, const RR& r) +{ + return l*r; +} template accelerator_inline auto outerProduct (const iVector& lhs,const iVector& rhs) -> iMatrix @@ -57,17 +67,6 @@ auto outerProduct (const iScalar& lhs,const iScalar& rhs) -> iScalar = 0> -accelerator_inline CC outerProduct(const CC &l, const CC& r) -{ - return l*conj(r); -} -template = 0> -accelerator_inline RR outerProduct(const RR &l, const RR& r) -{ - return l*r; -} - NAMESPACE_END(Grid); #endif diff --git a/Grid/threads/Accelerator.cc b/Grid/threads/Accelerator.cc index 6fc1afe2..9d9d851c 100644 --- a/Grid/threads/Accelerator.cc +++ b/Grid/threads/Accelerator.cc @@ -1,6 +1,7 @@ #include NAMESPACE_BEGIN(Grid); +int acceleratorAbortOnGpuError=1; uint32_t accelerator_threads=2; uint32_t acceleratorThreads(void) {return accelerator_threads;}; void acceleratorThreads(uint32_t t) {accelerator_threads = t;}; @@ -52,7 +53,6 @@ void acceleratorInit(void) prop = gpu_props[i]; totalDeviceMem = prop.totalGlobalMem; if ( world_rank == 0) { -#ifndef GRID_DEFAULT_GPU if ( i==rank ) { printf("AcceleratorCudaInit[%d]: ========================\n",rank); printf("AcceleratorCudaInit[%d]: Device Number : %d\n", rank,i); @@ -66,8 +66,8 @@ void acceleratorInit(void) GPU_PROP(warpSize); GPU_PROP(pciBusID); GPU_PROP(pciDeviceID); + printf("AcceleratorCudaInit[%d]: maxGridSize (%d,%d,%d)\n",rank,prop.maxGridSize[0],prop.maxGridSize[1],prop.maxGridSize[2]); } -#endif // GPU_PROP(unifiedAddressing); // GPU_PROP(l2CacheSize); // GPU_PROP(singleToDoublePrecisionPerfRatio); diff --git a/Grid/threads/Accelerator.h b/Grid/threads/Accelerator.h index 6232aea8..56b85c72 100644 --- a/Grid/threads/Accelerator.h +++ b/Grid/threads/Accelerator.h @@ -100,9 +100,11 @@ void acceleratorInit(void); #define accelerator __host__ __device__ #define accelerator_inline __host__ __device__ inline +extern int acceleratorAbortOnGpuError; + accelerator_inline int acceleratorSIMTlane(int Nsimd) { #ifdef GRID_SIMT - return threadIdx.z; + return threadIdx.x; #else return 0; #endif @@ -110,36 +112,77 @@ accelerator_inline int acceleratorSIMTlane(int Nsimd) { #define accelerator_for2dNB( iter1, num1, iter2, num2, nsimd, ... ) \ { \ + int nt=acceleratorThreads(); \ typedef uint64_t Iterator; \ auto lambda = [=] accelerator \ (Iterator iter1,Iterator iter2,Iterator lane) mutable { \ __VA_ARGS__; \ }; \ - int nt=acceleratorThreads(); \ - dim3 cu_threads(acceleratorThreads(),1,nsimd); \ + dim3 cu_threads(nsimd,acceleratorThreads(),1); \ dim3 cu_blocks ((num1+nt-1)/nt,num2,1); \ LambdaApply<<>>(num1,num2,nsimd,lambda); \ } +#define accelerator_for6dNB(iter1, num1, \ + iter2, num2, \ + iter3, num3, \ + iter4, num4, \ + iter5, num5, \ + iter6, num6, ... ) \ + { \ + typedef uint64_t Iterator; \ + auto lambda = [=] accelerator \ + (Iterator iter1,Iterator iter2, \ + Iterator iter3,Iterator iter4, \ + Iterator iter5,Iterator iter6) mutable { \ + __VA_ARGS__; \ + }; \ + dim3 cu_blocks (num1,num2,num3); \ + dim3 cu_threads(num4,num5,num6); \ + Lambda6Apply<<>>(num1,num2,num3,num4,num5,num6,lambda); \ + } + template __global__ void LambdaApply(uint64_t num1, uint64_t num2, uint64_t num3, lambda Lambda) { - uint64_t x = threadIdx.x + blockDim.x*blockIdx.x; - uint64_t y = threadIdx.y + blockDim.y*blockIdx.y; - uint64_t z = threadIdx.z; + // Weird permute is to make lane coalesce for large blocks + uint64_t x = threadIdx.y + blockDim.y*blockIdx.x; + uint64_t y = threadIdx.z + blockDim.z*blockIdx.y; + uint64_t z = threadIdx.x; if ( (x < num1) && (y __global__ +void Lambda6Apply(uint64_t num1, uint64_t num2, uint64_t num3, + uint64_t num4, uint64_t num5, uint64_t num6, + lambda Lambda) +{ + uint64_t iter1 = blockIdx.x; + uint64_t iter2 = blockIdx.y; + uint64_t iter3 = blockIdx.z; + uint64_t iter4 = threadIdx.x; + uint64_t iter5 = threadIdx.y; + uint64_t iter6 = threadIdx.z; + + if ( (iter1 < num1) && (iter2 global{unum1,unum2,nsimd}; \ cgh.parallel_for( \ cl::sycl::nd_range<3>(global,local), \ - [=] (cl::sycl::nd_item<3> item) mutable { \ + [=] (cl::sycl::nd_item<3> item) /*mutable*/ { \ auto iter1 = item.get_global_id(0); \ auto iter2 = item.get_global_id(1); \ auto lane = item.get_global_id(2); \ @@ -414,7 +457,7 @@ accelerator_inline void acceleratorSynchronise(void) __syncwarp(); #endif #ifdef GRID_SYCL - // No barrier call on SYCL?? // Option get __spir:: stuff to do warp barrier + cl::sycl::detail::workGroupBarrier(); #endif #ifdef GRID_HIP __syncthreads(); diff --git a/Grid/util/CompilerCompatible.h b/Grid/util/CompilerCompatible.h index 37331668..7c4a056d 100644 --- a/Grid/util/CompilerCompatible.h +++ b/Grid/util/CompilerCompatible.h @@ -1,5 +1,16 @@ #pragma once +#if defined(__NVCC__) + +#if (__CUDACC_VER_MAJOR__ == 11) && (__CUDACC_VER_MINOR__ == 0) +#error "NVCC version 11.0 breaks on Ampere, see Github issue 346" +#endif +#if (__CUDACC_VER_MAJOR__ == 11) && (__CUDACC_VER_MINOR__ == 1) +#error "NVCC version 11.1 breaks on Ampere, see Github issue 346" +#endif + +#endif + #if defined(__clang__) #if __clang_major__ < 3 diff --git a/Grid/util/Init.cc b/Grid/util/Init.cc index 9be39e94..bfbc464d 100644 --- a/Grid/util/Init.cc +++ b/Grid/util/Init.cc @@ -140,7 +140,7 @@ void GridCmdOptionCSL(std::string str,std::vector & vec) } template -void GridCmdOptionIntVector(std::string &str,VectorInt & vec) +void GridCmdOptionIntVector(const std::string &str,VectorInt & vec) { vec.resize(0); std::stringstream ss(str); @@ -153,6 +153,9 @@ void GridCmdOptionIntVector(std::string &str,VectorInt & vec) return; } +template void GridCmdOptionIntVector(const std::string &str,std::vector & vec); +template void GridCmdOptionIntVector(const std::string &str,Coordinate & vec); + void GridCmdOptionInt(std::string &str,int & val) { std::stringstream ss(str); diff --git a/Grid/util/Init.h b/Grid/util/Init.h index dad963a0..4eb8f06c 100644 --- a/Grid/util/Init.h +++ b/Grid/util/Init.h @@ -55,7 +55,7 @@ template std::string GridCmdVectorIntToString(const VectorInt & vec); void GridCmdOptionCSL(std::string str,std::vector & vec); template -void GridCmdOptionIntVector(std::string &str,VectorInt & vec); +void GridCmdOptionIntVector(const std::string &str,VectorInt & vec); void GridCmdOptionInt(std::string &str,int & val); diff --git a/HMC/Mobius2p1fRHMC.cc b/HMC/Mobius2p1fRHMC.cc index 82ca4d37..b958d548 100644 --- a/HMC/Mobius2p1fRHMC.cc +++ b/HMC/Mobius2p1fRHMC.cc @@ -56,12 +56,12 @@ int main(int argc, char **argv) { MD.trajL = 1.0; HMCparameters HMCparams; - HMCparams.StartTrajectory = 30; + HMCparams.StartTrajectory = 0; HMCparams.Trajectories = 200; HMCparams.NoMetropolisUntil= 0; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; - // HMCparams.StartingType =std::string("ColdStart"); - HMCparams.StartingType =std::string("CheckpointStart"); + HMCparams.StartingType =std::string("ColdStart"); + // HMCparams.StartingType =std::string("CheckpointStart"); HMCparams.MD = MD; HMCWrapper TheHMC(HMCparams); diff --git a/README.md b/README.md index 4cbae720..fff68dc6 100644 --- a/README.md +++ b/README.md @@ -149,7 +149,6 @@ If you want to build all the tests at once just use `make tests`. - `--enable-numa`: enable NUMA first touch optimisation - `--enable-simd=`: setup Grid for the SIMD target `` (default: `GEN`). A list of possible SIMD targets is detailed in a section below. - `--enable-gen-simd-width=`: select the size (in bytes) of the generic SIMD vector type (default: 32 bytes). -- `--enable-precision={single|double}`: set the default precision (default: `double`). **Deprecated option** - `--enable-comms=`: Use `` for message passing (default: `none`). A list of possible SIMD targets is detailed in a section below. - `--enable-rng={sitmo|ranlux48|mt19937}`: choose the RNG (default: `sitmo `). - `--disable-timers`: disable system dependent high-resolution timers. diff --git a/TODO b/TODO index f1175560..e23e040d 100644 --- a/TODO +++ b/TODO @@ -1,3 +1,6 @@ +-- comms threads issue?? +-- Part done: Staggered kernel performance on GPU + ========================================================= General ========================================================= @@ -5,28 +8,18 @@ General - Make representations code take Gimpl - Simplify the HMCand remove modules - Lattice_arith - are the mult, mac etc.. still needed after ET engine? -- Lattice_rng -- Lattice_transfer.h -- accelerate A2Autils -- off critical path for HMC +- Lattice_rng - faster local only loop in init +- Audit: accelerate A2Autils -- off critical path for HMC ========================================================= -GPU branch code item work list +GPU work list ========================================================= -* sum_cpu promote to double during summation for increased precisoin. +* sum_cpu promote to double during summation for increased precision. * Introduce sumD & ReduceD * GPU sum is probably better currently. - * Accelerate the cshift & benchmark -* 0) Single GPU -- 128 bit integer table load in GPU code. - - ImprovedStaggered accelerate & measure perf - - Gianluca's changes to Cayley into gpu-port - - Mobius kernel fusion. -- Gianluca? - - Lebesque order reintroduction. StencilView should have pointer to it - - Lebesgue reorder in all kernels - * 3) Comms/NVlink - OpenMP tasks to run comms threads. Experiment with it - Remove explicit openMP in staggered. @@ -35,14 +28,6 @@ GPU branch code item work list - Stencil gather ?? - SIMD dirs in stencil -* 4) ET enhancements -- eval -> scalar ops in ET engine -- coalescedRead, coalescedWrite in expressions. - -* 5) Misc -- Conserved current clean up. -- multLinkProp eliminate - 8) Merge develop and test HMC 9) Gamma tables on GPU; check this. Appear to work, but no idea why. Are these done on CPU? @@ -52,7 +37,7 @@ GPU branch code item work list - Audit NAMESPACE CHANGES - Audit changes ------ +--------- Gianluca's changes - Performance impact of construct in aligned allocator??? --------- @@ -62,6 +47,33 @@ Gianluca's changes ----------------------------- DONE: ----------------------------- +===== +-- Done: Remez X^-1/2 X^-1/2 X = 1 test. + Feed in MdagM^2 as a test and take its sqrt. + Automated test that MdagM invsqrt(MdagM)invsqrt(MdagM) = 1 in HMC for bounds satisfaction. + +-- Done: Sycl Kernels into develop. Compare to existing unroll and just use. +-- Done: sRNG into refresh functions +-- Done: Tuned decomposition on CUDA into develop +-- Done: Sycl friend accessor. Const view attempt via typedef?? + + +* Done 5) Misc +- Conserved current clean up. +- multLinkProp eliminate + +* Done 0) Single GPU +- 128 bit integer table load in GPU code. + - ImprovedStaggered accelerate & measure perf + - Gianluca's changes to Cayley into gpu-port + - Mobius kernel fusion. -- Gianluca? + - Lebesque order reintroduction. StencilView should have pointer to it + - Lebesgue reorder in all kernels + +* 4) ET enhancements +- Done eval -> scalar ops in ET engine +- Done coalescedRead, coalescedWrite in expressions. + ============================================================================================= AUDIT ContractWWVV with respect to develop -- DONE - GPU accelerate EOFA -- DONE @@ -125,23 +137,6 @@ AUDIT ContractWWVV with respect to develop -- DONE - - (4) omp parallel for collapse(n) - - Only (1) has a natural mirror in accelerator_loop - - Nested loop macros get cumbersome made a generic interface for N deep -- - Don't like thread_region and thread_loop_in_region -- - Could replace with - - thread_nested(1, - for { - - } - ); - thread_nested(2, - for (){ - for (){ - - } - } - ); - - and same "in_region". ----------------------------- diff --git a/benchmarks/Benchmark_dwf_fp32.cc b/benchmarks/Benchmark_dwf_fp32.cc index cb86177e..03f3ee61 100644 --- a/benchmarks/Benchmark_dwf_fp32.cc +++ b/benchmarks/Benchmark_dwf_fp32.cc @@ -53,7 +53,7 @@ int main (int argc, char ** argv) int threads = GridThread::GetThreads(); Coordinate latt4 = GridDefaultLatt(); - int Ls=8; + int Ls=16; for(int i=0;i> Ls; diff --git a/configure.ac b/configure.ac index f077ca93..4e5e33c8 100644 --- a/configure.ac +++ b/configure.ac @@ -7,7 +7,12 @@ AM_INIT_AUTOMAKE([subdir-objects 1.13]) AM_EXTRA_RECURSIVE_TARGETS([tests bench]) AC_CONFIG_MACRO_DIR([m4]) AC_CONFIG_SRCDIR([Grid/Grid.h]) -AC_CONFIG_HEADERS([Grid/Config.h],[sed -i 's|PACKAGE_|GRID_|' Grid/Config.h]) +AC_CONFIG_HEADERS([Grid/Config.h],[[$SED_INPLACE -e 's|PACKAGE_|GRID_|' -e 's|[[:space:]]PACKAGE[[:space:]]| GRID_PACKAGE |' -e 's|[[:space:]]VERSION[[:space:]]| GRID_PACKAGE_VERSION |' Grid/Config.h]], + [if test x"$host_os" == x"${host_os#darwin}" ; then] + [SED_INPLACE="sed -i"] + [else] + [SED_INPLACE="sed -i .bak"] + [fi]) m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])]) ################ Get git info @@ -125,7 +130,7 @@ esac ############### fermions AC_ARG_ENABLE([fermion-reps], - [AC_HELP_STRING([--fermion-reps=yes|no], [enable extra fermion representation support])], + [AC_HELP_STRING([--enable-fermion-reps=yes|no], [enable extra fermion representation support])], [ac_FERMION_REPS=${enable_fermion_reps}], [ac_FERMION_REPS=yes]) AM_CONDITIONAL(BUILD_FERMION_REPS, [ test "${ac_FERMION_REPS}X" == "yesX" ]) @@ -135,12 +140,23 @@ AC_ARG_ENABLE([gparity], [ac_GPARITY=${enable_gparity}], [ac_GPARITY=yes]) 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 yes) AC_DEFINE([ENABLE_FERMION_REPS],[1],[non QCD fermion reps]);; esac case ${ac_GPARITY} in yes) AC_DEFINE([ENABLE_GPARITY],[1],[fermion actions with GPARITY BCs]);; esac +case ${ac_ZMOBIUS} in + yes) AC_DEFINE([ENABLE_ZMOBIUS],[1],[Zmobius fermion actions]);; +esac ############### Nc AC_ARG_ENABLE([Nc], [AC_HELP_STRING([--enable-Nc=2|3|4], [enable number of colours])], @@ -428,7 +444,7 @@ case ${ax_cv_cxx_compiler_vendor} in SIMD_FLAGS='-mavx2 -mfma -mf16c';; AVX512) AC_DEFINE([AVX512],[1],[AVX512 intrinsics]) - SIMD_FLAGS='-mavx512f -mavx512pf -mavx512er -mavx512cd';; + SIMD_FLAGS='-mavx512f -mavx512cd';; SKL) AC_DEFINE([AVX512],[1],[AVX512 intrinsics for SkyLake Xeon]) SIMD_FLAGS='-march=skylake-avx512';; @@ -481,6 +497,9 @@ case ${ax_cv_cxx_compiler_vendor} in AC_DEFINE([AVX2],[1],[AVX2 intrinsics]) SIMD_FLAGS='-march=core-avx2 -xcore-avx2';; AVX512) + AC_DEFINE([AVX512],[1],[AVX512 intrinsics]) + SIMD_FLAGS='-xcommon-avx512';; + SKL) AC_DEFINE([AVX512],[1],[AVX512 intrinsics]) SIMD_FLAGS='-xcore-avx512';; KNC) diff --git a/tests/core/Test_main.cc b/tests/core/Test_main.cc index d7ed04ba..6e316aa6 100644 --- a/tests/core/Test_main.cc +++ b/tests/core/Test_main.cc @@ -231,6 +231,20 @@ int main(int argc, char **argv) { scalar = localInnerProduct(cVec, cVec); scalar = localNorm2(cVec); + std::cout << "Testing maxLocalNorm2" < shiftcoor = coor; shiftcoor[dir] = (shiftcoor[dir] + shift + latt_size[dir]) % - (latt_size[dir] / mpi_layout[dir]); + (latt_size[dir]); + // (latt_size[dir] / mpi_layout[dir]); std::vector rl(4); for (int dd = 0; dd < 4; dd++) { diff --git a/tests/core/Test_where.cc b/tests/core/Test_where.cc index 050b711b..deb29865 100644 --- a/tests/core/Test_where.cc +++ b/tests/core/Test_where.cc @@ -40,9 +40,9 @@ int main (int argc, char ** argv) int N=16; - std::vector latt_size ({N,4,4}); - std::vector simd_layout({vComplexD::Nsimd(),1,1}); - std::vector mpi_layout ({1,1,1}); + std::vector latt_size ({N,N,N,N}); + std::vector simd_layout({vComplexD::Nsimd(),1,1,1}); + std::vector mpi_layout ({1,1,1,1}); int vol = 1; int nd = latt_size.size(); @@ -69,7 +69,7 @@ int main (int argc, char ** argv) for(int t=0;t +Author: Peter Boyle + + 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. + + 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 */ +#include + +using namespace Grid; + ; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout< latt_size ({N,4,4}); + std::vector simd_layout({vComplexD::Nsimd(),1,1}); + std::vector mpi_layout ({1,1,1}); + + int vol = 1; + int nd = latt_size.size(); + for(int d=0;d({45,12,81,9})); + + std::cout< -void TestConserved(What & Ddwf, What & Ddwfrev, +void TestConserved(What & Ddwf, LatticeGaugeField &Umu, GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, RealD mass, RealD M5, GridParallelRNG *RNG4, - GridParallelRNG *RNG5); + GridParallelRNG *RNG5, + What *Ddwfrev=nullptr); Gamma::Algebra Gmu [] = { Gamma::Algebra::GammaX, @@ -102,10 +103,11 @@ int main (int argc, char ** argv) GridRedBlackCartesian * FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGridF); - std::vector seeds4({1,2,3,4}); std::vector seeds5({5,6,7,8}); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); - GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG RNG4(UGrid); + std::vector seeds4({1,2,3,4}); RNG4.SeedFixedIntegers(seeds4); + //const std::string seeds4{ "test-gauge-3000" }; RNG4.SeedUniqueString( seeds4 ); LatticeGaugeField Umu(UGrid); if( argc > 1 && argv[1][0] != '-' ) @@ -116,8 +118,8 @@ int main (int argc, char ** argv) } else { - std::cout<::ColdConfiguration(Umu); + std::cout<::ColdConfiguration(Umu); SU::HotConfiguration(RNG4,Umu); } @@ -127,7 +129,7 @@ int main (int argc, char ** argv) std::cout<(Ddwf,Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestConserved(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); RealD b=1.5;// Scale factor b+c=2, b-c=1 RealD c=0.5; @@ -137,13 +139,13 @@ int main (int argc, char ** argv) std::cout<(Dmob,Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestConserved(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dsham,Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestConserved(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(ZDmob,ZDmobrev,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestConserved(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5,&ZDmobrev); Grid_finalize(); } @@ -161,22 +162,17 @@ int main (int argc, char ** argv) template -void TestConserved(Action & Ddwf, - Action & Ddwfrev, +void TestConserved(Action & Ddwf, LatticeGaugeField &Umu, GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, RealD mass, RealD M5, GridParallelRNG *RNG4, - GridParallelRNG *RNG5) + GridParallelRNG *RNG5, + Action * Ddwfrev) { - int Ls=Ddwf.Ls; - - LatticePropagator phys_src(UGrid); - - std::vector U(4,UGrid); - - LatticePropagator seqsrc(FGrid); + LatticePropagator phys_src(UGrid); + LatticePropagator seqsrc(FGrid); LatticePropagator prop5(FGrid); LatticePropagator prop5rev(FGrid); LatticePropagator prop4(UGrid); @@ -194,9 +190,9 @@ void TestConserved(Action & Ddwf, phys_src=Zero(); pokeSite(kronecker,phys_src,coor); - MdagMLinearOperator HermOp(Ddwf); - MdagMLinearOperator HermOprev(Ddwfrev); ConjugateGradient CG(1.0e-16,100000); + SchurRedBlackDiagTwoSolve schur(CG); + ZeroGuesser zpg; for(int s=0;s(prop5,result5,s,c); LatticeFermion result4(UGrid); Ddwf.ExportPhysicalFermionSolution(result5,result4); FermToProp(prop4,result4,s,c); - Ddwfrev.ImportPhysicalFermionSource(src4,src5); - Ddwfrev.Mdag(src5,Mdagsrc5); - CG(HermOprev,Mdagsrc5,result5); + if( Ddwfrev ) { + Ddwfrev->ImportPhysicalFermionSource(src4,src5); + result5 = Zero(); + schur(*Ddwfrev,src5,result5,zpg); + } FermToProp(prop5rev,result5,s,c); } } @@ -251,11 +247,7 @@ void TestConserved(Action & Ddwf, PropToFerm(src5,seqsrc,s,c); LatticeFermion result5(FGrid); result5=Zero(); - - // CGNE - LatticeFermion Mdagsrc5 (FGrid); - Ddwf.Mdag(src5,Mdagsrc5); - CG(HermOp,Mdagsrc5,result5); + schur(Ddwf,src5,result5,zpg); LatticeFermion result4(UGrid); Ddwf.ExportPhysicalFermionSolution(result5,result4); @@ -276,10 +268,10 @@ void TestConserved(Action & Ddwf, Ddwf.ContractConservedCurrent(prop5rev,prop5,Vector_mu,phys_src,Current::Vector,Tdir); Ddwf.ContractJ5q(prop5,PJ5q); - PA = trace(g5*Axial_mu); - SV = trace(Vector_mu); - VV = trace(gT*Vector_mu); - PP = trace(adj(prop4)*prop4); + PA = trace(g5*Axial_mu); // Pseudoscalar-Axial conserved current + SV = trace(Vector_mu); // Scalar-Vector conserved current + VV = trace(gT*Vector_mu); // (local) Vector-Vector conserved current + PP = trace(adj(prop4)*prop4); // Pseudoscalar density // Spatial sum sliceSum(PA,sumPA,Tdir); @@ -288,15 +280,17 @@ void TestConserved(Action & Ddwf, sliceSum(PP,sumPP,Tdir); sliceSum(PJ5q,sumPJ5q,Tdir); - int Nt=sumPA.size(); + const int Nt{static_cast(sumPA.size())}; + std::cout< CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, true); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5 ); + RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_dwf_gpforce.cc b/tests/forces/Test_dwf_gpforce.cc index 28133cc6..1fa1c6e4 100644 --- a/tests/forces/Test_dwf_gpforce.cc +++ b/tests/forces/Test_dwf_gpforce.cc @@ -84,6 +84,13 @@ int main (int argc, char ** argv) GparityDomainWallFermionR::ImplParams params; params.twists = twists; + /* + params.boundary_phases[0] = 1.0; + params.boundary_phases[1] = 1.0; + params.boundary_phases[2] = 1.0; + params.boundary_phases[3] =- 1.0; + */ + GparityDomainWallFermionR Dw(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); Dw.M (phi,Mphi); @@ -96,6 +103,16 @@ int main (int argc, char ** argv) Dw.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; Dw.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp); + + // ***************************************************************************************** + // *** There is a funny negative sign in all derivatives. This is - UdSdU. *** + // *** *** + // *** Deriv in both Wilson gauge action and the TwoFlavour.h seems to miss a minus sign *** + // *** UdSdU is negated relative to what I think - call what is returned mUdSdU, *** + // *** and insert minus sign *** + // ***************************************************************************************** + + UdSdU = - UdSdU ; // Follow sign convention of actions in Grid. Seems crazy. FermionField Ftmp (FGrid); @@ -106,7 +123,7 @@ int main (int argc, char ** argv) RealD Hmom = 0.0; RealD Hmomprime = 0.0; LatticeColourMatrix mommu(UGrid); - LatticeColourMatrix forcemu(UGrid); + LatticeColourMatrix mUdSdUmu(UGrid); LatticeGaugeField mom(UGrid); LatticeGaugeField Uprime(UGrid); @@ -114,10 +131,20 @@ int main (int argc, char ** argv) SU::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg - Hmom -= real(sum(trace(mommu*mommu))); + // Momentum Hamiltonian is - trace(p^2)/HMC_MOM_DENOMINATOR + // + // Integrator.h: RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom // GaugeImplTypes.h: Hloc += trace(Pmu * Pmu); + // Sign comes from a sneaky multiply by "i" in GaussianFundemantalLie algebra + // P is i P^a_\mu T^a, not Pa Ta + // + // Integrator.h: H = Hmom + sum S(action) + Hmom -= real(sum(trace(mommu*mommu)))/ HMC_MOMENTUM_DENOMINATOR; PokeIndex(mom,mommu,mu); + // -- Drops factor of "i" in the U update: U' = e^{P dt} U [ _not_ e^{iPdt}U ]. P is anti hermitian already + // -- Udot = p U + // fourth order exponential approx autoView( mom_v, mom, CpuRead); autoView( U_v , U, CpuRead); @@ -134,8 +161,8 @@ int main (int argc, char ** argv) ; }); } - std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <(mom,mu); - std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<(UdSdU,mu); - std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<(UdSdU,mu); - mommu=Ta(mommu)*2.0; + mommu=Ta(mommu); // projectForce , GaugeImplTypes.h PokeIndex(UdSdU,mommu,mu); } for(int mu=0;mu(mom,mu); - std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<(UdSdU,mu); - std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<(UdSdU,mu); + mUdSdUmu= PeekIndex(UdSdU,mu); mommu = PeekIndex(mom,mu); - // Update PF action density - dS = dS+trace(mommu*forcemu)*dt; + // + // Derive HMC eom: + // + // Sdot = - 2 trace( p p^dot ) / D - trace( p [ mUdSdU - h.c. ] ) = 0 + // + // + // Sdot = 0 = - 2 trace( p p^dot ) / D - 2 trace( p Ta( mUdSdU ) = 0 + // + // EOM: + // + // pdot = - D Ta( mUdSdU ) -- source of sign is the "funny sign" above + // + // dSqcd_dt = - 2.0*trace(mommu* Ta(mUdSdU) )*dt -- i.e. mUdSdU with adjoint term -> force has a 2x implicit + // + // dH_mom/dt = - 2 trace (p pdot)/Denom + // + // dH_tot / dt = 0 <= pdot = - Denom * mUdSdU + // + // dH_mom/dt = 2 trace (p mUdSdU ) + // + // True Momentum delta H has a dt^2 piece + // + // dSmom = [ trace mom*mom - trace ( (mom-Denom*f*dt)(mom-Denom*f*dt) ) ] / Denom + // = 2*trace(mom*f) dt - Denom*dt*dt * trace(f*f). + // = dSmom + dSmom2 + // - dSmom = dSmom - trace(mommu*forcemu) * dt; - dSmom2 = dSmom2 - trace(forcemu*forcemu) *(0.25* dt*dt); + dS = dS - 2.0*trace(mommu*mUdSdUmu)*dt; // U and Udagger derivs hence 2x. - // Update mom action density - mommu = mommu + forcemu*(dt*0.5); + dSmom = dSmom + 2.0*trace(mommu*mUdSdUmu) * dt; // this 2.0 coms from derivative of p^2 + + dSmom2 = dSmom2 - trace(mUdSdUmu*mUdSdUmu) * dt*dt* HMC_MOMENTUM_DENOMINATOR; // Remnant - Hmomprime -= real(sum(trace(mommu*mommu))); + // Update mom action density . Verbatim update_P in Integrator.h + mommu = mommu - mUdSdUmu * dt* HMC_MOMENTUM_DENOMINATOR;; + + Hmomprime -= real(sum(trace(mommu*mommu))) / HMC_MOMENTUM_DENOMINATOR; } @@ -199,20 +233,25 @@ int main (int argc, char ** argv) ComplexD dSm = sum(dSmom); ComplexD dSm2 = sum(dSmom2); + std::cout << GridLogMessage <<"dSm "<< dSm< CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, true); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5); RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_laplacian_force.cc b/tests/forces/Test_laplacian_force.cc index 18508860..dbaf1cbd 100644 --- a/tests/forces/Test_laplacian_force.cc +++ b/tests/forces/Test_laplacian_force.cc @@ -46,6 +46,7 @@ int main (int argc, char ** argv) std::vector seeds({1,2,3,4}); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers({4,5,6,7}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(std::vector({15,91,21,3})); @@ -67,7 +68,7 @@ int main (int argc, char ** argv) LaplacianAdjointField Laplacian(&Grid, CG, LapPar, Kappa); GeneralisedMomenta LaplacianMomenta(&Grid, Laplacian); LaplacianMomenta.M.ImportGauge(U); - LaplacianMomenta.MomentaDistribution(pRNG);// fills the Momenta with the correct distr + LaplacianMomenta.MomentaDistribution(sRNG,pRNG);// fills the Momenta with the correct distr std::cout << std::setprecision(15); diff --git a/tests/forces/Test_mobius_force.cc b/tests/forces/Test_mobius_force.cc index ba7bc363..d2326a81 100644 --- a/tests/forces/Test_mobius_force.cc +++ b/tests/forces/Test_mobius_force.cc @@ -69,7 +69,14 @@ int main (int argc, char ** argv) RealD M5=1.8; RealD b=0.5; RealD c=0.5; - MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); + + WilsonImplParams p; + p.boundary_phases[0] = 1.0; + p.boundary_phases[1] = 1.0; + p.boundary_phases[2] = 1.0; + p.boundary_phases[3] =- 1.0; + + MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,p); Ddwf.M (phi,Mphi); ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p @@ -82,24 +89,44 @@ int main (int argc, char ** argv) Ddwf.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; Ddwf.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp); + // ***************************************************************************************** + // *** There is a funny negative sign in all derivatives. This is - UdSdU. *** + // *** *** + // *** Deriv in both Wilson gauge action and the TwoFlavour.h seems to miss a minus sign *** + // *** UdSdU is negated relative to what I think - call what is returned mUdSdU, *** + // *** and insert minus sign *** + // ***************************************************************************************** + + UdSdU = - UdSdU ; // Follow sign convention of actions in Grid. Seems crazy. + LatticeFermion Ftmp (FGrid); //////////////////////////////////// // Modify the gauge field a little //////////////////////////////////// - RealD dt = 0.0001; + RealD dt = 0.001; + RealD Hmom = 0.0; + RealD Hmomprime = 0.0; LatticeColourMatrix mommu(UGrid); - LatticeColourMatrix forcemu(UGrid); + LatticeColourMatrix mUdSdUmu(UGrid); LatticeGaugeField mom(UGrid); LatticeGaugeField Uprime(UGrid); for(int mu=0;mu::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg - PokeIndex(mom,mommu,mu); + // Momentum Hamiltonian is - trace(p^2)/HMC_MOM_DENOMINATOR + // + // Integrator.h: RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom // GaugeImplTypes.h: Hloc += trace(Pmu * Pmu); + // Sign comes from a sneaky multiply by "i" in GaussianFundemantalLie algebra + // P is i P^a_\mu T^a, not Pa Ta + // + // Integrator.h: H = Hmom + sum S(action) + Hmom -= real(sum(trace(mommu*mommu)))/ HMC_MOMENTUM_DENOMINATOR; + // fourth order exponential approx autoView( U_v , U, CpuRead); autoView( mom_v, mom, CpuRead); @@ -115,6 +142,7 @@ int main (int argc, char ** argv) ; }); } + std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <(UdSdU,mu); - mommu=Ta(mommu)*2.0; + mommu=Ta(mommu); PokeIndex(UdSdU,mommu,mu); } for(int mu=0;mu(UdSdU,mu); + + mUdSdUmu= PeekIndex(UdSdU,mu); mommu = PeekIndex(mom,mu); - // Update PF action density - dS = dS+trace(mommu*forcemu)*dt; + // + // Derive HMC eom: + // + // Sdot = - 2 trace( p p^dot ) / D - trace( p [ mUdSdU - h.c. ] ) = 0 + // + // + // Sdot = 0 = - 2 trace( p p^dot ) / D - 2 trace( p Ta( mUdSdU ) = 0 + // + // EOM: + // + // pdot = - D Ta( mUdSdU ) -- source of sign is the "funny sign" above + // + // dSqcd_dt = - 2.0*trace(mommu* Ta(mUdSdU) )*dt -- i.e. mUdSdU with adjoint term -> force has a 2x implicit + // + // dH_mom/dt = - 2 trace (p pdot)/Denom + // + // dH_tot / dt = 0 <= pdot = - Denom * mUdSdU + // + // dH_mom/dt = 2 trace (p mUdSdU ) + // + // True Momentum delta H has a dt^2 piece + // + // dSmom = [ trace mom*mom - trace ( (mom-Denom*f*dt)(mom-Denom*f*dt) ) ] / Denom + // = 2*trace(mom*f) dt - Denom*dt*dt * trace(f*f). + // = dSmom + dSmom2 + // + + dS = dS - 2.0*trace(mommu*mUdSdUmu)*dt; // U and Udagger derivs hence 2x. + + dSmom = dSmom + 2.0*trace(mommu*mUdSdUmu) * dt; // this 2.0 coms from derivative of p^2 + + dSmom2 = dSmom2 - trace(mUdSdUmu*mUdSdUmu) * dt*dt* HMC_MOMENTUM_DENOMINATOR; // Remnant + + mommu = mommu - mUdSdUmu * dt* HMC_MOMENTUM_DENOMINATOR;; + + Hmomprime -= real(sum(trace(mommu*mommu))) / HMC_MOMENTUM_DENOMINATOR; } ComplexD dSpred = sum(dS); + ComplexD dSm = sum(dSmom); + ComplexD dSm2 = sum(dSmom2); - std::cout << GridLogMessage << " -- S "< CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5 ); RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_mobius_gpforce_eofa.cc b/tests/forces/Test_mobius_gpforce_eofa.cc index 9c80b2aa..7f114615 100644 --- a/tests/forces/Test_mobius_gpforce_eofa.cc +++ b/tests/forces/Test_mobius_gpforce_eofa.cc @@ -93,7 +93,8 @@ int main (int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5 ); RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_momentum_filter.cc b/tests/forces/Test_momentum_filter.cc new file mode 100644 index 00000000..856ea0f2 --- /dev/null +++ b/tests/forces/Test_momentum_filter.cc @@ -0,0 +1,154 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_force.cc + + Copyright (C) 2015 + +Author: Christopher Kelly + + 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. + + 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 */ +#include + +using namespace std; +using namespace Grid; + +//Get the mu-directected links on the upper boundary and the bulk remainder +template +void getLinksBoundaryBulk(Field &bound, Field &bulk, Field &from, const Coordinate &latt_size){ + bound = Zero(); bulk = Zero(); + for(int mu=0;mu seeds({1,2,3,4}); + + GridParallelRNG pRNG(&Grid); + pRNG.SeedFixedIntegers(seeds); + + typedef PeriodicGimplR Gimpl; + typedef WilsonGaugeAction GaugeAction; + typedef NoHirep Representation; //fundamental + typedef NoSmearing Smearing; + typedef MinimumNorm2 Omelyan; + typedef Gimpl::Field Field; + typedef MomentumFilterApplyPhase Filter; + Filter filter(&Grid); + + //Setup a filter that disables link update on links passing through the global lattice boundary + typedef Filter::LatticeLorentzScalarType MaskType; + typedef Filter::LorentzScalarType MaskSiteType; + + MaskSiteType zero, one; + for(int mu=0;mu::HotConfiguration(pRNG,U); + + //Get the original links on the bulk and boundary for later use + Field Ubnd_orig(&Grid), Ubulk_orig(&Grid); + getLinksBoundaryBulk(Ubnd_orig, Ubulk_orig, U, latt_size); + + ActionSet actions(1); + double beta=6; + GaugeAction gauge_action(beta); + actions[0].push_back(&gauge_action); + + Smearing smear; + IntegratorParameters params(1,1.); //1 MD step + Omelyan integrator(&Grid, params, actions, smear); + + integrator.setMomentumFilter(filter); + + integrator.refresh(U, pRNG); //doesn't actually change the gauge field + + //Check the momentum is zero on the boundary + const auto &P = integrator.getMomentum(); + Field Pbnd(&Grid), Pbulk(&Grid); + getLinksBoundaryBulk(Pbnd, Pbulk, const_cast(P), latt_size); + + RealD Pbnd_nrm = norm2(Pbnd); //expect zero + std::cout << GridLogMessage << "After refresh, norm2 of mu-directed conjugate momentum on boundary is: " << Pbnd_nrm << " (expect 0)" << std::endl; + RealD Pbulk_nrm = norm2(Pbulk); //expect non-zero + std::cout << GridLogMessage << "After refresh, norm2 of bulk conjugate momentum is: " << Pbulk_nrm << " (expect non-zero)" << std::endl; + + //Evolve the gauge field + integrator.integrate(U); + + //Check momentum is still zero on boundary + getLinksBoundaryBulk(Pbnd, Pbulk, const_cast(P), latt_size); + + Pbnd_nrm = norm2(Pbnd); //expect zero + std::cout << GridLogMessage << "After integrate, norm2 of mu-directed conjugate momentum on boundary is: " << Pbnd_nrm << " (expect 0)" << std::endl; + Pbulk_nrm = norm2(Pbulk); //expect non-zero + std::cout << GridLogMessage << "After integrate, norm2 of bulk conjugate momentum is: " << Pbulk_nrm << " (expect non-zero)" << std::endl; + + //Get the new bulk and bound links + Field Ubnd_new(&Grid), Ubulk_new(&Grid); + getLinksBoundaryBulk(Ubnd_new, Ubulk_new, U, latt_size); + + Field Ubnd_diff = Ubnd_new - Ubnd_orig; + Field Ubulk_diff = Ubulk_new - Ubulk_orig; + + RealD Ubnd_change = norm2( Ubnd_diff ); + RealD Ubulk_change = norm2( Ubulk_diff ); + std::cout << GridLogMessage << "After integrate, norm2 of change in mu-directed boundary links is : " << Ubnd_change << " (expect 0)" << std::endl; + std::cout << GridLogMessage << "After integrate, norm2 of change in bulk links is : " << Ubulk_change << " (expect non-zero)" << std::endl; + + Grid_finalize(); +}