From bbb657da5c31fdbcfbb2b695533305a85039ed90 Mon Sep 17 00:00:00 2001
From: paboyle <paboyle@ph.ed.ac.uk>
Date: Fri, 12 Jan 2018 18:10:11 +0000
Subject: [PATCH] NAMESPACE and formatting

---
 lib/simd/Grid_generic.h | 724 ++++++++++++++++++++--------------------
 1 file changed, 363 insertions(+), 361 deletions(-)

diff --git a/lib/simd/Grid_generic.h b/lib/simd/Grid_generic.h
index e1d5f894..d555a672 100644
--- a/lib/simd/Grid_generic.h
+++ b/lib/simd/Grid_generic.h
@@ -1,4 +1,4 @@
-    /*************************************************************************************
+/*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
 
@@ -25,413 +25,413 @@ Author: Antonin Portelli <antonin.portelli@me.com>
     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
     See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
+*************************************************************************************/
+/*  END LEGAL */
 
 #include "Grid_generic_types.h"
 
-namespace Grid {
-namespace Optimization {
+NAMESPACE_BEGIN(Grid);
+NAMESPACE_BEGIN(Optimization);
   
-  struct Vsplat{
-    // Complex
-    template <typename T>
-    inline vec<T> operator()(T a, T b){
-      vec<T> out;
+struct Vsplat{
+  // Complex
+  template <typename T>
+  inline vec<T> operator()(T a, T b){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::r, 2)
+    VECTOR_FOR(i, W<T>::r, 2)
       {
         out.v[i]   = a;
         out.v[i+1] = b;
       }
 
-      return out;
-    }
+    return out;
+  }
     
-    // Real
-    template <typename T>
-    inline vec<T> operator()(T a){
-      vec<T> out;
+  // Real
+  template <typename T>
+  inline vec<T> operator()(T a){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::r, 1)
+    VECTOR_FOR(i, W<T>::r, 1)
       {
         out.v[i] = a;
       }
       
-      return out;
-    }
-  };
+    return out;
+  }
+};
 
-  struct Vstore{
-    // Real
-    template <typename T>
-    inline void operator()(vec<T> a, T *D){
-      *((vec<T> *)D) = a;
-    }
-  };
+struct Vstore{
+  // Real
+  template <typename T>
+  inline void operator()(vec<T> a, T *D){
+    *((vec<T> *)D) = a;
+  }
+};
 
-  struct Vstream{
-    // Real
-    template <typename T>
-    inline void operator()(T * a, vec<T> b){
-      *((vec<T> *)a) = b;
-    }
-  };
+struct Vstream{
+  // Real
+  template <typename T>
+  inline void operator()(T * a, vec<T> b){
+    *((vec<T> *)a) = b;
+  }
+};
 
-  struct Vset{
-    // Complex
-    template <typename T>
-    inline vec<T> operator()(std::complex<T> *a){
-      vec<T> out;
+struct Vset{
+  // Complex
+  template <typename T>
+  inline vec<T> operator()(std::complex<T> *a){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::c, 1)
+    VECTOR_FOR(i, W<T>::c, 1)
       {
         out.v[2*i]   = a[i].real();
         out.v[2*i+1] = a[i].imag();
       }
       
-      return out;
-    }
+    return out;
+  }
     
-    // Real
-    template <typename T>
-    inline vec<T> operator()(T *a){
-      vec<T> out;
+  // Real
+  template <typename T>
+  inline vec<T> operator()(T *a){
+    vec<T> out;
       
-      out = *((vec<T> *)a);
+    out = *((vec<T> *)a);
       
-      return out;
-    }
-  };
+    return out;
+  }
+};
 
-  /////////////////////////////////////////////////////
-  // Arithmetic operations
-  /////////////////////////////////////////////////////
-  struct Sum{
-    // Complex/Real
-    template <typename T>
-    inline vec<T> operator()(vec<T> a, vec<T> b){
-      vec<T> out;
+/////////////////////////////////////////////////////
+// Arithmetic operations
+/////////////////////////////////////////////////////
+struct Sum{
+  // Complex/Real
+  template <typename T>
+  inline vec<T> operator()(vec<T> a, vec<T> b){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::r, 1)
+    VECTOR_FOR(i, W<T>::r, 1)
       {
         out.v[i] = a.v[i] + b.v[i];
       }
       
-      return out;
-    }
-  };
+    return out;
+  }
+};
 
-  struct Sub{
-    // Complex/Real
-    template <typename T>
-    inline vec<T> operator()(vec<T> a, vec<T> b){
-      vec<T> out;
+struct Sub{
+  // Complex/Real
+  template <typename T>
+  inline vec<T> operator()(vec<T> a, vec<T> b){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::r, 1)
+    VECTOR_FOR(i, W<T>::r, 1)
       {
         out.v[i] = a.v[i] - b.v[i];
       }
       
-      return out;
-    }
-  };
+    return out;
+  }
+};
 
-  struct Mult{
-    // Real
-    template <typename T>
-    inline vec<T> operator()(vec<T> a, vec<T> b){
-      vec<T> out;
+struct Mult{
+  // Real
+  template <typename T>
+  inline vec<T> operator()(vec<T> a, vec<T> b){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::r, 1)
+    VECTOR_FOR(i, W<T>::r, 1)
       {
         out.v[i] = a.v[i]*b.v[i];
       }
       
-      return out;
-    }
-  };
+    return out;
+  }
+};
   
-  #define cmul(a, b, c, i)\
-  c[i]   = a[i]*b[i]   - a[i+1]*b[i+1];\
+#define cmul(a, b, c, i)			\
+  c[i]   = a[i]*b[i]   - a[i+1]*b[i+1];		\
   c[i+1] = a[i]*b[i+1] + a[i+1]*b[i];
 
-  struct MultRealPart{
-    template <typename T>
-    inline vec<T> operator()(vec<T> a, vec<T> b){
-      vec<T> out;
+struct MultRealPart{
+  template <typename T>
+  inline vec<T> operator()(vec<T> a, vec<T> b){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::c, 1)
+    VECTOR_FOR(i, W<T>::c, 1)
       {
-         out.v[2*i]   = a.v[2*i]*b.v[2*i];
-         out.v[2*i+1] = a.v[2*i]*b.v[2*i+1];
+	out.v[2*i]   = a.v[2*i]*b.v[2*i];
+	out.v[2*i+1] = a.v[2*i]*b.v[2*i+1];
       }      
-      return out;
-    }
-  };
+    return out;
+  }
+};
 
-  struct MaddRealPart{
-    template <typename T>
-    inline vec<T> operator()(vec<T> a, vec<T> b, vec<T> c){
-      vec<T> out;
+struct MaddRealPart{
+  template <typename T>
+  inline vec<T> operator()(vec<T> a, vec<T> b, vec<T> c){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::c, 1)
+    VECTOR_FOR(i, W<T>::c, 1)
       {
-         out.v[2*i]   = a.v[2*i]*b.v[2*i] + c.v[2*i];
-         out.v[2*i+1] = a.v[2*i]*b.v[2*i+1] + c.v[2*i+1];
+	out.v[2*i]   = a.v[2*i]*b.v[2*i] + c.v[2*i];
+	out.v[2*i+1] = a.v[2*i]*b.v[2*i+1] + c.v[2*i+1];
       }      
-      return out;
-    }
-  };
+    return out;
+  }
+};
   
-  struct MultComplex{
-    // Complex
-    template <typename T>
-    inline vec<T> operator()(vec<T> a, vec<T> b){
-      vec<T> out;
+struct MultComplex{
+  // Complex
+  template <typename T>
+  inline vec<T> operator()(vec<T> a, vec<T> b){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::c, 1)
+    VECTOR_FOR(i, W<T>::c, 1)
       {
         cmul(a.v, b.v, out.v, 2*i);
       }      
       
-      return out;
-    }
-  };
+    return out;
+  }
+};
   
-  #undef cmul
+#undef cmul
 
-  struct Div{
-    // Real
-    template <typename T>
-    inline vec<T> operator()(vec<T> a, vec<T> b){
-      vec<T> out;
+struct Div{
+  // Real
+  template <typename T>
+  inline vec<T> operator()(vec<T> a, vec<T> b){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::r, 1)
+    VECTOR_FOR(i, W<T>::r, 1)
       {
         out.v[i] = a.v[i]/b.v[i];
       }
       
-      return out;
-    }
-  };
+    return out;
+  }
+};
   
-  #define conj(a, b, i)\
-  b[i]   = a[i];\
+#define conj(a, b, i)				\
+  b[i]   = a[i];				\
   b[i+1] = -a[i+1];
   
-  struct Conj{
-    // Complex
-    template <typename T>
-    inline vec<T> operator()(vec<T> a){
-      vec<T> out;
+struct Conj{
+  // Complex
+  template <typename T>
+  inline vec<T> operator()(vec<T> a){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::c, 1)
+    VECTOR_FOR(i, W<T>::c, 1)
       {
         conj(a.v, out.v, 2*i);
       }
       
-      return out;
-    }
-  };
+    return out;
+  }
+};
   
-  #undef conj
+#undef conj
 
-  #define timesmi(a, b, i)\
-  b[i]   = a[i+1];\
+#define timesmi(a, b, i)			\
+  b[i]   = a[i+1];				\
   b[i+1] = -a[i];
   
-  struct TimesMinusI{
-    // Complex
-    template <typename T>
-    inline vec<T> operator()(vec<T> a, vec<T> b){
-      vec<T> out;
+struct TimesMinusI{
+  // Complex
+  template <typename T>
+  inline vec<T> operator()(vec<T> a, vec<T> b){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::c, 1)
+    VECTOR_FOR(i, W<T>::c, 1)
       {
         timesmi(a.v, out.v, 2*i);
       }
       
-      return out;
-    }
-  };
+    return out;
+  }
+};
 
-  #undef timesmi
+#undef timesmi
   
-  #define timesi(a, b, i)\
-  b[i]   = -a[i+1];\
+#define timesi(a, b, i)				\
+  b[i]   = -a[i+1];				\
   b[i+1] = a[i];
   
-  struct TimesI{
-    // Complex
-    template <typename T>
-    inline vec<T> operator()(vec<T> a, vec<T> b){
-      vec<T> out;
+struct TimesI{
+  // Complex
+  template <typename T>
+  inline vec<T> operator()(vec<T> a, vec<T> b){
+    vec<T> out;
       
-      VECTOR_FOR(i, W<T>::c, 1)
+    VECTOR_FOR(i, W<T>::c, 1)
       {
         timesi(a.v, out.v, 2*i);
       }
       
-      return out;
-    }
-  };
+    return out;
+  }
+};
   
-  #undef timesi
+#undef timesi
 
-  struct PrecisionChange {
-    static inline vech StoH (const vecf &a,const vecf &b) {
-      vech ret;
+struct PrecisionChange {
+  static inline vech StoH (const vecf &a,const vecf &b) {
+    vech ret;
 #ifdef USE_FP16
-      vech *ha = (vech *)&a;
-      vech *hb = (vech *)&b;
-      const int nf = W<float>::r;
-      //      VECTOR_FOR(i, nf,1){ ret.v[i]    = ( (uint16_t *) &a.v[i])[1] ; }
-      //      VECTOR_FOR(i, nf,1){ ret.v[i+nf] = ( (uint16_t *) &b.v[i])[1] ; }
-      VECTOR_FOR(i, nf,1){ ret.v[i]    = ha->v[2*i+1]; }
-      VECTOR_FOR(i, nf,1){ ret.v[i+nf] = hb->v[2*i+1]; }
+    vech *ha = (vech *)&a;
+    vech *hb = (vech *)&b;
+    const int nf = W<float>::r;
+    //      VECTOR_FOR(i, nf,1){ ret.v[i]    = ( (uint16_t *) &a.v[i])[1] ; }
+    //      VECTOR_FOR(i, nf,1){ ret.v[i+nf] = ( (uint16_t *) &b.v[i])[1] ; }
+    VECTOR_FOR(i, nf,1){ ret.v[i]    = ha->v[2*i+1]; }
+    VECTOR_FOR(i, nf,1){ ret.v[i+nf] = hb->v[2*i+1]; }
 #else
-      assert(0);
+    assert(0);
 #endif
-      return ret;
-    }
-    static inline void  HtoS (vech h,vecf &sa,vecf &sb) {
+    return ret;
+  }
+  static inline void  HtoS (vech h,vecf &sa,vecf &sb) {
 #ifdef USE_FP16
-      const int nf = W<float>::r;
-      const int nh = W<uint16_t>::r;
-      vech *ha = (vech *)&sa;
-      vech *hb = (vech *)&sb;
-      VECTOR_FOR(i, nf, 1){ sb.v[i]= sa.v[i] = 0; }
-      //      VECTOR_FOR(i, nf, 1){ ( (uint16_t *) (&sa.v[i]))[1] = h.v[i];}
-      //      VECTOR_FOR(i, nf, 1){ ( (uint16_t *) (&sb.v[i]))[1] = h.v[i+nf];}
-      VECTOR_FOR(i, nf, 1){ ha->v[2*i+1]=h.v[i]; }
-      VECTOR_FOR(i, nf, 1){ hb->v[2*i+1]=h.v[i+nf]; }
+    const int nf = W<float>::r;
+    const int nh = W<uint16_t>::r;
+    vech *ha = (vech *)&sa;
+    vech *hb = (vech *)&sb;
+    VECTOR_FOR(i, nf, 1){ sb.v[i]= sa.v[i] = 0; }
+    //      VECTOR_FOR(i, nf, 1){ ( (uint16_t *) (&sa.v[i]))[1] = h.v[i];}
+    //      VECTOR_FOR(i, nf, 1){ ( (uint16_t *) (&sb.v[i]))[1] = h.v[i+nf];}
+    VECTOR_FOR(i, nf, 1){ ha->v[2*i+1]=h.v[i]; }
+    VECTOR_FOR(i, nf, 1){ hb->v[2*i+1]=h.v[i+nf]; }
 #else
-      assert(0);
+    assert(0);
 #endif
-    }
-    static inline vecf DtoS (vecd a,vecd b) {
-      const int nd = W<double>::r;
-      const int nf = W<float>::r;
-      vecf ret;
-      VECTOR_FOR(i, nd,1){ ret.v[i]    = a.v[i] ; }
-      VECTOR_FOR(i, nd,1){ ret.v[i+nd] = b.v[i] ; }
-      return ret;
-    }
-    static inline void StoD (vecf s,vecd &a,vecd &b) {
-      const int nd = W<double>::r;
-      VECTOR_FOR(i, nd,1){ a.v[i] = s.v[i] ; }
-      VECTOR_FOR(i, nd,1){ b.v[i] = s.v[i+nd] ; }
-    }
-    static inline vech DtoH (vecd a,vecd b,vecd c,vecd d) {
-      vecf sa,sb;
-      sa = DtoS(a,b);
-      sb = DtoS(c,d);
-      return StoH(sa,sb);
-    }
-    static inline void HtoD (vech h,vecd &a,vecd &b,vecd &c,vecd &d) {
-      vecf sa,sb;
-      HtoS(h,sa,sb);
-      StoD(sa,a,b);
-      StoD(sb,c,d);
-    }
-  };
+  }
+  static inline vecf DtoS (vecd a,vecd b) {
+    const int nd = W<double>::r;
+    const int nf = W<float>::r;
+    vecf ret;
+    VECTOR_FOR(i, nd,1){ ret.v[i]    = a.v[i] ; }
+    VECTOR_FOR(i, nd,1){ ret.v[i+nd] = b.v[i] ; }
+    return ret;
+  }
+  static inline void StoD (vecf s,vecd &a,vecd &b) {
+    const int nd = W<double>::r;
+    VECTOR_FOR(i, nd,1){ a.v[i] = s.v[i] ; }
+    VECTOR_FOR(i, nd,1){ b.v[i] = s.v[i+nd] ; }
+  }
+  static inline vech DtoH (vecd a,vecd b,vecd c,vecd d) {
+    vecf sa,sb;
+    sa = DtoS(a,b);
+    sb = DtoS(c,d);
+    return StoH(sa,sb);
+  }
+  static inline void HtoD (vech h,vecd &a,vecd &b,vecd &c,vecd &d) {
+    vecf sa,sb;
+    HtoS(h,sa,sb);
+    StoD(sa,a,b);
+    StoD(sb,c,d);
+  }
+};
 
-  //////////////////////////////////////////////
-  // Exchange support
-  struct Exchange{
+//////////////////////////////////////////////
+// Exchange support
+struct Exchange{
 
-    template <typename T,int n>
+  template <typename T,int n>
     static inline void ExchangeN(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
-      const int w = W<T>::r;
-      unsigned int mask = w >> (n + 1);
-      //      std::cout << " Exchange "<<n<<" nsimd "<<w<<" mask 0x" <<std::hex<<mask<<std::dec<<std::endl;
-      VECTOR_FOR(i, w, 1) {	
-	int j1 = i&(~mask);
-	if  ( (i&mask) == 0 ) { out1.v[i]=in1.v[j1];}
-	else                  { out1.v[i]=in2.v[j1];}
-	int j2 = i|mask;
-	if  ( (i&mask) == 0 ) { out2.v[i]=in1.v[j2];}
-	else                  { out2.v[i]=in2.v[j2];}
-      }      
-    }
-    template <typename T>
-    static inline void Exchange0(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
-      ExchangeN<T,0>(out1,out2,in1,in2);
-    };
-    template <typename T>
-    static inline void Exchange1(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
-      ExchangeN<T,1>(out1,out2,in1,in2);
-    };
-    template <typename T>
-    static inline void Exchange2(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
-      ExchangeN<T,2>(out1,out2,in1,in2);
-    };
-    template <typename T>
-    static inline void Exchange3(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
-      ExchangeN<T,3>(out1,out2,in1,in2);
-    };
+    const int w = W<T>::r;
+    unsigned int mask = w >> (n + 1);
+    //      std::cout << " Exchange "<<n<<" nsimd "<<w<<" mask 0x" <<std::hex<<mask<<std::dec<<std::endl;
+    VECTOR_FOR(i, w, 1) {	
+      int j1 = i&(~mask);
+      if  ( (i&mask) == 0 ) { out1.v[i]=in1.v[j1];}
+      else                  { out1.v[i]=in2.v[j1];}
+      int j2 = i|mask;
+      if  ( (i&mask) == 0 ) { out2.v[i]=in1.v[j2];}
+      else                  { out2.v[i]=in2.v[j2];}
+    }      
+  }
+  template <typename T>
+  static inline void Exchange0(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
+    ExchangeN<T,0>(out1,out2,in1,in2);
   };
+  template <typename T>
+  static inline void Exchange1(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
+    ExchangeN<T,1>(out1,out2,in1,in2);
+  };
+  template <typename T>
+  static inline void Exchange2(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
+    ExchangeN<T,2>(out1,out2,in1,in2);
+  };
+  template <typename T>
+  static inline void Exchange3(vec<T> &out1,vec<T> &out2,vec<T> &in1,vec<T> &in2){
+    ExchangeN<T,3>(out1,out2,in1,in2);
+  };
+};
 
 
-  //////////////////////////////////////////////
-  // Some Template specialization
-  #define perm(a, b, n, w)\
-  unsigned int _mask = w >> (n + 1);\
-  VECTOR_FOR(i, w, 1)\
-  {\
-    b[i] = a[i^_mask];\
+//////////////////////////////////////////////
+// Some Template specialization
+#define perm(a, b, n, w)			\
+  unsigned int _mask = w >> (n + 1);		\
+  VECTOR_FOR(i, w, 1)				\
+  {						\
+    b[i] = a[i^_mask];				\
   }
   
-  #define DECL_PERMUTE_N(n)\
-  template <typename T>\
-  static inline vec<T> Permute##n(vec<T> in) {\
-    vec<T> out;\
-    perm(in.v, out.v, n, W<T>::r);\
-    return out;\
+#define DECL_PERMUTE_N(n)			\
+  template <typename T>				\
+  static inline vec<T> Permute##n(vec<T> in) {	\
+    vec<T> out;					\
+    perm(in.v, out.v, n, W<T>::r);		\
+    return out;					\
   }
   
-  struct Permute{
-    DECL_PERMUTE_N(0);
-    DECL_PERMUTE_N(1);
-    DECL_PERMUTE_N(2);
-    DECL_PERMUTE_N(3);
-  };
+struct Permute{
+  DECL_PERMUTE_N(0);
+  DECL_PERMUTE_N(1);
+  DECL_PERMUTE_N(2);
+  DECL_PERMUTE_N(3);
+};
   
-  #undef perm
-  #undef DECL_PERMUTE_N
+#undef perm
+#undef DECL_PERMUTE_N
   
-  #define rot(a, b, n, w)\
-  VECTOR_FOR(i, w, 1)\
-  {\
-    b[i] = a[(i + n)%w];\
+#define rot(a, b, n, w)				\
+  VECTOR_FOR(i, w, 1)				\
+  {						\
+    b[i] = a[(i + n)%w];			\
   }
   
-  struct Rotate{
+struct Rotate{
       
-    template <int n, typename T> static inline vec<T> tRotate(vec<T> in){
-      return rotate(in, n);
-    }
+  template <int n, typename T> static inline vec<T> tRotate(vec<T> in){
+    return rotate(in, n);
+  }
     
-    template <typename T>
-    static inline vec<T> rotate(vec<T> in, int n){
-      vec<T> out;
+  template <typename T>
+  static inline vec<T> rotate(vec<T> in, int n){
+    vec<T> out;
       
-      rot(in.v, out.v, n, W<T>::r);
+    rot(in.v, out.v, n, W<T>::r);
       
-      return out;
-    }
-  };
-
-  #undef rot
-  
-  #define acc(v, a, off, step, n)\
-  for (unsigned int i = off; i < n; i += step)\
-  {\
-    a += v[i];\
+    return out;
   }
+};
+
+#undef rot
   
-  template <typename Out_type, typename In_type>
+#define acc(v, a, off, step, n)			\
+  for (unsigned int i = off; i < n; i += step)	\
+    {						\
+      a += v[i];				\
+    }
+  
+template <typename Out_type, typename In_type>
   struct Reduce{
     //Need templated class to overload output type
     //General form must generate error if compiled
@@ -442,89 +442,91 @@ namespace Optimization {
     }
   };
   
-  //Complex float Reduce
-  template <>
-  inline Grid::ComplexF Reduce<Grid::ComplexF, vecf>::operator()(vecf in){
-    float a = 0.f, b = 0.f;
+//Complex float Reduce
+template <>
+inline Grid::ComplexF Reduce<Grid::ComplexF, vecf>::operator()(vecf in){
+  float a = 0.f, b = 0.f;
     
-    acc(in.v, a, 0, 2, W<float>::r);
-    acc(in.v, b, 1, 2, W<float>::r);
+  acc(in.v, a, 0, 2, W<float>::r);
+  acc(in.v, b, 1, 2, W<float>::r);
     
-    return Grid::ComplexF(a, b);
-  }
-  
-  //Real float Reduce
-  template<>
-  inline Grid::RealF Reduce<Grid::RealF, vecf>::operator()(vecf in){
-    float a = 0.;
-    
-    acc(in.v, a, 0, 1, W<float>::r);
-    
-    return a;
-  }
-  
-  //Complex double Reduce
-  template<>
-  inline Grid::ComplexD Reduce<Grid::ComplexD, vecd>::operator()(vecd in){
-    double a = 0., b = 0.;
-    
-    acc(in.v, a, 0, 2, W<double>::r);
-    acc(in.v, b, 1, 2, W<double>::r);
-    
-    return Grid::ComplexD(a, b);
-  }
-  
-  //Real double Reduce
-  template<>
-  inline Grid::RealD Reduce<Grid::RealD, vecd>::operator()(vecd in){
-    double a = 0.f;
-    
-    acc(in.v, a, 0, 1, W<double>::r);
-    
-    return a;
-  }
-
-  //Integer Reduce
-  template<>
-  inline Integer Reduce<Integer, veci>::operator()(veci in){
-    Integer a = 0;
-    
-    acc(in.v, a, 0, 1, W<Integer>::r);
-    
-    return a;
-  }
-
-  #undef acc  // EIGEN compatibility
+  return Grid::ComplexF(a, b);
 }
+  
+//Real float Reduce
+template<>
+inline Grid::RealF Reduce<Grid::RealF, vecf>::operator()(vecf in){
+  float a = 0.;
+    
+  acc(in.v, a, 0, 1, W<float>::r);
+    
+  return a;
+}
+  
+//Complex double Reduce
+template<>
+inline Grid::ComplexD Reduce<Grid::ComplexD, vecd>::operator()(vecd in){
+  double a = 0., b = 0.;
+    
+  acc(in.v, a, 0, 2, W<double>::r);
+  acc(in.v, b, 1, 2, W<double>::r);
+    
+  return Grid::ComplexD(a, b);
+}
+  
+//Real double Reduce
+template<>
+inline Grid::RealD Reduce<Grid::RealD, vecd>::operator()(vecd in){
+  double a = 0.f;
+    
+  acc(in.v, a, 0, 1, W<double>::r);
+    
+  return a;
+}
+
+//Integer Reduce
+template<>
+inline Integer Reduce<Integer, veci>::operator()(veci in){
+  Integer a = 0;
+    
+  acc(in.v, a, 0, 1, W<Integer>::r);
+    
+  return a;
+}
+
+#undef acc  // EIGEN compatibility
+NAMESPACE_END(Optimization)
 
 //////////////////////////////////////////////////////////////////////////////////////
 // Here assign types 
 
-  typedef Optimization::vech SIMD_Htype; // Reduced precision type
-  typedef Optimization::vecf SIMD_Ftype; // Single precision type
-  typedef Optimization::vecd SIMD_Dtype; // Double precision type
-  typedef Optimization::veci SIMD_Itype; // Integer type
+typedef Optimization::vech SIMD_Htype; // Reduced precision type
+typedef Optimization::vecf SIMD_Ftype; // Single precision type
+typedef Optimization::vecd SIMD_Dtype; // Double precision type
+typedef Optimization::veci SIMD_Itype; // Integer type
 
-  // prefetch utilities
-  inline void v_prefetch0(int size, const char *ptr){};
-  inline void prefetch_HINT_T0(const char *ptr){};
+// prefetch utilities
+inline void v_prefetch0(int size, const char *ptr){};
+inline void prefetch_HINT_T0(const char *ptr){};
 
-  // Function name aliases
-  typedef Optimization::Vsplat   VsplatSIMD;
-  typedef Optimization::Vstore   VstoreSIMD;
-  typedef Optimization::Vset     VsetSIMD;
-  typedef Optimization::Vstream  VstreamSIMD;
-  template <typename S, typename T> using ReduceSIMD = Optimization::Reduce<S,T>;
+// Function name aliases
+typedef Optimization::Vsplat   VsplatSIMD;
+typedef Optimization::Vstore   VstoreSIMD;
+typedef Optimization::Vset     VsetSIMD;
+typedef Optimization::Vstream  VstreamSIMD;
+template <typename S, typename T> using ReduceSIMD = Optimization::Reduce<S,T>;
+
+// Arithmetic operations
+typedef Optimization::Sum         SumSIMD;
+typedef Optimization::Sub         SubSIMD;
+typedef Optimization::Div         DivSIMD;
+typedef Optimization::Mult        MultSIMD;
+typedef Optimization::MultComplex MultComplexSIMD;
+typedef Optimization::MultRealPart MultRealPartSIMD;
+typedef Optimization::MaddRealPart MaddRealPartSIMD;
+typedef Optimization::Conj        ConjSIMD;
+typedef Optimization::TimesMinusI TimesMinusISIMD;
+typedef Optimization::TimesI      TimesISIMD;
+
+NAMESPACE_END(Grid)
 
-  // Arithmetic operations
-  typedef Optimization::Sum         SumSIMD;
-  typedef Optimization::Sub         SubSIMD;
-  typedef Optimization::Div         DivSIMD;
-  typedef Optimization::Mult        MultSIMD;
-  typedef Optimization::MultComplex MultComplexSIMD;
-  typedef Optimization::MultRealPart MultRealPartSIMD;
-  typedef Optimization::MaddRealPart MaddRealPartSIMD;
-  typedef Optimization::Conj        ConjSIMD;
-  typedef Optimization::TimesMinusI TimesMinusISIMD;
-  typedef Optimization::TimesI      TimesISIMD;
-}