Threading support rework.

Placed parallel pragmas as macros; implemented deterministic thread reduction in style of
BFM.

TODO

@@ -1,33 +1,35 @@
+================================================================
+*** Hacks and bug fixes to clean up and Audits
+================================================================
+* Base class to share common code between vRealF, VComplexF etc... 
+
+* Performance check on Guido's reimplementation strategy
+
+* Bug in SeedFixedIntegers gives same output on each site. -- Think I fixed but NOT checked for sure
+
+* FIXME audit
+* const audit
+
+* Replace vset with a call to merge.; 
+* care in Gmerge,Gextract over vset .
+* extract / merge extra implementation removal      
+
+* Strong test for norm2, conj and all primitive types. -- tests/Grid_simd.cc is almost there
+
+================================================================
+*** New Functionality
+================================================================
+
+* Implement where to take template scheme.
+
 * - BinaryWriter, TextWriter etc...
   - use protocol buffers? replace xmlReader/Writer ec..
   - Binary use htonll, htonl
 
-*** Hacks and bug fixes to clean up
-* Had to hack assignment to 1.0 in the tests/Grid_gamma test
-* norm2l is a hack. figure out syntax error and make this norm2 c.f. tests/Grid_gamma.cc
-
-* Reduce implemention is poor ; need threaded reductions; OMP isn't able to do it for generic objects.
-  
-* Bug in SeedFixedIntegers gives same output on each site.
-* Bug in RNG with complex numbers ; only filling real values; need helper function -- DONE
-* Conformable test in Cshift routines.
-
-*** Functionality
-
-* Implement where to take template scheme.
-
-* Command line args for geometry, simd, etc. layout. Is it necessary to have
-  user pass these? Is this a QCD specific?
-
-* Strong test for norm2, conj and all primitive types. -- Grid_simd test is almost there
-
 * Expression template engine:
 
-   - Audit
+   -- Audit
-   - Introduce base clase for Grid Tensors.
+   -- Norm2(expression) problem: introduce norm2 unary op, or Introduce conversion automatic from expression to Lattice<vobj>
-   - Introduce norm2 unary op.
-   - Introduce conversion automatic from expression to Lattice<vobj>
-
 
 * CovariantShift support                             -----Use a class to store gauge field? (parallel transport?)
 
@@ -48,12 +50,10 @@
  - I have collated into single location at least.
  - Need to use _mm_*insert/extract routines.
 
-
 * Flavour matrices?
 * Pauli, SU subgroup, etc.. 
 * su3 exponentiation & log etc.. [Jamie's code?]
 * TaProj
-
 * FFTnD ?
 
 * Parallel MPI2 IO
@@ -62,20 +62,23 @@
 * rb4d support for 5th dimension in Mobius.
 
 * Check for missing functionality                    - partially audited against QDP++ layout
-   // Base class to share common code between vRealF, VComplexF etc...
-   //
    // Unary functions
    // cos,sin, tan, acos, asin, cosh, acosh, tanh, sinh, // Scalar<vReal> only arg
    // exp, log, sqrt, fabs
-   //
    // transposeColor, transposeSpin,
    // adjColor, adjSpin,
-   //
    // copyMask.
-   //
    // localMaxAbs
-   //
    // Fourier transform equivalent.
+Actions
+* Fermion
+  - Wilson
+  - Clover
+  - DomainWall
+  - Mobius
+  - z-Mobius
+* Gauge
+  - Wilson, symanzik, iwasaki
 
 Algorithms
 * LinearOperator
@@ -83,23 +86,21 @@ Algorithms
 * Polynomial 
 * Eigen
 * Pcg
+* Adef2
+* DeflCG
 * fPcg
 * MCR
 * HDCG
+* HMC, Heatbath
 * etc..
 
-AUDITS:
+======================================================================================================
-
+FUNCTIONALITY: it pleases me to keep track of things I have done (keeps sane)
-* FIXME audit
-* const audit
-* Replace vset with a call to merge.; 
-* care in Gmerge,Gextract over vset .
-* extract / merge extra implementation removal      
-
-
 ======================================================================================================
 
-FUNCTIONALITY:
+* Command line args for geometry, simd, etc. layout. Is it necessary to have -- DONE
+  user pass these? Is this a QCD specific?
+
 * Stencil -- DONE
 * Test infrastructure -- DONE
 * Fourspin, two spin project --- DONE
@@ -115,6 +116,7 @@ FUNCTIONALITY:
 * How to do U[mu] ... lorentz part of type structure or not. more like chroma if not. -- DONE
 
 * Twospin/Fourspin/Gamma/Proj/Recon                  ----- DONE
+* norm2l is a hack. figure out syntax error and make this norm2 c.f. tests/Grid_gamma.cc -- DONE
 
 * subdirs lib, tests ??                              ----- DONE
   - lib/math        
@@ -149,3 +151,10 @@ FUNCTIONALITY:
 
 * Controling std::cout                              ------- DONE
 
+* Had to hack assignment to 1.0 in the tests/Grid_gamma test -- DONE
+* Reduce implemention is poor ; need threaded reductions; OMP isn't able to do it for generic objects. -- DONE
+* Bug in RNG with complex numbers ; only filling real values; need helper function -- DONE
+* Conformable test in Cshift routines. -- none needed ; there is only one 
+* Conformable testing in expression templates -- DONE (recursive)
+
+

lib/Grid.h

@@ -53,7 +53,6 @@
 #include <Grid_lattice.h>
 #include <Grid_comparison.h>
 #include <Grid_cshift.h>
-#include <Grid_where.h>
 #include <Grid_stencil.h>
 #include <qcd/Grid_qcd.h>
 #include <parallelIO/GridNerscIO.h>

lib/Grid_comparison.h

@@ -142,5 +142,6 @@ namespace Grid {
 }
 
 #include <lattice/Grid_lattice_comparison.h>
+#include <lattice/Grid_lattice_where.h>
 
 #endif

lib/Grid_simd.h

@@ -68,10 +68,14 @@ namespace Grid {
     inline ComplexF  adj(const ComplexF& r ){ return(conj(r)); }
     //conj already supported for complex
 
-    inline ComplexF timesI(const ComplexF r)     { return(r*ComplexF(0.0,1.0));}
+    inline ComplexF timesI(const ComplexF &r)     { return(r*ComplexF(0.0,1.0));}
-    inline ComplexD timesI(const ComplexD r)     { return(r*ComplexD(0.0,1.0));}
+    inline ComplexD timesI(const ComplexD &r)     { return(r*ComplexD(0.0,1.0));}
-    inline ComplexF timesMinusI(const ComplexF r){ return(r*ComplexF(0.0,-1.0));}
+    inline ComplexF timesMinusI(const ComplexF &r){ return(r*ComplexF(0.0,-1.0));}
-    inline ComplexD timesMinusI(const ComplexD r){ return(r*ComplexD(0.0,-1.0));}
+    inline ComplexD timesMinusI(const ComplexD &r){ return(r*ComplexD(0.0,-1.0));}
+    inline void timesI(ComplexF &ret,const ComplexF &r)     { ret = timesI(r);}
+    inline void timesI(ComplexD &ret,const ComplexD &r)     { ret = timesI(r);}
+    inline void timesMinusI(ComplexF &ret,const ComplexF &r){ ret = timesMinusI(r);}
+    inline void timesMinusI(ComplexD &ret,const ComplexD &r){ ret = timesMinusI(r);}
 
     inline void mac (RealD * __restrict__ y,const RealD * __restrict__ a,const RealD *__restrict__ x){  *y = (*a) * (*x)+(*y);}
     inline void mult(RealD * __restrict__ y,const RealD * __restrict__ l,const RealD *__restrict__ r){ *y = (*l) * (*r);}

lib/Makefile.am

@@ -35,8 +35,7 @@ include_HEADERS =\
 	Grid_lattice.h\
 	Grid_math.h\
 	Grid_simd.h\
-	Grid_stencil.h\
+	Grid_stencil.h
-	Grid_where.h
 
 nobase_include_HEADERS=\
 	cartesian/Grid_cartesian_base.h\
@@ -56,6 +55,7 @@ nobase_include_HEADERS=\
 	lattice/Grid_lattice_trace.h\
 	lattice/Grid_lattice_transfer.h\
 	lattice/Grid_lattice_transpose.h\
+	lattice/Grid_lattice_where.h\
 	math/Grid_math_arith.h\
 	math/Grid_math_arith_add.h\
 	math/Grid_math_arith_mac.h\

lib/cshift/Grid_cshift_common.h

@@ -28,7 +28,7 @@ Gather_plane_simple (const Lattice<vobj> &rhs,std::vector<cobj,alignedAllocator<
   int so  = plane*rhs._grid->_ostride[dimension]; // base offset for start of plane 
   int bo  = 0;                                      // offset in buffer
   
-#pragma omp parallel for collapse(2)
+PARALLEL_NESTED_LOOP(2)
   for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
     for(int b=0;b<rhs._grid->_slice_block[dimension];b++){
       int o = n*rhs._grid->_slice_stride[dimension];
@@ -57,7 +57,7 @@ Gather_plane_extract(const Lattice<vobj> &rhs,std::vector<typename cobj::scalar_
   int so  = plane*rhs._grid->_ostride[dimension]; // base offset for start of plane 
   int bo  = 0;                                      // offset in buffer
     
-#pragma omp parallel for collapse(2)
+PARALLEL_NESTED_LOOP(2)
   for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
     for(int b=0;b<rhs._grid->_slice_block[dimension];b++){
 
@@ -106,7 +106,7 @@ template<class vobj> void Scatter_plane_simple (Lattice<vobj> &rhs,std::vector<v
   int so  = plane*rhs._grid->_ostride[dimension]; // base offset for start of plane 
   int bo  = 0;                                      // offset in buffer
     
-#pragma omp parallel for collapse(2)
+PARALLEL_NESTED_LOOP(2)
   for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
     for(int b=0;b<rhs._grid->_slice_block[dimension];b++){
       int o=n*rhs._grid->_slice_stride[dimension];
@@ -131,7 +131,7 @@ template<class vobj> void Scatter_plane_simple (Lattice<vobj> &rhs,std::vector<v
 
   int so  = plane*rhs._grid->_ostride[dimension]; // base offset for start of plane 
     
-#pragma omp parallel for collapse(2)
+PARALLEL_NESTED_LOOP(2)
   for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
     for(int b=0;b<rhs._grid->_slice_block[dimension];b++){
 
@@ -160,7 +160,7 @@ template<class vobj> void Copy_plane(Lattice<vobj>& lhs,Lattice<vobj> &rhs, int
   int ro  = rplane*rhs._grid->_ostride[dimension]; // base offset for start of plane 
   int lo  = lplane*lhs._grid->_ostride[dimension]; // base offset for start of plane 
   
-#pragma omp parallel for collapse(2)
+PARALLEL_NESTED_LOOP(2)
   for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
     for(int b=0;b<rhs._grid->_slice_block[dimension];b++){
       
@@ -185,7 +185,7 @@ template<class vobj> void Copy_plane_permute(Lattice<vobj>& lhs,Lattice<vobj> &r
   int ro  = rplane*rhs._grid->_ostride[dimension]; // base offset for start of plane 
   int lo  = lplane*lhs._grid->_ostride[dimension]; // base offset for start of plane 
   
-#pragma omp parallel for collapse(2)
+PARALLEL_NESTED_LOOP(2)
   for(int n=0;n<rhs._grid->_slice_nblock[dimension];n++){
     for(int b=0;b<rhs._grid->_slice_block[dimension];b++){
       int o =n*rhs._grid->_slice_stride[dimension];

lib/lattice/Grid_lattice_arith.h

@@ -10,7 +10,7 @@ namespace Grid {
   template<class obj1,class obj2,class obj3>
     void mult(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const Lattice<obj3> &rhs){
     conformable(lhs,rhs);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
       obj1 tmp;
       mult(&tmp,&lhs._odata[ss],&rhs._odata[ss]);
@@ -21,7 +21,7 @@ namespace Grid {
   template<class obj1,class obj2,class obj3>
     void mac(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const Lattice<obj3> &rhs){
     conformable(lhs,rhs);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
       obj1 tmp;
       mac(&tmp,&lhs._odata[ss],&rhs._odata[ss]);
@@ -32,7 +32,7 @@ namespace Grid {
   template<class obj1,class obj2,class obj3>
     void sub(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const Lattice<obj3> &rhs){
     conformable(lhs,rhs);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
       obj1 tmp;
       sub(&tmp,&lhs._odata[ss],&rhs._odata[ss]);
@@ -42,7 +42,7 @@ namespace Grid {
   template<class obj1,class obj2,class obj3>
     void add(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const Lattice<obj3> &rhs){
     conformable(lhs,rhs);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
       obj1 tmp;
       add(&tmp,&lhs._odata[ss],&rhs._odata[ss]);
@@ -56,7 +56,7 @@ namespace Grid {
   template<class obj1,class obj2,class obj3>
     void mult(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const obj3 &rhs){
     conformable(lhs,ret);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
       obj1 tmp;
       mult(&tmp,&lhs._odata[ss],&rhs);
@@ -67,7 +67,7 @@ namespace Grid {
   template<class obj1,class obj2,class obj3>
     void mac(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const obj3 &rhs){
     conformable(lhs,ret);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
       obj1 tmp;
       mac(&tmp,&lhs._odata[ss],&rhs);
@@ -78,7 +78,7 @@ namespace Grid {
   template<class obj1,class obj2,class obj3>
     void sub(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const obj3 &rhs){
     conformable(lhs,ret);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
       obj1 tmp;
       sub(&tmp,&lhs._odata[ss],&rhs);
@@ -88,7 +88,7 @@ namespace Grid {
   template<class obj1,class obj2,class obj3>
     void add(Lattice<obj1> &ret,const Lattice<obj2> &lhs,const obj3 &rhs){
     conformable(lhs,ret);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
       obj1 tmp;
       add(&tmp,&lhs._odata[ss],&rhs);
@@ -102,7 +102,7 @@ namespace Grid {
     template<class obj1,class obj2,class obj3>
     void mult(Lattice<obj1> &ret,const obj2 &lhs,const Lattice<obj3> &rhs){
     conformable(ret,rhs);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<rhs._grid->oSites();ss++){
       obj1 tmp;
       mult(&tmp,&lhs,&rhs._odata[ss]);
@@ -113,7 +113,7 @@ namespace Grid {
   template<class obj1,class obj2,class obj3>
     void mac(Lattice<obj1> &ret,const obj2 &lhs,const Lattice<obj3> &rhs){
     conformable(ret,rhs);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<rhs._grid->oSites();ss++){
       obj1 tmp;
       mac(&tmp,&lhs,&rhs._odata[ss]);
@@ -124,7 +124,7 @@ namespace Grid {
   template<class obj1,class obj2,class obj3>
     void sub(Lattice<obj1> &ret,const obj2 &lhs,const Lattice<obj3> &rhs){
     conformable(ret,rhs);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<rhs._grid->oSites();ss++){
       obj1 tmp;
       sub(&tmp,&lhs,&rhs._odata[ss]);
@@ -134,7 +134,7 @@ namespace Grid {
   template<class obj1,class obj2,class obj3>
     void add(Lattice<obj1> &ret,const obj2 &lhs,const Lattice<obj3> &rhs){
     conformable(ret,rhs);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<rhs._grid->oSites();ss++){
       obj1 tmp;
       add(&tmp,&lhs,&rhs._odata[ss]);
@@ -145,7 +145,7 @@ namespace Grid {
   template<class sobj,class vobj>
   inline void axpy(Lattice<vobj> &ret,sobj a,const Lattice<vobj> &lhs,const Lattice<vobj> &rhs){
     conformable(lhs,rhs);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<lhs._grid->oSites();ss++){
       vobj tmp = a*lhs._odata[ss];
       vstream(ret._odata[ss],tmp+rhs._odata[ss]);

lib/lattice/Grid_lattice_base.h

@@ -64,7 +64,7 @@ public:
   ////////////////////////////////////////////////////////////////////////////////
   template <typename Op, typename T1>                         inline Lattice<vobj> & operator=(const LatticeUnaryExpression<Op,T1> &expr)
   {
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
       vobj tmp= eval(ss,expr);
       vstream(_odata[ss] ,tmp);
@@ -73,7 +73,7 @@ public:
   }
   template <typename Op, typename T1,typename T2>             inline Lattice<vobj> & operator=(const LatticeBinaryExpression<Op,T1,T2> &expr)
   {
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
       vobj tmp= eval(ss,expr);
       vstream(_odata[ss] ,tmp);
@@ -82,7 +82,7 @@ public:
   }
   template <typename Op, typename T1,typename T2,typename T3> inline Lattice<vobj> & operator=(const LatticeTrinaryExpression<Op,T1,T2,T3> &expr)
   {
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
       vobj tmp= eval(ss,expr);
       vstream(_odata[ss] ,tmp);
@@ -95,7 +95,7 @@ public:
     GridFromExpression(_grid,expr);
     assert(_grid!=nullptr);
     _odata.resize(_grid->oSites());
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
       _odata[ss] = eval(ss,expr);
     }
@@ -105,7 +105,7 @@ public:
     GridFromExpression(_grid,expr);
     assert(_grid!=nullptr);
     _odata.resize(_grid->oSites());
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
       _odata[ss] = eval(ss,expr);
     }
@@ -115,7 +115,7 @@ public:
     GridFromExpression(_grid,expr);
     assert(_grid!=nullptr);
     _odata.resize(_grid->oSites());
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<_grid->oSites();ss++){
       _odata[ss] = eval(ss,expr);
     }
@@ -133,7 +133,7 @@ public:
     }
 
     template<class sobj> inline Lattice<vobj> & operator = (const sobj & r){
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<_grid->oSites();ss++){
             this->_odata[ss]=r;
         }
@@ -142,7 +142,7 @@ public:
     template<class robj> inline Lattice<vobj> & operator = (const Lattice<robj> & r){
       conformable(*this,r);
       std::cout<<"Lattice operator ="<<std::endl;
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<_grid->oSites();ss++){
             this->_odata[ss]=r._odata[ss];
         }
@@ -167,7 +167,7 @@ public:
     inline friend Lattice<vobj> operator / (const Lattice<vobj> &lhs,const Lattice<vobj> &rhs){
         conformable(lhs,rhs);
         Lattice<vobj> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
             ret._odata[ss] = lhs._odata[ss]/rhs._odata[ss];
         }

lib/lattice/Grid_lattice_comparison.h

@@ -15,7 +15,7 @@ namespace Grid {
     inline Lattice<vInteger> LLComparison(vfunctor op,const Lattice<lobj> &lhs,const Lattice<robj> &rhs)
     {
       Lattice<vInteger> ret(rhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<rhs._grid->oSites(); ss++){
 	  ret._odata[ss]=op(lhs._odata[ss],rhs._odata[ss]);
         }
@@ -25,7 +25,7 @@ namespace Grid {
     inline Lattice<vInteger> LSComparison(vfunctor op,const Lattice<lobj> &lhs,const robj &rhs)
     {
       Lattice<vInteger> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites(); ss++){
 	  ret._odata[ss]=op(lhs._odata[ss],rhs);
         }
@@ -35,7 +35,7 @@ namespace Grid {
     inline Lattice<vInteger> SLComparison(vfunctor op,const lobj &lhs,const Lattice<robj> &rhs)
     {
       Lattice<vInteger> ret(rhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<rhs._grid->oSites(); ss++){
 	  ret._odata[ss]=op(lhs._odata[ss],rhs);
         }

lib/lattice/Grid_lattice_local.h

@@ -16,7 +16,7 @@ namespace Grid {
     inline auto localNorm2 (const Lattice<vobj> &rhs)-> Lattice<typename vobj::tensor_reduced>
     {
       Lattice<typename vobj::tensor_reduced> ret(rhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<rhs._grid->oSites(); ss++){
 	  ret._odata[ss]=innerProduct(rhs._odata[ss],rhs._odata[ss]);
         }
@@ -29,7 +29,7 @@ namespace Grid {
       -> Lattice<typename vobj::tensor_reduced>
     {
       Lattice<typename vobj::tensor_reduced> ret(rhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
       for(int ss=0;ss<rhs._grid->oSites(); ss++){
 	ret._odata[ss]=innerProduct(lhs._odata[ss],rhs._odata[ss]);
       }
@@ -42,7 +42,7 @@ namespace Grid {
     inline auto outerProduct (const Lattice<ll> &lhs,const Lattice<rr> &rhs) -> Lattice<decltype(outerProduct(lhs._odata[0],rhs._odata[0]))>
     {
         Lattice<decltype(outerProduct(lhs._odata[0],rhs._odata[0]))> ret(rhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<rhs._grid->oSites(); ss++){
             ret._odata[ss]=outerProduct(lhs._odata[ss],rhs._odata[ss]);
         }

lib/lattice/Grid_lattice_overload.h

@@ -10,7 +10,7 @@ namespace Grid {
   inline Lattice<vobj> operator -(const Lattice<vobj> &r)
   {
     Lattice<vobj> ret(r._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<r._grid->oSites();ss++){
       vstream(ret._odata[ss], -r._odata[ss]);
     }
@@ -47,7 +47,7 @@ namespace Grid {
   inline auto operator * (const left &lhs,const Lattice<right> &rhs) -> Lattice<decltype(lhs*rhs._odata[0])>
   {
     Lattice<decltype(lhs*rhs._odata[0])> ret(rhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<rhs._grid->oSites(); ss++){
       decltype(lhs*rhs._odata[0]) tmp=lhs*rhs._odata[ss]; 
       vstream(ret._odata[ss],tmp);
@@ -59,7 +59,7 @@ namespace Grid {
     inline auto operator + (const left &lhs,const Lattice<right> &rhs) -> Lattice<decltype(lhs+rhs._odata[0])>
     {
       Lattice<decltype(lhs+rhs._odata[0])> ret(rhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
       for(int ss=0;ss<rhs._grid->oSites(); ss++){
 	decltype(lhs+rhs._odata[0]) tmp =lhs-rhs._odata[ss];  
 	vstream(ret._odata[ss],tmp);
@@ -71,7 +71,7 @@ namespace Grid {
     inline auto operator - (const left &lhs,const Lattice<right> &rhs) -> Lattice<decltype(lhs-rhs._odata[0])>
   {
     Lattice<decltype(lhs-rhs._odata[0])> ret(rhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<rhs._grid->oSites(); ss++){
       decltype(lhs-rhs._odata[0]) tmp=lhs-rhs._odata[ss];  
       vstream(ret._odata[ss],tmp);
@@ -83,7 +83,7 @@ namespace Grid {
       inline auto operator * (const Lattice<left> &lhs,const right &rhs) -> Lattice<decltype(lhs._odata[0]*rhs)>
     {
       Lattice<decltype(lhs._odata[0]*rhs)> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
       for(int ss=0;ss<lhs._grid->oSites(); ss++){
 	decltype(lhs._odata[0]*rhs) tmp =lhs._odata[ss]*rhs;
 	vstream(ret._odata[ss],tmp);
@@ -95,7 +95,7 @@ namespace Grid {
       inline auto operator + (const Lattice<left> &lhs,const right &rhs) -> Lattice<decltype(lhs._odata[0]+rhs)>
     {
         Lattice<decltype(lhs._odata[0]+rhs)> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<rhs._grid->oSites(); ss++){
 	  decltype(lhs._odata[0]+rhs) tmp=lhs._odata[ss]+rhs; 
 	  vstream(ret._odata[ss],tmp);
@@ -107,7 +107,7 @@ namespace Grid {
       inline auto operator - (const Lattice<left> &lhs,const right &rhs) -> Lattice<decltype(lhs._odata[0]-rhs)>
     {
       Lattice<decltype(lhs._odata[0]-rhs)> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
       for(int ss=0;ss<rhs._grid->oSites(); ss++){
 	  decltype(lhs._odata[0]-rhs) tmp=lhs._odata[ss]-rhs;
 	  vstream(ret._odata[ss],tmp);

lib/lattice/Grid_lattice_peekpoke.h

@@ -15,7 +15,7 @@ namespace Grid {
       -> Lattice<decltype(peekIndex<Index>(lhs._odata[0]))>
     {
       Lattice<decltype(peekIndex<Index>(lhs._odata[0]))> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
             ret._odata[ss] = peekIndex<Index>(lhs._odata[ss]);
         }
@@ -26,7 +26,7 @@ namespace Grid {
       -> Lattice<decltype(peekIndex<Index>(lhs._odata[0],i))>
     {
       Lattice<decltype(peekIndex<Index>(lhs._odata[0],i))> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
 	  ret._odata[ss] = peekIndex<Index>(lhs._odata[ss],i);
         }
@@ -37,7 +37,7 @@ namespace Grid {
       -> Lattice<decltype(peekIndex<Index>(lhs._odata[0],i,j))>
     {
       Lattice<decltype(peekIndex<Index>(lhs._odata[0],i,j))> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
 	  ret._odata[ss] = peekIndex<Index>(lhs._odata[ss],i,j);
         }
@@ -50,7 +50,7 @@ namespace Grid {
     template<int Index,class vobj> inline
     void pokeIndex(Lattice<vobj> &lhs,const Lattice<decltype(peekIndex<Index>(lhs._odata[0]))> & rhs)
     {
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
 	  pokeIndex<Index>(lhs._odata[ss],rhs._odata[ss]);
 	}      
@@ -58,7 +58,7 @@ namespace Grid {
     template<int Index,class vobj> inline
     void pokeIndex(Lattice<vobj> &lhs,const Lattice<decltype(peekIndex<Index>(lhs._odata[0],0))> & rhs,int i)
     {
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
 	  pokeIndex<Index>(lhs._odata[ss],rhs._odata[ss],i);
 	}      
@@ -66,7 +66,7 @@ namespace Grid {
     template<int Index,class vobj> inline
     void pokeIndex(Lattice<vobj> &lhs,const Lattice<decltype(peekIndex<Index>(lhs._odata[0],0,0))> & rhs,int i,int j)
     {
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
 	  pokeIndex<Index>(lhs._odata[ss],rhs._odata[ss],i,j);
 	}      

lib/lattice/Grid_lattice_reality.h

@@ -11,7 +11,7 @@ namespace Grid {
 
     template<class vobj> inline Lattice<vobj> adj(const Lattice<vobj> &lhs){
         Lattice<vobj> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
             ret._odata[ss] = adj(lhs._odata[ss]);
         }
@@ -20,7 +20,7 @@ namespace Grid {
 
     template<class vobj> inline Lattice<vobj> conj(const Lattice<vobj> &lhs){
         Lattice<vobj> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
             ret._odata[ss] = conj(lhs._odata[ss]);
         }
@@ -30,7 +30,7 @@ namespace Grid {
     template<class vobj> inline auto real(const Lattice<vobj> &z) -> Lattice<decltype(real(z._odata[0]))>
     {
       Lattice<decltype(real(z._odata[0]))> ret(z._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<z._grid->oSites();ss++){
             ret._odata[ss] = real(z._odata[ss]);
         }
@@ -40,7 +40,7 @@ namespace Grid {
     template<class vobj> inline auto imag(const Lattice<vobj> &z) -> Lattice<decltype(imag(z._odata[0]))>
     {
       Lattice<decltype(imag(z._odata[0]))> ret(z._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<z._grid->oSites();ss++){
             ret._odata[ss] = imag(z._odata[ss]);
         }

lib/lattice/Grid_lattice_reduction.h

@@ -7,69 +7,85 @@ namespace Grid {
 #endif     
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////
-    // Reduction operations
+    // Deterministic Reduction operations
     ////////////////////////////////////////////////////////////////////////////////////////////////////
-    template<class vobj>
+  template<class vobj> inline RealD norm2(const Lattice<vobj> &arg){
-    inline RealD norm2(const Lattice<vobj> &arg){
+    ComplexD nrm = innerProduct(arg,arg);
-
+    return real(nrm); 
-      typedef typename vobj::scalar_type scalar;
+  }
-      typedef typename vobj::vector_type vector;
-      decltype(innerProduct(arg._odata[0],arg._odata[0])) vnrm; 
-      scalar nrm;
-      //FIXME make this loop parallelisable
-      vnrm=zero;
-      for(int ss=0;ss<arg._grid->oSites(); ss++){
-	vnrm = vnrm + innerProduct(arg._odata[ss],arg._odata[ss]);
-      }
-      vector vvnrm =TensorRemove(vnrm) ;
-      nrm = Reduce(vvnrm);
-      arg._grid->GlobalSum(nrm);
-      return real(nrm);
-    }
 
     template<class vobj>
     inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right) 
-      //    inline auto innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right) 
-      //->decltype(innerProduct(left._odata[0],right._odata[0]))
     {
-      typedef typename vobj::scalar_type scalar;
+      typedef typename vobj::scalar_type scalar_type;
-      decltype(innerProduct(left._odata[0],right._odata[0])) vnrm; 
+      typedef typename vobj::vector_type vector_type;
+      vector_type vnrm;
+      scalar_type  nrm;
 
-      scalar nrm;
+      GridBase *grid = left._grid;
-      //FIXME make this loop parallelisable
+
-      vnrm=zero;
+      std::vector<vector_type,alignedAllocator<vector_type> > sumarray(grid->SumArraySize());
-      for(int ss=0;ss<left._grid->oSites(); ss++){
+      for(int i=0;i<grid->SumArraySize();i++){
-	vnrm = vnrm + innerProduct(left._odata[ss],right._odata[ss]);
+	sumarray[i]=zero;
       }
-      nrm = Reduce(vnrm);
+
+PARALLEL_FOR_LOOP
+      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]);
+	}
+	sumarray[thr]=TensorRemove(vnrm) ;
+      }
+    
+      vector_type vvnrm; vvnrm=zero;  // sum across threads
+      for(int i=0;i<grid->SumArraySize();i++){
+	vvnrm = vvnrm+sumarray[i];
+      } 
+      nrm = Reduce(vvnrm);// sum across simd
       right._grid->GlobalSum(nrm);
       return nrm;
     }
 
     template<class vobj>
-      inline typename vobj::scalar_object sum(const Lattice<vobj> &arg){
+    inline typename vobj::scalar_object sum(const Lattice<vobj> &arg){
 
       GridBase *grid=arg._grid;
       int Nsimd = grid->Nsimd();
 
-      typedef typename vobj::scalar_object sobj;
+      std::vector<vobj,alignedAllocator<vobj> > sumarray(grid->SumArraySize());
-      typedef typename vobj::scalar_type   scalar_type;
+      for(int i=0;i<grid->SumArraySize();i++){
-
+	sumarray[i]=zero;
-      vobj vsum;
-      sobj ssum;
-
-      vsum=zero;
-      ssum=zero;
-      //FIXME make this loop parallelisable
-      for(int ss=0;ss<arg._grid->oSites(); ss++){
-	vsum = vsum + arg._odata[ss];
       }
-      
+
+PARALLEL_FOR_LOOP
+      for(int thr=0;thr<grid->SumArraySize();thr++){
+	int nwork, mywork, myoff;
+	GridThread::GetWork(grid->oSites(),thr,mywork,myoff);
+
+	vobj vvsum=zero;
+        for(int ss=myoff;ss<mywork+myoff; ss++){
+	  vvsum = vvsum + arg._odata[ss];
+	}
+	sumarray[thr]=vvsum;
+      }
+
+      vobj vsum=zero;  // sum across threads
+      for(int i=0;i<grid->SumArraySize();i++){
+	vsum = vsum+sumarray[i];
+      } 
+
+      typedef typename vobj::scalar_object sobj;
+      sobj ssum=zero;
+
       std::vector<sobj>               buf(Nsimd);
       extract(vsum,buf);
 
       for(int i=0;i<Nsimd;i++) ssum = ssum + buf[i];
-
       arg._grid->GlobalSum(ssum);
 
       return ssum;
@@ -77,13 +93,15 @@ namespace Grid {
 
 
 
-
 template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<typename vobj::scalar_object> &result,int orthogdim)
 {
   typedef typename vobj::scalar_object sobj;
-
   GridBase  *grid = Data._grid;
   assert(grid!=NULL);
+
+  // FIXME
+  std::cout<<"WARNING ! SliceSum is unthreaded "<<grid->SumArraySize()<<" threads "<<std::endl;
+
   const int    Nd = grid->_ndimension;
   const int Nsimd = grid->Nsimd();
 
@@ -94,10 +112,8 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
   int ld=grid->_ldimensions[orthogdim];
   int rd=grid->_rdimensions[orthogdim];
 
-  sobj szero; szero=zero;
-
   std::vector<vobj,alignedAllocator<vobj> > lvSum(rd); // will locally sum vectors first
-  std::vector<sobj> lsSum(ld,szero); // sum across these down to scalars
+  std::vector<sobj> lsSum(ld,zero); // sum across these down to scalars
   std::vector<sobj> extracted(Nsimd);     // splitting the SIMD
 
   result.resize(fd); // And then global sum to return the same vector to every node for IO to file
@@ -108,6 +124,7 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
   std::vector<int>  coor(Nd);  
 
   // sum over reduced dimension planes, breaking out orthog dir
+
   for(int ss=0;ss<grid->oSites();ss++){
     GridBase::CoorFromIndex(coor,ss,grid->_rdimensions);
     int r = coor[orthogdim];

lib/lattice/Grid_lattice_trace.h

@@ -15,7 +15,7 @@ namespace Grid {
       -> Lattice<decltype(trace(lhs._odata[0]))>
     {
       Lattice<decltype(trace(lhs._odata[0]))> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
             ret._odata[ss] = trace(lhs._odata[ss]);
         }
@@ -30,7 +30,7 @@ namespace Grid {
       -> Lattice<decltype(traceIndex<Index>(lhs._odata[0]))>
     {
       Lattice<decltype(traceIndex<Index>(lhs._odata[0]))> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
       for(int ss=0;ss<lhs._grid->oSites();ss++){
 	ret._odata[ss] = traceIndex<Index>(lhs._odata[ss]);
       }

lib/lattice/Grid_lattice_transfer.h

@@ -23,7 +23,7 @@ inline void subdivides(GridBase *coarse,GridBase *fine)
   template<class vobj> inline void pickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full){
     half.checkerboard = cb;
     int ssh=0;
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<full._grid->oSites();ss++){
       std::vector<int> coor;
       int cbos;
@@ -41,7 +41,7 @@ inline void subdivides(GridBase *coarse,GridBase *fine)
   template<class vobj> inline void setCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half){
     int cb = half.checkerboard;
     int ssh=0;
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
     for(int ss=0;ss<full._grid->oSites();ss++){
       std::vector<int> coor;
       int cbos;

lib/lattice/Grid_lattice_transpose.h

@@ -13,7 +13,7 @@ namespace Grid {
   template<class vobj>
     inline Lattice<vobj> transpose(const Lattice<vobj> &lhs){
         Lattice<vobj> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
             ret._odata[ss] = transpose(lhs._odata[ss]);
         }
@@ -28,7 +28,7 @@ namespace Grid {
       -> Lattice<decltype(transposeIndex<Index>(lhs._odata[0]))>
     {
       Lattice<decltype(transposeIndex<Index>(lhs._odata[0]))> ret(lhs._grid);
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
         for(int ss=0;ss<lhs._grid->oSites();ss++){
             ret._odata[ss] = transposeIndex<Index>(lhs._odata[ss]);
         }

lib/lattice/Grid_lattice_where.h

@@ -1,5 +1,5 @@
-#ifndef GRID_WHERE_H
+#ifndef GRID_LATTICE_WHERE_H
-#define GRID_WHERE_H
+#define GRID_LATTICE_WHERE_H
 namespace Grid {
 // Must implement the predicate gating the 
 // Must be able to reduce the predicate down to a single vInteger per site.
@@ -27,7 +27,7 @@ inline void where(Lattice<vobj> &ret,const Lattice<iobj> &predicate,Lattice<vobj
   std::vector<scalar_object> truevals (Nsimd);
   std::vector<scalar_object> falsevals(Nsimd);
 
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
   for(int ss=0;ss<iftrue._grid->oSites(); ss++){
 
     extract(iftrue._odata[ss]   ,truevals);

lib/math/Grid_math_reality.h

@@ -6,23 +6,20 @@ namespace Grid {
     /////////////////////////////////////////// CONJ         ///////////////////////////////////////////
     ///////////////////////////////////////////////////////////////////////////////////////////////////
 
-#ifdef GRID_WARN_SUBOPTIMAL
-#warning "Optimisation alert switch over to two argument form to avoid copy back in perf critical timesI "
-#endif     
 /////////////////////////////////////////////// 
 // multiply by I; make recursive.
 /////////////////////////////////////////////// 
 template<class vtype> inline iScalar<vtype> timesI(const iScalar<vtype>&r) 
 {
     iScalar<vtype> ret;
-    ret._internal = timesI(r._internal);
+    timesI(ret._internal,r._internal);
     return ret;
 }
 template<class vtype,int N> inline iVector<vtype,N> timesI(const iVector<vtype,N>&r) 
 {
   iVector<vtype,N> ret;
   for(int i=0;i<N;i++){
-    ret._internal[i] = timesI(r._internal[i]);
+    timesI(ret._internal[i],r._internal[i]);
   }
   return ret;
 }
@@ -31,22 +28,41 @@ template<class vtype,int N> inline iMatrix<vtype,N> timesI(const iMatrix<vtype,N
   iMatrix<vtype,N> ret;
   for(int i=0;i<N;i++){
   for(int j=0;j<N;j++){
-    ret._internal[i][j] = timesI(r._internal[i][j]);
+    timesI(ret._internal[i][j],r._internal[i][j]);
   }}
   return ret;
 }
 
+template<class vtype> inline void timesI(iScalar<vtype> &ret,const iScalar<vtype>&r) 
+{
+  timesI(ret._internal,r._internal);
+}
+template<class vtype,int N> inline void timesI(iVector<vtype,N> &ret,const iVector<vtype,N>&r) 
+{
+  for(int i=0;i<N;i++){
+    timesI(ret._internal[i],r._internal[i]);
+  }
+}
+template<class vtype,int N> inline void  timesI(iMatrix<vtype,N> &ret,const iMatrix<vtype,N>&r)
+{
+  for(int i=0;i<N;i++){
+  for(int j=0;j<N;j++){
+    timesI(ret._internal[i][j],r._internal[i][j]);
+  }}
+}
+
+
 template<class vtype> inline iScalar<vtype> timesMinusI(const iScalar<vtype>&r) 
 {
     iScalar<vtype> ret;
-    ret._internal = timesMinusI(r._internal);
+    timesMinusI(ret._internal,r._internal);
     return ret;
 }
 template<class vtype,int N> inline iVector<vtype,N> timesMinusI(const iVector<vtype,N>&r) 
 {
   iVector<vtype,N> ret;
   for(int i=0;i<N;i++){
-    ret._internal[i] = timesMinusI(r._internal[i]);
+    timesMinusI(ret._internal[i],r._internal[i]);
   }
   return ret;
 }
@@ -55,10 +71,30 @@ template<class vtype,int N> inline iMatrix<vtype,N> timesMinusI(const iMatrix<vt
   iMatrix<vtype,N> ret;
   for(int i=0;i<N;i++){
   for(int j=0;j<N;j++){
-    ret._internal[i][j] = timesMinusI(r._internal[i][j]);
+    timesMinusI(ret._internal[i][j],r._internal[i][j]);
   }}
   return ret;
 }
+
+template<class vtype>  inline void timesMinusI(iScalar<vtype> &ret,const iScalar<vtype>&r) 
+{
+  timesMinusI(ret._internal,r._internal);
+}
+template<class vtype,int N> inline void timesMinusI(iVector<vtype,N> &ret,const iVector<vtype,N>&r) 
+{
+  for(int i=0;i<N;i++){
+    timesMinusI(ret._internal[i],r._internal[i]);
+  }
+}
+template<class vtype,int N> inline void  timesMinusI(iMatrix<vtype,N> &ret,const iMatrix<vtype,N>&r)
+{
+  for(int i=0;i<N;i++){
+  for(int j=0;j<N;j++){
+    timesMinusI(ret._internal[i][j],r._internal[i][j]);
+  }}
+}
+
+
 /////////////////////////////////////////////// 
 // Conj function for scalar, vector, matrix
 /////////////////////////////////////////////// 

lib/qcd/Grid_qcd_2spinor.h

@@ -1,9 +1,7 @@
 #ifndef GRID_QCD_TWOSPIN_H
 #define GRID_QCD_TWOSPIN_H
 namespace Grid{
-
 namespace QCD {
-  
 
   //////////////////////////////////////////////////////////////////////////////////////////////////////
   // Normalisation alert; the g5   project is 1/2(1+-G5) 
@@ -1073,10 +1071,6 @@ namespace QCD {
       spProj5m(hspin._internal[i][j],fspin._internal[i][j]);
     }}
   }
-
-
-
-
 }   //namespace QCD
 } // Grid
 #endif

lib/qcd/Grid_qcd_wilson_dop.cc

@@ -103,7 +103,7 @@ void WilsonMatrix::Dhop(const LatticeFermion &in, LatticeFermion &out)
   vHalfSpinColourVector *chi_p;
   int offset,local,perm, ptype;
 
-#pragma omp parallel for
+PARALLEL_FOR_LOOP
   for(int sss=0;sss<grid->oSites();sss++){
 
     int ss = sss;

lib/simd/Grid_vComplexD.h

@@ -13,6 +13,9 @@ namespace Grid {
             vzero(*this);
             return (*this);
         }
+        vComplexD( Zero & z){
+            vzero(*this);
+        }
         vComplexD()=default;
         vComplexD(ComplexD a){
 	  vsplat(*this,a);
@@ -286,8 +289,8 @@ friend inline void vstore(const vComplexD &ret, ComplexD *a){
             return ret;
         }
 
-        friend inline vComplexD timesMinusI(const vComplexD &in){
+        friend inline void timesMinusI(vComplexD &ret,const vComplexD &in){
-	  vComplexD ret; vzero(ret);
+	  vzero(ret);
 	  vComplexD tmp;
 #if defined (AVX1)|| defined (AVX2)
 	  tmp.v    =_mm256_addsub_pd(ret.v,in.v); // r,-i
@@ -304,11 +307,10 @@ friend inline void vstore(const vComplexD &ret, ComplexD *a){
 #ifdef QPX
             assert(0);
 #endif
-	  return ret;
 	}
 
-        friend inline vComplexD timesI(const vComplexD &in){
+        friend inline void timesI(vComplexD &ret, const vComplexD &in){
-	  vComplexD ret; vzero(ret);
+	  vzero(ret);
 	  vComplexD tmp;
 #if defined (AVX1)|| defined (AVX2)
 	  tmp.v    =_mm256_shuffle_pd(in.v,in.v,0x5);
@@ -325,9 +327,21 @@ friend inline void vstore(const vComplexD &ret, ComplexD *a){
 #ifdef QPX
             assert(0);
 #endif
+	}
+
+        friend inline vComplexD timesMinusI(const vComplexD &in){
+	  vComplexD ret; 
+	  timesMinusI(ret,in);
 	  return ret;
 	}
 
+        friend inline vComplexD timesI(const vComplexD &in){
+	  vComplexD ret; 
+	  timesI(ret,in);
+	  return ret;
+	}
+
+
 // REDUCE FIXME must be a cleaner implementation
        friend inline ComplexD Reduce(const vComplexD & in)
        { 

lib/simd/Grid_vComplexF.h

@@ -28,6 +28,9 @@ namespace Grid {
             vzero(*this);
             return (*this);
         }
+        vComplexF( Zero & z){
+            vzero(*this);
+        }
         vComplexF()=default;
         vComplexF(ComplexF a){
 	  vsplat(*this,a);
@@ -363,8 +366,7 @@ namespace Grid {
 #endif
             return ret;
         }
-        friend inline vComplexF timesMinusI(const vComplexF &in){
+        friend inline void timesMinusI( vComplexF &ret,const vComplexF &in){
-	  vComplexF ret; 
 	  vzero(ret);
 #if defined (AVX1)|| defined (AVX2)
 	  cvec tmp =_mm256_addsub_ps(ret.v,in.v); // r,-i
@@ -381,10 +383,9 @@ namespace Grid {
 #ifdef QPX
             assert(0);
 #endif
-	  return ret;
 	}
-        friend inline vComplexF timesI(const vComplexF &in){
+        friend inline void timesI(vComplexF &ret,const vComplexF &in){
-	  vComplexF ret; vzero(ret);
+	  vzero(ret);
 #if defined (AVX1)|| defined (AVX2)
 	  cvec tmp =_mm256_shuffle_ps(in.v,in.v,_MM_SHUFFLE(2,3,0,1));//i,r
           ret.v    =_mm256_addsub_ps(ret.v,tmp);     //i,-r
@@ -400,8 +401,18 @@ namespace Grid {
 #ifdef QPX
             assert(0);
 #endif
+	}
+        friend inline vComplexF timesMinusI(const vComplexF &in){
+	  vComplexF ret; 
+	  timesMinusI(ret,in);
 	  return ret;
 	}
+        friend inline vComplexF timesI(const vComplexF &in){
+	  vComplexF ret; 
+	  timesI(ret,in);
+	  return ret;
+	}
+
         
         // Unary negation
         friend inline vComplexF operator -(const vComplexF &r) {

lib/simd/Grid_vRealD.h

@@ -17,6 +17,10 @@ namespace Grid {
         vRealD(Zero &zero){
 	  zeroit(*this);
 	}
+        vRealD & operator = ( Zero & z){
+	  vzero(*this);
+	  return (*this);
+        }
 
         friend inline void mult(vRealD * __restrict__ y,const vRealD * __restrict__ l,const vRealD *__restrict__ r) {*y = (*l) * (*r);}
         friend inline void sub (vRealD * __restrict__ y,const vRealD * __restrict__ l,const vRealD *__restrict__ r) {*y = (*l) - (*r);}

lib/simd/Grid_vRealF.h

@@ -18,6 +18,10 @@ namespace Grid {
         vRealF(Zero &zero){
 	  zeroit(*this);
 	}
+        vRealF & operator = ( Zero & z){
+	  vzero(*this);
+	  return (*this);
+        }
         ////////////////////////////////////
         // Arithmetic operator overloads +,-,*
         ////////////////////////////////////

tests/Grid_gamma.cc

@@ -50,9 +50,8 @@ int main (int argc, char ** argv)
   //  std::cout << " Is triv Scalar<double> " <<std::has_trivial_default_constructor<iScalar<double> >::value << std::endl;
   //  std::cout << " Is triv Scalar<vComplexD> "<<std::has_trivial_default_constructor<iScalar<vComplexD> >::value  << std::endl;
 
-  std::complex<float> c(1.0);
   for(int a=0;a<Ns;a++){
-    ident()(a,a) = c;
+    ident()(a,a) = ComplexF(1.0);
   }
 
   const Gamma::GammaMatrix *g = Gamma::GammaMatrices;