Improving the reduction to go through our on permute.

Must also do this for avx512

INSTALL

@@ -1 +1 @@
-/usr/share/automake-1.14/INSTALL
+/opt/local/share/automake-1.15/INSTALL

configure

@@ -2483,7 +2483,7 @@ test -n "$target_alias" &&
     NONENONEs,x,x, &&
   program_prefix=${target_alias}-
 
-am__api_version='1.14'
+am__api_version='1.15'
 
 # Find a good install program.  We prefer a C program (faster),
 # so one script is as good as another.  But avoid the broken or
@@ -2655,8 +2655,8 @@ test "$program_suffix" != NONE &&
 ac_script='s/[\\$]/&&/g;s/;s,x,x,$//'
 program_transform_name=`$as_echo "$program_transform_name" | sed "$ac_script"`
 
-# expand $ac_aux_dir to an absolute path
+# Expand $ac_aux_dir to an absolute path.
-am_aux_dir=`cd $ac_aux_dir && pwd`
+am_aux_dir=`cd "$ac_aux_dir" && pwd`
 
 if test x"${MISSING+set}" != xset; then
   case $am_aux_dir in
@@ -2675,7 +2675,7 @@ else
 $as_echo "$as_me: WARNING: 'missing' script is too old or missing" >&2;}
 fi
 
-if test x"${install_sh}" != xset; then
+if test x"${install_sh+set}" != xset; then
   case $am_aux_dir in
   *\ * | *\	*)
     install_sh="\${SHELL} '$am_aux_dir/install-sh'" ;;
@@ -3003,8 +3003,8 @@ MAKEINFO=${MAKEINFO-"${am_missing_run}makeinfo"}
 # <http://lists.gnu.org/archive/html/automake/2012-07/msg00014.html>
 mkdir_p='$(MKDIR_P)'
 
-# We need awk for the "check" target.  The system "awk" is bad on
+# We need awk for the "check" target (and possibly the TAP driver).  The
-# some platforms.
+# system "awk" is bad on some platforms.
 # Always define AMTAR for backward compatibility.  Yes, it's still used
 # in the wild :-(  We should find a proper way to deprecate it ...
 AMTAR='$${TAR-tar}'
@@ -3063,6 +3063,7 @@ END
 fi
 
 
+
 ac_config_headers="$ac_config_headers lib/Grid_config.h"
 
 # Check whether --enable-silent-rules was given.

lib/simd/Grid_avx.h

@@ -14,7 +14,23 @@
 #endif
 
 namespace Optimization {
-  
+
+  template<class vtype>
+  union uconv {
+    __m256 f;
+    vtype v;
+  };
+
+  union u256f {
+    __m256 v;
+    float f[8];
+  };
+
+  union u256d {
+    __m256d v;
+    double f[4];
+  };
+
   struct Vsplat{
     //Complex float
     inline __m256 operator()(float a, float b){
@@ -54,7 +70,6 @@ namespace Optimization {
 
   };
 
-
   struct Vstream{
     //Float
     inline void operator()(float * a, __m256 b){
@@ -68,8 +83,6 @@ namespace Optimization {
 
   };
 
-
-
   struct Vset{
     // Complex float 
     inline __m256 operator()(Grid::ComplexF *a){
@@ -92,7 +105,6 @@ namespace Optimization {
       return _mm256_set_epi32(a[7],a[6],a[5],a[4],a[3],a[2],a[1],a[0]);
     }
 
-
   };
 
   template <typename Out_type, typename In_type>
@@ -106,9 +118,6 @@ namespace Optimization {
     }
   };
 
-
- 
-
   /////////////////////////////////////////////////////
   // Arithmetic operations
   /////////////////////////////////////////////////////
@@ -170,7 +179,6 @@ namespace Optimization {
     }
   };
 
-
   struct MultComplex{
     // Complex float
     inline __m256 operator()(__m256 a, __m256 b){
@@ -207,7 +215,6 @@ namespace Optimization {
 	IF IMM0[3] = 0
 	THEN DEST[255:192]=SRC2[191:128] ELSE DEST[255:192]=SRC2[255:192] FI; // Ox5 r<->i   ; 0xC unchanged
       */
-      
       __m256d ymm0,ymm1,ymm2;
       ymm0 = _mm256_shuffle_pd(a,a,0x0); // ymm0 <- ar ar, ar,ar b'00,00
       ymm0 = _mm256_mul_pd(ymm0,b);      // ymm0 <- ar bi, ar br
@@ -247,7 +254,6 @@ namespace Optimization {
     }
   };
 
-
   struct Conj{
     // Complex single
     inline __m256 operator()(__m256 in){
@@ -292,18 +298,13 @@ namespace Optimization {
     }
   };
 
-
-  
-
-
   //////////////////////////////////////////////
   // Some Template specialization
+  //////////////////////////////////////////////
+
   template < typename vtype > 
-    void permute(vtype &a, vtype &b, int perm) {
+    void permute(vtype &a,vtype b, int perm) {
-    union { 
+    uconv<vtype> conv;
-      __m256 f;
-      vtype v;
-    } conv;
     conv.v = b;
     switch (perm){
       // 8x32 bits=>3 permutes
@@ -313,24 +314,20 @@ namespace Optimization {
     default: assert(0); break;
     }
     a = conv.v;
-    
   }
-  
+
   //Complex float Reduce
   template<>
     inline Grid::ComplexF Reduce<Grid::ComplexF, __m256>::operator()(__m256 in){
     __m256 v1,v2;
-    union { 
-      __m256 v;
-      float f[8];
-    } conv;
     Optimization::permute(v1,in,0); // sse 128; paired complex single
     v1 = _mm256_add_ps(v1,in);
     Optimization::permute(v2,v1,1); // avx 256; quad complex single
     v1 = _mm256_add_ps(v1,v2);
-    conv.v = v1;
+    u256f conv; conv.v = v1;
     return Grid::ComplexF(conv.f[0],conv.f[1]);
   }
+
   //Real float Reduce
   template<>
   inline Grid::RealF Reduce<Grid::RealF, __m256>::operator()(__m256 in){
@@ -341,7 +338,8 @@ namespace Optimization {
     v1 = _mm256_add_ps(v1,v2);
     Optimization::permute(v2,v1,2); 
     v1 = _mm256_add_ps(v1,v2);
-    return v1[0];
+    u256f conv; conv.v=v1;
+    return conv.f[0];
   }
   
   
@@ -351,7 +349,8 @@ namespace Optimization {
     __m256d v1;
     Optimization::permute(v1,in,0); // sse 128; paired complex single
     v1 = _mm256_add_pd(v1,in);
-    return Grid::ComplexD(v1[0],v1[1]);
+    u256d conv; conv.v = v1;
+    return Grid::ComplexD(conv.f[0],conv.f[1]);
   }
   
   //Real double Reduce
@@ -362,7 +361,8 @@ namespace Optimization {
     v1 = _mm256_add_pd(v1,in);
     Optimization::permute(v2,v1,1); 
     v1 = _mm256_add_pd(v1,v2);
-    return v1[0];
+    u256d conv; conv.v = v1;
+    return conv.f[0];
   }
 
   //Integer Reduce
@@ -390,22 +390,9 @@ namespace Grid {
       _mm_prefetch(ptr+i+512,_MM_HINT_T0);
     }
   }
-  
   template < typename VectorSIMD > 
     inline void Gpermute(VectorSIMD &y,const VectorSIMD &b, int perm ) {
-    union { 
+    Optimization::permute(y.v,b.v,perm);
-      __m256 f;
-      decltype(VectorSIMD::v) v;
-    } conv;
-    conv.v = b.v;
-    switch(perm){
-    case 3: break; //empty for AVX1/2
-    case 2: conv.f = _mm256_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(2,3,0,1)); break;
-    case 1: conv.f = _mm256_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(1,0,3,2));  break; 
-    case 0: conv.f = _mm256_permute2f128_ps(conv.f,conv.f,0x01); break;
-    default: assert(0); break;
-    }
-    y.v=conv.v;
   };
 
   // Function name aliases

lib/simd/Grid_sse4.h

@@ -10,6 +10,21 @@
 #include <pmmintrin.h>
 
 namespace Optimization {
+
+  template<class vtype>
+  union uconv {
+    __m128 f;
+    vtype v;
+  };
+
+  union u128f {
+    __m128 v;
+    float f[4];
+  };
+  union u128d {
+    __m128d v;
+    double f[2];
+  };
   
   struct Vsplat{
     //Complex float
@@ -50,7 +65,6 @@ namespace Optimization {
 
   };
 
-
   struct Vstream{
     //Float
     inline void operator()(float * a, __m128 b){
@@ -64,8 +78,6 @@ namespace Optimization {
 
   };
 
-
-
   struct Vset{
     // Complex float 
     inline __m128 operator()(Grid::ComplexF *a){
@@ -102,9 +114,6 @@ namespace Optimization {
     }
   };
 
-
- 
-
   /////////////////////////////////////////////////////
   // Arithmetic operations
   /////////////////////////////////////////////////////
@@ -138,7 +147,6 @@ namespace Optimization {
     }
   };
 
-
   struct MultComplex{
     // Complex float
     inline __m128 operator()(__m128 a, __m128 b){
@@ -177,7 +185,6 @@ namespace Optimization {
     }
   };
 
-
   struct Conj{
     // Complex single
     inline __m128 operator()(__m128 in){
@@ -216,57 +223,61 @@ namespace Optimization {
       __m128d tmp = _mm_shuffle_pd(in,in,0x1);
       return _mm_addsub_pd(_mm_setzero_pd(),tmp); // r,-i
     }
-
-
   };
 
-
-
-
-
   //////////////////////////////////////////////
   // Some Template specialization
-  
+  template < typename vtype > 
+    void permute(vtype &a, vtype b, int perm) {
+    uconv<vtype> conv; 
+    conv.v = b;
+    switch(perm){
+    case 3: break; //empty for SSE4
+    case 2: break; //empty for SSE4
+    case 1: conv.f = _mm_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(2,3,0,1)); break;
+    case 0: conv.f = _mm_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(1,0,3,2)); break;
+    default: assert(0); break;
+    }
+    a=conv.v;
+  }; 
+
   //Complex float Reduce
   template<>
   inline Grid::ComplexF Reduce<Grid::ComplexF, __m128>::operator()(__m128 in){
-    union {
+    __m128 v1; // two complex
-      __m128 v1;  
+    Optimization::permute(v1,in,0); 
-      float f[4]; 
+    v1= _mm_add_ps(v1,in);
-    } u128;
+    u128f conv;    conv.v=v1;
-    u128.v1 = _mm_add_ps(in, _mm_shuffle_ps(in,in, 0b01001110)); // FIXME Prefer to use _MM_SHUFFLE macros
+    return Grid::ComplexF(conv.f[0],conv.f[1]);
-    return Grid::ComplexF(u128.f[0], u128.f[1]);   
   }
   //Real float Reduce
   template<>
   inline Grid::RealF Reduce<Grid::RealF, __m128>::operator()(__m128 in){
-    // FIXME Hack
+    __m128 v1,v2; // quad single
-    const Grid::RealF * ptr = (const Grid::RealF *) &in;
+    Optimization::permute(v1,in,0); 
-    Grid::RealF ret = 0; 
+    v1= _mm_add_ps(v1,in);
-    for(int i=0;i< 4 ;i++){ // 4 number of simd lanes for float
+    Optimization::permute(v2,v1,1); 
-      ret = ret+ptr[i];
+    v1 = _mm_add_ps(v1,v2);
-    }
+    u128f conv; conv.v=v1;
-    return ret;
+    return conv.f[0];
   }
   
   
   //Complex double Reduce
   template<>
   inline Grid::ComplexD Reduce<Grid::ComplexD, __m128d>::operator()(__m128d in){
-    printf("Reduce : Missing good complex double implementation -> FIX\n");
+    u128d conv; conv.v = in;
-    return Grid::ComplexD(in[0], in[1]); // inefficient
+    return Grid::ComplexD(conv.f[0],conv.f[1]);
   }
   
   //Real double Reduce
   template<>
   inline Grid::RealD Reduce<Grid::RealD, __m128d>::operator()(__m128d in){
-    // FIXME Hack
+    __m128d v1;
-    const Grid::RealD * ptr =(const Grid::RealD *)  &in;
+    Optimization::permute(v1,in,0); // avx 256; quad double
-    Grid::RealD ret = 0; 
+    v1 = _mm_add_pd(v1,in);
-    for(int i=0;i< 2 ;i++){// 2 number of simd lanes for float
+    u128d conv; conv.v = v1;
-      ret = ret+ptr[i];
+    return conv.f[0];
-    }
-    return ret;
   }
 
   //Integer Reduce
@@ -276,12 +287,6 @@ namespace Optimization {
    printf("Reduce : Missing integer implementation -> FIX\n");
     assert(0);
   }
-  
-
-
-
-
-  
 }
 
 //////////////////////////////////////////////////////////////////////////////////////
@@ -292,27 +297,13 @@ namespace Grid {
   typedef __m128d SIMD_Dtype; // Double precision type
   typedef __m128i SIMD_Itype; // Integer type
 
-
   inline void v_prefetch0(int size, const char *ptr){};  // prefetch utilities
 
   // Gpermute function
   template < typename VectorSIMD > 
     inline void Gpermute(VectorSIMD &y,const VectorSIMD &b, int perm ) {
-    union { 
+    Optimization::permute(y.v,b.v,perm);
-      __m128 f;
+  }
-      decltype(VectorSIMD::v) v;
-    } conv;
-    conv.v = b.v;
-    switch(perm){
-    case 3: break; //empty for SSE4
-    case 2: break; //empty for SSE4
-    case 1: conv.f = _mm_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(2,3,0,1)); break;
-    case 0: conv.f = _mm_shuffle_ps(conv.f,conv.f,_MM_SHUFFLE(1,0,3,2)); break;
-    default: assert(0); break;
-    }
-    y.v=conv.v;
-  }; 
-  
 
 
   // Function name aliases