Cleaning up the CG reproduciblity test. More info reported

lib/Bitwise.cc

@@ -0,0 +1,48 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/Bitwise.cc
+
+Copyright (C) 2016
+
+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 */
+#ifndef GRID_BITWISE_H
+#define GRID_BITWISE_H
+
+#include <iostream>
+#include <Bitwise.h>
+#include <bitset>
+#include <climits>
+
+namespace Grid {
+
+void show_binaryrep(const unsigned char* a, size_t size) {
+  const unsigned char* beg = a;
+  const unsigned char* end = a + size;
+  while (beg != end) std::cout << std::bitset<CHAR_BIT>(*beg++) << ' ';
+  std::cout << '\n';
+}
+
+} // namespace 
+
+#endif

lib/Bitwise.h

@@ -0,0 +1,59 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/Bitwise.h
+
+Copyright (C) 2016
+
+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 */
+#ifndef GRID_BITWISE_H
+#define GRID_BITWISE_H
+
+#include <bitset>
+#include <climits>
+
+namespace Grid {
+
+void show_binaryrep(const unsigned char* a, size_t size);
+
+template <typename T>
+void show_binaryrep(const T& a) {
+  const char* beg = reinterpret_cast<const char*>(&a);
+  const char* end = beg + sizeof(a);
+  while (beg != end) std::cout << std::bitset<CHAR_BIT>(*beg++) << ' ';
+  std::cout << '\n';
+}
+
+template <typename T>
+void bitwise_xor(T& l, T& r, unsigned char* xors) {
+  assert(sizeof(l) == sizeof(r));
+  unsigned char* org = reinterpret_cast<unsigned char*>(&l);
+  unsigned char* cur = reinterpret_cast<unsigned char*>(&r);
+  int words = sizeof(l) / sizeof(*org);
+  unsigned char result = 0;
+  for (int w = 0; w < words; w++) xors[w] = (org[w] ^ cur[w]);
+}
+
+}; // namespace 
+
+#endif

lib/Grid.h

@@ -62,6 +62,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/serialisation/Serialisation.h>
 #include "Config.h"
 #include <Grid/Timer.h>
+#include <Grid/Bitwise.h>
 #include <Grid/PerfCount.h>
 #include <Grid/Log.h>
 #include <Grid/AlignedAllocator.h>

lib/algorithms/iterative/ConjugateGradient.h

@@ -37,10 +37,7 @@ struct CG_state {
   bool do_repro;
   std::vector<RealD> residuals;
 
-  CG_state() {
+  CG_state() {reset();}
-    do_repro = false;
-    residuals.clear();
-  }
 
   void reset(){
     do_repro = false;
@@ -71,7 +68,14 @@ class ConjugateGradient : public OperatorFunction<Field> {
       : Tolerance(tol),
         MaxIterations(maxit),
         ErrorOnNoConverge(err_on_no_conv),
-        ReproTest(ReproducibilityTest){};
+        ReproTest(ReproducibilityTest){
+        	if(ReproducibilityTest == true && err_on_no_conv == true){
+        		std::cout << GridLogMessage << "CG: Reproducibility test ON "<<
+        					"and error on convergence ON are incompatible options" << std::endl;
+        	exit(1);
+        	}
+        	
+        };
 
 
   void operator()(LinearOperatorBase<Field> &Linop, const Field &src,
@@ -210,6 +214,7 @@ class ConjugateGradient : public OperatorFunction<Field> {
                 ReprTest.do_check = true;
                 ReprTest.reset_counter();
                 this->operator()(Linop, src, psi_start);// run the repro test
+                std::cout << GridLogMessage << "Test passed" << std::endl;
         }
 
         // Clear state

lib/communicator/Communicator_base.cc

@@ -65,6 +65,7 @@ const std::vector<int> & CartesianCommunicator::ThisProcessorCoor(void) { return
 const std::vector<int> & CartesianCommunicator::ProcessorGrid(void)     { return _processors; };
 int                      CartesianCommunicator::ProcessorCount(void)    { return _Nprocessors; };
 
+
 ////////////////////////////////////////////////////////////////////////////////
 // very VERY rarely (Log, serial RNG) we need world without a grid
 ////////////////////////////////////////////////////////////////////////////////
@@ -89,11 +90,11 @@ void CartesianCommunicator::GlobalSumVector(ComplexD *c,int N)
 #if !defined( GRID_COMMS_MPI3) && !defined (GRID_COMMS_MPI3L)
 
 void CartesianCommunicator::StencilSendToRecvFromBegin(std::vector<CommsRequest_t> &list,
-						       void *xmit,
+                                                       void *xmit,
-						       int xmit_to_rank,
+                                                       int xmit_to_rank,
-						       void *recv,
+                                                       void *recv,
-						       int recv_from_rank,
+                                                       int recv_from_rank,
-						       int bytes)
+                                                       int bytes)
 {
   SendToRecvFromBegin(list,xmit,xmit_to_rank,recv,recv_from_rank,bytes);
 }

lib/communicator/Communicator_base.h

@@ -68,6 +68,8 @@ class CartesianCommunicator {
   static MPI_Comm communicator_world;
          MPI_Comm communicator;
   typedef MPI_Request CommsRequest_t;
+  static char name[MPI_MAX_PROCESSOR_NAME]; // processing node physical name
+  static int length; 
 #else 
   typedef int CommsRequest_t;
 #endif
@@ -149,6 +151,7 @@ class CartesianCommunicator {
   const std::vector<int> & ProcessorGrid(void)     ;
   int                      ProcessorCount(void)    ;
 
+  void        		 PrintRankInfo(void)     ;
   ////////////////////////////////////////////////////////////////////////////////
   // very VERY rarely (Log, serial RNG) we need world without a grid
   ////////////////////////////////////////////////////////////////////////////////

lib/communicator/Communicator_mpi.cc

@@ -35,6 +35,8 @@ namespace Grid {
 // Info that is setup once and indept of cartesian layout
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 MPI_Comm CartesianCommunicator::communicator_world;
+char CartesianCommunicator::name[MPI_MAX_PROCESSOR_NAME]; // processing node physical name
+int CartesianCommunicator::length; 
 
 // Should error check all MPI calls.
 void CartesianCommunicator::Init(int *argc, char ***argv) {
@@ -44,6 +46,8 @@ void CartesianCommunicator::Init(int *argc, char ***argv) {
     MPI_Init(argc,argv);
   }
   MPI_Comm_dup (MPI_COMM_WORLD,&communicator_world);
+
+  MPI_Get_processor_name(name, &length);
   ShmInitGeneric();
 }
 
@@ -206,5 +210,10 @@ void CartesianCommunicator::BroadcastWorld(int root,void* data, int bytes)
   assert(ierr==0);
 }
 
+void CartesianCommunicator::PrintRankInfo(){
+    std::cout << "Grid: Rank "<< _processor << "  -  Physical node name: " << name << std::endl;
 }
 
+
+}// end of namespace
+

lib/communicator/Communicator_mpi3.cc

@@ -576,5 +576,10 @@ void CartesianCommunicator::BroadcastWorld(int root,void* data, int bytes)
   assert(ierr==0);
 }
 
+  void CartesianCommunicator::PrintRankInfo(){
+    std::cout << "Grid: Rank "<< _processor << "  -  Physical node name: " << name << std::endl;
+  }
+
+  
 }
 

lib/communicator/Communicator_mpi3_leader.cc

@@ -869,6 +869,9 @@ void *CartesianCommunicator::ShmBufferTranslate(int rank,void * local_p) {
   return NULL;
 }
 
+  void CartesianCommunicator::PrintRankInfo(){
+    std::cout << "Grid: Rank "<< _processor << "  -  Physical node name: " << name << std::endl;
+  }
 
 };
 

lib/communicator/Communicator_none.cc

@@ -37,6 +37,11 @@ void CartesianCommunicator::Init(int *argc, char *** arv)
   ShmInitGeneric();
 }
 
+void CartesianCommunicator::PrintRankInfo(){
+    std::cout << GridLogMessage << "No Rank Info available" << std::endl;
+}
+
+
 CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors)
 {
   _processors = processors;
@@ -60,10 +65,10 @@ void CartesianCommunicator::GlobalSum(uint64_t &){}
 void CartesianCommunicator::GlobalSumVector(double *,int N){}
 
 void CartesianCommunicator::SendRecvPacket(void *xmit,
-					   void *recv,
+                                           void *recv,
-					   int xmit_to_rank,
+                                           int xmit_to_rank,
-					   int recv_from_rank,
+                                           int recv_from_rank,
-					   int bytes)
+                                           int bytes)
 {
   assert(0);
 }
@@ -71,19 +76,19 @@ void CartesianCommunicator::SendRecvPacket(void *xmit,
 
 // Basic Halo comms primitive -- should never call in single node
 void CartesianCommunicator::SendToRecvFrom(void *xmit,
-					   int dest,
+                                           int dest,
-					   void *recv,
+                                           void *recv,
-					   int from,
+                                           int from,
-					   int bytes)
+                                           int bytes)
 {
   assert(0);
 }
 void CartesianCommunicator::SendToRecvFromBegin(std::vector<CommsRequest_t> &list,
-						void *xmit,
+                                                void *xmit,
-						int dest,
+                                                int dest,
-						void *recv,
+                                                void *recv,
-						int from,
+                                                int from,
-						int bytes)
+                                                int bytes)
 {
   assert(0);
 }

lib/communicator/Communicator_shmem.cc

@@ -333,5 +333,9 @@ void CartesianCommunicator::BroadcastWorld(int root,void* data, int bytes)
   }
 }
 
+  void CartesianCommunicator::PrintRankInfo(){
+    std::cout << GridLogMessage << "SHMEM: Rank Info not implemented yet" << std::endl;
+  }
+
 }
 

lib/lattice/Lattice_reduction.h

@@ -33,26 +33,28 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 namespace Grid {
 
 template <class T>
-struct ReproducibilityState {
+class ReproducibilityState {
+ public:
   typedef typename T::vector_type vector_type;
+  typedef std::vector<std::vector<vector_type, alignedAllocator<vector_type> > > sum_type;
   unsigned int n_call;
   bool do_check;
   bool enable_reprocheck;
-  std::vector<std::vector<vector_type, alignedAllocator<vector_type> > >
+  sum_type th_states;
-      th_states;
 
+  void reset_counter() { n_call = 0; }
   void reset() {
     th_states.clear();
     do_check = false;
     enable_reprocheck = false;
     n_call = 0;
-  }
+  };
-
-  void reset_counter() { n_call = 0; }
 
   ReproducibilityState() { reset(); }
 };
 
+
+
 #ifdef GRID_WARN_SUBOPTIMAL
 #warning "Optimisation alert all these reduction loops are NOT threaded "
 #endif     
@@ -94,8 +96,6 @@ struct ReproducibilityState {
         sumarray[i]=zero;
       }
 
-      // accumulation done in the same precision ad vobj...
-      // may need to froce higher precision
       PARALLEL_FOR_LOOP_STATIC //request statically scheduled threads for reproducibility
       for(int thr=0;thr<grid->SumArraySize();thr++){
         int nwork, mywork, myoff;
@@ -103,30 +103,41 @@ struct ReproducibilityState {
         
         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]);// accumulate here in higher precision
+          vnrm = vnrm + innerProduct(left._odata[ss],right._odata[ss]);
         }
         sumarray[thr]=TensorRemove(vnrm) ;
       }
       
-
       ///////////////////////  Reproducibility section
       if (repr.enable_reprocheck) {
         if (repr.do_check) {
-          //std::cout << GridLogMessage << "Checking thread state for inner product. Call n. " << repr.n_call << std::endl;
+          std::cout << GridLogDebug << "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]) {
+            int words = sizeof(sumarray[thread])/sizeof(unsigned char);
-              std::cout << GridLogMessage << "Reproducibility failure on node " << grid->ThisRank() << std::endl;
+            unsigned char xors[words];
+            bitwise_xor(sumarray[thread], repr.th_states[repr.n_call][thread],xors);
+            // XOR all words
+            unsigned char res = 0;
+            for (int w = 0; w < words; w++) res = res ^ xors[w];
+            if ( res ) {
+              std::cout << GridLogMessage << "Reproducibility failure report" << std::endl;
+              grid->PrintRankInfo();
               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 << "Current partial sum: " << sumarray[thread] << std::endl;
-              std::cout << GridLogMessage << repr.th_states[repr.n_call][thread] << std::endl;
+              std::cout << GridLogMessage << "Saved partial sum  : " << repr.th_states[repr.n_call][thread] << std::endl;
-              //exit(1);
+              std::cout << GridLogMessage << "Saved state " << std::endl;
+              show_binaryrep(repr.th_states[repr.n_call][thread]);
+              std::cout << GridLogMessage << "Current state" << std::endl;
+              show_binaryrep(sumarray[thread]);
+              std::cout << GridLogMessage << "Xor result" << std::endl;
+              show_binaryrep(xors, words);                          
             }
           }
           repr.n_call++;
         } else 
         {
-          //std::cout << GridLogMessage << "Saving thread state for inner product. Call n. " << repr.n_call << std::endl;
+          std::cout << GridLogDebug << "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

tests/hmc/Test_hmc_WilsonFermionGauge.cc

@@ -64,9 +64,9 @@ class HmcRunner : public NerscHmcRunner {
     FermionAction FermOp(U, *FGrid, *FrbGrid, mass);
 
     // To enable the CG reproducibility tests use 
-    // ConjugateGradient<FermionField> CG(1.0e-8, 10000, true, true);
+    ConjugateGradient<FermionField> CG(1.0e-8, 10000, false, true);
     // This is the plain version
-    ConjugateGradient<FermionField> CG(1.0e-8, 10000);
+    //ConjugateGradient<FermionField> CG(1.0e-8, 10000);
 
     TwoFlavourPseudoFermionAction<ImplPolicy> Nf2(FermOp, CG, CG);