HMC bit repro across checkpoints. Fixed parallel RNG issue with threading.

Conclusion: c++11 distributions not thread safe and must us distinct dist as well as distinct engine
per site. Makes sense when you think of box muller. Also added a reset of dist on fill to ensure
repro across checkpoints.

lib/algorithms/iterative/ConjugateGradient.h

@@ -30,7 +30,7 @@ public:
       
       //Initial residual computation & set up
       RealD guess = norm2(psi);
-      assert(std::isnan(psi)==0);
+      assert(std::isnan(guess)==0);
 
       Linop.HermOpAndNorm(psi,mmp,d,b);
       

lib/lattice/Lattice_rng.h

@@ -22,13 +22,16 @@ namespace Grid {
       assert(fine->_processors[d]==1);
     }
 
-    // local and global volumes subdivide cleanly after SIMDization
     int multiplicity=1;
+    for(int d=0;d<lowerdims;d++){
+      multiplicity=multiplicity*fine->_rdimensions[d];
+    }
+    // local and global volumes subdivide cleanly after SIMDization
     for(int d=0;d<rngdims;d++){
       int fd= d+lowerdims;
       assert(coarse->_processors[d]  == fine->_processors[fd]);
       assert(coarse->_simd_layout[d] == fine->_simd_layout[fd]);
-      assert((fine->_rdimensions[fd] / coarse->_rdimensions[d])* coarse->_rdimensions[d]==fine->_rdimensions[fd]); 
+      assert(((fine->_rdimensions[fd] / coarse->_rdimensions[d])* coarse->_rdimensions[d])==fine->_rdimensions[fd]); 
 
       multiplicity = multiplicity *fine->_rdimensions[fd] / coarse->_rdimensions[d]; 
     }
@@ -76,8 +79,6 @@ namespace Grid {
 
   public:
 
-   GridRNGbase() : _uniform{0,1}, _gaussian(0.0,1.0) {};
-
     int _seeded;
     // One generator per site.
     // Uniform and Gaussian distributions from these generators.
@@ -91,13 +92,12 @@ namespace Grid {
     static const int     RngStateCount = std::mt19937::state_size;
 #endif
     std::vector<RngEngine>             _generators;
-    std::uniform_real_distribution<double> _uniform;
+    std::vector<std::uniform_real_distribution<RealD> > _uniform;
-    std::normal_distribution<double>       _gaussian;
+    std::vector<std::normal_distribution<RealD> >       _gaussian;
 
     void GetState(std::vector<RngStateType> & saved,int gen) {
       saved.resize(RngStateCount);
       std::stringstream ss;
-      //      std::cout << _generators[gen]<<std::endl;
       ss<<_generators[gen];
       ss.seekg(0,ss.beg);
       for(int i=0;i<RngStateCount;i++){
@@ -112,7 +112,6 @@ namespace Grid {
       }
       ss.seekg(0,ss.beg);
       ss>>_generators[gen];
-      //      std::cout << _generators[gen]<<std::endl;
     }
   };
 
@@ -132,12 +131,14 @@ namespace Grid {
 
     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) );
       _seeded=0;
     }
 
 
 
-    template <class sobj,class distribution> inline void fill(sobj &l,distribution &dist){
+    template <class sobj,class distribution> inline void fill(sobj &l,std::vector<distribution> &dist){
 
       typedef typename sobj::scalar_type scalar_type;
  
@@ -145,56 +146,65 @@ namespace Grid {
 
       scalar_type *buf = (scalar_type *) & l;
 
+      dist[0].reset();
       for(int idx=0;idx<words;idx++){
-	fillScalar(buf[idx],dist,_generators[0]);
+	fillScalar(buf[idx],dist[0],_generators[0]);
       }
       
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
 
     };
 
-    template <class distribution>  inline void fill(ComplexF &l,distribution &dist){
+    template <class distribution>  inline void fill(ComplexF &l,std::vector<distribution> &dist){
-      fillScalar(l,dist,_generators[0]);
+      dist[0].reset();
+      fillScalar(l,dist[0],_generators[0]);
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
     }
-    template <class distribution>  inline void fill(ComplexD &l,distribution &dist){
+    template <class distribution>  inline void fill(ComplexD &l,std::vector<distribution> &dist){
-      fillScalar(l,dist,_generators[0]);
+      dist[0].reset();
+      fillScalar(l,dist[0],_generators[0]);
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
     }
-    template <class distribution>  inline void fill(RealF &l,distribution &dist){
+    template <class distribution>  inline void fill(RealF &l,std::vector<distribution> &dist){
-      fillScalar(l,dist,_generators[0]);
+      dist[0].reset();
+      fillScalar(l,dist[0],_generators[0]);
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
     }
-    template <class distribution>  inline void fill(RealD &l,distribution &dist){
+    template <class distribution>  inline void fill(RealD &l,std::vector<distribution> &dist){
-      fillScalar(l,dist,_generators[0]);
+      dist[0].reset();
+      fillScalar(l,dist[0],_generators[0]);
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
     }
     // vector fill
-    template <class distribution>  inline void fill(vComplexF &l,distribution &dist){
+    template <class distribution>  inline void fill(vComplexF &l,std::vector<distribution> &dist){
       RealF *pointer=(RealF *)&l;
+      dist[0].reset();
       for(int i=0;i<2*vComplexF::Nsimd();i++){
-	fillScalar(pointer[i],dist,_generators[0]);
+	fillScalar(pointer[i],dist[0],_generators[0]);
       }
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
     }
-    template <class distribution>  inline void fill(vComplexD &l,distribution &dist){
+    template <class distribution>  inline void fill(vComplexD &l,std::vector<distribution> &dist){
       RealD *pointer=(RealD *)&l;
+      dist[0].reset();
       for(int i=0;i<2*vComplexD::Nsimd();i++){
-	fillScalar(pointer[i],dist,_generators[0]);
+	fillScalar(pointer[i],dist[0],_generators[0]);
       }
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
     }
-    template <class distribution>  inline void fill(vRealF &l,distribution &dist){
+    template <class distribution>  inline void fill(vRealF &l,std::vector<distribution> &dist){
       RealF *pointer=(RealF *)&l;
+      dist[0].reset();
       for(int i=0;i<vRealF::Nsimd();i++){
-	fillScalar(pointer[i],dist,_generators[0]);
+	fillScalar(pointer[i],dist[0],_generators[0]);
       }
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
     }
-    template <class distribution>  inline void fill(vRealD &l,distribution &dist){
+    template <class distribution>  inline void fill(vRealD &l,std::vector<distribution> &dist){
       RealD *pointer=(RealD *)&l;
+      dist[0].reset();
       for(int i=0;i<vRealD::Nsimd();i++){
-	fillScalar(pointer[i],dist,_generators[0]);
+	fillScalar(pointer[i],dist[0],_generators[0]);
       }
       CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
     }
@@ -213,10 +223,6 @@ namespace Grid {
 
   class GridParallelRNG : public GridRNGbase {
   public:
-    
-    // Uniform and Gaussian distributions from these generators.
-    std::uniform_real_distribution<double> _uniform;
-    std::normal_distribution<double>       _gaussian;
 
     GridBase *_grid;
     int _vol;
@@ -228,7 +234,10 @@ namespace Grid {
     GridParallelRNG(GridBase *grid) : GridRNGbase() {
       _grid=grid;
       _vol =_grid->iSites()*_grid->oSites();
+
       _generators.resize(_vol);
+      _uniform.resize(_vol,std::uniform_real_distribution<RealD>{0,1});
+      _gaussian.resize(_vol,std::normal_distribution<RealD>(0.0,1.0) );
       _seeded=0;
     }
 
@@ -277,7 +286,7 @@ namespace Grid {
     //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,distribution &dist){
+    template <class vobj,class distribution> inline void fill(Lattice<vobj> &l,std::vector<distribution> &dist){
 
       typedef typename vobj::scalar_object scalar_object;
       typedef typename vobj::scalar_type scalar_type;
@@ -301,8 +310,9 @@ PARALLEL_FOR_LOOP
 	  for(int si=0;si<Nsimd;si++){
 	    int gdx = generator_idx(ss,si); // index of generator state
 	    scalar_type *pointer = (scalar_type *)&buf[si];
+	    dist[gdx].reset();
 	    for(int idx=0;idx<words;idx++){
-	      fillScalar(pointer[idx],dist,_generators[gdx]);
+	      fillScalar(pointer[idx],dist[gdx],_generators[gdx]);
 	    }
 	  }
 

lib/qcd/action/gauge/WilsonGaugeAction.h

@@ -22,10 +22,10 @@ namespace Grid{
       
       virtual RealD S(const GaugeField &U) {
 	RealD plaq = WilsonLoops<GaugeField>::avgPlaquette(U);
-	std::cout<<GridLogMessage << "Plaq : "<<plaq << "\n";
+	//	std::cout<<GridLogMessage << "Plaq : "<<plaq << "\n";
 	RealD vol = U._grid->gSites();
 	RealD action=beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
-	std::cout << GridLogMessage << "WilsonGauge action "<<action<<std::endl;
+	// std::cout << GridLogMessage << "WilsonGauge action "<<action<<std::endl;
 	return action;
       };
 

lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h

@@ -61,7 +61,7 @@ namespace Grid{
 	FermionField tmp    (NumOp.FermionRedBlackGrid());
 
 	gaussian(pRNG,eta);
-	
+
 	pickCheckerboard(Even,etaEven,eta);
 	pickCheckerboard(Odd,etaOdd,eta);
 

lib/qcd/hmc/HMC.h

@@ -39,7 +39,25 @@ namespace Grid{
       virtual void TrajectoryComplete (int traj, GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG & pRNG )=0;
     };
 
-    
+    template<class GaugeField> 
+    class PlaquetteLogger : public HmcObservable<GaugeField> {
+    private:
+      std::string Stem;
+    public:
+      PlaquetteLogger(std::string cf) {
+        Stem  = cf;
+      };
+
+      void TrajectoryComplete(int traj, GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG & pRNG )
+      {
+	  std::string file;   { std::ostringstream os; os << Stem     <<"."<< traj; file = os.str(); }
+	  std::ofstream of(file);
+	  RealD plaq = WilsonLoops<GaugeField>::avgPlaquette(U);
+	  of << plaq <<std::endl;
+	  std::cout<< GridLogMessage<< "Plaquette for trajectory "<< traj << " is " << plaq <<std::endl;
+      }
+    };
+
     //    template <class GaugeField, class Integrator, class Smearer, class Boundary> 
     template <class GaugeField, class IntegratorType>
     class HybridMonteCarlo {
@@ -141,9 +159,6 @@ namespace Grid{
 	    Ucur = Ucopy;
 	  }
 
-	  plaq = WilsonLoops<GaugeField>::avgPlaquette(Ucur);
-	  std::cout << " Now gauge field has plaq = "<< plaq <<std::endl;
-
 	  for(int obs = 0;obs<Observables.size();obs++){
 	    Observables[obs]->TrajectoryComplete (traj+1,Ucur,sRNG,pRNG);
 	  }

lib/qcd/hmc/integrators/Integrator.h

@@ -165,7 +165,6 @@ namespace Grid{
 	    H += Hterm;
 	  }
 	}
-	std::cout<<GridLogMessage << "Total action H = "<< H << "\n";
 	
 	return H;
       }

tests/Test_hmc_EODWFRatio.cc

@@ -1,62 +1,158 @@
 #include "Grid.h"
-
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
+namespace Grid { 
+  namespace QCD { 
+
+class NerscHmcRunner {
+public:
+
+  enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart };
+
+  ActionSet<LatticeGaugeField> TheAction;
+
+  GridCartesian         * UGrid   ;
+  GridCartesian         * FGrid   ;
+  GridRedBlackCartesian * UrbGrid ;
+  GridRedBlackCartesian * FrbGrid ;
+
+  void BuildTheAction (int argc, char **argv)
+  {
+    const int Ls = 8;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  
+    FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+    FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(5.6);
+
+    Real mass=0.04;
+    Real pv  =1.0;
+    RealD M5=1.5;
+    DomainWallFermionR DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+    DomainWallFermionR NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5);
+  
+    ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
+    TwoFlavourEvenOddRatioPseudoFermionAction<WilsonImplR> Nf2(NumOp, DenOp,CG,CG);
+  
+    //Collect actions
+    ActionLevel<LatticeGaugeField> Level1;
+    Level1.push_back(&Nf2);
+    Level1.push_back(&Waction);
+    TheAction.push_back(Level1);
+
+    Run(argc,argv);
+  };
+  
+  void Run (int argc, char  **argv){
+
+    StartType_t StartType = HotStart;
+
+    std::string arg;
+
+    if( GridCmdOptionExists(argv,argv+argc,"--StartType") ){
+      arg = GridCmdOptionPayload(argv,argv+argc,"--StartType");
+      if ( arg == "HotStart" ) { StartType = HotStart; }
+      else if ( arg == "ColdStart" ) { StartType = ColdStart; }
+      else if ( arg == "TepidStart" ) { StartType = TepidStart; }
+      else if ( arg == "CheckpointStart" ) { StartType = CheckpointStart; }
+      else assert(0);
+    }
+
+    int StartTraj = 0;
+    if( GridCmdOptionExists(argv,argv+argc,"--StartTrajectory") ){
+      arg= GridCmdOptionPayload(argv,argv+argc,"--StartTrajectory");
+      std::vector<int> ivec(0);
+      GridCmdOptionIntVector(arg,ivec);
+      StartTraj = ivec[0];
+    }    
+
+    int NumTraj = 1;
+    if( GridCmdOptionExists(argv,argv+argc,"--Trajectories") ){
+      arg= GridCmdOptionPayload(argv,argv+argc,"--Trajectories");
+      std::vector<int> ivec(0);
+      GridCmdOptionIntVector(arg,ivec);
+      NumTraj = ivec[0];
+    }
+
+    // Create integrator
+    typedef MinimumNorm2<LatticeGaugeField>  IntegratorType;// change here to change the algorithm
+    IntegratorParameters MDpar(20);
+    IntegratorType MDynamics(UGrid,MDpar, TheAction);
+
+    // Checkpoint strategy
+    NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
+    PlaquetteLogger<LatticeGaugeField> PlaqLog(std::string("plaq"));
+
+    HMCparameters HMCpar;
+    HMCpar.StartTrajectory = StartTraj;
+    HMCpar.Trajectories    = NumTraj;
+    
+    GridSerialRNG    sRNG;
+    GridParallelRNG  pRNG(UGrid);
+    LatticeGaugeField  U(UGrid);
+
+    std::vector<int> SerSeed({1,2,3,4,5});
+    std::vector<int> ParSeed({6,7,8,9,10});
+
+    if ( StartType == HotStart ) {
+      // Hot start
+      HMCpar.NoMetropolisUntil =0;
+      HMCpar.MetropolisTest = true;
+      sRNG.SeedFixedIntegers(SerSeed);
+      pRNG.SeedFixedIntegers(ParSeed);
+      SU3::HotConfiguration(pRNG, U);
+    } else if ( StartType == ColdStart ) { 
+      // Cold start
+      HMCpar.NoMetropolisUntil =0;
+      HMCpar.MetropolisTest = true;
+      sRNG.SeedFixedIntegers(SerSeed);
+      pRNG.SeedFixedIntegers(ParSeed);
+      SU3::ColdConfiguration(pRNG, U);
+    } else if ( StartType == TepidStart ) {       
+      // Tepid start
+      HMCpar.NoMetropolisUntil =0;
+      HMCpar.MetropolisTest = true;
+      sRNG.SeedFixedIntegers(SerSeed);
+      pRNG.SeedFixedIntegers(ParSeed);
+      SU3::TepidConfiguration(pRNG, U);
+    } else if ( StartType == CheckpointStart ) { 
+      HMCpar.NoMetropolisUntil =0;
+      HMCpar.MetropolisTest = true;
+      // CheckpointRestart
+      Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
+    }
+
+    HybridMonteCarlo<LatticeGaugeField,IntegratorType>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
+    HMC.AddObservable(&Checkpoint);
+    HMC.AddObservable(&PlaqLog);
+    
+    // Run it
+    HMC.evolve();
+    
+  }
+  
+};
+
+}}
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  const int Ls = 8;
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
 
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  NerscHmcRunner TheHMC;
-  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-
-  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-
-  GridSerialRNG    sRNG;
-  GridParallelRNG  pRNG(UGrid);
-  sRNG.SeedRandomDevice();
-  pRNG.SeedRandomDevice();
-
-  LatticeLorentzColourMatrix     U(UGrid);
-
-  SU3::HotConfiguration(pRNG, U);
-
-  // simplify template declaration? Strip the lorentz from the second template
-  WilsonGaugeActionR Waction(5.6);
-
-  Real mass=0.04;
-  Real pv  =1.0;
-  RealD M5=1.5;
-
-  DomainWallFermionR DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  DomainWallFermionR NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5);
   
-  ConjugateGradient<LatticeFermion>  CG(1.0e-8,10000);
+  TheHMC.BuildTheAction(argc,argv);
-  TwoFlavourEvenOddRatioPseudoFermionAction<WilsonImplR> Nf2(NumOp, DenOp,CG,CG);
-  
-  //Collect actions
-  ActionLevel<LatticeGaugeField> Level1;
-  Level1.push_back(&Nf2);
-  Level1.push_back(&Waction);
-  ActionSet<LatticeGaugeField> FullSet;
-  FullSet.push_back(Level1);
-
-  // Create integrator
-  typedef MinimumNorm2<LatticeGaugeField>  IntegratorType;// change here to change the algorithm
-  IntegratorParameters MDpar(20);
-  IntegratorType MDynamics(UGrid,MDpar, FullSet);
-
-  // Create HMC
-  NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
-  HMCparameters HMCpar;
-  HybridMonteCarlo<LatticeGaugeField,IntegratorType>  HMC(HMCpar, MDynamics,sRNG,pRNG,U);
-  HMC.AddObservable(&Checkpoint);
-
-  // Run it
-  HMC.evolve();
 
 }