Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
SSE/PacketMath.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_PACKET_MATH_SSE_H
11 #define EIGEN_PACKET_MATH_SSE_H
12 
13 #include "../../InternalHeaderCheck.h"
14 
15 namespace Eigen {
16 
17 namespace internal {
18 
19 #ifndef EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
20 #define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 8
21 #endif
22 
23 #if !defined(EIGEN_VECTORIZE_AVX) && !defined(EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS)
24 // 32 bits => 8 registers
25 // 64 bits => 16 registers
26 #define EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS (2*sizeof(void*))
27 #endif
28 
29 #ifdef EIGEN_VECTORIZE_FMA
30 #ifndef EIGEN_HAS_SINGLE_INSTRUCTION_MADD
31 #define EIGEN_HAS_SINGLE_INSTRUCTION_MADD
32 #endif
33 #endif
34 
35 #if ((defined EIGEN_VECTORIZE_AVX) && (EIGEN_COMP_GNUC_STRICT || EIGEN_COMP_MINGW || EIGEN_COMP_LCC) && (__GXX_ABI_VERSION < 1004)) || EIGEN_OS_QNX
36 // With GCC's default ABI version, a __m128 or __m256 are the same types and therefore we cannot
37 // have overloads for both types without linking error.
38 // One solution is to increase ABI version using -fabi-version=4 (or greater).
39 // Otherwise, we workaround this inconvenience by wrapping 128bit types into the following helper
40 // structure:
41 typedef eigen_packet_wrapper<__m128> Packet4f;
42 typedef eigen_packet_wrapper<__m128d> Packet2d;
43 #else
44 typedef __m128 Packet4f;
45 typedef __m128d Packet2d;
46 #endif
47 
48 typedef eigen_packet_wrapper<__m128i, 0> Packet4i;
49 typedef eigen_packet_wrapper<__m128i, 1> Packet16b;
50 
51 template<> struct is_arithmetic<__m128> { enum { value = true }; };
52 template<> struct is_arithmetic<__m128i> { enum { value = true }; };
53 template<> struct is_arithmetic<__m128d> { enum { value = true }; };
54 template<> struct is_arithmetic<Packet4i> { enum { value = true }; };
55 template<> struct is_arithmetic<Packet16b> { enum { value = true }; };
56 
57 template<int p, int q, int r, int s>
58 struct shuffle_mask{
59  enum { mask = (s)<<6|(r)<<4|(q)<<2|(p) };
60 };
61 
62 // TODO: change the implementation of all swizzle* ops from macro to template,
63 #define vec4f_swizzle1(v,p,q,r,s) \
64  Packet4f(_mm_castsi128_ps(_mm_shuffle_epi32( _mm_castps_si128(v), (shuffle_mask<p,q,r,s>::mask))))
65 
66 #define vec4i_swizzle1(v,p,q,r,s) \
67  Packet4i(_mm_shuffle_epi32( v, (shuffle_mask<p,q,r,s>::mask)))
68 
69 #define vec2d_swizzle1(v,p,q) \
70  Packet2d(_mm_castsi128_pd(_mm_shuffle_epi32( _mm_castpd_si128(v), (shuffle_mask<2*p,2*p+1,2*q,2*q+1>::mask))))
71 
72 #define vec4f_swizzle2(a,b,p,q,r,s) \
73  Packet4f(_mm_shuffle_ps( (a), (b), (shuffle_mask<p,q,r,s>::mask)))
74 
75 #define vec4i_swizzle2(a,b,p,q,r,s) \
76  Packet4i(_mm_castps_si128( (_mm_shuffle_ps( _mm_castsi128_ps(a), _mm_castsi128_ps(b), (shuffle_mask<p,q,r,s>::mask)))))
77 
78 EIGEN_STRONG_INLINE Packet4f vec4f_movelh(const Packet4f& a, const Packet4f& b)
79 {
80  return Packet4f(_mm_movelh_ps(a,b));
81 }
82 EIGEN_STRONG_INLINE Packet4f vec4f_movehl(const Packet4f& a, const Packet4f& b)
83 {
84  return Packet4f(_mm_movehl_ps(a,b));
85 }
86 EIGEN_STRONG_INLINE Packet4f vec4f_unpacklo(const Packet4f& a, const Packet4f& b)
87 {
88  return Packet4f(_mm_unpacklo_ps(a,b));
89 }
90 EIGEN_STRONG_INLINE Packet4f vec4f_unpackhi(const Packet4f& a, const Packet4f& b)
91 {
92  return Packet4f(_mm_unpackhi_ps(a,b));
93 }
94 #define vec4f_duplane(a,p) \
95  vec4f_swizzle2(a,a,p,p,p,p)
96 
97 #define vec2d_swizzle2(a,b,mask) \
98  Packet2d(_mm_shuffle_pd(a,b,mask))
99 
100 EIGEN_STRONG_INLINE Packet2d vec2d_unpacklo(const Packet2d& a, const Packet2d& b)
101 {
102  return Packet2d(_mm_unpacklo_pd(a,b));
103 }
104 EIGEN_STRONG_INLINE Packet2d vec2d_unpackhi(const Packet2d& a, const Packet2d& b)
105 {
106  return Packet2d(_mm_unpackhi_pd(a,b));
107 }
108 #define vec2d_duplane(a,p) \
109  vec2d_swizzle2(a,a,(p<<1)|p)
110 
111 #define EIGEN_DECLARE_CONST_Packet4f(NAME,X) \
112  const Packet4f p4f_##NAME = pset1<Packet4f>(X)
113 
114 #define EIGEN_DECLARE_CONST_Packet2d(NAME,X) \
115  const Packet2d p2d_##NAME = pset1<Packet2d>(X)
116 
117 #define EIGEN_DECLARE_CONST_Packet4f_FROM_INT(NAME,X) \
118  const Packet4f p4f_##NAME = pset1frombits<Packet4f>(X)
119 
120 #define EIGEN_DECLARE_CONST_Packet4i(NAME,X) \
121  const Packet4i p4i_##NAME = pset1<Packet4i>(X)
122 
123 
124 // Use the packet_traits defined in AVX/PacketMath.h instead if we're going
125 // to leverage AVX instructions.
126 #ifndef EIGEN_VECTORIZE_AVX
127 template <>
128 struct packet_traits<float> : default_packet_traits {
129  typedef Packet4f type;
130  typedef Packet4f half;
131  enum {
132  Vectorizable = 1,
133  AlignedOnScalar = 1,
134  size = 4,
135  HasHalfPacket = 0,
136 
137  HasCmp = 1,
138  HasDiv = 1,
139  HasReciprocal = EIGEN_FAST_MATH,
140  HasSin = EIGEN_FAST_MATH,
141  HasCos = EIGEN_FAST_MATH,
142  HasLog = 1,
143  HasLog1p = 1,
144  HasExpm1 = 1,
145  HasNdtri = 1,
146  HasExp = 1,
147  HasBessel = 1,
148  HasSqrt = 1,
149  HasRsqrt = 1,
150  HasTanh = EIGEN_FAST_MATH,
151  HasErf = EIGEN_FAST_MATH,
152  HasBlend = 1,
153  HasCeil = 1,
154  HasFloor = 1,
155 #ifdef EIGEN_VECTORIZE_SSE4_1
156  HasRound = 1,
157 #endif
158  HasRint = 1
159  };
160 };
161 template <>
162 struct packet_traits<double> : default_packet_traits {
163  typedef Packet2d type;
164  typedef Packet2d half;
165  enum {
166  Vectorizable = 1,
167  AlignedOnScalar = 1,
168  size=2,
169  HasHalfPacket = 0,
170 
171  HasCmp = 1,
172  HasDiv = 1,
173  HasLog = 1,
174  HasExp = 1,
175  HasSqrt = 1,
176  HasRsqrt = 1,
177  HasBlend = 1,
178  HasFloor = 1,
179  HasCeil = 1,
180 #ifdef EIGEN_VECTORIZE_SSE4_1
181  HasRound = 1,
182 #endif
183  HasRint = 1
184  };
185 };
186 template<> struct packet_traits<int> : default_packet_traits
187 {
188  typedef Packet4i type;
189  typedef Packet4i half;
190  enum {
191  Vectorizable = 1,
192  AlignedOnScalar = 1,
193  size=4,
194 
195  HasShift = 1,
196  HasBlend = 1
197  };
198 };
199 #endif
200 template<> struct packet_traits<bool> : default_packet_traits
201 {
202  typedef Packet16b type;
203  typedef Packet16b half;
204  enum {
205  Vectorizable = 1,
206  AlignedOnScalar = 1,
207  HasHalfPacket = 0,
208  size=16,
209 
210  HasAdd = 1,
211  HasSub = 1,
212  HasShift = 0,
213  HasMul = 1,
214  HasNegate = 1,
215  HasAbs = 0,
216  HasAbs2 = 0,
217  HasMin = 0,
218  HasMax = 0,
219  HasConj = 0,
220  HasSqrt = 1
221  };
222 };
223 
224 template<> struct unpacket_traits<Packet4f> {
225  typedef float type;
226  typedef Packet4f half;
227  typedef Packet4i integer_packet;
228  enum {size=4, alignment=Aligned16, vectorizable=true, masked_load_available=false, masked_store_available=false};
229 };
230 template<> struct unpacket_traits<Packet2d> {
231  typedef double type;
232  typedef Packet2d half;
233  enum {size=2, alignment=Aligned16, vectorizable=true, masked_load_available=false, masked_store_available=false};
234 };
235 template<> struct unpacket_traits<Packet4i> {
236  typedef int type;
237  typedef Packet4i half;
238  enum {size=4, alignment=Aligned16, vectorizable=true, masked_load_available=false, masked_store_available=false};
239 };
240 template<> struct unpacket_traits<Packet16b> {
241  typedef bool type;
242  typedef Packet16b half;
243  enum {size=16, alignment=Aligned16, vectorizable=true, masked_load_available=false, masked_store_available=false};
244 };
245 
246 #ifndef EIGEN_VECTORIZE_AVX
247 template<> struct scalar_div_cost<float,true> { enum { value = 7 }; };
248 template<> struct scalar_div_cost<double,true> { enum { value = 8 }; };
249 #endif
250 
251 template<> EIGEN_STRONG_INLINE Packet4f pset1<Packet4f>(const float& from) { return _mm_set_ps1(from); }
252 template<> EIGEN_STRONG_INLINE Packet2d pset1<Packet2d>(const double& from) { return _mm_set1_pd(from); }
253 template<> EIGEN_STRONG_INLINE Packet4i pset1<Packet4i>(const int& from) { return _mm_set1_epi32(from); }
254 template<> EIGEN_STRONG_INLINE Packet16b pset1<Packet16b>(const bool& from) { return _mm_set1_epi8(static_cast<char>(from)); }
255 
256 template<> EIGEN_STRONG_INLINE Packet4f pset1frombits<Packet4f>(unsigned int from) { return _mm_castsi128_ps(pset1<Packet4i>(from)); }
257 template<> EIGEN_STRONG_INLINE Packet2d pset1frombits<Packet2d>(uint64_t from) { return _mm_castsi128_pd(_mm_set1_epi64x(from)); }
258 
259 template<> EIGEN_STRONG_INLINE Packet4f peven_mask(const Packet4f& /*a*/) { return _mm_castsi128_ps(_mm_set_epi32(0, -1, 0, -1)); }
260 template<> EIGEN_STRONG_INLINE Packet4i peven_mask(const Packet4i& /*a*/) { return _mm_set_epi32(0, -1, 0, -1); }
261 template<> EIGEN_STRONG_INLINE Packet2d peven_mask(const Packet2d& /*a*/) { return _mm_castsi128_pd(_mm_set_epi32(0, 0, -1, -1)); }
262 
263 template<> EIGEN_STRONG_INLINE Packet4f pzero(const Packet4f& /*a*/) { return _mm_setzero_ps(); }
264 template<> EIGEN_STRONG_INLINE Packet2d pzero(const Packet2d& /*a*/) { return _mm_setzero_pd(); }
265 template<> EIGEN_STRONG_INLINE Packet4i pzero(const Packet4i& /*a*/) { return _mm_setzero_si128(); }
266 
267 // GCC generates a shufps instruction for _mm_set1_ps/_mm_load1_ps instead of the more efficient pshufd instruction.
268 // However, using inrinsics for pset1 makes gcc to generate crappy code in some cases (see bug 203)
269 // Using inline assembly is also not an option because then gcc fails to reorder properly the instructions.
270 // Therefore, we introduced the pload1 functions to be used in product kernels for which bug 203 does not apply.
271 // Also note that with AVX, we want it to generate a vbroadcastss.
272 #if EIGEN_COMP_GNUC_STRICT && (!defined __AVX__)
273 template<> EIGEN_STRONG_INLINE Packet4f pload1<Packet4f>(const float *from) {
274  return vec4f_swizzle1(_mm_load_ss(from),0,0,0,0);
275 }
276 #endif
277 
278 template<> EIGEN_STRONG_INLINE Packet4f plset<Packet4f>(const float& a) { return _mm_add_ps(pset1<Packet4f>(a), _mm_set_ps(3,2,1,0)); }
279 template<> EIGEN_STRONG_INLINE Packet2d plset<Packet2d>(const double& a) { return _mm_add_pd(pset1<Packet2d>(a),_mm_set_pd(1,0)); }
280 template<> EIGEN_STRONG_INLINE Packet4i plset<Packet4i>(const int& a) { return _mm_add_epi32(pset1<Packet4i>(a),_mm_set_epi32(3,2,1,0)); }
281 
282 template<> EIGEN_STRONG_INLINE Packet4f padd<Packet4f>(const Packet4f& a, const Packet4f& b) { return _mm_add_ps(a,b); }
283 template<> EIGEN_STRONG_INLINE Packet2d padd<Packet2d>(const Packet2d& a, const Packet2d& b) { return _mm_add_pd(a,b); }
284 template<> EIGEN_STRONG_INLINE Packet4i padd<Packet4i>(const Packet4i& a, const Packet4i& b) { return _mm_add_epi32(a,b); }
285 
286 template<> EIGEN_STRONG_INLINE Packet16b padd<Packet16b>(const Packet16b& a, const Packet16b& b) { return _mm_or_si128(a,b); }
287 
288 template<> EIGEN_STRONG_INLINE Packet4f psub<Packet4f>(const Packet4f& a, const Packet4f& b) { return _mm_sub_ps(a,b); }
289 template<> EIGEN_STRONG_INLINE Packet2d psub<Packet2d>(const Packet2d& a, const Packet2d& b) { return _mm_sub_pd(a,b); }
290 template<> EIGEN_STRONG_INLINE Packet4i psub<Packet4i>(const Packet4i& a, const Packet4i& b) { return _mm_sub_epi32(a,b); }
291 template<> EIGEN_STRONG_INLINE Packet16b psub<Packet16b>(const Packet16b& a, const Packet16b& b) { return _mm_xor_si128(a,b); }
292 
293 template<> EIGEN_STRONG_INLINE Packet4f pxor<Packet4f>(const Packet4f& a, const Packet4f& b);
294 template<> EIGEN_STRONG_INLINE Packet4f paddsub<Packet4f>(const Packet4f& a, const Packet4f& b)
295 {
296 #ifdef EIGEN_VECTORIZE_SSE3
297  return _mm_addsub_ps(a,b);
298 #else
299  const Packet4f mask = _mm_castsi128_ps(_mm_setr_epi32(0x80000000,0x0,0x80000000,0x0));
300  return padd(a, pxor(mask, b));
301 #endif
302 }
303 
304 template<> EIGEN_STRONG_INLINE Packet2d pxor<Packet2d>(const Packet2d& , const Packet2d& );
305 template<> EIGEN_STRONG_INLINE Packet2d paddsub<Packet2d>(const Packet2d& a, const Packet2d& b)
306 {
307 #ifdef EIGEN_VECTORIZE_SSE3
308  return _mm_addsub_pd(a,b);
309 #else
310  const Packet2d mask = _mm_castsi128_pd(_mm_setr_epi32(0x0,0x80000000,0x0,0x0));
311  return padd(a, pxor(mask, b));
312 #endif
313 }
314 
315 template<> EIGEN_STRONG_INLINE Packet4f pnegate(const Packet4f& a)
316 {
317  const Packet4f mask = _mm_castsi128_ps(_mm_setr_epi32(0x80000000,0x80000000,0x80000000,0x80000000));
318  return _mm_xor_ps(a,mask);
319 }
320 template<> EIGEN_STRONG_INLINE Packet2d pnegate(const Packet2d& a)
321 {
322  const Packet2d mask = _mm_castsi128_pd(_mm_setr_epi32(0x0,0x80000000,0x0,0x80000000));
323  return _mm_xor_pd(a,mask);
324 }
325 template<> EIGEN_STRONG_INLINE Packet4i pnegate(const Packet4i& a)
326 {
327  return psub(Packet4i(_mm_setr_epi32(0,0,0,0)), a);
328 }
329 
330 template<> EIGEN_STRONG_INLINE Packet16b pnegate(const Packet16b& a)
331 {
332  return psub(pset1<Packet16b>(false), a);
333 }
334 
335 template<> EIGEN_STRONG_INLINE Packet4f pconj(const Packet4f& a) { return a; }
336 template<> EIGEN_STRONG_INLINE Packet2d pconj(const Packet2d& a) { return a; }
337 template<> EIGEN_STRONG_INLINE Packet4i pconj(const Packet4i& a) { return a; }
338 
339 template<> EIGEN_STRONG_INLINE Packet4f pmul<Packet4f>(const Packet4f& a, const Packet4f& b) { return _mm_mul_ps(a,b); }
340 template<> EIGEN_STRONG_INLINE Packet2d pmul<Packet2d>(const Packet2d& a, const Packet2d& b) { return _mm_mul_pd(a,b); }
341 template<> EIGEN_STRONG_INLINE Packet4i pmul<Packet4i>(const Packet4i& a, const Packet4i& b)
342 {
343 #ifdef EIGEN_VECTORIZE_SSE4_1
344  return _mm_mullo_epi32(a,b);
345 #else
346  // this version is slightly faster than 4 scalar products
347  return vec4i_swizzle1(
348  vec4i_swizzle2(
349  _mm_mul_epu32(a,b),
350  _mm_mul_epu32(vec4i_swizzle1(a,1,0,3,2),
351  vec4i_swizzle1(b,1,0,3,2)),
352  0,2,0,2),
353  0,2,1,3);
354 #endif
355 }
356 
357 template<> EIGEN_STRONG_INLINE Packet16b pmul<Packet16b>(const Packet16b& a, const Packet16b& b) { return _mm_and_si128(a,b); }
358 
359 template<> EIGEN_STRONG_INLINE Packet4f pdiv<Packet4f>(const Packet4f& a, const Packet4f& b) { return _mm_div_ps(a,b); }
360 template<> EIGEN_STRONG_INLINE Packet2d pdiv<Packet2d>(const Packet2d& a, const Packet2d& b) { return _mm_div_pd(a,b); }
361 
362 // for some weird raisons, it has to be overloaded for packet of integers
363 template<> EIGEN_STRONG_INLINE Packet4i pmadd(const Packet4i& a, const Packet4i& b, const Packet4i& c) { return padd(pmul(a,b), c); }
364 #ifdef EIGEN_VECTORIZE_FMA
365 template<> EIGEN_STRONG_INLINE Packet4f pmadd(const Packet4f& a, const Packet4f& b, const Packet4f& c) { return _mm_fmadd_ps(a,b,c); }
366 template<> EIGEN_STRONG_INLINE Packet2d pmadd(const Packet2d& a, const Packet2d& b, const Packet2d& c) { return _mm_fmadd_pd(a,b,c); }
367 template<> EIGEN_STRONG_INLINE Packet4f pmsub(const Packet4f& a, const Packet4f& b, const Packet4f& c) { return _mm_fmsub_ps(a,b,c); }
368 template<> EIGEN_STRONG_INLINE Packet2d pmsub(const Packet2d& a, const Packet2d& b, const Packet2d& c) { return _mm_fmsub_pd(a,b,c); }
369 template<> EIGEN_STRONG_INLINE Packet4f pnmadd(const Packet4f& a, const Packet4f& b, const Packet4f& c) { return _mm_fnmadd_ps(a,b,c); }
370 template<> EIGEN_STRONG_INLINE Packet2d pnmadd(const Packet2d& a, const Packet2d& b, const Packet2d& c) { return _mm_fnmadd_pd(a,b,c); }
371 template<> EIGEN_STRONG_INLINE Packet4f pnmsub(const Packet4f& a, const Packet4f& b, const Packet4f& c) { return _mm_fnmsub_ps(a,b,c); }
372 template<> EIGEN_STRONG_INLINE Packet2d pnmsub(const Packet2d& a, const Packet2d& b, const Packet2d& c) { return _mm_fnmsub_pd(a,b,c); }
373 #endif
374 
375 #ifdef EIGEN_VECTORIZE_SSE4_1
376 template<> EIGEN_DEVICE_FUNC inline Packet4f pselect(const Packet4f& mask, const Packet4f& a, const Packet4f& b) {
377  return _mm_blendv_ps(b,a,mask);
378 }
379 
380 template<> EIGEN_DEVICE_FUNC inline Packet4i pselect(const Packet4i& mask, const Packet4i& a, const Packet4i& b) {
381  return _mm_castps_si128(_mm_blendv_ps(_mm_castsi128_ps(b),_mm_castsi128_ps(a),_mm_castsi128_ps(mask)));
382 }
383 
384 template<> EIGEN_DEVICE_FUNC inline Packet2d pselect(const Packet2d& mask, const Packet2d& a, const Packet2d& b) { return _mm_blendv_pd(b,a,mask); }
385 
386 template<> EIGEN_DEVICE_FUNC inline Packet16b pselect(const Packet16b& mask, const Packet16b& a, const Packet16b& b) {
387  return _mm_blendv_epi8(b,a,mask);
388 }
389 #else
390 template<> EIGEN_DEVICE_FUNC inline Packet16b pselect(const Packet16b& mask, const Packet16b& a, const Packet16b& b) {
391  Packet16b a_part = _mm_and_si128(mask, a);
392  Packet16b b_part = _mm_andnot_si128(mask, b);
393  return _mm_or_si128(a_part, b_part);
394 }
395 #endif
396 
397 template<> EIGEN_STRONG_INLINE Packet4i ptrue<Packet4i>(const Packet4i& a) { return _mm_cmpeq_epi32(a, a); }
398 template<> EIGEN_STRONG_INLINE Packet16b ptrue<Packet16b>(const Packet16b& a) { return _mm_cmpeq_epi8(a, a); }
399 template<> EIGEN_STRONG_INLINE Packet4f
400 ptrue<Packet4f>(const Packet4f& a) {
401  Packet4i b = _mm_castps_si128(a);
402  return _mm_castsi128_ps(_mm_cmpeq_epi32(b, b));
403 }
404 template<> EIGEN_STRONG_INLINE Packet2d
405 ptrue<Packet2d>(const Packet2d& a) {
406  Packet4i b = _mm_castpd_si128(a);
407  return _mm_castsi128_pd(_mm_cmpeq_epi32(b, b));
408 }
409 
410 
411 template<> EIGEN_STRONG_INLINE Packet4f pand<Packet4f>(const Packet4f& a, const Packet4f& b) { return _mm_and_ps(a,b); }
412 template<> EIGEN_STRONG_INLINE Packet2d pand<Packet2d>(const Packet2d& a, const Packet2d& b) { return _mm_and_pd(a,b); }
413 template<> EIGEN_STRONG_INLINE Packet4i pand<Packet4i>(const Packet4i& a, const Packet4i& b) { return _mm_and_si128(a,b); }
414 template<> EIGEN_STRONG_INLINE Packet16b pand<Packet16b>(const Packet16b& a, const Packet16b& b) { return _mm_and_si128(a,b); }
415 
416 template<> EIGEN_STRONG_INLINE Packet4f por<Packet4f>(const Packet4f& a, const Packet4f& b) { return _mm_or_ps(a,b); }
417 template<> EIGEN_STRONG_INLINE Packet2d por<Packet2d>(const Packet2d& a, const Packet2d& b) { return _mm_or_pd(a,b); }
418 template<> EIGEN_STRONG_INLINE Packet4i por<Packet4i>(const Packet4i& a, const Packet4i& b) { return _mm_or_si128(a,b); }
419 template<> EIGEN_STRONG_INLINE Packet16b por<Packet16b>(const Packet16b& a, const Packet16b& b) { return _mm_or_si128(a,b); }
420 
421 template<> EIGEN_STRONG_INLINE Packet4f pxor<Packet4f>(const Packet4f& a, const Packet4f& b) { return _mm_xor_ps(a,b); }
422 template<> EIGEN_STRONG_INLINE Packet2d pxor<Packet2d>(const Packet2d& a, const Packet2d& b) { return _mm_xor_pd(a,b); }
423 template<> EIGEN_STRONG_INLINE Packet4i pxor<Packet4i>(const Packet4i& a, const Packet4i& b) { return _mm_xor_si128(a,b); }
424 template<> EIGEN_STRONG_INLINE Packet16b pxor<Packet16b>(const Packet16b& a, const Packet16b& b) { return _mm_xor_si128(a,b); }
425 
426 template<> EIGEN_STRONG_INLINE Packet4f pandnot<Packet4f>(const Packet4f& a, const Packet4f& b) { return _mm_andnot_ps(b,a); }
427 template<> EIGEN_STRONG_INLINE Packet2d pandnot<Packet2d>(const Packet2d& a, const Packet2d& b) { return _mm_andnot_pd(b,a); }
428 template<> EIGEN_STRONG_INLINE Packet4i pandnot<Packet4i>(const Packet4i& a, const Packet4i& b) { return _mm_andnot_si128(b,a); }
429 
430 template<> EIGEN_STRONG_INLINE Packet4f pcmp_le(const Packet4f& a, const Packet4f& b) { return _mm_cmple_ps(a,b); }
431 template<> EIGEN_STRONG_INLINE Packet4f pcmp_lt(const Packet4f& a, const Packet4f& b) { return _mm_cmplt_ps(a,b); }
432 template<> EIGEN_STRONG_INLINE Packet4f pcmp_lt_or_nan(const Packet4f& a, const Packet4f& b) { return _mm_cmpnge_ps(a,b); }
433 template<> EIGEN_STRONG_INLINE Packet4f pcmp_eq(const Packet4f& a, const Packet4f& b) { return _mm_cmpeq_ps(a,b); }
434 
435 template<> EIGEN_STRONG_INLINE Packet2d pcmp_le(const Packet2d& a, const Packet2d& b) { return _mm_cmple_pd(a,b); }
436 template<> EIGEN_STRONG_INLINE Packet2d pcmp_lt(const Packet2d& a, const Packet2d& b) { return _mm_cmplt_pd(a,b); }
437 template<> EIGEN_STRONG_INLINE Packet2d pcmp_lt_or_nan(const Packet2d& a, const Packet2d& b) { return _mm_cmpnge_pd(a,b); }
438 template<> EIGEN_STRONG_INLINE Packet2d pcmp_eq(const Packet2d& a, const Packet2d& b) { return _mm_cmpeq_pd(a,b); }
439 
440 template<> EIGEN_STRONG_INLINE Packet4i pcmp_lt(const Packet4i& a, const Packet4i& b) { return _mm_cmplt_epi32(a,b); }
441 template<> EIGEN_STRONG_INLINE Packet4i pcmp_eq(const Packet4i& a, const Packet4i& b) { return _mm_cmpeq_epi32(a,b); }
442 template<> EIGEN_STRONG_INLINE Packet16b pcmp_eq(const Packet16b& a, const Packet16b& b) { return _mm_cmpeq_epi8(a,b); }
443 template<> EIGEN_STRONG_INLINE Packet4i pcmp_le(const Packet4i& a, const Packet4i& b) { return por(pcmp_lt(a,b), pcmp_eq(a,b)); }
444 
445 template<> EIGEN_STRONG_INLINE Packet4f pmin<Packet4f>(const Packet4f& a, const Packet4f& b) {
446 #if EIGEN_COMP_GNUC_STRICT && EIGEN_COMP_GNUC < 63
447  // There appears to be a bug in GCC, by which the optimizer may
448  // flip the argument order in calls to _mm_min_ps, so we have to
449  // resort to inline ASM here. This is supposed to be fixed in gcc6.3,
450  // see also: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=72867
451  #ifdef EIGEN_VECTORIZE_AVX
452  Packet4f res;
453  asm("vminps %[a], %[b], %[res]" : [res] "=x" (res) : [a] "x" (a), [b] "x" (b));
454  #else
455  Packet4f res = b;
456  asm("minps %[a], %[res]" : [res] "+x" (res) : [a] "x" (a));
457  #endif
458  return res;
459 #else
460  // Arguments are reversed to match NaN propagation behavior of std::min.
461  return _mm_min_ps(b, a);
462 #endif
463 }
464 template<> EIGEN_STRONG_INLINE Packet2d pmin<Packet2d>(const Packet2d& a, const Packet2d& b) {
465 #if EIGEN_COMP_GNUC_STRICT && EIGEN_COMP_GNUC < 63
466  // There appears to be a bug in GCC, by which the optimizer may
467  // flip the argument order in calls to _mm_min_pd, so we have to
468  // resort to inline ASM here. This is supposed to be fixed in gcc6.3,
469  // see also: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=72867
470  #ifdef EIGEN_VECTORIZE_AVX
471  Packet2d res;
472  asm("vminpd %[a], %[b], %[res]" : [res] "=x" (res) : [a] "x" (a), [b] "x" (b));
473  #else
474  Packet2d res = b;
475  asm("minpd %[a], %[res]" : [res] "+x" (res) : [a] "x" (a));
476  #endif
477  return res;
478 #else
479  // Arguments are reversed to match NaN propagation behavior of std::min.
480  return _mm_min_pd(b, a);
481 #endif
482 }
483 template<> EIGEN_STRONG_INLINE Packet4i pmin<Packet4i>(const Packet4i& a, const Packet4i& b)
484 {
485 #ifdef EIGEN_VECTORIZE_SSE4_1
486  return _mm_min_epi32(a,b);
487 #else
488  // after some bench, this version *is* faster than a scalar implementation
489  Packet4i mask = _mm_cmplt_epi32(a,b);
490  return _mm_or_si128(_mm_and_si128(mask,a),_mm_andnot_si128(mask,b));
491 #endif
492 }
493 
494 
495 template<> EIGEN_STRONG_INLINE Packet4f pmax<Packet4f>(const Packet4f& a, const Packet4f& b) {
496 #if EIGEN_COMP_GNUC_STRICT && EIGEN_COMP_GNUC < 63
497  // There appears to be a bug in GCC, by which the optimizer may
498  // flip the argument order in calls to _mm_max_ps, so we have to
499  // resort to inline ASM here. This is supposed to be fixed in gcc6.3,
500  // see also: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=72867
501  #ifdef EIGEN_VECTORIZE_AVX
502  Packet4f res;
503  asm("vmaxps %[a], %[b], %[res]" : [res] "=x" (res) : [a] "x" (a), [b] "x" (b));
504  #else
505  Packet4f res = b;
506  asm("maxps %[a], %[res]" : [res] "+x" (res) : [a] "x" (a));
507  #endif
508  return res;
509 #else
510  // Arguments are reversed to match NaN propagation behavior of std::max.
511  return _mm_max_ps(b, a);
512 #endif
513 }
514 template<> EIGEN_STRONG_INLINE Packet2d pmax<Packet2d>(const Packet2d& a, const Packet2d& b) {
515 #if EIGEN_COMP_GNUC_STRICT && EIGEN_COMP_GNUC < 63
516  // There appears to be a bug in GCC, by which the optimizer may
517  // flip the argument order in calls to _mm_max_pd, so we have to
518  // resort to inline ASM here. This is supposed to be fixed in gcc6.3,
519  // see also: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=72867
520  #ifdef EIGEN_VECTORIZE_AVX
521  Packet2d res;
522  asm("vmaxpd %[a], %[b], %[res]" : [res] "=x" (res) : [a] "x" (a), [b] "x" (b));
523  #else
524  Packet2d res = b;
525  asm("maxpd %[a], %[res]" : [res] "+x" (res) : [a] "x" (a));
526  #endif
527  return res;
528 #else
529  // Arguments are reversed to match NaN propagation behavior of std::max.
530  return _mm_max_pd(b, a);
531 #endif
532 }
533 template<> EIGEN_STRONG_INLINE Packet4i pmax<Packet4i>(const Packet4i& a, const Packet4i& b)
534 {
535 #ifdef EIGEN_VECTORIZE_SSE4_1
536  return _mm_max_epi32(a,b);
537 #else
538  // after some bench, this version *is* faster than a scalar implementation
539  Packet4i mask = _mm_cmpgt_epi32(a,b);
540  return _mm_or_si128(_mm_and_si128(mask,a),_mm_andnot_si128(mask,b));
541 #endif
542 }
543 
544 template <typename Packet, typename Op>
545 EIGEN_STRONG_INLINE Packet pminmax_propagate_numbers(const Packet& a, const Packet& b, Op op) {
546  // In this implementation, we take advantage of the fact that pmin/pmax for SSE
547  // always return a if either a or b is NaN.
548  Packet not_nan_mask_a = pcmp_eq(a, a);
549  Packet m = op(a, b);
550  return pselect<Packet>(not_nan_mask_a, m, b);
551 }
552 
553 template <typename Packet, typename Op>
554 EIGEN_STRONG_INLINE Packet pminmax_propagate_nan(const Packet& a, const Packet& b, Op op) {
555  // In this implementation, we take advantage of the fact that pmin/pmax for SSE
556  // always return a if either a or b is NaN.
557  Packet not_nan_mask_a = pcmp_eq(a, a);
558  Packet m = op(b, a);
559  return pselect<Packet>(not_nan_mask_a, m, a);
560 }
561 
562 // Add specializations for min/max with prescribed NaN progation.
563 template<>
564 EIGEN_STRONG_INLINE Packet4f pmin<PropagateNumbers, Packet4f>(const Packet4f& a, const Packet4f& b) {
565  return pminmax_propagate_numbers(a, b, pmin<Packet4f>);
566 }
567 template<>
568 EIGEN_STRONG_INLINE Packet2d pmin<PropagateNumbers, Packet2d>(const Packet2d& a, const Packet2d& b) {
569  return pminmax_propagate_numbers(a, b, pmin<Packet2d>);
570 }
571 template<>
572 EIGEN_STRONG_INLINE Packet4f pmax<PropagateNumbers, Packet4f>(const Packet4f& a, const Packet4f& b) {
573  return pminmax_propagate_numbers(a, b, pmax<Packet4f>);
574 }
575 template<>
576 EIGEN_STRONG_INLINE Packet2d pmax<PropagateNumbers, Packet2d>(const Packet2d& a, const Packet2d& b) {
577  return pminmax_propagate_numbers(a, b, pmax<Packet2d>);
578 }
579 template<>
580 EIGEN_STRONG_INLINE Packet4f pmin<PropagateNaN, Packet4f>(const Packet4f& a, const Packet4f& b) {
581  return pminmax_propagate_nan(a, b, pmin<Packet4f>);
582 }
583 template<>
584 EIGEN_STRONG_INLINE Packet2d pmin<PropagateNaN, Packet2d>(const Packet2d& a, const Packet2d& b) {
585  return pminmax_propagate_nan(a, b, pmin<Packet2d>);
586 }
587 template<>
588 EIGEN_STRONG_INLINE Packet4f pmax<PropagateNaN, Packet4f>(const Packet4f& a, const Packet4f& b) {
589  return pminmax_propagate_nan(a, b, pmax<Packet4f>);
590 }
591 template<>
592 EIGEN_STRONG_INLINE Packet2d pmax<PropagateNaN, Packet2d>(const Packet2d& a, const Packet2d& b) {
593  return pminmax_propagate_nan(a, b, pmax<Packet2d>);
594 }
595 
596 template<int N> EIGEN_STRONG_INLINE Packet4i parithmetic_shift_right(const Packet4i& a) { return _mm_srai_epi32(a,N); }
597 template<int N> EIGEN_STRONG_INLINE Packet4i plogical_shift_right (const Packet4i& a) { return _mm_srli_epi32(a,N); }
598 template<int N> EIGEN_STRONG_INLINE Packet4i plogical_shift_left (const Packet4i& a) { return _mm_slli_epi32(a,N); }
599 
600 template<> EIGEN_STRONG_INLINE Packet4f pabs(const Packet4f& a)
601 {
602  const Packet4f mask = _mm_castsi128_ps(_mm_setr_epi32(0x7FFFFFFF,0x7FFFFFFF,0x7FFFFFFF,0x7FFFFFFF));
603  return _mm_and_ps(a,mask);
604 }
605 template<> EIGEN_STRONG_INLINE Packet2d pabs(const Packet2d& a)
606 {
607  const Packet2d mask = _mm_castsi128_pd(_mm_setr_epi32(0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF));
608  return _mm_and_pd(a,mask);
609 }
610 template<> EIGEN_STRONG_INLINE Packet4i pabs(const Packet4i& a)
611 {
612  #ifdef EIGEN_VECTORIZE_SSSE3
613  return _mm_abs_epi32(a);
614  #else
615  Packet4i aux = _mm_srai_epi32(a,31);
616  return _mm_sub_epi32(_mm_xor_si128(a,aux),aux);
617  #endif
618 }
619 
620 #ifdef EIGEN_VECTORIZE_SSE4_1
621 template<> EIGEN_STRONG_INLINE Packet4f pround<Packet4f>(const Packet4f& a)
622 {
623  // Unfortunately _mm_round_ps doesn't have a rounding mode to implement numext::round.
624  const Packet4f mask = pset1frombits<Packet4f>(0x80000000u);
625  const Packet4f prev0dot5 = pset1frombits<Packet4f>(0x3EFFFFFFu);
626  return _mm_round_ps(padd(por(pand(a, mask), prev0dot5), a), _MM_FROUND_TO_ZERO);
627 }
628 
629 template<> EIGEN_STRONG_INLINE Packet2d pround<Packet2d>(const Packet2d& a)
630 {
631  const Packet2d mask = _mm_castsi128_pd(_mm_set_epi64x(0x8000000000000000ull, 0x8000000000000000ull));
632  const Packet2d prev0dot5 = _mm_castsi128_pd(_mm_set_epi64x(0x3FDFFFFFFFFFFFFFull, 0x3FDFFFFFFFFFFFFFull));
633  return _mm_round_pd(padd(por(pand(a, mask), prev0dot5), a), _MM_FROUND_TO_ZERO);
634 }
635 
636 template<> EIGEN_STRONG_INLINE Packet4f print<Packet4f>(const Packet4f& a) { return _mm_round_ps(a, _MM_FROUND_CUR_DIRECTION); }
637 template<> EIGEN_STRONG_INLINE Packet2d print<Packet2d>(const Packet2d& a) { return _mm_round_pd(a, _MM_FROUND_CUR_DIRECTION); }
638 
639 template<> EIGEN_STRONG_INLINE Packet4f pceil<Packet4f>(const Packet4f& a) { return _mm_ceil_ps(a); }
640 template<> EIGEN_STRONG_INLINE Packet2d pceil<Packet2d>(const Packet2d& a) { return _mm_ceil_pd(a); }
641 
642 template<> EIGEN_STRONG_INLINE Packet4f pfloor<Packet4f>(const Packet4f& a) { return _mm_floor_ps(a); }
643 template<> EIGEN_STRONG_INLINE Packet2d pfloor<Packet2d>(const Packet2d& a) { return _mm_floor_pd(a); }
644 #else
645 template<> EIGEN_STRONG_INLINE Packet4f print(const Packet4f& a) {
646  // Adds and subtracts signum(a) * 2^23 to force rounding.
647  const Packet4f limit = pset1<Packet4f>(static_cast<float>(1<<23));
648  const Packet4f abs_a = pabs(a);
649  Packet4f r = padd(abs_a, limit);
650  // Don't compile-away addition and subtraction.
651  EIGEN_OPTIMIZATION_BARRIER(r);
652  r = psub(r, limit);
653  // If greater than limit, simply return a. Otherwise, account for sign.
654  r = pselect(pcmp_lt(abs_a, limit),
655  pselect(pcmp_lt(a, pzero(a)), pnegate(r), r), a);
656  return r;
657 }
658 
659 template<> EIGEN_STRONG_INLINE Packet2d print(const Packet2d& a) {
660  // Adds and subtracts signum(a) * 2^52 to force rounding.
661  const Packet2d limit = pset1<Packet2d>(static_cast<double>(1ull<<52));
662  const Packet2d abs_a = pabs(a);
663  Packet2d r = padd(abs_a, limit);
664  // Don't compile-away addition and subtraction.
665  EIGEN_OPTIMIZATION_BARRIER(r);
666  r = psub(r, limit);
667  // If greater than limit, simply return a. Otherwise, account for sign.
668  r = pselect(pcmp_lt(abs_a, limit),
669  pselect(pcmp_lt(a, pzero(a)), pnegate(r), r), a);
670  return r;
671 }
672 
673 template<> EIGEN_STRONG_INLINE Packet4f pfloor<Packet4f>(const Packet4f& a)
674 {
675  const Packet4f cst_1 = pset1<Packet4f>(1.0f);
676  Packet4f tmp = print<Packet4f>(a);
677  // If greater, subtract one.
678  Packet4f mask = _mm_cmpgt_ps(tmp, a);
679  mask = pand(mask, cst_1);
680  return psub(tmp, mask);
681 }
682 
683 template<> EIGEN_STRONG_INLINE Packet2d pfloor<Packet2d>(const Packet2d& a)
684 {
685  const Packet2d cst_1 = pset1<Packet2d>(1.0);
686  Packet2d tmp = print<Packet2d>(a);
687  // If greater, subtract one.
688  Packet2d mask = _mm_cmpgt_pd(tmp, a);
689  mask = pand(mask, cst_1);
690  return psub(tmp, mask);
691 }
692 
693 template<> EIGEN_STRONG_INLINE Packet4f pceil<Packet4f>(const Packet4f& a)
694 {
695  const Packet4f cst_1 = pset1<Packet4f>(1.0f);
696  Packet4f tmp = print<Packet4f>(a);
697  // If smaller, add one.
698  Packet4f mask = _mm_cmplt_ps(tmp, a);
699  mask = pand(mask, cst_1);
700  return padd(tmp, mask);
701 }
702 
703 template<> EIGEN_STRONG_INLINE Packet2d pceil<Packet2d>(const Packet2d& a)
704 {
705  const Packet2d cst_1 = pset1<Packet2d>(1.0);
706  Packet2d tmp = print<Packet2d>(a);
707  // If smaller, add one.
708  Packet2d mask = _mm_cmplt_pd(tmp, a);
709  mask = pand(mask, cst_1);
710  return padd(tmp, mask);
711 }
712 #endif
713 
714 template<> EIGEN_STRONG_INLINE Packet4f pload<Packet4f>(const float* from) { EIGEN_DEBUG_ALIGNED_LOAD return _mm_load_ps(from); }
715 template<> EIGEN_STRONG_INLINE Packet2d pload<Packet2d>(const double* from) { EIGEN_DEBUG_ALIGNED_LOAD return _mm_load_pd(from); }
716 template<> EIGEN_STRONG_INLINE Packet4i pload<Packet4i>(const int* from) { EIGEN_DEBUG_ALIGNED_LOAD return _mm_load_si128(reinterpret_cast<const __m128i*>(from)); }
717 template<> EIGEN_STRONG_INLINE Packet16b pload<Packet16b>(const bool* from) { EIGEN_DEBUG_ALIGNED_LOAD return _mm_load_si128(reinterpret_cast<const __m128i*>(from)); }
718 
719 #if EIGEN_COMP_MSVC
720  template<> EIGEN_STRONG_INLINE Packet4f ploadu<Packet4f>(const float* from) {
721  EIGEN_DEBUG_UNALIGNED_LOAD
722  return _mm_loadu_ps(from);
723  }
724 #else
725 // NOTE: with the code below, MSVC's compiler crashes!
726 
727 template<> EIGEN_STRONG_INLINE Packet4f ploadu<Packet4f>(const float* from)
728 {
729  EIGEN_DEBUG_UNALIGNED_LOAD
730  return _mm_loadu_ps(from);
731 }
732 #endif
733 
734 template<> EIGEN_STRONG_INLINE Packet2d ploadu<Packet2d>(const double* from)
735 {
736  EIGEN_DEBUG_UNALIGNED_LOAD
737  return _mm_loadu_pd(from);
738 }
739 template<> EIGEN_STRONG_INLINE Packet4i ploadu<Packet4i>(const int* from)
740 {
741  EIGEN_DEBUG_UNALIGNED_LOAD
742  return _mm_loadu_si128(reinterpret_cast<const __m128i*>(from));
743 }
744 template<> EIGEN_STRONG_INLINE Packet16b ploadu<Packet16b>(const bool* from) {
745  EIGEN_DEBUG_UNALIGNED_LOAD
746  return _mm_loadu_si128(reinterpret_cast<const __m128i*>(from));
747 }
748 
749 
750 template<> EIGEN_STRONG_INLINE Packet4f ploaddup<Packet4f>(const float* from)
751 {
752  return vec4f_swizzle1(_mm_castpd_ps(_mm_load_sd(reinterpret_cast<const double*>(from))), 0, 0, 1, 1);
753 }
754 template<> EIGEN_STRONG_INLINE Packet2d ploaddup<Packet2d>(const double* from)
755 { return pset1<Packet2d>(from[0]); }
756 template<> EIGEN_STRONG_INLINE Packet4i ploaddup<Packet4i>(const int* from)
757 {
758  Packet4i tmp;
759  tmp = _mm_loadl_epi64(reinterpret_cast<const __m128i*>(from));
760  return vec4i_swizzle1(tmp, 0, 0, 1, 1);
761 }
762 
763 // Loads 8 bools from memory and returns the packet
764 // {b0, b0, b1, b1, b2, b2, b3, b3, b4, b4, b5, b5, b6, b6, b7, b7}
765 template<> EIGEN_STRONG_INLINE Packet16b ploaddup<Packet16b>(const bool* from)
766 {
767  __m128i tmp = _mm_castpd_si128(pload1<Packet2d>(reinterpret_cast<const double*>(from)));
768  return _mm_unpacklo_epi8(tmp, tmp);
769 }
770 
771 // Loads 4 bools from memory and returns the packet
772 // {b0, b0 b0, b0, b1, b1, b1, b1, b2, b2, b2, b2, b3, b3, b3, b3}
773 template<> EIGEN_STRONG_INLINE Packet16b
774 ploadquad<Packet16b>(const bool* from) {
775  __m128i tmp = _mm_castps_si128(pload1<Packet4f>(reinterpret_cast<const float*>(from)));
776  tmp = _mm_unpacklo_epi8(tmp, tmp);
777  return _mm_unpacklo_epi16(tmp, tmp);
778 }
779 
780 template<> EIGEN_STRONG_INLINE void pstore<float>(float* to, const Packet4f& from) { EIGEN_DEBUG_ALIGNED_STORE _mm_store_ps(to, from); }
781 template<> EIGEN_STRONG_INLINE void pstore<double>(double* to, const Packet2d& from) { EIGEN_DEBUG_ALIGNED_STORE _mm_store_pd(to, from); }
782 template<> EIGEN_STRONG_INLINE void pstore<int>(int* to, const Packet4i& from) { EIGEN_DEBUG_ALIGNED_STORE _mm_store_si128(reinterpret_cast<__m128i*>(to), from); }
783 template<> EIGEN_STRONG_INLINE void pstore<bool>(bool* to, const Packet16b& from) { EIGEN_DEBUG_ALIGNED_STORE _mm_store_si128(reinterpret_cast<__m128i*>(to), from); }
784 
785 template<> EIGEN_STRONG_INLINE void pstoreu<double>(double* to, const Packet2d& from) { EIGEN_DEBUG_UNALIGNED_STORE _mm_storeu_pd(to, from); }
786 template<> EIGEN_STRONG_INLINE void pstoreu<float>(float* to, const Packet4f& from) { EIGEN_DEBUG_UNALIGNED_STORE _mm_storeu_ps(to, from); }
787 template<> EIGEN_STRONG_INLINE void pstoreu<int>(int* to, const Packet4i& from) { EIGEN_DEBUG_UNALIGNED_STORE _mm_storeu_si128(reinterpret_cast<__m128i*>(to), from); }
788 template<> EIGEN_STRONG_INLINE void pstoreu<bool>(bool* to, const Packet16b& from) { EIGEN_DEBUG_ALIGNED_STORE _mm_storeu_si128(reinterpret_cast<__m128i*>(to), from); }
789 
790 template<> EIGEN_DEVICE_FUNC inline Packet4f pgather<float, Packet4f>(const float* from, Index stride)
791 {
792  return _mm_set_ps(from[3*stride], from[2*stride], from[1*stride], from[0*stride]);
793 }
794 template<> EIGEN_DEVICE_FUNC inline Packet2d pgather<double, Packet2d>(const double* from, Index stride)
795 {
796  return _mm_set_pd(from[1*stride], from[0*stride]);
797 }
798 template<> EIGEN_DEVICE_FUNC inline Packet4i pgather<int, Packet4i>(const int* from, Index stride)
799 {
800  return _mm_set_epi32(from[3*stride], from[2*stride], from[1*stride], from[0*stride]);
801 }
802 
803 template<> EIGEN_DEVICE_FUNC inline Packet16b pgather<bool, Packet16b>(const bool* from, Index stride)
804 {
805  return _mm_set_epi8(from[15*stride], from[14*stride], from[13*stride], from[12*stride],
806  from[11*stride], from[10*stride], from[9*stride], from[8*stride],
807  from[7*stride], from[6*stride], from[5*stride], from[4*stride],
808  from[3*stride], from[2*stride], from[1*stride], from[0*stride]);
809 }
810 
811 template<> EIGEN_DEVICE_FUNC inline void pscatter<float, Packet4f>(float* to, const Packet4f& from, Index stride)
812 {
813  to[stride*0] = _mm_cvtss_f32(from);
814  to[stride*1] = _mm_cvtss_f32(_mm_shuffle_ps(from, from, 1));
815  to[stride*2] = _mm_cvtss_f32(_mm_shuffle_ps(from, from, 2));
816  to[stride*3] = _mm_cvtss_f32(_mm_shuffle_ps(from, from, 3));
817 }
818 template<> EIGEN_DEVICE_FUNC inline void pscatter<double, Packet2d>(double* to, const Packet2d& from, Index stride)
819 {
820  to[stride*0] = _mm_cvtsd_f64(from);
821  to[stride*1] = _mm_cvtsd_f64(_mm_shuffle_pd(from, from, 1));
822 }
823 template<> EIGEN_DEVICE_FUNC inline void pscatter<int, Packet4i>(int* to, const Packet4i& from, Index stride)
824 {
825  to[stride*0] = _mm_cvtsi128_si32(from);
826  to[stride*1] = _mm_cvtsi128_si32(_mm_shuffle_epi32(from, 1));
827  to[stride*2] = _mm_cvtsi128_si32(_mm_shuffle_epi32(from, 2));
828  to[stride*3] = _mm_cvtsi128_si32(_mm_shuffle_epi32(from, 3));
829 }
830 template<> EIGEN_DEVICE_FUNC inline void pscatter<bool, Packet16b>(bool* to, const Packet16b& from, Index stride)
831 {
832  to[4*stride*0] = _mm_cvtsi128_si32(from);
833  to[4*stride*1] = _mm_cvtsi128_si32(_mm_shuffle_epi32(from, 1));
834  to[4*stride*2] = _mm_cvtsi128_si32(_mm_shuffle_epi32(from, 2));
835  to[4*stride*3] = _mm_cvtsi128_si32(_mm_shuffle_epi32(from, 3));
836 }
837 
838 
839 // some compilers might be tempted to perform multiple moves instead of using a vector path.
840 template<> EIGEN_STRONG_INLINE void pstore1<Packet4f>(float* to, const float& a)
841 {
842  Packet4f pa = _mm_set_ss(a);
843  pstore(to, Packet4f(vec4f_swizzle1(pa,0,0,0,0)));
844 }
845 // some compilers might be tempted to perform multiple moves instead of using a vector path.
846 template<> EIGEN_STRONG_INLINE void pstore1<Packet2d>(double* to, const double& a)
847 {
848  Packet2d pa = _mm_set_sd(a);
849  pstore(to, Packet2d(vec2d_swizzle1(pa,0,0)));
850 }
851 
852 #if EIGEN_COMP_PGI && EIGEN_COMP_PGI < 1900
853 typedef const void * SsePrefetchPtrType;
854 #else
855 typedef const char * SsePrefetchPtrType;
856 #endif
857 
858 #ifndef EIGEN_VECTORIZE_AVX
859 template<> EIGEN_STRONG_INLINE void prefetch<float>(const float* addr) { _mm_prefetch((SsePrefetchPtrType)(addr), _MM_HINT_T0); }
860 template<> EIGEN_STRONG_INLINE void prefetch<double>(const double* addr) { _mm_prefetch((SsePrefetchPtrType)(addr), _MM_HINT_T0); }
861 template<> EIGEN_STRONG_INLINE void prefetch<int>(const int* addr) { _mm_prefetch((SsePrefetchPtrType)(addr), _MM_HINT_T0); }
862 #endif
863 
864 #if EIGEN_COMP_MSVC_STRICT && EIGEN_OS_WIN64
865 // The temporary variable fixes an internal compilation error in vs <= 2008 and a wrong-result bug in vs 2010
866 // Direct of the struct members fixed bug #62.
867 template<> EIGEN_STRONG_INLINE float pfirst<Packet4f>(const Packet4f& a) { return a.m128_f32[0]; }
868 template<> EIGEN_STRONG_INLINE double pfirst<Packet2d>(const Packet2d& a) { return a.m128d_f64[0]; }
869 template<> EIGEN_STRONG_INLINE int pfirst<Packet4i>(const Packet4i& a) { int x = _mm_cvtsi128_si32(a); return x; }
870 #elif EIGEN_COMP_MSVC_STRICT
871 // The temporary variable fixes an internal compilation error in vs <= 2008 and a wrong-result bug in vs 2010
872 template<> EIGEN_STRONG_INLINE float pfirst<Packet4f>(const Packet4f& a) { float x = _mm_cvtss_f32(a); return x; }
873 template<> EIGEN_STRONG_INLINE double pfirst<Packet2d>(const Packet2d& a) { double x = _mm_cvtsd_f64(a); return x; }
874 template<> EIGEN_STRONG_INLINE int pfirst<Packet4i>(const Packet4i& a) { int x = _mm_cvtsi128_si32(a); return x; }
875 #else
876 template<> EIGEN_STRONG_INLINE float pfirst<Packet4f>(const Packet4f& a) { return _mm_cvtss_f32(a); }
877 template<> EIGEN_STRONG_INLINE double pfirst<Packet2d>(const Packet2d& a) { return _mm_cvtsd_f64(a); }
878 template<> EIGEN_STRONG_INLINE int pfirst<Packet4i>(const Packet4i& a) { return _mm_cvtsi128_si32(a); }
879 #endif
880 template<> EIGEN_STRONG_INLINE bool pfirst<Packet16b>(const Packet16b& a) { int x = _mm_cvtsi128_si32(a); return static_cast<bool>(x & 1); }
881 
882 template<> EIGEN_STRONG_INLINE Packet4f preverse(const Packet4f& a) { return _mm_shuffle_ps(a,a,0x1B); }
883 template<> EIGEN_STRONG_INLINE Packet2d preverse(const Packet2d& a) { return _mm_shuffle_pd(a,a,0x1); }
884 template<> EIGEN_STRONG_INLINE Packet4i preverse(const Packet4i& a) { return _mm_shuffle_epi32(a,0x1B); }
885 template<> EIGEN_STRONG_INLINE Packet16b preverse(const Packet16b& a) {
886 #ifdef EIGEN_VECTORIZE_SSSE3
887  __m128i mask = _mm_set_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15);
888  return _mm_shuffle_epi8(a, mask);
889 #else
890  Packet16b tmp = _mm_shuffle_epi32(a, _MM_SHUFFLE(0, 1, 2, 3));
891  tmp = _mm_shufflehi_epi16(_mm_shufflelo_epi16(tmp, _MM_SHUFFLE(2, 3, 0, 1)), _MM_SHUFFLE(2, 3, 0, 1));
892  return _mm_or_si128(_mm_slli_epi16(tmp, 8), _mm_srli_epi16(tmp, 8));
893 #endif
894 }
895 
896 template<> EIGEN_STRONG_INLINE Packet4f pfrexp<Packet4f>(const Packet4f& a, Packet4f& exponent) {
897  return pfrexp_generic(a,exponent);
898 }
899 
900 // Extract exponent without existence of Packet2l.
901 template<>
902 EIGEN_STRONG_INLINE
903 Packet2d pfrexp_generic_get_biased_exponent(const Packet2d& a) {
904  const Packet2d cst_exp_mask = pset1frombits<Packet2d>(static_cast<uint64_t>(0x7ff0000000000000ull));
905  __m128i a_expo = _mm_srli_epi64(_mm_castpd_si128(pand(a, cst_exp_mask)), 52);
906  return _mm_cvtepi32_pd(vec4i_swizzle1(a_expo, 0, 2, 1, 3));
907 }
908 
909 template<> EIGEN_STRONG_INLINE Packet2d pfrexp<Packet2d>(const Packet2d& a, Packet2d& exponent) {
910  return pfrexp_generic(a, exponent);
911 }
912 
913 template<> EIGEN_STRONG_INLINE Packet4f pldexp<Packet4f>(const Packet4f& a, const Packet4f& exponent) {
914  return pldexp_generic(a,exponent);
915 }
916 
917 // We specialize pldexp here, since the generic implementation uses Packet2l, which is not well
918 // supported by SSE, and has more range than is needed for exponents.
919 template<> EIGEN_STRONG_INLINE Packet2d pldexp<Packet2d>(const Packet2d& a, const Packet2d& exponent) {
920  // Clamp exponent to [-2099, 2099]
921  const Packet2d max_exponent = pset1<Packet2d>(2099.0);
922  const Packet2d e = pmin(pmax(exponent, pnegate(max_exponent)), max_exponent);
923 
924  // Convert e to integer and swizzle to low-order bits.
925  const Packet4i ei = vec4i_swizzle1(_mm_cvtpd_epi32(e), 0, 3, 1, 3);
926 
927  // Split 2^e into four factors and multiply:
928  const Packet4i bias = _mm_set_epi32(0, 1023, 0, 1023);
929  Packet4i b = parithmetic_shift_right<2>(ei); // floor(e/4)
930  Packet2d c = _mm_castsi128_pd(_mm_slli_epi64(padd(b, bias), 52)); // 2^b
931  Packet2d out = pmul(pmul(pmul(a, c), c), c); // a * 2^(3b)
932  b = psub(psub(psub(ei, b), b), b); // e - 3b
933  c = _mm_castsi128_pd(_mm_slli_epi64(padd(b, bias), 52)); // 2^(e - 3b)
934  out = pmul(out, c); // a * 2^e
935  return out;
936 }
937 
938 // with AVX, the default implementations based on pload1 are faster
939 #ifndef __AVX__
940 template<> EIGEN_STRONG_INLINE void
941 pbroadcast4<Packet4f>(const float *a,
942  Packet4f& a0, Packet4f& a1, Packet4f& a2, Packet4f& a3)
943 {
944  a3 = pload<Packet4f>(a);
945  a0 = vec4f_swizzle1(a3, 0,0,0,0);
946  a1 = vec4f_swizzle1(a3, 1,1,1,1);
947  a2 = vec4f_swizzle1(a3, 2,2,2,2);
948  a3 = vec4f_swizzle1(a3, 3,3,3,3);
949 }
950 template<> EIGEN_STRONG_INLINE void
951 pbroadcast4<Packet2d>(const double *a,
952  Packet2d& a0, Packet2d& a1, Packet2d& a2, Packet2d& a3)
953 {
954 #ifdef EIGEN_VECTORIZE_SSE3
955  a0 = _mm_loaddup_pd(a+0);
956  a1 = _mm_loaddup_pd(a+1);
957  a2 = _mm_loaddup_pd(a+2);
958  a3 = _mm_loaddup_pd(a+3);
959 #else
960  a1 = pload<Packet2d>(a);
961  a0 = vec2d_swizzle1(a1, 0,0);
962  a1 = vec2d_swizzle1(a1, 1,1);
963  a3 = pload<Packet2d>(a+2);
964  a2 = vec2d_swizzle1(a3, 0,0);
965  a3 = vec2d_swizzle1(a3, 1,1);
966 #endif
967 }
968 #endif
969 
970 EIGEN_STRONG_INLINE void punpackp(Packet4f* vecs)
971 {
972  vecs[1] = _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(vecs[0]), 0x55));
973  vecs[2] = _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(vecs[0]), 0xAA));
974  vecs[3] = _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(vecs[0]), 0xFF));
975  vecs[0] = _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(vecs[0]), 0x00));
976 }
977 
978 template<> EIGEN_STRONG_INLINE float predux<Packet4f>(const Packet4f& a)
979 {
980  // Disable SSE3 _mm_hadd_pd that is extremely slow on all existing Intel's architectures
981  // (from Nehalem to Haswell)
982 // #ifdef EIGEN_VECTORIZE_SSE3
983 // Packet4f tmp = _mm_add_ps(a, vec4f_swizzle1(a,2,3,2,3));
984 // return pfirst<Packet4f>(_mm_hadd_ps(tmp, tmp));
985 // #else
986  Packet4f tmp = _mm_add_ps(a, _mm_movehl_ps(a,a));
987  return pfirst<Packet4f>(_mm_add_ss(tmp, _mm_shuffle_ps(tmp,tmp, 1)));
988 // #endif
989 }
990 
991 template<> EIGEN_STRONG_INLINE double predux<Packet2d>(const Packet2d& a)
992 {
993  // Disable SSE3 _mm_hadd_pd that is extremely slow on all existing Intel's architectures
994  // (from Nehalem to Haswell)
995 // #ifdef EIGEN_VECTORIZE_SSE3
996 // return pfirst<Packet2d>(_mm_hadd_pd(a, a));
997 // #else
998  return pfirst<Packet2d>(_mm_add_sd(a, _mm_unpackhi_pd(a,a)));
999 // #endif
1000 }
1001 
1002 #ifdef EIGEN_VECTORIZE_SSSE3
1003 template<> EIGEN_STRONG_INLINE int predux<Packet4i>(const Packet4i& a)
1004 {
1005  Packet4i tmp0 = _mm_hadd_epi32(a,a);
1006  return pfirst<Packet4i>(_mm_hadd_epi32(tmp0,tmp0));
1007 }
1008 
1009 #else
1010 template<> EIGEN_STRONG_INLINE int predux<Packet4i>(const Packet4i& a)
1011 {
1012  Packet4i tmp = _mm_add_epi32(a, _mm_unpackhi_epi64(a,a));
1013  return pfirst(tmp) + pfirst<Packet4i>(_mm_shuffle_epi32(tmp, 1));
1014 }
1015 #endif
1016 
1017 template<> EIGEN_STRONG_INLINE bool predux<Packet16b>(const Packet16b& a) {
1018  Packet4i tmp = _mm_or_si128(a, _mm_unpackhi_epi64(a,a));
1019  return (pfirst(tmp) != 0) || (pfirst<Packet4i>(_mm_shuffle_epi32(tmp, 1)) != 0);
1020 }
1021 
1022 // Other reduction functions:
1023 
1024 
1025 // mul
1026 template<> EIGEN_STRONG_INLINE float predux_mul<Packet4f>(const Packet4f& a)
1027 {
1028  Packet4f tmp = _mm_mul_ps(a, _mm_movehl_ps(a,a));
1029  return pfirst<Packet4f>(_mm_mul_ss(tmp, _mm_shuffle_ps(tmp,tmp, 1)));
1030 }
1031 template<> EIGEN_STRONG_INLINE double predux_mul<Packet2d>(const Packet2d& a)
1032 {
1033  return pfirst<Packet2d>(_mm_mul_sd(a, _mm_unpackhi_pd(a,a)));
1034 }
1035 template<> EIGEN_STRONG_INLINE int predux_mul<Packet4i>(const Packet4i& a)
1036 {
1037  // after some experiments, it is seems this is the fastest way to implement it
1038  // for GCC (eg., reusing pmul is very slow !)
1039  // TODO try to call _mm_mul_epu32 directly
1040  EIGEN_ALIGN16 int aux[4];
1041  pstore(aux, a);
1042  return (aux[0] * aux[1]) * (aux[2] * aux[3]);
1043 }
1044 
1045 template<> EIGEN_STRONG_INLINE bool predux_mul<Packet16b>(const Packet16b& a) {
1046  Packet4i tmp = _mm_and_si128(a, _mm_unpackhi_epi64(a,a));
1047  return ((pfirst<Packet4i>(tmp) == 0x01010101) &&
1048  (pfirst<Packet4i>(_mm_shuffle_epi32(tmp, 1)) == 0x01010101));
1049 }
1050 
1051 // min
1052 template<> EIGEN_STRONG_INLINE float predux_min<Packet4f>(const Packet4f& a)
1053 {
1054  Packet4f tmp = _mm_min_ps(a, _mm_movehl_ps(a,a));
1055  return pfirst<Packet4f>(_mm_min_ss(tmp, _mm_shuffle_ps(tmp,tmp, 1)));
1056 }
1057 template<> EIGEN_STRONG_INLINE double predux_min<Packet2d>(const Packet2d& a)
1058 {
1059  return pfirst<Packet2d>(_mm_min_sd(a, _mm_unpackhi_pd(a,a)));
1060 }
1061 template<> EIGEN_STRONG_INLINE int predux_min<Packet4i>(const Packet4i& a)
1062 {
1063 #ifdef EIGEN_VECTORIZE_SSE4_1
1064  Packet4i tmp = _mm_min_epi32(a, _mm_shuffle_epi32(a, _MM_SHUFFLE(0,0,3,2)));
1065  return pfirst<Packet4i>(_mm_min_epi32(tmp,_mm_shuffle_epi32(tmp, 1)));
1066 #else
1067  // after some experiments, it is seems this is the fastest way to implement it
1068  // for GCC (eg., it does not like using std::min after the pstore !!)
1069  EIGEN_ALIGN16 int aux[4];
1070  pstore(aux, a);
1071  int aux0 = aux[0]<aux[1] ? aux[0] : aux[1];
1072  int aux2 = aux[2]<aux[3] ? aux[2] : aux[3];
1073  return aux0<aux2 ? aux0 : aux2;
1074 #endif // EIGEN_VECTORIZE_SSE4_1
1075 }
1076 
1077 // max
1078 template<> EIGEN_STRONG_INLINE float predux_max<Packet4f>(const Packet4f& a)
1079 {
1080  Packet4f tmp = _mm_max_ps(a, _mm_movehl_ps(a,a));
1081  return pfirst<Packet4f>(_mm_max_ss(tmp, _mm_shuffle_ps(tmp,tmp, 1)));
1082 }
1083 template<> EIGEN_STRONG_INLINE double predux_max<Packet2d>(const Packet2d& a)
1084 {
1085  return pfirst<Packet2d>(_mm_max_sd(a, _mm_unpackhi_pd(a,a)));
1086 }
1087 template<> EIGEN_STRONG_INLINE int predux_max<Packet4i>(const Packet4i& a)
1088 {
1089 #ifdef EIGEN_VECTORIZE_SSE4_1
1090  Packet4i tmp = _mm_max_epi32(a, _mm_shuffle_epi32(a, _MM_SHUFFLE(0,0,3,2)));
1091  return pfirst<Packet4i>(_mm_max_epi32(tmp,_mm_shuffle_epi32(tmp, 1)));
1092 #else
1093  // after some experiments, it is seems this is the fastest way to implement it
1094  // for GCC (eg., it does not like using std::min after the pstore !!)
1095  EIGEN_ALIGN16 int aux[4];
1096  pstore(aux, a);
1097  int aux0 = aux[0]>aux[1] ? aux[0] : aux[1];
1098  int aux2 = aux[2]>aux[3] ? aux[2] : aux[3];
1099  return aux0>aux2 ? aux0 : aux2;
1100 #endif // EIGEN_VECTORIZE_SSE4_1
1101 }
1102 
1103 // not needed yet
1104 // template<> EIGEN_STRONG_INLINE bool predux_all(const Packet4f& x)
1105 // {
1106 // return _mm_movemask_ps(x) == 0xF;
1107 // }
1108 
1109 template<> EIGEN_STRONG_INLINE bool predux_any(const Packet4f& x)
1110 {
1111  return _mm_movemask_ps(x) != 0x0;
1112 }
1113 
1114 EIGEN_DEVICE_FUNC inline void
1115 ptranspose(PacketBlock<Packet4f,4>& kernel) {
1116  _MM_TRANSPOSE4_PS(kernel.packet[0], kernel.packet[1], kernel.packet[2], kernel.packet[3]);
1117 }
1118 
1119 EIGEN_DEVICE_FUNC inline void
1120 ptranspose(PacketBlock<Packet2d,2>& kernel) {
1121  __m128d tmp = _mm_unpackhi_pd(kernel.packet[0], kernel.packet[1]);
1122  kernel.packet[0] = _mm_unpacklo_pd(kernel.packet[0], kernel.packet[1]);
1123  kernel.packet[1] = tmp;
1124 }
1125 
1126 EIGEN_DEVICE_FUNC inline void
1127 ptranspose(PacketBlock<Packet4i,4>& kernel) {
1128  __m128i T0 = _mm_unpacklo_epi32(kernel.packet[0], kernel.packet[1]);
1129  __m128i T1 = _mm_unpacklo_epi32(kernel.packet[2], kernel.packet[3]);
1130  __m128i T2 = _mm_unpackhi_epi32(kernel.packet[0], kernel.packet[1]);
1131  __m128i T3 = _mm_unpackhi_epi32(kernel.packet[2], kernel.packet[3]);
1132 
1133  kernel.packet[0] = _mm_unpacklo_epi64(T0, T1);
1134  kernel.packet[1] = _mm_unpackhi_epi64(T0, T1);
1135  kernel.packet[2] = _mm_unpacklo_epi64(T2, T3);
1136  kernel.packet[3] = _mm_unpackhi_epi64(T2, T3);
1137 }
1138 
1139 EIGEN_DEVICE_FUNC inline void
1140 ptranspose(PacketBlock<Packet16b,4>& kernel) {
1141  __m128i T0 = _mm_unpacklo_epi8(kernel.packet[0], kernel.packet[1]);
1142  __m128i T1 = _mm_unpackhi_epi8(kernel.packet[0], kernel.packet[1]);
1143  __m128i T2 = _mm_unpacklo_epi8(kernel.packet[2], kernel.packet[3]);
1144  __m128i T3 = _mm_unpackhi_epi8(kernel.packet[2], kernel.packet[3]);
1145  kernel.packet[0] = _mm_unpacklo_epi16(T0, T2);
1146  kernel.packet[1] = _mm_unpackhi_epi16(T0, T2);
1147  kernel.packet[2] = _mm_unpacklo_epi16(T1, T3);
1148  kernel.packet[3] = _mm_unpackhi_epi16(T1, T3);
1149 }
1150 
1151 EIGEN_DEVICE_FUNC inline void
1152 ptranspose(PacketBlock<Packet16b,16>& kernel) {
1153  // If we number the elements in the input thus:
1154  // kernel.packet[ 0] = {00, 01, 02, 03, 04, 05, 06, 07, 08, 09, 0a, 0b, 0c, 0d, 0e, 0f}
1155  // kernel.packet[ 1] = {10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1a, 1b, 1c, 1d, 1e, 1f}
1156  // ...
1157  // kernel.packet[15] = {f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, fa, fb, fc, fd, fe, ff},
1158  //
1159  // the desired output is:
1160  // kernel.packet[ 0] = {00, 10, 20, 30, 40, 50, 60, 70, 80, 90, a0, b0, c0, d0, e0, f0}
1161  // kernel.packet[ 1] = {01, 11, 21, 31, 41, 51, 61, 71, 81, 91, a1, b1, c1, d1, e1, f1}
1162  // ...
1163  // kernel.packet[15] = {0f, 1f, 2f, 3f, 4f, 5f, 6f, 7f, 8f, 9f, af, bf, cf, df, ef, ff},
1164  __m128i t0 = _mm_unpacklo_epi8(kernel.packet[0], kernel.packet[1]); // 00 10 01 11 02 12 03 13 04 14 05 15 06 16 07 17
1165  __m128i t1 = _mm_unpackhi_epi8(kernel.packet[0], kernel.packet[1]); // 08 18 09 19 0a 1a 0b 1b 0c 1c 0d 1d 0e 1e 0f 1f
1166  __m128i t2 = _mm_unpacklo_epi8(kernel.packet[2], kernel.packet[3]); // 20 30 21 31 22 32 ... 27 37
1167  __m128i t3 = _mm_unpackhi_epi8(kernel.packet[2], kernel.packet[3]); // 28 38 29 39 2a 3a ... 2f 3f
1168  __m128i t4 = _mm_unpacklo_epi8(kernel.packet[4], kernel.packet[5]); // 40 50 41 51 42 52 47 57
1169  __m128i t5 = _mm_unpackhi_epi8(kernel.packet[4], kernel.packet[5]); // 48 58 49 59 4a 5a
1170  __m128i t6 = _mm_unpacklo_epi8(kernel.packet[6], kernel.packet[7]);
1171  __m128i t7 = _mm_unpackhi_epi8(kernel.packet[6], kernel.packet[7]);
1172  __m128i t8 = _mm_unpacklo_epi8(kernel.packet[8], kernel.packet[9]);
1173  __m128i t9 = _mm_unpackhi_epi8(kernel.packet[8], kernel.packet[9]);
1174  __m128i ta = _mm_unpacklo_epi8(kernel.packet[10], kernel.packet[11]);
1175  __m128i tb = _mm_unpackhi_epi8(kernel.packet[10], kernel.packet[11]);
1176  __m128i tc = _mm_unpacklo_epi8(kernel.packet[12], kernel.packet[13]);
1177  __m128i td = _mm_unpackhi_epi8(kernel.packet[12], kernel.packet[13]);
1178  __m128i te = _mm_unpacklo_epi8(kernel.packet[14], kernel.packet[15]);
1179  __m128i tf = _mm_unpackhi_epi8(kernel.packet[14], kernel.packet[15]);
1180 
1181  __m128i s0 = _mm_unpacklo_epi16(t0, t2); // 00 10 20 30 01 11 21 31 02 12 22 32 03 13 23 33
1182  __m128i s1 = _mm_unpackhi_epi16(t0, t2); // 04 14 24 34
1183  __m128i s2 = _mm_unpacklo_epi16(t1, t3); // 08 18 28 38 ...
1184  __m128i s3 = _mm_unpackhi_epi16(t1, t3); // 0c 1c 2c 3c ...
1185  __m128i s4 = _mm_unpacklo_epi16(t4, t6); // 40 50 60 70 41 51 61 71 42 52 62 72 43 53 63 73
1186  __m128i s5 = _mm_unpackhi_epi16(t4, t6); // 44 54 64 74 ...
1187  __m128i s6 = _mm_unpacklo_epi16(t5, t7);
1188  __m128i s7 = _mm_unpackhi_epi16(t5, t7);
1189  __m128i s8 = _mm_unpacklo_epi16(t8, ta);
1190  __m128i s9 = _mm_unpackhi_epi16(t8, ta);
1191  __m128i sa = _mm_unpacklo_epi16(t9, tb);
1192  __m128i sb = _mm_unpackhi_epi16(t9, tb);
1193  __m128i sc = _mm_unpacklo_epi16(tc, te);
1194  __m128i sd = _mm_unpackhi_epi16(tc, te);
1195  __m128i se = _mm_unpacklo_epi16(td, tf);
1196  __m128i sf = _mm_unpackhi_epi16(td, tf);
1197 
1198  __m128i u0 = _mm_unpacklo_epi32(s0, s4); // 00 10 20 30 40 50 60 70 01 11 21 31 41 51 61 71
1199  __m128i u1 = _mm_unpackhi_epi32(s0, s4); // 02 12 22 32 42 52 62 72 03 13 23 33 43 53 63 73
1200  __m128i u2 = _mm_unpacklo_epi32(s1, s5);
1201  __m128i u3 = _mm_unpackhi_epi32(s1, s5);
1202  __m128i u4 = _mm_unpacklo_epi32(s2, s6);
1203  __m128i u5 = _mm_unpackhi_epi32(s2, s6);
1204  __m128i u6 = _mm_unpacklo_epi32(s3, s7);
1205  __m128i u7 = _mm_unpackhi_epi32(s3, s7);
1206  __m128i u8 = _mm_unpacklo_epi32(s8, sc);
1207  __m128i u9 = _mm_unpackhi_epi32(s8, sc);
1208  __m128i ua = _mm_unpacklo_epi32(s9, sd);
1209  __m128i ub = _mm_unpackhi_epi32(s9, sd);
1210  __m128i uc = _mm_unpacklo_epi32(sa, se);
1211  __m128i ud = _mm_unpackhi_epi32(sa, se);
1212  __m128i ue = _mm_unpacklo_epi32(sb, sf);
1213  __m128i uf = _mm_unpackhi_epi32(sb, sf);
1214 
1215  kernel.packet[0] = _mm_unpacklo_epi64(u0, u8);
1216  kernel.packet[1] = _mm_unpackhi_epi64(u0, u8);
1217  kernel.packet[2] = _mm_unpacklo_epi64(u1, u9);
1218  kernel.packet[3] = _mm_unpackhi_epi64(u1, u9);
1219  kernel.packet[4] = _mm_unpacklo_epi64(u2, ua);
1220  kernel.packet[5] = _mm_unpackhi_epi64(u2, ua);
1221  kernel.packet[6] = _mm_unpacklo_epi64(u3, ub);
1222  kernel.packet[7] = _mm_unpackhi_epi64(u3, ub);
1223  kernel.packet[8] = _mm_unpacklo_epi64(u4, uc);
1224  kernel.packet[9] = _mm_unpackhi_epi64(u4, uc);
1225  kernel.packet[10] = _mm_unpacklo_epi64(u5, ud);
1226  kernel.packet[11] = _mm_unpackhi_epi64(u5, ud);
1227  kernel.packet[12] = _mm_unpacklo_epi64(u6, ue);
1228  kernel.packet[13] = _mm_unpackhi_epi64(u6, ue);
1229  kernel.packet[14] = _mm_unpacklo_epi64(u7, uf);
1230  kernel.packet[15] = _mm_unpackhi_epi64(u7, uf);
1231 }
1232 
1233 template<> EIGEN_STRONG_INLINE Packet4i pblend(const Selector<4>& ifPacket, const Packet4i& thenPacket, const Packet4i& elsePacket) {
1234  const __m128i zero = _mm_setzero_si128();
1235  const __m128i select = _mm_set_epi32(ifPacket.select[3], ifPacket.select[2], ifPacket.select[1], ifPacket.select[0]);
1236  __m128i false_mask = _mm_cmpeq_epi32(select, zero);
1237 #ifdef EIGEN_VECTORIZE_SSE4_1
1238  return _mm_blendv_epi8(thenPacket, elsePacket, false_mask);
1239 #else
1240  return _mm_or_si128(_mm_andnot_si128(false_mask, thenPacket), _mm_and_si128(false_mask, elsePacket));
1241 #endif
1242 }
1243 template<> EIGEN_STRONG_INLINE Packet4f pblend(const Selector<4>& ifPacket, const Packet4f& thenPacket, const Packet4f& elsePacket) {
1244  const __m128 zero = _mm_setzero_ps();
1245  const __m128 select = _mm_set_ps(ifPacket.select[3], ifPacket.select[2], ifPacket.select[1], ifPacket.select[0]);
1246  __m128 false_mask = _mm_cmpeq_ps(select, zero);
1247 #ifdef EIGEN_VECTORIZE_SSE4_1
1248  return _mm_blendv_ps(thenPacket, elsePacket, false_mask);
1249 #else
1250  return _mm_or_ps(_mm_andnot_ps(false_mask, thenPacket), _mm_and_ps(false_mask, elsePacket));
1251 #endif
1252 }
1253 template<> EIGEN_STRONG_INLINE Packet2d pblend(const Selector<2>& ifPacket, const Packet2d& thenPacket, const Packet2d& elsePacket) {
1254  const __m128d zero = _mm_setzero_pd();
1255  const __m128d select = _mm_set_pd(ifPacket.select[1], ifPacket.select[0]);
1256  __m128d false_mask = _mm_cmpeq_pd(select, zero);
1257 #ifdef EIGEN_VECTORIZE_SSE4_1
1258  return _mm_blendv_pd(thenPacket, elsePacket, false_mask);
1259 #else
1260  return _mm_or_pd(_mm_andnot_pd(false_mask, thenPacket), _mm_and_pd(false_mask, elsePacket));
1261 #endif
1262 }
1263 
1264 // Scalar path for pmadd with FMA to ensure consistency with vectorized path.
1265 #ifdef EIGEN_VECTORIZE_FMA
1266 template<> EIGEN_STRONG_INLINE float pmadd(const float& a, const float& b, const float& c) {
1267  return ::fmaf(a,b,c);
1268 }
1269 template<> EIGEN_STRONG_INLINE double pmadd(const double& a, const double& b, const double& c) {
1270  return ::fma(a,b,c);
1271 }
1272 template<> EIGEN_STRONG_INLINE float pmsub(const float& a, const float& b, const float& c) {
1273  return ::fmaf(a,b,-c);
1274 }
1275 template<> EIGEN_STRONG_INLINE double pmsub(const double& a, const double& b, const double& c) {
1276  return ::fma(a,b,-c);
1277 }
1278 template<> EIGEN_STRONG_INLINE float pnmadd(const float& a, const float& b, const float& c) {
1279  return ::fmaf(-a,b,c);
1280 }
1281 template<> EIGEN_STRONG_INLINE double pnmadd(const double& a, const double& b, const double& c) {
1282  return ::fma(-a,b,c);
1283 }
1284 template<> EIGEN_STRONG_INLINE float pnmsub(const float& a, const float& b, const float& c) {
1285  return ::fmaf(-a,b,-c);
1286 }
1287 template<> EIGEN_STRONG_INLINE double pnmsub(const double& a, const double& b, const double& c) {
1288  return ::fma(-a,b,-c);
1289 }
1290 #endif
1291 
1292 #ifdef EIGEN_VECTORIZE_SSE4_1
1293 // Helpers for half->float and float->half conversions.
1294 // Currently only used by the AVX code.
1295 EIGEN_STRONG_INLINE __m128i half2floatsse(__m128i h) {
1296  __m128i input = _mm_cvtepu16_epi32(h);
1297 
1298  // Direct vectorization of half_to_float, C parts in the comments.
1299  __m128i shifted_exp = _mm_set1_epi32(0x7c00 << 13);
1300  // o.u = (h.x & 0x7fff) << 13; // exponent/mantissa bits
1301  __m128i ou = _mm_slli_epi32(_mm_and_si128(input, _mm_set1_epi32(0x7fff)), 13);
1302  // exp = shifted_exp & o.u; // just the exponent
1303  __m128i exp = _mm_and_si128(ou, shifted_exp);
1304  // o.u += (127 - 15) << 23;
1305  ou = _mm_add_epi32(ou, _mm_set1_epi32((127 - 15) << 23));
1306 
1307  // Inf/NaN?
1308  __m128i naninf_mask = _mm_cmpeq_epi32(exp, shifted_exp);
1309  // Inf/NaN adjust
1310  __m128i naninf_adj =
1311  _mm_and_si128(_mm_set1_epi32((128 - 16) << 23), naninf_mask);
1312  // extra exp adjust for Inf/NaN
1313  ou = _mm_add_epi32(ou, naninf_adj);
1314 
1315  // Zero/Denormal?
1316  __m128i zeroden_mask = _mm_cmpeq_epi32(exp, _mm_setzero_si128());
1317  __m128i zeroden_adj = _mm_and_si128(zeroden_mask, _mm_set1_epi32(1 << 23));
1318  // o.u += 1 << 23;
1319  ou = _mm_add_epi32(ou, zeroden_adj);
1320  // magic.u = 113 << 23
1321  __m128i magic = _mm_and_si128(zeroden_mask, _mm_set1_epi32(113 << 23));
1322  // o.f -= magic.f
1323  ou = _mm_castps_si128(
1324  _mm_sub_ps(_mm_castsi128_ps(ou), _mm_castsi128_ps(magic)));
1325 
1326  __m128i sign =
1327  _mm_slli_epi32(_mm_and_si128(input, _mm_set1_epi32(0x8000)), 16);
1328  // o.u |= (h.x & 0x8000) << 16; // sign bit
1329  ou = _mm_or_si128(ou, sign);
1330  // return o.f;
1331  // We are actually returning uint version, to make
1332  // _mm256_insertf128_si256 work.
1333  return ou;
1334 }
1335 
1336 EIGEN_STRONG_INLINE __m128i float2half(__m128 f) {
1337  __m128i o = _mm_setzero_si128();
1338 
1339  // unsigned int sign_mask = 0x80000000u;
1340  __m128i sign = _mm_set1_epi32(0x80000000u);
1341  // unsigned int sign = f.u & sign_mask;
1342  sign = _mm_and_si128(sign, _mm_castps_si128(f));
1343  // f.u ^= sign;
1344  f = _mm_xor_ps(f, _mm_castsi128_ps(sign));
1345 
1346  __m128i fu = _mm_castps_si128(f);
1347 
1348  __m128i f16max = _mm_set1_epi32((127 + 16) << 23);
1349  __m128i f32infty = _mm_set1_epi32(255 << 23);
1350  // if (f.u >= f16max.u) // result is Inf or NaN (all exponent bits set)
1351  // there is no _mm_cmpge_epi32, so use lt and swap operands
1352  __m128i infnan_mask = _mm_cmplt_epi32(f16max, _mm_castps_si128(f));
1353  __m128i inf_mask = _mm_cmpgt_epi32(_mm_castps_si128(f), f32infty);
1354  __m128i nan_mask = _mm_andnot_si128(inf_mask, infnan_mask);
1355  __m128i inf_value = _mm_and_si128(inf_mask, _mm_set1_epi32(0x7e00));
1356  __m128i nan_value = _mm_and_si128(nan_mask, _mm_set1_epi32(0x7c00));
1357  // o.x = (f.u > f32infty.u) ? 0x7e00 : 0x7c00; // NaN->qNaN and Inf->Inf
1358  __m128i naninf_value = _mm_or_si128(inf_value, nan_value);
1359 
1360  __m128i denorm_magic = _mm_set1_epi32(((127 - 15) + (23 - 10) + 1) << 23);
1361  __m128i subnorm_mask =
1362  _mm_cmplt_epi32(_mm_castps_si128(f), _mm_set1_epi32(113 << 23));
1363  // f.f += denorm_magic.f;
1364  f = _mm_add_ps(f, _mm_castsi128_ps(denorm_magic));
1365  // f.u - denorm_magic.u
1366  o = _mm_sub_epi32(_mm_castps_si128(f), denorm_magic);
1367  o = _mm_and_si128(o, subnorm_mask);
1368  // Correct result for inf/nan/zero/subnormal, 0 otherwise
1369  o = _mm_or_si128(o, naninf_value);
1370 
1371  __m128i mask = _mm_or_si128(infnan_mask, subnorm_mask);
1372  o = _mm_and_si128(o, mask);
1373 
1374  // mant_odd = (f.u >> 13) & 1;
1375  __m128i mand_odd = _mm_and_si128(_mm_srli_epi32(fu, 13), _mm_set1_epi32(0x1));
1376  // f.u += 0xc8000fffU;
1377  fu = _mm_add_epi32(fu, _mm_set1_epi32(0xc8000fffU));
1378  // f.u += mant_odd;
1379  fu = _mm_add_epi32(fu, mand_odd);
1380  fu = _mm_andnot_si128(mask, fu);
1381  // f.u >> 13
1382  fu = _mm_srli_epi32(fu, 13);
1383  o = _mm_or_si128(fu, o);
1384 
1385  // o.x |= static_cast<numext::uint16_t>(sign >> 16);
1386  o = _mm_or_si128(o, _mm_srli_epi32(sign, 16));
1387 
1388  // 16 bit values
1389  return _mm_and_si128(o, _mm_set1_epi32(0xffff));
1390 }
1391 #endif
1392 
1393 // Packet math for Eigen::half
1394 // Disable the following code since it's broken on too many platforms / compilers.
1395 //#elif defined(EIGEN_VECTORIZE_SSE) && (!EIGEN_ARCH_x86_64) && (!EIGEN_COMP_MSVC)
1396 #if 0
1397 
1398 typedef struct {
1399  __m64 x;
1400 } Packet4h;
1401 
1402 
1403 template<> struct is_arithmetic<Packet4h> { enum { value = true }; };
1404 
1405 template <>
1406 struct packet_traits<Eigen::half> : default_packet_traits {
1407  typedef Packet4h type;
1408  // There is no half-size packet for Packet4h.
1409  typedef Packet4h half;
1410  enum {
1411  Vectorizable = 1,
1412  AlignedOnScalar = 1,
1413  size = 4,
1414  HasHalfPacket = 0,
1415  HasAdd = 1,
1416  HasSub = 1,
1417  HasMul = 1,
1418  HasDiv = 1,
1419  HasNegate = 0,
1420  HasAbs = 0,
1421  HasAbs2 = 0,
1422  HasMin = 0,
1423  HasMax = 0,
1424  HasConj = 0,
1425  HasSetLinear = 0,
1426  HasSqrt = 0,
1427  HasRsqrt = 0,
1428  HasExp = 0,
1429  HasLog = 0,
1430  HasBlend = 0
1431  };
1432 };
1433 
1434 
1435 template<> struct unpacket_traits<Packet4h> { typedef Eigen::half type; enum {size=4, alignment=Aligned16, vectorizable=true, masked_load_available=false, masked_store_available=false}; typedef Packet4h half; };
1436 
1437 template<> EIGEN_STRONG_INLINE Packet4h pset1<Packet4h>(const Eigen::half& from) {
1438  Packet4h result;
1439  result.x = _mm_set1_pi16(from.x);
1440  return result;
1441 }
1442 
1443 template<> EIGEN_STRONG_INLINE Eigen::half pfirst<Packet4h>(const Packet4h& from) {
1444  return half_impl::raw_uint16_to_half(static_cast<unsigned short>(_mm_cvtsi64_si32(from.x)));
1445 }
1446 
1447 template<> EIGEN_STRONG_INLINE Packet4h pconj(const Packet4h& a) { return a; }
1448 
1449 template<> EIGEN_STRONG_INLINE Packet4h padd<Packet4h>(const Packet4h& a, const Packet4h& b) {
1450  __int64_t a64 = _mm_cvtm64_si64(a.x);
1451  __int64_t b64 = _mm_cvtm64_si64(b.x);
1452 
1453  Eigen::half h[4];
1454 
1455  Eigen::half ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64));
1456  Eigen::half hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64));
1457  h[0] = ha + hb;
1458  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 16));
1459  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 16));
1460  h[1] = ha + hb;
1461  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 32));
1462  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 32));
1463  h[2] = ha + hb;
1464  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 48));
1465  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 48));
1466  h[3] = ha + hb;
1467  Packet4h result;
1468  result.x = _mm_set_pi16(h[3].x, h[2].x, h[1].x, h[0].x);
1469  return result;
1470 }
1471 
1472 template<> EIGEN_STRONG_INLINE Packet4h psub<Packet4h>(const Packet4h& a, const Packet4h& b) {
1473  __int64_t a64 = _mm_cvtm64_si64(a.x);
1474  __int64_t b64 = _mm_cvtm64_si64(b.x);
1475 
1476  Eigen::half h[4];
1477 
1478  Eigen::half ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64));
1479  Eigen::half hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64));
1480  h[0] = ha - hb;
1481  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 16));
1482  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 16));
1483  h[1] = ha - hb;
1484  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 32));
1485  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 32));
1486  h[2] = ha - hb;
1487  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 48));
1488  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 48));
1489  h[3] = ha - hb;
1490  Packet4h result;
1491  result.x = _mm_set_pi16(h[3].x, h[2].x, h[1].x, h[0].x);
1492  return result;
1493 }
1494 
1495 template<> EIGEN_STRONG_INLINE Packet4h pmul<Packet4h>(const Packet4h& a, const Packet4h& b) {
1496  __int64_t a64 = _mm_cvtm64_si64(a.x);
1497  __int64_t b64 = _mm_cvtm64_si64(b.x);
1498 
1499  Eigen::half h[4];
1500 
1501  Eigen::half ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64));
1502  Eigen::half hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64));
1503  h[0] = ha * hb;
1504  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 16));
1505  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 16));
1506  h[1] = ha * hb;
1507  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 32));
1508  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 32));
1509  h[2] = ha * hb;
1510  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 48));
1511  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 48));
1512  h[3] = ha * hb;
1513  Packet4h result;
1514  result.x = _mm_set_pi16(h[3].x, h[2].x, h[1].x, h[0].x);
1515  return result;
1516 }
1517 
1518 template<> EIGEN_STRONG_INLINE Packet4h pdiv<Packet4h>(const Packet4h& a, const Packet4h& b) {
1519  __int64_t a64 = _mm_cvtm64_si64(a.x);
1520  __int64_t b64 = _mm_cvtm64_si64(b.x);
1521 
1522  Eigen::half h[4];
1523 
1524  Eigen::half ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64));
1525  Eigen::half hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64));
1526  h[0] = ha / hb;
1527  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 16));
1528  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 16));
1529  h[1] = ha / hb;
1530  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 32));
1531  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 32));
1532  h[2] = ha / hb;
1533  ha = half_impl::raw_uint16_to_half(static_cast<unsigned short>(a64 >> 48));
1534  hb = half_impl::raw_uint16_to_half(static_cast<unsigned short>(b64 >> 48));
1535  h[3] = ha / hb;
1536  Packet4h result;
1537  result.x = _mm_set_pi16(h[3].x, h[2].x, h[1].x, h[0].x);
1538  return result;
1539 }
1540 
1541 template<> EIGEN_STRONG_INLINE Packet4h pload<Packet4h>(const Eigen::half* from) {
1542  Packet4h result;
1543  result.x = _mm_cvtsi64_m64(*reinterpret_cast<const __int64_t*>(from));
1544  return result;
1545 }
1546 
1547 template<> EIGEN_STRONG_INLINE Packet4h ploadu<Packet4h>(const Eigen::half* from) {
1548  Packet4h result;
1549  result.x = _mm_cvtsi64_m64(*reinterpret_cast<const __int64_t*>(from));
1550  return result;
1551 }
1552 
1553 template<> EIGEN_STRONG_INLINE void pstore<Eigen::half>(Eigen::half* to, const Packet4h& from) {
1554  __int64_t r = _mm_cvtm64_si64(from.x);
1555  *(reinterpret_cast<__int64_t*>(to)) = r;
1556 }
1557 
1558 template<> EIGEN_STRONG_INLINE void pstoreu<Eigen::half>(Eigen::half* to, const Packet4h& from) {
1559  __int64_t r = _mm_cvtm64_si64(from.x);
1560  *(reinterpret_cast<__int64_t*>(to)) = r;
1561 }
1562 
1563 template<> EIGEN_STRONG_INLINE Packet4h
1564 ploadquad<Packet4h>(const Eigen::half* from) {
1565  return pset1<Packet4h>(*from);
1566 }
1567 
1568 template<> EIGEN_STRONG_INLINE Packet4h pgather<Eigen::half, Packet4h>(const Eigen::half* from, Index stride)
1569 {
1570  Packet4h result;
1571  result.x = _mm_set_pi16(from[3*stride].x, from[2*stride].x, from[1*stride].x, from[0*stride].x);
1572  return result;
1573 }
1574 
1575 template<> EIGEN_STRONG_INLINE void pscatter<Eigen::half, Packet4h>(Eigen::half* to, const Packet4h& from, Index stride)
1576 {
1577  __int64_t a = _mm_cvtm64_si64(from.x);
1578  to[stride*0].x = static_cast<unsigned short>(a);
1579  to[stride*1].x = static_cast<unsigned short>(a >> 16);
1580  to[stride*2].x = static_cast<unsigned short>(a >> 32);
1581  to[stride*3].x = static_cast<unsigned short>(a >> 48);
1582 }
1583 
1584 EIGEN_STRONG_INLINE void
1585 ptranspose(PacketBlock<Packet4h,4>& kernel) {
1586  __m64 T0 = _mm_unpacklo_pi16(kernel.packet[0].x, kernel.packet[1].x);
1587  __m64 T1 = _mm_unpacklo_pi16(kernel.packet[2].x, kernel.packet[3].x);
1588  __m64 T2 = _mm_unpackhi_pi16(kernel.packet[0].x, kernel.packet[1].x);
1589  __m64 T3 = _mm_unpackhi_pi16(kernel.packet[2].x, kernel.packet[3].x);
1590 
1591  kernel.packet[0].x = _mm_unpacklo_pi32(T0, T1);
1592  kernel.packet[1].x = _mm_unpackhi_pi32(T0, T1);
1593  kernel.packet[2].x = _mm_unpacklo_pi32(T2, T3);
1594  kernel.packet[3].x = _mm_unpackhi_pi32(T2, T3);
1595 }
1596 
1597 #endif
1598 
1599 
1600 } // end namespace internal
1601 
1602 } // end namespace Eigen
1603 
1604 #if EIGEN_COMP_PGI && EIGEN_COMP_PGI < 1900
1605 // PGI++ does not define the following intrinsics in C++ mode.
1606 static inline __m128 _mm_castpd_ps (__m128d x) { return reinterpret_cast<__m128&>(x); }
1607 static inline __m128i _mm_castpd_si128(__m128d x) { return reinterpret_cast<__m128i&>(x); }
1608 static inline __m128d _mm_castps_pd (__m128 x) { return reinterpret_cast<__m128d&>(x); }
1609 static inline __m128i _mm_castps_si128(__m128 x) { return reinterpret_cast<__m128i&>(x); }
1610 static inline __m128 _mm_castsi128_ps(__m128i x) { return reinterpret_cast<__m128&>(x); }
1611 static inline __m128d _mm_castsi128_pd(__m128i x) { return reinterpret_cast<__m128d&>(x); }
1612 #endif
1613 
1614 #endif // EIGEN_PACKET_MATH_SSE_H
Eigen::Aligned16
@ Aligned16
Definition: Constants.h:237
Eigen
Namespace containing all symbols from the Eigen library.
Definition: Core:139
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Eigen::sign
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sign_op< typename Derived::Scalar >, const Derived > sign(const Eigen::ArrayBase< Derived > &x)
Eigen::exp
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_exp_op< typename Derived::Scalar >, const Derived > exp(const Eigen::ArrayBase< Derived > &x)