T-direction terms done

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

Grid/cartesian/CartesianCrossIcosahedron.h

@@ -41,8 +41,17 @@ enum NorthSouth {
   North = 1,
   South = 0
 };
+enum IcoshedralDirections {
+  IcosahedronPatchX = 0,
+  IcosahedronPatchY = 1,
+  IcosahedronPatchDiagonal=2,
+  NumIcosahedralPolarizations
+};
 
-const int num_icosahedron_tiles = 10;
+const int IcosahedralPatches = 10;
+const int HemiPatches=IcosahedralPatches/2;
+const int NorthernHemisphere = HemiPatches;
+const int SouthernHemisphere = 0;
 
 class GridCartesianCrossIcosahedron: public GridCartesian {
 
@@ -86,9 +95,12 @@ public:
 
     assert(simd_layout[0]==1); // Force simd into perpendicular dimensions
     assert(simd_layout[1]==1); // to avoid pole storage complexity interacting with SIMD.
-    assert(dimensions[_ndimension-1]==num_icosahedron_tiles);
+    assert(dimensions[_ndimension-1]==IcosahedralPatches);
     assert(processor_grid[_ndimension-1]<=2); // Keeps the patches that need a pole on the same node
 
+    // Save a copy of the basic cartesian initialisation volume
+    cartesianOsites = this->_osites;
+
     // allocate the pole storage if we are seeking vertex domain data
     if ( meshType == IcosahedralVertices ) {
       InitPoles();
@@ -106,14 +118,46 @@ public:
   int southPoleOsite;
   int northPoleOsites;
   int southPoleOsites;
+  int cartesianOsites;
 
-  virtual int Icosahedral(void)     override { return 1;}
+  virtual int isIcosahedral(void)           override { return 1;}
-  virtual int ownsNorthPole(void)   const override { return hasNorthPole; };
+  virtual int isIcosahedralVertex(void)     override { return meshType==IcosahedralVertices;}
+  virtual int isIcosahedralEdge  (void)     override { return meshType==IcosahedralEdges;}
   virtual int NorthPoleOsite(void)  const override { return northPoleOsite; };
   virtual int NorthPoleOsites(void) const override { return northPoleOsites; };
-  virtual int ownsSouthPole(void)   const override { return hasSouthPole; };
   virtual int SouthPoleOsite(void)  const override { return southPoleOsite; };
   virtual int SouthPoleOsites(void) const override { return southPoleOsites; };
+  virtual int ownsNorthPole(void)   const override { return hasNorthPole; };
+  virtual int ownsSouthPole(void)   const override { return hasSouthPole; };
+  virtual int CartesianOsites(void) const override { return cartesianOsites; };
+  virtual int64_t PoleIdxForOcoor(Coordinate &Coor) override
+  {
+    // Work out the pole_osite. Pick the higher dims
+    Coordinate rdims;
+    Coordinate ocoor;
+    int64_t pole_idx;
+    int Ndm1 = this->Nd()-1;
+    for(int d=2;d<Ndm1;d++){
+      int dd=d-2;
+      rdims.push_back(this->_rdimensions[d]);
+      ocoor.push_back(Coor[d]%this->_rdimensions[d]);
+    }
+    Lexicographic::IndexFromCoor(ocoor,pole_idx,rdims);
+    return pole_idx;
+  }
+  virtual int64_t PoleSiteForOcoor(Coordinate &Coor) override
+  {
+    int Ndm1 = this->Nd()-1;
+    int64_t pole_idx = this->PoleIdxForOcoor(Coor);
+    int64_t pole_osite;
+    if ( Coor[Ndm1] >= HemiPatches ) {
+      pole_osite = pole_idx + this->NorthPoleOsite();
+    } else {
+      pole_osite = pole_idx + this->SouthPoleOsite();
+    }
+    return pole_osite;
+  }
+
 
   void InitPoles(void)
   {
@@ -185,11 +229,11 @@ public:
       northPoleOsites=0;
       northPoleOsite=0;
     }
-    std::cout << "Icosahedral vertex field volume " << this->_osites<<std::endl;
+    std::cout << GridLogDebug<<"Icosahedral vertex field volume " << this->_osites<<std::endl;
-    std::cout << "Icosahedral south pole offset   " << this->southPoleOsite<<std::endl;
+    std::cout << GridLogDebug<<"Icosahedral south pole offset   " << this->southPoleOsite<<std::endl;
-    std::cout << "Icosahedral north pole offset   " << this->northPoleOsite<<std::endl;
+    std::cout << GridLogDebug<<"Icosahedral north pole offset   " << this->northPoleOsite<<std::endl;
-    std::cout << "Icosahedral south pole size     " << this->southPoleOsites<<std::endl;
+    std::cout << GridLogDebug<<"Icosahedral south pole size     " << this->southPoleOsites<<std::endl;
-    std::cout << "Icosahedral north pole size     " << this->northPoleOsites<<std::endl;
+    std::cout << GridLogDebug<<"Icosahedral north pole size     " << this->northPoleOsites<<std::endl;
   };
 
 };

Grid/cartesian/Cartesian_base.h

@@ -87,13 +87,18 @@ public:
 public:
 
   // Icosahedral decisions
-  virtual int Icosahedral(void) { return 0;}
+  virtual int isIcosahedral(void) { return 0;}
+  virtual int isIcosahedralVertex(void) { return 0;}
+  virtual int isIcosahedralEdge  (void) { return 0;}
   virtual int ownsNorthPole(void) const { return 0; };
   virtual int ownsSouthPole(void) const { return 0; };
   virtual int NorthPoleOsite(void) const { return 0; };
   virtual int SouthPoleOsite(void) const { return 0; };
   virtual int NorthPoleOsites(void) const { std::cout << "base osites" <<std::endl;return 0; };
   virtual int SouthPoleOsites(void) const { std::cout << "base osites" <<std::endl;return 0; };
+  virtual int CartesianOsites(void) const { return this->oSites(); };
+  virtual int64_t PoleIdxForOcoor(Coordinate &Coor) { return 0;};
+  virtual int64_t PoleSiteForOcoor(Coordinate &Coor){ return 0;}
 
   ////////////////////////////////////////////////////////////////
   // Checkerboarding interface is virtual and overridden by 

Grid/cshift/Cshift_mpi.h

@@ -34,7 +34,7 @@ NAMESPACE_BEGIN(Grid);
 const int Cshift_verbose=0;
 template<class vobj> Lattice<vobj> Cshift(const Lattice<vobj> &rhs,int dimension,int shift)
 {
-  assert(!rhs.Grid()->Icosahedral());
+  assert(!rhs.Grid()->isIcosahedral());
   
   typedef typename vobj::vector_type vector_type;
   typedef typename vobj::scalar_type scalar_type;

Grid/cshift/Cshift_none.h

@@ -30,7 +30,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 NAMESPACE_BEGIN(Grid);
 template<class vobj> Lattice<vobj> Cshift(const Lattice<vobj> &rhs,int dimension,int shift)
 {
-  assert(!rhs.Grid()->Icosahedral());
+  assert(!rhs.Grid()->isIcosahedral());
   Lattice<vobj> ret(rhs.Grid());
   ret.Checkerboard() = rhs.Grid()->CheckerBoardDestination(rhs.Checkerboard(),shift,dimension);
   Cshift_local(ret,rhs,dimension,shift);

Grid/lattice/Lattice_base.h

@@ -391,7 +391,7 @@ template<class vobj> std::ostream& operator<< (std::ostream& stream, const Latti
     stream<<"]\t";
     stream<<ss<<std::endl;
   }
-  if ( o.Grid()->Icosahedral() ) {
+  if ( o.Grid()->isIcosahedralVertex() ) {
     uint64_t psites=1;
     Coordinate perpdims;
     for(int d=2;d<o.Grid()->_ndimension-1;d++){

Grid/lattice/Lattice_coordinate.h

@@ -40,7 +40,7 @@ template<class iobj> inline void LatticeCoordinate(Lattice<iobj> &l,int mu)
   int Nsimd = grid->iSites();
 
   int cartesian_vol = grid->oSites();
-  if ( grid->Icosahedral() ) {
+  if ( grid->isIcosahedral() ) {
     cartesian_vol = cartesian_vol - grid->NorthPoleOsites()-grid->SouthPoleOsites();
   }
   {
@@ -58,7 +58,7 @@ template<class iobj> inline void LatticeCoordinate(Lattice<iobj> &l,int mu)
       });
   }
 
-  if (grid->Icosahedral()) {
+  if (grid->isIcosahedralVertex()) {
     uint64_t psites=1;
     Coordinate perpdims;
     typename iobj::scalar_object ss;
@@ -95,7 +95,9 @@ template<class iobj> inline void LatticePole(Lattice<iobj> &l,NorthSouth pole)
 
   l=Zero();
 
-  if (grid->Icosahedral()) {
+  assert(grid->isIcosahedralVertex());
+  
+  if (grid->isIcosahedralVertex()) {
     uint64_t psites=1;
     Coordinate perpdims;
     sobj ss;

Grid/lattice/Lattice_peekpoke.h

@@ -159,7 +159,8 @@ void peekPole(sobj &s,const Lattice<vobj> &l,const Coordinate &orthog,NorthSouth
   
   GridBase *grid=l.Grid();
 
-  assert(grid->Icosahedral());
+  assert(grid->isIcosahedral());
+  assert(grid->isIcosahedralVertex());
 
   int Nsimd = grid->Nsimd();
 
@@ -219,7 +220,8 @@ void pokePole(const sobj &s,Lattice<vobj> &l,const Coordinate &orthog,NorthSouth
 {
   GridBase *grid=l.Grid();
 
-  assert(grid->Icosahedral());
+  assert(grid->isIcosahedral());
+  assert(grid->isIcosahedralVertex());
 
   grid->Broadcast(grid->BossRank(),s);
 
@@ -279,6 +281,131 @@ void pokePole(const sobj &s,Lattice<vobj> &l,const Coordinate &orthog,NorthSouth
   return;
 };
 
+
+template<class vobj,class sobj>
+void peekLocalPole(sobj &s,const Lattice<vobj> &l,const Coordinate &orthog,NorthSouth isNorth)
+{
+  s=Zero();
+  
+  GridBase *grid=l.Grid();
+
+  assert(grid->isIcosahedral());
+  assert(grid->isIcosahedralVertex());
+
+  int Nsimd = grid->Nsimd();
+
+  int rank;
+
+  int Ndm1         = grid->_ndimension-1;
+  Coordinate pgrid = grid->ProcessorGrid();
+  const int xdim=0;
+  const int ydim=1;
+  const int pdim=Ndm1;
+
+  int64_t pole_osite;
+  int64_t pole_isite;
+  Coordinate rdims;
+  Coordinate idims;
+  Coordinate ocoor;
+  Coordinate icoor;
+  //  Coordinate pcoor(grid->_ndimension);
+  for(int d=2;d<Ndm1;d++){
+    int dd=d-2;
+    rdims.push_back(grid->_rdimensions[d]);
+    idims.push_back(grid->_simd_layout[d]);
+    icoor.push_back((orthog[dd]%grid->_ldimensions[d])/grid->_rdimensions[d]);
+    ocoor.push_back(orthog[dd]%grid->_rdimensions[d]);
+    //    pcoor[d] = orthog[dd]/grid->_ldimensions[d];
+  }
+  Lexicographic::IndexFromCoor(ocoor,pole_osite,rdims);
+  Lexicographic::IndexFromCoor(icoor,pole_isite,idims);
+  
+  int64_t osite;
+  if(isNorth == North){
+    //    pcoor[xdim] = 0;
+    //    pcoor[ydim] = pgrid[ydim]-1;
+    //    pcoor[Ndm1] = pgrid[Ndm1]-1;
+    osite = pole_osite + grid->NorthPoleOsite();
+    assert(grid->ownsNorthPole());
+  } else {
+    //    pcoor[xdim] = pgrid[xdim]-1;
+    //    pcoor[ydim] = 0;
+    //    pcoor[Ndm1] = 0;
+    osite = pole_osite + grid->SouthPoleOsite();
+    assert(grid->ownsSouthPole());
+  }
+
+  ExtractBuffer<sobj> buf(Nsimd);
+  autoView( l_v , l, CpuWrite);
+  extract(l_v[osite],buf);
+  s = buf[pole_isite];
+
+  return;
+};
+template<class vobj,class sobj>
+void pokeLocalPole(const sobj &s,Lattice<vobj> &l,const Coordinate &orthog,NorthSouth isNorth)
+{
+  GridBase *grid=l.Grid();
+
+  assert(grid->isIcosahedral());
+  assert(grid->isIcosahedralVertex());
+
+  int Nsimd = grid->Nsimd();
+  int rank;
+  int Ndm1         = grid->_ndimension-1;
+
+  const int xdim=0;
+  const int ydim=1;
+  const int pdim=Ndm1;
+
+  int64_t pole_osite;
+  int64_t pole_isite;
+  Coordinate rdims;
+  Coordinate idims;
+  Coordinate ocoor;
+  Coordinate icoor;
+  //  Coordinate pcoor(grid->_ndimension,0);
+  for(int d=2;d<Ndm1;d++){
+    int dd = d-2;
+    rdims.push_back(grid->_rdimensions[d]);
+    idims.push_back(grid->_simd_layout[d]);
+    icoor.push_back((orthog[dd]%grid->_ldimensions[d])/grid->_rdimensions[d]);
+    ocoor.push_back(orthog[dd]%grid->_rdimensions[d]);
+    //    pcoor[d] = orthog[dd]/grid->_ldimensions[d];
+
+    int o = orthog[dd];
+    int r = grid->_rdimensions[d];
+    int omr = o % r;
+  }
+  Lexicographic::IndexFromCoor(ocoor,pole_osite,rdims);
+  Lexicographic::IndexFromCoor(icoor,pole_isite,idims);
+  
+  int64_t osite;
+  int insert=0;
+  if(isNorth ==North){
+    //    pcoor[xdim] = 0;
+    //    pcoor[ydim] = pgrid[ydim]-1;
+    //    pcoor[Ndm1] = pgrid[Ndm1]-1;
+    osite = pole_osite + grid->NorthPoleOsite();
+    assert(grid->ownsNorthPole());
+  } else {
+    //    pcoor[xdim] = pgrid[xdim]-1;
+    //    pcoor[ydim] = 0;
+    //    pcoor[Ndm1] = 0;
+    osite = pole_osite + grid->SouthPoleOsite();
+    assert(grid->ownsSouthPole());
+  }
+
+  // extract-modify-merge cycle is easiest way and this is not perf critical
+  ExtractBuffer<sobj> buf(Nsimd);
+  autoView( l_v , l, CpuWrite);
+  extract(l_v[osite],buf);
+  buf[pole_isite] = s;
+  merge(l_v[osite],buf);
+  
+  return;
+};
+
 //////////////////////////////////////////////////////////
 // Peek a scalar object from the SIMD array
 //////////////////////////////////////////////////////////

Grid/lattice/Lattice_rng.h

@@ -47,6 +47,19 @@ NAMESPACE_BEGIN(Grid);
 // Allow the RNG state to be less dense than the fine grid
 //////////////////////////////////////////////////////////////
 inline int RNGfillable(GridBase *coarse,GridBase *fine)
+{
+  if ( coarse == fine ) return 1;
+
+  if ( coarse->isIcosahedral()) assert(coarse->isIcosahedralEdge());
+  
+  if ( fine->isIcosahedralVertex() && coarse->isIcosahedralEdge() ) {
+    assert(fine->Nd()==coarse->Nd());
+    for(int d=0;d<fine->Nd();d++){
+      assert(fine->LocalDimensions()[d] == coarse->LocalDimensions()[d]);
+    }
+    return 1;
+  }
+    
   {
     
     int rngdims = coarse->_ndimension;
@@ -74,12 +87,26 @@ inline int RNGfillable(GridBase *coarse,GridBase *fine)
     }
     return multiplicity;
   }
+}
 
   
 // merge of April 11 2017
 // this function is necessary for the LS vectorised field
 inline int RNGfillable_general(GridBase *coarse,GridBase *fine)
 {
+
+  if ( coarse == fine ) return 1;
+
+  if ( coarse->isIcosahedral()) assert(coarse->isIcosahedralEdge());
+  
+  if ( fine->isIcosahedralVertex() && coarse->isIcosahedralEdge() ) {
+    assert(fine->Nd()==coarse->Nd());
+    for(int d=0;d<fine->Nd();d++){
+      assert(fine->LocalDimensions()[d] == coarse->LocalDimensions()[d]);
+    }
+    return 1;
+  }
+
   int rngdims = coarse->_ndimension;
     
   // trivially extended in higher dims, with locality guaranteeing RNG state is local to node
@@ -352,12 +379,12 @@ private:
 public:
   GridBase *Grid(void) const { return _grid; }
   int generator_idx(int os,int is) {
-    return is*_grid->oSites()+os;
+    return (is*_grid->CartesianOsites()+os)%_grid->lSites(); // On the pole sites wrap back to normal generators; Icosahedral hack
   }
 
   GridParallelRNG(GridBase *grid) : GridRNGbase() {
     _grid = grid;
-    _vol  =_grid->iSites()*_grid->oSites();
+    _vol  =_grid->lSites();
 
     _generators.resize(_vol);
     _uniform.resize(_vol,std::uniform_real_distribution<RealD>{0,1});
@@ -381,7 +408,7 @@ public:
 
     int multiplicity = RNGfillable_general(_grid, l.Grid()); // l has finer or same grid
     int Nsimd  = _grid->Nsimd();  // guaranteed to be the same for l.Grid() too
-    int osites = _grid->oSites();  // guaranteed to be <= l.Grid()->oSites() by a factor multiplicity
+    int osites = _grid->CartesianOsites();  // guaranteed to be <= l.Grid()->oSites() by a factor multiplicity, except on Icosahedral
     int words  = sizeof(scalar_object) / sizeof(scalar_type);
 
     autoView(l_v, l, CpuWrite);
@@ -403,7 +430,26 @@ public:
 	merge(l_v[sm], buf);
       }
     });
-    //    });
+
+    /*
+     * Fill in the poles for an Icosahedral vertex mesh
+     */
+    if (l.Grid()->isIcosahedralVertex()) { 
+      int64_t pole_sites=l.Grid()->NorthPoleOsites()+l.Grid()->SouthPoleOsites();
+      int64_t pole_base =l.Grid()->CartesianOsites();
+
+      ExtractBuffer<scalar_object> buf(Nsimd);
+      for (int m = 0; m < pole_sites; m++) {  // Draw from same generator multiplicity times                                                                                                           
+        for (int si = 0; si < Nsimd; si++) {
+          int gdx = 0;
+	  scalar_type *pointer = (scalar_type *)&buf[si];
+          dist[gdx].reset();
+          for (int idx = 0; idx < words; idx++)
+            fillScalar(pointer[idx], dist[gdx], _generators[gdx]);
+        }
+        merge(l_v[pole_base+m], buf);
+      }      
+    }
 
     _time_counter += usecond()- inner_time_counter;
   }

Grid/simd/Simd.h

@@ -252,7 +252,7 @@ inline std::ostream& operator<< (std::ostream& stream, const vComplexF &o){
  
 inline std::ostream& operator<< (std::ostream& stream, const vComplexD &o){
   int nn=vComplexD::Nsimd();
-  std::vector<ComplexD,alignedAllocator<ComplexD> > buf(nn);
+  std::vector<ComplexD> buf(nn);
   vstore(o,&buf[0]);
   stream<<"<";
   for(int i=0;i<nn;i++){
@@ -272,7 +272,7 @@ inline std::ostream& operator<< (std::ostream& stream, const vComplexD2 &o){
 
 inline std::ostream& operator<< (std::ostream& stream, const vRealF &o){
   int nn=vRealF::Nsimd();
-  std::vector<RealF,alignedAllocator<RealF> > buf(nn);
+  std::vector<RealF> buf(nn);
   vstore(o,&buf[0]);
   stream<<"<";
   for(int i=0;i<nn;i++){

Grid/util/Init.cc

@@ -187,8 +187,9 @@ void GridParseLayout(char **argv,int argc,
 		     Coordinate &latt_c,
 		     Coordinate &mpi_c)
 {
-  auto mpi =std::vector<int>(1,Nd);
+  auto mpi =std::vector<int>(Nd,1);
-  auto latt=std::vector<int>(8,Nd);
+  auto latt=std::vector<int>(Nd,8);
+
 
   GridThread::SetMaxThreads();
 
@@ -228,6 +229,9 @@ void GridParseLayout(char **argv,int argc,
   }
   // Copy back into coordinate format
   int nd = mpi.size();
+  //  std::cout << "mpi.size() "<<nd<<std::endl;
+  //  std::cout << "latt.size() "<<latt.size()<<std::endl;
+  //  std::cout << "Nd "<<Nd<<std::endl;
   assert(latt.size()==nd);
   latt_c.resize(nd);
    mpi_c.resize(nd);

configure.ac

@@ -226,7 +226,7 @@ case ${ac_Nd} in
     4)
         AC_DEFINE([Config_Nd],[4],[Gauge field dimension Nd]);;
     *)
-      AC_MSG_ERROR(["Unsupport gauge group choice Nc = ${ac_Nc}"]);;
+      AC_MSG_ERROR(["Unsupport dimension Nd = ${ac_Nd}"]);;
 esac
 
 ############### Symplectic group
@@ -835,6 +835,7 @@ os (target)                 : $target_os
 compiler vendor             : ${ax_cv_cxx_compiler_vendor}
 compiler version            : ${ax_cv_gxx_version}
 ----- BUILD OPTIONS -----------------------------------
+Nd                          : ${ac_Nd}
 Nc                          : ${ac_Nc}
 SIMD                        : ${ac_SIMD}${SIMD_GEN_WIDTH_MSG}
 Threading                   : ${ac_openmp}

systems/mac-arm/config-command

@@ -0,0 +1,12 @@
+MPICXX=mpicxx CXXFLAGS=-I/opt/local/include LDFLAGS=-L/opt/local/lib/ CXX=clang++ ../../configure \
+	      --enable-simd=GEN \
+	      --enable-Nc=1 \
+	      --enable-debug \
+	      --enable-unified=yes \
+	      --prefix $HOME/QCD/GridInstall \
+	      --with-lime=/Users/peterboyle/QCD/SciDAC/install/ \
+	      --with-openssl=$BREW \
+	      --disable-fermion-reps \
+	      --disable-gparity \
+	      --enable-debug
+