Reproducibility checks for inner product

.gitignore

@@ -49,6 +49,7 @@ config.status
 .deps
 Make.inc
 eigen.inc
+Eigen.inc
 
 # http://www.gnu.org/software/autoconf #
 ########################################

lib/Threads.h

@@ -46,6 +46,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #endif
 #define PARALLEL_NESTED_LOOP2    _Pragma("omp parallel for collapse(2)")
 #define PARALLEL_REGION          _Pragma("omp parallel")
+#define PARALLEL_FOR_LOOP_STATIC _Pragma("omp parallel for schedule(static)")
 #else
 #define PARALLEL_FOR_LOOP
 #define PARALLEL_FOR_LOOP_INTERN

lib/algorithms/iterative/ConjugateGradient.h

@@ -62,6 +62,7 @@ class ConjugateGradient : public OperatorFunction<Field> {
   Integer MaxIterations;
   bool ReproTest;
   CG_state CGState;//to check reproducibility by repeating the CG
+  ReproducibilityState<typename Field::vector_object> ReprTest;
 
   ConjugateGradient(RealD tol, Integer maxit, bool err_on_no_conv = true,
         bool ReproducibilityTest = false)
@@ -93,6 +94,10 @@ class ConjugateGradient : public OperatorFunction<Field> {
     
     Linop.HermOpAndNorm(psi, mmp, d, b);
     
+    if(!ReprTest.do_check)
+    	ReprTest.reset();
+    ReprTest.enable_reprocheck=ReproTest;
+
 
     r = src - mmp;
     p = r;
@@ -146,7 +151,8 @@ class ConjugateGradient : public OperatorFunction<Field> {
       b_pred = a * (a * qq - d) / c;// a check
 
 
-      cp = axpy_norm(r, -a, mmp, r);// new residual r = r_old - a * Ap
+      axpy(r, -a, mmp, r);// new residual r = r_old - a * Ap
+      cp = norm2(r, ReprTest);// 
       if (ReproTest && !CGState.do_repro) {
         CGState.residuals.push_back(cp);  // save residuals state
                 std::cout << GridLogIterative << "ReproTest: Saving state" << std::endl;
@@ -199,11 +205,14 @@ class ConjugateGradient : public OperatorFunction<Field> {
 
         if (!CGState.do_repro && ReproTest){
                 CGState.do_repro = true;
+                ReprTest.do_check = true;
+                ReprTest.reset_counter();
                 this->operator()(Linop, src, psi_start);// run the repro test
         }
 
         // Clear state
         CGState.reset();
+        ReprTest.reset();
         return;
       }
     }

lib/lattice/Lattice_reduction.h

@@ -34,21 +34,24 @@ namespace Grid {
 
 template <class T>
 struct ReproducibilityState {
-  int n_call;
+  typedef typename T::vector_type vector_type;
+  unsigned int n_call;
   bool do_check;
-  std::vector<std::vector<T, alignedAllocator<T> > > th_states;
+  bool enable_reprocheck;
+  std::vector<std::vector<vector_type, alignedAllocator<vector_type> > >
+      th_states;
 
   void reset() {
     th_states.clear();
     do_check = false;
+    enable_reprocheck = false;
     n_call = 0;
   }
 
-  ReproducibilityState(){
+  void reset_counter() { n_call = 0; }
-    reset();
-  }
-};
 
+  ReproducibilityState() { reset(); }
+};
 
 #ifdef GRID_WARN_SUBOPTIMAL
 #warning "Optimisation alert all these reduction loops are NOT threaded "
@@ -58,15 +61,27 @@ struct ReproducibilityState {
     // Deterministic Reduction operations
     ////////////////////////////////////////////////////////////////////////////////////////////////////
   template<class vobj> inline RealD norm2(const Lattice<vobj> &arg){
-      ComplexD nrm = innerProduct(arg,arg);
+        ReproducibilityState<vobj> repr;
+        ComplexD nrm = innerProduct(arg, arg, repr);
         return std::real(nrm); 
     }
 
+    template <class vobj>
+    inline RealD norm2(const Lattice<vobj> &arg, ReproducibilityState<vobj>& rep) {
+      ComplexD nrm = innerProduct(arg, arg, rep);
+      return std::real(nrm);
+    }
+
+    template<class vobj>
+    inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right){
+      ReproducibilityState<vobj> repr;
+      return innerProduct(left, right, repr);
+    } 
+    
 
 
-
     template<class vobj>
-    inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right) 
+    inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right, ReproducibilityState<vobj>& repr) 
     {
       typedef typename vobj::scalar_type scalar_type;
       typedef typename vobj::vector_type vector_type;
@@ -81,30 +96,45 @@ struct ReproducibilityState {
 
       // accumulation done in the same precision ad vobj...
       // may need to froce higher precision
-      PARALLEL_FOR_LOOP
+      PARALLEL_FOR_LOOP_STATIC //request statically scheduled threads for reproducibility
       for(int thr=0;thr<grid->SumArraySize();thr++){
         int nwork, mywork, myoff;
         GridThread::GetWork(left._grid->oSites(),thr,mywork,myoff);
         
         decltype(innerProduct(left._odata[0],right._odata[0])) vnrm = zero; // private to thread; sub summation
         for(int ss=myoff;ss<mywork+myoff; ss++){
-          vnrm = vnrm + innerProduct(left._odata[ss],right._odata[ss]);
+          vnrm = vnrm + innerProduct(left._odata[ss],right._odata[ss]);// accumulate here in higher precision
         }
         sumarray[thr]=TensorRemove(vnrm) ;
       }
       
-/*
+
+      ///////////////////////  Reproducibility section
+      if (repr.enable_reprocheck) {
         if (repr.do_check) {
-        if (sumarray!=repr.th_states[n_call]){
+          //std::cout << GridLogMessage << "Checking thread state for inner product. Call n. " << repr.n_call << std::endl;
+          for (int thread = 0; thread < sumarray.size(); thread++) {
+            if (sumarray[thread] != repr.th_states[repr.n_call][thread]) {
               std::cout << GridLogMessage << "Reproducibility failure on node " << grid->ThisRank() << std::endl;
-          exit(1);
+              std::cout << GridLogMessage << "Call: "<< repr.n_call << " Thread: " << thread << std::endl;
+              std::cout << GridLogMessage << "Size of states: " << repr.th_states.size() << std::endl;
+              std::cout << GridLogMessage << sumarray[thread] << std::endl;
+              std::cout << GridLogMessage << repr.th_states[repr.n_call][thread] << std::endl;
+              //exit(1);
             }
           }
-      } else {
+          repr.n_call++;
-        repr.th_states.push_back(sumarray);//save threads state
+        } else 
-        repr.n_call +=1;
+        {
+          //std::cout << GridLogMessage << "Saving thread state for inner product. Call n. " << repr.n_call << std::endl;
+          repr.th_states.resize(repr.n_call+1);
+          repr.th_states[repr.n_call].resize(grid->SumArraySize());
+          repr.th_states[repr.n_call] = sumarray;  // save threads state
+          repr.n_call++;
         }
-      */
+      }
+      ////////////////////////////////////////////////////////
+
 
       vector_type vvnrm; vvnrm=zero;  // sum across threads
       for(int i=0;i<grid->SumArraySize();i++){

lib/simd/Grid_avx.h

@@ -66,6 +66,9 @@ namespace Optimization {
     double f[4];
   };
 
+
+
+
   struct Vsplat{
     //Complex float
     inline __m256 operator()(float a, float b){

lib/simd/Grid_vector_types.h

@@ -146,6 +146,23 @@ class Grid_simd {
   Grid_simd(const Grid_simd &rhs) : v(rhs.v){};  // compiles in movaps
   Grid_simd(const Grid_simd &&rhs) : v(rhs.v){};
 
+  /////////////////////////////
+  // Comparisons
+  /////////////////////////////
+  inline bool operator==(const Grid_simd& lhs){ 
+    Grid_simd::conv_t conv1;
+    Grid_simd::conv_t conv2;
+    conv1.v = v;
+    conv2.v = lhs.v;
+    return std::equal(std::begin(conv1.s), std::end(conv1.s), std::begin(conv2.s));
+  }
+
+  inline bool operator!=(const Grid_simd& lhs){ 
+    return !(*this == lhs);
+  }
+
+
+
   /////////////////////////////
   // Constructors
   /////////////////////////////