Got unpreconditioned conjugate gradient to run and converge on a random (uniform random,

not even SU(3) for now) gauge field. Convergence history is correctly indepdendent of decomposition
on 1,2,4,8,16 mpi tasks.
Found a couple of simd bugs which required fixed and enhanced the Grid_simd.cc test suite.
Implemented the Mdag, M, MdagM, Meooe Mooee schur type stuff in the wilson dop.

tests/Grid_main.cc

@@ -140,7 +140,7 @@ int main (int argc, char ** argv)
 
 //     -=,+=,*=,()
 //     add,+,sub,-,mult,mac,*
-//     adj,conj
+//     adj,conjugate
 //     real,imag
 //     transpose,transposeIndex  
 //     trace,traceIndex
@@ -437,7 +437,7 @@ int main (int argc, char ** argv)
         for(int r=0;r<3;r++){
         for(int c=0;c<3;c++){
             diff =shifted1()()(r,c)-shifted2()()(r,c);
-            nn=real(conj(diff)*diff);
+            nn=real(conjugate(diff)*diff);
             if ( nn > 0 )
                 cout<<"Shift fail (shifted1/shifted2-ref) "<<coor[0]<<coor[1]<<coor[2]<<coor[3] <<" "
                     <<shifted1()()(r,c)<<" "<<shifted2()()(r,c)
@@ -451,7 +451,7 @@ int main (int argc, char ** argv)
         for(int r=0;r<3;r++){
         for(int c=0;c<3;c++){
             diff =shifted3()()(r,c)-shifted2()()(r,c);
-            nn=real(conj(diff)*diff);
+            nn=real(conjugate(diff)*diff);
             if ( nn > 0 )
                 cout<<"Shift rb fail (shifted3/shifted2-ref) "<<coor[0]<<coor[1]<<coor[2]<<coor[3] <<" "
                 <<shifted3()()(r,c)<<" "<<shifted2()()(r,c)
@@ -468,7 +468,7 @@ int main (int argc, char ** argv)
         for(int r=0;r<Nc;r++){
         for(int c=0;c<Nc;c++){
             diff =foobar2()()(r,c)-foobar1()()(r,c);
-            nrm = nrm + real(conj(diff)*diff);
+            nrm = nrm + real(conjugate(diff)*diff);
         }}
     }}}}
 	if( Fine.IsBoss() ){

tests/Grid_rng_fixed.cc

@@ -0,0 +1,52 @@
+#include <Grid.h>
+#include <parallelIO/GridNerscIO.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(4,vComplexF::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+     
+  GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
+
+  std::vector<int> seeds({1,2,3,4});
+
+  GridSerialRNG            fsRNG;  fsRNG.SeedFixedIntegers(seeds);
+  GridParallelRNG          fpRNG(&Grid);  fpRNG.SeedFixedIntegers(seeds);
+
+  vComplexF tmp; random(fsRNG,tmp);
+  std::cout<<"Random vComplexF (fixed seed)\n"<< tmp<<std::endl;
+  std::cout<<"conjugate(tmp)\n"<< conjugate(tmp)<<std::endl;
+  std::cout<<"conjugate(tmp)*tmp\n"<< conjugate(tmp)*tmp<<std::endl;
+  std::cout<<"innerProduct"<< innerProduct(tmp,tmp)<<std::endl;
+  std::cout<<"Reduce(innerProduct)"<< Reduce(innerProduct(tmp,tmp))<<std::endl;
+
+  SpinMatrix rnd  ; 
+  random(fsRNG,rnd);
+  std::cout<<"Random Spin Matrix (fixed seed)\n"<< rnd<<std::endl;
+
+  SpinVector rv; 
+  random(fsRNG,rv);
+  std::cout<<"Random Spin Vector (fixed seed)\n"<< rv<<std::endl;
+
+  gaussian(fsRNG,rv);
+  std::cout<<"Gaussian Spin Vector (fixed seed)\n"<< rv<<std::endl;
+
+  LatticeColourVector lcv(&Grid);
+
+  LatticeFermion src(&Grid); random(fpRNG,src);
+  std::cout << "src norm : " << norm2(src)<<std::endl;
+  std::cout << "src " << src<<std::endl;
+
+  random(fpRNG,lcv);
+  std::cout<<"Random Lattice Colour Vector (fixed seed)\n"<< lcv<<std::endl;
+
+
+  Grid_finalize();
+}

tests/Grid_simd.cc

@@ -26,7 +26,7 @@ public:
 class funcConj {
 public:
   funcConj() {};
-  template<class vec> void operator()(vec &rr,vec &i1,vec &i2) const { rr = conj(i1);}
+  template<class vec> void operator()(vec &rr,vec &i1,vec &i2) const { rr = conjugate(i1);}
   std::string name(void) const { return std::string("Conj"); }
 };
 class funcAdj {
@@ -49,6 +49,34 @@ public:
   template<class vec> void operator()(vec &rr,vec &i1,vec &i2) const { rr = timesMinusI(i1);}
   std::string name(void) const { return std::string("timesMinusI"); }
 };
+class funcInnerProduct {
+public:
+  funcInnerProduct() {};
+  template<class vec> void operator()(vec &rr,vec &i1,vec &i2) const { rr = innerProduct(i1,i2);}
+  std::string name(void) const { return std::string("innerProduct"); }
+};
+
+// FIXME still to test:
+//
+//  innerProduct,
+//  norm2, 
+//  Reduce,
+//
+//  mac,mult,sub,add, vone,vzero,vcomplex_i, =zero,
+//  vset,vsplat,vstore,vstream,vload, scalar*vec, vec*scalar
+//  unary -,
+//  *= , -=, +=
+//  outerproduct, 
+//  zeroit
+//  permute
+
+class funcReduce {
+public:
+  funcReduce() {};
+template<class reduce,class vec>    void vfunc(reduce &rr,vec &i1,vec &i2)   const { rr = Reduce(i1);}
+template<class reduce,class scal>   void sfunc(reduce &rr,scal &i1,scal &i2) const { rr = i1;}
+  std::string name(void) const { return std::string("Reduce"); }
+};
 
 template<class scal, class vec,class functor > 
 void Tester(const functor &func)
@@ -96,8 +124,59 @@ void Tester(const functor &func)
       ok++;
     }
   }
-  if ( ok==0 ) std::cout << " OK!" <<std::endl;
+  if ( ok==0 ) {
+    std::cout << " OK!" <<std::endl;
+  }
+  assert(ok==0);
+}
 
+
+template<class reduced,class scal, class vec,class functor > 
+void ReductionTester(const functor &func)
+{
+  GridSerialRNG          sRNG;
+  sRNG.SeedRandomDevice();
+  
+  int Nsimd = vec::Nsimd();
+
+  std::vector<scal> input1(Nsimd);
+  std::vector<scal> input2(Nsimd);
+  reduced result(0);
+  reduced reference(0);
+  reduced tmp;
+
+  std::vector<vec,alignedAllocator<vec> > buf(3);
+  vec & v_input1 = buf[0];
+  vec & v_input2 = buf[1];
+
+
+  for(int i=0;i<Nsimd;i++){
+    random(sRNG,input1[i]);
+    random(sRNG,input2[i]);
+  }
+
+  merge<vec,scal>(v_input1,input1);
+  merge<vec,scal>(v_input2,input2);
+
+  func.template vfunc<reduced,vec>(result,v_input1,v_input2);
+
+  for(int i=0;i<Nsimd;i++) {
+    func.template sfunc<reduced,scal>(tmp,input1[i],input2[i]);
+    reference+=tmp;
+  }
+
+  std::cout << " " << func.name()<<std::endl;
+
+  int ok=0;
+  if ( abs(reference-result)/abs(reference) > 1.0e-6 ){ // rounding is possible for reduce order
+    std::cout<< "*****" << std::endl;
+    std::cout<< abs(reference-result) << " " <<reference<< " " << result<<std::endl;
+    ok++;
+  }
+  if ( ok==0 ) {
+    std::cout << " OK!" <<std::endl;
+  }
+  assert(ok==0);
 }
 
 
@@ -127,6 +206,8 @@ int main (int argc, char ** argv)
   Tester<ComplexF,vComplexF>(funcTimes());
   Tester<ComplexF,vComplexF>(funcConj());
   Tester<ComplexF,vComplexF>(funcAdj());
+  Tester<ComplexF,vComplexF>(funcInnerProduct());
+  ReductionTester<ComplexF,ComplexF,vComplexF>(funcReduce());
 
   std::cout << "==================================="<<  std::endl;
   std::cout << "Testing vComplexD "<<std::endl;
@@ -140,6 +221,8 @@ int main (int argc, char ** argv)
   Tester<ComplexD,vComplexD>(funcTimes());
   Tester<ComplexD,vComplexD>(funcConj());
   Tester<ComplexD,vComplexD>(funcAdj());
+  Tester<ComplexD,vComplexD>(funcInnerProduct());
+  ReductionTester<ComplexD,ComplexD,vComplexD>(funcReduce());
 
   std::cout << "==================================="<<  std::endl;
   std::cout << "Testing vRealF "<<std::endl;
@@ -150,6 +233,9 @@ int main (int argc, char ** argv)
   Tester<RealF,vRealF>(funcMinus());
   Tester<RealF,vRealF>(funcTimes());
   Tester<RealF,vRealF>(funcAdj());
+  Tester<RealF,vRealF>(funcConj());
+  Tester<RealF,vRealF>(funcInnerProduct());
+  ReductionTester<RealF,RealF,vRealF>(funcReduce());
 
   std::cout << "==================================="<<  std::endl;
   std::cout << "Testing vRealD "<<std::endl;
@@ -159,6 +245,9 @@ int main (int argc, char ** argv)
   Tester<RealD,vRealD>(funcMinus());
   Tester<RealD,vRealD>(funcTimes());
   Tester<RealD,vRealD>(funcAdj());
+  Tester<RealD,vRealD>(funcConj());
+  Tester<RealD,vRealD>(funcInnerProduct());
+  ReductionTester<RealD,RealD,vRealD>(funcReduce());
 
   Grid_finalize();
 }

tests/Grid_stencil.cc

@@ -93,7 +93,7 @@ int main (int argc, char ** argv)
 	  for(int r=0;r<3;r++){
 	  for(int c=0;c<3;c++){
             diff =check()()(r,c)-bar()()(r,c);
-            double nn=real(conj(diff)*diff);
+            double nn=real(conjugate(diff)*diff);
             if ( nn > 0){
 	      printf("Coor (%d %d %d %d) \t rc %d%d \t %le (%le,%le) %le\n",
 		     coor[0],coor[1],coor[2],coor[3],r,c,
@@ -103,8 +103,8 @@ int main (int argc, char ** argv)
 		     real(bar()()(r,c))
 		     );
 	    }
-	    snrmC=snrmC+real(conj(check()()(r,c))*check()()(r,c));
-	    snrmB=snrmB+real(conj(bar()()(r,c))*bar()()(r,c));
+	    snrmC=snrmC+real(conjugate(check()()(r,c))*check()()(r,c));
+	    snrmB=snrmB+real(conjugate(bar()()(r,c))*bar()()(r,c));
 	    snrm=snrm+nn;
 	  }}
 	 

tests/Makefile.am

@@ -5,7 +5,7 @@ AM_LDFLAGS = -L$(top_builddir)/lib
 #
 # Test code
 #
-bin_PROGRAMS = Grid_main Grid_stencil Grid_nersc_io Grid_cshift Grid_gamma  Grid_simd Grid_rng Grid_remez
+bin_PROGRAMS = Grid_main Grid_stencil Grid_nersc_io Grid_cshift Grid_gamma  Grid_simd Grid_rng Grid_remez Grid_rng_fixed
 
 Grid_main_SOURCES = Grid_main.cc
 Grid_main_LDADD = -lGrid
@@ -13,6 +13,9 @@ Grid_main_LDADD = -lGrid
 Grid_rng_SOURCES = Grid_rng.cc
 Grid_rng_LDADD = -lGrid
 
+Grid_rng_fixed_SOURCES = Grid_rng_fixed.cc
+Grid_rng_fixed_LDADD = -lGrid
+
 Grid_remez_SOURCES = Grid_remez.cc
 Grid_remez_LDADD = -lGrid
 

tests/test_Grid_jacobi.cc

@@ -186,7 +186,7 @@ int main (int argc, char ** argv)
 	  for(int r=0;r<3;r++){
 	  for(int c=0;c<3;c++){
             diff =check._internal._internal[r][c]-bar._internal._internal[r][c];
-            double nn=real(conj(diff)*diff);
+            double nn=real(conjugate(diff)*diff);
             if ( nn > 0 ){
 	      printf("Coor (%d %d %d %d) \t rc %d%d \t %le %le %le\n",
 		     coor[0],coor[1],coor[2],coor[3],r,c,