From c6e88d9a110499bb5ed318829dccfa24808d6ba2 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Tue, 11 Nov 2025 08:36:57 -0500 Subject: [PATCH] T-direction terms done --- Grid/stencil/IcosahedralStencil.h | 9 ++- tests/debug/Test_icosahedron.cc | 101 +++++++++++++++++++++++++----- 2 files changed, 91 insertions(+), 19 deletions(-) diff --git a/Grid/stencil/IcosahedralStencil.h b/Grid/stencil/IcosahedralStencil.h index 4df63960..d67c6e3f 100644 --- a/Grid/stencil/IcosahedralStencil.h +++ b/Grid/stencil/IcosahedralStencil.h @@ -809,6 +809,7 @@ public: SE._offset = grid->oIndex(TpCoor); SE._permute = permuteTp; acceleratorPut(this->_entries[lexTp],SE); + SE._offset = grid->oIndex(TmCoor); SE._permute = permuteTm; acceleratorPut(this->_entries[lexTm],SE); @@ -846,9 +847,8 @@ public: int permute; int64_t nbr_pole_site; - tCoor = Coor; - lex = pole_site*np+6;// tp + tCoor = Coor; tCoor[tdim] = periAdd(tCoor[tdim],1,Rt,permute); nbr_pole_site = vertexgrid->PoleSiteForOcoor(tCoor); SE._offset = nbr_pole_site; @@ -856,7 +856,9 @@ public: acceleratorPut(this->_entries[lex],SE); lex = pole_site*np+7;// tm + tCoor = Coor; tCoor[tdim] = periAdd(tCoor[tdim],-1,Rt,permute); + // std::cout << " pole_site "<PoleSiteForOcoor(tCoor); SE._offset = nbr_pole_site; SE._permute= permute; @@ -893,9 +895,9 @@ public: int permute; int64_t nbr_pole_site; - tCoor = Coor; lex = pole_site*np+6;// tp + tCoor = Coor; tCoor[tdim] = periAdd(tCoor[tdim],1,Rt,permute); nbr_pole_site = vertexgrid->PoleSiteForOcoor(tCoor); SE._offset = nbr_pole_site; @@ -904,6 +906,7 @@ public: // std::cout << " Put nbr "<PoleSiteForOcoor(tCoor); SE._offset = nbr_pole_site; diff --git a/tests/debug/Test_icosahedron.cc b/tests/debug/Test_icosahedron.cc index 046f6879..3ff9cd6d 100644 --- a/tests/debug/Test_icosahedron.cc +++ b/tests/debug/Test_icosahedron.cc @@ -282,6 +282,8 @@ public: const int ent_Xm = 3; const int ent_Ym = 4; const int ent_Dm = 5; + const int ent_Tp = 6; + const int ent_Tm = 7; accelerator_for(ss,VertexGrid->oSites(),vComplex::Nsimd(),{ @@ -303,6 +305,14 @@ public: SE = stencil_v.GetEntry(ent_Dm,ss); uint64_t dm_idx = SE->_offset; int missingLink = SE->_missing_link; + + SE = stencil_v.GetEntry(ent_Tm,ss); + uint64_t tm_idx = SE->_offset; + int permuteTm = SE->_permute; + + SE = stencil_v.GetEntry(ent_Tp,ss); + uint64_t tp_idx = SE->_offset; + int permuteTp = SE->_permute; auto i = in_v(ss)(); auto inxp = in_v(xp_idx)(); @@ -312,16 +322,22 @@ public: auto inym = in_v(ym_idx)(); auto indm = in_v(dm_idx)(); - auto o = out_v(ss)(); + auto intm = in_v(tm_idx)(); + auto intp = in_v(tp_idx)(); + const int Nsimdm1 = vComplex::Nsimd()-1; + if ( permuteTp ) rotate(intp,intp,1); + if ( permuteTm ) rotate(intm,intm,Nsimdm1); + + auto o = i; if ( missingLink ) { - o = (1.0/5.0)*(inxp+inyp+indp+inxm+inym)-o; + o = (1.0/5.0)*(inxp+inyp+indp+inxm+inym)-i; } else { - o = (1.0/6.0)*(inxp+inyp+indp+inxm+inym+indm)-o; + o = (1.0/6.0)*(inxp+inyp+indp+inxm+inym+indm)-i; } - + o = o + (1.0/2.0)*(intp+intm)-i; // Coefficients need to be supplied; relative time/space weight coalescedWrite(out_v[ss](),o); - }); + }); } /* * Temporarily compute all without summing p and m directions, and XYD for the temporal links @@ -887,8 +903,25 @@ public: indp = U_v(ss)(2)*indp; inxm = U_v(ss)(3)*inxm; inym = U_v(ss)(4)*inym; + + const int Nsimdm1 = vComplex::Nsimd()-1; + + SE = stencil_v.GetEntry(ent_Tp,ss); + uint64_t tp_idx = SE->_offset; + auto permuteTp = SE->_permute; + + SE = stencil_v.GetEntry(ent_Tm,ss); + uint64_t tm_idx = SE->_offset; + auto permuteTm = SE->_permute; - auto o = i; + auto intm = in_v(tm_idx)(); + auto intp = in_v(tp_idx)(); + if ( permuteTp ) rotate(intp,intp,1); + if ( permuteTm ) rotate(intm,intm,Nsimdm1); + intp = U_v(ss)(6)*intp; + intm = U_v(ss)(7)*intm; + + auto o = i; if ( missingLink ) { o = (1.0/5.0)*(inxp+inyp+indp+inxm+inym)-i; @@ -896,17 +929,19 @@ public: indm = U_v(ss)(5)*indm; o = (1.0/6.0)*(inxp+inyp+indp+inxm+inym+indm)-i; } + o = o + (1.0/2.0)*(intp+intm)-i; // Coefficients need to be supplied; relative time/space weight coalescedWrite(out_v[ss](),o); }); } - void DoubleStore(GaugeField &U,DoubledGaugeField &Uds) + void DoubleStore(GaugeField &U,GaugeLinkField &Ut,DoubledGaugeField &Uds) { assert(U.Grid()==EdgeGrid); assert(Uds.Grid()==VertexGrid); autoView(Uds_v,Uds,AcceleratorWrite); autoView(U_v,U,AcceleratorRead); + autoView(Ut_v,Ut,AcceleratorRead); autoView(stencil_v,NNvv,AcceleratorRead); // Vertex result @@ -984,8 +1019,20 @@ public: coalescedWrite(Uds_v[ss](5),Link); } } - coalescedWrite(Uds_v[ss](6),Link); - coalescedWrite(Uds_v[ss](7),Link); + { + auto Lt = Ut_v(ss)(); + coalescedWrite(Uds_v[ss](6),Lt); + } + { + const int Nsimdm1 = vComplex::Nsimd()-1; + auto SE = stencil_v.GetEntry(ent_Tm,ss); + auto s = SE->_offset; + auto p = SE->_permute; + auto Lt_at_tm = Ut_v(s)(); + if(p) rotate(Lt_at_tm,Lt_at_tm,Nsimdm1); + Lt_at_tm = adj(Lt_at_tm); + coalescedWrite(Uds_v[ss](7),Lt_at_tm); + } }); auto pole_sites = VertexGrid->oSites() - VertexGrid->CartesianOsites(); auto pole_offset= VertexGrid->CartesianOsites(); @@ -1001,8 +1048,20 @@ public: auto Link = adj(U_v(0)(0)); Link=Zero(); coalescedWrite(Uds_v[pole_offset+ss](5),Link); - coalescedWrite(Uds_v[pole_offset+ss](6),Link); - coalescedWrite(Uds_v[pole_offset+ss](7),Link); + { + auto Lt = Ut_v(pole_offset+ss)(); + coalescedWrite(Uds_v[pole_offset+ss](6),Lt); + } + { + const int Nsimdm1 = vComplex::Nsimd()-1; + auto SE = stencil_v.GetEntry(ent_Tm,pole_offset+ss); + auto s = SE->_offset; + auto p = SE->_permute; + auto Link = Ut_v(s)(); + if ( p ) rotate(Link,Link,Nsimdm1); + Link = adj(Link); + coalescedWrite(Uds_v[pole_offset+ss](7),Link); + } }); } void GaugeTransform(GaugeLinkField >, GaugeField &Umu, GaugeLinkField &Ut) @@ -1304,8 +1363,11 @@ int main (int argc, char ** argv) std::cout << GridLogMessage<< "Calling double storing gauge field" <