/Users/mourrain/Devel/mmx/numerix/include/numerix/sse.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : sse.hpp
* DESCRIPTION: Wrapper for SSE instructions
* COPYRIGHT  : (C) 2008  Joris van der Hoeven and Gregoire Lecerf
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_SSE_HPP
#define __MMX_SSE_HPP
#include <numerix/simd.hpp>

#if defined (NUMERIX_ENABLE_SIMD) && defined (__SSE2__)
#include <stdint.h>
#ifdef __SSE2__
#include <emmintrin.h>
#endif
#ifdef __SSE3__
#include <pmmintrin.h>
#endif
#ifdef __SSSE3__
#include <tmmintrin.h>
#endif
#ifdef __SSE4A__
#include <ammintrin.h>
#endif

#include <basix/compound.hpp>
#include <basix/identifiers.hpp>
#include <basix/syntactic.hpp>
#include <numerix/complex.hpp>

namespace mmx {

/******************************************************************************
* Generic vectorial routines
******************************************************************************/
  
template<typename C> inline typename Simd_type (C)
simd_load_aligned (const C* v) {
  return (typename Simd_type (C))
    _mm_load_si128 ((const __m128i*) v); }

template<typename C> inline void
simd_save_aligned (C* v, const typename Simd_type (C)& x) {
  _mm_store_si128 ((__m128i*) v, (const __m128i) x); }

template<typename C> inline void
simd_save (C* v, const typename Simd_type(C)& x) {
  _mm_storeu_si128 ((__m128i*) v, (const __m128i) x); }

template<typename C> inline typename Simd_type (C)
simd_load (const C* v0, const C* v1) {
  return simd_set (*v0, *v1); }

template<typename C> inline void
simd_save (C* v0, C* v1, const typename Simd_type (C)& x) {
  static C v[Simd_size (C)]; // ensures alignment
  simd_save_aligned (v, x);
  *v0 = v[0]; *v1 = v[1]; }

template<typename C> inline typename Simd_type (C)
simd_load (const C* v0, const C* v1, const C* v2, const C* v3) {
  return simd_set (*v0, *v1, *v2, *v3); }

template<typename C> inline void
simd_save (C* v0, C* v1, C* v2, C* v3, const typename Simd_type (C)& x) {
  static C v[Simd_size (C)]; // ensures alignment
  simd_save_aligned (v, x);
  *v0 = v[0]; *v1 = v[1]; *v2 = v[2]; *v3 = v[3]; }

template<typename C> inline typename Simd_type (C)
simd_load (const C* v0, const C* v1, const C* v2, const C* v3,
           const C* v4, const C* v5, const C* v6, const C* v7) {
  return simd_set (*v0, *v1, *v2, *v3, *v4, *v5, *v6, *v7); }

template<typename C> inline void
simd_save (C* v0, C* v1, C* v2, C* v3, C* v4, C* v5, C* v6, C* v7,
           const typename Simd_type (C)& x) {
  static C v[Simd_size (C)]; // ensures alignment
  simd_save_aligned (v, x);
  *v0 = v[0]; *v1 = v[1]; *v2 = v[2]; *v3 = v[3];
  *v4 = v[4]; *v5 = v[5]; *v6 = v[6]; *v7 = v[7]; }

template<typename C> inline typename Simd_type (C)
simd_load (const C* v0, const C* v1, const C* v2, const C* v3,
           const C* v4, const C* v5, const C* v6, const C* v7,
           const C* v8, const C* v9, const C* v10, const C* v11,
           const C* v12, const C* v13, const C* v14, const C* v15) {
  return simd_set (*v0, *v1, *v2, *v3, *v4, *v5, *v6, *v7,
                   *v8, *v9, *v10, *v11, *v12, *v13, *v14, *v15); }

template<typename C> inline void
simd_save (C* v0, C* v1, C* v2, C* v3, C* v4, C* v5, C* v6, C* v7,
           C* v8, C* v9, C* v10, C* v11, C* v12, C* v13, C* v14, C* v15,
           const typename Simd_type (C)& x) {
  static C v[Simd_size (C)]; // ensures alignment
  simd_save_aligned (v, x);
  *v0 = v[0]; *v1 = v[1]; *v2 = v[2]; *v3 = v[3];
  *v4 = v[4]; *v5 = v[5]; *v6 = v[6]; *v7 = v[7];
  *v8 = v[8]; *v9 = v[9]; *v10 = v[10]; *v11 = v[11];
  *v12 = v[12]; *v13 = v[13]; *v14 = v[14]; *v15 = v[15]; }

template<typename C> inline C
simd_big_add (const typename Simd_type (C)& x) {
  C r = 0;
  for (nat i = 0; i < Simd_size (C); i++)
    r += ((C*) &x) [i];
  return r; }
 
/******************************************************************************
* Vector light interface
******************************************************************************/

// Note that vectorial +, -, *, /, unary minus, ^, |, &, ~ operations
// are supposed to be furnished by the compiler

template<typename V> inline syntactic
simd_flatten (const V& x) {
  typedef typename Simd_base_type(V) C;
  static const nat size = Simd_size(C);
  C* v = mmx_new<C> (size);
  simd_save_aligned (v, x);
  vector<syntactic> w = fill <syntactic> (size);;
  for (nat i = 0; i < size; i++)
    w[i] = flatten (v[i]);
  mmx_delete<C> (v, size);
  return apply (GEN_SQTUPLE, w); 
}

#define SIMD_SUGAR(C,V)                                                 \
  inline syntactic flatten (const V& x) {                               \
    return simd_flatten (x); }                                          \
  inline bool equal (const V& x, const V& y) {                          \
    return _mm_movemask_epi8 (                                          \
      (__m128i) _mm_cmpeq_epi32 ((const __m128i) x,                     \
                                 (const __m128i) y)) == 131071; }       \
  inline bool unequal (const V& x, const V& y) {                        \
    return ! equal (x, y); }                                            \
  STMPL inline void clear (V& x) {                                      \
    x = simd_set_duplicate (C (0)); }                                   \
  STMPL inline void mul (V& x, const V& y1, const C& y2) {              \
    x = y1 * simd_set_duplicate (y2); }                                 \
  STMPL inline void mul_add (V& x, const V& y1, const C& y2) {          \
    x += y1 * simd_set_duplicate (y2); }

/******************************************************************************
* Vectors of two doubles
******************************************************************************/

typedef double __attribute__((vector_size(16))) sse_double;

template<>
struct simd_helper<double> {
  typedef sse_double type;
  static const nat size = 2; };

template<>
struct simd_base_helper<sse_double> {
  typedef double type; };

inline sse_double
simd_load (const double* v) {
  return _mm_loadu_pd (v); }

inline sse_double
simd_set_duplicate (double x) {
  return _mm_set1_pd (x); }

inline sse_double
simd_set (double v0, double v1) {
  return _mm_set_pd (v1, v0); }

#ifdef __SSE3__
STMPL inline double
simd_big_add (const sse_double& x) {
  double r;
  sse_double y = _mm_hadd_pd (x, simd_set_duplicate((double) 0));
  _mm_storel_pd (&r, y);
  return r;
}
#endif

// Comparisons
inline sse_double
simd_equal (const sse_double& x, const sse_double& y) {
  return _mm_cmpeq_pd (x, y); }
  
inline sse_double
simd_unequal (const sse_double& x, const sse_double& y) {
  return _mm_cmpneq_pd (x, y); }

inline sse_double
simd_less (const sse_double& x, const sse_double& y) {
  return _mm_cmplt_pd (x, y); }

inline sse_double
simd_gtr (const sse_double& x, const sse_double& y) {
  return _mm_cmpgt_pd (x, y); }

inline sse_double
simd_lesseq (const sse_double& x, const sse_double& y) {
  return _mm_cmple_pd (x, y); }

inline sse_double
simd_gtreq (const sse_double& x, const sse_double& y) {
  return _mm_cmpge_pd (x, y); }

// Min, max
inline sse_double
min (const sse_double& x, const sse_double& y) { 
  return _mm_min_pd (x, y); }

inline sse_double
max (const sse_double& x, const sse_double& y) {
  return _mm_max_pd (x, y); }

inline sse_double
simd_shuffle (const sse_double& x, const sse_double& y, int i) {
  return _mm_shuffle_pd (x, y, i); }

// Specific
inline sse_double
simd_load_duplicate (const double* v) {
  return _mm_load1_pd (v); }

inline sse_double
simd_load (const double* v0, const double* v1) {
  return _mm_loadh_pd (_mm_load1_pd (v0), v1); }

inline void
simd_save (double* v0, double* v1, const sse_double& x) {
  _mm_storel_pd (v0, x); _mm_storeh_pd (v1, x); }

inline sse_double
simd_swap (const sse_double& x) {
  return _mm_shuffle_pd (x, x, 1); }

// Printing and equalities
SIMD_SUGAR (double, sse_double)

/******************************************************************************
* Vectors of two complexified doubles
******************************************************************************/

typedef complex<    double>      complex_double;
typedef complex<sse_double>  sse_complex_double;

template<>
struct simd_helper<complex_double> {
  typedef sse_complex_double type;
  static const nat size = 2; };

inline sse_complex_double
simd_set_duplicate (const complex_double& z) {
  return sse_complex_double (simd_set_duplicate (Re (z)),
                             simd_set_duplicate (Im (z))); }
inline sse_complex_double
simd_set (const complex_double& z0, const complex_double& z1) {
  return sse_complex_double (simd_set (Re (z0), Re (z1)),
                             simd_set (Im (z0), Im (z1))); }
inline sse_complex_double
simd_load_duplicate (const complex_double* v) {
  const double* w= (double*) ((void*) v);
  return sse_complex_double (simd_load_duplicate (w),
                             simd_load_duplicate (w + 1)); }

template<> inline syntactic
flatten (const sse_complex_double& z) {
  return flatten (Re (z)) + flatten (Im (z)) * Imaginary (syntactic); }

/******************************************************************************
* Vectors of int64_t
******************************************************************************/

typedef int64_t __attribute__((vector_size(16))) sse_int64_t;

template<>
struct simd_helper<int64_t> {
  typedef sse_int64_t type;
  static const nat size = 2; };

template<>
struct simd_base_helper<sse_int64_t> {
  typedef int64_t type; };

inline sse_int64_t
simd_set_duplicate (int64_t x) {
  return sse_int64_t(_mm_set1_epi64 ((__m64) x)); }

inline sse_int64_t
simd_set (int64_t x0, int64_t x1) {
  return sse_int64_t(_mm_set_epi64 ((__m64) x1, (__m64) x0)); }

// Shifts
inline sse_int64_t
simd_sll (const sse_int64_t& x, int i) {
  return sse_int64_t(_mm_slli_epi64 ((__m128i) x, i)); }

inline sse_int64_t
simd_srl (const sse_int64_t& x, int i) {
  return sse_int64_t(_mm_srli_epi64 ((__m128i) x, i)); }

// Printing and equalities
SIMD_SUGAR (int64_t, sse_int64_t)

/******************************************************************************
* Vectors of uint64_t
******************************************************************************/

typedef uint64_t __attribute__((vector_size(16))) sse_uint64_t;

template<>
struct simd_helper<uint64_t> {
  typedef sse_uint64_t type;
  static const nat size = 2; };

template<>
struct simd_base_helper<sse_uint64_t> {
  typedef uint64_t type; };

inline sse_uint64_t
simd_set_duplicate (uint64_t x) {
  return sse_uint64_t(_mm_set1_epi64 ((__m64) x)); }

inline sse_uint64_t
simd_set (uint64_t x0, uint64_t x1) {
  return sse_uint64_t(_mm_set_epi64 ((__m64) x1, (__m64) x0)); }

// Shifts
inline sse_uint64_t
simd_sll (const sse_uint64_t& x, int i) {
  return sse_uint64_t(_mm_slli_epi64 ((__m128i) x, i)); }

inline sse_uint64_t
simd_srl (const sse_uint64_t& x, int i) {
  return sse_uint64_t(_mm_srli_epi64 ((__m128i) x, i)); }

// Printing and equalities
SIMD_SUGAR (uint64_t, sse_uint64_t)

/******************************************************************************
* Vectors of int32_t
******************************************************************************/

typedef int32_t __attribute__((vector_size(16))) sse_int32_t;

template<>
struct simd_helper<int32_t> {
  typedef sse_int32_t type;
  static const nat size = 4; };

template<>
struct simd_base_helper<sse_int32_t> {
  typedef int32_t type; };

inline sse_int32_t
simd_set_duplicate (int32_t x) {
  return sse_int32_t(_mm_set1_epi32 (x)); }

inline sse_int32_t
simd_set (int32_t x0, int32_t x1, int32_t x2, int32_t x3) {
  return sse_int32_t(_mm_set_epi32 (x3, x2, x1, x0)); }

// Comparisons
inline sse_int32_t
simd_equal (const sse_int32_t& x, const sse_int32_t& y) {
  return sse_int32_t(_mm_cmpeq_epi32 ((__m128i) x, (__m128i) y)); }
  
inline sse_int32_t
simd_less (const sse_int32_t& x, const sse_int32_t& y) {
  return sse_int32_t(_mm_cmplt_epi32 ((__m128i) x, (__m128i) y)); }

inline sse_int32_t
simd_gtr (const sse_int32_t& x, const sse_int32_t& y) {
  return sse_int32_t(_mm_cmpgt_epi32 ((__m128i) x, (__m128i) y)); }

// Shifts
inline sse_int32_t
simd_sll (const sse_int32_t& x, int i) {
  return sse_int32_t(_mm_slli_epi32 ((__m128i) x, i)); }

inline sse_int32_t
simd_srl (const sse_int32_t& x, int i) {
  return sse_int32_t(_mm_srli_epi32 ((__m128i) x, i)); }

inline sse_int32_t
simd_sra (const sse_int32_t& x, int i) {
  return sse_int32_t(_mm_srai_epi32 ((__m128i) x, i)); }

// Printing and equalities
SIMD_SUGAR (int32_t, sse_int32_t)

/******************************************************************************
* Vectors of uint32_t
******************************************************************************/

typedef uint32_t __attribute__((vector_size(16))) sse_uint32_t;

template<>
struct simd_helper<uint32_t> {
  typedef sse_uint32_t type;
  static const nat size = 4; };

template<>
struct simd_base_helper<sse_uint32_t> {
  typedef uint32_t type; };

inline sse_uint32_t
simd_set_duplicate (uint32_t x) {
  return sse_uint32_t(_mm_set1_epi32 (x)); }

inline sse_uint32_t
simd_set (uint32_t x0, uint32_t x1, uint32_t x2, uint32_t x3) {
  return sse_uint32_t(_mm_set_epi32 (x3, x2, x1, x0)); }

// Comparisons
inline sse_uint32_t
simd_equal (const sse_uint32_t& x, const sse_uint32_t& y) {
  return sse_uint32_t(_mm_cmpeq_epi32 ((__m128i) x, (__m128i) y)); }
  
// Shifts
inline sse_uint32_t
simd_sll (const sse_uint32_t& x, int i) {
  return sse_uint32_t(_mm_slli_epi32 ((__m128i) x, i)); }

inline sse_uint32_t
simd_srl (const sse_uint32_t& x, int i) {
  return sse_uint32_t(_mm_srli_epi32 ((__m128i) x, i)); }

// Printing and equalities
SIMD_SUGAR (uint32_t, sse_uint32_t)

/******************************************************************************
* Vectors of eight int16_t
******************************************************************************/

typedef int16_t __attribute__((vector_size(16))) sse_int16_t;

template<>
struct simd_helper<int16_t> {
  typedef sse_int16_t type;
  static const nat size = 8; };

template<>
struct simd_base_helper<sse_int16_t> {
  typedef int16_t type; };

inline sse_int16_t
simd_set_duplicate (int16_t x) {
  return sse_int16_t(_mm_set1_epi16 (x)); }

inline sse_int16_t
simd_set (int16_t x0, int16_t x1, int16_t x2, int16_t x3,
          int16_t x4, int16_t x5, int16_t x6, int16_t x7) {
  return sse_int16_t(_mm_set_epi16 (x7, x6, x5, x4, x3, x2, x1, x0)); }

// Comparisons
inline sse_int16_t
simd_equal (const sse_int16_t& x, const sse_int16_t& y) {
  return sse_int16_t(_mm_cmpeq_epi16 ((__m128i) x, (__m128i) y)); }
  
inline sse_int16_t
simd_less (const sse_int16_t& x, const sse_int16_t& y) {
  return sse_int16_t(_mm_cmplt_epi16 ((__m128i) x, (__m128i) y)); }

inline sse_int16_t
simd_gtr (const sse_int16_t& x, const sse_int16_t& y) {
  return sse_int16_t(_mm_cmpgt_epi16 ((__m128i) x, (__m128i) y)); }

// Min, max
inline sse_int16_t
min (const sse_int16_t& x, const sse_int16_t& y) { 
  return sse_int16_t(_mm_min_epi16 ((__m128i) x, (__m128i) y)); }

inline sse_int16_t
max (const sse_int16_t& x, const sse_int16_t& y) {
  return sse_int16_t(_mm_max_epi16 ((__m128i) x, (__m128i) y)); }

// Shifts
inline sse_int16_t
simd_sll (const sse_int16_t& x, int i) {
  return sse_int16_t(_mm_slli_epi16 ((__m128i) x, i)); }

inline sse_int16_t
simd_sra (const sse_int16_t& x, int i) {
  return sse_int16_t(_mm_srai_epi16 ((__m128i) x, i)); }

inline sse_int16_t
simd_srl (const sse_int16_t& x, int i) {
  return sse_int16_t(_mm_srli_epi16 ((__m128i) x, i)); }

// Printing and equalities
SIMD_SUGAR (int16_t, sse_int16_t)

/******************************************************************************
* Vectors of eight uint16_t
******************************************************************************/

typedef uint16_t __attribute__((vector_size(16))) sse_uint16_t;

template<>
struct simd_helper<uint16_t> {
  typedef sse_uint16_t type;
  static const nat size = 8; };

template<>
struct simd_base_helper<sse_uint16_t> {
  typedef uint16_t type; };

inline sse_uint16_t
simd_set_duplicate (uint16_t x) {
  return sse_uint16_t(_mm_set1_epi16 ((short) x)); }

inline sse_uint16_t
simd_set (uint16_t x0, uint16_t x1, uint16_t x2, uint16_t x3,
          uint16_t x4, uint16_t x5, uint16_t x6, uint16_t x7) {
  return sse_uint16_t(_mm_set_epi16 ((short) x7, (short) x6, (short) x5,
        (short) x4, (short) x3, (short) x2, (short)  x1, (short) x0)); }

// Comparisons
inline sse_uint16_t
simd_equal (const sse_uint16_t& x, const sse_uint16_t& y) {
  return sse_uint16_t(_mm_cmpeq_epi16 ((__m128i) x, (__m128i) y)); }
  
// Shifts
inline sse_uint16_t
simd_sll (const sse_uint16_t& x, int i) {
  return sse_uint16_t(_mm_slli_epi16 ((__m128i) x, i)); }

inline sse_uint16_t
simd_srl (const sse_uint16_t& x, int i) {
  return sse_uint16_t(_mm_srli_epi16 ((__m128i) x, i)); }

// Printing and equalities
SIMD_SUGAR (uint16_t, sse_uint16_t)

/******************************************************************************
* Vectors of sixteen int8_t
******************************************************************************/

typedef int8_t __attribute__((vector_size(16))) sse_int8_t;

template<>
struct simd_helper<int8_t> {
  typedef sse_int8_t type;
  static const nat size = 16;
};

template<>
struct simd_base_helper<sse_int8_t> {
  typedef int8_t type; };

inline sse_int8_t
simd_set_duplicate (int8_t x) {
  return sse_int8_t(_mm_set1_epi8 ((char) x)); }

inline sse_int8_t 
simd_set (int8_t x0, int8_t x1,  int8_t x2,  int8_t x3,
          int8_t x4, int8_t x5,  int8_t x6,  int8_t x7,
          int8_t x8, int8_t x9,  int8_t x10, int8_t x11,
          int8_t x12, int8_t x13, int8_t x14, int8_t x15) {
  return sse_int8_t(_mm_set_epi8 (x15, x14,  x13,  x12,  x11,  x10,  x9,  x8,
                       x7, x6, x5, x4, x3, x2, x1, x0)); }

// Comparisons
inline sse_int8_t
simd_equal (const sse_int8_t& x, const sse_int8_t& y) {
  return sse_int8_t(_mm_cmpeq_epi8 ((__m128i) x, (__m128i) y)); }
  
inline sse_int8_t
simd_less (const sse_int8_t& x, const sse_int8_t& y) {
  return sse_int8_t(_mm_cmplt_epi8 ((__m128i) x, (__m128i) y)); }

inline sse_int8_t
simd_gtr (const sse_int8_t& x, const sse_int8_t& y) {
  return sse_int8_t(_mm_cmpgt_epi8 ((__m128i) x, (__m128i) y)); }

// Printing and equalities
SIMD_SUGAR (int8_t, sse_int8_t)

/******************************************************************************
* Vectors of sixteen uint8_t
******************************************************************************/

typedef uint8_t __attribute__((vector_size(16))) sse_uint8_t;

template<>
struct simd_helper<uint8_t> {
  typedef sse_uint8_t type;
  static const nat size = 16;
};

template<>
struct simd_base_helper<sse_uint8_t> {
  typedef uint8_t type; };

inline sse_uint8_t
simd_set_duplicate (uint8_t x) {
  return sse_uint8_t(_mm_set1_epi8 (x)); }

inline sse_uint8_t 
simd_set (uint8_t x0, uint8_t x1,  uint8_t x2,  uint8_t x3,
          uint8_t x4, uint8_t x5,  uint8_t x6,  uint8_t x7,
          uint8_t x8, uint8_t x9,  uint8_t x10, uint8_t x11,
          uint8_t x12, uint8_t x13, uint8_t x14, uint8_t x15) {
  return sse_uint8_t(_mm_set_epi8 (x15, x14,  x13,  x12,  x11,  x10,  x9,  x8,
                       x7, x6, x5, x4, x3, x2, x1, x0)); }

// Comparisons
inline sse_uint8_t
simd_equal (const sse_uint8_t& x, const sse_uint8_t& y) {
  return sse_uint8_t(_mm_cmpeq_epi8 ((__m128i) x, (__m128i) y)); }
  
// Min, max
inline sse_uint8_t
min (const sse_uint8_t& x, const sse_uint8_t& y) { 
  return sse_uint8_t(_mm_min_epu8 ((__m128i) x, (__m128i) y)); }

inline sse_uint8_t
max (const sse_uint8_t& x, const sse_uint8_t& y) {
  return sse_uint8_t(_mm_max_epu8 ((__m128i) x, (__m128i) y)); }

// Printing and equalities
SIMD_SUGAR (uint8_t, sse_uint8_t)

#undef SIMD_SUGAR
} // namespace mmx

#endif // NUMERIX_ENABLE_SIMD
#endif // __MMX_SSE_HPP