First time it compiles, but fails to pass the tests.

2 files changed, 312 insertions, 71 deletions

diff --git a/Eigen/src/Core/arch/AltiVec/Complex.h b/Eigen/src/Core/arch/AltiVec/Complex.h
index 064341a3b..d90cfe0de 100644
--- a/Eigen/src/Core/arch/AltiVec/Complex.h
+++ b/Eigen/src/Core/arch/AltiVec/Complex.h
@@ -18,7 +18,7 @@ namespace internal {
 static Packet4ui  p4ui_CONJ_XOR = vec_mergeh((Packet4ui)p4i_ZERO, (Packet4ui)p4f_ZERO_);//{ 0x00000000, 0x80000000, 0x00000000, 0x80000000 };
 #else
 static Packet4ui  p4ui_CONJ_XOR = vec_mergeh((Packet4ui)p4f_ZERO_, (Packet4ui)p4i_ZERO);//{ 0x80000000, 0x00000000, 0x80000000, 0x00000000 };
-static Packet2ul  p2ul_CONJ_XOR = (Packet2ul) vec_sld((Packet4ui)p2d_ZERO_, (Packet4ui)p2l_ZERO, 8);//{ 0x8000000000000000, 0x0000000000000000 };
+static Packet2ul  p2ul_CONJ_XOR = (Packet2ul) vec_sld((Packet4ui) p2d_ZERO_, (Packet4ui)p2l_ZERO, 8);//{ 0x8000000000000000, 0x0000000000000000 };
 #endif
 
 static Packet16uc p16uc_COMPLEX_RE   = 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 };
@@ -276,28 +276,58 @@ template<> struct packet_traits<std::complex<double> >  : default_packet_traits
 
 template<> struct unpacket_traits<Packet1cd> { typedef std::complex<double> type; enum {size=1}; typedef Packet1cd half; };
 
+template<> EIGEN_STRONG_INLINE Packet1cd pload <Packet1cd>(const std::complex<double>* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet1cd(pload<Packet2d>((const double*)from)); }
+template<> EIGEN_STRONG_INLINE Packet1cd ploadu<Packet1cd>(const std::complex<double>* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet1cd(ploadu<Packet2d>((const double*)from)); }
+template<> EIGEN_STRONG_INLINE void pstore <std::complex<double> >(std::complex<double> *   to, const Packet1cd& from) { EIGEN_DEBUG_ALIGNED_STORE pstore((double*)to, from.v); }
+template<> EIGEN_STRONG_INLINE void pstoreu<std::complex<double> >(std::complex<double> *   to, const Packet1cd& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, from.v); }
+
+template<> EIGEN_STRONG_INLINE Packet1cd pset1<Packet1cd>(const std::complex<double>&  from)
+{
+  Packet1cd res;
+  /* On AltiVec we cannot load 64-bit registers, so wa have to take care of alignment */
+  if((ptrdiff_t(&from) % 16) == 0)
+    res.v = pload<Packet2d>((const double *)&from);
+  else
+    res.v = ploadu<Packet2d>((const double *)&from);
+  res.v = vec_perm(res.v, res.v, p16uc_PSET_HI);
+  return res;
+}
+template<> EIGEN_DEVICE_FUNC inline Packet1cd pgather<std::complex<double>, Packet1cd>(const std::complex<double>* from, DenseIndex stride)
+{
+  std::complex<double> EIGEN_ALIGN16 af[2];
+  af[0] = from[0*stride];
+  af[1] = from[1*stride];
+  return pload<Packet1cd>(af);
+}
+template<> EIGEN_DEVICE_FUNC inline void pscatter<std::complex<double>, Packet1cd>(std::complex<double>* to, const Packet1cd& from, DenseIndex stride)
+{
+  std::complex<double> EIGEN_ALIGN16 af[2];
+  pstore<std::complex<double> >(af, from);
+  to[0*stride] = af[0];
+  to[1*stride] = af[1];
+}
+
 template<> EIGEN_STRONG_INLINE Packet1cd padd<Packet1cd>(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_add(a.v,b.v)); }
 template<> EIGEN_STRONG_INLINE Packet1cd psub<Packet1cd>(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_sub(a.v,b.v)); }
 template<> EIGEN_STRONG_INLINE Packet1cd pnegate(const Packet1cd& a) { return Packet1cd(pnegate(Packet2d(a.v))); }
-template<> EIGEN_STRONG_INLINE Packet1cd pconj(const Packet1cd& a) { return Packet1cd((Packet2d)vec_xor((Packet2ul)a.v, p2ul_CONJ_XOR)); }
+template<> EIGEN_STRONG_INLINE Packet1cd pconj(const Packet1cd& a) { return Packet1cd((Packet2d)vec_xor((Packet2d)a.v, (Packet2d)p2ul_CONJ_XOR)); }
 
 template<> EIGEN_STRONG_INLINE Packet1cd pmul<Packet1cd>(const Packet1cd& a, const Packet1cd& b)
 {
-  Packet2d v1, v2;
+  Packet2d a_re, a_im, v1, v2;
 
   // Permute and multiply the real parts of a and b
-  v1 = vec_perm(a.v, a.v, p16uc_COMPLEX_RE);
+  a_re = vec_perm(a.v, a.v, p16uc_PSET_HI);
   // Get the imaginary parts of a
-  v2 = vec_perm(a.v, a.v, p16uc_COMPLEX_IM);
+  a_im = vec_perm(a.v, a.v, p16uc_PSET_HI);
   // multiply a_re * b
-  v1 = vec_madd(v1, b.v, p4f_ZERO);
+  v1 = vec_madd(a_re, b.v, p2d_ZERO_);
   // multiply a_im * b and get the conjugate result
-  v2 = vec_madd(v2, b.v, p4f_ZERO);
-  v2 = (Packet4f) vec_xor((Packet4ui)v2, p4ui_CONJ_XOR);
-  // permute back to a proper order
-  v2 = vec_perm(v2, v2, p16uc_COMPLEX_REV);
+  v2 = vec_madd(a_im, b.v, p2d_ZERO_);
+  v2 = (Packet2d) vec_sld((Packet4ui)v2, (Packet4ui)v2, 8);
+  v2 = (Packet2d) vec_xor((Packet2d)v2, (Packet2d) p2ul_CONJ_XOR);
 
-  return Packet2cf(vec_add(v1, v2));
+  return Packet1cd(vec_add(v1, v2));
 }
 
 template<> EIGEN_STRONG_INLINE Packet1cd pand   <Packet1cd>(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_and(a.v,b.v)); }
@@ -305,23 +335,17 @@ template<> EIGEN_STRONG_INLINE Packet1cd por    <Packet1cd>(const Packet1cd& a,
 template<> EIGEN_STRONG_INLINE Packet1cd pxor   <Packet1cd>(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_xor(a.v,b.v)); }
 template<> EIGEN_STRONG_INLINE Packet1cd pandnot<Packet1cd>(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_and(a.v, vec_nor(b.v,b.v))); }
 
-template<> EIGEN_STRONG_INLINE Packet1cd pload <Packet1cd>(const std::complex<double>* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet1cd(pload<Packet4f>((const double*)from)); }
-template<> EIGEN_STRONG_INLINE Packet1cd ploadu<Packet1cd>(const std::complex<double>* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet1cd(ploadu<Packet4f>((const double*)from)); }
-
 template<> EIGEN_STRONG_INLINE Packet1cd ploaddup<Packet1cd>(const std::complex<double>*     from)
 {
   return pset1<Packet1cd>(*from);
 }
 
-template<> EIGEN_STRONG_INLINE void pstore <std::complex<double> >(std::complex<double> *   to, const Packet1cd& from) { EIGEN_DEBUG_ALIGNED_STORE pstore((double*)to, from.v); }
-template<> EIGEN_STRONG_INLINE void pstoreu<std::complex<double> >(std::complex<double> *   to, const Packet1cd& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, from.v); }
-
-template<> EIGEN_STRONG_INLINE void prefetch<std::complex<double> >(const std::complex<double> *   addr) { vec_dstt((double *)addr, DST_CTRL(2,2,32), DST_CHAN); }
+template<> EIGEN_STRONG_INLINE void prefetch<std::complex<double> >(const std::complex<double> *   addr) { vec_dstt((long *)addr, DST_CTRL(2,2,32), DST_CHAN); }
 
 template<> EIGEN_STRONG_INLINE std::complex<double>  pfirst<Packet1cd>(const Packet1cd& a)
 {
   std::complex<double> EIGEN_ALIGN16 res[2];
-  pstore((float *)&res, a.v);
+  pstore<std::complex<double> >(res, a);
 
   return res[0];
 }
@@ -336,7 +360,7 @@ template<> EIGEN_STRONG_INLINE Packet1cd preverse(const Packet1cd& a)
 template<> EIGEN_STRONG_INLINE std::complex<double> predux<Packet1cd>(const Packet1cd& a)
 {
   Packet2d b;
-  b = (Packet2d) vec_sld(a.v, a.v, 8);
+  b = (Packet2d) vec_sld((Packet4ui) a.v, (Packet4ui) a.v, 8);
   b = padd(a.v, b);
   return pfirst(Packet1cd(b));
 }
@@ -345,9 +369,9 @@ template<> EIGEN_STRONG_INLINE Packet1cd preduxp<Packet1cd>(const Packet1cd* vec
 {
   Packet2d b1, b2;
  
-  b1 = (Packet2d) vec_sld(vecs[0].v, vecs[1].v, 8);
-  b2 = (Packet2d) vec_sld(vecs[1].v, vecs[0].v, 8);
-  b2 = (Packet2d) vec_sld(b2, b2, 8);
+  b1 = (Packet2d) vec_sld((Packet4ui) vecs[0].v, (Packet4ui) vecs[1].v, 8);
+  b2 = (Packet2d) vec_sld((Packet4ui) vecs[1].v, (Packet4ui) vecs[0].v, 8);
+  b2 = (Packet2d) vec_sld((Packet4ui) b2, (Packet4ui) b2, 8);
   b2 = padd(b1, b2);
 
   return Packet1cd(b2);
@@ -357,7 +381,7 @@ template<> EIGEN_STRONG_INLINE std::complex<double> predux_mul<Packet1cd>(const
 {
   Packet2d b;
   Packet1cd prod;
-  b = (Packet2d) vec_sld(a.v, a.v, 8);
+  b = (Packet2d) vec_sld((Packet4ui) a.v, (Packet4ui) a.v, 8);
   prod = pmul(a, Packet1cd(b));
 
   return pfirst(prod);
@@ -407,7 +431,7 @@ template<> EIGEN_STRONG_INLINE Packet1cd pdiv<Packet1cd>(const Packet1cd& a, con
 {
   // TODO optimize it for AltiVec
   Packet1cd res = conj_helper<Packet1cd,Packet1cd,false,true>().pmul(a,b);
-  Packet2d s = vec_madd(b.v, b.v, p4f_ZERO);
+  Packet2d s = vec_madd(b.v, b.v, p2d_ZERO_);
   return Packet1cd(pdiv(res.v, vec_add(s,vec_perm(s, s, p16uc_COMPLEX_REV))));
 }
 
diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h
index e70039ae2..0489ad886 100755
--- a/Eigen/src/Core/arch/AltiVec/PacketMath.h
+++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h
@@ -31,12 +31,6 @@ namespace internal {
 #define EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS 16
 #endif
 
-#ifdef __VSX__
-typedef __vector double         Packet2d;
-typedef __vector unsigned long  Packet2ul;
-typedef __vector long           Packet2l;
-#endif // __VSX__
-
 typedef __vector float          Packet4f;
 typedef __vector int            Packet4i;
 typedef __vector unsigned int   Packet4ui;
@@ -82,8 +76,6 @@ static Packet4i p4i_COUNTDOWN = { 3, 2, 1, 0 };
 static Packet16uc p16uc_REVERSE = { 0,1,2,3, 4,5,6,7, 8,9,10,11, 12,13,14,15 };
 static Packet16uc p16uc_FORWARD = { 12,13,14,15, 8,9,10,11, 4,5,6,7, 0,1,2,3}; 
 static Packet16uc p16uc_DUPLICATE = { 4,5,6,7, 4,5,6,7, 0,1,2,3, 0,1,2,3 };
-static Packet2d p2d_ZERO_ = (Packet2d) { 0x8000000000000000, 0x8000000000000000 };
-static Packet2l p2l_ZERO = (Packet2l) { 0x0, 0x0 };
 #endif // _BIG_ENDIAN
 
 // These constants are endian-agnostic
@@ -114,24 +106,6 @@ template<> struct packet_traits<float>  : default_packet_traits
     HasSqrt = 0
   };
 };
-#ifdef __VSX__
-template<> struct packet_traits<double> : default_packet_traits
-{
-  typedef Packet2d type;
-  typedef Packet2d half;
-  enum {
-    Vectorizable = 1,
-    AlignedOnScalar = 1,
-    size=2,
-    HasHalfPacket = 0,
-
-    HasDiv  = 1,
-    HasExp  = 0,
-    HasSqrt = 0
-  };
-};
-#endif // __VSX__
-
 template<> struct packet_traits<int>    : default_packet_traits
 {
   typedef Packet4i type;
@@ -144,9 +118,6 @@ template<> struct packet_traits<int>    : default_packet_traits
   };
 };
 
-#ifdef __VSX__
-template<> struct unpacket_traits<Packet2d> { typedef double type; enum {size=2}; typedef Packet2d half; };
-#endif // __VSX__
 
 template<> struct unpacket_traits<Packet4f> { typedef float  type; enum {size=4}; typedef Packet4f half; };
 template<> struct unpacket_traits<Packet4i> { typedef int    type; enum {size=4}; typedef Packet4i half; };
@@ -196,11 +167,18 @@ inline std::ostream & operator <<(std::ostream & s, const Packetbi & v)
 }
 */
 
+// Need to define them first or we get specialization after instantiation errors
+template<> EIGEN_STRONG_INLINE Packet4f pload<Packet4f>(const float* from) { EIGEN_DEBUG_ALIGNED_LOAD return vec_ld(0, from); }
+template<> EIGEN_STRONG_INLINE Packet4i pload<Packet4i>(const int*     from) { EIGEN_DEBUG_ALIGNED_LOAD return vec_ld(0, from); }
+
+template<> EIGEN_STRONG_INLINE void pstore<float>(float*   to, const Packet4f& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st(from, 0, to); }
+template<> EIGEN_STRONG_INLINE void pstore<int>(int*       to, const Packet4i& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st(from, 0, to); }
+
 template<> EIGEN_STRONG_INLINE Packet4f pset1<Packet4f>(const float&  from) {
   // Taken from http://developer.apple.com/hardwaredrivers/ve/alignment.html
   float EIGEN_ALIGN16 af[4];
   af[0] = from;
-  Packet4f vc = vec_ld(0, af);
+  Packet4f vc = pload<Packet4f>(af);
   vc = vec_splat(vc, 0);
   return vc;
 }
@@ -208,17 +186,15 @@ template<> EIGEN_STRONG_INLINE Packet4f pset1<Packet4f>(const float&  from) {
 template<> EIGEN_STRONG_INLINE Packet4i pset1<Packet4i>(const int&    from)   {
   int EIGEN_ALIGN16 ai[4];
   ai[0] = from;
-  Packet4i vc = vec_ld(0, ai);
+  Packet4i vc = pload<Packet4i>(ai);
   vc = vec_splat(vc, 0);
   return vc;
 }
-
-
 template<> EIGEN_STRONG_INLINE void
 pbroadcast4<Packet4f>(const float *a,
                       Packet4f& a0, Packet4f& a1, Packet4f& a2, Packet4f& a3)
 {
-  a3 = vec_ld(0,a);
+  a3 = pload<Packet4f>(a);
   a0 = vec_splat(a3, 0);
   a1 = vec_splat(a3, 1);
   a2 = vec_splat(a3, 2);
@@ -228,7 +204,7 @@ template<> EIGEN_STRONG_INLINE void
 pbroadcast4<Packet4i>(const int *a,
                       Packet4i& a0, Packet4i& a1, Packet4i& a2, Packet4i& a3)
 {
-  a3 = vec_ld(0,a);
+  a3 = pload<Packet4i>(a);
   a0 = vec_splat(a3, 0);
   a1 = vec_splat(a3, 1);
   a2 = vec_splat(a3, 2);
@@ -242,7 +218,7 @@ template<> EIGEN_DEVICE_FUNC inline Packet4f pgather<float, Packet4f>(const floa
   af[1] = from[1*stride];
   af[2] = from[2*stride];
   af[3] = from[3*stride];
- return vec_ld(0, af);
+ return pload<Packet4f>(af);
 }
 template<> EIGEN_DEVICE_FUNC inline Packet4i pgather<int, Packet4i>(const int* from, DenseIndex stride)
 {
@@ -251,12 +227,12 @@ template<> EIGEN_DEVICE_FUNC inline Packet4i pgather<int, Packet4i>(const int* f
   ai[1] = from[1*stride];
   ai[2] = from[2*stride];
   ai[3] = from[3*stride];
- return vec_ld(0, ai);
+ return pload<Packet4i>(ai);
 }
 template<> EIGEN_DEVICE_FUNC inline void pscatter<float, Packet4f>(float* to, const Packet4f& from, DenseIndex stride)
 {
   float EIGEN_ALIGN16 af[4];
-  vec_st(from, 0, af);
+  pstore<float>(af, from);
   to[0*stride] = af[0];
   to[1*stride] = af[1];
   to[2*stride] = af[2];
@@ -265,7 +241,7 @@ template<> EIGEN_DEVICE_FUNC inline void pscatter<float, Packet4f>(float* to, co
 template<> EIGEN_DEVICE_FUNC inline void pscatter<int, Packet4i>(int* to, const Packet4i& from, DenseIndex stride)
 {
   int EIGEN_ALIGN16 ai[4];
-  vec_st(from, 0, ai);
+  pstore<int>((int *)ai, from);
   to[0*stride] = ai[0];
   to[1*stride] = ai[1];
   to[2*stride] = ai[2];
@@ -366,10 +342,7 @@ template<> EIGEN_STRONG_INLINE Packet4i pxor<Packet4i>(const Packet4i& a, const
 template<> EIGEN_STRONG_INLINE Packet4f pandnot<Packet4f>(const Packet4f& a, const Packet4f& b) { return vec_and(a, vec_nor(b, b)); }
 template<> EIGEN_STRONG_INLINE Packet4i pandnot<Packet4i>(const Packet4i& a, const Packet4i& b) { return vec_and(a, vec_nor(b, b)); }
 
-template<> EIGEN_STRONG_INLINE Packet4f pload<Packet4f>(const float* from) { EIGEN_DEBUG_ALIGNED_LOAD return vec_ld(0, from); }
-template<> EIGEN_STRONG_INLINE Packet4i pload<Packet4i>(const int*     from) { EIGEN_DEBUG_ALIGNED_LOAD return vec_ld(0, from); }
-
-#ifndef __VSX__
+#ifdef _BIG_ENDIAN
 template<> EIGEN_STRONG_INLINE Packet4f ploadu<Packet4f>(const float* from)
 {
   EIGEN_DEBUG_ALIGNED_LOAD
@@ -393,6 +366,7 @@ template<> EIGEN_STRONG_INLINE Packet4i ploadu<Packet4i>(const int* from)
   return (Packet4i) vec_perm(MSQ, LSQ, mask);    // align the data
 }
 #else
+// We also need ot redefine little endian loading of Packet4i/Packet4f using VSX
 template<> EIGEN_STRONG_INLINE Packet4i ploadu<Packet4i>(const int* from)
 {
   EIGEN_DEBUG_ALIGNED_LOAD
@@ -420,9 +394,6 @@ template<> EIGEN_STRONG_INLINE Packet4i ploaddup<Packet4i>(const int*     from)
   return vec_perm(p, p, p16uc_DUPLICATE);
 }
 
-template<> EIGEN_STRONG_INLINE void pstore<float>(float*   to, const Packet4f& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st(from, 0, to); }
-template<> EIGEN_STRONG_INLINE void pstore<int>(int*       to, const Packet4i& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st(from, 0, to); }
-
 template<> EIGEN_STRONG_INLINE void pstoreu<float>(float*  to, const Packet4f& from)
 {
   EIGEN_DEBUG_UNALIGNED_STORE
@@ -642,6 +613,252 @@ ptranspose(PacketBlock<Packet4i,4>& kernel) {
   kernel.packet[3] = vec_mergel(t1, t3);
 }
 
+
+//---------- double ----------
+#ifdef __VSX__
+typedef __vector double              Packet2d;
+typedef __vector unsigned long long  Packet2ul;
+typedef __vector long long           Packet2l;
+
+static Packet2l p2l_ZERO = { 0x0, 0x0 };
+static Packet2d p2l_ONE = { 1, 1 };
+static Packet2d p2d_ONE = { 1.0, 1.0 }; 
+static Packet2d p2d_ZERO = { 0.0, 0.0 };
+static Packet2d p2d_ZERO_ = { (double) 0x8000000000000000, (double) 0x8000000000000000 };
+static Packet2d p2d_COUNTDOWN = { 1.0, 0.0 };
+
+#ifdef _BIG_ENDIAN
+static Packet16uc p16uc_SPLATQ1 = { 0,1,2,3, 4,5,6,7, 0,1,2,3, 4,5,6,7 }; 
+static Packet16uc p16uc_SPLATQ2 = { 8,9,10,11, 12,13,14,15, 8,9,10,11, 12,13,14,15 }; 
+#else
+static Packet16uc p16uc_SPLATQ1 = { 4,5,6,7, 0,1,2,3, 4,5,6,7, 0,1,3,4 }; 
+static Packet16uc p16uc_SPLATQ2 = { 12,13,14,15, 8,9,10,11, 12,13,14,15, 8,9,10,11 }; 
+#endif
+
+
+static EIGEN_STRONG_INLINE Packet2d vec_splat_dbl(Packet2d& a, int index)
+{
+  switch (index) {
+  case 0:
+    return (Packet2d) vec_perm(a, a, p16uc_SPLATQ1);
+  case 1:
+    return (Packet2d) vec_perm(a, a, p16uc_SPLATQ2);
+  }
+  return a;
+}
+
+template<> struct packet_traits<double> : default_packet_traits
+{
+  typedef Packet2d type;
+  typedef Packet2d half;
+  enum {
+    Vectorizable = 1,
+    AlignedOnScalar = 1,
+    size=2,
+    HasHalfPacket = 0,
+
+    HasDiv  = 1,
+    HasExp  = 0,
+    HasSqrt = 0
+  };
+};
+
+template<> struct unpacket_traits<Packet2d> { typedef double type; enum {size=2}; typedef Packet2d half; };
+
+/*
+inline std::ostream & operator <<(std::ostream & s, const Packet2d & v)
+{
+  union {
+    Packetdf   v;
+    double n[2];
+  } vt;
+  vt.v = v;
+  s << vt.n[0] << ", " << vt.n[1];
+  return s;
+}*/
+
+// Need to define them first or we get specialization after instantiation errors
+template<> EIGEN_STRONG_INLINE Packet2d pload<Packet2d>(const double* from) { EIGEN_DEBUG_ALIGNED_LOAD return (Packet2d) vec_ld(0, (const float *) from); } //FIXME
+
+template<> EIGEN_STRONG_INLINE void pstore<double>(double*   to, const Packet2d& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st((Packet4f)from, 0, (float *)to); }
+
+template<> EIGEN_STRONG_INLINE Packet2d pset1<Packet2d>(const double&  from) {
+  double EIGEN_ALIGN16 af[2];
+  af[0] = from;
+  Packet2d vc = pload<Packet2d>(af);
+  vc = vec_splat_dbl(vc, 0);
+  return vc;
+}
+template<> EIGEN_STRONG_INLINE void
+pbroadcast4<Packet2d>(const double *a,
+                      Packet2d& a0, Packet2d& a1, Packet2d& a2, Packet2d& a3)
+{
+  a1 = pload<Packet2d>(a);
+  a0 = vec_splat_dbl(a1, 0);
+  a1 = vec_splat_dbl(a1, 1);
+  a3 = pload<Packet2d>(a+2);
+  a2 = vec_splat_dbl(a3, 0);
+  a3 = vec_splat_dbl(a1, 1);
+}
+template<> EIGEN_DEVICE_FUNC inline Packet2d pgather<double, Packet2d>(const double* from, DenseIndex stride)
+{
+  double EIGEN_ALIGN16 af[2];
+  af[0] = from[0*stride];
+  af[1] = from[1*stride];
+ return pload<Packet2d>(af);
+}
+template<> EIGEN_DEVICE_FUNC inline void pscatter<double, Packet2d>(double* to, const Packet2d& from, DenseIndex stride)
+{
+  double EIGEN_ALIGN16 af[2];
+  pstore<double>(af, from);
+  to[0*stride] = af[0];
+  to[1*stride] = af[1];
+}
+template<> EIGEN_STRONG_INLINE Packet2d plset<double>(const double& a) { return vec_add(pset1<Packet2d>(a), p2d_COUNTDOWN); }
+
+template<> EIGEN_STRONG_INLINE Packet2d padd<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_add(a,b); }
+
+template<> EIGEN_STRONG_INLINE Packet2d psub<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_sub(a,b); }
+
+template<> EIGEN_STRONG_INLINE Packet2d pnegate(const Packet2d& a) { return psub<Packet2d>(p2d_ZERO, a); }
+
+template<> EIGEN_STRONG_INLINE Packet2d pconj(const Packet2d& a) { return a; }
+
+template<> EIGEN_STRONG_INLINE Packet2d pmul<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_madd(a,b,p2d_ZERO); }
+template<> EIGEN_STRONG_INLINE Packet2d pdiv<Packet2d>(const Packet2d& a, const Packet2d& b)
+{
+  Packet2d t, y_0, y_1, res;
+
+  // Altivec does not offer a divide instruction, we have to do a reciprocal approximation
+  y_0 = vec_re(b);
+
+  // Do one Newton-Raphson iteration to get the needed accuracy
+  t   = vec_nmsub(y_0, b, p2d_ONE);
+  y_1 = vec_madd(y_0, t, y_0);
+
+  res = vec_madd(a, y_1, p2d_ZERO);
+  return res;
+}
+
+// for some weird raisons, it has to be overloaded for packet of integers
+template<> EIGEN_STRONG_INLINE Packet2d pmadd(const Packet2d& a, const Packet2d& b, const Packet2d& c) { return vec_madd(a, b, c); }
+
+template<> EIGEN_STRONG_INLINE Packet2d pmin<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_min(a, b); }
+
+template<> EIGEN_STRONG_INLINE Packet2d pmax<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_max(a, b); }
+
+template<> EIGEN_STRONG_INLINE Packet2d pand<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_and(a, b); }
+
+template<> EIGEN_STRONG_INLINE Packet2d por<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_or(a, b); }
+
+template<> EIGEN_STRONG_INLINE Packet2d pxor<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_xor(a, b); }
+
+template<> EIGEN_STRONG_INLINE Packet2d pandnot<Packet2d>(const Packet2d& a, const Packet2d& b) { return vec_and(a, vec_nor(b, b)); }
+
+
+template<> EIGEN_STRONG_INLINE Packet2d ploadu<Packet2d>(const double* from)
+{
+  EIGEN_DEBUG_ALIGNED_LOAD
+  return (Packet2d) vec_vsx_ld((int)((size_t)(from) & 15), (const float *)from);      // align the data
+}
+template<> EIGEN_STRONG_INLINE Packet2d ploaddup<Packet2d>(const double*   from)
+{
+  Packet2d p;
+  if((ptrdiff_t(from) % 16) == 0)  p = pload<Packet2d>(from);
+  else                             p = ploadu<Packet2d>(from);
+  return vec_perm(p, p, p16uc_DUPLICATE);
+}
+
+template<> EIGEN_STRONG_INLINE void pstoreu<double>(double*  to, const Packet2d& from)
+{
+  EIGEN_DEBUG_UNALIGNED_STORE
+  // Taken from http://developer.apple.com/hardwaredrivers/ve/alignment.html
+  // Warning: not thread safe!
+  Packet16uc MSQ, LSQ, edges;
+  Packet16uc edgeAlign, align;
+
+  MSQ = vec_ld(0, (unsigned char *)to);                     // most significant quadword
+  LSQ = vec_ld(15, (unsigned char *)to);                    // least significant quadword
+  edgeAlign = vec_lvsl(0, to);                              // permute map to extract edges
+  edges=vec_perm(LSQ,MSQ,edgeAlign);                        // extract the edges
+  align = vec_lvsr( 0, to );                                // permute map to misalign data
+  MSQ = vec_perm(edges,(Packet16uc)from,align);             // misalign the data (MSQ)
+  LSQ = vec_perm((Packet16uc)from,edges,align);             // misalign the data (LSQ)
+  vec_st( LSQ, 15, (unsigned char *)to );                   // Store the LSQ part first
+  vec_st( MSQ, 0, (unsigned char *)to );                    // Store the MSQ part
+}
+template<> EIGEN_STRONG_INLINE void prefetch<double>(const double* addr) { vec_dstt((const float *) addr, DST_CTRL(2,2,32), DST_CHAN); }
+
+template<> EIGEN_STRONG_INLINE double  pfirst<Packet2d>(const Packet2d& a) { double EIGEN_ALIGN16 x[2]; pstore(x, a); return x[0]; }
+
+template<> EIGEN_STRONG_INLINE Packet2d preverse(const Packet2d& a) { return (Packet2d)vec_perm((Packet16uc)a,(Packet16uc)a, p16uc_REVERSE); }
+
+template<> EIGEN_STRONG_INLINE Packet2d pabs(const Packet2d& a) { return vec_abs(a); }
+
+template<> EIGEN_STRONG_INLINE double predux<Packet2d>(const Packet2d& a)
+{
+  Packet2d b, sum;
+  b   = (Packet2d) vec_sld((Packet4ui) a, (Packet4ui)a, 8);
+  sum = vec_add(a, b);
+  return pfirst(sum);
+}
+
+template<> EIGEN_STRONG_INLINE Packet2d preduxp<Packet2d>(const Packet2d* vecs)
+{
+  Packet2d v[2], sum;
+
+  v[0] = (Packet2d) vec_sld((Packet4ui) vecs[0], (Packet4ui) vecs[1], 8);
+  v[1] = (Packet2d) vec_sld((Packet4ui) vecs[1], (Packet4ui) vecs[0], 8);
+  sum = vec_add(v[0], v[1]);
+
+  return sum;
+}
+// Other reduction functions:
+// mul
+template<> EIGEN_STRONG_INLINE double predux_mul<Packet2d>(const Packet2d& a)
+{
+  Packet2d prod;
+  prod = pmul(a, (Packet2d)vec_sld((Packet4ui) a, (Packet4ui)a, 8));
+  return pfirst(pmul(prod, (Packet2d)vec_sld((Packet4ui) prod, (Packet4ui) prod, 4)));
+}
+
+// min
+template<> EIGEN_STRONG_INLINE double predux_min<Packet2d>(const Packet2d& a)
+{
+  Packet2d b, res;
+  b = vec_min(a, (Packet2d) vec_sld((Packet4ui) a, (Packet4ui) a, 8));
+  res = vec_min(b, (Packet2d) vec_sld((Packet4ui) b, (Packet4ui) b, 4));
+  return pfirst(res);
+}
+
+// max
+template<> EIGEN_STRONG_INLINE double predux_max<Packet2d>(const Packet2d& a)
+{
+  Packet2d res;
+  res = vec_max(a, (Packet2d) vec_sld((Packet4ui) a, (Packet4ui) a, 8));
+  return pfirst(res);
+}
+
+template<int Offset>
+struct palign_impl<Offset,Packet2d>
+{
+  static EIGEN_STRONG_INLINE void run(Packet2d& first, const Packet2d& second)
+  {
+    if (Offset!=0)
+      first = (Packet2d) vec_sld((Packet4ui) first, (Packet4ui) second, Offset*8);
+  }
+};
+
+EIGEN_DEVICE_FUNC inline void
+ptranspose(PacketBlock<Packet2d,2>& kernel) {
+  Packet2d t0, t1;
+  t0 = (Packet2d) vec_sld((Packet4ui) kernel.packet[0], (Packet4ui) kernel.packet[1], 8);
+  t1 = (Packet2d) vec_sld((Packet4ui) kernel.packet[1], (Packet4ui) kernel.packet[0], 8);
+  kernel.packet[0] = t0;
+  kernel.packet[1] = t1;
+}
+
+#endif // __VSX__
 } // end namespace internal
 
 } // end namespace Eigen