diff --git a/Grid/cartesian/CartesianCrossIcosahedron.h b/Grid/cartesian/CartesianCrossIcosahedron.h index 1be63d2e..776b8265 100644 --- a/Grid/cartesian/CartesianCrossIcosahedron.h +++ b/Grid/cartesian/CartesianCrossIcosahedron.h @@ -41,6 +41,12 @@ enum NorthSouth { North = 1, South = 0 }; +enum IcoshedralDirections { + IcosahedronPatchX = 0, + IcosahedronPatchY = 1, + IcosahedronPatchDiagonal=2, + NumIcosahedralPolarizations +}; const int IcosahedralPatches = 10; const int HemiPatches=IcosahedralPatches/2; diff --git a/Grid/stencil/IcosahedralStencil.h b/Grid/stencil/IcosahedralStencil.h index 1fd3b36f..4df63960 100644 --- a/Grid/stencil/IcosahedralStencil.h +++ b/Grid/stencil/IcosahedralStencil.h @@ -35,16 +35,17 @@ struct IcosahedralStencilEntry { uint8_t _adjoint; // is this with our lattice array (else in a comms buf) uint8_t _polarisation; // which lorentz index on the neighbours patch uint8_t _missing_link; // + uint8_t _permute; // did this wrap (in T-direction) }; -enum IcoshedralDirections { - IcosahedronPatchX = 0, - IcosahedronPatchY = 1, - IcosahedronPatchDiagonal=2, - IcosahedronTime=3 -}; inline int periAdd(int A,int inc,int L) { return (A+inc+L)%L ; } +inline int periAdd(int A,int inc,int L,int & wrap) { + int r = (A+inc+L)%L; + if ( r != (A+inc) ) wrap = 1; + else wrap =0; + return r; +} class IcosahedralStencilView { public: @@ -86,6 +87,7 @@ public: // If needing edge mesh "neigbours" to assemble loops we must find the mapping of a forward link // to a corresponding "backward link" on the pole deviceVector _entries; + void GetNbrForPlusX(GridBase *grid,Coordinate &Coor,Coordinate &NbrCoor, int &isPole) { @@ -559,16 +561,6 @@ public: assert(dmxy_count + dmxy_count_special + num_missing == triangle_ref); assert(dmyx_count + dmyx_count_special + num_missing == triangle_ref); - std::cout << GridLogMessage<< "------------------------------------"<_grid = grid; this->_vertexgrid = vertexgrid; + assert(grid->ProcessorCount() ==1); // Loop over L^2 x T x npatch and the assert(grid->isIcosahedral()); assert(grid->isIcosahedral()); @@ -598,16 +591,21 @@ public: GridBase * vertexgrid = this->_vertexgrid; uint64_t cart_sites = grid->CartesianOsites(); - uint64_t Npole_sites = grid->NorthPoleOsites(); - uint64_t Spole_sites = grid->SouthPoleOsites(); + uint64_t Npole_sites = vertexgrid->NorthPoleOsites(); + uint64_t Spole_sites = vertexgrid->SouthPoleOsites(); + Coordinate pcoor = grid->ThisProcessorCoor(); Coordinate pgrid = grid->ProcessorGrid(); /* * resize the stencil entries array and set npoints */ - const int np=6; + const int np=8; this->_npoints=np; // Move to template param? - this->_entries.resize(this->_npoints * cart_sites); + if ( vertexOutputs ) { + this->_entries.resize(this->_npoints * (cart_sites+Npole_sites+Spole_sites)); + } else { + this->_entries.resize(this->_npoints * cart_sites); + } this->_entries_p = &_entries[0]; int nd = grid->Nd(); @@ -618,13 +616,14 @@ public: Coordinate Coor; Coordinate NbrCoor; - Integer lexXp = site*np ; Integer lexYp = site*np+1; Integer lexDp = site*np+2; Integer lexXm = site*np+3; Integer lexYm = site*np+4; Integer lexDm = site*np+5; + Integer lexTp = site*np+6; + Integer lexTm = site*np+7; IcosahedralStencilEntry SE; @@ -659,6 +658,8 @@ public: Coordinate XmCoor; Coordinate YmCoor; Coordinate DmCoor; + Coordinate TpCoor; + Coordinate TmCoor; GetNbrForPlusDiagonal(grid,Coor,DpCoor); GetNbrForPlusX(grid,Coor,XpCoor,isPoleX); @@ -668,7 +669,13 @@ public: GetNbrForMinusX(grid,Coor,XmCoor); GetNbrForMinusY(grid,Coor,YmCoor); - int DpPatch = DpCoor[nd-1]; + int tdim = 2; int delta = 1; + int permuteTp; + int permuteTm; + GetOrthogNbr(grid,Coor,TpCoor,tdim, delta,permuteTp); + GetOrthogNbr(grid,Coor,TmCoor,tdim,-delta,permuteTm); + + int DpPatch = DpCoor[nd-1]; int DpHemiPatch = DpCoor[nd-1]%HemiPatches; int DpHemisphere = DpCoor[nd-1]/HemiPatches; @@ -682,10 +689,11 @@ public: int YmHemiPatch = YmCoor[nd-1]%HemiPatches; int YmHemisphere = YmCoor[nd-1]/HemiPatches; - int DmPatch = DmCoor[nd-1]; + int DmPatch = DmCoor[nd-1]; int DmHemiPatch = DmCoor[nd-1]%HemiPatches; int DmHemisphere = DmCoor[nd-1]/HemiPatches; + SE._permute=0; if ( vertexInputs ) {// Neighbour will live on poles and peer point //////////////////////////////////////////////// // XpCoor stencil entry; consider isPole case @@ -794,6 +802,17 @@ public: } SE._missing_link = missingLink; acceleratorPut(this->_entries[lexDm],SE); + ///////////////////////////////////////////////////////////////////// + // for Tp/mCoor + ///////////////////////////////////////////////////////////////////// + SE._polarisation = 0; + 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); + SE._permute = 0; } if ( vertexOutputs ) { @@ -802,23 +821,47 @@ public: IcosahedralStencilEntry SE; for(uint64_t site=0;siteoCoorFromOindex(Coor,site); int north = Coor[ndm1]/HemiPatches; int south = 1-north; if( (Coor[0]==(L-1))&&(Coor[1]==0) &&south ) { int64_t pole_site = vertexgrid->PoleSiteForOcoor(Coor); int64_t lex = pole_site*np+(Coor[ndm1]%HemiPatches); - // std::cout << "Coor "<_entries[lex],SE); int64_t lex5 = pole_site*np+5; // We miss the backwards link SE._missing_link = true; acceleratorPut(this->_entries[lex5],SE); + + int tdim = 2; + int Rt = vertexgrid->_rdimensions[tdim]; + int permute; + int64_t nbr_pole_site; + + tCoor = Coor; + + lex = pole_site*np+6;// tp + tCoor[tdim] = periAdd(tCoor[tdim],1,Rt,permute); + nbr_pole_site = vertexgrid->PoleSiteForOcoor(tCoor); + SE._offset = nbr_pole_site; + SE._permute= permute; + acceleratorPut(this->_entries[lex],SE); + + lex = pole_site*np+7;// tm + tCoor[tdim] = periAdd(tCoor[tdim],-1,Rt,permute); + nbr_pole_site = vertexgrid->PoleSiteForOcoor(tCoor); + SE._offset = nbr_pole_site; + SE._permute= permute; + acceleratorPut(this->_entries[lex],SE); + } } } @@ -826,6 +869,7 @@ public: IcosahedralStencilEntry SE; for(uint64_t site=0;siteoCoorFromOindex(Coor,site); int north = Coor[ndm1]/HemiPatches; if( (Coor[0]==0)&&(Coor[1]==(L-1))&&north ) { @@ -837,11 +881,36 @@ public: SE._polarisation = IcosahedronPatchY; // ignored SE._adjoint = false; // ignored SE._missing_link = false; + SE._permute=0; acceleratorPut(this->_entries[lex],SE); int64_t lex5 = pole_site*np+5; // We miss the backwards link SE._missing_link = true; acceleratorPut(this->_entries[lex5],SE); + + int tdim = 2; + int Rt = vertexgrid->_rdimensions[tdim]; + int permute; + int64_t nbr_pole_site; + + tCoor = Coor; + + lex = pole_site*np+6;// tp + tCoor[tdim] = periAdd(tCoor[tdim],1,Rt,permute); + nbr_pole_site = vertexgrid->PoleSiteForOcoor(tCoor); + SE._offset = nbr_pole_site; + SE._permute= permute; + acceleratorPut(this->_entries[lex],SE); + // std::cout << " Put nbr "<PoleSiteForOcoor(tCoor); + SE._offset = nbr_pole_site; + SE._permute= permute; + acceleratorPut(this->_entries[lex],SE); + // std::cout << " Put nbr "<Nd(); int L = grid->LocalDimensions()[0]; @@ -887,7 +955,6 @@ public: // Outer index of neighbour Offset calculation //////////////////////////////////////////////// grid->oCoorFromOindex(Coor,site); - NbrCoor = Coor; assert( grid->LocalDimensions()[1]==grid->LocalDimensions()[0]); assert( grid->_simd_layout[0]==1); // Cannot vectorise in these dims assert( grid->_simd_layout[1]==1); @@ -949,12 +1016,14 @@ public: SE._polarisation = IcosahedronPatchX; SE._adjoint = true; SE._missing_link = false; + SE._permute=0; } else { SE._offset = grid->oIndex(XpCoor); SE._is_local = true; SE._polarisation = IcosahedronPatchY; SE._adjoint = false; SE._missing_link = false; + SE._permute=0; } //////////////////////////////////////////////// @@ -973,12 +1042,14 @@ public: SE._polarisation = IcosahedronPatchY; SE._adjoint = true; SE._missing_link = false; + SE._permute=0; } else { SE._offset = grid->oIndex(YpCoor); SE._is_local = true; SE._polarisation = IcosahedronPatchX; SE._adjoint = false; SE._missing_link = false; + SE._permute=0; } //////////////////////////////////////////////// // Store in look up table @@ -986,24 +1057,459 @@ public: acceleratorPut(this->_entries[lexYX],SE); }; } + // + // Orthogonal direction support + // + // Plaquettes: + // Must be able to get the sites +T, and +X, +Y, +D + // + // Staples + // Must be able to get the sites (+T, and +X, +Y, +D) (-T, and -T+X, -T+Y, -T+D) + // + // Laplacian: + // Must be able to get the sites +-T + // + void GetOrthogNbr(GridBase *grid,Coordinate &Coor,Coordinate &NbrCoor,int dim,int delta, int & permute) + { + assert(delta==1 || delta==-1); + int L = grid->_rdimensions[dim]; + + NbrCoor = Coor; + NbrCoor[dim] = periAdd(Coor[dim],delta,L,permute); + } + void TemporalPlaquetteStencil(void) + { + GridBase * grid = this->_grid; // the edge grid + GridBase * vertexgrid = this->_vertexgrid; + + int osites = grid->oSites(); + + uint64_t cart_sites = grid->CartesianOsites(); + uint64_t Npole_sites = grid->NorthPoleOsites(); + uint64_t Spole_sites = grid->SouthPoleOsites(); + Coordinate pcoor = grid->ThisProcessorCoor(); + Coordinate pgrid = grid->ProcessorGrid(); + + /* + * resize the stencil entries array and set npoints + */ + this->_npoints=4; + this->_entries.resize(this->_npoints * cart_sites); + this->_entries_p = &_entries[0]; + + for(uint64_t site=0;siteNd(); + int L = grid->LocalDimensions()[0]; + + Coordinate Coor; + IcosahedralStencilEntry SE; + + Integer lexT = site*4+0; + Integer lexX = site*4+1; + Integer lexY = site*4+2; + Integer lexD = site*4+3; + + //////////////////////////////////////////////// + // Outer index of neighbour Offset calculation + //////////////////////////////////////////////// + grid->oCoorFromOindex(Coor,site); + assert( grid->LocalDimensions()[1]==grid->LocalDimensions()[0]); + assert( grid->_simd_layout[0]==1); // Cannot vectorise in these dims + assert( grid->_simd_layout[1]==1); + assert( grid->_processors[0]==1); // Cannot mpi distribute in these dims + assert( grid->_processors[1]==1); + + int Patch = Coor[nd-1]; + int HemiPatch = Patch%HemiPatches; + int Hemisphere= Patch/HemiPatches; + int north = Patch/HemiPatches; + int south = 1-north; + int isPoleY; + int isPoleX; + + assert(PatchoIndex(TpCoor); + SE._permute=permute; + // std::cout << " Plaq stencil "<_entries[lexT],SE); + SE._permute=0; + //////////////////////////////////////////////// + // X+ direction + //////////////////////////////////////////////// + if ( isPoleX ) { + SE._offset = vertexgrid->PoleSiteForOcoor(Coor); + } else { + SE._offset = vertexgrid->oIndex(XpCoor); + } + acceleratorPut(this->_entries[lexX],SE); + //////////////////////////////////////////////// + // Y+ direction + //////////////////////////////////////////////// + if ( isPoleY ) { + SE._offset = vertexgrid->PoleSiteForOcoor(Coor); + } else { + SE._offset = vertexgrid->oIndex(YpCoor); + } + acceleratorPut(this->_entries[lexY],SE); + //////////////////////////////////////////////// + // D+ direction + //////////////////////////////////////////////// + SE._offset = vertexgrid->oIndex(DpCoor); + // std::cout << "Coor "<_entries[lexD],SE); + }; + } + + /* + * enough to build staples in ico-T dir + */ /* * For gauge action derivative implementation * Staple * * Case1: I x T loops - * - * Need: DirP this site; no entry - * Tp @ dir++ - * DirP @ t++ - * Tp @ t-- - * DirP @ t-- - * Tp @ t--, dir++ - * * There is no complex rotation of links on other site * * Case2: I x I loops * Just use a general 6 point stencil and cherry pick terms */ + void TemporalStapleStencil(void) + { + GridBase * grid = this->_grid; // the edge grid + GridBase * vertexgrid = this->_vertexgrid; + + int osites = grid->oSites(); + + uint64_t cart_sites = grid->CartesianOsites(); + uint64_t Npole_sites = grid->NorthPoleOsites(); + uint64_t Spole_sites = grid->SouthPoleOsites(); + Coordinate pcoor = grid->ThisProcessorCoor(); + Coordinate pgrid = grid->ProcessorGrid(); + + /* + * resize the stencil entries array and set npoints + */ + this->_npoints=14; + this->_entries.resize(this->_npoints * cart_sites); + this->_entries_p = &_entries[0]; + int np= this->_npoints; + + for(uint64_t site=0;siteNd(); + int L = grid->LocalDimensions()[0]; + + Coordinate Coor; + IcosahedralStencilEntry SE; + + Integer lexTp = site*np+0; + Integer lexXp = site*np+1; + Integer lexYp = site*np+2; + Integer lexDp = site*np+3; + Integer lexTm = site*np+4; + Integer lexTmXp = site*np+5; + Integer lexTmYp = site*np+6; + Integer lexTmDp = site*np+7; + + Integer lexXm = site*np+8; + Integer lexYm = site*np+9; + Integer lexDm = site*np+10; // If !missingLink + Integer lexXmTp = site*np+11; + Integer lexYmTp = site*np+12; + Integer lexDmTp = site*np+13; // If !missingLink + + //////////////////////////////////////////////// + // Outer index of neighbour Offset calculation + //////////////////////////////////////////////// + grid->oCoorFromOindex(Coor,site); + assert( grid->LocalDimensions()[1]==grid->LocalDimensions()[0]); + assert( grid->_simd_layout[0]==1); // Cannot vectorise in these dims + assert( grid->_simd_layout[1]==1); + assert( grid->_processors[0]==1); // Cannot mpi distribute in these dims + assert( grid->_processors[1]==1); + + int Patch = Coor[nd-1]; + int HemiPatch = Coor[nd-1]%HemiPatches; + int Hemisphere = Coor[nd-1]/HemiPatches; + int north = Patch/HemiPatches; + int south = 1-north; + int isPoleY; + int isPoleX; + + assert(PatchoIndex(TpCoor); + SE._permute= permute; + acceleratorPut(this->_entries[lexTp],SE); + SE._permute=0; + //////////////////////////////////////////////// + // X+ direction + //////////////////////////////////////////////// + if ( isPoleX ) { + SE._offset = vertexgrid->PoleSiteForOcoor(Coor); + } else { + SE._offset = vertexgrid->oIndex(XpCoor); + } + acceleratorPut(this->_entries[lexXp],SE); + //////////////////////////////////////////////// + // Y+ direction + //////////////////////////////////////////////// + if ( isPoleY ) { + SE._offset = vertexgrid->PoleSiteForOcoor(Coor); + } else { + SE._offset = vertexgrid->oIndex(YpCoor); + } + acceleratorPut(this->_entries[lexYp],SE); + //////////////////////////////////////////////// + // D+ direction + //////////////////////////////////////////////// + SE._offset = vertexgrid->oIndex(DpCoor); + acceleratorPut(this->_entries[lexDp],SE); + ////////////////////////////////////////////////// + // Backward one site in time direction + ////////////////////////////////////////////////// + GetOrthogNbr(grid,Coor,TmCoor,tdim,-delta,permute); + SE._permute = permute; + SE._offset = vertexgrid->oIndex(TmCoor); + acceleratorPut(this->_entries[lexTm],SE); + //////////////////////////////////////////////// + // T-X+ hop + // T-Y+ hop + // T-D+ hop + //////////////////////////////////////////////// + GetNbrForPlusX(grid,TmCoor,TmXpCoor,isPoleX); + GetNbrForPlusY(grid,TmCoor,TmYpCoor,isPoleY); + GetNbrForPlusDiagonal(grid,TmCoor,TmDpCoor); + if ( isPoleX ) { + SE._offset = vertexgrid->PoleSiteForOcoor(TmCoor); + } else { + SE._offset = vertexgrid->oIndex(TmXpCoor); + } + acceleratorPut(this->_entries[lexTmXp],SE); + + if ( isPoleY ) { + SE._offset = vertexgrid->PoleSiteForOcoor(TmCoor); + } else { + SE._offset = vertexgrid->oIndex(TmYpCoor); + } + acceleratorPut(this->_entries[lexTmYp],SE); + + SE._offset = vertexgrid->oIndex(TmDpCoor); + acceleratorPut(this->_entries[lexTmDp],SE); + + //////////////////////////////////////////////// + // Links for the negative XYD staple on the forward T link + //////////////////////////////////////////////// + SE._permute = 0; + SE._offset = vertexgrid->oIndex(XmCoor); + SE._polarisation = IcosahedronPatchX; + if ( XmHemiPatch != HemiPatch && north ) { + SE._polarisation = IcosahedronPatchDiagonal; + } + acceleratorPut(this->_entries[lexXm],SE); + GetOrthogNbr(grid,XmCoor,XmTpCoor,tdim,delta,permute); + SE._permute = permute; + SE._offset = vertexgrid->oIndex(XmTpCoor); + acceleratorPut(this->_entries[lexXmTp],SE); + + SE._permute = 0; + SE._offset = vertexgrid->oIndex(YmCoor); + SE._polarisation = IcosahedronPatchY; + if ( YmHemiPatch != HemiPatch && south ) { + SE._polarisation = IcosahedronPatchDiagonal; + } + acceleratorPut(this->_entries[lexYm],SE); + GetOrthogNbr(grid,YmCoor,YmTpCoor,tdim,delta,permute); + SE._permute = permute; + SE._offset = vertexgrid->oIndex(YmTpCoor); + acceleratorPut(this->_entries[lexYmTp],SE); + + if ( ! missingLink ) { + + SE._polarisation = IcosahedronPatchDiagonal; // Basis rotates + if ( (DmHemiPatch != HemiPatch) && (DmHemisphere==Hemisphere) && south ) { + SE._polarisation = IcosahedronPatchX; // Basis rotates + } + if ( (DmHemiPatch != HemiPatch) && (DmHemisphere==Hemisphere) && north ) { + SE._polarisation = IcosahedronPatchY; // Basis rotates + } + + SE._permute = 0; + SE._offset = vertexgrid->oIndex(DmCoor); + acceleratorPut(this->_entries[lexDm],SE); + + GetOrthogNbr(grid,DmCoor,DmTpCoor,tdim,delta,permute); + SE._permute = permute; + SE._offset = vertexgrid->oIndex(DmTpCoor); + acceleratorPut(this->_entries[lexDmTp],SE); + + } else { + + SE._permute = 0; + SE._offset = 0; + SE._missing_link = 1; + acceleratorPut(this->_entries[lexDm],SE); + acceleratorPut(this->_entries[lexDmTp],SE); + + } + }; + + int ndm1 = grid->Nd()-1; + int L = grid->LocalDimensions()[0]; + if ( vertexgrid->ownsSouthPole() ) { + IcosahedralStencilEntry SE; + for(uint64_t site=0;siteoCoorFromOindex(Coor,site); + int north = Coor[ndm1]/HemiPatches; + int south = 1-north; + int tdim = 2; + + if( (Coor[0]==(L-1))&&(Coor[1]==0)&&south ) { + + int64_t pole_site = vertexgrid->PoleSiteForOcoor(Coor); + int64_t lex = pole_site*np+(Coor[ndm1]%HemiPatches)*2; + + SE._offset = site; + SE._is_local = true; + SE._polarisation = IcosahedronPatchX; // ignored + SE._adjoint = false; // ignored + SE._missing_link = false; + SE._permute=0; + acceleratorPut(this->_entries[lex],SE); + + lex = pole_site*np+(Coor[ndm1]%HemiPatches)*2+1; + + int Rt = vertexgrid->_rdimensions[tdim]; + int permute; + tCoor = Coor; + tCoor[tdim] = periAdd(tCoor[tdim],1,Rt,permute); + SE._offset = vertexgrid->oIndex(tCoor); + SE._permute= permute; + acceleratorPut(this->_entries[lex],SE); + + } + } + } + if ( vertexgrid->ownsNorthPole() ) { + IcosahedralStencilEntry SE; + for(uint64_t site=0;siteoCoorFromOindex(Coor,site); + int north = Coor[ndm1]/HemiPatches; + int tdim = 2; + + if( (Coor[0]==0)&&(Coor[1]==(L-1))&&north ) { + + int64_t pole_site = vertexgrid->PoleSiteForOcoor(Coor); + int64_t lex = pole_site*np+(Coor[ndm1]%HemiPatches)*2; + + SE._offset = site; + SE._is_local = true; + SE._polarisation = IcosahedronPatchY; // ignored + SE._adjoint = false; // ignored + SE._missing_link = false; + SE._permute=0; + acceleratorPut(this->_entries[lex],SE); + + lex = pole_site*np+(Coor[ndm1]%HemiPatches)*2+1; + + int Rt = vertexgrid->_rdimensions[tdim]; + int permute; + tCoor = Coor; + tCoor[tdim] = periAdd(tCoor[tdim],1,Rt,permute); + SE._offset = vertexgrid->oIndex(tCoor); + SE._permute= permute; + acceleratorPut(this->_entries[lex],SE); + + } + } + } + } }; NAMESPACE_END(Grid); diff --git a/systems/mac-arm/config-command b/systems/mac-arm/config-command index 5b75cfb5..dbb83abd 100644 --- a/systems/mac-arm/config-command +++ b/systems/mac-arm/config-command @@ -1,7 +1,7 @@ MPICXX=mpicxx CXXFLAGS=-I/opt/local/include LDFLAGS=-L/opt/local/lib/ CXX=clang++ ../../configure \ --enable-simd=GEN \ --enable-Nc=1 \ - --enable-comms=mpi-auto \ + --enable-debug \ --enable-unified=yes \ --prefix $HOME/QCD/GridInstall \ --with-lime=/Users/peterboyle/QCD/SciDAC/install/ \ diff --git a/tests/debug/Test_icosahedron.cc b/tests/debug/Test_icosahedron.cc index 3af8d588..046f6879 100644 --- a/tests/debug/Test_icosahedron.cc +++ b/tests/debug/Test_icosahedron.cc @@ -31,14 +31,29 @@ Author: Peter Boyle using namespace std; using namespace Grid; -const int MyNd=3; + +// FIXME -- the ROTATE code is CPU vector only +// Need coalescedReadPermute to support rotations + +/////////////////////////////////////////////////////////// +// 3 links per edge site in icos-plane; 3.10.L^2xLt +// 1 temporal per vertex perp to plane; (10.L^2+2)xLt +// +// This makes the HMC complicated -- the field type integrated is not simply U but a collection of fields +// and momenta with different support. Easiest to view as a "bunch of links", but breaks current code designs +// unless we have a way to split a bigger gauge field with "empty" or unused slots; probably the cleanest. +/////////////////////////////////////////////////////////// +template using iIcosahedralLorentzColourMatrix = iVector >, NumIcosahedralPolarizations> ; +typedef iIcosahedralLorentzColourMatrix IcosahedralLorentzColourMatrix; +typedef iIcosahedralLorentzColourMatrix vIcosahedralLorentzColourMatrix; +typedef Lattice LatticeIcosahedralLorentzColourMatrix; class IcosahedralGimpl { public: - typedef LatticeLorentzColourMatrix GaugeField; - typedef LatticeColourMatrix GaugeLinkField; - typedef LatticeDoubledGaugeField DoubledGaugeField; + typedef LatticeIcosahedralLorentzColourMatrix GaugeField; + typedef LatticeDoubledGaugeField DoubledGaugeField; + typedef LatticeColourMatrix GaugeLinkField; typedef LatticeComplex ComplexField; }; @@ -59,6 +74,9 @@ public: IcosahedralStencil NNev; // vertex neighbours but in edge domain IcosahedralStencil NNvv; // vertex neighbours with vertex domain IcosahedralStencil NNve; // edge neighbours with vertex domain + + IcosahedralStencil TimePlaqStencil; + IcosahedralStencil TimeStapleStencil; IcosahedralSupport(GridBase *_VertexGrid,GridBase *_EdgeGrid) : FaceStencil (_EdgeGrid,_VertexGrid), @@ -66,21 +84,25 @@ public: NNev(_EdgeGrid,_VertexGrid), NNve(_EdgeGrid,_VertexGrid), NNvv(_EdgeGrid,_VertexGrid), + TimePlaqStencil(_EdgeGrid,_VertexGrid), + TimeStapleStencil(_EdgeGrid,_VertexGrid), VertexGrid(_VertexGrid), EdgeGrid(_EdgeGrid) { FaceStencil.FaceStencil(); - std::cout << "NNee"<oSites(),vComplex::Nsimd(),{ - // Three forward links from this site + // Three forward links from this site auto Lx = Umu_v(ss)(IcosahedronPatchX); auto Ly = Umu_v(ss)(IcosahedronPatchY); auto Ld = Umu_v(ss)(IcosahedronPatchDiagonal); @@ -301,7 +323,515 @@ public: coalescedWrite(out_v[ss](),o); }); } + /* + * Temporarily compute all without summing p and m directions, and XYD for the temporal links + * Do this so we can compute and check all link x staple loops + void TemporalStaplesOld(GaugeLinkField &Ut, + GaugeField &Umu, + GaugeLinkField &stapleXTp, + GaugeLinkField &stapleYTp, + GaugeLinkField &stapleDTp, + GaugeLinkField &stapleXTm, + GaugeLinkField &stapleYTm, + GaugeLinkField &stapleDTm, + GaugeLinkField &stapleTXp, // these 5 are on vertex graph + GaugeLinkField &stapleTYp, // and contain the north south pole temporal link staples + GaugeLinkField &stapleTDp, + GaugeLinkField &stapleTXm, + GaugeLinkField &stapleTYm, + GaugeLinkField &stapleTDm) + { + + autoView(Umu_v,Umu,AcceleratorRead); + autoView(Ut_v,Ut,AcceleratorRead); + autoView(stapleXTp_v,stapleXTp,AcceleratorWrite); + autoView(stapleYTp_v,stapleYTp,AcceleratorWrite); + autoView(stapleDTp_v,stapleDTp,AcceleratorWrite); + + autoView(stapleXTm_v,stapleXTm,AcceleratorWrite); + autoView(stapleYTm_v,stapleYTm,AcceleratorWrite); + autoView(stapleDTm_v,stapleDTm,AcceleratorWrite); + + autoView(stapleTXp_v,stapleTXp,AcceleratorWrite); + autoView(stapleTYp_v,stapleTYp,AcceleratorWrite); + autoView(stapleTDp_v,stapleTDp,AcceleratorWrite); + + autoView(stapleTXm_v,stapleTXm,AcceleratorWrite); + autoView(stapleTYm_v,stapleTYm,AcceleratorWrite); + autoView(stapleTDm_v,stapleTDm,AcceleratorWrite); + + autoView(stencil_v,TimeStapleStencil,AcceleratorRead); + + const int np = TimeStapleStencil._npoints; + + Integer ent_Tp = 0; + Integer ent_Xp = 1; + Integer ent_Yp = 2; + Integer ent_Dp = 3; + Integer ent_Tm = 4; + Integer ent_TmXp = 5; + Integer ent_TmYp = 6; + Integer ent_TmDp = 7; + Integer ent_Xm = 8; + Integer ent_Ym = 9; + Integer ent_Dm = 10; + Integer ent_XmTp = 11; + Integer ent_YmTp = 12; + Integer ent_DmTp = 13; + + accelerator_for(ss,EdgeGrid->oSites(),vComplex::Nsimd(),{ + + // Four forward links from this site + auto Lx = Umu_v(ss)(IcosahedronPatchX); + auto Ly = Umu_v(ss)(IcosahedronPatchY); + auto Ld = Umu_v(ss)(IcosahedronPatchDiagonal); + auto Lt = Ut_v(ss)(); + IcosahedralStencilEntry *SE; + int64_t s; + int p; + /////////////////////////////////////////////////////////////////// + // Load all links required and multiple round loops + /////////////////////////////////////////////////////////////////// + SE = stencil_v.GetEntry(ent_Tp,ss); + s = SE->_offset; + p = SE->_permute; + auto Lx_at_tp = Umu_v(s)(IcosahedronPatchX); + auto Ly_at_tp = Umu_v(s)(IcosahedronPatchY); + auto Ld_at_tp = Umu_v(s)(IcosahedronPatchDiagonal); + if ( p ) { + #ifdef GRID_SIMT + #error GPU not supported yet + #endif + rotate(Lx_at_tp,Lx_at_tp,1); + rotate(Ly_at_tp,Ly_at_tp,1); + rotate(Ld_at_tp,Ld_at_tp,1); + } + SE = stencil_v.GetEntry(ent_Xp,ss); + s = SE->_offset; + auto Lt_at_xp = Ut_v(s)(); + + SE = stencil_v.GetEntry(ent_Yp,ss); + s = SE->_offset; + auto Lt_at_yp = Ut_v(s)(); + + SE = stencil_v.GetEntry(ent_Dp,ss); + s = SE->_offset; + auto Lt_at_dp = Ut_v(s)(); + + Coordinate Coor; + EdgeGrid->oCoorFromOindex(Coor,ss); + + // Lines start at origin and progress to top of link + // close as gauge invariant if multiplied from left by link + coalescedWrite(stapleXTp_v[ss](),Lt_at_xp*adj(Lx_at_tp)*adj(Lt)); + coalescedWrite(stapleYTp_v[ss](),Lt_at_yp*adj(Ly_at_tp)*adj(Lt)); + coalescedWrite(stapleDTp_v[ss](),Lt_at_dp*adj(Ld_at_tp)*adj(Lt)); + + SE = stencil_v.GetEntry(ent_Tm,ss); + s = SE->_offset; + p = SE->_permute; + auto Lt_at_tm = Ut_v(s)(); + auto Lx_at_tm = Umu_v(s)(IcosahedronPatchX); + auto Ly_at_tm = Umu_v(s)(IcosahedronPatchY); + auto Ld_at_tm = Umu_v(s)(IcosahedronPatchDiagonal); + const int Nsimdm1 = vComplex::Nsimd()-1; + + if ( p ) { + rotate(Lt_at_tm,Lt_at_tm,Nsimdm1); + rotate(Lx_at_tm,Lx_at_tm,Nsimdm1); + rotate(Ly_at_tm,Ly_at_tm,Nsimdm1); + rotate(Ld_at_tm,Ld_at_tm,Nsimdm1); + } + + SE = stencil_v.GetEntry(ent_TmXp,ss); + s = SE->_offset; + p = SE->_permute; + auto Lt_at_tmxp = Ut_v(s)(); + + SE = stencil_v.GetEntry(ent_TmYp,ss); + s = SE->_offset; + auto Lt_at_tmyp = Ut_v(s)(); + + SE = stencil_v.GetEntry(ent_TmDp,ss); + s = SE->_offset; + auto Lt_at_tmdp = Ut_v(s)(); + + if ( p ) { + rotate(Lt_at_tmxp,Lt_at_tmxp,Nsimdm1); + rotate(Lt_at_tmyp,Lt_at_tmyp,Nsimdm1); + rotate(Lt_at_tmdp,Lt_at_tmdp,Nsimdm1); + } + + coalescedWrite(stapleXTm_v[ss](),adj(Lt_at_tmxp)*adj(Lx_at_tm)*Lt_at_tm); + coalescedWrite(stapleYTm_v[ss](),adj(Lt_at_tmyp)*adj(Ly_at_tm)*Lt_at_tm); + coalescedWrite(stapleDTm_v[ss](),adj(Lt_at_tmdp)*adj(Ld_at_tm)*Lt_at_tm); + + coalescedWrite(stapleTXp_v[ss](),Lx_at_tp*adj(Lt_at_xp)*adj(Lx)); + coalescedWrite(stapleTYp_v[ss](),Ly_at_tp*adj(Lt_at_yp)*adj(Ly)); + coalescedWrite(stapleTDp_v[ss](),Ld_at_tp*adj(Lt_at_dp)*adj(Ld)); + + int pol; + SE = stencil_v.GetEntry(ent_Xm,ss); + s = SE->_offset; + pol= SE->_polarisation; + auto Lt_at_xm = Ut_v(s)(); + auto Lx_at_xm = Umu_v(s)(pol); + + SE = stencil_v.GetEntry(ent_XmTp,ss); + s = SE->_offset; + p = SE->_permute; + pol= SE->_polarisation; + auto Lx_at_xmtp = Umu_v(s)(pol); + + SE = stencil_v.GetEntry(ent_Ym,ss); + s = SE->_offset; + pol= SE->_polarisation; + auto Lt_at_ym = Ut_v(s)(); + auto Ly_at_ym = Umu_v(s)(pol); + + SE = stencil_v.GetEntry(ent_YmTp,ss); + s = SE->_offset; + p = SE->_permute; + pol= SE->_polarisation; + auto Ly_at_ymtp = Umu_v(s)(pol); + + if ( p ) { + rotate(Ly_at_ymtp,Ly_at_ymtp,1); + rotate(Lx_at_xmtp,Lx_at_xmtp,1); + } + + coalescedWrite(stapleTXm_v[ss](),adj(Lx_at_xmtp)*adj(Lt_at_xm)*Lx_at_xm); + coalescedWrite(stapleTYm_v[ss](),adj(Ly_at_ymtp)*adj(Lt_at_ym)*Ly_at_ym); + + SE = stencil_v.GetEntry(ent_Dm,ss); + int missingLink = SE->_missing_link; + s = SE->_offset; + pol= SE->_polarisation; + if ( missingLink ) { + auto zero = Lx; + zero = Zero(); + coalescedWrite(stapleTDm_v[ss](),zero); + } else { + auto Lt_at_dm = Ut_v(s)(); + auto Ld_at_dm = Umu_v(s)(pol); + SE = stencil_v.GetEntry(ent_DmTp,ss); + s = SE->_offset; + auto Ld_at_dmtp = Umu_v(s)(pol); + if ( p ) { + rotate(Ld_at_dmtp,Ld_at_dmtp,1); + } + coalescedWrite(stapleTDm_v[ss](),adj(Ld_at_dmtp)*adj(Lt_at_dm)*Ld_at_dm); + } + + }); + + const int ent_h0 = 0; + const int ent_h0tp = 1; + const int ent_h1 = 2; + const int ent_h1tp = 3; + const int ent_h2 = 4; + const int ent_h2tp = 5; + const int ent_h3 = 6; + const int ent_h3tp = 7; + const int ent_h4 = 8; + const int ent_h4tp = 9; + + uint64_t cart_sites = EdgeGrid->oSites(); + uint64_t pole_sites = VertexGrid->NorthPoleOsites() + + VertexGrid->SouthPoleOsites(); + +#define LINK(ENT_0,ENT_1) \ + {SE = stencil_v.GetEntry(ENT_0,ss);\ + s = SE->_offset; \ + pol = SE->_polarisation; \ + \ + auto Lt = Ut_v(ss)(); \ + auto Li = Umu_v(s)(pol); \ + auto Lth = Ut_v(s)(); \ + \ + SE=stencil_v.GetEntry(ENT_1,ss); \ + s = SE->_offset; \ + p = SE->_permute; \ + auto Litp = Umu_v(s)(pol); \ + if ( p ) { \ + rotate(Litp,Litp,1); \ + } \ + \ + auto stap = adj(Litp)*adj(Lth)*Li; \ + std::cout << ss<<" "<oSites(),vComplex::Nsimd(),{ + + // Four forward links from this site + auto Lx = Umu_v(ss)(IcosahedronPatchX); + auto Ly = Umu_v(ss)(IcosahedronPatchY); + auto Ld = Umu_v(ss)(IcosahedronPatchDiagonal); + auto Lt = Ut_v(ss)(); + IcosahedralStencilEntry *SE; + int64_t s; + int p; + /////////////////////////////////////////////////////////////////// + // Load all links required and multiple round loops + /////////////////////////////////////////////////////////////////// + SE = stencil_v.GetEntry(ent_Tp,ss); + s = SE->_offset; + p = SE->_permute; + auto Lx_at_tp = Umu_v(s)(IcosahedronPatchX); + auto Ly_at_tp = Umu_v(s)(IcosahedronPatchY); + auto Ld_at_tp = Umu_v(s)(IcosahedronPatchDiagonal); + if ( p ) { + rotate(Lx_at_tp,Lx_at_tp,1); + rotate(Ly_at_tp,Ly_at_tp,1); + rotate(Ld_at_tp,Ld_at_tp,1); + } + SE = stencil_v.GetEntry(ent_Xp,ss); + s = SE->_offset; + auto Lt_at_xp = Ut_v(s)(); + + SE = stencil_v.GetEntry(ent_Yp,ss); + s = SE->_offset; + auto Lt_at_yp = Ut_v(s)(); + + SE = stencil_v.GetEntry(ent_Dp,ss); + s = SE->_offset; + auto Lt_at_dp = Ut_v(s)(); + + Coordinate Coor; + EdgeGrid->oCoorFromOindex(Coor,ss); + + + SE = stencil_v.GetEntry(ent_Tm,ss); + s = SE->_offset; + p = SE->_permute; + auto Lt_at_tm = Ut_v(s)(); + auto Lx_at_tm = Umu_v(s)(IcosahedronPatchX); + auto Ly_at_tm = Umu_v(s)(IcosahedronPatchY); + auto Ld_at_tm = Umu_v(s)(IcosahedronPatchDiagonal); + const int Nsimdm1 = vComplex::Nsimd()-1; + + if ( p ) { + rotate(Lt_at_tm,Lt_at_tm,Nsimdm1); + rotate(Lx_at_tm,Lx_at_tm,Nsimdm1); + rotate(Ly_at_tm,Ly_at_tm,Nsimdm1); + rotate(Ld_at_tm,Ld_at_tm,Nsimdm1); + } + + SE = stencil_v.GetEntry(ent_TmXp,ss); + s = SE->_offset; + p = SE->_permute; + auto Lt_at_tmxp = Ut_v(s)(); + + SE = stencil_v.GetEntry(ent_TmYp,ss); + s = SE->_offset; + auto Lt_at_tmyp = Ut_v(s)(); + + SE = stencil_v.GetEntry(ent_TmDp,ss); + s = SE->_offset; + auto Lt_at_tmdp = Ut_v(s)(); + + if ( p ) { + rotate(Lt_at_tmxp,Lt_at_tmxp,Nsimdm1); + rotate(Lt_at_tmyp,Lt_at_tmyp,Nsimdm1); + rotate(Lt_at_tmdp,Lt_at_tmdp,Nsimdm1); + } + + // Lines start at origin and progress to top of link + // close as gauge invariant if multiplied from left by link + + coalescedWrite(stapleX_v[ss](),Lt_at_xp*adj(Lx_at_tp)*adj(Lt) + adj(Lt_at_tmxp)*adj(Lx_at_tm)*Lt_at_tm); + coalescedWrite(stapleY_v[ss](),Lt_at_yp*adj(Ly_at_tp)*adj(Lt) + adj(Lt_at_tmyp)*adj(Ly_at_tm)*Lt_at_tm); + coalescedWrite(stapleD_v[ss](),Lt_at_dp*adj(Ld_at_tp)*adj(Lt) + adj(Lt_at_tmdp)*adj(Ld_at_tm)*Lt_at_tm); + + int pol; + SE = stencil_v.GetEntry(ent_Xm,ss); + s = SE->_offset; + pol= SE->_polarisation; + auto Lt_at_xm = Ut_v(s)(); + auto Lx_at_xm = Umu_v(s)(pol); + + SE = stencil_v.GetEntry(ent_XmTp,ss); + s = SE->_offset; + p = SE->_permute; + pol= SE->_polarisation; + auto Lx_at_xmtp = Umu_v(s)(pol); + + SE = stencil_v.GetEntry(ent_Ym,ss); + s = SE->_offset; + pol= SE->_polarisation; + auto Lt_at_ym = Ut_v(s)(); + auto Ly_at_ym = Umu_v(s)(pol); + + SE = stencil_v.GetEntry(ent_YmTp,ss); + s = SE->_offset; + p = SE->_permute; + pol= SE->_polarisation; + auto Ly_at_ymtp = Umu_v(s)(pol); + + if ( p ) { + rotate(Ly_at_ymtp,Ly_at_ymtp,1); + rotate(Lx_at_xmtp,Lx_at_xmtp,1); + } + + auto stapT = Lx; + stapT= Lx_at_tp*adj(Lt_at_xp)*adj(Lx); + stapT= stapT + Ly_at_tp*adj(Lt_at_yp)*adj(Ly); + stapT= stapT + Ld_at_tp*adj(Lt_at_dp)*adj(Ld); + stapT= stapT + adj(Lx_at_xmtp)*adj(Lt_at_xm)*Lx_at_xm; + stapT= stapT + adj(Ly_at_ymtp)*adj(Lt_at_ym)*Ly_at_ym; + + SE = stencil_v.GetEntry(ent_Dm,ss); + int missingLink = SE->_missing_link; + s = SE->_offset; + pol= SE->_polarisation; + if ( !missingLink ) { + + auto Lt_at_dm = Ut_v(s)(); + auto Ld_at_dm = Umu_v(s)(pol); + SE = stencil_v.GetEntry(ent_DmTp,ss); + s = SE->_offset; + auto Ld_at_dmtp = Umu_v(s)(pol); + if ( p ) { + rotate(Ld_at_dmtp,Ld_at_dmtp,1); + } + stapT = stapT + adj(Ld_at_dmtp)*adj(Lt_at_dm)*Ld_at_dm; + } + + coalescedWrite(stapleT_v[ss](),stapT); + + }); + + const int ent_h0 = 0; + const int ent_h0tp = 1; + const int ent_h1 = 2; + const int ent_h1tp = 3; + const int ent_h2 = 4; + const int ent_h2tp = 5; + const int ent_h3 = 6; + const int ent_h3tp = 7; + const int ent_h4 = 8; + const int ent_h4tp = 9; + + uint64_t cart_sites = EdgeGrid->oSites(); + uint64_t pole_sites = VertexGrid->NorthPoleOsites() + + VertexGrid->SouthPoleOsites(); + + +#define LINK(ENT_0,ENT_1) \ + SE = stencil_v.GetEntry(ENT_0,ss); \ + s = SE->_offset; \ + pol = SE->_polarisation; \ + \ + Li = Umu_v(s)(pol); \ + Lth = Ut_v(s)(); \ + \ + SE=stencil_v.GetEntry(ENT_1,ss); \ + s = SE->_offset; \ + p = SE->_permute; \ + Litp = Umu_v(s)(pol); \ + if ( p ) { \ + rotate(Litp,Litp,1); \ + } + +#define LINK0(ENT_0,ENT_1) LINK(ENT_0,ENT_1); auto stapT = adj(Litp)*adj(Lth)*Li; +#define LINKn(ENT_0,ENT_1) LINK(ENT_0,ENT_1); stapT = stapT + adj(Litp)*adj(Lth)*Li; + + accelerator_for(sss,pole_sites,vComplex::Nsimd(),{ + + IcosahedralStencilEntry *SE; + int64_t s; + int p; + int pol; + + uint64_t ss = sss+cart_sites; + + auto Lt = Ut_v(ss)(); + auto Lth = Lt; + auto Li = Lt; + auto Litp= Lt; + + LINK0(ent_h0,ent_h0tp); + LINKn(ent_h1,ent_h1tp); + LINKn(ent_h2,ent_h2tp); + LINKn(ent_h3,ent_h3tp); + LINKn(ent_h4,ent_h4tp); + + coalescedWrite(stapleT_v[ss](),stapT); + + }); +#undef LINK +#undef LINK0 +#undef LINKn + } + + + template + void CopyNonPoles(Field &in,Field &out) + { + autoView(out_v,out,AcceleratorWrite); + autoView(in_v,in,AcceleratorRead); + + accelerator_for(ss,EdgeGrid->oSites(),vComplex::Nsimd(),{ + + auto i = in_v(ss)(); + coalescedWrite(out_v[ss](),i); + }); + } template void CovariantLaplacian(MatterField &in,MatterField &out,DoubledGaugeField &Uds) @@ -319,6 +849,9 @@ 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(),{ @@ -385,6 +918,9 @@ public: const int ent_Ym = 4; const int ent_Dm = 5; + const int ent_Tp = 6; + const int ent_Tm = 7; + accelerator_for(ss,VertexGrid->CartesianOsites(),vComplex::Nsimd(),{ // Three local links @@ -432,6 +968,8 @@ public: coalescedWrite(Uds_v[ss](4),Ly_at_ym ); } int missingLink; + auto Link = adj(U_v(ss)(0)); + Link=Zero(); { auto SE = stencil_v.GetEntry(ent_Dm,ss); auto s = SE->_offset; @@ -442,8 +980,12 @@ public: auto Ld_at_dm = U_v(s)(pol); Ld_at_dm = adj(Ld_at_dm); coalescedWrite(Uds_v[ss](5),Ld_at_dm ); + } else { + coalescedWrite(Uds_v[ss](5),Link); } } + coalescedWrite(Uds_v[ss](6),Link); + coalescedWrite(Uds_v[ss](7),Link); }); auto pole_sites = VertexGrid->oSites() - VertexGrid->CartesianOsites(); auto pole_offset= VertexGrid->CartesianOsites(); @@ -455,16 +997,21 @@ public: auto pol= SE->_polarisation; auto Link = adj(U_v(s)(pol)); coalescedWrite(Uds_v[pole_offset+ss](p),Link); - } + } + 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); }); } - void GaugeTransform(GaugeLinkField >, GaugeField &Umu) + void GaugeTransform(GaugeLinkField >, GaugeField &Umu, GaugeLinkField &Ut) { autoView(Umu_v,Umu,AcceleratorWrite); + autoView(Ut_v,Ut,AcceleratorWrite); autoView(g_v,gt,AcceleratorRead); - autoView(stencil_v,NNev,AcceleratorRead); - - const int np = NNev._npoints; + autoView(stencil_v ,NNev,AcceleratorRead); + autoView(stencil_vv,NNvv,AcceleratorRead); const int ent_Xp = 0; const int ent_Yp = 1; @@ -503,15 +1050,106 @@ public: coalescedWrite(Umu_v[ss](IcosahedronPatchY),ly); coalescedWrite(Umu_v[ss](IcosahedronPatchDiagonal),ld); }); - } + const int ent_Tp = 6; + accelerator_for(ss,VertexGrid->oSites(),vComplex::Nsimd(),{ + auto lt = Ut_v(ss)(); + auto g = g_v(ss)(); + auto SE = stencil_vv.GetEntry(ent_Tp,ss); + uint64_t tp_idx = SE->_offset; + auto permute = SE->_permute; + auto gtp = g_v(tp_idx)(); + // std::cout << ss<<" tp "<_permute<<" : g"<oSites(),vComplex::Nsimd(),{ + + // Three forward links from this site + auto Lx = Uj_v(ss)(IcosahedronPatchX); + auto Ly = Uj_v(ss)(IcosahedronPatchY); + auto Ld = Uj_v(ss)(IcosahedronPatchDiagonal); + auto Lt = Ut_v(ss)(); + + // Three forward links from this site at Tp + auto SE = stencil_v.GetEntry(ent_Tp,ss); + uint64_t tp_idx = SE->_offset; + uint64_t permute= SE->_permute; + + auto Ltp_x = Uj_v(tp_idx)(IcosahedronPatchX); + auto Ltp_y = Uj_v(tp_idx)(IcosahedronPatchY); + auto Ltp_d = Uj_v(tp_idx)(IcosahedronPatchDiagonal); + + if ( permute ) { + rotate(Ltp_x,Ltp_x,1); // does in=out work? + rotate(Ltp_y,Ltp_y,1); + rotate(Ltp_d,Ltp_d,1); + } + + SE = stencil_v.GetEntry(ent_Xp,ss); + uint64_t xp_idx = SE->_offset; + + SE = stencil_v.GetEntry(ent_Yp,ss); + uint64_t yp_idx = SE->_offset; + + SE = stencil_v.GetEntry(ent_Dp,ss); + uint64_t dp_idx = SE->_offset; + + auto Lxp_t = Ut_v(xp_idx)(); + auto Lyp_t = Ut_v(yp_idx)(); + auto Ldp_t = Ut_v(dp_idx)(); + + coalescedWrite(plaq1_v[ss](),trace(Lx*Lxp_t*adj(Ltp_x)*adj(Lt))); + coalescedWrite(plaq2_v[ss](),trace(Ly*Lyp_t*adj(Ltp_y)*adj(Lt))); + coalescedWrite(plaq3_v[ss](),trace(Ld*Ldp_t*adj(Ltp_d)*adj(Lt))); + + /* + std::cout << ss <<"----------------"< seeds({1,2,3,4}); GridParallelRNG vRNG(&EdgeGrid); vRNG.SeedFixedIntegers(seeds); - LatticeColourMatrix g(&VertexGrid); + GaugeLinkField g(&VertexGrid); LatticeReal gr(&VertexGrid); - LatticeComplex gc(&VertexGrid); + ComplexField gc(&VertexGrid); // gr = Zero(); gaussian(vRNG,gr); Complex ci(0.0,1.0); @@ -599,19 +1245,16 @@ int main (int argc, char ** argv) std::cout << GridLogMessage << " Check plaquette is gauge invariant "<