portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2026-02-11 01:13:27 +00:00
Code Releases Activity

Gauge staples for temporal direction added (ico-T staples and T-ico staples).

Browse Source 
Passes gauge covariance test, requiring the link x its staples = 1 on a random gauge transform.

Still to do:

temporal link double store
temporal laplacian + covariant laplacian terms

Make the "rotate" -> coalescedReadRotate as presently it is hardwired CPU domain in a couple of places in Test_icosahedron
This commit is contained in:
Peter Boyle
2025-11-06 13:54:39 -05:00
parent d3ca16c76d
commit d5c0d54f89
4 changed files with 1358 additions and 72 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
Grid/cartesian/CartesianCrossIcosahedron.h
View File

@@ -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;

568
Grid/stencil/IcosahedralStencil.h
View File

@@ -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:
@@ -87,6 +88,7 @@ public:
// to a corresponding "backward link" on the pole
deviceVector<IcosahedralStencilEntry> _entries;
void GetNbrForPlusX(GridBase *grid,Coordinate &Coor,Coordinate &NbrCoor, int &isPole)
{
int nd = grid->Nd();
@@ -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<< "------------------------------------"<<std::endl;
std::cout << GridLogMessage<< "NOT testing diag -x-y = identity "<<std::endl;
std::cout << GridLogMessage<< "NOT testing diag -y-x = identity"<<std::endl;
std::cout << GridLogMessage<< "------------------------------------"<<std::endl;
std::cout << GridLogMessage<< "------------------------------------"<<std::endl;
std::cout << GridLogMessage<< "NOT testing -diag = -x-y "<<std::endl;
std::cout << GridLogMessage<< "NOT testing -diag = -y-x "<<std::endl;
std::cout << GridLogMessage<< "------------------------------------"<<std::endl;
std::cout << GridLogMessage<< "*************************************"<<std::endl;
std::cout << GridLogMessage<< " Icosahedral Stencil Geometry Test Complete"<<std::endl;
std::cout << GridLogMessage<< "*************************************"<<std::endl;
@@ -577,6 +569,7 @@ public:
{
this->_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?
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,6 +669,12 @@ public:
GetNbrForMinusX(grid,Coor,XmCoor);
GetNbrForMinusY(grid,Coor,YmCoor);
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;
@@ -686,6 +693,7 @@ public:
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;site<cart_sites; site ++) { // loops over volume
Coordinate Coor;
Coordinate tCoor;
vertexgrid->oCoorFromOindex(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 "<<Coor<<" connects to south pole_site "<<pole_site<<std::endl;
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);
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;site<cart_sites; site ++) {
Coordinate Coor;
Coordinate tCoor;
vertexgrid->oCoorFromOindex(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 "<<SE._offset<<" for north site "<<lex<<std::endl;
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);
// std::cout << " Put nbr "<<SE._offset<<" for north site "<<lex<<std::endl;
}
}
}
@@ -873,7 +942,6 @@ public:
for(uint64_t site=0;site<cart_sites; site ++) {
Coordinate Coor;
Coordinate NbrCoor;
int nd = grid->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;site<cart_sites; site ++) {
int nd = grid->Nd();
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(Patch<IcosahedralPatches);
assert((north==1)||(south==1));
Coordinate XpCoor;
Coordinate YpCoor;
Coordinate DpCoor;
Coordinate TpCoor;
GetNbrForPlusDiagonal(grid,Coor,DpCoor);
GetNbrForPlusX(grid,Coor,XpCoor,isPoleX);
GetNbrForPlusY(grid,Coor,YpCoor,isPoleY);
int DpPatch = DpCoor[nd-1];
int DpHemiPatch = DpCoor[nd-1]%HemiPatches;
int DpHemisphere = DpCoor[nd-1]/HemiPatches;
SE._is_local = true;
SE._polarisation = 0;
SE._adjoint = false;
SE._missing_link = false;
SE._permute=0;
//////////////////////////////////////////////////
// Forward one site in time direction
//////////////////////////////////////////////////
int tdim = 2;
int delta = 1;
int permute;
GetOrthogNbr(grid,Coor,TpCoor,tdim,delta,permute);
SE._offset = vertexgrid->oIndex(TpCoor);
SE._permute=permute;
// std::cout << " Plaq stencil "<<Coor<<" Tp "<<TpCoor<<" perm "<<permute<<std::endl;
acceleratorPut(this->_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 "<<Coor<<" DpCoor "<<DpCoor<<" site "<<SE._offset<<std::endl;
acceleratorPut(this->_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;site<cart_sites; site ++) {
int nd = grid->Nd();
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(Patch<IcosahedralPatches);
assert((north==1)||(south==1));
Coordinate XpCoor;
Coordinate YpCoor;
Coordinate DpCoor;
Coordinate TpCoor;
Coordinate TmCoor;
Coordinate TmXpCoor;
Coordinate TmYpCoor;
Coordinate TmDpCoor;
int missingLink;
Coordinate XmCoor;
Coordinate YmCoor;
Coordinate DmCoor;
Coordinate XmTpCoor;
Coordinate YmTpCoor;
Coordinate DmTpCoor;
GetNbrForPlusDiagonal(grid,Coor,DpCoor);
GetNbrForPlusX(grid,Coor,XpCoor,isPoleX);
GetNbrForPlusY(grid,Coor,YpCoor,isPoleY);
GetNbrForMinusDiagonal(grid,Coor,DmCoor,missingLink);
GetNbrForMinusX(grid,Coor,XmCoor);
GetNbrForMinusY(grid,Coor,YmCoor);
int DpPatch = DpCoor[nd-1];
int DpHemiPatch = DpCoor[nd-1]%HemiPatches;
int DpHemisphere = DpCoor[nd-1]/HemiPatches;
int XmHemiPatch = XmCoor[nd-1]%HemiPatches;
int XmHemisphere = XmCoor[nd-1]/HemiPatches;
int YmHemiPatch = YmCoor[nd-1]%HemiPatches;
int YmHemisphere = YmCoor[nd-1]/HemiPatches;
int DmPatch = DmCoor[nd-1];
int DmHemiPatch = DmCoor[nd-1]%HemiPatches;
int DmHemisphere = DmCoor[nd-1]/HemiPatches;
SE._is_local = true;
SE._polarisation = 0;
SE._adjoint = false;
SE._missing_link = false;
SE._permute=0;
//////////////////////////////////////////////////
// Forward one site in time direction
//////////////////////////////////////////////////
int tdim = 2;
int delta = 1;
int permute;
GetOrthogNbr(grid,Coor,TpCoor,tdim,delta,permute);
SE._offset = vertexgrid->oIndex(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;site<cart_sites; site ++) { // loops over volume
Coordinate Coor;
Coordinate tCoor;
vertexgrid->oCoorFromOindex(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;site<cart_sites; site ++) {
Coordinate Coor;
Coordinate tCoor;
vertexgrid->oCoorFromOindex(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);

2
systems/mac-arm/config-command
View File

@@ -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/ \

838
tests/debug/Test_icosahedron.cc
View File

@@ -31,14 +31,29 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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<typename vtype> using iIcosahedralLorentzColourMatrix = iVector<iScalar<iMatrix<vtype, Nc> >, NumIcosahedralPolarizations> ;
typedef iIcosahedralLorentzColourMatrix<Complex> IcosahedralLorentzColourMatrix;
typedef iIcosahedralLorentzColourMatrix<vComplex > vIcosahedralLorentzColourMatrix;
typedef Lattice<vIcosahedralLorentzColourMatrix> LatticeIcosahedralLorentzColourMatrix;
class IcosahedralGimpl
{
public:
typedef LatticeLorentzColourMatrix GaugeField;
typedef LatticeColourMatrix GaugeLinkField;
typedef LatticeIcosahedralLorentzColourMatrix GaugeField;
typedef LatticeDoubledGaugeField DoubledGaugeField;
typedef LatticeColourMatrix GaugeLinkField;
typedef LatticeComplex ComplexField;
};
@@ -60,27 +75,34 @@ public:
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),
NNee(_EdgeGrid,_VertexGrid),
NNev(_EdgeGrid,_VertexGrid),
NNve(_EdgeGrid,_VertexGrid),
NNvv(_EdgeGrid,_VertexGrid),
TimePlaqStencil(_EdgeGrid,_VertexGrid),
TimeStapleStencil(_EdgeGrid,_VertexGrid),
VertexGrid(_VertexGrid), EdgeGrid(_EdgeGrid)
{
FaceStencil.FaceStencil();
std::cout << "NNee"<<std::endl;
// true/false is "vertexInput, VertexOutput"
NNee.NearestNeighbourStencil(false,false);// edge input + output ; used by face stencil
NNev.NearestNeighbourStencil(true,false); // vertex input, edge ouput ; used by gauge transform
NNvv.NearestNeighbourStencil(true,true); // vertex input + output ; used by Laplacian
NNve.NearestNeighbourStencil(false,true); // edge input, vertex output; used by double store
TimePlaqStencil.TemporalPlaquetteStencil();
TimeStapleStencil.TemporalStapleStencil();
}
////////////////////////////////////////////////////////////////////////////////////
// Gauge Link field GT is the gauge transform and lives on the VERTEX field
////////////////////////////////////////////////////////////////////////////////////
void ForwardTriangles(GaugeField &Umu,LatticeComplex &plaq1,LatticeComplex &plaq2)
void ForwardTriangles(GaugeField &Umu,ComplexField &plaq1,ComplexField &plaq2)
{
{
autoView(Umu_v,Umu,AcceleratorRead);
@@ -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<<" "<<ENT_0<<" plaq "<<Lt*stap<<std::endl;}
accelerator_for(sss,pole_sites,vComplex::Nsimd(),{
IcosahedralStencilEntry *SE;
int64_t s;
int p;
int pol;
uint64_t ss = sss+cart_sites;
LINK(ent_h0,ent_h0tp);
LINK(ent_h1,ent_h1tp);
LINK(ent_h2,ent_h2tp);
LINK(ent_h3,ent_h3tp);
LINK(ent_h4,ent_h4tp);
});
#undef LINK
}
*/
void TemporalStaples(GaugeLinkField &Ut,
GaugeField &Umu,
GaugeLinkField &stapleX,
GaugeLinkField &stapleY,
GaugeLinkField &stapleD,
GaugeLinkField &stapleT)
{
autoView(Umu_v,Umu,AcceleratorRead);
autoView(Ut_v,Ut,AcceleratorRead);
autoView(stapleX_v,stapleX,AcceleratorWrite);
autoView(stapleY_v,stapleY,AcceleratorWrite);
autoView(stapleD_v,stapleD,AcceleratorWrite);
autoView(stapleT_v,stapleT,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 ) {
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<class Field>
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<class MatterField>
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();
@@ -456,15 +998,20 @@ public:
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 &gt, GaugeField &Umu)
void GaugeTransform(GaugeLinkField &gt, 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 "<<tp_idx<<" permute "<<(int)SE->_permute<<" : g"<<g<<" l "<<lt<<" gtp "<<gtp<<std::endl;
if ( permute ) rotate(gtp,gtp,1);
lt = g*lt*adj(gtp);
coalescedWrite(Ut_v[ss](),lt);
});
}
void TemporalPlaquette(GaugeLinkField &Ut, GaugeField &Uj,ComplexField &plaq1,ComplexField &plaq2,ComplexField &plaq3)
{
autoView(Uj_v,Uj,AcceleratorRead);
autoView(Ut_v,Ut,AcceleratorRead);
autoView(plaq1_v,plaq1,AcceleratorWrite);
autoView(plaq2_v,plaq2,AcceleratorWrite);
autoView(plaq3_v,plaq3,AcceleratorWrite);
autoView(stencil_v,TimePlaqStencil,AcceleratorRead);
const int np = TimePlaqStencil._npoints;
const int ent_Tp = 0;
const int ent_Xp = 1;
const int ent_Yp = 2;
const int ent_Dp = 3;
accelerator_for(ss,EdgeGrid->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 <<"----------------"<<std::endl;
std::cout << ss <<" 1 "<<plaq1_v[ss]<<std::endl;
std::cout << ss <<" 2 "<<plaq2_v[ss]<<std::endl;
std::cout << ss <<" 3 "<<plaq3_v[ss]<<std::endl;
std::cout << ss <<" xp_idx "<<xp_idx<<std::endl;
std::cout << ss <<" yp_idx "<<yp_idx<<std::endl;
std::cout << ss <<" dp_idx "<<dp_idx<<std::endl;
std::cout << ss <<" tp_idx "<<dp_idx<<std::endl;
std::cout << ss <<" Lx "<<Lx<<std::endl;
std::cout << ss <<" Ly "<<Ly<<std::endl;
std::cout << ss <<" Ld "<<Ld<<std::endl;
std::cout << ss <<" Lt "<<Lt<<std::endl;
std::cout << ss <<" Ldp_t "<<Ldp_t<<std::endl;
std::cout << ss <<" Lxp_t "<<Lxp_t<<std::endl;
std::cout << ss <<" Lyp_t "<<Lyp_t<<std::endl;
std::cout << ss <<" Ltp_x "<<Ltp_x<<std::endl;
std::cout << ss <<" Ltp_y "<<Ltp_y<<std::endl;
std::cout << ss <<" Ltp_d "<<Ltp_d<<std::endl;
*/
});
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int L=8;
// const int L=8;
const int Npatch = IcosahedralPatches;
// Put SIMD all in time direction
@@ -527,9 +1165,15 @@ int main (int argc, char ** argv)
GridCartesianCrossIcosahedron VertexGrid(latt_size,simd_layout,mpi_layout,IcosahedralVertices);
std::cout << GridLogMessage << " Created vertex grid "<<std::endl;
LatticeLorentzColourMatrix Umu(&EdgeGrid);
LatticeLorentzColourMatrix Umuck(&EdgeGrid);
LatticeComplex Phi(&VertexGrid);
typedef IcosahedralGimpl::ComplexField ComplexField;
typedef IcosahedralGimpl::GaugeLinkField GaugeLinkField;
typedef IcosahedralGimpl::GaugeField GaugeField;
typedef IcosahedralGimpl::DoubledGaugeField DoubledGaugeField;
GaugeLinkField Ut(&VertexGrid);
GaugeLinkField Utnp(&EdgeGrid);
GaugeField Umu(&EdgeGrid);
GaugeField Umuck(&EdgeGrid);
ComplexField Phi(&VertexGrid);
std::cout << GridLogMessage << " Created two fields "<<std::endl;
Phi = Zero();
@@ -539,6 +1183,7 @@ int main (int argc, char ** argv)
Complex one (1.0);
Phi = one;
Umu = one;
Ut = one;
std::cout << GridLogMessage << " V = "<<norm2(Phi)<<std::endl;
std::cout << GridLogMessage << " Expect "<<latt_size[0]*latt_size[1]*latt_size[2]*10+2*latt_size[2]<<std::endl;
@@ -568,10 +1213,11 @@ int main (int argc, char ** argv)
Umu=one;
LatticeComplex plaq1(&EdgeGrid);
LatticeComplex plaq2(&EdgeGrid);
ComplexField plaq1(&EdgeGrid);
ComplexField plaq2(&EdgeGrid);
ComplexField plaq3(&EdgeGrid);
LatticeComplex plaq_ref(&EdgeGrid);
ComplexField plaq_ref(&EdgeGrid);
plaq_ref=1.0;
Support.ForwardTriangles(Umu,plaq1,plaq2);
@@ -585,9 +1231,9 @@ int main (int argc, char ** argv)
std::vector<int> 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 "<<std::endl;
std::cout << GridLogMessage << "****************************************"<<std::endl;
std::cout << GridLogMessage << " applying gauge transform"<<std::endl;
Support.GaugeTransform (g,Umu);
Support.GaugeTransform (g,Umu,Ut);
std::cout << GridLogMessage << " applied gauge transform "<<std::endl;
std::cout << GridLogMessage << " recalculating plaquette "<<std::endl;
Support.ForwardTriangles(Umu,plaq1,plaq2);
std::cout << GridLogMessage << " plaq1 "<< norm2(plaq1)<<std::endl;
std::cout << GridLogMessage << " plaq2 "<< norm2(plaq2)<<std::endl;
std::cout << GridLogMessage << " lower plaq1 "<< norm2(plaq1)<<std::endl;
std::cout << GridLogMessage << " upper plaq2 "<< norm2(plaq2)<<std::endl;
std::cout << GridLogMessage << " plaq1 err "<< norm2(plaq1-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " plaq2 err "<< norm2(plaq2-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " lower plaq1 err "<< norm2(plaq1-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " upper plaq2 err "<< norm2(plaq2-plaq_ref)<<std::endl;
typedef IcosahedralGimpl::GaugeLinkField GaugeLinkField;
typedef IcosahedralGimpl::GaugeField GaugeField;
GaugeLinkField stapleXY(&EdgeGrid);
GaugeLinkField stapleYX(&EdgeGrid);
@@ -635,11 +1278,9 @@ int main (int argc, char ** argv)
linkY = peekLorentz(Umu,IcosahedronPatchY);
linkD = peekLorentz(Umu,IcosahedronPatchDiagonal);
// OK
std::cout << GridLogMessage << " trace D*StapleXY "<<norm2(trace(linkD * stapleXY))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(linkD * stapleXY)-plaq_ref)<<std::endl;
// BAD
std::cout << GridLogMessage << " trace D*StapleYX "<<norm2(trace(linkD * stapleYX))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(linkD * stapleYX)-plaq_ref)<<std::endl;
@@ -662,7 +1303,7 @@ int main (int argc, char ** argv)
Support.Laplacian(in,out);
std::cout << GridLogMessage<< "Calling double storing gauge field" <<std::endl;
LatticeDoubledGaugeField Uds(&VertexGrid);
DoubledGaugeField Uds(&VertexGrid);
Support.DoubleStore(Umu,Uds);
Support.CovariantLaplacian(in,out,Uds);
@@ -677,16 +1318,149 @@ int main (int argc, char ** argv)
gout = g*out;
gin = g*in;
Support.GaugeTransform(g,Umu);
Support.GaugeTransform(g,Umu,Ut);
Support.DoubleStore(Umu,Uds);
Support.CovariantLaplacian(gin,out,Uds);
std::cout << GridLogMessage<< "Applied gauge transformed covariant laplacian to transformed vector !" <<std::endl;
auto ipgt = innerProduct(out, gin);
std::cout << "Testing D[U_gt](gF) = g D[U] F : defect is "<<norm2(out-gout)<<std::endl;
std::cout <<GridLogMessage<< "Testing D[U_gt](gF) = g D[U] F : defect is "<<norm2(out-gout)<<std::endl;
ip = ip - ipgt;
std::cout << "Testing F D[U](F) = (gF) D[U_gt] gF : defect is "<<ip<<std::endl;
std::cout <<GridLogMessage<< "Testing F D[U](F) = (gF) D[U_gt] gF : defect is "<<ip<<std::endl;
/*
* Temporal plaquettes
*/
std::cout << GridLogMessage<< " Computing temporal plaquettes"<<std::endl;
Support.TemporalPlaquette(Ut, Umu,plaq1,plaq2,plaq3);
std::cout << GridLogMessage<<" Computed temporal plaquettes"<<std::endl;
std::cout << GridLogMessage<<" XT plaq1 "<<norm2(plaq1)<<std::endl;
std::cout << GridLogMessage<<" YT plaq2 "<<norm2(plaq2)<<std::endl;
std::cout << GridLogMessage<<" DT plaq3 "<<norm2(plaq3)<<std::endl;
std::cout << GridLogMessage << " XT plaq1 err "<< norm2(plaq1-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " YT plaq2 err "<< norm2(plaq2-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " DT plaq3 err "<< norm2(plaq3-plaq_ref)<<std::endl;
GaugeLinkField stapleX(&EdgeGrid);
GaugeLinkField stapleY(&EdgeGrid);
GaugeLinkField stapleD(&EdgeGrid);
GaugeLinkField stapleT(&VertexGrid); // Note live on different Grids
std::cout << GridLogMessage<<" Computing temporal staples"<<std::endl;
Support.TemporalStaples(Ut,Umu,stapleX,stapleY,stapleD,stapleT);
std::cout << GridLogMessage<<" Computed temporal staples"<<std::endl;
linkX = peekLorentz(Umu,IcosahedronPatchX);
linkY = peekLorentz(Umu,IcosahedronPatchY);
linkD = peekLorentz(Umu,IcosahedronPatchDiagonal);
std::cout << GridLogMessage <<"------------------------------------------------"<<std::endl;
std::cout << GridLogMessage << " trace X*StapleX "<<norm2(trace(linkX * stapleX))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(linkX * stapleX)-2.0*plaq_ref)<<std::endl;
std::cout << GridLogMessage << " trace Y*StapleY "<<norm2(trace(linkY * stapleY))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(linkY * stapleY)-2.0*plaq_ref)<<std::endl;
std::cout << GridLogMessage << " trace D*StapleD "<<norm2(trace(linkD * stapleD))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(linkD * stapleD)-2.0*plaq_ref)<<std::endl;
std::cout << GridLogMessage <<"------------------------------------------------"<<std::endl;
int L = latt_size[0];
int T = latt_size[2];
double T_staple_norm = 10*(L*L-1)*T*36.0 + 12*T*25.0;
double tr = norm2(trace(Ut * stapleT));
std::cout << GridLogMessage << " trace T*StapleT "<<tr<<std::endl;
std::cout << GridLogMessage << " got " << tr<<" expect "<<T_staple_norm<<" diff "<<tr-T_staple_norm <<std::endl;
/*
GaugeLinkField stapleXTp(&EdgeGrid);
GaugeLinkField stapleYTp(&EdgeGrid);
GaugeLinkField stapleDTp(&EdgeGrid);
GaugeLinkField stapleXTm(&EdgeGrid);
GaugeLinkField stapleYTm(&EdgeGrid);
GaugeLinkField stapleDTm(&EdgeGrid);
GaugeLinkField stapleTXp(&EdgeGrid);
GaugeLinkField stapleTYp(&EdgeGrid);
GaugeLinkField stapleTDp(&EdgeGrid);
GaugeLinkField stapleTXm(&EdgeGrid);
GaugeLinkField stapleTYm(&EdgeGrid);
GaugeLinkField stapleTDm(&EdgeGrid);
std::cout << GridLogMessage<<" Computing temporal staples"<<std::endl;
Support.TemporalStaplesOld(Ut,Umu,
stapleXTp,stapleYTp,stapleDTp,
stapleXTm,stapleYTm,stapleDTm,
stapleTXp,stapleTYp,stapleTDp,
stapleTXm,stapleTYm,stapleTDm);
std::cout << GridLogMessage<<" Computed temporal staples"<<std::endl;
linkX = peekLorentz(Umu,IcosahedronPatchX);
linkY = peekLorentz(Umu,IcosahedronPatchY);
linkD = peekLorentz(Umu,IcosahedronPatchDiagonal);
Support.CopyNonPoles(Ut,Utnp);
std::cout << GridLogMessage <<"------------------------------------------------"<<std::endl;
std::cout << GridLogMessage << " trace X*StapleXTp "<<norm2(trace(linkX * stapleXTp))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(linkX * stapleXTp)-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " trace Y*StapleYTp "<<norm2(trace(linkY * stapleYTp))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(linkY * stapleYTp)-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " trace D*StapleDTp "<<norm2(trace(linkD * stapleDTp))<<std::endl;
std::cout << GridLogMessage << " norm StapleDTp "<<norm2(stapleDTp)<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(linkD * stapleDTp)-plaq_ref)<<std::endl;
std::cout << GridLogMessage <<"------------------------------------------------"<<std::endl;
std::cout << GridLogMessage << " trace T*StapleTXp "<<norm2(trace(Utnp * stapleTXp))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(Utnp * stapleTXp)-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " trace T*StapleTYp "<<norm2(trace(Utnp * stapleTYp))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(Utnp * stapleTYp)-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " trace T*StapleTDp "<<norm2(trace(Utnp * stapleTDp))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(Utnp * stapleTDp)-plaq_ref)<<std::endl;
std::cout << GridLogMessage <<"------------------------------------------------"<<std::endl;
std::cout << GridLogMessage << " trace X*StapleXTm "<<norm2(trace(linkX * stapleXTm))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(linkX * stapleXTm)-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " trace Y*StapleYTm "<<norm2(trace(linkY * stapleYTm))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(linkY * stapleYTm)-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " trace D*StapleDTm "<<norm2(trace(linkD * stapleDTm))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(linkD * stapleDTm)-plaq_ref)<<std::endl;
std::cout << GridLogMessage <<"------------------------------------------------"<<std::endl;
std::cout << GridLogMessage << " trace T*StapleTXm "<<norm2(trace(Utnp * stapleTXm))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(Utnp * stapleTXm)-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " trace T*StapleTYm "<<norm2(trace(Utnp * stapleTYm))<<std::endl;
std::cout << GridLogMessage << " err " << norm2(trace(Utnp * stapleTYm)-plaq_ref)<<std::endl;
std::cout << GridLogMessage << " trace T*StapleTDm "<<norm2(trace(Utnp * stapleTDm))<<std::endl;
std::cout << GridLogMessage << " diff expected to be number of exceptional points " << norm2(trace(Utnp * stapleTDm)-plaq_ref)<<std::endl;
std::cout << GridLogMessage <<"------------------------------------------------"<<std::endl;
// std::cout << closure(Utnp*stapleTXm) << std::endl;
*/
Grid_finalize();
}