Debugged the real() and imag() functions and added tests to Test_Simd

lib/Simd.h

@@ -24,7 +24,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 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
+See the full license in the file "LICENSE" in the top level distribution
+directory
 *************************************************************************************/
 /*  END LEGAL */
 #ifndef GRID_SIMD_H
@@ -118,6 +119,14 @@ namespace Grid {
   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 ComplexD timesMinusI(const ComplexD &r){ return(r*ComplexD(0.0,-1.0));}
+
+  // define projections to real and imaginay parts
+  inline ComplexF projReal(const ComplexF &r){return( ComplexF(std::real(r), 0.0));}
+  inline ComplexD projReal(const ComplexD &r){return( ComplexD(std::real(r), 0.0));}
+  inline ComplexF projImag(const ComplexF &r){return (ComplexF(std::imag(r), 0.0 ));}
+  inline ComplexD projImag(const ComplexD &r){return (ComplexD(std::imag(r), 0.0));}
+
+  // define auxiliary functions for complex computations
   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);}

lib/lattice/Lattice_reduction.h

@@ -40,7 +40,7 @@ namespace Grid {
     ////////////////////////////////////////////////////////////////////////////////////////////////////
   template<class vobj> inline RealD norm2(const Lattice<vobj> &arg){
     ComplexD nrm = innerProduct(arg,arg);
-    return real(nrm); 
+    return std::real(nrm); 
   }
 
     template<class vobj>

lib/qcd/smearing/StoutSmearing.h

@@ -18,12 +18,14 @@
   INHERIT_GIMPL_TYPES(Gimpl)
 
   Smear_Stout(Smear<Gimpl>* base) : SmearBase(base) {
-  				static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3");
+    static_assert(Nc == 3,
+                  "Stout smearing currently implemented only for Nc==3");
   }
 
   /*! Default constructor */
   Smear_Stout(double rho = 1.0) : SmearBase(new Smear_APE<Gimpl>(rho)) {
-  				static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3");
+    static_assert(Nc == 3,
+                  "Stout smearing currently implemented only for Nc==3");
   }
 
   ~Smear_Stout() {}  // delete SmearBase...
@@ -37,32 +39,28 @@
     // Smear the configurations
     SmearBase->smear(C, U);
 
-  				for (int mu = 0; mu<Nd; mu++)
+    for (int mu = 0; mu < Nd; mu++) {
-  				{
       tmp = peekLorentz(C, mu);
       Umu = peekLorentz(U, mu);
-		  			iq_mu = Ta(tmp * adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper
+      iq_mu = Ta(
+          tmp *
+          adj(Umu));  // iq_mu = Ta(Omega_mu) to match the signs with the paper
       exponentiate_iQ(tmp, iq_mu);
       pokeLorentz(u_smr, tmp * Umu, mu);  // u_smr = exp(iQ_mu)*U_mu
     }
     std::cout << GridLogDebug << "Stout smearing completed\n";
   };
 
-
+  void derivative(GaugeField& SigmaTerm, const GaugeField& iLambda,
-			void derivative(GaugeField& SigmaTerm,
-				const GaugeField& iLambda,
                   const GaugeField& Gauge) const {
     SmearBase->derivative(SigmaTerm, iLambda, Gauge);
   };
 
-
+  void BaseSmear(GaugeField& C, const GaugeField& U) const {
-			void BaseSmear(GaugeField& C,
-				const GaugeField& U) const{
     SmearBase->smear(C, U);
   };
 
-			void exponentiate_iQ(GaugeLinkField& e_iQ,
+  void exponentiate_iQ(GaugeLinkField& e_iQ, const GaugeLinkField& iQ) const {
-				const GaugeLinkField& iQ) const{
     // Put this outside
     // only valid for SU(3) matrices
 
@@ -88,13 +86,10 @@
     set_fj(f0, f1, f2, u, w);
 
     e_iQ = f0 * unity + timesMinusI(f1) * iQ - f2 * iQ2;
-
-
   };
 
-
+  void set_uw(LatticeComplex& u, LatticeComplex& w, GaugeLinkField& iQ2,
-			void set_uw(LatticeComplex& u, LatticeComplex& w,
+              GaugeLinkField& iQ3) const {
-				GaugeLinkField& iQ2, GaugeLinkField& iQ3) const{
     Complex one_over_three = 1.0 / 3.0;
     Complex one_over_two = 1.0 / 2.0;
 
@@ -102,21 +97,21 @@
     LatticeComplex c0(grid), c1(grid), tmp(grid), c0max(grid), theta(grid);
 
     // sign in c0 from the conventions on the Ta
-				c0    = - real(timesMinusI(trace(iQ3))) * one_over_three; //temporary hack 
+    c0 = -imag(trace(iQ3)) * one_over_three;  
     c1 = -real(trace(iQ2)) * one_over_two;
 
     // Cayley Hamilton checks to machine precision, tested
     tmp = c1 * one_over_three;
     c0max = 2.0 * pow(tmp, 1.5);
 
-				theta = acos(c0/c0max)*one_over_three; // divide by three here, now leave as it is
+    theta = acos(c0 / c0max) *
+            one_over_three;  // divide by three here, now leave as it is
     u = sqrt(tmp) * cos(theta);
     w = sqrt(c1) * sin(theta);
   }
 
   void set_fj(LatticeComplex& f0, LatticeComplex& f1, LatticeComplex& f2,
               const LatticeComplex& u, const LatticeComplex& w) const {
-
     GridBase* grid = u._grid;
     LatticeComplex xi0(grid), u2(grid), w2(grid), cosw(grid);
     LatticeComplex fden(grid);
@@ -134,7 +129,8 @@
     emiu = cos(u) - timesI(sin(u));
     e2iu = cos(2.0 * u) + timesI(sin(2.0 * u));
 
-				h0    = e2iu * (u2 - w2) + emiu * ( (8.0*u2*cosw) + (2.0*u*(3.0*u2 + w2)*ixi0));
+    h0 = e2iu * (u2 - w2) +
+         emiu * ((8.0 * u2 * cosw) + (2.0 * u * (3.0 * u2 + w2) * ixi0));
     h1 = e2iu * (2.0 * u) - emiu * ((2.0 * u * cosw) - (3.0 * u2 - w2) * ixi0);
     h2 = e2iu - emiu * (cosw + (3.0 * u) * ixi0);
 
@@ -144,23 +140,20 @@
     f2 = h2 * fden;
   }
 
-
-
-
   LatticeComplex func_xi0(const LatticeComplex& w) const {
     // Define a function to do the check
-	//if( w < 1e-4 ) std::cout << GridLogWarning<< "[Smear_stout] w too small: "<< w <<"\n";
+    // if( w < 1e-4 ) std::cout << GridLogWarning<< "[Smear_stout] w too small:
+    // "<< w <<"\n";
     return sin(w) / w;
   }
 
   LatticeComplex func_xi1(const LatticeComplex& w) const {
     // Define a function to do the check
-	//if( w < 1e-4 ) std::cout << GridLogWarning << "[Smear_stout] w too small: "<< w <<"\n";
+    // if( w < 1e-4 ) std::cout << GridLogWarning << "[Smear_stout] w too small:
+    // "<< w <<"\n";
     return cos(w) / (w * w) - sin(w) / (w * w * w);
   }
-
 };
-
 }
 }
 

lib/simd/Grid_vector_types.h

@@ -25,7 +25,8 @@ Author: neo <cossu@post.kek.jp>
 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
+See the full license in the file "LICENSE" in the top level distribution
+directory
 *************************************************************************************/
 /*  END LEGAL */
 //---------------------------------------------------------------------------
@@ -64,10 +65,12 @@ namespace Grid {
 //////////////////////////////////////
 // To take the floating point type of real/complex type
 //////////////////////////////////////
-  template <typename T> struct RealPart {
+template <typename T>
+struct RealPart {
   typedef T type;
 };
-  template <typename T> struct RealPart< std::complex<T> >{
+template <typename T>
+struct RealPart<std::complex<T> > {
   typedef T type;
 };
 
@@ -76,24 +79,36 @@ namespace Grid {
 //////////////////////////////////////
 
 // type alias used to simplify the syntax of std::enable_if
-  template <typename T> using Invoke                                  =  typename T::type;
+template <typename T>
-  template <typename Condition, typename ReturnType> using EnableIf   =  Invoke<std::enable_if<Condition::value, ReturnType> >;
+using Invoke = typename T::type;
-  template <typename Condition, typename ReturnType> using NotEnableIf=  Invoke<std::enable_if<!Condition::value, ReturnType> >;
+template <typename Condition, typename ReturnType>
-
+using EnableIf = Invoke<std::enable_if<Condition::value, ReturnType> >;
+template <typename Condition, typename ReturnType>
+using NotEnableIf = Invoke<std::enable_if<!Condition::value, ReturnType> >;
 
 ////////////////////////////////////////////////////////
 // Check for complexity with type traits
-  template <typename T>   struct is_complex : public std::false_type {};
+template <typename T>
-  template <> struct is_complex<std::complex<double> >: public std::true_type {};
+struct is_complex : public std::false_type {};
-  template <> struct is_complex<std::complex<float> > : public std::true_type {};
+template <>
+struct is_complex<std::complex<double> > : public std::true_type {};
+template <>
+struct is_complex<std::complex<float> > : public std::true_type {};
 
-  template <typename T> using IfReal    = Invoke<std::enable_if<std::is_floating_point<T>::value,int> > ;
+template <typename T>
-  template <typename T> using IfComplex = Invoke<std::enable_if<is_complex<T>::value,int> > ;
+using IfReal = Invoke<std::enable_if<std::is_floating_point<T>::value, int> >;
-  template <typename T> using IfInteger = Invoke<std::enable_if<std::is_integral<T>::value,int> > ;
+template <typename T>
+using IfComplex = Invoke<std::enable_if<is_complex<T>::value, int> >;
+template <typename T>
+using IfInteger = Invoke<std::enable_if<std::is_integral<T>::value, int> >;
 
-  template <typename T> using IfNotReal    = Invoke<std::enable_if<!std::is_floating_point<T>::value,int> > ;
+template <typename T>
-  template <typename T> using IfNotComplex = Invoke<std::enable_if<!is_complex<T>::value,int> > ;
+using IfNotReal =
-  template <typename T> using IfNotInteger = Invoke<std::enable_if<!std::is_integral<T>::value,int> > ;
+    Invoke<std::enable_if<!std::is_floating_point<T>::value, int> >;
+template <typename T>
+using IfNotComplex = Invoke<std::enable_if<!is_complex<T>::value, int> >;
+template <typename T>
+using IfNotInteger = Invoke<std::enable_if<!std::is_integral<T>::value, int> >;
 
 ////////////////////////////////////////////////////////
 // Define the operation templates functors
@@ -111,33 +126,36 @@ namespace Grid {
 }
 ///////////////////////////////////////////////
 
-
-
 /*
   @brief Grid_simd class for the SIMD vector type operations
  */
 template <class Scalar_type, class Vector_type>
 class Grid_simd {
-    
  public:
   typedef typename RealPart<Scalar_type>::type Real;
   typedef Vector_type vector_type;
   typedef Scalar_type scalar_type;
 
-
   typedef union conv_t_union {
     Vector_type v;
     Scalar_type s[sizeof(Vector_type) / sizeof(Scalar_type)];
     conv_t_union(){};
   } conv_t;
 
-   
   Vector_type v;
 
-    static inline int Nsimd(void) { return sizeof(Vector_type)/sizeof(Scalar_type);}
+  static inline int Nsimd(void) {
+    return sizeof(Vector_type) / sizeof(Scalar_type);
+  }
 
-    Grid_simd& operator=(const Grid_simd&& rhs){v=rhs.v;return *this;};
+  Grid_simd &operator=(const Grid_simd &&rhs) {
-    Grid_simd& operator=(const Grid_simd& rhs){v=rhs.v;return *this;}; //faster than not declaring it and leaving to the compiler
+    v = rhs.v;
+    return *this;
+  };
+  Grid_simd &operator=(const Grid_simd &rhs) {
+    v = rhs.v;
+    return *this;
+  };  // faster than not declaring it and leaving to the compiler
   Grid_simd() = default;
   Grid_simd(const Grid_simd &rhs) : v(rhs.v){};  // compiles in movaps
   Grid_simd(const Grid_simd &&rhs) : v(rhs.v){};
@@ -156,31 +174,76 @@ namespace Grid {
     vsplat(*this, a);
   };
 
-    Grid_simd(const Real a){
+  Grid_simd(const Real a) { vsplat(*this, Scalar_type(a)); };
-      vsplat(*this,Scalar_type(a));
-    };
 
   ///////////////////////////////////////////////
   // mac, mult, sub, add, adj
   ///////////////////////////////////////////////
 
   // FIXME -- alias this to an inline MAC struct.
-    friend inline void mac (Grid_simd * __restrict__ y,const Grid_simd * __restrict__ a,const Grid_simd *__restrict__ x){ *y = (*a)*(*x)+(*y); };
+  friend inline void mac(Grid_simd *__restrict__ y,
+                         const Grid_simd *__restrict__ a,
+                         const Grid_simd *__restrict__ x) {
+    *y = (*a) * (*x) + (*y);
+  };
 
+  friend inline void mult(Grid_simd *__restrict__ y,
+                          const Grid_simd *__restrict__ l,
+                          const Grid_simd *__restrict__ r) {
+    *y = (*l) * (*r);
+  }
+  friend inline void sub(Grid_simd *__restrict__ y,
+                         const Grid_simd *__restrict__ l,
+                         const Grid_simd *__restrict__ r) {
+    *y = (*l) - (*r);
+  }
+  friend inline void add(Grid_simd *__restrict__ y,
+                         const Grid_simd *__restrict__ l,
+                         const Grid_simd *__restrict__ r) {
+    *y = (*l) + (*r);
+  }
 
-    friend inline void mult(Grid_simd * __restrict__ y,const Grid_simd * __restrict__ l,const Grid_simd *__restrict__ r){ *y = (*l) * (*r); }
+  friend inline void mac(Grid_simd *__restrict__ y,
-    friend inline void sub (Grid_simd * __restrict__ y,const Grid_simd * __restrict__ l,const Grid_simd *__restrict__ r){ *y = (*l) - (*r); }
+                         const Scalar_type *__restrict__ a,
-    friend inline void add (Grid_simd * __restrict__ y,const Grid_simd * __restrict__ l,const Grid_simd *__restrict__ r){ *y = (*l) + (*r); }
+                         const Grid_simd *__restrict__ x) {
+    *y = (*a) * (*x) + (*y);
+  };
+  friend inline void mult(Grid_simd *__restrict__ y,
+                          const Scalar_type *__restrict__ l,
+                          const Grid_simd *__restrict__ r) {
+    *y = (*l) * (*r);
+  }
+  friend inline void sub(Grid_simd *__restrict__ y,
+                         const Scalar_type *__restrict__ l,
+                         const Grid_simd *__restrict__ r) {
+    *y = (*l) - (*r);
+  }
+  friend inline void add(Grid_simd *__restrict__ y,
+                         const Scalar_type *__restrict__ l,
+                         const Grid_simd *__restrict__ r) {
+    *y = (*l) + (*r);
+  }
 
-    friend inline void mac (Grid_simd *__restrict__ y,const Scalar_type *__restrict__ a,const Grid_simd   *__restrict__ x){ *y = (*a)*(*x)+(*y); };
+  friend inline void mac(Grid_simd *__restrict__ y,
-    friend inline void mult(Grid_simd *__restrict__ y,const Scalar_type *__restrict__ l,const Grid_simd   *__restrict__ r){ *y = (*l) * (*r); }
+                         const Grid_simd *__restrict__ a,
-    friend inline void sub (Grid_simd *__restrict__ y,const Scalar_type *__restrict__ l,const Grid_simd   *__restrict__ r){ *y = (*l) - (*r); }
+                         const Scalar_type *__restrict__ x) {
-    friend inline void add (Grid_simd *__restrict__ y,const Scalar_type *__restrict__ l,const Grid_simd   *__restrict__ r){ *y = (*l) + (*r); }
+    *y = (*a) * (*x) + (*y);
-
+  };
-    friend inline void mac (Grid_simd *__restrict__ y,const Grid_simd   *__restrict__ a,const Scalar_type *__restrict__ x){ *y = (*a)*(*x)+(*y); };
+  friend inline void mult(Grid_simd *__restrict__ y,
-    friend inline void mult(Grid_simd *__restrict__ y,const Grid_simd   *__restrict__ l,const Scalar_type *__restrict__ r){ *y = (*l) * (*r); }
+                          const Grid_simd *__restrict__ l,
-    friend inline void sub (Grid_simd *__restrict__ y,const Grid_simd   *__restrict__ l,const Scalar_type *__restrict__ r){ *y = (*l) - (*r); }
+                          const Scalar_type *__restrict__ r) {
-    friend inline void add (Grid_simd *__restrict__ y,const Grid_simd   *__restrict__ l,const Scalar_type *__restrict__ r){ *y = (*l) + (*r); }
+    *y = (*l) * (*r);
+  }
+  friend inline void sub(Grid_simd *__restrict__ y,
+                         const Grid_simd *__restrict__ l,
+                         const Scalar_type *__restrict__ r) {
+    *y = (*l) - (*r);
+  }
+  friend inline void add(Grid_simd *__restrict__ y,
+                         const Grid_simd *__restrict__ l,
+                         const Scalar_type *__restrict__ r) {
+    *y = (*l) + (*r);
+  }
 
   ////////////////////////////////////////////////////////////////////////
   // FIXME:  gonna remove these load/store, get, set, prefetch
@@ -199,16 +262,14 @@ namespace Grid {
   ///////////////////////
   // Vprefetch
   ///////////////////////
-    friend inline void vprefetch(const Grid_simd &v)
+  friend inline void vprefetch(const Grid_simd &v) {
-    {
     prefetch_HINT_T0((const char *)&v.v);
   }
 
   ///////////////////////
   // Reduce
   ///////////////////////
-    friend inline Scalar_type Reduce(const Grid_simd & in)
+  friend inline Scalar_type Reduce(const Grid_simd &in) {
-    {
     return unary<Scalar_type>(in.v, ReduceSIMD<Scalar_type, Vector_type>());
   }
 
@@ -255,7 +316,8 @@ namespace Grid {
   // provides support
   ///////////////////////////////////////
 
-    template<class functor> friend inline Grid_simd SimdApply (const functor &func,const Grid_simd &v) {
+  template <class functor>
+  friend inline Grid_simd SimdApply(const functor &func, const Grid_simd &v) {
     Grid_simd ret;
     Grid_simd::conv_t conv;
 
@@ -266,7 +328,10 @@ namespace Grid {
     ret.v = conv.v;
     return ret;
   }
-    template<class functor> friend inline Grid_simd SimdApplyBinop (const functor &func,const Grid_simd &x,const Grid_simd &y) {
+  template <class functor>
+  friend inline Grid_simd SimdApplyBinop(const functor &func,
+                                         const Grid_simd &x,
+                                         const Grid_simd &y) {
     Grid_simd ret;
     Grid_simd::conv_t cx;
     Grid_simd::conv_t cy;
@@ -297,19 +362,27 @@ namespace Grid {
   friend inline void permute3(Grid_simd &y, Grid_simd b) {
     y.v = Optimization::Permute::Permute3(b.v);
   }
-    friend inline void permute(Grid_simd &y,Grid_simd b,int perm)
+  friend inline void permute(Grid_simd &y, Grid_simd b, int perm) {
-    {
     if (perm & RotateBit) {
       int dist = perm & 0xF;
       y = rotate(b, dist);
       return;
     }
     switch (perm) {
-      case 3: permute3(y,b); break;
+      case 3:
-      case 2: permute2(y,b); break;
+        permute3(y, b);
-      case 1: permute1(y,b); break;
+        break;
-      case 0: permute0(y,b); break;
+      case 2:
-      default: assert(0);
+        permute2(y, b);
+        break;
+      case 1:
+        permute1(y, b);
+        break;
+      case 0:
+        permute0(y, b);
+        break;
+      default:
+        assert(0);
     }
   }
 
@@ -319,8 +392,7 @@ namespace Grid {
 // General rotate
 ////////////////////////////////////////////////////////////////////
 template <class S, class V, IfNotComplex<S> = 0>
-  inline Grid_simd<S,V> rotate(Grid_simd<S,V> b,int nrot)
+inline Grid_simd<S, V> rotate(Grid_simd<S, V> b, int nrot) {
-  {
   nrot = nrot % Grid_simd<S, V>::Nsimd();
   Grid_simd<S, V> ret;
   //    std::cout << "Rotate Real by "<<nrot<<std::endl;
@@ -328,8 +400,7 @@ namespace Grid {
   return ret;
 }
 template <class S, class V, IfComplex<S> = 0>
-  inline Grid_simd<S,V> rotate(Grid_simd<S,V> b,int nrot)
+inline Grid_simd<S, V> rotate(Grid_simd<S, V> b, int nrot) {
-  {
   nrot = nrot % Grid_simd<S, V>::Nsimd();
   Grid_simd<S, V> ret;
   //    std::cout << "Rotate Complex by "<<nrot<<std::endl;
@@ -348,11 +419,13 @@ namespace Grid {
 }
 
 // overload if complex
-  template <class S,class V> inline void vsplat(Grid_simd<S,V> &ret, EnableIf<is_complex < S >, S> c) {
+template <class S, class V>
+inline void vsplat(Grid_simd<S, V> &ret, EnableIf<is_complex<S>, S> c) {
   vsplat(ret, real(c), imag(c));
 }
 
-  //if real fill with a, if complex fill with a in the real part (first function above)
+// if real fill with a, if complex fill with a in the real part (first function
+// above)
 template <class S, class V>
 inline void vsplat(Grid_simd<S, V> &ret, NotEnableIf<is_complex<S>, S> a) {
   ret.v = unary<V>(a, VsplatSIMD());
@@ -363,24 +436,60 @@ namespace Grid {
 // Initialise to 1,0,i for the correct types
 ///////////////////////////////////////////////
 // For complex types
-  template <class S,class V, IfComplex<S> = 0 > inline void vone(Grid_simd<S,V>  &ret)     { vsplat(ret,S(1.0,0.0)); }
+template <class S, class V, IfComplex<S> = 0>
-  template <class S,class V, IfComplex<S> = 0 > inline void vzero(Grid_simd<S,V> &ret)     { vsplat(ret,S(0.0,0.0)); }// use xor?
+inline void vone(Grid_simd<S, V> &ret) {
-  template <class S,class V, IfComplex<S> = 0 > inline void vcomplex_i(Grid_simd<S,V> &ret){ vsplat(ret,S(0.0,1.0));} 
+  vsplat(ret, S(1.0, 0.0));
+}
+template <class S, class V, IfComplex<S> = 0>
+inline void vzero(Grid_simd<S, V> &ret) {
+  vsplat(ret, S(0.0, 0.0));
+}  // use xor?
+template <class S, class V, IfComplex<S> = 0>
+inline void vcomplex_i(Grid_simd<S, V> &ret) {
+  vsplat(ret, S(0.0, 1.0));
+}
 
-  template <class S,class V, IfComplex<S> = 0 > inline void visign(Grid_simd<S,V> &ret){ vsplat(ret,S(1.0,-1.0));} 
+template <class S, class V, IfComplex<S> = 0>
-  template <class S,class V, IfComplex<S> = 0 > inline void vrsign(Grid_simd<S,V> &ret){ vsplat(ret,S(-1.0,1.0));} 
+inline void visign(Grid_simd<S, V> &ret) {
+  vsplat(ret, S(1.0, -1.0));
+}
+template <class S, class V, IfComplex<S> = 0>
+inline void vrsign(Grid_simd<S, V> &ret) {
+  vsplat(ret, S(-1.0, 1.0));
+}
 
 // if not complex overload here
-  template <class S,class V, IfReal<S> = 0 > inline void vone (Grid_simd<S,V> &ret){ vsplat(ret,S(1.0)); }
+template <class S, class V, IfReal<S> = 0>
-  template <class S,class V, IfReal<S> = 0 > inline void vzero(Grid_simd<S,V> &ret){ vsplat(ret,S(0.0)); }
+inline void vone(Grid_simd<S, V> &ret) {
+  vsplat(ret, S(1.0));
+}
+template <class S, class V, IfReal<S> = 0>
+inline void vzero(Grid_simd<S, V> &ret) {
+  vsplat(ret, S(0.0));
+}
 
 // For integral types
-  template <class S,class V,IfInteger<S> = 0 > inline void vone(Grid_simd<S,V> &ret)  {vsplat(ret,1); }
+template <class S, class V, IfInteger<S> = 0>
-  template <class S,class V,IfInteger<S> = 0 > inline void vzero(Grid_simd<S,V> &ret) {vsplat(ret,0); }
+inline void vone(Grid_simd<S, V> &ret) {
-  template <class S,class V,IfInteger<S> = 0 > inline void vtrue (Grid_simd<S,V> &ret){vsplat(ret,0xFFFFFFFF);}
+  vsplat(ret, 1);
-  template <class S,class V,IfInteger<S> = 0 > inline void vfalse(Grid_simd<S,V> &ret){vsplat(ret,0);}
+}
+template <class S, class V, IfInteger<S> = 0>
+inline void vzero(Grid_simd<S, V> &ret) {
+  vsplat(ret, 0);
+}
+template <class S, class V, IfInteger<S> = 0>
+inline void vtrue(Grid_simd<S, V> &ret) {
+  vsplat(ret, 0xFFFFFFFF);
+}
+template <class S, class V, IfInteger<S> = 0>
+inline void vfalse(Grid_simd<S, V> &ret) {
+  vsplat(ret, 0);
+}
 
-  template<class S,class V> inline void zeroit(Grid_simd<S,V> &z){ vzero(z);}
+template <class S, class V>
+inline void zeroit(Grid_simd<S, V> &z) {
+  vzero(z);
+}
 
 ///////////////////////
 // Vstream
@@ -403,33 +512,36 @@ namespace Grid {
 ////////////////////////////////////
 // Arithmetic operator overloads +,-,*
 ////////////////////////////////////
-  template<class S,class V> inline Grid_simd<S,V> operator + (Grid_simd<S,V> a, Grid_simd<S,V> b) {
+template <class S, class V>
+inline Grid_simd<S, V> operator+(Grid_simd<S, V> a, Grid_simd<S, V> b) {
   Grid_simd<S, V> ret;
   ret.v = binary<V>(a.v, b.v, SumSIMD());
   return ret;
 };
 
-  template<class S,class V> inline Grid_simd<S,V> operator - (Grid_simd<S,V> a, Grid_simd<S,V> b) {
+template <class S, class V>
+inline Grid_simd<S, V> operator-(Grid_simd<S, V> a, Grid_simd<S, V> b) {
   Grid_simd<S, V> ret;
   ret.v = binary<V>(a.v, b.v, SubSIMD());
   return ret;
 };
 
 // Distinguish between complex types and others
-  template<class S,class V, IfComplex<S> = 0 > inline Grid_simd<S,V> operator * (Grid_simd<S,V> a, Grid_simd<S,V> b) {
+template <class S, class V, IfComplex<S> = 0>
+inline Grid_simd<S, V> operator*(Grid_simd<S, V> a, Grid_simd<S, V> b) {
   Grid_simd<S, V> ret;
   ret.v = binary<V>(a.v, b.v, MultComplexSIMD());
   return ret;
 };
 
 // Real/Integer types
-  template<class S,class V, IfNotComplex<S> = 0 > inline Grid_simd<S,V> operator * (Grid_simd<S,V> a, Grid_simd<S,V> b) {
+template <class S, class V, IfNotComplex<S> = 0>
+inline Grid_simd<S, V> operator*(Grid_simd<S, V> a, Grid_simd<S, V> b) {
   Grid_simd<S, V> ret;
   ret.v = binary<V>(a.v, b.v, MultSIMD());
   return ret;
 };
 
-
 ///////////////////////
 // Conjugate
 ///////////////////////
@@ -439,13 +551,16 @@ namespace Grid {
   ret.v = unary<V>(in.v, ConjSIMD());
   return ret;
 }
-  template <class S,class V, IfNotComplex<S> = 0 > inline Grid_simd<S,V> conjugate(const Grid_simd<S,V>  &in){
+template <class S, class V, IfNotComplex<S> = 0>
+inline Grid_simd<S, V> conjugate(const Grid_simd<S, V> &in) {
   return in;  // for real objects
 }
 
 // Suppress adj for integer types... // odd; why conjugate above but not adj??
 template <class S, class V, IfNotInteger<S> = 0>
-    inline Grid_simd<S,V> adj(const Grid_simd<S,V> &in){ return conjugate(in); }
+inline Grid_simd<S, V> adj(const Grid_simd<S, V> &in) {
+  return conjugate(in);
+}
 
 ///////////////////////
 // timesMinusI
@@ -492,14 +607,14 @@ namespace Grid {
 /////////////////////
 
 template <class S, class V>
-    inline Grid_simd< S, V>  innerProduct(const Grid_simd< S, V> & l, const Grid_simd< S, V> & r) 
+inline Grid_simd<S, V> innerProduct(const Grid_simd<S, V> &l,
-  {
+                                    const Grid_simd<S, V> &r) {
   return conjugate(l) * r;
 }
 
 template <class S, class V>
-  inline Grid_simd< S, V> outerProduct(const Grid_simd< S, V> &l, const Grid_simd< S, V> & r)
+inline Grid_simd<S, V> outerProduct(const Grid_simd<S, V> &l,
-  {
+                                    const Grid_simd<S, V> &r) {
   return l * conjugate(r);
 }
 
@@ -508,7 +623,6 @@ namespace Grid {
   return arg;
 }
 
-
 ////////////////////////////////////////////////////////////
 // copy/splat complex real parts into real;
 // insert real into complex and zero imag;
@@ -516,8 +630,7 @@ namespace Grid {
 
 // real = toReal( complex )
 template <class S, class V, IfReal<S> = 0>
-  inline Grid_simd<S,V> toReal(const Grid_simd<std::complex<S>,V> &in)
+inline Grid_simd<S, V> toReal(const Grid_simd<std::complex<S>, V> &in) {
-  {
   typedef Grid_simd<S, V> simd;
   simd ret;
   typename simd::conv_t conv;
@@ -531,16 +644,17 @@ namespace Grid {
 
 // complex = toComplex( real )
 template <class R, class V, IfReal<R> = 0>  // must be a real arg
-  inline Grid_simd<std::complex<R>,V> toComplex (const Grid_simd<R,V> &in)
+inline Grid_simd<std::complex<R>, V> toComplex(const Grid_simd<R, V> &in) {
-  {
   typedef Grid_simd<R, V> Rsimd;
   typedef Grid_simd<std::complex<R>, V> Csimd;
   typename Rsimd::conv_t conv;  // address as real
 
   conv.v = in.v;
   for (int i = 0; i < Rsimd::Nsimd(); i += 2) {
-      assert(conv.s[i+1]==conv.s[i]); // trap any cases where real was not duplicated 
+    assert(conv.s[i + 1] ==
-      // indicating the SIMD grids of real and imag assignment did not correctly match
+           conv.s[i]);  // trap any cases where real was not duplicated
+    // indicating the SIMD grids of real and imag assignment did not correctly
+    // match
     conv.s[i + 1] = 0.0;  // zero imaginary parts
   }
   Csimd ret;
@@ -548,9 +662,6 @@ namespace Grid {
   return ret;
 }
 
-
-
-  
 ///////////////////////////////
 // Define available types
 ///////////////////////////////
@@ -563,16 +674,23 @@ namespace Grid {
 /////////////////////////////////////////
 // Some traits to recognise the types
 /////////////////////////////////////////
-  template <typename T> struct is_simd : public std::false_type{};
+template <typename T>
-  template <> struct is_simd<vRealF>   : public std::true_type {};
+struct is_simd : public std::false_type {};
-  template <> struct is_simd<vRealD>   : public std::true_type {};
+template <>
-  template <> struct is_simd<vComplexF>: public std::true_type {};
+struct is_simd<vRealF> : public std::true_type {};
-  template <> struct is_simd<vComplexD>: public std::true_type {};
+template <>
-  template <> struct is_simd<vInteger> : public std::true_type {};
+struct is_simd<vRealD> : public std::true_type {};
-
+template <>
-  template <typename T> using IfSimd     = Invoke<std::enable_if< is_simd<T>::value,int> > ;
+struct is_simd<vComplexF> : public std::true_type {};
-  template <typename T> using IfNotSimd  = Invoke<std::enable_if<!is_simd<T>::value,unsigned> > ;
+template <>
+struct is_simd<vComplexD> : public std::true_type {};
+template <>
+struct is_simd<vInteger> : public std::true_type {};
 
+template <typename T>
+using IfSimd = Invoke<std::enable_if<is_simd<T>::value, int> >;
+template <typename T>
+using IfNotSimd = Invoke<std::enable_if<!is_simd<T>::value, unsigned> >;
 }
 
 #endif

lib/simd/Grid_vector_unops.h

@@ -25,7 +25,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 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
+See the full license in the file "LICENSE" in the top level distribution
+directory
 *************************************************************************************/
 /*  END LEGAL */
 #ifndef GRID_VECTOR_UNOPS
@@ -35,97 +36,84 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 namespace Grid {
 
-  template<class scalar> struct SqrtRealFunctor {
+template <class scalar>
-    scalar operator()(const scalar &a) const {
+struct SqrtRealFunctor {
-      return sqrt(real(a));
+  scalar operator()(const scalar &a) const { return sqrt(real(a)); }
-    }
 };
 
-  template<class scalar> struct RSqrtRealFunctor {
+template <class scalar>
+struct RSqrtRealFunctor {
   scalar operator()(const scalar &a) const {
     return scalar(1.0 / sqrt(real(a)));
   }
 };
 
-  template<class scalar> struct CosRealFunctor {
+template <class scalar>
-    scalar operator()(const scalar &a)  const {
+struct CosRealFunctor {
-      return cos(real(a));
+  scalar operator()(const scalar &a) const { return cos(real(a)); }
-    }
 };
 
-  template<class scalar> struct SinRealFunctor {
+template <class scalar>
-    scalar operator()(const scalar &a)  const {
+struct SinRealFunctor {
-      return sin(real(a));
+  scalar operator()(const scalar &a) const { return sin(real(a)); }
-    }
 };
 
-  template<class scalar> struct AcosRealFunctor {
+template <class scalar>
-    scalar operator()(const scalar &a)  const {
+struct AcosRealFunctor {
-      return acos(real(a));
+  scalar operator()(const scalar &a) const { return acos(real(a)); }
-    }
 };
 
-  template<class scalar> struct AsinRealFunctor {
+template <class scalar>
-    scalar operator()(const scalar &a)  const {
+struct AsinRealFunctor {
-      return asin(real(a));
+  scalar operator()(const scalar &a) const { return asin(real(a)); }
-    }
 };
 
-  template<class scalar> struct LogRealFunctor {
+template <class scalar>
-    scalar operator()(const scalar &a)  const {
+struct LogRealFunctor {
-      return log(real(a));
+  scalar operator()(const scalar &a) const { return log(real(a)); }
-    }
 };
 
-  template<class scalar> struct ExpRealFunctor {
+template <class scalar>
-    scalar operator()(const scalar &a)  const {
+struct ExpRealFunctor {
-      return exp(real(a));
+  scalar operator()(const scalar &a) const { return exp(real(a)); }
-    }
 };
-  template<class scalar> struct NotFunctor {
+template <class scalar>
-    scalar operator()(const scalar &a)  const {
+struct NotFunctor {
-      return (!a);
+  scalar operator()(const scalar &a) const { return (!a); }
-    }
 };
-  template<class scalar> struct AbsRealFunctor {
+template <class scalar>
-    scalar operator()(const scalar &a)  const {
+struct AbsRealFunctor {
-      return std::abs(real(a));
+  scalar operator()(const scalar &a) const { return std::abs(real(a)); }
-    }
 };
 
-  template<class scalar> struct PowRealFunctor {
+template <class scalar>
+struct PowRealFunctor {
   double y;
   PowRealFunctor(double _y) : y(_y){};
-    scalar operator()(const scalar &a)  const {
+  scalar operator()(const scalar &a) const { return pow(real(a), y); }
-      return pow(real(a),y);
-    }
 };
 
-  template<class scalar> struct ModIntFunctor {
+template <class scalar>
+struct ModIntFunctor {
   Integer y;
   ModIntFunctor(Integer _y) : y(_y){};
-    scalar operator()(const scalar &a)  const {
+  scalar operator()(const scalar &a) const { return Integer(a) % y; }
-      return Integer(a)%y;
-    }
 };
 
-  template<class scalar> struct DivIntFunctor {
+template <class scalar>
+struct DivIntFunctor {
   Integer y;
   DivIntFunctor(Integer _y) : y(_y){};
-    scalar operator()(const scalar &a)  const {
+  scalar operator()(const scalar &a) const { return Integer(a) / y; }
-      return Integer(a)/y;
-    }
 };
 
-  template<class scalar> struct RealFunctor {
+template <class scalar>
-    scalar operator()(const std::complex<scalar> &a)  const {
+struct RealFunctor {
-      return real(a);
+  scalar operator()(const scalar &a) const { return std::real(a); }
-    }
 };
-  template<class scalar> struct ImagFunctor {
+template <class scalar>
-    scalar operator()(const std::complex<scalar> &a)  const {
+struct ImagFunctor {
-      return imag(a);
+  scalar operator()(const scalar &a) const { return std::imag(a); }
-    }
 };
 template <class S, class V>
 inline Grid_simd<S, V> real(const Grid_simd<S, V> &r) {
@@ -135,7 +123,6 @@ namespace Grid {
 inline Grid_simd<S, V> imag(const Grid_simd<S, V> &r) {
   return SimdApply(ImagFunctor<S>(), r);
 }
-
 template <class S, class V>
 inline Grid_simd<S, V> sqrt(const Grid_simd<S, V> &r) {
   return SimdApply(SqrtRealFunctor<S>(), r);
@@ -197,51 +184,51 @@ namespace Grid {
 // Allows us to assign into **conformable** real vectors from complex
 ////////////////////////////////////////////////////////////////////////////
 //  template < class S, class V >
-  //  inline auto ComplexRemove(const Grid_simd<S,V> &c) -> Grid_simd<Grid_simd<S,V>::Real,V> {
+//  inline auto ComplexRemove(const Grid_simd<S,V> &c) ->
+//  Grid_simd<Grid_simd<S,V>::Real,V> {
 //    Grid_simd<Grid_simd<S,V>::Real,V> ret;
 //    ret.v = c.v;
 //    return ret;
 //  }
-  template<class scalar> struct AndFunctor {
+template <class scalar>
-    scalar operator()(const scalar &x, const scalar &y)  const {
+struct AndFunctor {
-      return x & y;
+  scalar operator()(const scalar &x, const scalar &y) const { return x & y; }
-    }
 };
-  template<class scalar> struct OrFunctor {
+template <class scalar>
-    scalar operator()(const scalar &x, const scalar &y)  const {
+struct OrFunctor {
-      return x | y;
+  scalar operator()(const scalar &x, const scalar &y) const { return x | y; }
-    }
 };
-  template<class scalar> struct AndAndFunctor {
+template <class scalar>
-    scalar operator()(const scalar &x, const scalar &y)  const {
+struct AndAndFunctor {
-      return x && y;
+  scalar operator()(const scalar &x, const scalar &y) const { return x && y; }
-    }
 };
-  template<class scalar> struct OrOrFunctor {
+template <class scalar>
-    scalar operator()(const scalar &x, const scalar &y)  const {
+struct OrOrFunctor {
-      return x || y;
+  scalar operator()(const scalar &x, const scalar &y) const { return x || y; }
-    }
 };
 
 ////////////////////////////////
 // Calls to simd binop functors
 ////////////////////////////////
 template <class S, class V>
-  inline Grid_simd<S,V> operator &(const Grid_simd<S,V> &x,const Grid_simd<S,V> &y) {
+inline Grid_simd<S, V> operator&(const Grid_simd<S, V> &x,
+                                 const Grid_simd<S, V> &y) {
   return SimdApplyBinop(AndFunctor<S>(), x, y);
 }
 template <class S, class V>
-  inline Grid_simd<S,V> operator &&(const Grid_simd<S,V> &x,const Grid_simd<S,V> &y) {
+inline Grid_simd<S, V> operator&&(const Grid_simd<S, V> &x,
+                                  const Grid_simd<S, V> &y) {
   return SimdApplyBinop(AndAndFunctor<S>(), x, y);
 }
 template <class S, class V>
-  inline Grid_simd<S,V> operator |(const Grid_simd<S,V> &x,const Grid_simd<S,V> &y) {
+inline Grid_simd<S, V> operator|(const Grid_simd<S, V> &x,
+                                 const Grid_simd<S, V> &y) {
   return SimdApplyBinop(OrFunctor<S>(), x, y);
 }
 template <class S, class V>
-  inline Grid_simd<S,V> operator ||(const Grid_simd<S,V> &x,const Grid_simd<S,V> &y) {
+inline Grid_simd<S, V> operator||(const Grid_simd<S, V> &x,
+                                  const Grid_simd<S, V> &y) {
   return SimdApplyBinop(OrOrFunctor<S>(), x, y);
 }
-
 }
 #endif

tests/Test_simd.cc

@@ -23,7 +23,8 @@ Author: neo <cossu@post.kek.jp>
 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
+See the full license in the file "LICENSE" in the top level distribution
+directory
 *************************************************************************************/
 /*  END LEGAL */
 #include <Grid.h>
@@ -62,6 +63,18 @@ public:
   template<class vec> void operator()(vec &rr,vec &i1,vec &i2) const { rr = adj(i1);}
   std::string name(void) const { return std::string("Adj"); }
 };
+class funcImag {
+public:
+  funcImag() {};
+  template<class vec> void operator()(vec &rr,vec &i1,vec &i2) const { rr = imag(i1);}
+  std::string name(void) const { return std::string("imag"); }
+};
+class funcReal {
+public:
+  funcReal() {};
+  template<class vec> void operator()(vec &rr,vec &i1,vec &i2) const { rr = real(i1);}
+  std::string name(void) const { return std::string("real"); }
+};
 
 class funcTimesI {
 public:
@@ -141,8 +154,14 @@ void Tester(const functor &func)
   }
 
   extract<vec,scal>(v_result,result);
+
   std::cout << GridLogMessage << " " << func.name() << std::endl;
 
+  std::cout << GridLogDebug << v_input1 << std::endl;
+  std::cout << GridLogDebug << v_result << std::endl;
+
+
+
   int ok=0;
   for(int i=0;i<Nsimd;i++){
     if ( abs(reference[i]-result[i])>1.0e-7){
@@ -389,6 +408,8 @@ int main (int argc, char ** argv)
   Tester<ComplexF,vComplexF>(funcTimes());
   Tester<ComplexF,vComplexF>(funcConj());
   Tester<ComplexF,vComplexF>(funcAdj());
+  Tester<ComplexF,vComplexF>(funcReal());
+  Tester<ComplexF,vComplexF>(funcImag());
   Tester<ComplexF,vComplexF>(funcInnerProduct());
   ReductionTester<ComplexF,ComplexF,vComplexF>(funcReduce());
 
@@ -421,13 +442,17 @@ int main (int argc, char ** argv)
   Tester<ComplexD,vComplexD>(funcTimes());
   Tester<ComplexD,vComplexD>(funcConj());
   Tester<ComplexD,vComplexD>(funcAdj());
+  Tester<ComplexD, vComplexD>(funcReal());
+  Tester<ComplexD, vComplexD>(funcImag());
+
   Tester<ComplexD, vComplexD>(funcInnerProduct());
   ReductionTester<ComplexD, ComplexD, vComplexD>(funcReduce());
 
-
+  std::cout << GridLogMessage
-  std::cout<<GridLogMessage << "==================================="<<  std::endl;
+            << "===================================" << std::endl;
   std::cout << GridLogMessage << "Testing vComplexD permutes " << std::endl;
-  std::cout<<GridLogMessage << "==================================="<<  std::endl;
+  std::cout << GridLogMessage
+            << "===================================" << std::endl;
 
   // Log2 iteration
   for (int i = 0; (1 << i) < vComplexD::Nsimd(); i++) {