Completed testing of the HMC for Ls vectorised version (on AVX2)

lib/cartesian/Cartesian_base.h

@@ -8,6 +8,7 @@
 
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Guido Cossu <guido.cossu@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

lib/lattice/Lattice_rng.h

@@ -148,7 +148,7 @@ namespace Grid {
       ss<<_generators[gen];
       ss.seekg(0,ss.beg);
       for(int i=0;i<RngStateCount;i++){
-      	ss>>saved[i];
+        ss>>saved[i];
       }
     }
     void SetState(std::vector<RngStateType> & saved,int gen){
@@ -315,26 +315,26 @@ namespace Grid {
 
       for(int gidx=0;gidx<gsites;gidx++){
 
-	int rank,o_idx,i_idx;
-	_grid->GlobalIndexToGlobalCoor(gidx,gcoor);
-	_grid->GlobalCoorToRankIndex(rank,o_idx,i_idx,gcoor);
+        int rank,o_idx,i_idx;
+        _grid->GlobalIndexToGlobalCoor(gidx,gcoor);
+        _grid->GlobalCoorToRankIndex(rank,o_idx,i_idx,gcoor);
         
-	int l_idx=generator_idx(o_idx,i_idx);
+        int l_idx=generator_idx(o_idx,i_idx);
 
-	const int num_rand_seed=16;
-	std::vector<int> site_seeds(num_rand_seed);
-	for(int i=0;i<site_seeds.size();i++){
-	  site_seeds[i]= ui(pseeder);
-	}
+        const int num_rand_seed=16;
+        std::vector<int> site_seeds(num_rand_seed);
+        for(int i=0;i<site_seeds.size();i++){
+          site_seeds[i]= ui(pseeder);
+        }
 
 
-      	_grid->Broadcast(0,(void *)&site_seeds[0],sizeof(int)*site_seeds.size());
+        _grid->Broadcast(0,(void *)&site_seeds[0],sizeof(int)*site_seeds.size());
 
-      	if( rank == _grid->ThisRank() ){
+        if( rank == _grid->ThisRank() ){
           fixedSeed ssrc(site_seeds);
           typename source::result_type sinit = ssrc();
           _generators[l_idx] = RngEngine(sinit);
-      	}
+        }
       }
       _seeded=1;
     }    
@@ -373,7 +373,7 @@ namespace Grid {
             for (int idx = 0; idx < words; idx++) 
               fillScalar(pointer[idx], dist[gdx], _generators[gdx]);
           }
-          // merge into SIMD lanes
+          // merge into SIMD lanes, FIXME suboptimal implementation
           merge(l._odata[sm], buf);
         }
       }

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -304,6 +304,8 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
     
     // FIXME 
     // Current implementation works, thread safe, probably suboptimal
+    // Passing through the local coordinate for grid transformation
+    // the force grid is in general very different from the Ls vectorized grid
 
     PARALLEL_FOR_LOOP
     for (int so = 0; so < grid->oSites(); so++) {
@@ -314,9 +316,6 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
         std::vector<int> local_coor;      
         std::vector<int> icoor; grid->iCoorFromIindex(icoor,si);
         grid->InOutCoorToLocalCoor(ocoor, icoor, local_coor);
-        //for (int i = 0; i < dimU; i++) local_coor[i] = ocoor[i] + grid->_rdimensions[i]*icoor[i];
-        //std::cout << "so: " << so << "  si: "<< si << "  local_coor: " << local_coor << std::endl;
-
         for (int s = 0; s < LLs; s++) {
           std::vector<int> slocal_coor(dimF);
           slocal_coor[0] = s;

tests/hmc/Test_hmc_EODWFRatioLsVectorised_Binary.cc

@@ -0,0 +1,170 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+
+Copyright (C) 2015
+
+Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
+Author: neo <cossu@post.kek.jp>
+Author: Guido Cossu <guido.cossu@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;
+
+namespace Grid {
+namespace QCD {
+
+  class HMCRunnerParameters : Serializable {
+  public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(HMCRunnerParameters,
+      double, beta,
+      double, mass,
+      int, MaxCGIterations,
+      double, StoppingCondition,
+      bool, smearedAction,
+      int, SaveInterval,
+      std::string, format,
+      std::string, conf_prefix,
+      std::string, rng_prefix,
+      double, rho,
+      int, SmearingLevels,
+      );
+
+    HMCRunnerParameters() {}
+  };
+
+// Derive from the BinaryHmcRunner (templated for gauge fields)
+class HmcRunner : public BinaryHmcRunner {
+public:
+  void BuildTheAction(int argc, char **argv)
+
+  {
+    typedef DomainWallVec5dImplR ImplPolicy;
+    typedef ScaledShamirFermion<ImplPolicy> FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
+    const int Ls = 8;
+
+    UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+
+    GridCartesian* sUGrid   = SpaceTimeGrid::makeFourDimDWFGrid(GridDefaultLatt(),GridDefaultMpi());
+    GridRedBlackCartesian* sUrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(sUGrid);
+
+    
+    FGrid   = SpaceTimeGrid::makeFiveDimDWFGrid(Ls,UGrid);
+    FrbGrid = SpaceTimeGrid::makeFiveDimDWFRedBlackGrid(Ls,UGrid);
+
+    // temporarily need a gauge field
+    LatticeGaugeField  U(UGrid);
+
+    // Gauge action
+    double beta = 4.0;
+    WilsonGaugeActionR Waction(beta);
+
+    Real mass = 0.04;
+    Real pv   = 1.0;
+    RealD M5  = 1.5;
+    RealD scale = 2.0;
+    FermionAction DenOp(U,*FGrid,*FrbGrid,*sUGrid,*sUrbGrid,mass,M5,scale);
+    FermionAction NumOp(U,*FGrid,*FrbGrid,*sUGrid,*sUrbGrid,pv,M5,scale);
+
+    double StoppingCondition = 1.0e-8;
+    double MaxCGIterations = 10000;
+    ConjugateGradient<FermionField>  CG(StoppingCondition,MaxCGIterations);
+    TwoFlavourEvenOddRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
+
+    // Set smearing (true/false), default: false
+    Nf2.is_smeared = true;
+
+    // Collect actions
+    // here an example of 2 level integration
+    ActionLevel<Field> Level1(1);
+    Level1.push_back(&Nf2);
+
+    // this level will integrate with a
+    // step that is 4 times finer
+    // than the previous level
+    ActionLevel<Field> Level2(4);
+    Level2.push_back(&Waction);
+
+    TheAction.push_back(Level1);
+    TheAction.push_back(Level2);
+
+    // Add observables
+    int SaveInterval = 1;
+    std::string format = std::string("IEEE64BIG");
+    std::string conf_prefix = std::string("DWF_ckpoint_lat");
+    std::string rng_prefix = std::string("DWF_ckpoint_rng");
+    BinaryHmcCheckpointer<BinaryHmcRunner::ImplPolicy> Checkpoint(conf_prefix, rng_prefix, SaveInterval, format);
+    // Can implement also a specific function in the hmcrunner
+    // AddCheckpoint (...) that takes the same parameters + a string/tag
+    // defining the type of the checkpointer
+    // with tags can be implemented by overloading and no ifs
+    // Then force all checkpoint to have few common functions
+    // return an object that is then passed to the Run function
+
+    PlaquetteLogger<BinaryHmcRunner::ImplPolicy> PlaqLog(std::string("Plaquette"));
+    ObservablesList.push_back(&PlaqLog);
+    ObservablesList.push_back(&Checkpoint);
+
+    // Smearing section, omit if not needed
+    double rho = 0.1;  // smearing parameter
+    int Nsmear = 2;    // number of smearing levels
+    Smear_Stout<BinaryHmcRunner::ImplPolicy> Stout(rho);
+    SmearedConfiguration<BinaryHmcRunner::ImplPolicy> SmearingPolicy(
+        UGrid, Nsmear, Stout);
+    ///////////////////
+
+    Run(argc, argv, Checkpoint, SmearingPolicy); 
+    //Run(argc, argv, Checkpoint);  // no smearing
+};
+};
+}
+}
+
+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;
+
+  HmcRunner TheHMC;
+
+
+  // Seeds for the random number generators
+  std::vector<int> SerSeed({1, 2, 3, 4, 5});
+  std::vector<int> ParSeed({6, 7, 8, 9, 5});
+  TheHMC.RNGSeeds(SerSeed, ParSeed);
+
+  TheHMC.MDparameters.set(20, 1.0);// MDsteps, traj length
+
+  TheHMC.BuildTheAction(argc, argv);
+
+  Grid_finalize();
+}

tests/hmc/Test_hmc_EODWFRatio_Binary.cc

@@ -62,8 +62,7 @@ public:
   void BuildTheAction(int argc, char **argv)
 
   {
-    //typedef WilsonImplR ImplPolicy;
-    typedef DomainWallVec5dImplR ImplPolicy;
+    typedef WilsonImplR ImplPolicy;
     typedef ScaledShamirFermion<ImplPolicy> FermionAction;
     typedef typename FermionAction::FermionField FermionField;
 
@@ -71,37 +70,23 @@ public:
 
     UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
     UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-
-
-    GridCartesian* sUGrid   = SpaceTimeGrid::makeFourDimDWFGrid(GridDefaultLatt(),GridDefaultMpi());
-    GridRedBlackCartesian* sUrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(sUGrid);
-
     
-    FGrid   = SpaceTimeGrid::makeFiveDimDWFGrid(Ls,UGrid);
-    FrbGrid = SpaceTimeGrid::makeFiveDimDWFRedBlackGrid(Ls,UGrid);
-    
-    /*
     FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
     FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-    */
-
 
     // temporarily need a gauge field
     LatticeGaugeField  U(UGrid);
 
     // Gauge action
     double beta = 4.0;
-    WilsonGaugeActionR Iaction(beta);
+    WilsonGaugeActionR Waction(beta);
 
     Real mass = 0.04;
     Real pv   = 1.0;
     RealD M5  = 1.5;
     RealD scale = 2.0;
-    FermionAction DenOp(U,*FGrid,*FrbGrid,*sUGrid,*sUrbGrid,mass,M5,scale);
-    FermionAction NumOp(U,*FGrid,*FrbGrid,*sUGrid,*sUrbGrid,pv,M5,scale);
-
-    std::cout << GridLogMessage << "Frb Osites: " << FrbGrid->oSites() << std::endl;
-    std::cout << GridLogMessage << "sUGrid Osites: " << sUGrid->oSites() << std::endl;
+    FermionAction DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,scale);
+    FermionAction NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5,scale);
 
     double StoppingCondition = 1.0e-8;
     double MaxCGIterations = 10000;
@@ -109,7 +94,7 @@ public:
     TwoFlavourEvenOddRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
 
     // Set smearing (true/false), default: false
-    Nf2.is_smeared = false;
+    Nf2.is_smeared = true;
 
     // Collect actions
     // here an example of 2 level integration
@@ -120,7 +105,7 @@ public:
     // step that is 4 times finer
     // than the previous level
     ActionLevel<Field> Level2(4);
-    Level2.push_back(&Iaction);
+    Level2.push_back(&Waction);
 
     TheAction.push_back(Level1);
     TheAction.push_back(Level2);