Nodes on padded cell

Grid/lattice/Lattice_transfer.h

@@ -707,9 +707,9 @@ void localCopyRegion(const Lattice<vobj> &From,Lattice<vobj> & To,Coordinate Fro
   Coordinate ist = Tg->_istride;
   Coordinate ost = Tg->_ostride;
 
-  autoView( t_v , To, AcceleratorWrite);
-  autoView( f_v , From, AcceleratorRead);
-  accelerator_for(idx,Fg->lSites(),1,{
+  autoView( t_v , To, CpuWrite);
+  autoView( f_v , From, CpuRead);
+  thread_for(idx,Fg->lSites(),{
     sobj s;
     Coordinate Fcoor(nd);
     Coordinate Tcoor(nd);
@@ -722,15 +722,20 @@ void localCopyRegion(const Lattice<vobj> &From,Lattice<vobj> & To,Coordinate Fro
       Tcoor[d] = ToLowerLeft[d]+ Fcoor[d]-FromLowerLeft[d];
     }
     if (in_region) {
-      Integer idx_f = 0; for(int d=0;d<nd;d++) idx_f+=isf[d]*(Fcoor[d]/rdf[d]);
-      Integer idx_t = 0; for(int d=0;d<nd;d++) idx_t+=ist[d]*(Tcoor[d]/rdt[d]);
-      Integer odx_f = 0; for(int d=0;d<nd;d++) odx_f+=osf[d]*(Fcoor[d]%rdf[d]);
-      Integer odx_t = 0; for(int d=0;d<nd;d++) odx_t+=ost[d]*(Tcoor[d]%rdt[d]);
-      vector_type * fp = (vector_type *)&f_v[odx_f];
-      vector_type * tp = (vector_type *)&t_v[odx_t];
+#if 0      
+      Integer idx_f = 0; for(int d=0;d<nd;d++) idx_f+=isf[d]*(Fcoor[d]/rdf[d]); // inner index from
+      Integer idx_t = 0; for(int d=0;d<nd;d++) idx_t+=ist[d]*(Tcoor[d]/rdt[d]); // inner index to
+      Integer odx_f = 0; for(int d=0;d<nd;d++) odx_f+=osf[d]*(Fcoor[d]%rdf[d]); // outer index from
+      Integer odx_t = 0; for(int d=0;d<nd;d++) odx_t+=ost[d]*(Tcoor[d]%rdt[d]); // outer index to
+      scalar_type * fp = (scalar_type *)&f_v[odx_f];
+      scalar_type * tp = (scalar_type *)&t_v[odx_t];
       for(int w=0;w<words;w++){
 	tp[w].putlane(fp[w].getlane(idx_f),idx_t);
       }
+#else
+    peekLocalSite(s,f_v,Fcoor);
+    pokeLocalSite(s,t_v,Tcoor);
+#endif
     }
   });
 }
@@ -841,9 +846,9 @@ void InsertSliceLocal(const Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int
 
   for(int d=0;d<nh;d++){
     if ( d!=orthog ) {
-    assert(lg->_processors[d]  == hg->_processors[d]);
-    assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
-  }
+      assert(lg->_processors[d]  == hg->_processors[d]);
+      assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
+    }
   }
 
   // the above should guarantee that the operations are local

Grid/lattice/PaddedCell.h

@@ -0,0 +1,136 @@
+/*************************************************************************************
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/lattice/PaddedCell.h
+
+    Copyright (C) 2019
+
+Author: Peter Boyle pboyle@bnl.gov
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+class PaddedCell {
+public:
+  GridCartesian * unpadded_grid;
+  int dims;
+  int depth;
+  std::vector<GridCartesian *> grids;
+  ~PaddedCell()
+  {
+    DeleteGrids();
+  }
+  PaddedCell(int _depth,GridCartesian *_grid)
+  {
+    unpadded_grid = _grid;
+    depth=_depth;
+    dims=_grid->Nd();
+    AllocateGrids();
+    Coordinate local     =unpadded_grid->LocalDimensions();
+    for(int d=0;d<dims;d++){
+      assert(local[d]>=depth);
+    }
+  }
+  void DeleteGrids(void)
+  {
+    for(int d=0;d<grids.size();d++){
+      delete grids[d];
+    }
+    grids.resize(0);
+  };
+  void AllocateGrids(void)
+  {
+    Coordinate local     =unpadded_grid->LocalDimensions();
+    Coordinate simd      =unpadded_grid->_simd_layout;
+    Coordinate processors=unpadded_grid->_processors;
+    Coordinate plocal    =unpadded_grid->LocalDimensions();
+    Coordinate global(dims);
+
+    // expand up one dim at a time
+    for(int d=0;d<dims;d++){
+
+      plocal[d] += 2*depth; 
+
+      for(int d=0;d<dims;d++){
+	global[d] = plocal[d]*processors[d];
+      }
+
+      grids.push_back(new GridCartesian(global,simd,processors));
+    }
+  };
+  template<class vobj>
+  inline Lattice<vobj> Extract(Lattice<vobj> &in)
+  {
+    Lattice<vobj> out(unpadded_grid);
+
+    Coordinate local     =unpadded_grid->LocalDimensions();
+    Coordinate fll(dims,depth); // depends on the MPI spread
+    Coordinate tll(dims,0); // depends on the MPI spread
+    localCopyRegion(in,out,fll,tll,local);
+    return out;
+  }
+  template<class vobj>
+  inline Lattice<vobj> Exchange(Lattice<vobj> &in)
+  {
+    GridBase *old_grid = in.Grid();
+    int dims = old_grid->Nd();
+    Lattice<vobj> tmp = in;
+    for(int d=0;d<dims;d++){
+      tmp = Expand(d,tmp); // rvalue && assignment
+    }
+    return tmp;
+  }
+  // expand up one dim at a time
+  template<class vobj>
+  inline Lattice<vobj> Expand(int dim,Lattice<vobj> &in)
+  {
+    GridBase *old_grid = in.Grid();
+    GridCartesian *new_grid = grids[dim];//These are new grids
+    Lattice<vobj>  padded(new_grid);
+    Lattice<vobj> shifted(old_grid);    
+    Coordinate local     =old_grid->LocalDimensions();
+    Coordinate plocal    =new_grid->LocalDimensions();
+    if(dim==0) conformable(old_grid,unpadded_grid);
+    else       conformable(old_grid,grids[dim-1]);
+
+    std::cout << " dim "<<dim<<" local "<<local << " padding to "<<plocal<<std::endl;
+    // Middle bit
+    for(int x=0;x<local[dim];x++){
+      InsertSliceLocal(in,padded,x,depth+x,dim);
+    }
+    // High bit
+    shifted = Cshift(in,dim,depth);
+    for(int x=0;x<depth;x++){
+      InsertSliceLocal(shifted,padded,local[dim]-depth+x,depth+local[dim]+x,dim);
+    }
+    // Low bit
+    shifted = Cshift(in,dim,-depth);
+    for(int x=0;x<depth;x++){
+      InsertSliceLocal(shifted,padded,x,x,dim);
+    }
+    return padded;
+  }
+
+};
+ 
+
+NAMESPACE_END(Grid);
+

Grid/qcd/QCD.h

@@ -104,6 +104,7 @@ template<typename vtype> using iSpinMatrix                = iScalar<iMatrix<iSca
 template<typename vtype> using iColourMatrix              = iScalar<iScalar<iMatrix<vtype, Nc> > > ;
 template<typename vtype> using iSpinColourMatrix          = iScalar<iMatrix<iMatrix<vtype, Nc>, Ns> >;
 template<typename vtype> using iLorentzColourMatrix       = iVector<iScalar<iMatrix<vtype, Nc> >, Nd > ;
+template<typename vtype> using iLorentzComplex            = iVector<iScalar<iScalar<vtype> >, Nd > ;
 template<typename vtype> using iDoubleStoredColourMatrix  = iVector<iScalar<iMatrix<vtype, Nc> >, Nds > ;
 template<typename vtype> using iSpinVector                = iScalar<iVector<iScalar<vtype>, Ns> >;
 template<typename vtype> using iColourVector              = iScalar<iScalar<iVector<vtype, Nc> > >;
@@ -178,6 +179,15 @@ typedef iLorentzColourMatrix<vComplexF>  vLorentzColourMatrixF;
 typedef iLorentzColourMatrix<vComplexD>  vLorentzColourMatrixD;
 typedef iLorentzColourMatrix<vComplexD2> vLorentzColourMatrixD2;
 
+// LorentzComplex
+typedef iLorentzComplex<Complex  > LorentzComplex;
+typedef iLorentzComplex<ComplexF > LorentzComplexF;
+typedef iLorentzComplex<ComplexD > LorentzComplexD;
+
+typedef iLorentzComplex<vComplex > vLorentzComplex;
+typedef iLorentzComplex<vComplexF> vLorentzComplexF;
+typedef iLorentzComplex<vComplexD> vLorentzComplexD;
+
 // DoubleStored gauge field
 typedef iDoubleStoredColourMatrix<Complex  > DoubleStoredColourMatrix;
 typedef iDoubleStoredColourMatrix<ComplexF > DoubleStoredColourMatrixF;
@@ -307,6 +317,10 @@ typedef Lattice<vLorentzColourMatrixF>  LatticeLorentzColourMatrixF;
 typedef Lattice<vLorentzColourMatrixD>  LatticeLorentzColourMatrixD;
 typedef Lattice<vLorentzColourMatrixD2> LatticeLorentzColourMatrixD2;
 
+typedef Lattice<vLorentzComplex>  LatticeLorentzComplex;
+typedef Lattice<vLorentzComplexF> LatticeLorentzComplexF;
+typedef Lattice<vLorentzComplexD> LatticeLorentzComplexD;
+
 // DoubleStored gauge field
 typedef Lattice<vDoubleStoredColourMatrix>   LatticeDoubleStoredColourMatrix;
 typedef Lattice<vDoubleStoredColourMatrixF>  LatticeDoubleStoredColourMatrixF;

Grid/stencil/GeneralLocalStencil.h

@@ -123,7 +123,7 @@ public:
 	  }
 	  if ( permute_slice ) {
 	    int ptype       =grid->PermuteType(d);
-	    uint8_t mask    =grid->Nsimd() >> (ptype + 1);		
+	    uint8_t mask    =0x1<<ptype;
 	    SE._permute    |= mask;
 	  }
 	}	

documentation/David_notes.txt

@@ -0,0 +1,90 @@
+Branch: develop
+
+Files:
+
+Grid/lattice/PaddedCell.h -- Halo exchange
+tests/Test_general_stencil.cc -- test local off axis stencil addressing
+tests/debug/Test_padded_cell.cc -- test PaddedCell halo exchange and the General local stencil  by computing ALL plaquettes on lattice
+
+Functionality:
+
+-- extend a lattice field:
+Grid/lattice/PaddedCell.h
+
+// Constructor
+  PaddedCell(int _depth,GridCartesian *_grid)
+
+// Expand a field "in" to depth "d"
+  template<class vobj>
+  inline Lattice<vobj> Exchange(Lattice<vobj> &in)
+  
+// Take the "apple core" of in to a smaller local volume
+  template<class vobj>
+  inline Lattice<vobj> Extract(Lattice<vobj> &in)
+
+-- Plaquette test:
+tests/debug/Test_padded_cell.cc
+  /////////////////////////////////////////////////
+  // Create a padded cell of extra padding depth=1
+  /////////////////////////////////////////////////
+  int depth = 1;
+  PaddedCell Ghost(depth,&GRID);
+  LatticeGaugeField Ughost = Ghost.Exchange(Umu);
+
+///// Array for the site plaquette
+  GridBase *GhostGrid = Ughost.Grid();
+  LatticeComplex gplaq(GhostGrid); 
+
+  std::vector<Coordinate> shifts;
+  for(int mu=0;mu<Nd;mu++){
+    for(int nu=mu+1;nu<Nd;nu++){
+  
+      //    Umu(x) Unu(x+mu) Umu^dag(x+nu) Unu^dag(x)
+      Coordinate shift_0(Nd,0);
+      Coordinate shift_mu(Nd,0); shift_mu[mu]=1;
+      Coordinate shift_nu(Nd,0); shift_nu[nu]=1;
+      shifts.push_back(shift_0);
+      shifts.push_back(shift_mu);
+      shifts.push_back(shift_nu);
+      shifts.push_back(shift_0);
+    }
+  }
+  GeneralLocalStencil gStencil(GhostGrid,shifts);
+
+  gplaq=Zero();
+  {
+    autoView( gp_v , gplaq, CpuWrite);
+    autoView( t_v , trplaq, CpuRead);
+    autoView( U_v , Ughost, CpuRead);
+    for(int ss=0;ss<gp_v.size();ss++){
+      int s=0;
+      for(int mu=0;mu<Nd;mu++){
+	for(int nu=mu+1;nu<Nd;nu++){
+
+	  auto SE0 = gStencil.GetEntry(s+0,ss);
+	  auto SE1 = gStencil.GetEntry(s+1,ss);
+	  auto SE2 = gStencil.GetEntry(s+2,ss);
+	  auto SE3 = gStencil.GetEntry(s+3,ss);
+	
+	  int o0 = SE0->_offset;
+	  int o1 = SE1->_offset;
+	  int o2 = SE2->_offset;
+	  int o3 = SE3->_offset;
+	  
+	  auto U0 = U_v[o0](mu);
+	  auto U1 = U_v[o1](nu);
+	  auto U2 = adj(U_v[o2](mu));
+	  auto U3 = adj(U_v[o3](nu));
+
+	  gpermute(U0,SE0->_permute);
+	  gpermute(U1,SE1->_permute);
+	  gpermute(U2,SE2->_permute);
+	  gpermute(U3,SE3->_permute);
+	  
+	  gp_v[ss]() =gp_v[ss]() + trace( U0*U1*U2*U3 );
+	  s=s+4;
+	}
+      }
+    }
+  }
+  cplaq = Ghost.Extract(gplaq);

systems/mac-arm/config-command-mpi

@@ -1,2 +1,2 @@
-CXX=mpicxx-openmpi-mp CXXFLAGS=-I/opt/local/include/ LDFLAGS=-L/opt/local/lib/ ../../configure --enable-simd=GEN --enable-debug --enable-comms=mpi --enable-unified=yes
-
+BREW=/opt/local/
+MPICXX=mpicxx CXX=c++-12 ../../configure --enable-simd=GEN --enable-comms=mpi-auto --enable-unified=yes --prefix $HOME/QCD/GridInstallOpt --with-lime=/Users/peterboyle/QCD/SciDAC/install/ --with-openssl=$BREW

tests/Test_general_stencil.cc

@@ -115,6 +115,7 @@ int main(int argc, char ** argv)
 	  if (SE->_permute & 0x2 ) { permute(check[i],tmp,1); tmp=check[i];}
 	  if (SE->_permute & 0x4 ) { permute(check[i],tmp,2); tmp=check[i];}
 	  if (SE->_permute & 0x8 ) { permute(check[i],tmp,3); tmp=check[i];}
+	  //	  std::cout<<GridLogMessage<<"stencil["<<i<<"] "<< check[i]<< " perm "<<(uint32_t)SE->_permute <<std::endl;
 	}
 
 	Real nrmC = norm2(Check);
@@ -138,18 +139,17 @@ int main(int argc, char ** argv)
 	  ddiff = check -bar;
 	  diff =norm2(ddiff);
 	  if ( diff > 0){
-	    std::cout <<"Coor (" << coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]
-		      <<") " <<check<<" vs "<<bar<<std::endl;
+	    std::cout <<"Diff at Coor (" << coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]
+		      <<") stencil " <<check<<" vs cshift "<<bar<<std::endl;
 	  }
 
-
 	}}}}
 
 	if (nrm > 1.0e-4) {
 	  autoView( check , Check, CpuRead);
 	  autoView(   bar ,   Bar, CpuRead);
 	  for(int i=0;i<check.size();i++){
-	    std::cout << i<<" Check "<<check[i]<< "\n"<<i<<" Bar "<<bar[i]<<std::endl;
+	    std::cout << i<<" ERROR Check \n"<<check[i]<< "\n"<<i<<" Bar \n"<<bar[i]<<std::endl;
 	  }
 	}
 	if (nrm > 1.0e-4) exit(-1);

tests/core/Test_lie_generators.cc

@@ -63,7 +63,9 @@ int main(int argc, char** argv) {
   std::cout << "Dimension of adjoint representation: "<< SU2Adjoint::Dimension << std::endl;
 
   // guard as this code fails to compile for Nc != 3
-#if (Nc == 3)
+#if 1
+
+  std::cout << " Printing  Adjoint Generators"<< std::endl;
     
   SU2Adjoint::printGenerators();
   SU2::testGenerators();
@@ -148,10 +150,33 @@ int main(int argc, char** argv) {
     typename AdjointRep<Nc>::LatticeMatrix Vrmu = peekLorentz(Vr,mu);
     pokeLorentz(UrVr,Urmu*Vrmu, mu);
   }
-    
+
+  typedef typename SU_Adjoint<Nc>::AMatrix AdjointMatrix;
   typename AdjointRep<Nc>::LatticeField Diff_check = UVr - UrVr;
   std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Adjoint representation) : " << norm2(Diff_check) << std::endl;
-    
+
+  std::cout << GridLogMessage << "****************************************** " << std::endl;
+  std::cout << GridLogMessage << " MAP BETWEEN FUNDAMENTAL AND ADJOINT CHECK " << std::endl;
+  std::cout << GridLogMessage << "****************************************** " << std::endl;
+  for(int a=0;a<Nc*Nc-1;a++){
+  for(int b=0;b<Nc*Nc-1;b++){
+  for(int c=0;c<Nc*Nc-1;c++){
+    ColourMatrix Ta;
+    ColourMatrix Tb;
+    ColourMatrix Tc;
+    SU3::generator(a, Ta);
+    SU3::generator(b, Tb);
+    SU3::generator(c, Tc);
+    AdjointMatrix TRa;
+    SU3Adjoint::generator(a,TRa);
+    Complex tr1 = trace ( Tc * ( Ta*Tb-Tb*Ta)); // i/2 fabc
+    Complex tr2 = TRa()()(b,c) * Complex(0,1);
+    std::cout << " 2 Tr( Tc[Ta,Tb]) " << 2.0*tr1<<std::endl;
+    std::cout << " - TRa_bc " << tr2<<std::endl;
+    assert(abs( (2.0*tr1-tr2) ) < 1.0e-7);
+    std::cout << "------------------"<<std::endl;
+  }}}
+  
   // Check correspondence of algebra and group transformations
   // Create a random vector
   SU3::LatticeAlgebraVector h_adj(grid);

tests/debug/Test_padded_cell.cc

@@ -0,0 +1,184 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_padded_cell.cc
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#include <Grid/Grid.h>
+#include <Grid/lattice/PaddedCell.h>
+#include <Grid/stencil/GeneralLocalStencil.h>
+
+using namespace std;
+using namespace Grid;
+
+template<class vobj> void gpermute(vobj & inout,int perm){
+  vobj tmp=inout;
+  if (perm & 0x1 ) { permute(inout,tmp,0); tmp=inout;}
+  if (perm & 0x2 ) { permute(inout,tmp,1); tmp=inout;}
+  if (perm & 0x4 ) { permute(inout,tmp,2); tmp=inout;}
+  if (perm & 0x8 ) { permute(inout,tmp,3); tmp=inout;}
+}
+  
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  Coordinate latt_size  = GridDefaultLatt();
+  Coordinate simd_layout= GridDefaultSimd(Nd,vComplexD::Nsimd());
+  Coordinate mpi_layout = GridDefaultMpi();
+  std::cout << " mpi "<<mpi_layout<<std::endl;
+  std::cout << " simd "<<simd_layout<<std::endl;
+  std::cout << " latt "<<latt_size<<std::endl;
+  GridCartesian GRID(latt_size,simd_layout,mpi_layout);
+
+  GridParallelRNG   pRNG(&GRID);
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+  LatticeGaugeField Umu(&GRID);
+
+  SU<Nc>::HotConfiguration(pRNG,Umu);
+
+  Real plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
+  LatticeComplex trplaq(&GRID);
+
+  std::vector<LatticeColourMatrix> U(Nd, Umu.Grid());
+  for (int mu = 0; mu < Nd; mu++) {
+    U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
+  }
+
+  std::cout << GridLogMessage << " Average plaquette "<<plaq<<std::endl;
+
+  LatticeComplex cplaq(&GRID); cplaq=Zero();
+
+  /////////////////////////////////////////////////
+  // Create a padded cell of extra padding depth=1
+  /////////////////////////////////////////////////
+  int depth = 1;
+  PaddedCell Ghost(depth,&GRID);
+  LatticeGaugeField Ughost = Ghost.Exchange(Umu);
+
+  ///////////////////////////////////////////////////////////////////
+  // Temporary debug Hack for single rank sim:
+  // Check the contents of the cell are periodcally replicated
+  // In future ONLY pad those dimensions that are not local to node
+  ///////////////////////////////////////////////////////////////////
+#if 0
+  {
+    double diff=0;
+    double n=0;
+  {
+    autoView( Ug_v , Ughost, CpuRead);
+    autoView( Ul_v , Umu   , CpuRead);
+  for(int x=0;x<latt_size[0]+2;x++){
+  for(int y=0;y<latt_size[1]+2;y++){
+  for(int z=0;z<latt_size[2]+2;z++){
+  for(int t=0;t<latt_size[3]+2;t++){
+    int lx=(x-1+latt_size[0])%latt_size[0];
+    int ly=(y-1+latt_size[1])%latt_size[1];
+    int lz=(z-1+latt_size[2])%latt_size[2];
+    int lt=(t-1+latt_size[3])%latt_size[3];
+    Coordinate gcoor({x,y,z,t});
+    Coordinate lcoor({lx,ly,lz,lt});
+    LorentzColourMatrix g;
+    LorentzColourMatrix l;
+    peekLocalSite(g,Ug_v,gcoor);
+    peekLocalSite(l,Ul_v,lcoor);
+    g=g-l;
+    assert(norm2(g)==0);
+    diff = diff + norm2(g);
+    n = n + norm2(l);
+  }}}}
+  }
+  std::cout << "padded field check diff "<< diff <<" / "<< n<<std::endl;
+  std::cout << norm2(Ughost)<< " " << norm2(Umu)<<std::endl;
+  }
+#endif
+
+  ///// Array for the site plaquette
+  GridBase *GhostGrid = Ughost.Grid();
+  LatticeComplex gplaq(GhostGrid); 
+  
+  std::vector<Coordinate> shifts;
+  for(int mu=0;mu<Nd;mu++){
+    for(int nu=mu+1;nu<Nd;nu++){
+  
+      //    Umu(x) Unu(x+mu) Umu^dag(x+nu) Unu^dag(x)
+      Coordinate shift_0(Nd,0);
+      Coordinate shift_mu(Nd,0); shift_mu[mu]=1;
+      Coordinate shift_nu(Nd,0); shift_nu[nu]=1;
+      shifts.push_back(shift_0);
+      shifts.push_back(shift_mu);
+      shifts.push_back(shift_nu);
+      shifts.push_back(shift_0);
+    }
+  }
+  GeneralLocalStencil gStencil(GhostGrid,shifts);
+
+  gplaq=Zero();
+  {
+    autoView( gp_v , gplaq, CpuWrite);
+    autoView( t_v , trplaq, CpuRead);
+    autoView( U_v , Ughost, CpuRead);
+    for(int ss=0;ss<gp_v.size();ss++){
+      int s=0;
+      for(int mu=0;mu<Nd;mu++){
+	for(int nu=mu+1;nu<Nd;nu++){
+
+	  auto SE0 = gStencil.GetEntry(s+0,ss);
+	  auto SE1 = gStencil.GetEntry(s+1,ss);
+	  auto SE2 = gStencil.GetEntry(s+2,ss);
+	  auto SE3 = gStencil.GetEntry(s+3,ss);
+	
+	  int o0 = SE0->_offset;
+	  int o1 = SE1->_offset;
+	  int o2 = SE2->_offset;
+	  int o3 = SE3->_offset;
+	  
+	  auto U0 = U_v[o0](mu);
+	  auto U1 = U_v[o1](nu);
+	  auto U2 = adj(U_v[o2](mu));
+	  auto U3 = adj(U_v[o3](nu));
+
+	  gpermute(U0,SE0->_permute);
+	  gpermute(U1,SE1->_permute);
+	  gpermute(U2,SE2->_permute);
+	  gpermute(U3,SE3->_permute);
+	  
+	  gp_v[ss]() =gp_v[ss]() + trace( U0*U1*U2*U3 );
+	  s=s+4;
+	}
+      }
+    }
+  }
+  cplaq = Ghost.Extract(gplaq);
+  RealD vol = cplaq.Grid()->gSites();
+  RealD faces = (Nd * (Nd-1))/2;
+  auto p = TensorRemove(sum(cplaq));
+  auto result = p.real()/vol/faces/Nc;
+
+  std::cout << GridLogMessage << " Average plaquette via padded cell "<<result<<std::endl;
+  std::cout << GridLogMessage << " Diff "<<result-plaq<<std::endl;
+  
+  assert(fabs(result-plaq)<1.0e-8);
+  Grid_finalize();
+}