From 7d3512ab219747f443e6416c7222d16d14765b07 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 14 Aug 2015 00:01:04 +0100 Subject: [PATCH] Gparity valence test now working. Interface in FermionOperator will change a lot in future --- lib/Stencil.h | 16 +- lib/algorithms/CoarsenedMatrix.h | 20 +- lib/lattice/Lattice_peekpoke.h | 2 + lib/qcd/QCD.h | 2 +- lib/qcd/action/Actions.h | 11 + lib/qcd/action/fermion/FermionOperator.h | 66 ++-- lib/qcd/action/fermion/WilsonCompressor.h | 75 ++-- lib/qcd/action/fermion/WilsonFermion.cc | 6 +- lib/qcd/action/fermion/WilsonFermion5D.cc | 11 +- lib/qcd/action/fermion/WilsonKernels.cc | 399 ++++++++++------------ lib/stencil/Stencil_common.cc | 72 ++-- tests/Make.inc | 6 +- tests/Test_gparity.cc | 237 +++++++++++-- tests/Test_stencil.cc | 20 +- 14 files changed, 572 insertions(+), 371 deletions(-) diff --git a/lib/Stencil.h b/lib/Stencil.h index 8529e73a..c6e51059 100644 --- a/lib/Stencil.h +++ b/lib/Stencil.h @@ -41,6 +41,12 @@ namespace Grid { + struct StencilEntry { + int _offset; + int _is_local; + int _permute; + int _around_the_world; + }; class CartesianStencil { // Stencil runs along coordinate axes only; NO diagonal fill in. public: @@ -58,9 +64,9 @@ namespace Grid { std::vector _permute_type; // npoints x Osites() of these - std::vector > _offsets; - std::vector > _is_local; - std::vector > _permute; + std::vector > _entries; + + inline StencilEntry * GetEntry(int &ptype,int point,int osite) { ptype = _permute_type[point]; return & _entries[point][osite]; } int _unified_buffer_size; int _request_count; @@ -77,8 +83,8 @@ namespace Grid { // Can this be avoided with simpler coding of comms? void Local (int point, int dimension,int shift,int cbmask); void Comms (int point, int dimension,int shift,int cbmask); - void CopyPlane(int point, int dimension,int lplane,int rplane,int cbmask,int permute); - void ScatterPlane (int point,int dimension,int plane,int cbmask,int offset); + void CopyPlane(int point, int dimension,int lplane,int rplane,int cbmask,int permute,int wrap); + void ScatterPlane (int point,int dimension,int plane,int cbmask,int offset,int wrap); // Could allow a functional munging of the halo to another type during the comms. // this could implement the 16bit/32bit/64bit compression. diff --git a/lib/algorithms/CoarsenedMatrix.h b/lib/algorithms/CoarsenedMatrix.h index c4b039f1..a854bba7 100644 --- a/lib/algorithms/CoarsenedMatrix.h +++ b/lib/algorithms/CoarsenedMatrix.h @@ -193,20 +193,18 @@ PARALLEL_FOR_LOOP for(int ss=0;ssoSites();ss++){ siteVector res = zero; siteVector nbr; - int offset,local,perm,ptype; - + int ptype; + StencilEntry *SE; for(int point=0;point_is_local&&SE->_permute) { + permute(nbr,in._odata[SE->_offset],ptype); + } else if(SE->_is_local) { + nbr = in._odata[SE->_offset]; } else { - nbr = comm_buf[offset]; + nbr = comm_buf[SE->_offset]; } res = res + A[point]._odata[ss]*nbr; } diff --git a/lib/lattice/Lattice_peekpoke.h b/lib/lattice/Lattice_peekpoke.h index e60618c5..8ac952ea 100644 --- a/lib/lattice/Lattice_peekpoke.h +++ b/lib/lattice/Lattice_peekpoke.h @@ -14,6 +14,7 @@ namespace Grid { auto PeekIndex(const Lattice &lhs,int i) -> Lattice(lhs._odata[0],i))> { Lattice(lhs._odata[0],i))> ret(lhs._grid); + ret.checkerboard=lhs.checkerboard; PARALLEL_FOR_LOOP for(int ss=0;ssoSites();ss++){ ret._odata[ss] = peekIndex(lhs._odata[ss],i); @@ -24,6 +25,7 @@ PARALLEL_FOR_LOOP auto PeekIndex(const Lattice &lhs,int i,int j) -> Lattice(lhs._odata[0],i,j))> { Lattice(lhs._odata[0],i,j))> ret(lhs._grid); + ret.checkerboard=lhs.checkerboard; PARALLEL_FOR_LOOP for(int ss=0;ssoSites();ss++){ ret._odata[ss] = peekIndex(lhs._odata[ss],i,j); diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index 16920f2b..1286c6bd 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -36,7 +36,7 @@ namespace QCD { template struct isSpinor { static const bool value = (SpinorIndex==T::TensorLevel); }; - template using IfSpinor = Invoke::value,int> > ; + template using IfSpinor = Invoke::value,int> > ; template using IfNotSpinor = Invoke::value,int> > ; // ChrisK very keen to add extra space for Gparity doubling. diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index 9b119858..9976c670 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -46,6 +46,8 @@ //////////////////////////////////////////////////////////////////////////////////////////////////// #define FermOpTemplateInstantiate(A) \ + template class A; \ + template class A; \ template class A; \ template class A; @@ -131,6 +133,15 @@ typedef OverlapWilsonPartialFractionTanhFermion OverlapWilsonPartia typedef OverlapWilsonPartialFractionZolotarevFermion OverlapWilsonPartialFractionZolotarevFermionR; typedef OverlapWilsonPartialFractionZolotarevFermion OverlapWilsonPartialFractionZolotarevFermionF; typedef OverlapWilsonPartialFractionZolotarevFermion OverlapWilsonPartialFractionZolotarevFermionD; + +// Gparity cases; partial list until tested +typedef WilsonFermion GparityWilsonFermionR; +typedef WilsonFermion GparityWilsonFermionF; +typedef WilsonFermion GparityWilsonFermionD; +typedef DomainWallFermion GparityDomainWallFermionR; +typedef DomainWallFermion GparityDomainWallFermionF; +typedef DomainWallFermion GparityDomainWallFermionD; + }} /////////////////////////////////////////////////////////////////////////////// // G5 herm -- this has to live in QCD since dirac matrix is not in the broader sector of code diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 75641b83..e1f65a27 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -37,12 +37,10 @@ namespace Grid { typedef WilsonCompressor Compressor; - // provide the multiply by link that is differentiated between Gparity (with flavour index) and - // non-Gparity - static inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu){ + // provide the multiply by link that is differentiated between Gparity (with flavour index) and non-Gparity + static inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu,StencilEntry *SE){ mult(&phi(),&U(mu),&chi()); } - static inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) { conformable(Uds._grid,GaugeGrid); @@ -55,6 +53,11 @@ namespace Grid { PokeIndex(Uds,U,mu+4); } } + static inline void InsertForce(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu){ + GaugeLinkField link(mat._grid); + link = TraceIndex(outerProduct(Btilde,A)); + PokeIndex(mat,link,mu); + } }; @@ -62,8 +65,6 @@ namespace Grid { typedef WilsonImpl WilsonImplF; // Float typedef WilsonImpl WilsonImplD; // Double - - template class GparityWilsonImpl { public: @@ -74,7 +75,7 @@ namespace Grid { template using iImplHalfSpinor = iVector, Nhs>, Ngp >; template using iImplGaugeField = iVector >, Nd >; - template using iImplGaugeLink = iScalar > > >; + template using iImplGaugeLink = iScalar > >; template using iImplDoubledGaugeField = iVector >, Nds >, Ngp >; typedef iImplSpinor SiteSpinor; @@ -88,16 +89,26 @@ namespace Grid { typedef Lattice GaugeField; typedef Lattice DoubledGaugeField; - typedef GparityWilsonCompressor Compressor; + // typedef GparityWilsonCompressor Compressor; + typedef WilsonCompressor Compressor; // provide the multiply by link that is differentiated between Gparity (with flavour index) and // non-Gparity - static inline void multLink(SiteHalfSpinor &phi,const SiteDoubledGaugeField &U,const SiteHalfSpinor &chi,int mu){ - for(int f=0;f_around_the_world && ((mu==Xp)||(mu==Xm)) ) { + mult(&phi(0),&U(0)(mu),&chi(1)); + mult(&phi(1),&U(1)(mu),&chi(0)); + } else { + mult(&phi(0),&U(0)(mu),&chi(0)); + mult(&phi(1),&U(1)(mu),&chi(1)); } + } + static inline void InsertForce(GaugeField &mat,const FermionField &Btilde,const FermionField &A,int mu){ + // Fixme + return; + } static inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) { conformable(Uds._grid,GaugeGrid); @@ -116,44 +127,51 @@ namespace Grid { LatticeCoordinate(coor,mu); U = PeekIndex(Umu,mu); - Uconj = conj(U); + Uconj = conjugate(U); int neglink = GaugeGrid->GlobalDimensions()[mu]-1; if ( gpdirs[mu] ) { Uconj = where(coor==neglink,-Uconj,Uconj); } + PARALLEL_FOR_LOOP for(auto ss=U.begin();ss WilsonImplR; // Real.. whichever prec - typedef WilsonImpl WilsonImplF; // Float - typedef WilsonImpl WilsonImplD; // Double + typedef GparityWilsonImpl GparityWilsonImplR; // Real.. whichever prec + typedef GparityWilsonImpl GparityWilsonImplF; // Float + typedef GparityWilsonImpl GparityWilsonImplD; // Double diff --git a/lib/qcd/action/fermion/WilsonCompressor.h b/lib/qcd/action/fermion/WilsonCompressor.h index c879e9f1..f2aae8aa 100644 --- a/lib/qcd/action/fermion/WilsonCompressor.h +++ b/lib/qcd/action/fermion/WilsonCompressor.h @@ -63,6 +63,7 @@ namespace QCD { } }; + /* template class GparityWilsonCompressor : public WilsonCompressor{ @@ -77,55 +78,65 @@ namespace QCD { const int Nsimd = grid->Nsimd(); - int ocb=grid->CheckerBoardFromOindex(osite); + int checkered=grid->CheckerBoarded(dim); + int ocb=grid->CheckerBoardFromOindex(osite); - SiteHalfSpinor tmp = spinproject(in); // spin projected + SiteHalfSpinor tmp = this->spinproject(in); // spin projected ////////////////////////////////////////////////////////////// // Check whether we must flavour flip // do this if source is plane 0 on processor 0 in dimension dim ////////////////////////////////////////////////////////////// - if ( (this->mu==Xp)||(this->mu==Yp)||(this->mu==Zp)||(this->mu==Tp) ) { + int do_flip = 0; + int flipicoor= 0; + if(Gbcs[this->mu]){ - if ( Gbcs[this->mu] && (grid->_processor_coor[dim] == 0) && (plane==0) && (ocb==0) ) { + std::cout << "Applying Gparity BC's in direction "<mu< flat(Nsimd); - - extract(tmp,flat); - - for(int i=0;i coor; - - grid->iCoorFromIindex(coor,i); - - if ( coor[dim]==0 ) { - scalar_object stmp; - stmp(0) = flat[i](1); - stmp(1) = flat[i](0); - flat[i] = stmp; - } + if ( (this->mu==Xp)||(this->mu==Yp)||(this->mu==Zp)||(this->mu==Tp) ) { + if ( (grid->_processor_coor[dim] == 0) + && (plane==0) + && ( (!checkered)||(ocb==0) ) ) { + do_flip=1; + flipicoor=0; } + } + if ( (this->mu==Xm)||(this->mu==Ym)||(this->mu==Zm)||(this->mu==Tm) ) { + if ( (grid->_processor_coor[dim] == (grid->_processors[dim]-1) ) + && (plane==grid->_rdimensions[dim]-1) + && ( (!checkered)||(ocb==1) ) ) { + do_flip=1; + flipicoor=grid->_simd_layout[dim]-1; + } + } + } + if ( do_flip ) { - merge(tmp,flat); + std::cout << "Applying Gparity BC's in direction "<mu<< " osite " << osite << " plane "< flat(Nsimd); + std::vector coor; + + extract(tmp,flat); + for(int i=0;iiCoorFromIindex(coor,i); + if ( coor[dim]==flipicoor ) { + scalar_object stmp; + stmp(0) = flat[i](1); + stmp(1) = flat[i](0); + flat[i] = stmp; + } + } + merge(tmp,flat); + + } return tmp; } - SiteHalfSpinor flavourflip(const SiteHalfSpinor &in) { - - SiteHalfSpinor ret; - for(int f=0;f(B,comm_buf,compressor); @@ -141,9 +140,8 @@ PARALLEL_FOR_LOOP ////////////////////////////////////////////////// // spin trace outer product ////////////////////////////////////////////////// - tmp = TraceIndex(outerProduct(Btilde,A)); - PokeIndex(mat,tmp,mu); - + Impl::InsertForce(mat,Btilde,A,mu); + } } diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index f324db8d..10caad13 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -137,7 +137,6 @@ void WilsonFermion5D::DerivInternal(CartesianStencil & st, //////////////////////// // Call the single hop //////////////////////// - tmp = zero; PARALLEL_FOR_LOOP for(int sss=0;sssoSites();sss++){ @@ -154,14 +153,18 @@ PARALLEL_FOR_LOOP // spin trace outer product //////////////////////////// - tmp[sU] = tmp[sU]+ traceIndex( - outerProduct(Btilde[sF],Atilde[sF])); // ordering here - } } + Impl::InsertForce(mat,Btilde,A,mu); + /* + tmp = zero; + for(int sss=0;sssoSites();sss++){ + tmp[sU] = tmp[sU]+ traceIndex(outerProduct(Btilde[sF],Atilde[sF])); // ordering here + } PokeIndex(mat,tmp,mu); + */ } } diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 6d558e03..d098d138 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -9,143 +9,116 @@ void WilsonKernels::DiracOptDhopSite(CartesianStencil &st,DoubledGaugeFiel { SiteHalfSpinor tmp; SiteHalfSpinor chi; - SiteSpinor result; SiteHalfSpinor Uchi; - int offset,local,perm, ptype; + SiteSpinor result; + StencilEntry *SE; + int ptype; // Xp - int ss = sF; - offset = st._offsets [Xp][ss]; - local = st._is_local[Xp][ss]; - perm = st._permute[Xp][ss]; - - ptype = st._permute_type[Xp]; - if ( local && perm ) { - spProjXp(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Xp,sF); + if ( SE->_is_local && SE->_permute ) { + spProjXp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjXp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjXp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Xp); + Impl::multLink(Uchi,U._odata[sU],chi,Xp,SE); spReconXp(result,Uchi); // Yp - offset = st._offsets [Yp][ss]; - local = st._is_local[Yp][ss]; - perm = st._permute[Yp][ss]; - ptype = st._permute_type[Yp]; - if ( local && perm ) { - spProjYp(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Yp,sF); + if ( SE->_is_local && SE->_permute ) { + spProjYp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjYp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjYp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Yp); + Impl::multLink(Uchi,U._odata[sU],chi,Yp,SE); accumReconYp(result,Uchi); // Zp - offset = st._offsets [Zp][ss]; - local = st._is_local[Zp][ss]; - perm = st._permute[Zp][ss]; - ptype = st._permute_type[Zp]; - if ( local && perm ) { - spProjZp(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Zp,sF); + if ( SE->_is_local && SE->_permute ) { + spProjZp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjZp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjZp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Zp); + Impl::multLink(Uchi,U._odata[sU],chi,Zp,SE); accumReconZp(result,Uchi); // Tp - offset = st._offsets [Tp][ss]; - local = st._is_local[Tp][ss]; - perm = st._permute[Tp][ss]; - ptype = st._permute_type[Tp]; - if ( local && perm ) { - spProjTp(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Tp,sF); + if ( SE->_is_local && SE->_permute ) { + spProjTp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjTp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjTp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Tp); + Impl::multLink(Uchi,U._odata[sU],chi,Tp,SE); accumReconTp(result,Uchi); // Xm - offset = st._offsets [Xm][ss]; - local = st._is_local[Xm][ss]; - perm = st._permute[Xm][ss]; - ptype = st._permute_type[Xm]; - - if ( local && perm ) { - spProjXm(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Xm,sF); + if ( SE->_is_local && SE->_permute ) { + spProjXm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjXm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjXm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Xm); + Impl::multLink(Uchi,U._odata[sU],chi,Xm,SE); accumReconXm(result,Uchi); // Ym - offset = st._offsets [Ym][ss]; - local = st._is_local[Ym][ss]; - perm = st._permute[Ym][ss]; - ptype = st._permute_type[Ym]; - - if ( local && perm ) { - spProjYm(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Ym,sF); + if ( SE->_is_local && SE->_permute ) { + spProjYm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjYm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjYm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Ym); + Impl::multLink(Uchi,U._odata[sU],chi,Ym,SE); accumReconYm(result,Uchi); // Zm - offset = st._offsets [Zm][ss]; - local = st._is_local[Zm][ss]; - perm = st._permute[Zm][ss]; - ptype = st._permute_type[Zm]; - if ( local && perm ) { - spProjZm(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Zm,sF); + if ( SE->_is_local && SE->_permute ) { + spProjZm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjZm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjZm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Zm); + Impl::multLink(Uchi,U._odata[sU],chi,Zm,SE); accumReconZm(result,Uchi); // Tm - offset = st._offsets [Tm][ss]; - local = st._is_local[Tm][ss]; - perm = st._permute[Tm][ss]; - ptype = st._permute_type[Tm]; - if ( local && perm ) { - spProjTm(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Tm,sF); + if ( SE->_is_local && SE->_permute ) { + spProjTm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjTm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjTm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Tm); + Impl::multLink(Uchi,U._odata[sU],chi,Tm,SE); accumReconTm(result,Uchi); - vstream(out._odata[ss],result*(-0.5)); + vstream(out._odata[sF],result*(-0.5)); }; template @@ -157,141 +130,114 @@ void WilsonKernels::DiracOptDhopSiteDag(CartesianStencil &st,DoubledGaugeF SiteHalfSpinor chi; SiteSpinor result; SiteHalfSpinor Uchi; - int offset,local,perm, ptype; + StencilEntry *SE; + int ptype; // Xp - int ss=sF; - offset = st._offsets [Xm][ss]; - local = st._is_local[Xm][ss]; - perm = st._permute[Xm][ss]; - - ptype = st._permute_type[Xm]; - if ( local && perm ) { - spProjXp(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Xm,sF); + if ( SE->_is_local && SE->_permute ) { + spProjXp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjXp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjXp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Xm); + Impl::multLink(Uchi,U._odata[sU],chi,Xm,SE); spReconXp(result,Uchi); // Yp - offset = st._offsets [Ym][ss]; - local = st._is_local[Ym][ss]; - perm = st._permute[Ym][ss]; - ptype = st._permute_type[Ym]; - if ( local && perm ) { - spProjYp(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Ym,sF); + if ( SE->_is_local && SE->_permute ) { + spProjYp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjYp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjYp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Ym); + Impl::multLink(Uchi,U._odata[sU],chi,Ym,SE); accumReconYp(result,Uchi); // Zp - offset = st._offsets [Zm][ss]; - local = st._is_local[Zm][ss]; - perm = st._permute[Zm][ss]; - ptype = st._permute_type[Zm]; - if ( local && perm ) { - spProjZp(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Zm,sF); + if ( SE->_is_local && SE->_permute ) { + spProjZp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjZp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjZp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Zm); + Impl::multLink(Uchi,U._odata[sU],chi,Zm,SE); accumReconZp(result,Uchi); // Tp - offset = st._offsets [Tm][ss]; - local = st._is_local[Tm][ss]; - perm = st._permute[Tm][ss]; - ptype = st._permute_type[Tm]; - if ( local && perm ) { - spProjTp(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Tm,sF); + if ( SE->_is_local && SE->_permute ) { + spProjTp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjTp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjTp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Tm); + Impl::multLink(Uchi,U._odata[sU],chi,Tm,SE); accumReconTp(result,Uchi); // Xm - offset = st._offsets [Xp][ss]; - local = st._is_local[Xp][ss]; - perm = st._permute[Xp][ss]; - ptype = st._permute_type[Xp]; - - if ( local && perm ) { - spProjXm(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Xp,sF); + if ( SE->_is_local && SE->_permute ) { + spProjXm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjXm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjXm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Xp); + Impl::multLink(Uchi,U._odata[sU],chi,Xp,SE); accumReconXm(result,Uchi); // Ym - offset = st._offsets [Yp][ss]; - local = st._is_local[Yp][ss]; - perm = st._permute[Yp][ss]; - ptype = st._permute_type[Yp]; - - if ( local && perm ) { - spProjYm(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Yp,sF); + if ( SE->_is_local && SE->_permute ) { + spProjYm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjYm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjYm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Yp); + Impl::multLink(Uchi,U._odata[sU],chi,Yp,SE); accumReconYm(result,Uchi); // Zm - offset = st._offsets [Zp][ss]; - local = st._is_local[Zp][ss]; - perm = st._permute[Zp][ss]; - ptype = st._permute_type[Zp]; - if ( local && perm ) { - spProjZm(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Zp,sF); + if ( SE->_is_local && SE->_permute ) { + spProjZm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjZm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjZm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Zp); + Impl::multLink(Uchi,U._odata[sU],chi,Zp,SE); accumReconZm(result,Uchi); // Tm - offset = st._offsets [Tp][ss]; - local = st._is_local[Tp][ss]; - perm = st._permute[Tp][ss]; - ptype = st._permute_type[Tp]; - if ( local && perm ) { - spProjTm(tmp,in._odata[offset]); + SE=st.GetEntry(ptype,Tp,sF); + if ( SE->_is_local && SE->_permute ) { + spProjTm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjTm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjTm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,Tp); + Impl::multLink(Uchi,U._odata[sU],chi,Tp,SE); accumReconTm(result,Uchi); - vstream(out._odata[ss],result*(-0.5)); + vstream(out._odata[sF],result*(-0.5)); } template @@ -303,127 +249,124 @@ void WilsonKernels::DiracOptDhopDir(CartesianStencil &st,DoubledGaugeField SiteHalfSpinor chi; SiteSpinor result; SiteHalfSpinor Uchi; - int offset,local,perm, ptype; - int ss=sF; + StencilEntry *SE; + int ptype; - offset = st._offsets [dir][ss]; - local = st._is_local[dir][ss]; - perm = st._permute[dir][ss]; - ptype = st._permute_type[dir]; + SE=st.GetEntry(ptype,dir,sF); // Xp if(gamma==Xp){ - if ( local && perm ) { - spProjXp(tmp,in._odata[offset]); + if ( SE->_is_local && SE->_permute ) { + spProjXp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjXp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjXp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,dir); + Impl::multLink(Uchi,U._odata[sU],chi,dir,SE); spReconXp(result,Uchi); } // Yp if ( gamma==Yp ){ - if ( local && perm ) { - spProjYp(tmp,in._odata[offset]); + if ( SE->_is_local && SE->_permute ) { + spProjYp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjYp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjYp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,dir); + Impl::multLink(Uchi,U._odata[sU],chi,dir,SE); spReconYp(result,Uchi); } // Zp if ( gamma ==Zp ){ - if ( local && perm ) { - spProjZp(tmp,in._odata[offset]); + if ( SE->_is_local && SE->_permute ) { + spProjZp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjZp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjZp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,dir); + Impl::multLink(Uchi,U._odata[sU],chi,dir,SE); spReconZp(result,Uchi); } // Tp if ( gamma ==Tp ){ - if ( local && perm ) { - spProjTp(tmp,in._odata[offset]); + if ( SE->_is_local && SE->_permute ) { + spProjTp(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjTp(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjTp(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,dir); + Impl::multLink(Uchi,U._odata[sU],chi,dir,SE); spReconTp(result,Uchi); } // Xm if ( gamma==Xm ){ - if ( local && perm ) { - spProjXm(tmp,in._odata[offset]); + if ( SE->_is_local && SE->_permute ) { + spProjXm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjXm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjXm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,dir); + Impl::multLink(Uchi,U._odata[sU],chi,dir,SE); spReconXm(result,Uchi); } // Ym if ( gamma == Ym ){ - if ( local && perm ) { - spProjYm(tmp,in._odata[offset]); + if ( SE->_is_local && SE->_permute ) { + spProjYm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjYm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjYm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,dir); + Impl::multLink(Uchi,U._odata[sU],chi,dir,SE); spReconYm(result,Uchi); } // Zm if ( gamma == Zm ){ - if ( local && perm ) { - spProjZm(tmp,in._odata[offset]); + if ( SE->_is_local && SE->_permute ) { + spProjZm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjZm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjZm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,dir); + Impl::multLink(Uchi,U._odata[sU],chi,dir,SE); spReconZm(result,Uchi); } // Tm if ( gamma==Tm ) { - if ( local && perm ) { - spProjTm(tmp,in._odata[offset]); + if ( SE->_is_local && SE->_permute ) { + spProjTm(tmp,in._odata[SE->_offset]); permute(chi,tmp,ptype); - } else if ( local ) { - spProjTm(chi,in._odata[offset]); + } else if ( SE->_is_local ) { + spProjTm(chi,in._odata[SE->_offset]); } else { - chi=buf[offset]; + chi=buf[SE->_offset]; } - Impl::multLink(Uchi,U._odata[sU],chi,dir); + Impl::multLink(Uchi,U._odata[sU],chi,dir,SE); spReconTm(result,Uchi); } - vstream(out._odata[ss],result*(-0.5)); + vstream(out._odata[sF],result*(-0.5)); } FermOpTemplateInstantiate(WilsonKernels); diff --git a/lib/stencil/Stencil_common.cc b/lib/stencil/Stencil_common.cc index 7f894faf..8a8b3a54 100644 --- a/lib/stencil/Stencil_common.cc +++ b/lib/stencil/Stencil_common.cc @@ -8,11 +8,7 @@ namespace Grid { int checkerboard, const std::vector &directions, const std::vector &distances) - : _offsets(npoints), - _is_local(npoints), - _comm_buf_size(npoints), - _permute_type(npoints), - _permute(npoints) + : _entries(npoints), _permute_type(npoints) { _npoints = npoints; _grid = grid; @@ -27,9 +23,7 @@ namespace Grid { int point = i; - _offsets[i].resize( osites); - _is_local[i].resize(osites); - _permute[i].resize( osites); + _entries[i].resize( osites); int dimension = directions[i]; int displacement = distances[i]; @@ -76,7 +70,7 @@ namespace Grid { } - void CartesianStencil::Local (int point, int dimension,int shift,int cbmask) + void CartesianStencil::Local (int point, int dimension,int shiftpm,int cbmask) { int fd = _grid->_fdimensions[dimension]; int rd = _grid->_rdimensions[dimension]; @@ -84,7 +78,7 @@ namespace Grid { int gd = _grid->_gdimensions[dimension]; // Map to always positive shift modulo global full dimension. - shift = (shift+fd)%fd; + int shift = (shiftpm+fd)%fd; // the permute type int permute_dim =_grid->PermuteDim(dimension); @@ -98,6 +92,15 @@ namespace Grid { int sshift = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,cb); int sx = (x+sshift)%rd; + + int wraparound=0; + if ( (shiftpm==-1) && (sx>x) ) { + wraparound = 1; + } + if ( (shiftpm== 1) && (sx=0); assert(shiftx) && (grid->_processor_coor[dimension]==0) ) { + wraparound = 1; + } + if ( (shiftpm== 1) && (sx_processor_coor[dimension]==grid->_processors[dimension]-1) ) { + wraparound = 1; + } if (!offnode) { int permute_slice=0; - CopyPlane(point,dimension,x,sx,cbmask,permute_slice); + CopyPlane(point,dimension,x,sx,cbmask,permute_slice,wraparound); } else { @@ -161,13 +171,13 @@ namespace Grid { int unified_buffer_offset = _unified_buffer_size; _unified_buffer_size += words; - ScatterPlane(point,dimension,x,cbmask,unified_buffer_offset); // permute/extract/merge is done in comms phase + ScatterPlane(point,dimension,x,cbmask,unified_buffer_offset,wraparound); // permute/extract/merge is done in comms phase } } } // Routine builds up integer table for each site in _offsets, _is_local, _permute - void CartesianStencil::CopyPlane(int point, int dimension,int lplane,int rplane,int cbmask,int permute) + void CartesianStencil::CopyPlane(int point, int dimension,int lplane,int rplane,int cbmask,int permute,int wrap) { int rd = _grid->_rdimensions[dimension]; @@ -180,9 +190,10 @@ namespace Grid { // Simple block stride gather of SIMD objects for(int n=0;n<_grid->_slice_nblock[dimension];n++){ for(int b=0;b<_grid->_slice_block[dimension];b++){ - _offsets [point][lo+o+b]=ro+o+b; - _is_local[point][lo+o+b]=1; - _permute [point][lo+o+b]=permute; + _entries[point][lo+o+b]._offset =ro+o+b; + _entries[point][lo+o+b]._is_local=1; + _entries[point][lo+o+b]._permute=permute; + _entries[point][lo+o+b]._around_the_world=wrap; } o +=_grid->_slice_stride[dimension]; } @@ -199,9 +210,10 @@ namespace Grid { int ocb=1<<_grid->CheckerBoardFromOindex(o+b); if ( ocb&cbmask ) { - _offsets [point][lo+o+b]=ro+o+b; - _is_local[point][lo+o+b]=1; - _permute [point][lo+o+b]=permute; + _entries[point][lo+o+b]._offset =ro+o+b; + _entries[point][lo+o+b]._is_local=1; + _entries[point][lo+o+b]._permute=permute; + _entries[point][lo+o+b]._around_the_world=wrap; } } @@ -211,7 +223,7 @@ namespace Grid { } } // Routine builds up integer table for each site in _offsets, _is_local, _permute - void CartesianStencil::ScatterPlane (int point,int dimension,int plane,int cbmask,int offset) + void CartesianStencil::ScatterPlane (int point,int dimension,int plane,int cbmask,int offset, int wrap) { int rd = _grid->_rdimensions[dimension]; @@ -224,9 +236,10 @@ namespace Grid { // Simple block stride gather of SIMD objects for(int n=0;n<_grid->_slice_nblock[dimension];n++){ for(int b=0;b<_grid->_slice_block[dimension];b++){ - _offsets [point][so+o+b]=offset+(bo++); - _is_local[point][so+o+b]=0; - _permute [point][so+o+b]=0; + _entries[point][so+o+b]._offset =offset+(bo++); + _entries[point][so+o+b]._is_local=0; + _entries[point][so+o+b]._permute=0; + _entries[point][so+o+b]._around_the_world=wrap; } o +=_grid->_slice_stride[dimension]; } @@ -242,9 +255,10 @@ namespace Grid { int ocb=1<<_grid->CheckerBoardFromOindex(o+b);// Could easily be a table lookup if ( ocb & cbmask ) { - _offsets [point][so+o+b]=offset+(bo++); - _is_local[point][so+o+b]=0; - _permute [point][so+o+b]=0; + _entries[point][so+o+b]._offset =offset+(bo++); + _entries[point][so+o+b]._is_local=0; + _entries[point][so+o+b]._permute =0; + _entries[point][so+o+b]._around_the_world=wrap; } } o +=_grid->_slice_stride[dimension]; diff --git a/tests/Make.inc b/tests/Make.inc index 111eed35..c765cdf7 100644 --- a/tests/Make.inc +++ b/tests/Make.inc @@ -1,5 +1,5 @@ -bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_gparity Test_hmc_EOWilsonFermionGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi +bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_gparity Test_gpwilson_even_odd Test_hmc_EOWilsonFermionGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_GaugeAction_SOURCES=Test_GaugeAction.cc @@ -94,6 +94,10 @@ Test_gparity_SOURCES=Test_gparity.cc Test_gparity_LDADD=-lGrid +Test_gpwilson_even_odd_SOURCES=Test_gpwilson_even_odd.cc +Test_gpwilson_even_odd_LDADD=-lGrid + + Test_hmc_EOWilsonFermionGauge_SOURCES=Test_hmc_EOWilsonFermionGauge.cc Test_hmc_EOWilsonFermionGauge_LDADD=-lGrid diff --git a/tests/Test_gparity.cc b/tests/Test_gparity.cc index 46c019c6..52bf5e27 100644 --- a/tests/Test_gparity.cc +++ b/tests/Test_gparity.cc @@ -4,26 +4,54 @@ using namespace std; using namespace Grid; using namespace Grid::QCD; -template -struct scal { - d internal; -}; +typedef typename GparityDomainWallFermionR::FermionField FermionField; - Gamma::GammaMatrix Gmu [] = { - Gamma::GammaX, - Gamma::GammaY, - Gamma::GammaZ, - Gamma::GammaT - }; +template +void Replicate(Lattice &coarse,Lattice & fine) +{ + typedef typename vobj::scalar_object sobj; + + GridBase *cg = coarse._grid; + GridBase *fg = fine._grid; + + int nd = cg->_ndimension; + + subdivides(cg,fg); + + assert(cg->_ndimension==fg->_ndimension); + + std::vector ratio(cg->_ndimension); + + for(int d=0;d_ndimension;d++){ + ratio[d] = fg->_fdimensions[d]/cg->_fdimensions[d]; + } + + std::vector fcoor(nd); + std::vector ccoor(nd); + for(int g=0;ggSites();g++){ + + fg->GlobalIndexToGlobalCoor(g,fcoor); + for(int d=0;d_gdimensions[d]; + } + + sobj tmp; + peekSite(tmp,coarse,ccoor); + pokeSite(tmp,fine,fcoor); + } + +} int main (int argc, char ** argv) { + const int nu = 0; + Grid_init(&argc,&argv); const int Ls=4; - const int L=4; + const int L =4; std::vector latt_2f(Nd,L); - std::vector latt_1f(Nd,L); latt_1f[0] = 2L; + std::vector latt_1f(Nd,L); latt_1f[nu] = 2*L; std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); std::vector mpi_layout = GridDefaultMpi(); //node layout @@ -33,30 +61,64 @@ int main (int argc, char ** argv) GridCartesian * FGrid_1f = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_1f); GridRedBlackCartesian * FrbGrid_1f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_1f); + + GridCartesian * UGrid_2f = SpaceTimeGrid::makeFourDimGrid(latt_2f, simd_layout, mpi_layout); + GridRedBlackCartesian * UrbGrid_2f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_2f); + GridCartesian * FGrid_2f = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_2f); + GridRedBlackCartesian * FrbGrid_2f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_2f); + std::vector seeds4({1,2,3,4}); std::vector seeds5({5,6,7,8}); - GridParallelRNG RNG5_1f(FGrid_1f); RNG5_1f.SeedFixedIntegers(seeds5); - GridParallelRNG RNG4_1f(UGrid_1f); RNG4_1f.SeedFixedIntegers(seeds4); + GridParallelRNG RNG5_2f(FGrid_2f); RNG5_2f.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4_2f(UGrid_2f); RNG4_2f.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu_2f(UGrid_2f); + SU3::HotConfiguration(RNG4_2f,Umu_2f); + + LatticeFermion src (FGrid_2f); + LatticeFermion tmpsrc(FGrid_2f); + FermionField src_2f(FGrid_2f); + LatticeFermion src_1f(FGrid_1f); + + // Replicate fermion source + random(RNG5_2f,src); + PokeIndex<0>(src_2f,src,0); + tmpsrc=src*2.0; + PokeIndex<0>(src_2f,tmpsrc,1); - LatticeFermion src_1f(FGrid_1f); random(RNG5_1f,src_1f); LatticeFermion result_1f(FGrid_1f); result_1f=zero; - LatticeGaugeField Umu_1f(UGrid_1f); random(RNG4_1f,Umu_1f); + LatticeGaugeField Umu_1f(UGrid_1f); + Replicate(Umu_2f,Umu_1f); //Coordinate grid for reference LatticeInteger xcoor_1f(UGrid_1f); - LatticeCoordinate(xcoor_1f,0); + LatticeCoordinate(xcoor_1f,nu); //Copy-conjugate the gauge field //First C-shift the lattice by Lx/2 { - LatticeGaugeField Umu_shift = conjugate( Cshift(Umu_1f,0,L) ); + LatticeGaugeField Umu_shift = conjugate( Cshift(Umu_1f,nu,L) ); Umu_1f = where( xcoor_1f >= Integer(L), Umu_shift, Umu_1f ); + + // hack test to check the same + Replicate(Umu_2f,Umu_shift); + Umu_shift=Umu_shift-Umu_1f; + cout << GridLogMessage << "Umu diff " << norm2(Umu_shift)<(Umu_1f,nu); + Unu = where(xcoor_1f == Integer(2*L-1), -Unu, Unu); + PokeIndex(Umu_1f,Unu,nu); } - //Make the gauge field antiperiodic in x-direction - Umu_1f = where(xcoor_1f == Integer(L-1), -Umu_1f, Umu_1f); - - RealD mass=0.1; + //Coordinate grid for reference + LatticeInteger xcoor_1f5(FGrid_1f); + LatticeCoordinate(xcoor_1f5,1+nu); + Replicate(src,src_1f); + src_1f = where( xcoor_1f5 >= Integer(L), 2.0*src_1f,src_1f ); + + RealD mass=0.0; RealD M5=1.8; DomainWallFermionR Ddwf(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5); @@ -69,5 +131,136 @@ int main (int argc, char ** argv) ConjugateGradient CG(1.0e-8,10000); CG(HermOpEO,src_o_1f,result_o_1f); + GparityDomainWallFermionR GPDdwf(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f,mass,M5); + + for(int disp=-1;disp<=1;disp+=2) + for(int mu=0;mu<5;mu++) + { + FermionField Dsrc_2f(FGrid_2f); + + LatticeFermion Dsrc_1f(FGrid_1f); + LatticeFermion Dsrc_2freplica(FGrid_1f); + LatticeFermion Dsrc_2freplica0(FGrid_1f); + LatticeFermion Dsrc_2freplica1(FGrid_1f); + + if ( mu ==0 ) { + std::cout << GridLogMessage<< " Cross checking entire hopping term"<(Dsrc_2f,0); + LatticeFermion Dsrc_2f1(FGrid_2f); Dsrc_2f1 = PeekIndex<0>(Dsrc_2f,1); + + // Dsrc_2f1 = Dsrc_2f1 - Dsrc_2f0; + // std::cout << GridLogMessage << " Cross check two halves " <= Integer(L), Dsrc_2freplica1,Dsrc_2freplica0 ); + Dsrc_2freplica = Dsrc_2freplica - Dsrc_1f ; + std::cout << GridLogMessage << " Cross check against doubled latt " < CG2f(1.0e-8,10000); + SchurDiagMooeeOperator HermOpEO2f(GPDdwf); + CG2f(HermOpEO2f,src_o_2f,result_o_2f); + + std::cout << "2f cb "<(result_o_2f,0); + res1o = PeekIndex<0>(result_o_2f,1); + + std::cout << "res cb "<= Integer(L), replica1,replica0 ); + + replica0 = zero; + setCheckerboard(replica0,result_o_1f); + + std::cout << "Norm2 solutions is " <oSites();i++){ - - int offset = myStencil._offsets [0][i]; - int local = myStencil._is_local[0][i]; - int permute_type = myStencil._permute_type[0]; - int perm =myStencil._permute[0][i]; - if ( local && perm ) - permute(Check._odata[i],Foo._odata[offset],permute_type); - else if (local) - Check._odata[i] = Foo._odata[offset]; + + int permute_type; + StencilEntry *SE; + SE = myStencil.GetEntry(permute_type,0,i); + + if ( SE->_is_local && SE->_permute ) + permute(Check._odata[i],Foo._odata[SE->_offset],permute_type); + else if (SE->_is_local) + Check._odata[i] = Foo._odata[SE->_offset]; else - Check._odata[i] = comm_buf[offset]; + Check._odata[i] = comm_buf[SE->_offset]; } Real nrmC = norm2(Check);