Integrated Lebesgue code and been playing with alternate implementations of the wilson dop without

any particular success in increasing the performance.

lib/stencil/Grid_stencil_common.cc

@@ -2,6 +2,96 @@
 
 namespace Grid {
 
+
+
+void CartesianStencil::LebesgueOrder(void) 
+{
+  _LebesgueReorder.resize(0);
+  
+  // Align up dimensions to power of two.
+  const StencilInteger one=1;
+  StencilInteger ND = _grid->_ndimension;
+  std::vector<StencilInteger> dims(ND);
+  std::vector<StencilInteger> adims(ND);
+  std::vector<std::vector<StencilInteger> > bitlist(ND);
+
+
+  for(StencilInteger mu=0;mu<ND;mu++){
+    dims[mu] = _grid->_rdimensions[mu];
+    assert ( dims[mu] != 0 );
+    adims[mu] = alignup(dims[mu]);
+  }
+
+  // List which bits of padded volume coordinate contribute; this strategy 
+  // i) avoids recursion 
+  // ii) has loop lengths at most the width of a 32 bit word.
+  int sitebit=0;
+  int split=24;
+  for(int mu=0;mu<ND;mu++){   // mu 0 takes bit 0; mu 1 takes bit 1 etc...
+    for(int bit=0;bit<split;bit++){
+    StencilInteger mask = one<<bit;
+      if ( mask&(adims[mu]-1) ){
+	bitlist[mu].push_back(sitebit);
+	sitebit++;
+      }
+    }
+  }
+  for(int bit=split;bit<32;bit++){
+    StencilInteger mask = one<<bit;
+    for(int mu=0;mu<ND;mu++){   // mu 0 takes bit 0; mu 1 takes bit 1 etc...
+      if ( mask&(adims[mu]-1) ){
+	bitlist[mu].push_back(sitebit);
+	sitebit++;
+      }
+    }
+  }
+
+  // Work out padded and unpadded volumes
+  StencilInteger avol = 1;
+  for(int mu=0;mu<ND;mu++) avol = avol * adims[mu];
+
+  StencilInteger vol = 1;
+  for(int mu=0;mu<ND;mu++) vol = vol * dims[mu];
+  
+  // Loop over padded volume, following Lebesgue curve
+  // We interleave the bits from sequential "mu".
+  std::vector<StencilInteger> ax(ND);
+  
+  for(StencilInteger asite=0;asite<avol;asite++){
+
+    // Start with zero and collect bits
+    for(int mu=0;mu<ND;mu++) ax[mu] = 0;
+
+    int contained = 1;
+    for(int mu=0;mu<ND;mu++){
+
+      // Build the coordinate on the aligned volume
+      for(int bit=0;bit<bitlist[mu].size();bit++){
+	int sbit=bitlist[mu][bit];
+
+	if(asite&(one<<sbit)){
+	  ax[mu]|=one<<bit;
+	}
+      }
+
+      // Is it contained in original box
+      if ( ax[mu]>dims[mu]-1 ) contained = 0;
+
+    }
+
+    if ( contained ) {
+      int site = ax[0]
+	+        dims[0]*ax[1]
+        +dims[0]*dims[1]*ax[2]
+        +dims[0]*dims[1]*dims[2]*ax[3];
+
+      _LebesgueReorder.push_back(site);
+    }
+  }
+
+  assert( _LebesgueReorder.size() == vol );
+}
+
   CartesianStencil::CartesianStencil(GridBase *grid,
 				     int npoints,
 				     int checkerboard,
@@ -20,6 +110,8 @@ namespace Grid {
       _unified_buffer_size=0;
       _request_count =0;
 
+      LebesgueOrder();
+
       int osites  = _grid->oSites();
 
       for(int i=0;i<npoints;i++){