Creation mode better implementation

UID fix
Merge branch 'feature/sitmo-skipahead' into develop
2025-06-13 20:57:06 +01:00 · 2017-04-05 02:35:34 +09:00 · 2017-04-05 02:18:15 +09:00 · 2017-04-05 01:24:49 +09:00 · 2017-04-02 23:13:48 +09:00 · 2017-04-02 12:10:51 +09:00
50 changed files with 1345 additions and 311 deletions
--- a/benchmarks/Benchmark_memory_asynch.cc
+++ b/benchmarks/Benchmark_memory_asynch.cc
@ -66,7 +66,8 @@ int main (int argc, char ** argv)

    Vec tsum; tsum = zero;

-    GridParallelRNG          pRNG(&Grid);      pRNG.SeedRandomDevice();
+    GridParallelRNG          pRNG(&Grid);      
+    pRNG.SeedFixedIntegers(std::vector<int>({56,17,89,101}));

    std::vector<double> stop(threads);
    Vector<Vec> sum(threads);
--- a/benchmarks/Benchmark_memory_bandwidth.cc
+++ b/benchmarks/Benchmark_memory_bandwidth.cc
@ -65,7 +65,7 @@ int main (int argc, char ** argv)

      uint64_t Nloop=NLOOP;

-      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedRandomDevice();
+      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});

      LatticeVec z(&Grid); //random(pRNG,z);
      LatticeVec x(&Grid); //random(pRNG,x);
@ -100,7 +100,7 @@ int main (int argc, char ** argv)
      int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);

-      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedRandomDevice();
+      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});

      LatticeVec z(&Grid); //random(pRNG,z);
      LatticeVec x(&Grid); //random(pRNG,x);
@ -138,7 +138,7 @@ int main (int argc, char ** argv)

      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);

-      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedRandomDevice();
+      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});

      LatticeVec z(&Grid); //random(pRNG,z);
      LatticeVec x(&Grid); //random(pRNG,x);
@ -173,7 +173,7 @@ int main (int argc, char ** argv)
      uint64_t Nloop=NLOOP;
      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);

-      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedRandomDevice();
+      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});
      LatticeVec z(&Grid); //random(pRNG,z);
      LatticeVec x(&Grid); //random(pRNG,x);
      LatticeVec y(&Grid); //random(pRNG,y);
--- a/benchmarks/Benchmark_staggered.cc
+++ b/benchmarks/Benchmark_staggered.cc
@ -51,7 +51,7 @@ int main (int argc, char ** argv)
  std::vector<int> seeds({1,2,3,4});
  GridParallelRNG          pRNG(&Grid);
  pRNG.SeedFixedIntegers(seeds);
-  //  pRNG.SeedRandomDevice();
+  //  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});

  typedef typename ImprovedStaggeredFermionR::FermionField FermionField; 
  typename ImprovedStaggeredFermionR::ImplParams params; 
--- a/benchmarks/Benchmark_su3.cc
+++ b/benchmarks/Benchmark_su3.cc
@ -55,7 +55,7 @@ int main (int argc, char ** argv)
      std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
      int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedRandomDevice();
+      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});

      LatticeColourMatrix z(&Grid);// random(pRNG,z);
      LatticeColourMatrix x(&Grid);// random(pRNG,x);
@ -88,7 +88,7 @@ int main (int argc, char ** argv)
      int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];

      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedRandomDevice();
+      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});

      LatticeColourMatrix z(&Grid); //random(pRNG,z);
      LatticeColourMatrix x(&Grid); //random(pRNG,x);
@ -119,7 +119,7 @@ int main (int argc, char ** argv)
      int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];

      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedRandomDevice();
+      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});

      LatticeColourMatrix z(&Grid); //random(pRNG,z);
      LatticeColourMatrix x(&Grid); //random(pRNG,x);
@ -150,7 +150,7 @@ int main (int argc, char ** argv)
      int vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];

      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedRandomDevice();
+      //      GridParallelRNG          pRNG(&Grid);      pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});

      LatticeColourMatrix z(&Grid); //random(pRNG,z);
      LatticeColourMatrix x(&Grid); //random(pRNG,x);
--- a/benchmarks/Benchmark_wilson.cc
+++ b/benchmarks/Benchmark_wilson.cc
@ -69,7 +69,7 @@ int main (int argc, char ** argv)
  std::vector<int> seeds({1,2,3,4});
  GridParallelRNG          pRNG(&Grid);
  pRNG.SeedFixedIntegers(seeds);
-  //  pRNG.SeedRandomDevice();
+  //  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});

  LatticeFermion src   (&Grid); random(pRNG,src);
  LatticeFermion result(&Grid); result=zero;
--- a/configure.ac
+++ b/configure.ac
@ -321,7 +321,7 @@ AM_CONDITIONAL(BUILD_COMMS_NONE,  [ test "${comms_type}X" == "noneX" ])
 ############### RNG selection
 AC_ARG_ENABLE([rng],[AC_HELP_STRING([--enable-rng=ranlux48|mt19937|sitmo],\
 	            [Select Random Number Generator to be used])],\
-	            [ac_RNG=${enable_rng}],[ac_RNG=ranlux48])
+	            [ac_RNG=${enable_rng}],[ac_RNG=sitmo])

 case ${ac_RNG} in
     ranlux48)
@ -401,6 +401,7 @@ AC_CONFIG_FILES(tests/hadrons/Makefile)
 AC_CONFIG_FILES(tests/hmc/Makefile)
 AC_CONFIG_FILES(tests/solver/Makefile)
 AC_CONFIG_FILES(tests/qdpxx/Makefile)
+AC_CONFIG_FILES(tests/testu01/Makefile)
 AC_CONFIG_FILES(benchmarks/Makefile)
 AC_CONFIG_FILES(extras/Makefile)
 AC_CONFIG_FILES(extras/Hadrons/Makefile)
--- a/lib/algorithms/CoarsenedMatrix.h
+++ b/lib/algorithms/CoarsenedMatrix.h
@ -425,7 +425,7 @@ namespace Grid {
 	A[p]=zero;
      }

-      GridParallelRNG  RNG(Grid()); RNG.SeedRandomDevice();
+      GridParallelRNG  RNG(Grid()); RNG.SeedFixedIntegers(std::vector<int>({55,72,19,17,34}));
      Lattice<iScalar<CComplex> > val(Grid()); random(RNG,val);

      Complex one(1.0);
--- a/lib/cartesian/Cartesian_base.h
+++ b/lib/cartesian/Cartesian_base.h
@ -177,9 +177,11 @@ public:
    // Global addressing
    ////////////////////////////////////////////////////////////////
    void GlobalIndexToGlobalCoor(int gidx,std::vector<int> &gcoor){
+      assert(gidx< gSites());
      Lexicographic::CoorFromIndex(gcoor,gidx,_gdimensions);
    }
    void LocalIndexToLocalCoor(int lidx,std::vector<int> &lcoor){
+      assert(lidx<lSites());
      Lexicographic::CoorFromIndex(lcoor,lidx,_ldimensions);
    }
    void GlobalCoorToGlobalIndex(const std::vector<int> & gcoor,int & gidx){
--- a/lib/communicator/Communicator_mpi3.cc
+++ b/lib/communicator/Communicator_mpi3.cc
@ -206,7 +206,7 @@ void CartesianCommunicator::Init(int *argc, char ***argv) {
      sprintf(shm_name,"/Grid_mpi3_shm_%d_%d",GroupRank,r);

      shm_unlink(shm_name);
-      int fd=shm_open(shm_name,O_RDWR|O_CREAT,0660);
+      int fd=shm_open(shm_name,O_RDWR|O_CREAT,0666);
      if ( fd < 0 ) {	perror("failed shm_open");	assert(0);      }
      ftruncate(fd, size);

@ -226,7 +226,7 @@ void CartesianCommunicator::Init(int *argc, char ***argv) {
    
      sprintf(shm_name,"/Grid_mpi3_shm_%d_%d",GroupRank,r);

-      int fd=shm_open(shm_name,O_RDWR,0660);
+      int fd=shm_open(shm_name,O_RDWR,0666);
      if ( fd<0 ) {	perror("failed shm_open");	assert(0);      }

      void * ptr =  mmap(NULL,size, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0);
--- a/lib/lattice/Lattice_rng.h
+++ b/lib/lattice/Lattice_rng.h
@ -30,12 +30,19 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #define GRID_LATTICE_RNG_H

 #include <random>
+
+#ifdef RNG_SITMO
 #include <Grid/sitmo_rng/sitmo_prng_engine.hpp>
+#endif 
+
+#if defined(RNG_SITMO)
+#define RNG_FAST_DISCARD
+#else 
+#undef  RNG_FAST_DISCARD
+#endif

 namespace Grid {

-  //http://nvlpubs.nist.gov/nistpubs/SpecialPublications/NIST.SP.800-90Ar1.pdf ?
-
  //////////////////////////////////////////////////////////////
  // Allow the RNG state to be less dense than the fine grid
  //////////////////////////////////////////////////////////////
@ -65,120 +72,139 @@ namespace Grid {

      multiplicity = multiplicity *fine->_rdimensions[fd] / coarse->_rdimensions[d]; 
    }
-
    return multiplicity;
  }

-  // Wrap seed_seq to give common interface with random_device
-  // Should rather wrap random_device and have a generate
-  class fixedSeed {
-  public:
-
-    typedef std::seed_seq::result_type result_type;
-
-    std::seed_seq src;
-    
-    template<class int_type> fixedSeed(const std::vector<int_type> &seeds) : src(seeds.begin(),seeds.end()) {};
-
-    template< class RandomIt > void generate( RandomIt begin, RandomIt end ) {
-      src.generate(begin,end);
-    }
-
-  };
-
-
-  class deviceSeed {
-  public:
-
-    std::random_device rd;
-
-    typedef std::random_device::result_type result_type;
-    
-    deviceSeed(void) : rd(){};
-
-    template< class RandomIt > void generate( RandomIt begin, RandomIt end ) {
-      for(RandomIt it=begin; it!=end;it++){
-	*it = rd();
-      }
-    }
-  };
-
  // real scalars are one component
-  template<class scalar,class distribution,class generator> void fillScalar(scalar &s,distribution &dist,generator & gen)
+  template<class scalar,class distribution,class generator> 
+  void fillScalar(scalar &s,distribution &dist,generator & gen)
  {
    s=dist(gen);
  }
-  template<class distribution,class generator> void fillScalar(ComplexF &s,distribution &dist, generator &gen)
+  template<class distribution,class generator> 
+  void fillScalar(ComplexF &s,distribution &dist, generator &gen)
  {
    s=ComplexF(dist(gen),dist(gen));
  }
-  template<class distribution,class generator> void fillScalar(ComplexD &s,distribution &dist,generator &gen)
+  template<class distribution,class generator> 
+  void fillScalar(ComplexD &s,distribution &dist,generator &gen)
  {
    s=ComplexD(dist(gen),dist(gen));
  }
  
  class GridRNGbase {
-
  public:
-
-    int _seeded;
    // One generator per site.
    // Uniform and Gaussian distributions from these generators.
 #ifdef RNG_RANLUX
-    typedef uint64_t      RngStateType;
    typedef std::ranlux48 RngEngine;
+    typedef uint64_t      RngStateType;
    static const int RngStateCount = 15;
-#elif RNG_MT19937 
+#endif 
+#ifdef RNG_MT19937 
    typedef std::mt19937 RngEngine;
    typedef uint32_t     RngStateType;
    static const int     RngStateCount = std::mt19937::state_size;
-#elif RNG_SITMO
+#endif
+#ifdef RNG_SITMO
    typedef sitmo::prng_engine 	RngEngine;
    typedef uint64_t    	RngStateType;
    static const int    	RngStateCount = 4;
 #endif
-    std::vector<RngEngine>                             _generators;
-    std::vector<std::uniform_real_distribution<RealD>> _uniform;
-    std::vector<std::normal_distribution<RealD>>       _gaussian;
-    std::vector<std::discrete_distribution<int32_t>>   _bernoulli;

-    void GetState(std::vector<RngStateType> & saved,int gen) {
+    std::vector<RngEngine>                             _generators;
+    std::vector<std::uniform_real_distribution<RealD> > _uniform;
+    std::vector<std::normal_distribution<RealD> >       _gaussian;
+    std::vector<std::discrete_distribution<int32_t> >   _bernoulli;
+    std::vector<std::uniform_int_distribution<uint32_t> > _uid;
+
+    ///////////////////////
+    // support for parallel init
+    ///////////////////////
+#ifdef RNG_FAST_DISCARD
+    static void Skip(RngEngine &eng)
+    {
+      /////////////////////////////////////////////////////////////////////////////////////
+      // Skip by 2^40 elements between successive lattice sites
+      // This goes by 10^12.
+      // Consider quenched updating; likely never exceeding rate of 1000 sweeps
+      // per second on any machine. This gives us of order 10^9 seconds, or 100 years
+      // skip ahead.
+      // For HMC unlikely to go at faster than a solve per second, and 
+      // tens of seconds per trajectory so this is clean in all reasonable cases,
+      // and margin of safety is orders of magnitude.
+      // We could hack Sitmo to skip in the higher order words of state if necessary
+      /////////////////////////////////////////////////////////////////////////////////////
+      uint64_t skip = 0x1; skip = skip<<40;
+      eng.discard(skip);
+    } 
+#endif
+    static RngEngine Reseed(RngEngine &eng)
+    {
+      std::vector<uint32_t> newseed;
+      std::uniform_int_distribution<uint32_t> uid;
+      return Reseed(eng,newseed,uid);
+    }
+    static RngEngine Reseed(RngEngine &eng,std::vector<uint32_t> & newseed,
+			    std::uniform_int_distribution<uint32_t> &uid)
+    {
+      const int reseeds=4;
+      
+      newseed.resize(reseeds);
+      for(int i=0;i<reseeds;i++){
+	newseed[i] = uid(eng);
+      }
+      std::seed_seq sseq(newseed.begin(),newseed.end());
+      return RngEngine(sseq);
+    }    
+
+    void GetState(std::vector<RngStateType> & saved,RngEngine &eng) {
      saved.resize(RngStateCount);
      std::stringstream ss;
-      ss<<_generators[gen];
+      ss<<eng;
      ss.seekg(0,ss.beg);
      for(int i=0;i<RngStateCount;i++){
 	ss>>saved[i];
      }
    }
-    void SetState(std::vector<RngStateType> & saved,int gen){
+    void GetState(std::vector<RngStateType> & saved,int gen) {
+      GetState(saved,_generators[gen]);
+    }
+    void SetState(std::vector<RngStateType> & saved,RngEngine &eng){
      assert(saved.size()==RngStateCount);
      std::stringstream ss;
      for(int i=0;i<RngStateCount;i++){
 	ss<< saved[i]<<" ";
      }
      ss.seekg(0,ss.beg);
-      ss>>_generators[gen];
+      ss>>eng;
    }
+    void SetState(std::vector<RngStateType> & saved,int gen){
+      SetState(saved,_generators[gen]);
+    }
+    void SetEngine(RngEngine &Eng, int gen){
+      _generators[gen]=Eng;
+    }
+    void GetEngine(RngEngine &Eng, int gen){
+      Eng=_generators[gen];
+    }
+    template<class source> void Seed(source &src, int gen)
+    {
+      _generators[gen] = RngEngine(src);
+    }    
  };

  class GridSerialRNG : public GridRNGbase {
  public:

-    // FIXME ... do we require lockstep draws of randoms 
-    // from all nodes keeping seeds consistent.
-    // place a barrier/broadcast in the fill routine
-
    GridSerialRNG() : GridRNGbase() {
      _generators.resize(1);
      _uniform.resize(1,std::uniform_real_distribution<RealD>{0,1});
      _gaussian.resize(1,std::normal_distribution<RealD>(0.0,1.0) );
      _bernoulli.resize(1,std::discrete_distribution<int32_t>{1,1});
-      _seeded=0;
+      _uid.resize(1,std::uniform_int_distribution<uint32_t>() );
    }

-
-
    template <class sobj,class distribution> inline void fill(sobj &l,std::vector<distribution> &dist){

      typedef typename sobj::scalar_type scalar_type;
@ -191,7 +217,7 @@ namespace Grid {
      for(int idx=0;idx<words;idx++){
 	fillScalar(buf[idx],dist[0],_generators[0]);
      }
-      
+
      CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));

    };
@ -250,28 +276,18 @@ namespace Grid {
      CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
    }

-    template<class source> void Seed(source &src)
-    {
-      _generators[0] = RngEngine(src);
-      _seeded=1;
-    }    
-    void SeedRandomDevice(void){
-      deviceSeed src;
-      Seed(src);
-    }
    void SeedFixedIntegers(const std::vector<int> &seeds){
      CartesianCommunicator::BroadcastWorld(0,(void *)&seeds[0],sizeof(int)*seeds.size());
-      fixedSeed src(seeds);
-      Seed(src);
+      std::seed_seq src(seeds.begin(),seeds.end());
+      Seed(src,0);
    }
-
  };

  class GridParallelRNG : public GridRNGbase {
  public:
-
    GridBase *_grid;
    int _vol;
+  public:

    int generator_idx(int os,int is){
      return is*_grid->oSites()+os;
@ -285,15 +301,9 @@ namespace Grid {
      _uniform.resize(_vol,std::uniform_real_distribution<RealD>{0,1});
      _gaussian.resize(_vol,std::normal_distribution<RealD>(0.0,1.0) );
      _bernoulli.resize(_vol,std::discrete_distribution<int32_t>{1,1});
-      _seeded=0;
+      _uid.resize(_vol,std::uniform_int_distribution<uint32_t>() );
    }

-
-
-    //FIXME implement generic IO and create state save/restore
-    //void SaveState(const std::string<char> &file);
-    //void LoadState(const std::string<char> &file);
-
    template <class vobj,class distribution> inline void fill(Lattice<vobj> &l,std::vector<distribution> &dist){

      typedef typename vobj::scalar_object scalar_object;
@ -306,7 +316,6 @@ namespace Grid {
      int     osites=_grid->oSites();
      int words=sizeof(scalar_object)/sizeof(scalar_type);

-
      parallel_for(int ss=0;ss<osites;ss++){

 	std::vector<scalar_object> buf(Nsimd);
@ -329,104 +338,114 @@ namespace Grid {
      }
    };

-    // This loop could be made faster to avoid the Ahmdahl by
-    // i)  seed generators on each timeslice, for x=y=z=0;
-    // ii) seed generators on each z for x=y=0
-    // iii)seed generators on each y,z for x=0
-    // iv) seed generators on each y,z,x 
-    // made possible by physical indexing.
-    template<class source> void Seed(source &src)
-    {
+    void SeedFixedIntegers(const std::vector<int> &seeds){

-      typedef typename source::result_type seed_t;
-      std::uniform_int_distribution<seed_t> uid;
+      // Everyone generates the same seed_seq based on input seeds
+      CartesianCommunicator::BroadcastWorld(0,(void *)&seeds[0],sizeof(int)*seeds.size());

-      int numseed=4;
-      int gsites = _grid->_gsites;
-      std::vector<seed_t> site_init(numseed);
+      std::seed_seq source(seeds.begin(),seeds.end());
+
+      RngEngine master_engine(source);
+
+#ifdef RNG_FAST_DISCARD
+      ////////////////////////////////////////////////
+      // Skip ahead through a single stream.
+      // Applicable to SITMO and other has based/crypto RNGs
+      // Should be applicable to Mersenne Twister, but the C++11
+      // MT implementation does not implement fast discard even though
+      // in principle this is possible
+      ////////////////////////////////////////////////
      std::vector<int> gcoor;
+      int rank,o_idx,i_idx;

+      // Everybody loops over global volume.
+      for(int gidx=0;gidx<_grid->_gsites;gidx++){

-      // Master RngEngine
-      std::vector<seed_t> master_init(numseed);  src.generate(master_init.begin(),master_init.end());
-      _grid->Broadcast(0,(void *)&master_init[0],sizeof(seed_t)*numseed);
-      fixedSeed master_seed(master_init);
-      RngEngine master_engine(master_seed);
-
-      // Per node RngEngine
-      std::vector<seed_t> node_init(numseed);
-      for(int r=0;r<_grid->ProcessorCount();r++) {
-
-	std::vector<seed_t> rank_init(numseed);
-	for(int i=0;i<numseed;i++) rank_init[i] = uid(master_engine);
-
-	std::cout << GridLogMessage << "SeedSeq for rank "<<r;
-	for(int i=0;i<numseed;i++) std::cout<<" "<<rank_init[i];
-	std::cout <<std::endl;
-
-	if ( r==_grid->ThisRank() ) { 
-	  for(int i=0;i<numseed;i++) node_init[i] = rank_init[i];
-	}
-
-      }
-
-      ////////////////////////////////////////////////////
-      // Set up a seed_seq wrapper with these 8 words
-      // and draw for each site within node.
-      ////////////////////////////////////////////////////
-      fixedSeed node_seed(node_init);
-      RngEngine node_engine(node_seed);
-
-      for(int gidx=0;gidx<gsites;gidx++){
-	int rank,o_idx,i_idx;
+	Skip(master_engine); // Skip to next RNG sequence

+	// Where is it?
 	_grid->GlobalIndexToGlobalCoor(gidx,gcoor);
 	_grid->GlobalCoorToRankIndex(rank,o_idx,i_idx,gcoor);

+	// If this is one of mine we take it
 	if( rank == _grid->ThisRank() ){
 	  int l_idx=generator_idx(o_idx,i_idx);
-	  for(int i=0;i<numseed;i++)  site_init[i] = uid(node_engine);
-	  fixedSeed site_seed(site_init);
-	  _generators[l_idx] = RngEngine(site_seed);
+	  _generators[l_idx] = master_engine;
+	}
+
+      }
+#else 
+      ////////////////////////////////////////////////////////////////
+      // Machine and thread decomposition dependent seeding is efficient
+      // and maximally parallel; but NOT reproducible from machine to machine. 
+      // Not ideal, but fastest way to reseed all nodes.
+      ////////////////////////////////////////////////////////////////
+      {
+	// Obtain one Reseed per processor
+	int Nproc = _grid->ProcessorCount();
+	std::vector<RngEngine> seeders(Nproc);
+	int me= _grid->ThisRank();
+	for(int p=0;p<Nproc;p++){
+	  seeders[p] = Reseed(master_engine);
+	}
+	master_engine = seeders[me];
+      }
+
+      {
+	// Obtain one reseeded generator per thread
+	int Nthread = GridThread::GetThreads();
+	std::vector<RngEngine> seeders(Nthread);
+	for(int t=0;t<Nthread;t++){
+	  seeders[t] = Reseed(master_engine);
+	}
+
+	parallel_for(int t=0;t<Nthread;t++) {
+	  // set up one per local site in threaded fashion
+	  std::vector<uint32_t> newseeds;
+	  std::uniform_int_distribution<uint32_t> uid;	
+	  for(int l=0;l<_grid->lSites();l++) {
+	    if ( (l%Nthread)==t ) {
+	      _generators[l] = Reseed(seeders[t],newseeds,uid);
+	    }
+	  }
 	}
      }
-      _seeded=1;
-    }    
-    void SeedRandomDevice(void){
-      deviceSeed src;
-      Seed(src);
+#endif
    }
-    void SeedFixedIntegers(const std::vector<int> &seeds){
-      CartesianCommunicator::BroadcastWorld(0,(void *)&seeds[0],sizeof(int)*seeds.size());
-      fixedSeed src(seeds);
-      Seed(src);
+    ////////////////////////////////////////////////////////////////////////
+    // Support for rigorous test of RNG's
+    // Return uniform random uint32_t from requested site generator
+    ////////////////////////////////////////////////////////////////////////
+    uint32_t GlobalU01(int gsite){
+
+      uint32_t the_number;
+
+      // who
+      std::vector<int> gcoor;
+      int rank,o_idx,i_idx;
+      _grid->GlobalIndexToGlobalCoor(gsite,gcoor);
+      _grid->GlobalCoorToRankIndex(rank,o_idx,i_idx,gcoor);
+
+      // draw
+      int l_idx=generator_idx(o_idx,i_idx);
+      if( rank == _grid->ThisRank() ){
+	the_number = _uid[l_idx](_generators[l_idx]);
+      }
+      
+      // share & return
+      _grid->Broadcast(rank,(void *)&the_number,sizeof(the_number));
+      return the_number;
    }

  };

-  template <class vobj> inline void random(GridParallelRNG &rng,Lattice<vobj> &l){
-    rng.fill(l,rng._uniform);
-  }
+  template <class vobj> inline void random(GridParallelRNG &rng,Lattice<vobj> &l)   { rng.fill(l,rng._uniform);  }
+  template <class vobj> inline void gaussian(GridParallelRNG &rng,Lattice<vobj> &l) { rng.fill(l,rng._gaussian); }
+  template <class vobj> inline void bernoulli(GridParallelRNG &rng,Lattice<vobj> &l){ rng.fill(l,rng._bernoulli);}

-  template <class vobj> inline void gaussian(GridParallelRNG &rng,Lattice<vobj> &l){
-    rng.fill(l,rng._gaussian);
-  }
-  
-  template <class vobj> inline void bernoulli(GridParallelRNG &rng,Lattice<vobj> &l){
-    rng.fill(l,rng._bernoulli);
-  }
-
-  template <class sobj> inline void random(GridSerialRNG &rng,sobj &l){
-    rng.fill(l,rng._uniform);
-  }
-  
-  template <class sobj> inline void gaussian(GridSerialRNG &rng,sobj &l){
-    rng.fill(l,rng._gaussian);
-  }
-  
-  template <class sobj> inline void bernoulli(GridSerialRNG &rng,sobj &l){
-    rng.fill(l,rng._bernoulli);
-  }
+  template <class sobj> inline void random(GridSerialRNG &rng,sobj &l)   { rng.fill(l,rng._uniform  ); }
+  template <class sobj> inline void gaussian(GridSerialRNG &rng,sobj &l) { rng.fill(l,rng._gaussian ); }
+  template <class sobj> inline void bernoulli(GridSerialRNG &rng,sobj &l){ rng.fill(l,rng._bernoulli); }

 }
 #endif
--- a/lib/parallelIO/NerscIO.h
+++ b/lib/parallelIO/NerscIO.h
@ -491,10 +491,15 @@ static inline void writeRNGState(GridSerialRNG &serial,GridParallelRNG &parallel
 #ifdef RNG_RANLUX
    header.floating_point = std::string("UINT64");
    header.data_type      = std::string("RANLUX48");
-#else
+#endif
+#ifdef RNG_MT19937
    header.floating_point = std::string("UINT32");
    header.data_type      = std::string("MT19937");
 #endif
+#ifdef RNG_SITMO
+    header.floating_point = std::string("UINT64");
+    header.data_type      = std::string("SITMO");
+#endif

  truncate(file);
  offset = writeHeader(header,file);
@ -522,10 +527,15 @@ static inline void readRNGState(GridSerialRNG &serial,GridParallelRNG & parallel
 #ifdef RNG_RANLUX
  assert(format == std::string("UINT64"));
  assert(data_type == std::string("RANLUX48"));
-#else
+#endif
+#ifdef RNG_MT19937
  assert(format == std::string("UINT32"));
  assert(data_type == std::string("MT19937"));
 #endif
+#ifdef RNG_SITMO
+  assert(format == std::string("UINT64"));
+  assert(data_type == std::string("SITMO"));
+#endif

  // depending on datatype, set up munger;
  // munger is a function of <floating point, Real, data_type>
--- a/lib/qcd/action/fermion/CayleyFermion5D.cc
+++ b/lib/qcd/action/fermion/CayleyFermion5D.cc
@ -170,7 +170,6 @@ void CayleyFermion5D<Impl>::Mooee       (const FermionField &psi, FermionField &
  lower[0]   =-mass*lower[0];
  M5D(psi,psi,chi,lower,diag,upper);
 }
-
 template<class Impl>
 void CayleyFermion5D<Impl>::MooeeDag    (const FermionField &psi, FermionField &chi)
 {
@ -192,7 +191,7 @@ void CayleyFermion5D<Impl>::MooeeDag    (const FermionField &psi, FermionField &
      lower[s]=-cee[s-1];
    }
  }
-  // Conjugate the terms ?
+  // Conjugate the terms 
  for (int s=0;s<Ls;s++){
    diag[s] =conjugate(diag[s]);
    upper[s]=conjugate(upper[s]);
@ -219,14 +218,22 @@ void CayleyFermion5D<Impl>::MeooeDag5D    (const FermionField &psi, FermionField
  int Ls=this->Ls;
  std::vector<Coeff_t> diag =bs;
  std::vector<Coeff_t> upper=cs;
-  std::vector<Coeff_t> lower=cs;
-  upper[Ls-1]=-mass*upper[Ls-1];
-  lower[0]   =-mass*lower[0];
-  // Conjugate the terms ?
+  std::vector<Coeff_t> lower=cs; 
+
  for (int s=0;s<Ls;s++){
-    diag[s] =conjugate(diag[s]);
-    upper[s]=conjugate(upper[s]);
-    lower[s]=conjugate(lower[s]);
+    if ( s== 0 ) {
+      upper[s] = cs[s+1];
+      lower[s] =-mass*cs[Ls-1];
+    } else if ( s==(Ls-1) ) { 
+      upper[s] =-mass*cs[0];
+      lower[s] = cs[s-1];
+    } else { 
+      upper[s] = cs[s+1];
+      lower[s] = cs[s-1];
+    }
+    upper[s] = conjugate(upper[s]);
+    lower[s] = conjugate(lower[s]);
+    diag[s]  = conjugate(diag[s]);
  }
  M5Ddag(psi,psi,Din,lower,diag,upper);
 }
@ -313,7 +320,7 @@ void CayleyFermion5D<Impl>::MDeriv  (GaugeField &mat,const FermionField &U,const
    this->DhopDeriv(mat,U,Din,dag);
  } else {
    //      U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
-    Meooe5D(U,Din);
+    MeooeDag5D(U,Din);
    this->DhopDeriv(mat,Din,V,dag);
  }
 };
@ -328,7 +335,7 @@ void CayleyFermion5D<Impl>::MoeDeriv(GaugeField &mat,const FermionField &U,const
    this->DhopDerivOE(mat,U,Din,dag);
  } else {
    //      U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
-      Meooe5D(U,Din);
+      MeooeDag5D(U,Din);
      this->DhopDerivOE(mat,Din,V,dag);
  }
 };
@ -343,7 +350,7 @@ void CayleyFermion5D<Impl>::MeoDeriv(GaugeField &mat,const FermionField &U,const
    this->DhopDerivEO(mat,U,Din,dag);
  } else {
    //      U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
-    Meooe5D(U,Din);
+    MeooeDag5D(U,Din);
    this->DhopDerivEO(mat,Din,V,dag);
  }
 };
--- a/lib/qcd/action/fermion/CayleyFermion5D.h
+++ b/lib/qcd/action/fermion/CayleyFermion5D.h
@ -194,7 +194,9 @@ template void CayleyFermion5D< A >::M5Ddag(const FermionField &psi,const Fermion
 template void CayleyFermion5D< A >::MooeeInv    (const FermionField &psi, FermionField &chi); \
 template void CayleyFermion5D< A >::MooeeInvDag (const FermionField &psi, FermionField &chi);

-#define CAYLEY_DPERP_CACHE
+#undef  CAYLEY_DPERP_DENSE
+#define  CAYLEY_DPERP_CACHE
 #undef  CAYLEY_DPERP_LINALG
+#define CAYLEY_DPERP_VEC

 #endif
--- a/lib/qcd/action/fermion/CayleyFermion5Dcache.cc
+++ b/lib/qcd/action/fermion/CayleyFermion5Dcache.cc
@ -181,6 +181,18 @@ void CayleyFermion5D<Impl>::MooeeInvDag (const FermionField &psi, FermionField &
  assert(psi.checkerboard == psi.checkerboard);
  chi.checkerboard=psi.checkerboard;

+  std::vector<Coeff_t> ueec(Ls);
+  std::vector<Coeff_t> deec(Ls);
+  std::vector<Coeff_t> leec(Ls);
+  std::vector<Coeff_t> ueemc(Ls);
+  std::vector<Coeff_t> leemc(Ls);
+  for(int s=0;s<ueec.size();s++){
+    ueec[s] = conjugate(uee[s]);
+    deec[s] = conjugate(dee[s]);
+    leec[s] = conjugate(lee[s]);
+    ueemc[s]= conjugate(ueem[s]);
+    leemc[s]= conjugate(leem[s]);
+  }
  MooeeInvCalls++;
  MooeeInvTime-=usecond();

@ -192,25 +204,25 @@ void CayleyFermion5D<Impl>::MooeeInvDag (const FermionField &psi, FermionField &
    chi[ss]=psi[ss];
    for (int s=1;s<Ls;s++){
                            spProj5m(tmp,chi[ss+s-1]);
-      chi[ss+s] = psi[ss+s]-uee[s-1]*tmp;
+      chi[ss+s] = psi[ss+s]-ueec[s-1]*tmp;
    }
    // U_m^{-\dagger} 
    for (int s=0;s<Ls-1;s++){
                                   spProj5p(tmp,chi[ss+s]);
-      chi[ss+Ls-1] = chi[ss+Ls-1] - ueem[s]*tmp;
+      chi[ss+Ls-1] = chi[ss+Ls-1] - ueemc[s]*tmp;
    }

    // L_m^{-\dagger} D^{-dagger}
    for (int s=0;s<Ls-1;s++){
      spProj5m(tmp,chi[ss+Ls-1]);
-      chi[ss+s] = (1.0/dee[s])*chi[ss+s]-(leem[s]/dee[Ls-1])*tmp;
+      chi[ss+s] = (1.0/deec[s])*chi[ss+s]-(leemc[s]/deec[Ls-1])*tmp;
    }	
-    chi[ss+Ls-1]= (1.0/dee[Ls-1])*chi[ss+Ls-1];
+    chi[ss+Ls-1]= (1.0/deec[Ls-1])*chi[ss+Ls-1];
  
    // Apply L^{-dagger}
    for (int s=Ls-2;s>=0;s--){
      spProj5p(tmp,chi[ss+s+1]);
-      chi[ss+s] = chi[ss+s] - lee[s]*tmp;
+      chi[ss+s] = chi[ss+s] - leec[s]*tmp;
    }
  }

--- a/lib/qcd/action/fermion/CayleyFermion5Ddense.cc
+++ b/lib/qcd/action/fermion/CayleyFermion5Ddense.cc
@ -39,20 +39,17 @@ namespace QCD {
  /*
   * Dense matrix versions of routines
   */
-
-  /*
 template<class Impl>
 void CayleyFermion5D<Impl>::MooeeInvDag (const FermionField &psi, FermionField &chi)
 {
  this->MooeeInternal(psi,chi,DaggerYes,InverseYes);
 }
-  
 template<class Impl>
 void CayleyFermion5D<Impl>::MooeeInv(const FermionField &psi, FermionField &chi)
 {
  this->MooeeInternal(psi,chi,DaggerNo,InverseYes);
 }
-  */
+
 template<class Impl>
 void CayleyFermion5D<Impl>::MooeeInternal(const FermionField &psi, FermionField &chi,int dag, int inv)
 {
@ -126,9 +123,20 @@ void CayleyFermion5D<Impl>::MooeeInternal(const FermionField &psi, FermionField
  }
 }

+#ifdef CAYLEY_DPERP_DENSE
+INSTANTIATE_DPERP(GparityWilsonImplF);
+INSTANTIATE_DPERP(GparityWilsonImplD);
+INSTANTIATE_DPERP(WilsonImplF);
+INSTANTIATE_DPERP(WilsonImplD);
+INSTANTIATE_DPERP(ZWilsonImplF);
+INSTANTIATE_DPERP(ZWilsonImplD);
+
 template void CayleyFermion5D<GparityWilsonImplF>::MooeeInternal(const FermionField &psi, FermionField &chi,int dag, int inv);
 template void CayleyFermion5D<GparityWilsonImplD>::MooeeInternal(const FermionField &psi, FermionField &chi,int dag, int inv);
 template void CayleyFermion5D<WilsonImplF>::MooeeInternal(const FermionField &psi, FermionField &chi,int dag, int inv);
 template void CayleyFermion5D<WilsonImplD>::MooeeInternal(const FermionField &psi, FermionField &chi,int dag, int inv);
+template void CayleyFermion5D<ZWilsonImplF>::MooeeInternal(const FermionField &psi, FermionField &chi,int dag, int inv);
+template void CayleyFermion5D<ZWilsonImplD>::MooeeInternal(const FermionField &psi, FermionField &chi,int dag, int inv);
+#endif

 }}
--- a/lib/qcd/action/fermion/CayleyFermion5Dssp.cc
+++ b/lib/qcd/action/fermion/CayleyFermion5Dssp.cc
@ -48,17 +48,18 @@ void CayleyFermion5D<Impl>::M5D(const FermionField &psi,
 				std::vector<Coeff_t> &diag,
 				std::vector<Coeff_t> &upper)
 {
+  Coeff_t one(1.0);
  int Ls=this->Ls;
  for(int s=0;s<Ls;s++){
    if ( s==0 ) {
      axpby_ssp_pminus(chi,diag[s],phi,upper[s],psi,s,s+1);
-      axpby_ssp_pplus (chi,1.0,chi,lower[s],psi,s,Ls-1);
+      axpby_ssp_pplus (chi,one,chi,lower[s],psi,s,Ls-1);
    } else if ( s==(Ls-1)) { 
      axpby_ssp_pminus(chi,diag[s],phi,upper[s],psi,s,0);
-      axpby_ssp_pplus (chi,1.0,chi,lower[s],psi,s,s-1);
+      axpby_ssp_pplus (chi,one,chi,lower[s],psi,s,s-1);
    } else {
      axpby_ssp_pminus(chi,diag[s],phi,upper[s],psi,s,s+1);
-      axpby_ssp_pplus(chi,1.0,chi,lower[s],psi,s,s-1);
+      axpby_ssp_pplus(chi,one,chi,lower[s],psi,s,s-1);
    }
  }
 }
@ -70,17 +71,18 @@ void CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi,
 				   std::vector<Coeff_t> &diag,
 				   std::vector<Coeff_t> &upper)
 {
+  Coeff_t one(1.0);
  int Ls=this->Ls;
  for(int s=0;s<Ls;s++){
    if ( s==0 ) {
      axpby_ssp_pplus (chi,diag[s],phi,upper[s],psi,s,s+1);
-      axpby_ssp_pminus(chi,1.0,chi,lower[s],psi,s,Ls-1);
+      axpby_ssp_pminus(chi,one,chi,lower[s],psi,s,Ls-1);
    } else if ( s==(Ls-1)) { 
      axpby_ssp_pplus (chi,diag[s],phi,upper[s],psi,s,0);
-      axpby_ssp_pminus(chi,1.0,chi,lower[s],psi,s,s-1);
+      axpby_ssp_pminus(chi,one,chi,lower[s],psi,s,s-1);
    } else {
      axpby_ssp_pplus (chi,diag[s],phi,upper[s],psi,s,s+1);
-      axpby_ssp_pminus(chi,1.0,chi,lower[s],psi,s,s-1);
+      axpby_ssp_pminus(chi,one,chi,lower[s],psi,s,s-1);
    }
  }
 }
@ -88,62 +90,68 @@ void CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi,
 template<class Impl>
 void CayleyFermion5D<Impl>::MooeeInv    (const FermionField &psi, FermionField &chi)
 {
+  Coeff_t one(1.0);
+  Coeff_t czero(0.0);
  chi.checkerboard=psi.checkerboard;
  int Ls=this->Ls;
  // Apply (L^{\prime})^{-1}
-  axpby_ssp (chi,1.0,psi,     0.0,psi,0,0);      // chi[0]=psi[0]
+  axpby_ssp (chi,one,psi,     czero,psi,0,0);      // chi[0]=psi[0]
  for (int s=1;s<Ls;s++){
-    axpby_ssp_pplus(chi,1.0,psi,-lee[s-1],chi,s,s-1);// recursion Psi[s] -lee P_+ chi[s-1]
+    axpby_ssp_pplus(chi,one,psi,-lee[s-1],chi,s,s-1);// recursion Psi[s] -lee P_+ chi[s-1]
  }
  // L_m^{-1} 
  for (int s=0;s<Ls-1;s++){ // Chi[ee] = 1 - sum[s<Ls-1] -leem[s]P_- chi
-    axpby_ssp_pminus(chi,1.0,chi,-leem[s],chi,Ls-1,s);
+    axpby_ssp_pminus(chi,one,chi,-leem[s],chi,Ls-1,s);
  }
  // U_m^{-1} D^{-1}
  for (int s=0;s<Ls-1;s++){
    // Chi[s] + 1/d chi[s] 
-    axpby_ssp_pplus(chi,1.0/dee[s],chi,-ueem[s]/dee[Ls-1],chi,s,Ls-1);
+    axpby_ssp_pplus(chi,one/dee[s],chi,-ueem[s]/dee[Ls-1],chi,s,Ls-1);
  }	
-  axpby_ssp(chi,1.0/dee[Ls-1],chi,0.0,chi,Ls-1,Ls-1); // Modest avoidable 
+  axpby_ssp(chi,one/dee[Ls-1],chi,czero,chi,Ls-1,Ls-1); // Modest avoidable 
  
  // Apply U^{-1}
  for (int s=Ls-2;s>=0;s--){
-    axpby_ssp_pminus (chi,1.0,chi,-uee[s],chi,s,s+1);  // chi[Ls]
+    axpby_ssp_pminus (chi,one,chi,-uee[s],chi,s,s+1);  // chi[Ls]
  }
 }

 template<class Impl>
 void CayleyFermion5D<Impl>::MooeeInvDag (const FermionField &psi, FermionField &chi)
 {
+  Coeff_t one(1.0);
+  Coeff_t czero(0.0);
  chi.checkerboard=psi.checkerboard;
  int Ls=this->Ls;
  // Apply (U^{\prime})^{-dagger}
-  axpby_ssp (chi,1.0,psi,     0.0,psi,0,0);      // chi[0]=psi[0]
+  axpby_ssp (chi,one,psi,     czero,psi,0,0);      // chi[0]=psi[0]
  for (int s=1;s<Ls;s++){
-    axpby_ssp_pminus(chi,1.0,psi,-uee[s-1],chi,s,s-1);
+    axpby_ssp_pminus(chi,one,psi,-conjugate(uee[s-1]),chi,s,s-1);
  }
  // U_m^{-\dagger} 
  for (int s=0;s<Ls-1;s++){
-    axpby_ssp_pplus(chi,1.0,chi,-ueem[s],chi,Ls-1,s);
+    axpby_ssp_pplus(chi,one,chi,-conjugate(ueem[s]),chi,Ls-1,s);
  }
  // L_m^{-\dagger} D^{-dagger}
  for (int s=0;s<Ls-1;s++){
-    axpby_ssp_pminus(chi,1.0/dee[s],chi,-leem[s]/dee[Ls-1],chi,s,Ls-1);
+    axpby_ssp_pminus(chi,one/conjugate(dee[s]),chi,-conjugate(leem[s]/dee[Ls-1]),chi,s,Ls-1);
  }	
-  axpby_ssp(chi,1.0/dee[Ls-1],chi,0.0,chi,Ls-1,Ls-1); // Modest avoidable 
+  axpby_ssp(chi,one/conjugate(dee[Ls-1]),chi,czero,chi,Ls-1,Ls-1); // Modest avoidable 
  
  // Apply L^{-dagger}
  for (int s=Ls-2;s>=0;s--){
-    axpby_ssp_pplus (chi,1.0,chi,-lee[s],chi,s,s+1);  // chi[Ls]
+    axpby_ssp_pplus (chi,one,chi,-conjugate(lee[s]),chi,s,s+1);  // chi[Ls]
  }
 }


 #ifdef CAYLEY_DPERP_LINALG
-  INSTANTIATE(WilsonImplF);
-  INSTANTIATE(WilsonImplD);
-  INSTANTIATE(GparityWilsonImplF);
-  INSTANTIATE(GparityWilsonImplD);
+  INSTANTIATE_DPERP(WilsonImplF);
+  INSTANTIATE_DPERP(WilsonImplD);
+  INSTANTIATE_DPERP(GparityWilsonImplF);
+  INSTANTIATE_DPERP(GparityWilsonImplD);
+  INSTANTIATE_DPERP(ZWilsonImplF);
+  INSTANTIATE_DPERP(ZWilsonImplD);
 #endif

 }
--- a/lib/qcd/action/fermion/CayleyFermion5Dvec.cc
+++ b/lib/qcd/action/fermion/CayleyFermion5Dvec.cc
@ -35,7 +35,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>


 namespace Grid {
-namespace QCD {  /*
+namespace QCD {  
+  /*
   * Dense matrix versions of routines
   */
 template<class Impl>
--- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc
+++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc
@ -40,10 +40,10 @@ ImprovedStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3,
 // Constructor and gauge import
 /////////////////////////////////

+
 template <class Impl>
-ImprovedStaggeredFermion<Impl>::ImprovedStaggeredFermion(GaugeField &_Uthin, GaugeField &_Ufat, GridCartesian &Fgrid,
-							 GridRedBlackCartesian &Hgrid, RealD _mass,
-							 RealD _c1, RealD _c2,RealD _u0,
+ImprovedStaggeredFermion<Impl>::ImprovedStaggeredFermion(GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, 
+							 RealD _mass,
 							 const ImplParams &p)
    : Kernels(p),
      _grid(&Fgrid),
@ -52,9 +52,6 @@ ImprovedStaggeredFermion<Impl>::ImprovedStaggeredFermion(GaugeField &_Uthin, Gau
      StencilEven(&Hgrid, npoint, Even, directions, displacements),  // source is Even
      StencilOdd(&Hgrid, npoint, Odd, directions, displacements),  // source is Odd
      mass(_mass),
-      c1(_c1),
-      c2(_c2),
-      u0(_u0),
      Lebesgue(_grid),
      LebesgueEvenOdd(_cbgrid),
      Umu(&Fgrid),
@ -65,9 +62,29 @@ ImprovedStaggeredFermion<Impl>::ImprovedStaggeredFermion(GaugeField &_Uthin, Gau
      UUUmuOdd(&Hgrid) ,
      _tmp(&Hgrid)
 {
-  // Allocate the required comms buffer
+}
+
+template <class Impl>
+ImprovedStaggeredFermion<Impl>::ImprovedStaggeredFermion(GaugeField &_Uthin, GaugeField &_Ufat, GridCartesian &Fgrid,
+							 GridRedBlackCartesian &Hgrid, RealD _mass,
+							 RealD _c1, RealD _c2,RealD _u0,
+							 const ImplParams &p)
+  : ImprovedStaggeredFermion(Fgrid,Hgrid,_mass,p)
+{
+  c1=_c1;
+  c2=_c2;
+  u0=_u0;
  ImportGauge(_Uthin,_Ufat);
 }
+template <class Impl>
+ImprovedStaggeredFermion<Impl>::ImprovedStaggeredFermion(GaugeField &_Uthin,GaugeField &_Utriple, GaugeField &_Ufat, GridCartesian &Fgrid,
+							 GridRedBlackCartesian &Hgrid, RealD _mass,
+							 const ImplParams &p)
+  : ImprovedStaggeredFermion(Fgrid,Hgrid,_mass,p)
+{
+  ImportGaugeSimple(_Utriple,_Ufat);
+}
+

  ////////////////////////////////////////////////////////////
  // Momentum space propagator should be 
@ -86,6 +103,34 @@ void ImprovedStaggeredFermion<Impl>::ImportGauge(const GaugeField &_Uthin)
  ImportGauge(_Uthin,_Uthin);
 };
 template <class Impl>
+void ImprovedStaggeredFermion<Impl>::ImportGaugeSimple(const GaugeField &_Utriple,const GaugeField &_Ufat) 
+{
+  /////////////////////////////////////////////////////////////////
+  // Trivial import; phases and fattening and such like preapplied
+  /////////////////////////////////////////////////////////////////
+  GaugeLinkField U(GaugeGrid());
+
+  for (int mu = 0; mu < Nd; mu++) {
+
+    U = PeekIndex<LorentzIndex>(_Utriple, mu);
+    PokeIndex<LorentzIndex>(UUUmu, U, mu );
+
+    U = adj( Cshift(U, mu, -3));
+    PokeIndex<LorentzIndex>(UUUmu, -U, mu+4 );
+
+    U = PeekIndex<LorentzIndex>(_Ufat, mu);
+    PokeIndex<LorentzIndex>(Umu, U, mu);
+
+    U = adj( Cshift(U, mu, -1));
+    PokeIndex<LorentzIndex>(Umu, -U, mu+4);
+
+  }
+  pickCheckerboard(Even, UmuEven,  Umu);
+  pickCheckerboard(Odd,  UmuOdd ,  Umu);
+  pickCheckerboard(Even, UUUmuEven,UUUmu);
+  pickCheckerboard(Odd,  UUUmuOdd, UUUmu);
+}
+template <class Impl>
 void ImprovedStaggeredFermion<Impl>::ImportGauge(const GaugeField &_Uthin,const GaugeField &_Ufat) 
 {
  GaugeLinkField U(GaugeGrid());
@ -115,6 +160,8 @@ void ImprovedStaggeredFermion<Impl>::ImportGauge(const GaugeField &_Uthin,const
    PokeIndex<LorentzIndex>(UUUmu, U*(-0.5*c2/u0/u0/u0), mu+4);
  }

+  std::cout << " Umu " << Umu._odata[0]<<std::endl;
+  std::cout << " UUUmu " << UUUmu._odata[0]<<std::endl;
  pickCheckerboard(Even, UmuEven, Umu);
  pickCheckerboard(Odd,  UmuOdd , Umu);
  pickCheckerboard(Even, UUUmuEven, UUUmu);
--- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h
+++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h
@ -112,7 +112,16 @@ class ImprovedStaggeredFermion : public StaggeredKernels<Impl>, public ImprovedS
 			   RealD _c1=9.0/8.0, RealD _c2=-1.0/24.0,RealD _u0=1.0,
 			   const ImplParams &p = ImplParams());

+  ImprovedStaggeredFermion(GaugeField &_Uthin, GaugeField &_Utriple, GaugeField &_Ufat, GridCartesian &Fgrid,
+			   GridRedBlackCartesian &Hgrid, RealD _mass,
+			   const ImplParams &p = ImplParams());
+
+  ImprovedStaggeredFermion(GridCartesian &Fgrid, GridRedBlackCartesian &Hgrid, RealD _mass,
+			   const ImplParams &p = ImplParams());
+
+
  // DoubleStore impl dependent
+  void ImportGaugeSimple(const GaugeField &_Utriple, const GaugeField &_Ufat);
  void ImportGauge(const GaugeField &_Uthin, const GaugeField &_Ufat);
  void ImportGauge(const GaugeField &_Uthin);

--- a/lib/qcd/action/fermion/StaggeredKernelsAsm.cc
+++ b/lib/qcd/action/fermion/StaggeredKernelsAsm.cc
@ -27,8 +27,11 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
    *************************************************************************************/
    /*  END LEGAL */
 #include <Grid.h>
+
+#ifdef AVX512
 #include <simd/Intel512common.h>
 #include <simd/Intel512avx.h>
+#endif

 // Interleave operations from two directions
 // This looks just like a 2 spin multiply and reuse same sequence from the Wilson
@ -302,7 +305,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
  VRDUP(Chi_00,T0)           VIDUP(Chi_00,Chi_00)	          \
   VRDUP(Chi_10,T1)           VIDUP(Chi_10,Chi_10)		  \
   VMUL(Z00,Chi_00,Z1)        VMUL(Z10,Chi_10,Z2)		  \
-   VSHUFMEM(3,%r8,Z00)	      VSHUFMEM(3,%r9,Z10)		  \    
+   VSHUFMEM(3,%r8,Z00)	      VSHUFMEM(3,%r9,Z10)		  \
   VMUL(Z00,Chi_00,Z3)        VMUL(Z10,Chi_10,Z4)		  \
   VSHUFMEM(6,%r8,Z00)	      VSHUFMEM(6,%r9,Z10)		  \
   VMUL(Z00,Chi_00,Z5)        VMUL(Z10,Chi_10,Z6)		  \
@ -584,7 +587,6 @@ void StaggeredKernels<Impl>::DhopSiteAsm(StencilImpl &st, LebesgueOrder &lo,
 					 int sU, const FermionField &in, FermionField &out) 
 {
  assert(0);
-
 };


@ -902,9 +904,17 @@ template <> void StaggeredKernels<StaggeredImplD>::DhopSiteAsm(StencilImpl &st,
 #endif
 }

+#define KERNEL_INSTANTIATE(CLASS,FUNC,IMPL)			    \
+  template void CLASS<IMPL>::FUNC(StencilImpl &st, LebesgueOrder &lo,	\
+				  DoubledGaugeField &U,			\
+				  DoubledGaugeField &UUU,		\
+				  SiteSpinor *buf, int LLs,		\
+				  int sU, const FermionField &in, FermionField &out);

-FermOpStaggeredTemplateInstantiate(StaggeredKernels);
-FermOpStaggeredVec5dTemplateInstantiate(StaggeredKernels);
+KERNEL_INSTANTIATE(StaggeredKernels,DhopSiteAsm,StaggeredImplD);
+KERNEL_INSTANTIATE(StaggeredKernels,DhopSiteAsm,StaggeredImplF);
+KERNEL_INSTANTIATE(StaggeredKernels,DhopSiteAsm,StaggeredVec5dImplD);
+KERNEL_INSTANTIATE(StaggeredKernels,DhopSiteAsm,StaggeredVec5dImplF);

 }}

--- a/lib/qcd/action/fermion/StaggeredKernelsHand.cc
+++ b/lib/qcd/action/fermion/StaggeredKernelsHand.cc
@ -299,7 +299,24 @@ void StaggeredKernels<Impl>::DhopSiteDepthHand(StencilImpl &st, LebesgueOrder &l

 }

-FermOpStaggeredTemplateInstantiate(StaggeredKernels);
-FermOpStaggeredVec5dTemplateInstantiate(StaggeredKernels);
+#define DHOP_SITE_HAND_INSTANTIATE(IMPL)				\
+  template void StaggeredKernels<IMPL>::DhopSiteHand(StencilImpl &st, LebesgueOrder &lo, \
+						     DoubledGaugeField &U,DoubledGaugeField &UUU, \
+						     SiteSpinor *buf, int LLs, \
+						     int sU, const FermionField &in, FermionField &out, int dag);
+
+#define DHOP_SITE_DEPTH_HAND_INSTANTIATE(IMPL)				\
+  template void StaggeredKernels<IMPL>::DhopSiteDepthHand(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, \
+							  SiteSpinor *buf, int sF, \
+							  int sU, const FermionField &in, SiteSpinor &out,int threeLink) ;
+DHOP_SITE_HAND_INSTANTIATE(StaggeredImplD);
+DHOP_SITE_HAND_INSTANTIATE(StaggeredImplF);
+DHOP_SITE_HAND_INSTANTIATE(StaggeredVec5dImplD);
+DHOP_SITE_HAND_INSTANTIATE(StaggeredVec5dImplF);
+
+DHOP_SITE_DEPTH_HAND_INSTANTIATE(StaggeredImplD);
+DHOP_SITE_DEPTH_HAND_INSTANTIATE(StaggeredImplF);
+DHOP_SITE_DEPTH_HAND_INSTANTIATE(StaggeredVec5dImplD);
+DHOP_SITE_DEPTH_HAND_INSTANTIATE(StaggeredVec5dImplF);

 }}
--- a/lib/qcd/hmc/HmcRunner.h
+++ b/lib/qcd/hmc/HmcRunner.h
@ -114,8 +114,8 @@ class NerscHmcRunnerTemplate {
    */
    //////////////
    NoSmearing<Gimpl> SmearingPolicy;
-    typedef MinimumNorm2<GaugeField, NoSmearing<Gimpl>, RepresentationsPolicy >
-        IntegratorType;  // change here to change the algorithm
+    // change here to change the algorithm
+    typedef MinimumNorm2<GaugeField, NoSmearing<Gimpl>, RepresentationsPolicy >  IntegratorType;  
    IntegratorParameters MDpar(40, 1.0);
    IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);

--- a/scripts/zmobius.sh
+++ b/scripts/zmobius.sh
@ -0,0 +1,35 @@
+#!/bin/bash
+fn=$1
+
+grep "double zmobius_" $fn |
+awk 'BEGIN{ m["zmobius_b_coeff"]=0; m["zmobius_c_coeff"]=1; }{ val[m[substr($2,0,15)]][substr($2,17)+0]=$4; }END{
+
+    ls=length(val[0])/2;
+
+    print "ls = " ls
+
+    bmc=-111;
+
+    for (s=0;s<ls;s++) {
+      br[s] = val[0][2*s + 0];
+      bi[s] = val[0][2*s + 1];
+      cr[s] = val[1][2*s + 0];
+      ci[s] = val[1][2*s + 1];
+
+      t=br[s] - cr[s];
+      if (bmc == -111)
+        bmc=t;
+      else if (bmc != t)
+        print "Warning: b-c is not constant!";
+
+      omegar[s] = (-1.0 + 2.0* br[s])/(4.0*bi[s]**2.0 + (1.0 - 2.0* br[s])**2);
+      omegai[s] = - 2.0* bi[s]/(4.0*bi[s]**2.0 + (1.0 - 2.0* br[s])**2);
+    }
+
+    print "b-c = " bmc
+
+    for (s=0;s<ls;s++) {
+      printf( "omega.push_back( std::complex<double>(%.15g,%.15g) );\n",omegar[s],omegai[s]);
+    }
+
+}'
--- a/tests/IO/Test_nersc_io.cc
+++ b/tests/IO/Test_nersc_io.cc
@ -54,8 +54,8 @@ int main (int argc, char ** argv)
  GridSerialRNG     sRNGa;
  GridSerialRNG     sRNGb;

-  pRNGa.SeedRandomDevice();
-  sRNGa.SeedRandomDevice();
+  pRNGa.SeedFixedIntegers(std::vector<int>({45,12,81,9});
+  sRNGa.SeedFixedIntegers(std::vector<int>({45,12,81,9});
  
  std::string rfile("./ckpoint_rng.4000");
  NerscIO::writeRNGState(sRNGa,pRNGa,rfile);
--- a/tests/Test_cshift.cc
+++ b/tests/Test_cshift.cc
@ -41,7 +41,7 @@ int main (int argc, char ** argv)

  GridCartesian        Fine(latt_size,simd_layout,mpi_layout);

-  GridParallelRNG      FineRNG(&Fine);  FineRNG.SeedRandomDevice();
+  GridParallelRNG      FineRNG(&Fine);  FineRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  LatticeComplex U(&Fine);
  LatticeComplex ShiftU(&Fine);
--- a/tests/Test_simd.cc
+++ b/tests/Test_simd.cc
@ -125,7 +125,7 @@ template<class scal, class vec,class functor >
 void Tester(const functor &func)
 {
  GridSerialRNG          sRNG;
-  sRNG.SeedRandomDevice();
+  sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
  
  int Nsimd = vec::Nsimd();

@ -184,7 +184,7 @@ void IntTester(const functor &func)
  typedef Integer  scal;
  typedef vInteger vec;
  GridSerialRNG          sRNG;
-  sRNG.SeedRandomDevice();
+  sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  int Nsimd = vec::Nsimd();

@ -242,7 +242,7 @@ template<class reduced,class scal, class vec,class functor >
 void ReductionTester(const functor &func)
 {
  GridSerialRNG          sRNG;
-  sRNG.SeedRandomDevice();
+  sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
  
  int Nsimd = vec::Nsimd();

@ -343,7 +343,7 @@ template<class scal, class vec,class functor >
 void PermTester(const functor &func)
 {
  GridSerialRNG          sRNG;
-  sRNG.SeedRandomDevice();
+  sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
  
  int Nsimd = vec::Nsimd();

@ -409,7 +409,7 @@ template<class scal, class vec,class functor >
 void ExchangeTester(const functor &func)
 {
  GridSerialRNG          sRNG;
-  sRNG.SeedRandomDevice();
+  sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
  
  int Nsimd = vec::Nsimd();

--- a/tests/Test_stencil.cc
+++ b/tests/Test_stencil.cc
@ -52,7 +52,7 @@ int main (int argc, char ** argv)
  GridRedBlackCartesian rbFine(latt_size,simd_layout,mpi_layout);
  GridParallelRNG       fRNG(&Fine);

-  //  fRNG.SeedRandomDevice();
+  //  fRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});
  std::vector<int> seeds({1,2,3,4});
  fRNG.SeedFixedIntegers(seeds);
  
--- a/tests/core/Test_cshift_red_black.cc
+++ b/tests/core/Test_cshift_red_black.cc
@ -49,7 +49,7 @@ int main (int argc, char ** argv)
  GridCartesian         Fine  (latt_size,simd_layout,mpi_layout);
  GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout,mask,1);

-  GridParallelRNG      FineRNG(&Fine);  FineRNG.SeedRandomDevice();
+  GridParallelRNG      FineRNG(&Fine);  FineRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  LatticeComplex U(&Fine);
  LatticeComplex ShiftU(&Fine);
--- a/tests/core/Test_cshift_red_black_rotate.cc
+++ b/tests/core/Test_cshift_red_black_rotate.cc
@ -49,7 +49,7 @@ int main (int argc, char ** argv)
  GridCartesian         Fine  (latt_size,simd_layout,mpi_layout);
  GridRedBlackCartesian RBFine(latt_size,simd_layout,mpi_layout,mask,1);

-  GridParallelRNG      FineRNG(&Fine);  FineRNG.SeedRandomDevice();
+  GridParallelRNG      FineRNG(&Fine);  FineRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  LatticeComplex err(&Fine);
  LatticeComplex U(&Fine);
--- a/tests/core/Test_cshift_rotate.cc
+++ b/tests/core/Test_cshift_rotate.cc
@ -41,7 +41,7 @@ int main (int argc, char ** argv)

  GridCartesian        Fine(latt_size,simd_layout,mpi_layout);

-  GridParallelRNG      FineRNG(&Fine);  FineRNG.SeedRandomDevice();
+  GridParallelRNG      FineRNG(&Fine);  FineRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  LatticeComplex U(&Fine);
  LatticeComplex ShiftU(&Fine);
--- a/tests/core/Test_gamma.cc
+++ b/tests/core/Test_gamma.cc
@ -245,7 +245,7 @@ int main(int argc, char *argv[])
  GridCartesian Grid(latt_size,simd_layout,mpi_layout);
  GridSerialRNG sRNG;
  
-  sRNG.SeedRandomDevice();
+  sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
  
  std::cout << GridLogMessage << "======== Test algebra" << std::endl;
  createTestAlgebra();
--- a/tests/core/Test_gpwilson_even_odd.cc
+++ b/tests/core/Test_gpwilson_even_odd.cc
@ -50,7 +50,7 @@ int main (int argc, char ** argv)
  GridParallelRNG          pRNG(&Grid);
  //  std::vector<int> seeds({1,2,3,4});
  //  pRNG.SeedFixedIntegers(seeds);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  typedef typename GparityWilsonFermionR::FermionField FermionField;

--- a/tests/core/Test_lie_generators.cc
+++ b/tests/core/Test_lie_generators.cc
@ -86,7 +86,7 @@ int main(int argc, char** argv) {

  // Projectors 
  GridParallelRNG gridRNG(grid);
-  gridRNG.SeedRandomDevice();
+  gridRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
  SU3Adjoint::LatticeAdjMatrix Gauss(grid);
  SU3::LatticeAlgebraVector ha(grid);
  SU3::LatticeAlgebraVector hb(grid);
--- a/tests/core/Test_main.cc
+++ b/tests/core/Test_main.cc
@ -89,8 +89,8 @@ int main(int argc, char **argv) {
      GridSerialRNG SerialRNG;
      GridSerialRNG SerialRNG1;

-      FineRNG.SeedRandomDevice();
-      SerialRNG.SeedRandomDevice();
+      FineRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+      SerialRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

      std::cout << "SerialRNG" << SerialRNG._generators[0] << std::endl;

--- a/tests/core/Test_rng.cc
+++ b/tests/core/Test_rng.cc
@ -43,10 +43,10 @@ int main (int argc, char ** argv)

  std::vector<int> seeds({1,2,3,4});

-  GridSerialRNG             sRNG;   sRNG.SeedRandomDevice();
+  GridSerialRNG             sRNG;   sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
  GridSerialRNG            fsRNG;  fsRNG.SeedFixedIntegers(seeds);

-  GridParallelRNG           pRNG(&Grid);   pRNG.SeedRandomDevice();
+  GridParallelRNG           pRNG(&Grid);   pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
  GridParallelRNG          fpRNG(&Grid);  fpRNG.SeedFixedIntegers(seeds);

  SpinMatrix rnd  ; 
--- a/tests/core/Test_staggered.cc
+++ b/tests/core/Test_staggered.cc
@ -51,7 +51,7 @@ int main (int argc, char ** argv)
  std::vector<int> seeds({1,2,3,4});
  GridParallelRNG          pRNG(&Grid);
  pRNG.SeedFixedIntegers(seeds);
-  //  pRNG.SeedRandomDevice();
+  //  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});

  typedef typename ImprovedStaggeredFermionR::FermionField FermionField; 
  typedef typename ImprovedStaggeredFermionR::ComplexField ComplexField; 
--- a/tests/core/Test_wilson_even_odd.cc
+++ b/tests/core/Test_wilson_even_odd.cc
@ -62,7 +62,7 @@ int main (int argc, char ** argv)
  GridParallelRNG          pRNG(&Grid);
  //  std::vector<int> seeds({1,2,3,4});
  //  pRNG.SeedFixedIntegers(seeds);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  LatticeFermion src   (&Grid); random(pRNG,src);
  LatticeFermion phi   (&Grid); random(pRNG,phi);
--- a/tests/core/Test_wilson_tm_even_odd.cc
+++ b/tests/core/Test_wilson_tm_even_odd.cc
@ -61,7 +61,7 @@ int main (int argc, char ** argv)
  GridParallelRNG          pRNG(&Grid);
  //  std::vector<int> seeds({1,2,3,4});
  //  pRNG.SeedFixedIntegers(seeds);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  LatticeFermion src   (&Grid); random(pRNG,src);
  LatticeFermion phi   (&Grid); random(pRNG,phi);
--- a/tests/core/Test_zmobius_even_odd.cc
+++ b/tests/core/Test_zmobius_even_odd.cc
@ -0,0 +1,287 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_dwf_even_odd.cc
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: paboyle <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>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+template<class d>
+struct scal {
+  d internal;
+};
+
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
+  };
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+
+  const int Ls=10;
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  std::vector<int> seeds4({1,2,3,4});
+  std::vector<int> seeds5({5,6,7,8});
+
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+  GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
+
+  LatticeFermion src   (FGrid); random(RNG5,src);
+  LatticeFermion phi   (FGrid); random(RNG5,phi);
+  LatticeFermion chi   (FGrid); random(RNG5,chi);
+  LatticeFermion result(FGrid); result=zero;
+  LatticeFermion    ref(FGrid);    ref=zero;
+  LatticeFermion    tmp(FGrid);    tmp=zero;
+  LatticeFermion    err(FGrid);    tmp=zero;
+  LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
+  std::vector<LatticeColourMatrix> U(4,UGrid);
+
+  // Only one non-zero (y)
+  Umu=zero;
+  for(int nn=0;nn<Nd;nn++){
+    random(RNG4,U[nn]);
+    if ( nn>0 ) 
+      U[nn]=zero;
+    PokeIndex<LorentzIndex>(Umu,U[nn],nn);
+  }
+
+  RealD mass=0.1;
+  RealD M5  =1.8;
+  std::vector < std::complex<double>  > omegas;
+#if 0
+  for(int i=0;i<Ls;i++){
+    double imag = 0.;
+    if (i==0) imag=1.;
+    if (i==Ls-1) imag=-1.;
+    std::complex<double> temp (0.25+0.01*i, imag*0.01);
+    omegas.push_back(temp);
+  }
+#else
+  omegas.push_back( std::complex<double>(1.45806438985048,-0) );
+  omegas.push_back( std::complex<double>(1.18231318389348,-0) );
+  omegas.push_back( std::complex<double>(0.830951166685955,-0) );
+  omegas.push_back( std::complex<double>(0.542352409156791,-0) );
+  omegas.push_back( std::complex<double>(0.341985020453729,-0) );
+  omegas.push_back( std::complex<double>(0.21137902619029,-0) );
+  omegas.push_back( std::complex<double>(0.126074299502912,-0) );
+  omegas.push_back( std::complex<double>(0.0990136651962626,-0) );
+  omegas.push_back( std::complex<double>(0.0686324988446592,0.0550658530827402) );
+  omegas.push_back( std::complex<double>(0.0686324988446592,-0.0550658530827402) );
+#endif
+
+  ZMobiusFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, omegas,1.,0.);
+//  DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+
+  LatticeFermion src_e (FrbGrid);
+  LatticeFermion src_o (FrbGrid);
+  LatticeFermion r_e   (FrbGrid);
+  LatticeFermion r_o   (FrbGrid);
+  LatticeFermion r_eo  (FGrid);
+  LatticeFermion r_eeoo(FGrid);
+
+  std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
+  std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
+
+  pickCheckerboard(Even,src_e,src);
+  pickCheckerboard(Odd,src_o,src);
+
+  Ddwf.Meooe(src_e,r_o);  std::cout<<GridLogMessage<<"Applied Meo"<<std::endl;
+  Ddwf.Meooe(src_o,r_e);  std::cout<<GridLogMessage<<"Applied Moe"<<std::endl;
+  setCheckerboard(r_eo,r_o);
+  setCheckerboard(r_eo,r_e);
+
+  Ddwf.Mooee(src_e,r_e);  std::cout<<GridLogMessage<<"Applied Mee"<<std::endl;
+  Ddwf.Mooee(src_o,r_o);  std::cout<<GridLogMessage<<"Applied Moo"<<std::endl;
+  setCheckerboard(r_eeoo,r_e);
+  setCheckerboard(r_eeoo,r_o);
+
+  r_eo=r_eo+r_eeoo;
+  Ddwf.M(src,ref);  
+
+  //  std::cout<<GridLogMessage << r_eo<<std::endl;
+  //  std::cout<<GridLogMessage << ref <<std::endl;
+
+  err= ref - r_eo;
+  std::cout<<GridLogMessage << "EO norm diff   "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
+    
+  LatticeComplex cerr(FGrid);
+  cerr = localInnerProduct(err,err);
+  //  std::cout<<GridLogMessage << cerr<<std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Test MooeeDagger is the dagger of Mooee by requiring                "<<std::endl;
+  std::cout<<GridLogMessage<<"=  < phi | Deo | chi > * = < chi | Deo^dag| phi>  "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+
+  LatticeFermion chi_e   (FrbGrid);
+  LatticeFermion chi_o   (FrbGrid);
+
+  LatticeFermion dchi_e  (FrbGrid);
+  LatticeFermion dchi_o  (FrbGrid);
+
+  LatticeFermion phi_e   (FrbGrid);
+  LatticeFermion phi_o   (FrbGrid);
+
+  LatticeFermion dphi_e  (FrbGrid);
+  LatticeFermion dphi_o  (FrbGrid);
+
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+  pickCheckerboard(Even,phi_e,phi);
+  pickCheckerboard(Odd ,phi_o,phi);
+
+  Ddwf.Mooee(chi_e,dchi_o);
+  Ddwf.Mooee(chi_o,dchi_e);
+  Ddwf.MooeeDag(phi_e,dphi_o);
+  Ddwf.MooeeDag(phi_o,dphi_e);
+
+  ComplexD pDce = innerProduct(phi_e,dchi_e);
+  ComplexD pDco = innerProduct(phi_o,dchi_o);
+  ComplexD cDpe = innerProduct(chi_e,dphi_e);
+  ComplexD cDpo = innerProduct(chi_o,dphi_o);
+
+
+  std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
+  std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
+
+  std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
+  std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring                "<<std::endl;
+  std::cout<<GridLogMessage<<"=  < phi | Deo | chi > * = < chi | Deo^dag| phi>  "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  
+
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+  pickCheckerboard(Even,phi_e,phi);
+  pickCheckerboard(Odd ,phi_o,phi);
+
+  Ddwf.Meooe(chi_e,dchi_o);
+  Ddwf.Meooe(chi_o,dchi_e);
+  Ddwf.MeooeDag(phi_e,dphi_o);
+  Ddwf.MeooeDag(phi_o,dphi_e);
+
+  pDce = innerProduct(phi_e,dchi_e);
+  pDco = innerProduct(phi_o,dchi_o);
+  cDpe = innerProduct(chi_e,dphi_e);
+  cDpo = innerProduct(chi_o,dphi_o);
+
+  std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
+  std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
+
+  std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
+  std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1                                         "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+
+  Ddwf.Mooee(chi_e,src_e);
+  Ddwf.MooeeInv(src_e,phi_e);
+
+  Ddwf.Mooee(chi_o,src_o);
+  Ddwf.MooeeInv(src_o,phi_o);
+  
+  setCheckerboard(phi,phi_e);
+  setCheckerboard(phi,phi_o);
+
+  err = phi-chi;
+  std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<< std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1                                   "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+
+  Ddwf.MooeeDag(chi_e,src_e);
+  Ddwf.MooeeInvDag(src_e,phi_e);
+
+  Ddwf.MooeeDag(chi_o,src_o);
+  Ddwf.MooeeInvDag(src_o,phi_o);
+  
+  setCheckerboard(phi,phi_e);
+  setCheckerboard(phi,phi_o);
+
+  err = phi-chi;
+  std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<< std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian              "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  
+  random(RNG5,phi);
+  random(RNG5,chi);
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+  pickCheckerboard(Even,phi_e,phi);
+  pickCheckerboard(Odd ,phi_o,phi);
+  RealD t1,t2;
+
+
+  SchurDiagMooeeOperator<ZMobiusFermionR,LatticeFermion> HermOpEO(Ddwf);
+  HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
+  HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
+
+  HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
+  HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
+
+  pDce = innerProduct(phi_e,dchi_e);
+  pDco = innerProduct(phi_o,dchi_o);
+  cDpe = innerProduct(chi_e,dphi_e);
+  cDpo = innerProduct(chi_o,dphi_o);
+
+  std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
+  std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
+
+  std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
+  std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
+  
+  Grid_finalize();
+}
--- a/tests/forces/Test_dwf_gpforce.cc
+++ b/tests/forces/Test_dwf_gpforce.cc
@ -54,8 +54,8 @@ int main (int argc, char ** argv)

  std::vector<int> seeds({1,2,3,4});

-  GridParallelRNG          RNG5(FGrid);  RNG5.SeedRandomDevice();
-  GridParallelRNG          RNG4(UGrid);  RNG4.SeedRandomDevice();
+  GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
  
  FermionField phi        (FGrid); gaussian(RNG5,phi);
  FermionField Mphi       (FGrid); 
--- a/tests/forces/Test_gp_rect_force.cc
+++ b/tests/forces/Test_gp_rect_force.cc
@ -50,7 +50,7 @@ int main (int argc, char ** argv)
  std::vector<int> seeds({1,2,3,4});

  GridParallelRNG          pRNG(&Grid);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  LatticeGaugeField U(&Grid);

--- a/tests/forces/Test_rect_force.cc
+++ b/tests/forces/Test_rect_force.cc
@ -50,7 +50,7 @@ int main (int argc, char ** argv)
  std::vector<int> seeds({1,2,3,4});

  GridParallelRNG          pRNG(&Grid);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  LatticeGaugeField U(&Grid);

--- a/tests/forces/Test_wilson_force.cc
+++ b/tests/forces/Test_wilson_force.cc
@ -50,7 +50,7 @@ int main (int argc, char ** argv)
  std::vector<int> seeds({1,2,3,4});

  GridParallelRNG          pRNG(&Grid);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  LatticeFermion phi        (&Grid); gaussian(pRNG,phi);
  LatticeFermion Mphi       (&Grid); 
--- a/tests/forces/Test_wilson_force_phiMdagMphi.cc
+++ b/tests/forces/Test_wilson_force_phiMdagMphi.cc
@ -50,7 +50,7 @@ int main (int argc, char ** argv)
  std::vector<int> seeds({1,2,3,4});

  GridParallelRNG          pRNG(&Grid);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  LatticeFermion phi        (&Grid); gaussian(pRNG,phi);
  LatticeFermion Mphi       (&Grid); 
--- a/tests/forces/Test_wilson_force_phiMphi.cc
+++ b/tests/forces/Test_wilson_force_phiMphi.cc
@ -50,7 +50,7 @@ int main (int argc, char ** argv)
  std::vector<int> seeds({1,2,3,4});

  GridParallelRNG          pRNG(&Grid);
-  pRNG.SeedRandomDevice();
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));

  LatticeFermion phi        (&Grid); gaussian(pRNG,phi);
  LatticeFermion Mphi       (&Grid); 
--- a/tests/qdpxx/Test_qdpxx_loops_staples.cc
+++ b/tests/qdpxx/Test_qdpxx_loops_staples.cc
@ -282,8 +282,8 @@ double calc_grid_p(Grid::QCD::LatticeGaugeField & Umu)
  Grid::QCD::LatticeColourMatrix tmp(UGrid); 
  tmp = Grid::zero;

-  Grid::QCD::PokeIndex<Grid::QCD::LorentzIndex>(Umu,tmp,2);
-  Grid::QCD::PokeIndex<Grid::QCD::LorentzIndex>(Umu,tmp,3);
+  Grid::QCD::PokeIndex<LorentzIndex>(Umu,tmp,2);
+  Grid::QCD::PokeIndex<LorentzIndex>(Umu,tmp,3);

  Grid::QCD::WilsonGaugeActionR Wilson(beta); // Just take beta = 1.0
  
@ -311,7 +311,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)

  std::vector<Grid::QCD::LatticeColourMatrix> U(4,UGrid);
  for(int mu=0;mu<Nd;mu++){
-    U[mu] = Grid::PeekIndex<Grid::QCD::LorentzIndex>(Umu,mu);
+    U[mu] = Grid::PeekIndex<LorentzIndex>(Umu,mu);
  }

  Grid::QCD::LatticeComplex rect(UGrid);
@ -322,7 +322,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
  for(int nu=0;nu<Grid::QCD::Nd;nu++){
    if ( mu!=nu ) {

-      Grid::QCD::WilsonLoops<Grid::QCD::LatticeGaugeField>::traceDirRectangle(rect,U,mu,nu);
+      Grid::QCD::ColourWilsonLoops::traceDirRectangle(rect,U,mu,nu);
      trect = Grid::sum(rect);
      crect = Grid::TensorRemove(trect);
      std::cout<< "mu/nu = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/2.0/3.0<<std::endl;
@ -344,10 +344,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	//           __ ___ 
 	//          |    __ |
 	Stap = 
-	  Grid::Cshift(Grid::QCD::CovShiftForward (U[mu],mu,
-		       Grid::QCD::CovShiftForward (U[nu],nu,
-		       Grid::QCD::CovShiftBackward(U[mu],mu,
-                       Grid::QCD::CovShiftBackward(U[mu],mu,
+	  Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[mu],mu,
+		       Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu,
+		       Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
+                       Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
 		       Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);

 	TrStap = Grid::trace (U[mu]*Stap);
@ -361,10 +361,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	//              __ 
 	//          |__ __ |

-	Stap = Grid::Cshift(Grid::QCD::CovShiftForward (U[mu],mu,
-		            Grid::QCD::CovShiftBackward(U[nu],nu,
-   		            Grid::QCD::CovShiftBackward(U[mu],mu,
-                            Grid::QCD::CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1);
+	Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[mu],mu,
+		            Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
+   		            Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
+                            Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1);

 	TrStap = Grid::trace (U[mu]*Stap);

@ -375,10 +375,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	//           __ 
 	//          |__ __ |

-	Stap = Grid::Cshift(Grid::QCD::CovShiftBackward(U[nu],nu,
-		            Grid::QCD::CovShiftBackward(U[mu],mu,
-                            Grid::QCD::CovShiftBackward(U[mu],mu,
-   		            Grid::QCD::CovShiftForward(U[nu],nu,U[mu])))) , mu, 1);
+	Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
+		            Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
+                            Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
+   		            Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[mu])))) , mu, 1);

 	TrStap = Grid::trace (U[mu]*Stap);

@ -390,10 +390,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	//           __ ___ 
 	//          |__    |

-	Stap = Grid::Cshift(Grid::QCD::CovShiftForward (U[nu],nu,
-		            Grid::QCD::CovShiftBackward(U[mu],mu,
-                            Grid::QCD::CovShiftBackward(U[mu],mu,
-                            Grid::QCD::CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1);
+	Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu,
+		            Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
+                            Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
+                            Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1);


 	TrStap = Grid::trace (U[mu]*Stap);
@ -412,12 +412,12 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	 * Make staple for loops centered at coor of link ; this one is ok.     //     |
 	 */
 	//	Stap = 
-	//	  Grid::Cshift(Grid::QCD::CovShiftForward(U[nu],nu,U[nu]),mu,1)* // ->||
-	//	  Grid::adj(Grid::QCD::CovShiftForward(U[nu],nu,Grid::QCD::CovShiftForward(U[nu],nu,U[mu]))) ;
-	Stap = Grid::Cshift(Grid::QCD::CovShiftForward(U[nu],nu,
-		            Grid::QCD::CovShiftForward(U[nu],nu,
-                            Grid::QCD::CovShiftBackward(U[mu],mu,
-                            Grid::QCD::CovShiftBackward(U[nu],nu,  Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
+	//	  Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]),mu,1)* // ->||
+	//	  Grid::adj(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[mu]))) ;
+	Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,
+		            Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,
+                            Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
+                            Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,  Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
 	  
 	TrStap = Grid::trace (U[mu]*Stap);
 	SumTrStap += TrStap;
@ -433,10 +433,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	//      |  | 
 	//       -- 

-	Stap = Grid::Cshift(Grid::QCD::CovShiftBackward(U[nu],nu,
-		            Grid::QCD::CovShiftBackward(U[nu],nu,
-                            Grid::QCD::CovShiftBackward(U[mu],mu,
-                            Grid::QCD::CovShiftForward (U[nu],nu,U[nu])))) , mu, 1);
+	Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
+		            Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
+                            Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
+                            Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu,U[nu])))) , mu, 1);

 	TrStap = Grid::trace (U[mu]*Stap);
 	trect = Grid::sum(TrStap);
@ -460,10 +460,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	Grid::QCD::LatticeColourMatrix tmp(UGrid);
 	
 	// 2 (mu)x1(nu)
-	left_2=  Grid::QCD::CovShiftForward(U[mu],mu,U[mu]);   // Umu(x) Umu(x+mu)
+	left_2=  Grid::QCD::PeriodicBC::CovShiftForward(U[mu],mu,U[mu]);   // Umu(x) Umu(x+mu)
 	tmp=Grid::Cshift(U[nu],mu,2);                          // Unu(x+2mu)

-	upper_l=  Grid::QCD::CovShiftForward(tmp,nu,Grid::adj(left_2)); //  Unu(x+2mu) Umu^dag(x+mu+nu) Umu^dag(x+nu) 
+	upper_l=  Grid::QCD::PeriodicBC::CovShiftForward(tmp,nu,Grid::adj(left_2)); //  Unu(x+2mu) Umu^dag(x+mu+nu) Umu^dag(x+nu) 
 	//                 __ __ 
 	//              =       |
 	
@ -533,9 +533,9 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	//   _
 	//  | |
 	//  | |
-	Grid::QCD::LatticeColourMatrix up2= Grid::QCD::CovShiftForward(U[nu],nu,U[nu]);
+	Grid::QCD::LatticeColourMatrix up2= Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]);

-	upper_l= Grid::QCD::CovShiftForward(Grid::Cshift(up2,mu,1),nu,Grid::Cshift(adj(U[mu]),nu,1));
+	upper_l= Grid::QCD::PeriodicBC::CovShiftForward(Grid::Cshift(up2,mu,1),nu,Grid::Cshift(adj(U[mu]),nu,1));
 	ds_U= upper_l*Grid::adj(up2);

 	RectPlaq_d = Grid::trace(U[mu]*ds_U);
@ -555,7 +555,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
   downer_l=           |  
               (x)<----V                 
 */    
-	down_l= Grid::adj(Grid::QCD::CovShiftForward(U[mu],mu,up2)); //downer_l
+	down_l= Grid::adj(Grid::QCD::PeriodicBC::CovShiftForward(U[mu],mu,up2)); //downer_l
 /*
                     ^     |
   down_staple  =    |     V 
@ -616,9 +616,9 @@ void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu)
    // Vol as for each site
    Grid::RealD RectScale(1.0/vol/12.0/6.0/3.0); 

-    Grid::QCD::WilsonLoops<Grid::QCD::LatticeGaugeField>::RectStaple(staple,Umu,mu);
+    Grid::QCD::ColourWilsonLoops::RectStaple(staple,Umu,mu);
    
-    link = Grid::QCD::PeekIndex<Grid::QCD::LorentzIndex>(Umu,mu);
+    link = Grid::QCD::PeekIndex<LorentzIndex>(Umu,mu);

    Traced = Grid::trace( link*staple) * RectScale;
    Grid::QCD::TComplex Tp = Grid::sum(Traced);
@ -655,9 +655,9 @@ void check_grid_p_staple(Grid::QCD::LatticeGaugeField & Umu)
    // Vol as for each site
    Grid::RealD Scale(1.0/vol/12.0/2.0/3.0); 

-    Grid::QCD::WilsonLoops<Grid::QCD::LatticeGaugeField>::Staple(staple,Umu,mu);
+    Grid::QCD::ColourWilsonLoops::Staple(staple,Umu,mu);
    
-    link = Grid::QCD::PeekIndex<Grid::QCD::LorentzIndex>(Umu,mu);
+    link = Grid::QCD::PeekIndex<LorentzIndex>(Umu,mu);

    Traced = Grid::trace( link*staple) * Scale;
    Grid::QCD::TComplex Tp = Grid::sum(Traced);
--- a/tests/qdpxx/Test_qdpxx_stag.cc
+++ b/tests/qdpxx/Test_qdpxx_stag.cc
@ -0,0 +1,364 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/qdpxx/Test_qdpxx_munprec.cc
+
+    Copyright (C) 2015
+
+Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+Author: paboyle <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>
+
+double mq=0.1;
+
+typedef Grid::QCD::StaggeredImplR::FermionField FermionField;
+typedef Grid::QCD::LatticeGaugeField GaugeField;
+
+void make_gauge     (GaugeField & lat, FermionField &src);
+void calc_grid      (GaugeField & lat, GaugeField & uthin,GaugeField & ufat, FermionField &src, FermionField &res,int dag);
+void calc_chroma    (GaugeField & lat,GaugeField & uthin,GaugeField & ufat, FermionField &src, FermionField &res,int dag);
+
+#include <chroma.h>
+#include <actions/ferm/invert/syssolver_linop_cg_array.h>
+#include <actions/ferm/invert/syssolver_linop_aggregate.h>
+
+namespace Chroma { 
+
+
+class ChromaWrapper {
+public:
+  
+  typedef multi1d<LatticeColorMatrix> U;
+  typedef LatticeStaggeredFermion T4;
+  
+  static void ImportGauge(GaugeField & gr,
+			  QDP::multi1d<QDP::LatticeColorMatrix> & ch) 
+  {
+    Grid::QCD::LorentzColourMatrix LCM;
+    Grid::Complex cc;
+    QDP::ColorMatrix cm;
+    QDP::Complex c;
+
+    std::vector<int> x(4);
+    QDP::multi1d<int> cx(4);
+    std::vector<int> gd= gr._grid->GlobalDimensions();
+
+    for (x[0]=0;x[0]<gd[0];x[0]++){
+    for (x[1]=0;x[1]<gd[1];x[1]++){
+    for (x[2]=0;x[2]<gd[2];x[2]++){
+    for (x[3]=0;x[3]<gd[3];x[3]++){
+      cx[0] = x[0];
+      cx[1] = x[1];
+      cx[2] = x[2];
+      cx[3] = x[3];
+      Grid::peekSite(LCM,gr,x);
+
+      for(int mu=0;mu<4;mu++){
+	for(int i=0;i<3;i++){
+	for(int j=0;j<3;j++){
+	  cc = LCM(mu)()(i,j);
+	  c = QDP::cmplx(QDP::Real(real(cc)),QDP::Real(imag(cc)));
+	  QDP::pokeColor(cm,c,i,j);
+	}}
+	QDP::pokeSite(ch[mu],cm,cx);
+      }
+
+    }}}}
+  }
+
+  static void ExportGauge(GaugeField & gr,
+			  QDP::multi1d<QDP::LatticeColorMatrix> & ch) 
+  {
+    Grid::QCD::LorentzColourMatrix LCM;
+    Grid::Complex cc;
+    QDP::ColorMatrix cm;
+    QDP::Complex c;
+
+    std::vector<int> x(4);
+    QDP::multi1d<int> cx(4);
+    std::vector<int> gd= gr._grid->GlobalDimensions();
+
+    for (x[0]=0;x[0]<gd[0];x[0]++){
+    for (x[1]=0;x[1]<gd[1];x[1]++){
+    for (x[2]=0;x[2]<gd[2];x[2]++){
+    for (x[3]=0;x[3]<gd[3];x[3]++){
+      cx[0] = x[0];
+      cx[1] = x[1];
+      cx[2] = x[2];
+      cx[3] = x[3];
+
+      for(int mu=0;mu<4;mu++){
+	for(int i=0;i<3;i++){
+	for(int j=0;j<3;j++){
+	  cm = QDP::peekSite(ch[mu],cx);
+	  c  = QDP::peekColor(cm,i,j);
+	  cc = Grid::Complex(toDouble(real(c)),toDouble(imag(c)));
+	  LCM(mu)()(i,j)= cc;
+	}}
+      }
+      Grid::pokeSite(LCM,gr,x);
+
+    }}}}
+  }
+
+  
+  static void ImportFermion(FermionField & gr,
+			    QDP::LatticeStaggeredFermion & ch  ) 
+  {
+    Grid::QCD::ColourVector F;
+    Grid::Complex c;
+
+
+    std::vector<int> x(5);
+    QDP::multi1d<int> cx(4);
+    std::vector<int> gd= gr._grid->GlobalDimensions();
+
+    for (x[0]=0;x[0]<gd[0];x[0]++){
+    for (x[1]=0;x[1]<gd[1];x[1]++){
+    for (x[2]=0;x[2]<gd[2];x[2]++){
+    for (x[3]=0;x[3]<gd[3];x[3]++){
+      cx[0] = x[0];
+      cx[1] = x[1];
+      cx[2] = x[2];
+      cx[3] = x[3];
+
+      Grid::peekSite(F,gr,x);
+      QDP::ColorVector cv;
+      for(int j=0;j<3;j++){
+	QDP::Complex cc;
+	c  = F()()(j) ;
+	cc = QDP::cmplx(QDP::Real(real(c)),QDP::Real(imag(c)));
+	pokeColor(cv,cc,j);
+      }
+      QDP::StaggeredFermion cF;
+      pokeSpin(cF,cv,0);
+      QDP::pokeSite(ch,cF,cx);
+    }}}}
+  }
+  static void ExportFermion(FermionField & gr,
+			    QDP::LatticeStaggeredFermion & ch  ) 
+  {
+    Grid::QCD::ColourVector F;
+    Grid::Complex c;
+
+    std::vector<int> x(5);
+    QDP::multi1d<int> cx(4);
+    std::vector<int> gd= gr._grid->GlobalDimensions();
+
+    for (x[0]=0;x[0]<gd[0];x[0]++){
+    for (x[1]=0;x[1]<gd[1];x[1]++){
+    for (x[2]=0;x[2]<gd[2];x[2]++){
+    for (x[3]=0;x[3]<gd[3];x[3]++){
+      cx[0] = x[0];
+      cx[1] = x[1];
+      cx[2] = x[2];
+      cx[3] = x[3];
+
+      QDP::StaggeredFermion cF = QDP::peekSite(ch,cx);
+      for(int j=0;j<3;j++){
+	QDP::ColorVector cS=QDP::peekSpin(cF,0);
+	QDP::Complex cc=QDP::peekColor(cS,j);
+	c = Grid::Complex(QDP::toDouble(QDP::real(cc)), 
+			  QDP::toDouble(QDP::imag(cc)));
+	F()()(j) = c;
+      }
+      Grid::pokeSite(F,gr,x);
+    }}}}
+  }
+
+  static Handle< Chroma::EvenOddLinearOperator<T4,U,U> >  GetLinOp (U &u,U &u_fat,U &u_triple)
+  {
+    QDP::Real _mq(mq);
+    QDP::multi1d<int> bcs(QDP::Nd);
+
+    bcs[0] = bcs[1] = bcs[2] = bcs[3] = 1;
+
+    Chroma::AsqtadFermActParams p; 
+    p.Mass = _mq; 
+    p.u0 = Real(1.0);
+
+
+    Chroma::Handle<Chroma::FermBC<T4,U,U> > fbc(new Chroma::SimpleFermBC< T4, U, U >(bcs));
+    Chroma::Handle<Chroma::CreateFermState<T4,U,U> > cfs( new Chroma::CreateSimpleFermState<T4,U,U>(fbc));
+    Chroma::AsqtadFermAct S_f(cfs,p);
+    Chroma::Handle< Chroma::FermState<T4,U,U> >  ffs(  S_f.createState(u) );
+    u_fat   =ffs.cast<AsqtadConnectStateBase>()->getFatLinks();
+    u_triple=ffs.cast<AsqtadConnectStateBase>()->getTripleLinks();
+    return S_f.linOp(ffs);
+  }
+
+};
+}
+
+int main (int argc,char **argv )
+{
+
+  /********************************************************
+   * Setup QDP
+   *********************************************************/
+  Chroma::initialize(&argc,&argv);
+  Chroma::WilsonTypeFermActs4DEnv::registerAll(); 
+
+  /********************************************************
+   * Setup Grid
+   *********************************************************/
+  Grid::Grid_init(&argc,&argv);
+  Grid::GridCartesian * UGrid   = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(), 
+									    Grid::GridDefaultSimd(Grid::QCD::Nd,Grid::vComplex::Nsimd()),
+									    Grid::GridDefaultMpi());
+  
+  std::vector<int> gd = UGrid->GlobalDimensions();
+  QDP::multi1d<int> nrow(QDP::Nd);
+  for(int mu=0;mu<4;mu++) nrow[mu] = gd[mu];
+
+  QDP::Layout::setLattSize(nrow);
+  QDP::Layout::create();
+
+  GaugeField uthin  (UGrid);
+  GaugeField ufat   (UGrid);
+  GaugeField utriple(UGrid);
+  FermionField    src(UGrid);
+  FermionField    res_chroma(UGrid);
+  FermionField    res_grid  (UGrid);
+  
+
+  {
+
+    std::cout << "*****************************"<<std::endl;
+    std::cout << "Staggered Action "            <<std::endl;
+    std::cout << "*****************************"<<std::endl;
+
+    make_gauge(uthin,src);
+
+    for(int dag=0;dag<2;dag++) {
+
+      std::cout << "Dag =  "<<dag<<std::endl;
+      
+      calc_chroma(uthin,utriple,ufat,src,res_chroma,dag);
+
+      // Remove the normalisation of Chroma Gauge links ??
+      std::cout << "Norm of chroma Asqtad multiply "<<Grid::norm2(res_chroma)<<std::endl;
+      calc_grid  (uthin,utriple,ufat,src,res_grid,dag);
+
+      std::cout << "Norm of thin gauge "<< Grid::norm2(uthin) <<std::endl;
+      std::cout << "Norm of fat  gauge "<< Grid::norm2(ufat) <<std::endl;
+
+      std::cout << "Norm of Grid Asqtad multiply "<<Grid::norm2(res_grid)<<std::endl;
+      
+      /*
+      std::cout << " site 0 of Uthin  "<<uthin._odata[0] <<std::endl;
+      std::cout << " site 0 of Utriple"<<utriple._odata[0] <<std::endl;
+      std::cout << " site 0 of Ufat   "<<ufat._odata[0] <<std::endl;
+
+      std::cout << " site 0 of Grid   "<<res_grid._odata[0] <<std::endl;
+      std::cout << " site 0 of Chroma "<<res_chroma._odata[0] <<std::endl;
+      */
+
+      res_chroma=res_chroma - res_grid;
+      std::cout << "Norm of difference "<<Grid::norm2(res_chroma)<<std::endl;
+    }
+  }
+
+  std::cout << "Finished test "<<std::endl;
+
+  Chroma::finalize();
+}
+
+void calc_chroma(GaugeField & lat, GaugeField &uthin, GaugeField &ufat, FermionField &src, FermionField &res,int dag)
+{
+  typedef QDP::LatticeStaggeredFermion T;
+  typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
+  
+  U u(4);
+  U ut(4);
+  U uf(4);
+
+  //  Chroma::HotSt(u);
+  Chroma::ChromaWrapper::ImportGauge(lat,u) ;
+
+  QDP::LatticeStaggeredFermion  check;
+  QDP::LatticeStaggeredFermion  result;
+  QDP::LatticeStaggeredFermion  tmp;
+  QDP::LatticeStaggeredFermion  psi;
+
+  Chroma::ChromaWrapper::ImportFermion(src,psi);
+
+  auto linop =Chroma::ChromaWrapper::GetLinOp(u,uf,ut);
+
+  Chroma::ChromaWrapper::ExportGauge(uthin,ut) ;
+  Chroma::ChromaWrapper::ExportGauge(ufat ,uf) ;
+
+  enum Chroma::PlusMinus isign;
+  if ( dag ) {
+    isign=Chroma::MINUS;
+  } else {
+    isign=Chroma::PLUS;
+  }
+
+  std::cout << "Calling Chroma Linop "<< std::endl;
+  linop->evenEvenLinOp(tmp,psi,isign); check[rb[0]] = tmp;
+  linop->oddOddLinOp  (tmp,psi,isign); check[rb[1]] = tmp;
+  linop->evenOddLinOp(tmp,psi,isign) ; check[rb[0]]+= tmp;
+  linop->oddEvenLinOp(tmp,psi,isign) ; check[rb[1]]+= tmp;
+
+  Chroma::ChromaWrapper::ExportFermion(res,check) ;
+}
+
+
+void make_gauge(GaugeField & Umu,FermionField &src)
+{
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  std::vector<int> seeds4({1,2,3,4});
+
+  Grid::GridCartesian         * UGrid   = (Grid::GridCartesian *) Umu._grid;
+  Grid::GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+  Grid::QCD::SU3::HotConfiguration(RNG4,Umu);
+  Grid::gaussian(RNG4,src);
+}
+
+void calc_grid(GaugeField & Uthin, GaugeField & Utriple, GaugeField & Ufat, FermionField &src, FermionField &res,int dag)
+{
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  Grid::GridCartesian         * UGrid   = (Grid::GridCartesian *) Uthin._grid;
+  Grid::GridRedBlackCartesian * UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+  Grid::QCD::ImprovedStaggeredFermionR Dstag(Uthin,Utriple,Ufat,*UGrid,*UrbGrid,mq*2.0);
+
+  std::cout << Grid::GridLogMessage <<" Calling Grid staggered multiply "<<std::endl;
+
+  if ( dag ) 
+    Dstag.Mdag(src,res);  
+  else 
+    Dstag.M(src,res);  
+
+  res = res ; // Convention mismatch to Chroma
+  return;
+} 
+
+
+
+
+
--- a/tests/solver/Test_zmobius_cg_prec.cc
+++ b/tests/solver/Test_zmobius_cg_prec.cc
@ -43,7 +43,7 @@ Gamma::Algebra Gmu[] = {Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::A
 int main(int argc, char** argv) {
  Grid_init(&argc, &argv);

-  const int Ls = 16;
+  const int Ls = 10;

  GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
      GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
@ -80,11 +80,27 @@ int main(int argc, char** argv) {
  RealD mass = 0.01;
  RealD M5 = 1.8;
  std::vector < std::complex<double>  > omegas;
+#if 0
  for(int i=0;i<Ls;i++){
-  	std::complex<double> temp (0.25+0.00*i, 0.0+0.00*i);
- 	 omegas.push_back(temp);
+    double imag = 0.;
+    if (i==0) imag=1.;
+    if (i==Ls-1) imag=-1.;
+    std::complex<double> temp (0.25+0.01*i, imag*0.01);
+    omegas.push_back(temp);
  }
-//  DomainWallFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5);
+#else
+  omegas.push_back( std::complex<double>(1.45806438985048,-0) );
+  omegas.push_back( std::complex<double>(1.18231318389348,-0) );
+  omegas.push_back( std::complex<double>(0.830951166685955,-0) );
+  omegas.push_back( std::complex<double>(0.542352409156791,-0) );
+  omegas.push_back( std::complex<double>(0.341985020453729,-0) );
+  omegas.push_back( std::complex<double>(0.21137902619029,-0) );
+  omegas.push_back( std::complex<double>(0.126074299502912,-0) );
+  omegas.push_back( std::complex<double>(0.0990136651962626,-0) );
+  omegas.push_back( std::complex<double>(0.0686324988446592,0.0550658530827402) );
+  omegas.push_back( std::complex<double>(0.0686324988446592,-0.0550658530827402) );
+#endif
+
  ZMobiusFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, omegas,1.,0.);

  LatticeFermion src_o(FrbGrid);
--- a/tests/testu01/Makefile.am
+++ b/tests/testu01/Makefile.am
@ -0,0 +1,3 @@
+AM_LDFLAGS += -L$(LIBRARY_PATH) -ltestu01 -lprobdist -lmylib -lm
+AM_CXXFLAGS += -I$(C_INCLUDE_PATH)
+include Make.inc
--- a/tests/testu01/Test_smallcrush.cc
+++ b/tests/testu01/Test_smallcrush.cc
@ -0,0 +1,175 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_smallcrush.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>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+// Wrap Grid's parallel RNG for testU01
+#undef BIG_CRUSH             // Big crush enable (long running)
+#define MIDDLE_CRUSH             // Big crush enable (long running)
+#undef SMALL_CRUSH             // Big crush enable (long running)
+#undef TEST_RNG_STANDALONE   // Test serial RNGs in isolation
+
+extern "C" { 
+#include "TestU01.h"
+}
+
+std::vector<std::ranlux48>      EngineRanlux;
+std::vector<std::mt19937>       EngineMT;
+
+#include <Grid/sitmo_rng/sitmo_prng_engine.hpp>
+std::vector<sitmo::prng_engine> EngineSitmo;
+
+std::uniform_int_distribution<uint32_t> uid;
+
+uint32_t GetU01Ranlux(void) {
+  return uid(EngineRanlux[0]);
+};
+uint32_t GetU01MT(void) {
+  return uid(EngineMT[0]);
+};
+uint32_t GetU01Sitmo(void) {
+  return uid(EngineSitmo[0]);
+};
+
+typedef Grid::GridRNGbase::RngEngine RngEngine;
+
+struct TestRNG { 
+public:
+  static GridParallelRNG *pRNG;
+  static GridSerialRNG *sRNG;
+  static GridBase *_grid;
+  static RngEngine Eng;
+  static uint64_t site;
+  static uint64_t gsites;
+  static char *name;
+
+  static void Init(GridParallelRNG *_pRNG,GridSerialRNG *_sRNG,GridBase *grid) {
+    pRNG = _pRNG;
+    sRNG = _sRNG;
+    _grid= grid;
+    gsites= grid->_gsites;
+    site = 0;
+  }
+  static uint32_t GetU01(void) { 
+    uint32_t ret_val;
+    ret_val = pRNG->GlobalU01(site);
+    site=(site+1)%gsites;
+    return ret_val;
+  }
+};
+
+GridParallelRNG *TestRNG::pRNG;
+GridSerialRNG   *TestRNG::sRNG;
+GridBase        *TestRNG::_grid;
+RngEngine        TestRNG::Eng;
+uint64_t         TestRNG::site;
+uint64_t         TestRNG::gsites;
+
+#ifdef RNG_SITMO
+char * TestRNG::name = (char *)"Grid_Sitmo";
+#endif
+#ifdef RNG_RANLUX
+char * TestRNG::name = (char *)"Grid_ranlux48";
+#endif
+#ifdef RNG_MT19937
+char * TestRNG::name = (char *)"Grid_mt19937";
+#endif
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+     
+  GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
+
+  std::vector<int> seeds({1,2,3,4});
+  std::seed_seq seq(seeds.begin(),seeds.end());
+
+  EngineRanlux.push_back(std::ranlux48(seq));
+  EngineMT.push_back(std::mt19937(seq));
+  EngineSitmo.push_back(sitmo::prng_engine(seq));
+
+  std::cout << GridLogMessage<< "Initialising Grid RNGs "<<std::endl; 
+  GridParallelRNG           pRNG(&Grid);   
+  pRNG.SeedFixedIntegers(std::vector<int>({43,12,7019,9}));
+  GridSerialRNG           sRNG;
+  sRNG.SeedFixedIntegers(std::vector<int>({102,12,99,15}));
+  std::cout << GridLogMessage<< "Initialised Grid RNGs "<<std::endl; 
+
+  TestRNG::Init(&pRNG,&sRNG,&Grid);
+  std::cout << GridLogMessage<< "Grid RNG's are "<< std::string(TestRNG::name) <<std::endl; 
+
+  unif01_Gen * gen;
+
+#ifdef TEST_RNG_STANDALONE
+  std::cout << GridLogMessage<< "Testing Standalone Ranlux" <<std::endl; 
+  gen = unif01_CreateExternGenBits ((char *)"GridRanlux",GetU01Ranlux);
+  bbattery_SmallCrush (gen);
+  unif01_DeleteExternGenBits(gen);
+  std::cout << GridLogMessage<< "Testing Standalone Ranlux is complete" <<std::endl; 
+
+  std::cout << GridLogMessage<< "Testing Standalone Mersenne Twister" <<std::endl; 
+  gen = unif01_CreateExternGenBits ((char *)"GridMT",GetU01MT);
+  bbattery_SmallCrush (gen);
+  unif01_DeleteExternGenBits(gen);
+  std::cout << GridLogMessage<< "Testing Standalone Mersenne Twister is complete" <<std::endl; 
+
+  std::cout << GridLogMessage<< "Testing Standalone Sitmo" <<std::endl; 
+  gen = unif01_CreateExternGenBits ((char *)"GridSitmo",GetU01Sitmo);
+  bbattery_SmallCrush (gen);
+  unif01_DeleteExternGenBits(gen);
+  std::cout << GridLogMessage<< "Testing Standalone Sitmo is complete" <<std::endl; 
+#endif
+
+#ifdef BIG_CRUSH
+  std::cout << GridLogMessage<< "Testing Grid BigCrush for "<< std::string(TestRNG::name) <<std::endl; 
+  gen = unif01_CreateExternGenBits(TestRNG::name,TestRNG::GetU01);
+  bbattery_BigCrush (gen);
+  std::cout << GridLogMessage<< "Testing Grid BigCrush "<< std::string(TestRNG::name)<<" is complete" <<std::endl; 
+#endif
+#ifdef MIDDLE_CRUSH
+  std::cout << GridLogMessage<< "Testing Grid Crush for "<< std::string(TestRNG::name) <<std::endl; 
+  gen = unif01_CreateExternGenBits(TestRNG::name,TestRNG::GetU01);
+  bbattery_Crush (gen);
+  std::cout << GridLogMessage<< "Testing Grid Crush "<< std::string(TestRNG::name)<<" is complete" <<std::endl; 
+#endif
+#ifdef SMALL_CRUSH
+  std::cout << GridLogMessage<< "Testing Grid SmallCrush for "<< std::string(TestRNG::name) <<std::endl; 
+  gen = unif01_CreateExternGenBits(TestRNG::name,TestRNG::GetU01);
+  bbattery_SmallCrush (gen);
+  std::cout << GridLogMessage<< "Testing Grid SmallCrush "<< std::string(TestRNG::name)<<" is complete" <<std::endl; 
+#endif
+  Grid_finalize();
+}
+
Author	SHA1	Message	Date
paboyle	5592f7b8c1	Creation mode better implementation	2017-04-05 02:35:34 +09:00
paboyle	35da4ece0b	UID fix	2017-04-05 02:18:15 +09:00
paboyle	061b15b9e9	Merge branch 'feature/sitmo-skipahead' into develop	2017-04-05 01:24:49 +09:00
paboyle	561426f6eb	Clean up	2017-04-02 23:13:48 +09:00
paboyle	83f6fab8fa	Big/Small crush test, and fast SITMO rng init, faster but not ideal MT and Ranlux init.	2017-04-02 12:10:51 +09:00
paboyle	0fade84ab2	No random device	2017-04-02 00:29:40 +09:00
paboyle	9dc7ca4c3b	Sitmo fast init	2017-04-02 00:28:22 +09:00
paboyle	935d82f5b1	sanity checks	2017-04-02 00:27:28 +09:00
paboyle	9cbcdd65d7	No random device seed	2017-04-02 00:26:57 +09:00
paboyle	f18f5ed926	Drop random device	2017-04-02 00:26:26 +09:00
paboyle	d1d63a4f2d	sitmo default	2017-04-02 00:26:05 +09:00
paboyle	7e5faa0f34	Multiple RNGs	2017-04-02 00:25:44 +09:00
paboyle	6af459cae4	Christoph's coefficients.	2017-03-31 17:07:43 +09:00
paboyle	1c4bc7ed38	Debugged staggered conventions	2017-03-31 14:41:48 +09:00
paboyle	93ea5d9468	Pretty code	2017-03-30 15:00:03 +09:00
paboyle	1ec5d32369	Chulwoo's test to zmobius helped me shake out	2017-03-30 13:45:13 +09:00
paboyle	9fd23faadf	Pretty layout	2017-03-30 13:44:45 +09:00
paboyle	10e4fa0dc8	Template instantiation improvements	2017-03-30 13:44:25 +09:00
paboyle	c4aca1dde4	Conjugate coefficients on adjoint	2017-03-30 13:44:05 +09:00
paboyle	b9e8ea3aaa	conjugate coefficient on the dagger	2017-03-30 13:43:13 +09:00
paboyle	077aa728b9	Fix the ZMobius (I think)	2017-03-30 13:42:09 +09:00
paboyle	a8d83d886e	Macro controls	2017-03-30 13:31:34 +09:00
paboyle	7fd46eeec4	Trailing whitespace removal	2017-03-30 13:31:10 +09:00
paboyle	e0c4eeb3ec	Compiles again	2017-03-30 13:30:45 +09:00
paboyle	cb9a297a0a	Chulwoo's Zmobius test	2017-03-30 13:30:25 +09:00
paboyle	2b115929dc	Small AVX512 asm ifdef patch	2017-03-29 18:51:23 +09:00
paboyle	5c6571dab1	Merge branch 'feature/bgq-asm' into develop	2017-03-29 18:48:55 +09:00