Home | History | Annotate | Download | only in ZVector
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2016 Konstantinos Margaritis <markos (at) freevec.org>
      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_ZVECTOR_H
     11 #define EIGEN_PACKET_MATH_ZVECTOR_H
     12 
     13 #include <stdint.h>
     14 
     15 namespace Eigen {
     16 
     17 namespace internal {
     18 
     19 #ifndef EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
     20 #define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 4
     21 #endif
     22 
     23 #ifndef EIGEN_HAS_SINGLE_INSTRUCTION_MADD
     24 #define EIGEN_HAS_SINGLE_INSTRUCTION_MADD
     25 #endif
     26 
     27 #ifndef EIGEN_HAS_SINGLE_INSTRUCTION_CJMADD
     28 #define EIGEN_HAS_SINGLE_INSTRUCTION_CJMADD
     29 #endif
     30 
     31 #ifndef EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS
     32 #define EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS  16
     33 #endif
     34 
     35 typedef __vector int                 Packet4i;
     36 typedef __vector unsigned int        Packet4ui;
     37 typedef __vector __bool int          Packet4bi;
     38 typedef __vector short int           Packet8i;
     39 typedef __vector unsigned char       Packet16uc;
     40 typedef __vector double              Packet2d;
     41 typedef __vector unsigned long long  Packet2ul;
     42 typedef __vector long long           Packet2l;
     43 
     44 typedef struct {
     45 	Packet2d  v4f[2];
     46 } Packet4f;
     47 
     48 typedef union {
     49   int32_t   i[4];
     50   uint32_t ui[4];
     51   int64_t   l[2];
     52   uint64_t ul[2];
     53   double    d[2];
     54   Packet4i  v4i;
     55   Packet4ui v4ui;
     56   Packet2l  v2l;
     57   Packet2ul v2ul;
     58   Packet2d  v2d;
     59 } Packet;
     60 
     61 // We don't want to write the same code all the time, but we need to reuse the constants
     62 // and it doesn't really work to declare them global, so we define macros instead
     63 
     64 #define _EIGEN_DECLARE_CONST_FAST_Packet4i(NAME,X) \
     65   Packet4i p4i_##NAME = reinterpret_cast<Packet4i>(vec_splat_s32(X))
     66 
     67 #define _EIGEN_DECLARE_CONST_FAST_Packet2d(NAME,X) \
     68   Packet2d p2d_##NAME = reinterpret_cast<Packet2d>(vec_splat_s64(X))
     69 
     70 #define _EIGEN_DECLARE_CONST_FAST_Packet2l(NAME,X) \
     71   Packet2l p2l_##NAME = reinterpret_cast<Packet2l>(vec_splat_s64(X))
     72 
     73 #define _EIGEN_DECLARE_CONST_Packet4i(NAME,X) \
     74   Packet4i p4i_##NAME = pset1<Packet4i>(X)
     75 
     76 #define _EIGEN_DECLARE_CONST_Packet2d(NAME,X) \
     77   Packet2d p2d_##NAME = pset1<Packet2d>(X)
     78 
     79 #define _EIGEN_DECLARE_CONST_Packet2l(NAME,X) \
     80   Packet2l p2l_##NAME = pset1<Packet2l>(X)
     81 
     82 // These constants are endian-agnostic
     83 //static _EIGEN_DECLARE_CONST_FAST_Packet4i(ZERO, 0); //{ 0, 0, 0, 0,}
     84 static _EIGEN_DECLARE_CONST_FAST_Packet4i(ONE, 1); //{ 1, 1, 1, 1}
     85 
     86 static _EIGEN_DECLARE_CONST_FAST_Packet2d(ZERO, 0);
     87 static _EIGEN_DECLARE_CONST_FAST_Packet2l(ZERO, 0);
     88 static _EIGEN_DECLARE_CONST_FAST_Packet2l(ONE, 1);
     89 
     90 static Packet2d p2d_ONE = { 1.0, 1.0 };
     91 static Packet2d p2d_ZERO_ = { -0.0, -0.0 };
     92 
     93 static Packet4i p4i_COUNTDOWN = { 0, 1, 2, 3 };
     94 static Packet4f p4f_COUNTDOWN = { 0.0, 1.0, 2.0, 3.0 };
     95 static Packet2d p2d_COUNTDOWN = reinterpret_cast<Packet2d>(vec_sld(reinterpret_cast<Packet16uc>(p2d_ZERO), reinterpret_cast<Packet16uc>(p2d_ONE), 8));
     96 
     97 static Packet16uc p16uc_PSET64_HI = { 0,1,2,3, 4,5,6,7, 0,1,2,3, 4,5,6,7 };
     98 static Packet16uc p16uc_DUPLICATE32_HI = { 0,1,2,3, 0,1,2,3, 4,5,6,7, 4,5,6,7 };
     99 
    100 // Mask alignment
    101 #define _EIGEN_MASK_ALIGNMENT	0xfffffffffffffff0
    102 
    103 #define _EIGEN_ALIGNED_PTR(x)	((std::ptrdiff_t)(x) & _EIGEN_MASK_ALIGNMENT)
    104 
    105 // Handle endianness properly while loading constants
    106 // Define global static constants:
    107 
    108 static Packet16uc p16uc_FORWARD =   { 0,1,2,3, 4,5,6,7, 8,9,10,11, 12,13,14,15 };
    109 static Packet16uc p16uc_REVERSE32 = { 12,13,14,15, 8,9,10,11, 4,5,6,7, 0,1,2,3 };
    110 static Packet16uc p16uc_REVERSE64 = { 8,9,10,11, 12,13,14,15, 0,1,2,3, 4,5,6,7 };
    111 
    112 static Packet16uc p16uc_PSET32_WODD   = vec_sld((Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 0), (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 2), 8);//{ 0,1,2,3, 0,1,2,3, 8,9,10,11, 8,9,10,11 };
    113 static Packet16uc p16uc_PSET32_WEVEN  = vec_sld(p16uc_DUPLICATE32_HI, (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 3), 8);//{ 4,5,6,7, 4,5,6,7, 12,13,14,15, 12,13,14,15 };
    114 /*static Packet16uc p16uc_HALF64_0_16 = vec_sld((Packet16uc)p4i_ZERO, vec_splat((Packet16uc) vec_abs(p4i_MINUS16), 3), 8);      //{ 0,0,0,0, 0,0,0,0, 16,16,16,16, 16,16,16,16};
    115 
    116 static Packet16uc p16uc_PSET64_HI = (Packet16uc) vec_mergeh((Packet4ui)p16uc_PSET32_WODD, (Packet4ui)p16uc_PSET32_WEVEN);     //{ 0,1,2,3, 4,5,6,7, 0,1,2,3, 4,5,6,7 };*/
    117 static Packet16uc p16uc_PSET64_LO = (Packet16uc) vec_mergel((Packet4ui)p16uc_PSET32_WODD, (Packet4ui)p16uc_PSET32_WEVEN);     //{ 8,9,10,11, 12,13,14,15, 8,9,10,11, 12,13,14,15 };
    118 /*static Packet16uc p16uc_TRANSPOSE64_HI = vec_add(p16uc_PSET64_HI, p16uc_HALF64_0_16);                                         //{ 0,1,2,3, 4,5,6,7, 16,17,18,19, 20,21,22,23};
    119 static Packet16uc p16uc_TRANSPOSE64_LO = vec_add(p16uc_PSET64_LO, p16uc_HALF64_0_16);                                         //{ 8,9,10,11, 12,13,14,15, 24,25,26,27, 28,29,30,31};*/
    120 static Packet16uc p16uc_TRANSPOSE64_HI = { 0,1,2,3, 4,5,6,7, 16,17,18,19, 20,21,22,23};
    121 static Packet16uc p16uc_TRANSPOSE64_LO = { 8,9,10,11, 12,13,14,15, 24,25,26,27, 28,29,30,31};
    122 
    123 //static Packet16uc p16uc_COMPLEX32_REV = vec_sld(p16uc_REVERSE32, p16uc_REVERSE32, 8);                                         //{ 4,5,6,7, 0,1,2,3, 12,13,14,15, 8,9,10,11 };
    124 
    125 //static Packet16uc p16uc_COMPLEX32_REV2 = vec_sld(p16uc_FORWARD, p16uc_FORWARD, 8);                                            //{ 8,9,10,11, 12,13,14,15, 0,1,2,3, 4,5,6,7 };
    126 
    127 
    128 #if EIGEN_HAS_BUILTIN(__builtin_prefetch) || EIGEN_COMP_GNUC
    129   #define EIGEN_ZVECTOR_PREFETCH(ADDR) __builtin_prefetch(ADDR);
    130 #else
    131   #define EIGEN_ZVECTOR_PREFETCH(ADDR) asm( "   pfd [%[addr]]\n" :: [addr] "r" (ADDR) : "cc" );
    132 #endif
    133 
    134 template<> struct packet_traits<int>    : default_packet_traits
    135 {
    136   typedef Packet4i type;
    137   typedef Packet4i half;
    138   enum {
    139     Vectorizable = 1,
    140     AlignedOnScalar = 1,
    141     size = 4,
    142     HasHalfPacket = 0,
    143 
    144     HasAdd  = 1,
    145     HasSub  = 1,
    146     HasMul  = 1,
    147     HasDiv  = 1,
    148     HasBlend = 1
    149   };
    150 };
    151 
    152 template<> struct packet_traits<float> : default_packet_traits
    153 {
    154   typedef Packet4f type;
    155   typedef Packet4f half;
    156   enum {
    157     Vectorizable = 1,
    158     AlignedOnScalar = 1,
    159     size=4,
    160     HasHalfPacket = 0,
    161 
    162     HasAdd  = 1,
    163     HasSub  = 1,
    164     HasMul  = 1,
    165     HasDiv  = 1,
    166     HasMin  = 1,
    167     HasMax  = 1,
    168     HasAbs  = 1,
    169     HasSin  = 0,
    170     HasCos  = 0,
    171     HasLog  = 0,
    172     HasExp  = 1,
    173     HasSqrt = 1,
    174     HasRsqrt = 1,
    175     HasRound = 1,
    176     HasFloor = 1,
    177     HasCeil = 1,
    178     HasNegate = 1,
    179     HasBlend = 1
    180   };
    181 };
    182 
    183 template<> struct packet_traits<double> : default_packet_traits
    184 {
    185   typedef Packet2d type;
    186   typedef Packet2d half;
    187   enum {
    188     Vectorizable = 1,
    189     AlignedOnScalar = 1,
    190     size=2,
    191     HasHalfPacket = 1,
    192 
    193     HasAdd  = 1,
    194     HasSub  = 1,
    195     HasMul  = 1,
    196     HasDiv  = 1,
    197     HasMin  = 1,
    198     HasMax  = 1,
    199     HasAbs  = 1,
    200     HasSin  = 0,
    201     HasCos  = 0,
    202     HasLog  = 0,
    203     HasExp  = 1,
    204     HasSqrt = 1,
    205     HasRsqrt = 1,
    206     HasRound = 1,
    207     HasFloor = 1,
    208     HasCeil = 1,
    209     HasNegate = 1,
    210     HasBlend = 1
    211   };
    212 };
    213 
    214 template<> struct unpacket_traits<Packet4i> { typedef int    type; enum {size=4, alignment=Aligned16}; typedef Packet4i half; };
    215 template<> struct unpacket_traits<Packet4f> { typedef float  type; enum {size=4, alignment=Aligned16}; typedef Packet4f half; };
    216 template<> struct unpacket_traits<Packet2d> { typedef double type; enum {size=2, alignment=Aligned16}; typedef Packet2d half; };
    217 
    218 /* Forward declaration */
    219 EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock<Packet4f,4>& kernel);
    220 
    221 inline std::ostream & operator <<(std::ostream & s, const Packet4i & v)
    222 {
    223   Packet vt;
    224   vt.v4i = v;
    225   s << vt.i[0] << ", " << vt.i[1] << ", " << vt.i[2] << ", " << vt.i[3];
    226   return s;
    227 }
    228 
    229 inline std::ostream & operator <<(std::ostream & s, const Packet4ui & v)
    230 {
    231   Packet vt;
    232   vt.v4ui = v;
    233   s << vt.ui[0] << ", " << vt.ui[1] << ", " << vt.ui[2] << ", " << vt.ui[3];
    234   return s;
    235 }
    236 
    237 inline std::ostream & operator <<(std::ostream & s, const Packet2l & v)
    238 {
    239   Packet vt;
    240   vt.v2l = v;
    241   s << vt.l[0] << ", " << vt.l[1];
    242   return s;
    243 }
    244 
    245 inline std::ostream & operator <<(std::ostream & s, const Packet2ul & v)
    246 {
    247   Packet vt;
    248   vt.v2ul = v;
    249   s << vt.ul[0] << ", " << vt.ul[1] ;
    250   return s;
    251 }
    252 
    253 inline std::ostream & operator <<(std::ostream & s, const Packet2d & v)
    254 {
    255   Packet vt;
    256   vt.v2d = v;
    257   s << vt.d[0] << ", " << vt.d[1];
    258   return s;
    259 }
    260 
    261 /* Helper function to simulate a vec_splat_packet4f
    262  */
    263 template<int element> EIGEN_STRONG_INLINE Packet4f vec_splat_packet4f(const Packet4f&   from)
    264 {
    265   Packet4f splat;
    266   switch (element) {
    267   case 0:
    268     splat.v4f[0] = vec_splat(from.v4f[0], 0);
    269     splat.v4f[1] = splat.v4f[0];
    270     break;
    271   case 1:
    272     splat.v4f[0] = vec_splat(from.v4f[0], 1);
    273     splat.v4f[1] = splat.v4f[0];
    274     break;
    275   case 2:
    276     splat.v4f[0] = vec_splat(from.v4f[1], 0);
    277     splat.v4f[1] = splat.v4f[0];
    278     break;
    279   case 3:
    280     splat.v4f[0] = vec_splat(from.v4f[1], 1);
    281     splat.v4f[1] = splat.v4f[0];
    282     break;
    283   }
    284   return splat;
    285 }
    286 
    287 template<int Offset>
    288 struct palign_impl<Offset,Packet4i>
    289 {
    290   static EIGEN_STRONG_INLINE void run(Packet4i& first, const Packet4i& second)
    291   {
    292     switch (Offset % 4) {
    293     case 1:
    294       first = vec_sld(first, second, 4); break;
    295     case 2:
    296       first = vec_sld(first, second, 8); break;
    297     case 3:
    298       first = vec_sld(first, second, 12); break;
    299     }
    300   }
    301 };
    302 
    303 /* This is a tricky one, we have to translate float alignment to vector elements of sizeof double
    304  */
    305 template<int Offset>
    306 struct palign_impl<Offset,Packet4f>
    307 {
    308   static EIGEN_STRONG_INLINE void run(Packet4f& first, const Packet4f& second)
    309   {
    310     switch (Offset % 4) {
    311     case 1:
    312       first.v4f[0] = vec_sld(first.v4f[0], first.v4f[1], 8);
    313       first.v4f[1] = vec_sld(first.v4f[1], second.v4f[0], 8);
    314       break;
    315     case 2:
    316       first.v4f[0] = first.v4f[1];
    317       first.v4f[1] = second.v4f[0];
    318       break;
    319     case 3:
    320       first.v4f[0] = vec_sld(first.v4f[1],  second.v4f[0], 8);
    321       first.v4f[1] = vec_sld(second.v4f[0], second.v4f[1], 8);
    322       break;
    323     }
    324   }
    325 };
    326 
    327 
    328 template<int Offset>
    329 struct palign_impl<Offset,Packet2d>
    330 {
    331   static EIGEN_STRONG_INLINE void run(Packet2d& first, const Packet2d& second)
    332   {
    333     if (Offset == 1)
    334       first = reinterpret_cast<Packet2d>(vec_sld(reinterpret_cast<Packet4i>(first), reinterpret_cast<Packet4i>(second), 8));
    335   }
    336 };
    337 
    338 template<> EIGEN_STRONG_INLINE Packet4i pload<Packet4i>(const int*     from)
    339 {
    340   // FIXME: No intrinsic yet
    341   EIGEN_DEBUG_ALIGNED_LOAD
    342   Packet *vfrom;
    343   vfrom = (Packet *) from;
    344   return vfrom->v4i;
    345 }
    346 
    347 template<> EIGEN_STRONG_INLINE Packet4f pload<Packet4f>(const float*   from)
    348 {
    349   // FIXME: No intrinsic yet
    350   EIGEN_DEBUG_ALIGNED_LOAD
    351   Packet4f vfrom;
    352   vfrom.v4f[0] = vec_ld2f(&from[0]);
    353   vfrom.v4f[1] = vec_ld2f(&from[2]);
    354   return vfrom;
    355 }
    356 
    357 template<> EIGEN_STRONG_INLINE Packet2d pload<Packet2d>(const double* from)
    358 {
    359   // FIXME: No intrinsic yet
    360   EIGEN_DEBUG_ALIGNED_LOAD
    361   Packet *vfrom;
    362   vfrom = (Packet *) from;
    363   return vfrom->v2d;
    364 }
    365 
    366 template<> EIGEN_STRONG_INLINE void pstore<int>(int*       to, const Packet4i& from)
    367 {
    368   // FIXME: No intrinsic yet
    369   EIGEN_DEBUG_ALIGNED_STORE
    370   Packet *vto;
    371   vto = (Packet *) to;
    372   vto->v4i = from;
    373 }
    374 
    375 template<> EIGEN_STRONG_INLINE void pstore<float>(float*   to, const Packet4f& from)
    376 {
    377   // FIXME: No intrinsic yet
    378   EIGEN_DEBUG_ALIGNED_STORE
    379   vec_st2f(from.v4f[0], &to[0]);
    380   vec_st2f(from.v4f[1], &to[2]);
    381 }
    382 
    383 
    384 template<> EIGEN_STRONG_INLINE void pstore<double>(double*   to, const Packet2d& from)
    385 {
    386   // FIXME: No intrinsic yet
    387   EIGEN_DEBUG_ALIGNED_STORE
    388   Packet *vto;
    389   vto = (Packet *) to;
    390   vto->v2d = from;
    391 }
    392 
    393 template<> EIGEN_STRONG_INLINE Packet4i pset1<Packet4i>(const int&    from)
    394 {
    395   return vec_splats(from);
    396 }
    397 template<> EIGEN_STRONG_INLINE Packet2d pset1<Packet2d>(const double& from) {
    398   return vec_splats(from);
    399 }
    400 template<> EIGEN_STRONG_INLINE Packet4f pset1<Packet4f>(const float&    from)
    401 {
    402   Packet4f to;
    403   to.v4f[0] = pset1<Packet2d>(static_cast<const double&>(from));
    404   to.v4f[1] = to.v4f[0];
    405   return to;
    406 }
    407 
    408 template<> EIGEN_STRONG_INLINE void
    409 pbroadcast4<Packet4i>(const int *a,
    410                       Packet4i& a0, Packet4i& a1, Packet4i& a2, Packet4i& a3)
    411 {
    412   a3 = pload<Packet4i>(a);
    413   a0 = vec_splat(a3, 0);
    414   a1 = vec_splat(a3, 1);
    415   a2 = vec_splat(a3, 2);
    416   a3 = vec_splat(a3, 3);
    417 }
    418 
    419 template<> EIGEN_STRONG_INLINE void
    420 pbroadcast4<Packet4f>(const float *a,
    421                       Packet4f& a0, Packet4f& a1, Packet4f& a2, Packet4f& a3)
    422 {
    423   a3 = pload<Packet4f>(a);
    424   a0 = vec_splat_packet4f<0>(a3);
    425   a1 = vec_splat_packet4f<1>(a3);
    426   a2 = vec_splat_packet4f<2>(a3);
    427   a3 = vec_splat_packet4f<3>(a3);
    428 }
    429 
    430 template<> EIGEN_STRONG_INLINE void
    431 pbroadcast4<Packet2d>(const double *a,
    432                       Packet2d& a0, Packet2d& a1, Packet2d& a2, Packet2d& a3)
    433 {
    434   a1 = pload<Packet2d>(a);
    435   a0 = vec_splat(a1, 0);
    436   a1 = vec_splat(a1, 1);
    437   a3 = pload<Packet2d>(a+2);
    438   a2 = vec_splat(a3, 0);
    439   a3 = vec_splat(a3, 1);
    440 }
    441 
    442 template<> EIGEN_DEVICE_FUNC inline Packet4i pgather<int, Packet4i>(const int* from, Index stride)
    443 {
    444   int EIGEN_ALIGN16 ai[4];
    445   ai[0] = from[0*stride];
    446   ai[1] = from[1*stride];
    447   ai[2] = from[2*stride];
    448   ai[3] = from[3*stride];
    449  return pload<Packet4i>(ai);
    450 }
    451 
    452 template<> EIGEN_DEVICE_FUNC inline Packet4f pgather<float, Packet4f>(const float* from, Index stride)
    453 {
    454   float EIGEN_ALIGN16 ai[4];
    455   ai[0] = from[0*stride];
    456   ai[1] = from[1*stride];
    457   ai[2] = from[2*stride];
    458   ai[3] = from[3*stride];
    459  return pload<Packet4f>(ai);
    460 }
    461 
    462 template<> EIGEN_DEVICE_FUNC inline Packet2d pgather<double, Packet2d>(const double* from, Index stride)
    463 {
    464   double EIGEN_ALIGN16 af[2];
    465   af[0] = from[0*stride];
    466   af[1] = from[1*stride];
    467  return pload<Packet2d>(af);
    468 }
    469 
    470 template<> EIGEN_DEVICE_FUNC inline void pscatter<int, Packet4i>(int* to, const Packet4i& from, Index stride)
    471 {
    472   int EIGEN_ALIGN16 ai[4];
    473   pstore<int>((int *)ai, from);
    474   to[0*stride] = ai[0];
    475   to[1*stride] = ai[1];
    476   to[2*stride] = ai[2];
    477   to[3*stride] = ai[3];
    478 }
    479 
    480 template<> EIGEN_DEVICE_FUNC inline void pscatter<float, Packet4f>(float* to, const Packet4f& from, Index stride)
    481 {
    482   float EIGEN_ALIGN16 ai[4];
    483   pstore<float>((float *)ai, from);
    484   to[0*stride] = ai[0];
    485   to[1*stride] = ai[1];
    486   to[2*stride] = ai[2];
    487   to[3*stride] = ai[3];
    488 }
    489 
    490 template<> EIGEN_DEVICE_FUNC inline void pscatter<double, Packet2d>(double* to, const Packet2d& from, Index stride)
    491 {
    492   double EIGEN_ALIGN16 af[2];
    493   pstore<double>(af, from);
    494   to[0*stride] = af[0];
    495   to[1*stride] = af[1];
    496 }
    497 
    498 template<> EIGEN_STRONG_INLINE Packet4i padd<Packet4i>(const Packet4i& a, const Packet4i& b) { return (a + b); }
    499 template<> EIGEN_STRONG_INLINE Packet4f padd<Packet4f>(const Packet4f& a, const Packet4f& b)
    500 {
    501   Packet4f c;
    502   c.v4f[0] = a.v4f[0] + b.v4f[0];
    503   c.v4f[1] = a.v4f[1] + b.v4f[1];
    504   return c;
    505 }
    506 template<> EIGEN_STRONG_INLINE Packet2d padd<Packet2d>(const Packet2d& a, const Packet2d& b) { return (a + b); }
    507 
    508 template<> EIGEN_STRONG_INLINE Packet4i psub<Packet4i>(const Packet4i& a, const Packet4i& b) { return (a - b); }
    509 template<> EIGEN_STRONG_INLINE Packet4f psub<Packet4f>(const Packet4f& a, const Packet4f& b)
    510 {
    511   Packet4f c;
    512   c.v4f[0] = a.v4f[0] - b.v4f[0];
    513   c.v4f[1] = a.v4f[1] - b.v4f[1];
    514   return c;
    515 }
    516 template<> EIGEN_STRONG_INLINE Packet2d psub<Packet2d>(const Packet2d& a, const Packet2d& b) { return (a - b); }
    517 
    518 template<> EIGEN_STRONG_INLINE Packet4i pmul<Packet4i>(const Packet4i& a, const Packet4i& b) { return (a * b); }
    519 template<> EIGEN_STRONG_INLINE Packet4f pmul<Packet4f>(const Packet4f& a, const Packet4f& b)
    520 {
    521   Packet4f c;
    522   c.v4f[0] = a.v4f[0] * b.v4f[0];
    523   c.v4f[1] = a.v4f[1] * b.v4f[1];
    524   return c;
    525 }
    526 template<> EIGEN_STRONG_INLINE Packet2d pmul<Packet2d>(const Packet2d& a, const Packet2d& b) { return (a * b); }
    527 
    528 template<> EIGEN_STRONG_INLINE Packet4i pdiv<Packet4i>(const Packet4i& a, const Packet4i& b) { return (a / b); }
    529 template<> EIGEN_STRONG_INLINE Packet4f pdiv<Packet4f>(const Packet4f& a, const Packet4f& b)
    530 {
    531   Packet4f c;
    532   c.v4f[0] = a.v4f[0] / b.v4f[0];
    533   c.v4f[1] = a.v4f[1] / b.v4f[1];
    534   return c;
    535 }
    536 template<> EIGEN_STRONG_INLINE Packet2d pdiv<Packet2d>(const Packet2d& a, const Packet2d& b) { return (a / b); }
    537 
    538 template<> EIGEN_STRONG_INLINE Packet4i pnegate(const Packet4i& a) { return (-a); }
    539 template<> EIGEN_STRONG_INLINE Packet4f pnegate(const Packet4f& a)
    540 {
    541   Packet4f c;
    542   c.v4f[0] = -a.v4f[0];
    543   c.v4f[1] = -a.v4f[1];
    544   return c;
    545 }
    546 template<> EIGEN_STRONG_INLINE Packet2d pnegate(const Packet2d& a) { return (-a); }
    547 
    548 template<> EIGEN_STRONG_INLINE Packet4i pconj(const Packet4i& a) { return a; }
    549 template<> EIGEN_STRONG_INLINE Packet4f pconj(const Packet4f& a) { return a; }
    550 template<> EIGEN_STRONG_INLINE Packet2d pconj(const Packet2d& a) { return a; }
    551 
    552 template<> EIGEN_STRONG_INLINE Packet4i pmadd(const Packet4i& a, const Packet4i& b, const Packet4i& c) { return padd<Packet4i>(pmul<Packet4i>(a, b), c); }
    553 template<> EIGEN_STRONG_INLINE Packet4f pmadd(const Packet4f& a, const Packet4f& b, const Packet4f& c)
    554 {
    555   Packet4f res;
    556   res.v4f[0] = vec_madd(a.v4f[0], b.v4f[0], c.v4f[0]);
    557   res.v4f[1] = vec_madd(a.v4f[1], b.v4f[1], c.v4f[1]);
    558   return res;
    559 }
    560 template<> EIGEN_STRONG_INLINE Packet2d pmadd(const Packet2d& a, const Packet2d& b, const Packet2d& c) { return vec_madd(a, b, c); }
    561 
    562 template<> EIGEN_STRONG_INLINE Packet4i plset<Packet4i>(const int& a)    { return padd<Packet4i>(pset1<Packet4i>(a), p4i_COUNTDOWN); }
    563 template<> EIGEN_STRONG_INLINE Packet4f plset<Packet4f>(const float& a)  { return padd<Packet4f>(pset1<Packet4f>(a), p4f_COUNTDOWN); }
    564 template<> EIGEN_STRONG_INLINE Packet2d plset<Packet2d>(const double& a) { return padd<Packet2d>(pset1<Packet2d>(a), p2d_COUNTDOWN); }
    565 
    566 template<> EIGEN_STRONG_INLINE Packet4i pmin<Packet4i>(const Packet4i& a, const Packet4i& b) { return vec_min(a, b); }
    567 template<> EIGEN_STRONG_INLINE Packet2d pmin<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_min(a, b); }
    568 template<> EIGEN_STRONG_INLINE Packet4f pmin<Packet4f>(const Packet4f& a, const Packet4f& b)
    569 {
    570   Packet4f res;
    571   res.v4f[0] = pmin(a.v4f[0], b.v4f[0]);
    572   res.v4f[1] = pmin(a.v4f[1], b.v4f[1]);
    573   return res;
    574 }
    575 
    576 template<> EIGEN_STRONG_INLINE Packet4i pmax<Packet4i>(const Packet4i& a, const Packet4i& b) { return vec_max(a, b); }
    577 template<> EIGEN_STRONG_INLINE Packet2d pmax<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_max(a, b); }
    578 template<> EIGEN_STRONG_INLINE Packet4f pmax<Packet4f>(const Packet4f& a, const Packet4f& b)
    579 {
    580   Packet4f res;
    581   res.v4f[0] = pmax(a.v4f[0], b.v4f[0]);
    582   res.v4f[1] = pmax(a.v4f[1], b.v4f[1]);
    583   return res;
    584 }
    585 
    586 template<> EIGEN_STRONG_INLINE Packet4i pand<Packet4i>(const Packet4i& a, const Packet4i& b) { return vec_and(a, b); }
    587 template<> EIGEN_STRONG_INLINE Packet2d pand<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_and(a, b); }
    588 template<> EIGEN_STRONG_INLINE Packet4f pand<Packet4f>(const Packet4f& a, const Packet4f& b)
    589 {
    590   Packet4f res;
    591   res.v4f[0] = pand(a.v4f[0], b.v4f[0]);
    592   res.v4f[1] = pand(a.v4f[1], b.v4f[1]);
    593   return res;
    594 }
    595 
    596 template<> EIGEN_STRONG_INLINE Packet4i por<Packet4i>(const Packet4i& a, const Packet4i& b) { return vec_or(a, b); }
    597 template<> EIGEN_STRONG_INLINE Packet2d por<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_or(a, b); }
    598 template<> EIGEN_STRONG_INLINE Packet4f por<Packet4f>(const Packet4f& a, const Packet4f& b)
    599 {
    600   Packet4f res;
    601   res.v4f[0] = pand(a.v4f[0], b.v4f[0]);
    602   res.v4f[1] = pand(a.v4f[1], b.v4f[1]);
    603   return res;
    604 }
    605 
    606 template<> EIGEN_STRONG_INLINE Packet4i pxor<Packet4i>(const Packet4i& a, const Packet4i& b) { return vec_xor(a, b); }
    607 template<> EIGEN_STRONG_INLINE Packet2d pxor<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_xor(a, b); }
    608 template<> EIGEN_STRONG_INLINE Packet4f pxor<Packet4f>(const Packet4f& a, const Packet4f& b)
    609 {
    610   Packet4f res;
    611   res.v4f[0] = pand(a.v4f[0], b.v4f[0]);
    612   res.v4f[1] = pand(a.v4f[1], b.v4f[1]);
    613   return res;
    614 }
    615 
    616 template<> EIGEN_STRONG_INLINE Packet4i pandnot<Packet4i>(const Packet4i& a, const Packet4i& b) { return pand<Packet4i>(a, vec_nor(b, b)); }
    617 template<> EIGEN_STRONG_INLINE Packet2d pandnot<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_and(a, vec_nor(b, b)); }
    618 template<> EIGEN_STRONG_INLINE Packet4f pandnot<Packet4f>(const Packet4f& a, const Packet4f& b)
    619 {
    620   Packet4f res;
    621   res.v4f[0] = pandnot(a.v4f[0], b.v4f[0]);
    622   res.v4f[1] = pandnot(a.v4f[1], b.v4f[1]);
    623   return res;
    624 }
    625 
    626 template<> EIGEN_STRONG_INLINE Packet4f pround<Packet4f>(const Packet4f& a)
    627 {
    628   Packet4f res;
    629   res.v4f[0] = vec_round(a.v4f[0]);
    630   res.v4f[1] = vec_round(a.v4f[1]);
    631   return res;
    632 }
    633 template<> EIGEN_STRONG_INLINE Packet2d pround<Packet2d>(const Packet2d& a) { return vec_round(a); }
    634 template<> EIGEN_STRONG_INLINE Packet4f pceil<Packet4f>(const  Packet4f& a)
    635 {
    636   Packet4f res;
    637   res.v4f[0] = vec_ceil(a.v4f[0]);
    638   res.v4f[1] = vec_ceil(a.v4f[1]);
    639   return res;
    640 }
    641 template<> EIGEN_STRONG_INLINE Packet2d pceil<Packet2d>(const  Packet2d& a) { return vec_ceil(a); }
    642 template<> EIGEN_STRONG_INLINE Packet4f pfloor<Packet4f>(const Packet4f& a)
    643 {
    644   Packet4f res;
    645   res.v4f[0] = vec_floor(a.v4f[0]);
    646   res.v4f[1] = vec_floor(a.v4f[1]);
    647   return res;
    648 }
    649 template<> EIGEN_STRONG_INLINE Packet2d pfloor<Packet2d>(const Packet2d& a) { return vec_floor(a); }
    650 
    651 template<> EIGEN_STRONG_INLINE Packet4i ploadu<Packet4i>(const int*       from) { return pload<Packet4i>(from); }
    652 template<> EIGEN_STRONG_INLINE Packet4f ploadu<Packet4f>(const float*     from) { return pload<Packet4f>(from); }
    653 template<> EIGEN_STRONG_INLINE Packet2d ploadu<Packet2d>(const double*    from) { return pload<Packet2d>(from); }
    654 
    655 
    656 template<> EIGEN_STRONG_INLINE Packet4i ploaddup<Packet4i>(const int*     from)
    657 {
    658   Packet4i p = pload<Packet4i>(from);
    659   return vec_perm(p, p, p16uc_DUPLICATE32_HI);
    660 }
    661 
    662 template<> EIGEN_STRONG_INLINE Packet4f ploaddup<Packet4f>(const float*    from)
    663 {
    664   Packet4f p = pload<Packet4f>(from);
    665   p.v4f[1] = vec_splat(p.v4f[0], 1);
    666   p.v4f[0] = vec_splat(p.v4f[0], 0);
    667   return p;
    668 }
    669 
    670 template<> EIGEN_STRONG_INLINE Packet2d ploaddup<Packet2d>(const double*   from)
    671 {
    672   Packet2d p = pload<Packet2d>(from);
    673   return vec_perm(p, p, p16uc_PSET64_HI);
    674 }
    675 
    676 template<> EIGEN_STRONG_INLINE void pstoreu<int>(int*        to, const Packet4i& from) { pstore<int>(to, from); }
    677 template<> EIGEN_STRONG_INLINE void pstoreu<float>(float*    to, const Packet4f& from) { pstore<float>(to, from); }
    678 template<> EIGEN_STRONG_INLINE void pstoreu<double>(double*  to, const Packet2d& from) { pstore<double>(to, from); }
    679 
    680 template<> EIGEN_STRONG_INLINE void prefetch<int>(const int*       addr) { EIGEN_ZVECTOR_PREFETCH(addr); }
    681 template<> EIGEN_STRONG_INLINE void prefetch<float>(const float*   addr) { EIGEN_ZVECTOR_PREFETCH(addr); }
    682 template<> EIGEN_STRONG_INLINE void prefetch<double>(const double* addr) { EIGEN_ZVECTOR_PREFETCH(addr); }
    683 
    684 template<> EIGEN_STRONG_INLINE int    pfirst<Packet4i>(const Packet4i& a) { int    EIGEN_ALIGN16 x[4]; pstore(x, a); return x[0]; }
    685 template<> EIGEN_STRONG_INLINE float  pfirst<Packet4f>(const Packet4f& a) { float  EIGEN_ALIGN16 x[2]; vec_st2f(a.v4f[0], &x[0]); return x[0]; }
    686 template<> EIGEN_STRONG_INLINE double pfirst<Packet2d>(const Packet2d& a) { double EIGEN_ALIGN16 x[2]; pstore(x, a); return x[0]; }
    687 
    688 template<> EIGEN_STRONG_INLINE Packet4i preverse(const Packet4i& a)
    689 {
    690   return reinterpret_cast<Packet4i>(vec_perm(reinterpret_cast<Packet16uc>(a), reinterpret_cast<Packet16uc>(a), p16uc_REVERSE32));
    691 }
    692 
    693 template<> EIGEN_STRONG_INLINE Packet2d preverse(const Packet2d& a)
    694 {
    695   return reinterpret_cast<Packet2d>(vec_perm(reinterpret_cast<Packet16uc>(a), reinterpret_cast<Packet16uc>(a), p16uc_REVERSE64));
    696 }
    697 
    698 template<> EIGEN_STRONG_INLINE Packet4f preverse(const Packet4f& a)
    699 {
    700   Packet4f rev;
    701   rev.v4f[0] = preverse<Packet2d>(a.v4f[1]);
    702   rev.v4f[1] = preverse<Packet2d>(a.v4f[0]);
    703   return rev;
    704 }
    705 
    706 template<> EIGEN_STRONG_INLINE Packet4i pabs<Packet4i>(const Packet4i& a) { return vec_abs(a); }
    707 template<> EIGEN_STRONG_INLINE Packet2d pabs<Packet2d>(const Packet2d& a) { return vec_abs(a); }
    708 template<> EIGEN_STRONG_INLINE Packet4f pabs<Packet4f>(const Packet4f& a)
    709 {
    710   Packet4f res;
    711   res.v4f[0] = pabs(a.v4f[0]);
    712   res.v4f[1] = pabs(a.v4f[1]);
    713   return res;
    714 }
    715 
    716 template<> EIGEN_STRONG_INLINE int predux<Packet4i>(const Packet4i& a)
    717 {
    718   Packet4i b, sum;
    719   b   = vec_sld(a, a, 8);
    720   sum = padd<Packet4i>(a, b);
    721   b   = vec_sld(sum, sum, 4);
    722   sum = padd<Packet4i>(sum, b);
    723   return pfirst(sum);
    724 }
    725 
    726 template<> EIGEN_STRONG_INLINE double predux<Packet2d>(const Packet2d& a)
    727 {
    728   Packet2d b, sum;
    729   b   = reinterpret_cast<Packet2d>(vec_sld(reinterpret_cast<Packet4i>(a), reinterpret_cast<Packet4i>(a), 8));
    730   sum = padd<Packet2d>(a, b);
    731   return pfirst(sum);
    732 }
    733 template<> EIGEN_STRONG_INLINE float predux<Packet4f>(const Packet4f& a)
    734 {
    735   Packet2d sum;
    736   sum = padd<Packet2d>(a.v4f[0], a.v4f[1]);
    737   double first = predux<Packet2d>(sum);
    738   return static_cast<float>(first);
    739 }
    740 
    741 template<> EIGEN_STRONG_INLINE Packet4i preduxp<Packet4i>(const Packet4i* vecs)
    742 {
    743   Packet4i v[4], sum[4];
    744 
    745   // It's easier and faster to transpose then add as columns
    746   // Check: http://www.freevec.org/function/matrix_4x4_transpose_floats for explanation
    747   // Do the transpose, first set of moves
    748   v[0] = vec_mergeh(vecs[0], vecs[2]);
    749   v[1] = vec_mergel(vecs[0], vecs[2]);
    750   v[2] = vec_mergeh(vecs[1], vecs[3]);
    751   v[3] = vec_mergel(vecs[1], vecs[3]);
    752   // Get the resulting vectors
    753   sum[0] = vec_mergeh(v[0], v[2]);
    754   sum[1] = vec_mergel(v[0], v[2]);
    755   sum[2] = vec_mergeh(v[1], v[3]);
    756   sum[3] = vec_mergel(v[1], v[3]);
    757 
    758   // Now do the summation:
    759   // Lines 0+1
    760   sum[0] = padd<Packet4i>(sum[0], sum[1]);
    761   // Lines 2+3
    762   sum[1] = padd<Packet4i>(sum[2], sum[3]);
    763   // Add the results
    764   sum[0] = padd<Packet4i>(sum[0], sum[1]);
    765 
    766   return sum[0];
    767 }
    768 
    769 template<> EIGEN_STRONG_INLINE Packet2d preduxp<Packet2d>(const Packet2d* vecs)
    770 {
    771   Packet2d v[2], sum;
    772   v[0] = padd<Packet2d>(vecs[0], reinterpret_cast<Packet2d>(vec_sld(reinterpret_cast<Packet4ui>(vecs[0]), reinterpret_cast<Packet4ui>(vecs[0]), 8)));
    773   v[1] = padd<Packet2d>(vecs[1], reinterpret_cast<Packet2d>(vec_sld(reinterpret_cast<Packet4ui>(vecs[1]), reinterpret_cast<Packet4ui>(vecs[1]), 8)));
    774 
    775   sum = reinterpret_cast<Packet2d>(vec_sld(reinterpret_cast<Packet4ui>(v[0]), reinterpret_cast<Packet4ui>(v[1]), 8));
    776 
    777   return sum;
    778 }
    779 
    780 template<> EIGEN_STRONG_INLINE Packet4f preduxp<Packet4f>(const Packet4f* vecs)
    781 {
    782   PacketBlock<Packet4f,4> transpose;
    783   transpose.packet[0] = vecs[0];
    784   transpose.packet[1] = vecs[1];
    785   transpose.packet[2] = vecs[2];
    786   transpose.packet[3] = vecs[3];
    787   ptranspose(transpose);
    788 
    789   Packet4f sum = padd(transpose.packet[0], transpose.packet[1]);
    790   sum = padd(sum, transpose.packet[2]);
    791   sum = padd(sum, transpose.packet[3]);
    792   return sum;
    793 }
    794 
    795 // Other reduction functions:
    796 // mul
    797 template<> EIGEN_STRONG_INLINE int predux_mul<Packet4i>(const Packet4i& a)
    798 {
    799   EIGEN_ALIGN16 int aux[4];
    800   pstore(aux, a);
    801   return aux[0] * aux[1] * aux[2] * aux[3];
    802 }
    803 
    804 template<> EIGEN_STRONG_INLINE double predux_mul<Packet2d>(const Packet2d& a)
    805 {
    806   return pfirst(pmul(a, reinterpret_cast<Packet2d>(vec_sld(reinterpret_cast<Packet4i>(a), reinterpret_cast<Packet4i>(a), 8))));
    807 }
    808 
    809 template<> EIGEN_STRONG_INLINE float predux_mul<Packet4f>(const Packet4f& a)
    810 {
    811   // Return predux_mul<Packet2d> of the subvectors product
    812   return static_cast<float>(pfirst(predux_mul(pmul(a.v4f[0], a.v4f[1]))));
    813 }
    814 
    815 // min
    816 template<> EIGEN_STRONG_INLINE int predux_min<Packet4i>(const Packet4i& a)
    817 {
    818   Packet4i b, res;
    819   b   = pmin<Packet4i>(a, vec_sld(a, a, 8));
    820   res = pmin<Packet4i>(b, vec_sld(b, b, 4));
    821   return pfirst(res);
    822 }
    823 
    824 template<> EIGEN_STRONG_INLINE double predux_min<Packet2d>(const Packet2d& a)
    825 {
    826   return pfirst(pmin<Packet2d>(a, reinterpret_cast<Packet2d>(vec_sld(reinterpret_cast<Packet4i>(a), reinterpret_cast<Packet4i>(a), 8))));
    827 }
    828 
    829 template<> EIGEN_STRONG_INLINE float predux_min<Packet4f>(const Packet4f& a)
    830 {
    831   Packet2d b, res;
    832   b   = pmin<Packet2d>(a.v4f[0], a.v4f[1]);
    833   res = pmin<Packet2d>(b, reinterpret_cast<Packet2d>(vec_sld(reinterpret_cast<Packet4i>(b), reinterpret_cast<Packet4i>(b), 8)));
    834   return static_cast<float>(pfirst(res));
    835 }
    836 
    837 // max
    838 template<> EIGEN_STRONG_INLINE int predux_max<Packet4i>(const Packet4i& a)
    839 {
    840   Packet4i b, res;
    841   b = pmax<Packet4i>(a, vec_sld(a, a, 8));
    842   res = pmax<Packet4i>(b, vec_sld(b, b, 4));
    843   return pfirst(res);
    844 }
    845 
    846 // max
    847 template<> EIGEN_STRONG_INLINE double predux_max<Packet2d>(const Packet2d& a)
    848 {
    849   return pfirst(pmax<Packet2d>(a, reinterpret_cast<Packet2d>(vec_sld(reinterpret_cast<Packet4i>(a), reinterpret_cast<Packet4i>(a), 8))));
    850 }
    851 
    852 template<> EIGEN_STRONG_INLINE float predux_max<Packet4f>(const Packet4f& a)
    853 {
    854   Packet2d b, res;
    855   b   = pmax<Packet2d>(a.v4f[0], a.v4f[1]);
    856   res = pmax<Packet2d>(b, reinterpret_cast<Packet2d>(vec_sld(reinterpret_cast<Packet4i>(b), reinterpret_cast<Packet4i>(b), 8)));
    857   return static_cast<float>(pfirst(res));
    858 }
    859 
    860 EIGEN_DEVICE_FUNC inline void
    861 ptranspose(PacketBlock<Packet4i,4>& kernel) {
    862   Packet4i t0 = vec_mergeh(kernel.packet[0], kernel.packet[2]);
    863   Packet4i t1 = vec_mergel(kernel.packet[0], kernel.packet[2]);
    864   Packet4i t2 = vec_mergeh(kernel.packet[1], kernel.packet[3]);
    865   Packet4i t3 = vec_mergel(kernel.packet[1], kernel.packet[3]);
    866   kernel.packet[0] = vec_mergeh(t0, t2);
    867   kernel.packet[1] = vec_mergel(t0, t2);
    868   kernel.packet[2] = vec_mergeh(t1, t3);
    869   kernel.packet[3] = vec_mergel(t1, t3);
    870 }
    871 
    872 EIGEN_DEVICE_FUNC inline void
    873 ptranspose(PacketBlock<Packet2d,2>& kernel) {
    874   Packet2d t0 = vec_perm(kernel.packet[0], kernel.packet[1], p16uc_TRANSPOSE64_HI);
    875   Packet2d t1 = vec_perm(kernel.packet[0], kernel.packet[1], p16uc_TRANSPOSE64_LO);
    876   kernel.packet[0] = t0;
    877   kernel.packet[1] = t1;
    878 }
    879 
    880 /* Split the Packet4f PacketBlock into 4 Packet2d PacketBlocks and transpose each one
    881  */
    882 EIGEN_DEVICE_FUNC inline void
    883 ptranspose(PacketBlock<Packet4f,4>& kernel) {
    884   PacketBlock<Packet2d,2> t0,t1,t2,t3;
    885   // copy top-left 2x2 Packet2d block
    886   t0.packet[0] = kernel.packet[0].v4f[0];
    887   t0.packet[1] = kernel.packet[1].v4f[0];
    888 
    889   // copy top-right 2x2 Packet2d block
    890   t1.packet[0] = kernel.packet[0].v4f[1];
    891   t1.packet[1] = kernel.packet[1].v4f[1];
    892 
    893   // copy bottom-left 2x2 Packet2d block
    894   t2.packet[0] = kernel.packet[2].v4f[0];
    895   t2.packet[1] = kernel.packet[3].v4f[0];
    896 
    897   // copy bottom-right 2x2 Packet2d block
    898   t3.packet[0] = kernel.packet[2].v4f[1];
    899   t3.packet[1] = kernel.packet[3].v4f[1];
    900 
    901   // Transpose all 2x2 blocks
    902   ptranspose(t0);
    903   ptranspose(t1);
    904   ptranspose(t2);
    905   ptranspose(t3);
    906 
    907   // Copy back transposed blocks, but exchange t1 and t2 due to transposition
    908   kernel.packet[0].v4f[0] = t0.packet[0];
    909   kernel.packet[0].v4f[1] = t2.packet[0];
    910   kernel.packet[1].v4f[0] = t0.packet[1];
    911   kernel.packet[1].v4f[1] = t2.packet[1];
    912   kernel.packet[2].v4f[0] = t1.packet[0];
    913   kernel.packet[2].v4f[1] = t3.packet[0];
    914   kernel.packet[3].v4f[0] = t1.packet[1];
    915   kernel.packet[3].v4f[1] = t3.packet[1];
    916 }
    917 
    918 template<> EIGEN_STRONG_INLINE Packet4i pblend(const Selector<4>& ifPacket, const Packet4i& thenPacket, const Packet4i& elsePacket) {
    919   Packet4ui select = { ifPacket.select[0], ifPacket.select[1], ifPacket.select[2], ifPacket.select[3] };
    920   Packet4ui mask = vec_cmpeq(select, reinterpret_cast<Packet4ui>(p4i_ONE));
    921   return vec_sel(elsePacket, thenPacket, mask);
    922 }
    923 
    924 template<> EIGEN_STRONG_INLINE Packet4f pblend(const Selector<4>& ifPacket, const Packet4f& thenPacket, const Packet4f& elsePacket) {
    925   Packet2ul select_hi = { ifPacket.select[0], ifPacket.select[1] };
    926   Packet2ul select_lo = { ifPacket.select[2], ifPacket.select[3] };
    927   Packet2ul mask_hi = vec_cmpeq(select_hi, reinterpret_cast<Packet2ul>(p2l_ONE));
    928   Packet2ul mask_lo = vec_cmpeq(select_lo, reinterpret_cast<Packet2ul>(p2l_ONE));
    929   Packet4f result;
    930   result.v4f[0] = vec_sel(elsePacket.v4f[0], thenPacket.v4f[0], mask_hi);
    931   result.v4f[1] = vec_sel(elsePacket.v4f[1], thenPacket.v4f[1], mask_lo);
    932   return result;
    933 }
    934 
    935 template<> EIGEN_STRONG_INLINE Packet2d pblend(const Selector<2>& ifPacket, const Packet2d& thenPacket, const Packet2d& elsePacket) {
    936   Packet2ul select = { ifPacket.select[0], ifPacket.select[1] };
    937   Packet2ul mask = vec_cmpeq(select, reinterpret_cast<Packet2ul>(p2l_ONE));
    938   return vec_sel(elsePacket, thenPacket, mask);
    939 }
    940 
    941 } // end namespace internal
    942 
    943 } // end namespace Eigen
    944 
    945 #endif // EIGEN_PACKET_MATH_ZVECTOR_H
    946