From 7ff266e3ce592ec1a6284cf16811965eec775f25 Mon Sep 17 00:00:00 2001 From: Konstantinos Margaritis Date: Fri, 29 Aug 2014 20:03:49 +0000 Subject: [PATCH 1/8] Initial VSX commit --- CMakeLists.txt | 6 + Eigen/Core | 17 +- Eigen/src/Core/arch/AltiVec/Complex.h | 186 ++++++++++++++++++ Eigen/src/Core/arch/AltiVec/PacketMath.h | 64 +++++- .../Core/products/GeneralBlockPanelKernel.h | 4 +- Eigen/src/Core/util/Constants.h | 6 + bench/btl/libs/eigen2/eigen2_interface.hh | 2 +- cmake/EigenTesting.cmake | 6 + test/main.h | 2 +- test/packetmath.cpp | 2 +- 10 files changed, 284 insertions(+), 11 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 96d6c8701..ea42cc8db 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -219,6 +219,12 @@ if(NOT MSVC) message(STATUS "Enabling AltiVec in tests/examples") endif() + option(EIGEN_TEST_VSX "Enable/Disable VSX in tests/examples" OFF) + if(EIGEN_TEST_VSX) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -m64 -mvsx") + message(STATUS "Enabling VSX in tests/examples") + endif() + option(EIGEN_TEST_NEON "Enable/Disable Neon in tests/examples" OFF) if(EIGEN_TEST_NEON) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mfpu=neon -mcpu=cortex-a8") diff --git a/Eigen/Core b/Eigen/Core index 776b7faf3..8ea165d5b 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -171,6 +171,15 @@ #undef bool #undef vector #undef pixel + #elif defined __VSX__ + #define EIGEN_VECTORIZE + #define EIGEN_VECTORIZE_VSX + #include + // We need to #undef all these ugly tokens defined in + // => use __vector instead of vector + #undef bool + #undef vector + #undef pixel #elif defined __ARM_NEON__ #define EIGEN_VECTORIZE #define EIGEN_VECTORIZE_NEON @@ -235,6 +244,8 @@ inline static const char *SimdInstructionSetsInUse(void) { return "SSE, SSE2"; #elif defined(EIGEN_VECTORIZE_ALTIVEC) return "AltiVec"; +#elif defined(EIGEN_VECTORIZE_VSX) + return "VSX"; #elif defined(EIGEN_VECTORIZE_NEON) return "ARM NEON"; #else @@ -286,9 +297,13 @@ using std::ptrdiff_t; #include "src/Core/arch/SSE/PacketMath.h" #include "src/Core/arch/SSE/MathFunctions.h" #include "src/Core/arch/SSE/Complex.h" -#elif defined EIGEN_VECTORIZE_ALTIVEC +#elif defined(EIGEN_VECTORIZE_ALTIVEC) #include "src/Core/arch/AltiVec/PacketMath.h" #include "src/Core/arch/AltiVec/Complex.h" +#elif defined(EIGEN_VECTORIZE_ALTIVEC) + #include "src/Core/arch/AltiVec/PacketMath.h" + #include "src/Core/arch/AltiVec/VSX.h" + #include "src/Core/arch/AltiVec/Complex.h" #elif defined EIGEN_VECTORIZE_NEON #include "src/Core/arch/NEON/PacketMath.h" #include "src/Core/arch/NEON/Complex.h" diff --git a/Eigen/src/Core/arch/AltiVec/Complex.h b/Eigen/src/Core/arch/AltiVec/Complex.h index 13b874d0c..064341a3b 100644 --- a/Eigen/src/Core/arch/AltiVec/Complex.h +++ b/Eigen/src/Core/arch/AltiVec/Complex.h @@ -14,7 +14,13 @@ namespace Eigen { namespace internal { +#ifdef _BIG_ENDIAN 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 }; +#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 }; static Packet16uc p16uc_COMPLEX_IM = vec_sld(p16uc_DUPLICATE, (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 3), 8);//{ 4,5,6,7, 4,5,6,7, 12,13,14,15, 12,13,14,15 }; static Packet16uc p16uc_COMPLEX_REV = vec_sld(p16uc_REVERSE, p16uc_REVERSE, 8);//{ 4,5,6,7, 0,1,2,3, 12,13,14,15, 8,9,10,11 }; @@ -237,6 +243,186 @@ EIGEN_STRONG_INLINE void ptranspose(PacketBlock& kernel) kernel.packet[0].v = tmp; } +//---------- double ---------- +struct Packet1cd +{ + EIGEN_STRONG_INLINE Packet1cd() {} + EIGEN_STRONG_INLINE explicit Packet1cd(const Packet2d& a) : v(a) {} + Packet2d v; +}; + +template<> struct packet_traits > : default_packet_traits +{ + typedef Packet1cd type; + typedef Packet1cd half; + enum { + Vectorizable = 1, + AlignedOnScalar = 0, + size = 1, + HasHalfPacket = 0, + + HasAdd = 1, + HasSub = 1, + HasMul = 1, + HasDiv = 1, + HasNegate = 1, + HasAbs = 0, + HasAbs2 = 0, + HasMin = 0, + HasMax = 0, + HasSetLinear = 0 + }; +}; + +template<> struct unpacket_traits { typedef std::complex type; enum {size=1}; typedef Packet1cd half; }; + +template<> EIGEN_STRONG_INLINE Packet1cd padd(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_add(a.v,b.v)); } +template<> EIGEN_STRONG_INLINE Packet1cd psub(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 pmul(const Packet1cd& a, const Packet1cd& b) +{ + Packet2d v1, v2; + + // Permute and multiply the real parts of a and b + v1 = vec_perm(a.v, a.v, p16uc_COMPLEX_RE); + // Get the imaginary parts of a + v2 = vec_perm(a.v, a.v, p16uc_COMPLEX_IM); + // multiply a_re * b + v1 = vec_madd(v1, b.v, p4f_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); + + return Packet2cf(vec_add(v1, v2)); +} + +template<> EIGEN_STRONG_INLINE Packet1cd pand (const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_and(a.v,b.v)); } +template<> EIGEN_STRONG_INLINE Packet1cd por (const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_or(a.v,b.v)); } +template<> EIGEN_STRONG_INLINE Packet1cd pxor (const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_xor(a.v,b.v)); } +template<> EIGEN_STRONG_INLINE Packet1cd pandnot(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_and(a.v, vec_nor(b.v,b.v))); } + +template<> EIGEN_STRONG_INLINE Packet1cd pload (const std::complex* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet1cd(pload((const double*)from)); } +template<> EIGEN_STRONG_INLINE Packet1cd ploadu(const std::complex* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet1cd(ploadu((const double*)from)); } + +template<> EIGEN_STRONG_INLINE Packet1cd ploaddup(const std::complex* from) +{ + return pset1(*from); +} + +template<> EIGEN_STRONG_INLINE void pstore >(std::complex * to, const Packet1cd& from) { EIGEN_DEBUG_ALIGNED_STORE pstore((double*)to, from.v); } +template<> EIGEN_STRONG_INLINE void pstoreu >(std::complex * to, const Packet1cd& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, from.v); } + +template<> EIGEN_STRONG_INLINE void prefetch >(const std::complex * addr) { vec_dstt((double *)addr, DST_CTRL(2,2,32), DST_CHAN); } + +template<> EIGEN_STRONG_INLINE std::complex pfirst(const Packet1cd& a) +{ + std::complex EIGEN_ALIGN16 res[2]; + pstore((float *)&res, a.v); + + return res[0]; +} + +template<> EIGEN_STRONG_INLINE Packet1cd preverse(const Packet1cd& a) +{ + Packet2d rev_a; + rev_a = vec_perm(a.v, a.v, p16uc_COMPLEX_REV2); + return Packet1cd(rev_a); +} + +template<> EIGEN_STRONG_INLINE std::complex predux(const Packet1cd& a) +{ + Packet2d b; + b = (Packet2d) vec_sld(a.v, a.v, 8); + b = padd(a.v, b); + return pfirst(Packet1cd(b)); +} + +template<> EIGEN_STRONG_INLINE Packet1cd preduxp(const Packet1cd* vecs) +{ + 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); + b2 = padd(b1, b2); + + return Packet1cd(b2); +} + +template<> EIGEN_STRONG_INLINE std::complex predux_mul(const Packet1cd& a) +{ + Packet2d b; + Packet1cd prod; + b = (Packet2d) vec_sld(a.v, a.v, 8); + prod = pmul(a, Packet1cd(b)); + + return pfirst(prod); +} + +template +struct palign_impl +{ + static EIGEN_STRONG_INLINE void run(Packet1cd& first, const Packet1cd& second) + { } +}; + +template<> struct conj_helper +{ + EIGEN_STRONG_INLINE Packet1cd pmadd(const Packet1cd& x, const Packet1cd& y, const Packet1cd& c) const + { return padd(pmul(x,y),c); } + + EIGEN_STRONG_INLINE Packet1cd pmul(const Packet1cd& a, const Packet1cd& b) const + { + return internal::pmul(a, pconj(b)); + } +}; + +template<> struct conj_helper +{ + EIGEN_STRONG_INLINE Packet1cd pmadd(const Packet1cd& x, const Packet1cd& y, const Packet1cd& c) const + { return padd(pmul(x,y),c); } + + EIGEN_STRONG_INLINE Packet1cd pmul(const Packet1cd& a, const Packet1cd& b) const + { + return internal::pmul(pconj(a), b); + } +}; + +template<> struct conj_helper +{ + EIGEN_STRONG_INLINE Packet1cd pmadd(const Packet1cd& x, const Packet1cd& y, const Packet1cd& c) const + { return padd(pmul(x,y),c); } + + EIGEN_STRONG_INLINE Packet1cd pmul(const Packet1cd& a, const Packet1cd& b) const + { + return pconj(internal::pmul(a, b)); + } +}; + +template<> EIGEN_STRONG_INLINE Packet1cd pdiv(const Packet1cd& a, const Packet1cd& b) +{ + // TODO optimize it for AltiVec + Packet1cd res = conj_helper().pmul(a,b); + Packet2d s = vec_madd(b.v, b.v, p4f_ZERO); + return Packet1cd(pdiv(res.v, vec_add(s,vec_perm(s, s, p16uc_COMPLEX_REV)))); +} + +template<> EIGEN_STRONG_INLINE Packet1cd pcplxflip(const Packet1cd& x) +{ + return Packet1cd(vec_perm(x.v, x.v, p16uc_COMPLEX_REV)); +} + +EIGEN_STRONG_INLINE void ptranspose(PacketBlock& kernel) +{ + Packet2d tmp = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_COMPLEX_TRANSPOSE_0); + kernel.packet[1].v = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_COMPLEX_TRANSPOSE_1); + kernel.packet[0].v = tmp; +} + } // end namespace internal } // end namespace Eigen diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h index b43e8ace3..e70039ae2 100755 --- a/Eigen/src/Core/arch/AltiVec/PacketMath.h +++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h @@ -31,6 +31,12 @@ 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; @@ -50,22 +56,37 @@ typedef __vector unsigned char Packet16uc; #define _EIGEN_DECLARE_CONST_Packet4f(NAME,X) \ Packet4f p4f_##NAME = pset1(X) -#define _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(NAME,X) \ - Packet4f p4f_##NAME = vreinterpretq_f32_u32(pset1(X)) - #define _EIGEN_DECLARE_CONST_Packet4i(NAME,X) \ Packet4i p4i_##NAME = pset1(X) +#define _EIGEN_DECLARE_CONST_Packet2d(NAME,X) \ + Packet2d p2d_##NAME = pset1(X) + +#define _EIGEN_DECLARE_CONST_Packet2l(NAME,X) \ + Packet2l p2l_##NAME = pset1(X) + #define DST_CHAN 1 #define DST_CTRL(size, count, stride) (((size) << 24) | ((count) << 16) | (stride)) +// Handle endianness properly while loading constants // Define global static constants: +#ifdef _BIG_ENDIAN static Packet4f p4f_COUNTDOWN = { 0.0, 1.0, 2.0, 3.0 }; static Packet4i p4i_COUNTDOWN = { 0, 1, 2, 3 }; static Packet16uc p16uc_REVERSE = { 12,13,14,15, 8,9,10,11, 4,5,6,7, 0,1,2,3}; static Packet16uc p16uc_FORWARD = vec_lvsl(0, (float*)0); //{ 0,1,2,3, 4,5,6,7, 8,9,10,11, 12,13,14,15} static Packet16uc p16uc_DUPLICATE = { 0,1,2,3, 0,1,2,3, 4,5,6,7, 4,5,6,7}; +#else +static Packet4f p4f_COUNTDOWN = { 3.0, 2.0, 1.0, 0.0 }; +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 static _EIGEN_DECLARE_CONST_FAST_Packet4f(ZERO, 0); //{ 0.0, 0.0, 0.0, 0.0} static _EIGEN_DECLARE_CONST_FAST_Packet4i(ZERO, 0); //{ 0, 0, 0, 0,} static _EIGEN_DECLARE_CONST_FAST_Packet4i(ONE,1); //{ 1, 1, 1, 1} @@ -93,6 +114,24 @@ template<> struct packet_traits : default_packet_traits HasSqrt = 0 }; }; +#ifdef __VSX__ +template<> struct packet_traits : 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 : default_packet_traits { typedef Packet4i type; @@ -105,6 +144,10 @@ template<> struct packet_traits : default_packet_traits }; }; +#ifdef __VSX__ +template<> struct unpacket_traits { typedef double type; enum {size=2}; typedef Packet2d half; }; +#endif // __VSX__ + template<> struct unpacket_traits { typedef float type; enum {size=4}; typedef Packet4f half; }; template<> struct unpacket_traits { typedef int type; enum {size=4}; typedef Packet4i half; }; /* @@ -311,7 +354,6 @@ template<> EIGEN_STRONG_INLINE Packet4i pmin(const Packet4i& a, const template<> EIGEN_STRONG_INLINE Packet4f pmax(const Packet4f& a, const Packet4f& b) { return vec_max(a, b); } template<> EIGEN_STRONG_INLINE Packet4i pmax(const Packet4i& a, const Packet4i& b) { return vec_max(a, b); } -// Logical Operations are not supported for float, so we have to reinterpret casts using NEON intrinsics template<> EIGEN_STRONG_INLINE Packet4f pand(const Packet4f& a, const Packet4f& b) { return vec_and(a, b); } template<> EIGEN_STRONG_INLINE Packet4i pand(const Packet4i& a, const Packet4i& b) { return vec_and(a, b); } @@ -327,10 +369,10 @@ template<> EIGEN_STRONG_INLINE Packet4i pandnot(const Packet4i& a, con template<> EIGEN_STRONG_INLINE Packet4f pload(const float* from) { EIGEN_DEBUG_ALIGNED_LOAD return vec_ld(0, from); } template<> EIGEN_STRONG_INLINE Packet4i pload(const int* from) { EIGEN_DEBUG_ALIGNED_LOAD return vec_ld(0, from); } +#ifndef __VSX__ template<> EIGEN_STRONG_INLINE Packet4f ploadu(const float* from) { EIGEN_DEBUG_ALIGNED_LOAD - // Taken from http://developer.apple.com/hardwaredrivers/ve/alignment.html Packet16uc MSQ, LSQ; Packet16uc mask; MSQ = vec_ld(0, (unsigned char *)from); // most significant quadword @@ -350,6 +392,18 @@ template<> EIGEN_STRONG_INLINE Packet4i ploadu(const int* from) mask = vec_lvsl(0, from); // create the permute mask return (Packet4i) vec_perm(MSQ, LSQ, mask); // align the data } +#else +template<> EIGEN_STRONG_INLINE Packet4i ploadu(const int* from) +{ + EIGEN_DEBUG_ALIGNED_LOAD + return (Packet4i) vec_vsx_ld((long)from & 15, from); // align the data +} +template<> EIGEN_STRONG_INLINE Packet4f ploadu(const float* from) +{ + EIGEN_DEBUG_ALIGNED_LOAD + return (Packet4f) vec_vsx_ld((long)from & 15, from); // align the data +} +#endif template<> EIGEN_STRONG_INLINE Packet4f ploaddup(const float* from) { diff --git a/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/Eigen/src/Core/products/GeneralBlockPanelKernel.h index 7da52c2e8..7b2ed6728 100644 --- a/Eigen/src/Core/products/GeneralBlockPanelKernel.h +++ b/Eigen/src/Core/products/GeneralBlockPanelKernel.h @@ -182,7 +182,7 @@ public: nr = 4, // register block size along the M direction (currently, this one cannot be modified) -#if defined(EIGEN_HAS_FUSED_MADD) && !defined(EIGEN_VECTORIZE_ALTIVEC) +#if defined(EIGEN_HAS_FUSED_MADD) && !defined(EIGEN_VECTORIZE_ALTIVEC) && !defined(EIGEN_VECTORIZE_VSX) // we assume 16 registers mr = 3*LhsPacketSize, #else @@ -290,7 +290,7 @@ public: NumberOfRegisters = EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS, nr = 4, -#if defined(EIGEN_HAS_FUSED_MADD) && !defined(EIGEN_VECTORIZE_ALTIVEC) +#if defined(EIGEN_HAS_FUSED_MADD) && !defined(EIGEN_VECTORIZE_ALTIVEC) && !defined(EIGEN_VECTORIZE_VSX) // we assume 16 registers mr = 3*LhsPacketSize, #else diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 05107fdfe..31073b990 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -413,10 +413,16 @@ namespace Architecture Generic = 0x0, SSE = 0x1, AltiVec = 0x2, + VSX = 0x3, + NEON = 0x4, #if defined EIGEN_VECTORIZE_SSE Target = SSE #elif defined EIGEN_VECTORIZE_ALTIVEC Target = AltiVec +#elif defined EIGEN_VECTORIZE_VSX + Target = VSX +#elif defined EIGEN_VECTORIZE_NEON + Target = NEON #else Target = Generic #endif diff --git a/bench/btl/libs/eigen2/eigen2_interface.hh b/bench/btl/libs/eigen2/eigen2_interface.hh index 47fe58135..1deabdae2 100644 --- a/bench/btl/libs/eigen2/eigen2_interface.hh +++ b/bench/btl/libs/eigen2/eigen2_interface.hh @@ -47,7 +47,7 @@ public : { #if defined(EIGEN_VECTORIZE_SSE) if (SIZE==Dynamic) return "eigen2"; else return "tiny_eigen2"; - #elif defined(EIGEN_VECTORIZE_ALTIVEC) + #elif defined(EIGEN_VECTORIZE_ALTIVEC) || defined(EIGEN_VECTORIZE_VSX) if (SIZE==Dynamic) return "eigen2"; else return "tiny_eigen2"; #else if (SIZE==Dynamic) return "eigen2_novec"; else return "tiny_eigen2_novec"; diff --git a/cmake/EigenTesting.cmake b/cmake/EigenTesting.cmake index 9b9776894..6383c6eac 100644 --- a/cmake/EigenTesting.cmake +++ b/cmake/EigenTesting.cmake @@ -282,6 +282,12 @@ macro(ei_testing_print_summary) message(STATUS "Altivec: Using architecture defaults") endif() + if(EIGEN_TEST_VSX) + message(STATUS "VSX: ON") + else() + message(STATUS "VSX: Using architecture defaults") + endif() + if(EIGEN_TEST_NEON) message(STATUS "ARM NEON: ON") else() diff --git a/test/main.h b/test/main.h index 773873a0d..7667eaa18 100644 --- a/test/main.h +++ b/test/main.h @@ -76,7 +76,7 @@ #endif // bounds integer values for AltiVec -#ifdef __ALTIVEC__ +#if defined(__ALTIVEC__) || defined(__VSX__) #define EIGEN_MAKING_DOCS #endif diff --git a/test/packetmath.cpp b/test/packetmath.cpp index e5dc473c2..e716d6d9a 100644 --- a/test/packetmath.cpp +++ b/test/packetmath.cpp @@ -156,7 +156,7 @@ template void packetmath() CHECK_CWISE2(REF_ADD, internal::padd); CHECK_CWISE2(REF_SUB, internal::psub); CHECK_CWISE2(REF_MUL, internal::pmul); - #ifndef EIGEN_VECTORIZE_ALTIVEC + #if !defined(EIGEN_VECTORIZE_ALTIVEC) && !defined(EIGEN_VECTORIZE_VSX) if (!internal::is_same::value) CHECK_CWISE2(REF_DIV, internal::pdiv); #endif From 470aa15c35539a7b79dd053c74578ac9e769f6e5 Mon Sep 17 00:00:00 2001 From: Konstantinos Margaritis Date: Tue, 9 Sep 2014 16:58:48 +0000 Subject: [PATCH 2/8] First time it compiles, but fails to pass the tests. --- Eigen/Core | 6 +- Eigen/src/Core/arch/AltiVec/Complex.h | 74 ++++-- Eigen/src/Core/arch/AltiVec/PacketMath.h | 309 +++++++++++++++++++---- 3 files changed, 313 insertions(+), 76 deletions(-) diff --git a/Eigen/Core b/Eigen/Core index 8ea165d5b..ac3cbd0c7 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -297,13 +297,9 @@ using std::ptrdiff_t; #include "src/Core/arch/SSE/PacketMath.h" #include "src/Core/arch/SSE/MathFunctions.h" #include "src/Core/arch/SSE/Complex.h" -#elif defined(EIGEN_VECTORIZE_ALTIVEC) +#elif defined(EIGEN_VECTORIZE_ALTIVEC) || defined(EIGEN_VECTORIZE_VSX) #include "src/Core/arch/AltiVec/PacketMath.h" #include "src/Core/arch/AltiVec/Complex.h" -#elif defined(EIGEN_VECTORIZE_ALTIVEC) - #include "src/Core/arch/AltiVec/PacketMath.h" - #include "src/Core/arch/AltiVec/VSX.h" - #include "src/Core/arch/AltiVec/Complex.h" #elif defined EIGEN_VECTORIZE_NEON #include "src/Core/arch/NEON/PacketMath.h" #include "src/Core/arch/NEON/Complex.h" 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 > : default_packet_traits template<> struct unpacket_traits { typedef std::complex type; enum {size=1}; typedef Packet1cd half; }; +template<> EIGEN_STRONG_INLINE Packet1cd pload (const std::complex* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet1cd(pload((const double*)from)); } +template<> EIGEN_STRONG_INLINE Packet1cd ploadu(const std::complex* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet1cd(ploadu((const double*)from)); } +template<> EIGEN_STRONG_INLINE void pstore >(std::complex * to, const Packet1cd& from) { EIGEN_DEBUG_ALIGNED_STORE pstore((double*)to, from.v); } +template<> EIGEN_STRONG_INLINE void pstoreu >(std::complex * to, const Packet1cd& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, from.v); } + +template<> EIGEN_STRONG_INLINE Packet1cd pset1(const std::complex& 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((const double *)&from); + else + res.v = ploadu((const double *)&from); + res.v = vec_perm(res.v, res.v, p16uc_PSET_HI); + return res; +} +template<> EIGEN_DEVICE_FUNC inline Packet1cd pgather, Packet1cd>(const std::complex* from, DenseIndex stride) +{ + std::complex EIGEN_ALIGN16 af[2]; + af[0] = from[0*stride]; + af[1] = from[1*stride]; + return pload(af); +} +template<> EIGEN_DEVICE_FUNC inline void pscatter, Packet1cd>(std::complex* to, const Packet1cd& from, DenseIndex stride) +{ + std::complex EIGEN_ALIGN16 af[2]; + pstore >(af, from); + to[0*stride] = af[0]; + to[1*stride] = af[1]; +} + template<> EIGEN_STRONG_INLINE Packet1cd padd(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_add(a.v,b.v)); } template<> EIGEN_STRONG_INLINE Packet1cd psub(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(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 (const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_and(a.v,b.v)); } @@ -305,23 +335,17 @@ template<> EIGEN_STRONG_INLINE Packet1cd por (const Packet1cd& a, template<> EIGEN_STRONG_INLINE Packet1cd pxor (const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_xor(a.v,b.v)); } template<> EIGEN_STRONG_INLINE Packet1cd pandnot(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_and(a.v, vec_nor(b.v,b.v))); } -template<> EIGEN_STRONG_INLINE Packet1cd pload (const std::complex* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet1cd(pload((const double*)from)); } -template<> EIGEN_STRONG_INLINE Packet1cd ploadu(const std::complex* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet1cd(ploadu((const double*)from)); } - template<> EIGEN_STRONG_INLINE Packet1cd ploaddup(const std::complex* from) { return pset1(*from); } -template<> EIGEN_STRONG_INLINE void pstore >(std::complex * to, const Packet1cd& from) { EIGEN_DEBUG_ALIGNED_STORE pstore((double*)to, from.v); } -template<> EIGEN_STRONG_INLINE void pstoreu >(std::complex * to, const Packet1cd& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, from.v); } - -template<> EIGEN_STRONG_INLINE void prefetch >(const std::complex * addr) { vec_dstt((double *)addr, DST_CTRL(2,2,32), DST_CHAN); } +template<> EIGEN_STRONG_INLINE void prefetch >(const std::complex * addr) { vec_dstt((long *)addr, DST_CTRL(2,2,32), DST_CHAN); } template<> EIGEN_STRONG_INLINE std::complex pfirst(const Packet1cd& a) { std::complex EIGEN_ALIGN16 res[2]; - pstore((float *)&res, a.v); + pstore >(res, a); return res[0]; } @@ -336,7 +360,7 @@ template<> EIGEN_STRONG_INLINE Packet1cd preverse(const Packet1cd& a) template<> EIGEN_STRONG_INLINE std::complex predux(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(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 predux_mul(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(const Packet1cd& a, con { // TODO optimize it for AltiVec Packet1cd res = conj_helper().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 : default_packet_traits HasSqrt = 0 }; }; -#ifdef __VSX__ -template<> struct packet_traits : 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 : default_packet_traits { typedef Packet4i type; @@ -144,9 +118,6 @@ template<> struct packet_traits : default_packet_traits }; }; -#ifdef __VSX__ -template<> struct unpacket_traits { typedef double type; enum {size=2}; typedef Packet2d half; }; -#endif // __VSX__ template<> struct unpacket_traits { typedef float type; enum {size=4}; typedef Packet4f half; }; template<> struct unpacket_traits { 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(const float* from) { EIGEN_DEBUG_ALIGNED_LOAD return vec_ld(0, from); } +template<> EIGEN_STRONG_INLINE Packet4i pload(const int* from) { EIGEN_DEBUG_ALIGNED_LOAD return vec_ld(0, from); } + +template<> EIGEN_STRONG_INLINE void pstore(float* to, const Packet4f& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st(from, 0, to); } +template<> EIGEN_STRONG_INLINE void pstore(int* to, const Packet4i& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st(from, 0, to); } + template<> EIGEN_STRONG_INLINE Packet4f pset1(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(af); vc = vec_splat(vc, 0); return vc; } @@ -208,17 +186,15 @@ template<> EIGEN_STRONG_INLINE Packet4f pset1(const float& from) { template<> EIGEN_STRONG_INLINE Packet4i pset1(const int& from) { int EIGEN_ALIGN16 ai[4]; ai[0] = from; - Packet4i vc = vec_ld(0, ai); + Packet4i vc = pload(ai); vc = vec_splat(vc, 0); return vc; } - - template<> EIGEN_STRONG_INLINE void pbroadcast4(const float *a, Packet4f& a0, Packet4f& a1, Packet4f& a2, Packet4f& a3) { - a3 = vec_ld(0,a); + a3 = pload(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(const int *a, Packet4i& a0, Packet4i& a1, Packet4i& a2, Packet4i& a3) { - a3 = vec_ld(0,a); + a3 = pload(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(const floa af[1] = from[1*stride]; af[2] = from[2*stride]; af[3] = from[3*stride]; - return vec_ld(0, af); + return pload(af); } template<> EIGEN_DEVICE_FUNC inline Packet4i pgather(const int* from, DenseIndex stride) { @@ -251,12 +227,12 @@ template<> EIGEN_DEVICE_FUNC inline Packet4i pgather(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(ai); } template<> EIGEN_DEVICE_FUNC inline void pscatter(float* to, const Packet4f& from, DenseIndex stride) { float EIGEN_ALIGN16 af[4]; - vec_st(from, 0, af); + pstore(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* to, co template<> EIGEN_DEVICE_FUNC inline void pscatter(int* to, const Packet4i& from, DenseIndex stride) { int EIGEN_ALIGN16 ai[4]; - vec_st(from, 0, ai); + pstore((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(const Packet4i& a, const template<> EIGEN_STRONG_INLINE Packet4f pandnot(const Packet4f& a, const Packet4f& b) { return vec_and(a, vec_nor(b, b)); } template<> EIGEN_STRONG_INLINE Packet4i pandnot(const Packet4i& a, const Packet4i& b) { return vec_and(a, vec_nor(b, b)); } -template<> EIGEN_STRONG_INLINE Packet4f pload(const float* from) { EIGEN_DEBUG_ALIGNED_LOAD return vec_ld(0, from); } -template<> EIGEN_STRONG_INLINE Packet4i pload(const int* from) { EIGEN_DEBUG_ALIGNED_LOAD return vec_ld(0, from); } - -#ifndef __VSX__ +#ifdef _BIG_ENDIAN template<> EIGEN_STRONG_INLINE Packet4f ploadu(const float* from) { EIGEN_DEBUG_ALIGNED_LOAD @@ -393,6 +366,7 @@ template<> EIGEN_STRONG_INLINE Packet4i ploadu(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(const int* from) { EIGEN_DEBUG_ALIGNED_LOAD @@ -420,9 +394,6 @@ template<> EIGEN_STRONG_INLINE Packet4i ploaddup(const int* from) return vec_perm(p, p, p16uc_DUPLICATE); } -template<> EIGEN_STRONG_INLINE void pstore(float* to, const Packet4f& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st(from, 0, to); } -template<> EIGEN_STRONG_INLINE void pstore(int* to, const Packet4i& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st(from, 0, to); } - template<> EIGEN_STRONG_INLINE void pstoreu(float* to, const Packet4f& from) { EIGEN_DEBUG_UNALIGNED_STORE @@ -642,6 +613,252 @@ ptranspose(PacketBlock& 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 : 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 { 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(const double* from) { EIGEN_DEBUG_ALIGNED_LOAD return (Packet2d) vec_ld(0, (const float *) from); } //FIXME + +template<> EIGEN_STRONG_INLINE void pstore(double* to, const Packet2d& from) { EIGEN_DEBUG_ALIGNED_STORE vec_st((Packet4f)from, 0, (float *)to); } + +template<> EIGEN_STRONG_INLINE Packet2d pset1(const double& from) { + double EIGEN_ALIGN16 af[2]; + af[0] = from; + Packet2d vc = pload(af); + vc = vec_splat_dbl(vc, 0); + return vc; +} +template<> EIGEN_STRONG_INLINE void +pbroadcast4(const double *a, + Packet2d& a0, Packet2d& a1, Packet2d& a2, Packet2d& a3) +{ + a1 = pload(a); + a0 = vec_splat_dbl(a1, 0); + a1 = vec_splat_dbl(a1, 1); + a3 = pload(a+2); + a2 = vec_splat_dbl(a3, 0); + a3 = vec_splat_dbl(a1, 1); +} +template<> EIGEN_DEVICE_FUNC inline Packet2d pgather(const double* from, DenseIndex stride) +{ + double EIGEN_ALIGN16 af[2]; + af[0] = from[0*stride]; + af[1] = from[1*stride]; + return pload(af); +} +template<> EIGEN_DEVICE_FUNC inline void pscatter(double* to, const Packet2d& from, DenseIndex stride) +{ + double EIGEN_ALIGN16 af[2]; + pstore(af, from); + to[0*stride] = af[0]; + to[1*stride] = af[1]; +} +template<> EIGEN_STRONG_INLINE Packet2d plset(const double& a) { return vec_add(pset1(a), p2d_COUNTDOWN); } + +template<> EIGEN_STRONG_INLINE Packet2d padd(const Packet2d& a, const Packet2d& b) { return vec_add(a,b); } + +template<> EIGEN_STRONG_INLINE Packet2d psub(const Packet2d& a, const Packet2d& b) { return vec_sub(a,b); } + +template<> EIGEN_STRONG_INLINE Packet2d pnegate(const Packet2d& a) { return psub(p2d_ZERO, a); } + +template<> EIGEN_STRONG_INLINE Packet2d pconj(const Packet2d& a) { return a; } + +template<> EIGEN_STRONG_INLINE Packet2d pmul(const Packet2d& a, const Packet2d& b) { return vec_madd(a,b,p2d_ZERO); } +template<> EIGEN_STRONG_INLINE Packet2d pdiv(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(const Packet2d& a, const Packet2d& b) { return vec_min(a, b); } + +template<> EIGEN_STRONG_INLINE Packet2d pmax(const Packet2d& a, const Packet2d& b) { return vec_max(a, b); } + +template<> EIGEN_STRONG_INLINE Packet2d pand(const Packet2d& a, const Packet2d& b) { return vec_and(a, b); } + +template<> EIGEN_STRONG_INLINE Packet2d por(const Packet2d& a, const Packet2d& b) { return vec_or(a, b); } + +template<> EIGEN_STRONG_INLINE Packet2d pxor(const Packet2d& a, const Packet2d& b) { return vec_xor(a, b); } + +template<> EIGEN_STRONG_INLINE Packet2d pandnot(const Packet2d& a, const Packet2d& b) { return vec_and(a, vec_nor(b, b)); } + + +template<> EIGEN_STRONG_INLINE Packet2d ploadu(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(const double* from) +{ + Packet2d p; + if((ptrdiff_t(from) % 16) == 0) p = pload(from); + else p = ploadu(from); + return vec_perm(p, p, p16uc_DUPLICATE); +} + +template<> EIGEN_STRONG_INLINE void pstoreu(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(const double* addr) { vec_dstt((const float *) addr, DST_CTRL(2,2,32), DST_CHAN); } + +template<> EIGEN_STRONG_INLINE double pfirst(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(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(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(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(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(const Packet2d& a) +{ + Packet2d res; + res = vec_max(a, (Packet2d) vec_sld((Packet4ui) a, (Packet4ui) a, 8)); + return pfirst(res); +} + +template +struct palign_impl +{ + 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& 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 From 60663a510a9d1c1dfb3518e95b495170ab4d66ba Mon Sep 17 00:00:00 2001 From: Konstantinos Margaritis Date: Fri, 19 Sep 2014 21:05:01 +0000 Subject: [PATCH 3/8] 32-bit floats/ints, 64-bit doubles pass packetmath tests, complex 32/64-bit remaining --- Eigen/src/Core/arch/AltiVec/Complex.h | 61 ++++--- Eigen/src/Core/arch/AltiVec/PacketMath.h | 221 ++++++++++++++--------- 2 files changed, 171 insertions(+), 111 deletions(-) diff --git a/Eigen/src/Core/arch/AltiVec/Complex.h b/Eigen/src/Core/arch/AltiVec/Complex.h index d90cfe0de..94ecec2de 100644 --- a/Eigen/src/Core/arch/AltiVec/Complex.h +++ b/Eigen/src/Core/arch/AltiVec/Complex.h @@ -7,8 +7,8 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. -#ifndef EIGEN_COMPLEX_ALTIVEC_H -#define EIGEN_COMPLEX_ALTIVEC_H +#ifndef EIGEN_COMPLEX32_ALTIVEC_H +#define EIGEN_COMPLEX32_ALTIVEC_H namespace Eigen { @@ -21,16 +21,6 @@ static Packet4ui p4ui_CONJ_XOR = vec_mergeh((Packet4ui)p4f_ZERO_, (Packet4ui)p4 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 }; -static Packet16uc p16uc_COMPLEX_IM = vec_sld(p16uc_DUPLICATE, (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 3), 8);//{ 4,5,6,7, 4,5,6,7, 12,13,14,15, 12,13,14,15 }; -static Packet16uc p16uc_COMPLEX_REV = vec_sld(p16uc_REVERSE, p16uc_REVERSE, 8);//{ 4,5,6,7, 0,1,2,3, 12,13,14,15, 8,9,10,11 }; -static Packet16uc p16uc_COMPLEX_REV2 = vec_sld(p16uc_FORWARD, p16uc_FORWARD, 8);//{ 8,9,10,11, 12,13,14,15, 0,1,2,3, 4,5,6,7 }; -static Packet16uc p16uc_PSET_HI = (Packet16uc) vec_mergeh((Packet4ui)p16uc_COMPLEX_RE, (Packet4ui)p16uc_COMPLEX_IM);//{ 0,1,2,3, 4,5,6,7, 0,1,2,3, 4,5,6,7 }; -static Packet16uc p16uc_PSET_LO = (Packet16uc) vec_mergel((Packet4ui)p16uc_COMPLEX_RE, (Packet4ui)p16uc_COMPLEX_IM);//{ 8,9,10,11, 12,13,14,15, 8,9,10,11, 12,13,14,15 }; -static Packet16uc p16uc_COMPLEX_MASK16 = 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}; -static Packet16uc p16uc_COMPLEX_TRANSPOSE_0 = vec_add(p16uc_PSET_HI, p16uc_COMPLEX_MASK16);//{ 0,1,2,3, 4,5,6,7, 16,17,18,19, 20,21,22,23}; -static Packet16uc p16uc_COMPLEX_TRANSPOSE_1 = vec_add(p16uc_PSET_LO, p16uc_COMPLEX_MASK16);//{ 8,9,10,11, 12,13,14,15, 24,25,26,27, 28,29,30,31}; - //---------- float ---------- struct Packet2cf { @@ -71,7 +61,7 @@ template<> EIGEN_STRONG_INLINE Packet2cf pset1(const std::complex((const float *)&from); else res.v = ploadu((const float *)&from); - res.v = vec_perm(res.v, res.v, p16uc_PSET_HI); + res.v = vec_perm(res.v, res.v, p16uc_PSET64_HI); return res; } @@ -100,17 +90,24 @@ template<> EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, con { Packet4f v1, v2; + std::cout << "a: " << a.v << ", b: " << b.v << std::endl; // Permute and multiply the real parts of a and b - v1 = vec_perm(a.v, a.v, p16uc_COMPLEX_RE); + v1 = vec_perm(a.v, a.v, p16uc_PSET32_WODD); + std::cout << "v1 = a(re): " << v1 << std::endl; // Get the imaginary parts of a - v2 = vec_perm(a.v, a.v, p16uc_COMPLEX_IM); + v2 = vec_perm(a.v, a.v, p16uc_PSET32_WEVEN); + std::cout << "v2 = a(im): " << v2 << std::endl; // multiply a_re * b v1 = vec_madd(v1, b.v, p4f_ZERO); + std::cout << "v1+ b : " << v1 << std::endl; // multiply a_im * b and get the conjugate result v2 = vec_madd(v2, b.v, p4f_ZERO); + std::cout << "v2+ b : " << v2 << std::endl; v2 = (Packet4f) vec_xor((Packet4ui)v2, p4ui_CONJ_XOR); + std::cout << "~v2 : " << v2 << std::endl; // permute back to a proper order - v2 = vec_perm(v2, v2, p16uc_COMPLEX_REV); + v2 = vec_perm(v2, v2, p16uc_COMPLEX32_REV); + std::cout << "v2 : " << v2 << std::endl; return Packet2cf(vec_add(v1, v2)); } @@ -144,7 +141,7 @@ template<> EIGEN_STRONG_INLINE std::complex pfirst(const Pack template<> EIGEN_STRONG_INLINE Packet2cf preverse(const Packet2cf& a) { Packet4f rev_a; - rev_a = vec_perm(a.v, a.v, p16uc_COMPLEX_REV2); + rev_a = vec_perm(a.v, a.v, p16uc_COMPLEX32_REV2); return Packet2cf(rev_a); } @@ -185,7 +182,11 @@ struct palign_impl { if (Offset==1) { +#ifdef _BIG_ENDIAN first.v = vec_sld(first.v, second.v, 8); +#else + first.v = vec_sld(second.v, first.v, 8); +#endif } } }; @@ -228,18 +229,18 @@ template<> EIGEN_STRONG_INLINE Packet2cf pdiv(const Packet2cf& a, con // TODO optimize it for AltiVec Packet2cf res = conj_helper().pmul(a,b); Packet4f s = vec_madd(b.v, b.v, p4f_ZERO); - return Packet2cf(pdiv(res.v, vec_add(s,vec_perm(s, s, p16uc_COMPLEX_REV)))); + return Packet2cf(pdiv(res.v, vec_add(s,vec_perm(s, s, p16uc_COMPLEX32_REV)))); } template<> EIGEN_STRONG_INLINE Packet2cf pcplxflip(const Packet2cf& x) { - return Packet2cf(vec_perm(x.v, x.v, p16uc_COMPLEX_REV)); + return Packet2cf(vec_perm(x.v, x.v, p16uc_COMPLEX32_REV)); } EIGEN_STRONG_INLINE void ptranspose(PacketBlock& kernel) { - Packet4f tmp = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_COMPLEX_TRANSPOSE_0); - kernel.packet[1].v = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_COMPLEX_TRANSPOSE_1); + Packet4f tmp = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_TRANSPOSE64_HI); + kernel.packet[1].v = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_TRANSPOSE64_LO); kernel.packet[0].v = tmp; } @@ -289,7 +290,7 @@ template<> EIGEN_STRONG_INLINE Packet1cd pset1(const std::complex((const double *)&from); else res.v = ploadu((const double *)&from); - res.v = vec_perm(res.v, res.v, p16uc_PSET_HI); + res.v = vec_perm(res.v, res.v, p16uc_PSET64_HI); return res; } template<> EIGEN_DEVICE_FUNC inline Packet1cd pgather, Packet1cd>(const std::complex* from, DenseIndex stride) @@ -317,9 +318,9 @@ template<> EIGEN_STRONG_INLINE Packet1cd pmul(const Packet1cd& a, con Packet2d a_re, a_im, v1, v2; // Permute and multiply the real parts of a and b - a_re = vec_perm(a.v, a.v, p16uc_PSET_HI); + a_re = vec_perm(a.v, a.v, p16uc_PSET64_HI); // Get the imaginary parts of a - a_im = vec_perm(a.v, a.v, p16uc_PSET_HI); + a_im = vec_perm(a.v, a.v, p16uc_PSET64_HI); // multiply a_re * b v1 = vec_madd(a_re, b.v, p2d_ZERO_); // multiply a_im * b and get the conjugate result @@ -353,7 +354,7 @@ template<> EIGEN_STRONG_INLINE std::complex pfirst(const Pac template<> EIGEN_STRONG_INLINE Packet1cd preverse(const Packet1cd& a) { Packet2d rev_a; - rev_a = vec_perm(a.v, a.v, p16uc_COMPLEX_REV2); + rev_a = vec_perm(a.v, a.v, p16uc_COMPLEX32_REV2); return Packet1cd(rev_a); } @@ -432,18 +433,18 @@ template<> EIGEN_STRONG_INLINE Packet1cd pdiv(const Packet1cd& a, con // TODO optimize it for AltiVec Packet1cd res = conj_helper().pmul(a,b); 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)))); + return Packet1cd(pdiv(res.v, vec_add(s,vec_perm(s, s, p16uc_COMPLEX32_REV)))); } template<> EIGEN_STRONG_INLINE Packet1cd pcplxflip(const Packet1cd& x) { - return Packet1cd(vec_perm(x.v, x.v, p16uc_COMPLEX_REV)); + return Packet1cd(vec_perm(x.v, x.v, p16uc_COMPLEX32_REV)); } EIGEN_STRONG_INLINE void ptranspose(PacketBlock& kernel) { - Packet2d tmp = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_COMPLEX_TRANSPOSE_0); - kernel.packet[1].v = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_COMPLEX_TRANSPOSE_1); + Packet2d tmp = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_TRANSPOSE64_HI); + kernel.packet[1].v = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_TRANSPOSE64_LO); kernel.packet[0].v = tmp; } @@ -451,4 +452,4 @@ EIGEN_STRONG_INLINE void ptranspose(PacketBlock& kernel) } // end namespace Eigen -#endif // EIGEN_COMPLEX_ALTIVEC_H +#endif // EIGEN_COMPLEX32_ALTIVEC_H diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h index 0489ad886..aea3e79dc 100755 --- a/Eigen/src/Core/arch/AltiVec/PacketMath.h +++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h @@ -62,21 +62,6 @@ typedef __vector unsigned char Packet16uc; #define DST_CHAN 1 #define DST_CTRL(size, count, stride) (((size) << 24) | ((count) << 16) | (stride)) -// Handle endianness properly while loading constants -// Define global static constants: -#ifdef _BIG_ENDIAN -static Packet4f p4f_COUNTDOWN = { 0.0, 1.0, 2.0, 3.0 }; -static Packet4i p4i_COUNTDOWN = { 0, 1, 2, 3 }; -static Packet16uc p16uc_REVERSE = { 12,13,14,15, 8,9,10,11, 4,5,6,7, 0,1,2,3}; -static Packet16uc p16uc_FORWARD = vec_lvsl(0, (float*)0); //{ 0,1,2,3, 4,5,6,7, 8,9,10,11, 12,13,14,15} -static Packet16uc p16uc_DUPLICATE = { 0,1,2,3, 0,1,2,3, 4,5,6,7, 4,5,6,7}; -#else -static Packet4f p4f_COUNTDOWN = { 3.0, 2.0, 1.0, 0.0 }; -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 }; -#endif // _BIG_ENDIAN // These constants are endian-agnostic static _EIGEN_DECLARE_CONST_FAST_Packet4f(ZERO, 0); //{ 0.0, 0.0, 0.0, 0.0} @@ -87,6 +72,40 @@ static _EIGEN_DECLARE_CONST_FAST_Packet4i(MINUS1,-1); //{ -1, -1, -1, -1} static Packet4f p4f_ONE = vec_ctf(p4i_ONE, 0); //{ 1.0, 1.0, 1.0, 1.0} static Packet4f p4f_ZERO_ = (Packet4f) vec_sl((Packet4ui)p4i_MINUS1, (Packet4ui)p4i_MINUS1); //{ 0x80000000, 0x80000000, 0x80000000, 0x80000000} +static Packet4f p4f_COUNTDOWN = { 0.0, 1.0, 2.0, 3.0 }; +static Packet4i p4i_COUNTDOWN = { 0, 1, 2, 3 }; + +static Packet16uc p16uc_REVERSE32 = { 12,13,14,15, 8,9,10,11, 4,5,6,7, 0,1,2,3 }; +static Packet16uc p16uc_REVERSE64 = { 8,9,10,11, 12,13,14,15, 0,1,2,3, 4,5,6,7 }; +static Packet16uc p16uc_DUPLICATE32_HI = { 0,1,2,3, 0,1,2,3, 4,5,6,7, 4,5,6,7 }; + +// Handle endianness properly while loading constants +// Define global static constants: +#ifdef _BIG_ENDIAN +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 }; +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 }; +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}; +#else +// Mask alignment +#define _EIGEN_MASK_ALIGNMENT 0xfffffffffffffff0 +#define _EIGEN_ALIGNED_PTR(x) ((ptrdiff_t)(x) & _EIGEN_MASK_ALIGNMENT) + +static Packet16uc p16uc_FORWARD = p16uc_REVERSE32; +static Packet16uc p16uc_PSET32_WODD = vec_sld((Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 1), (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 3), 8);//{ 0,1,2,3, 0,1,2,3, 8,9,10,11, 8,9,10,11 }; +static Packet16uc p16uc_PSET32_WEVEN = vec_sld((Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 0), (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 2), 8);//{ 4,5,6,7, 4,5,6,7, 12,13,14,15, 12,13,14,15 }; +static Packet16uc p16uc_HALF64_0_16 = vec_sld(vec_splat((Packet16uc) vec_abs(p4i_MINUS16), 0), (Packet16uc)p4i_ZERO, 8); //{ 0,0,0,0, 0,0,0,0, 16,16,16,16, 16,16,16,16}; + +#endif // _BIG_ENDIAN + +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 }; +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 }; + +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 }; +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 }; + +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}; +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}; + template<> struct packet_traits : default_packet_traits { typedef Packet4f type; @@ -121,7 +140,19 @@ template<> struct packet_traits : default_packet_traits template<> struct unpacket_traits { typedef float type; enum {size=4}; typedef Packet4f half; }; template<> struct unpacket_traits { typedef int type; enum {size=4}; typedef Packet4i half; }; -/* + +inline std::ostream & operator <<(std::ostream & s, const Packet16uc & v) +{ + union { + Packet16uc v; + unsigned char n[16]; + } vt; + vt.v = v; + for (int i=0; i< 16; i++) + s << (int)vt.n[i] << ", "; + return s; +} + inline std::ostream & operator <<(std::ostream & s, const Packet4f & v) { union { @@ -154,7 +185,7 @@ inline std::ostream & operator <<(std::ostream & s, const Packet4ui & v) s << vt.n[0] << ", " << vt.n[1] << ", " << vt.n[2] << ", " << vt.n[3]; return s; } - +/* inline std::ostream & operator <<(std::ostream & s, const Packetbi & v) { union { @@ -164,8 +195,8 @@ inline std::ostream & operator <<(std::ostream & s, const Packetbi & v) vt.v = v; s << vt.n[0] << ", " << vt.n[1] << ", " << vt.n[2] << ", " << vt.n[3]; return s; -} -*/ +}*/ + // Need to define them first or we get specialization after instantiation errors template<> EIGEN_STRONG_INLINE Packet4f pload(const float* from) { EIGEN_DEBUG_ALIGNED_LOAD return vec_ld(0, from); } @@ -370,12 +401,12 @@ template<> EIGEN_STRONG_INLINE Packet4i ploadu(const int* from) template<> EIGEN_STRONG_INLINE Packet4i ploadu(const int* from) { EIGEN_DEBUG_ALIGNED_LOAD - return (Packet4i) vec_vsx_ld((long)from & 15, from); // align the data + return (Packet4i) vec_vsx_ld((long)from & 15, (const int*) _EIGEN_ALIGNED_PTR(from)); } template<> EIGEN_STRONG_INLINE Packet4f ploadu(const float* from) { EIGEN_DEBUG_ALIGNED_LOAD - return (Packet4f) vec_vsx_ld((long)from & 15, from); // align the data + return (Packet4f) vec_vsx_ld((long)from & 15, (const float*) _EIGEN_ALIGNED_PTR(from)); } #endif @@ -384,16 +415,17 @@ template<> EIGEN_STRONG_INLINE Packet4f ploaddup(const float* from) Packet4f p; if((ptrdiff_t(from) % 16) == 0) p = pload(from); else p = ploadu(from); - return vec_perm(p, p, p16uc_DUPLICATE); + return vec_perm(p, p, p16uc_DUPLICATE32_HI); } template<> EIGEN_STRONG_INLINE Packet4i ploaddup(const int* from) { Packet4i p; if((ptrdiff_t(from) % 16) == 0) p = pload(from); else p = ploadu(from); - return vec_perm(p, p, p16uc_DUPLICATE); + return vec_perm(p, p, p16uc_DUPLICATE32_HI); } +#ifdef _BIG_ENDIAN template<> EIGEN_STRONG_INLINE void pstoreu(float* to, const Packet4f& from) { EIGEN_DEBUG_UNALIGNED_STORE @@ -430,6 +462,19 @@ template<> EIGEN_STRONG_INLINE void pstoreu(int* to, const Packet4i& f vec_st( LSQ, 15, (unsigned char *)to ); // Store the LSQ part first vec_st( MSQ, 0, (unsigned char *)to ); // Store the MSQ part } +#else +// We also need ot redefine little endian loading of Packet4i/Packet4f using VSX +template<> EIGEN_STRONG_INLINE void pstoreu(int* to, const Packet4i& from) +{ + EIGEN_DEBUG_ALIGNED_STORE + vec_vsx_st(from, (long)to & 15, (int*) _EIGEN_ALIGNED_PTR(to)); +} +template<> EIGEN_STRONG_INLINE void pstoreu(float* to, const Packet4f& from) +{ + EIGEN_DEBUG_ALIGNED_STORE + vec_vsx_st(from, (long)to & 15, (float*) _EIGEN_ALIGNED_PTR(to)); +} +#endif template<> EIGEN_STRONG_INLINE void prefetch(const float* addr) { vec_dstt(addr, DST_CTRL(2,2,32), DST_CHAN); } template<> EIGEN_STRONG_INLINE void prefetch(const int* addr) { vec_dstt(addr, DST_CTRL(2,2,32), DST_CHAN); } @@ -437,8 +482,8 @@ template<> EIGEN_STRONG_INLINE void prefetch(const int* addr) { vec_dst template<> EIGEN_STRONG_INLINE float pfirst(const Packet4f& a) { float EIGEN_ALIGN16 x[4]; vec_st(a, 0, x); return x[0]; } template<> EIGEN_STRONG_INLINE int pfirst(const Packet4i& a) { int EIGEN_ALIGN16 x[4]; vec_st(a, 0, x); return x[0]; } -template<> EIGEN_STRONG_INLINE Packet4f preverse(const Packet4f& a) { return (Packet4f)vec_perm((Packet16uc)a,(Packet16uc)a, p16uc_REVERSE); } -template<> EIGEN_STRONG_INLINE Packet4i preverse(const Packet4i& a) { return (Packet4i)vec_perm((Packet16uc)a,(Packet16uc)a, p16uc_REVERSE); } +template<> EIGEN_STRONG_INLINE Packet4f preverse(const Packet4f& a) { return (Packet4f)vec_perm((Packet16uc)a,(Packet16uc)a, p16uc_REVERSE32); } +template<> EIGEN_STRONG_INLINE Packet4i preverse(const Packet4i& a) { return (Packet4i)vec_perm((Packet16uc)a,(Packet16uc)a, p16uc_REVERSE32); } template<> EIGEN_STRONG_INLINE Packet4f pabs(const Packet4f& a) { return vec_abs(a); } template<> EIGEN_STRONG_INLINE Packet4i pabs(const Packet4i& a) { return vec_abs(a); } @@ -485,7 +530,11 @@ template<> EIGEN_STRONG_INLINE int predux(const Packet4i& a) { Packet4i sum; sum = vec_sums(a, p4i_ZERO); +#ifdef _BIG_ENDIAN sum = vec_sld(sum, p4i_ZERO, 12); +#else + sum = vec_sld(p4i_ZERO, sum, 4); +#endif return pfirst(sum); } @@ -572,8 +621,25 @@ struct palign_impl { static EIGEN_STRONG_INLINE void run(Packet4f& first, const Packet4f& second) { - if (Offset!=0) - first = vec_sld(first, second, Offset*4); +#ifdef _BIG_ENDIAN + switch (Offset % 4) { + case 1: + first = vec_sld(first, second, 4); break; + case 2: + first = vec_sld(first, second, 8); break; + case 3: + first = vec_sld(first, second, 12); break; + } +#else + switch (Offset % 4) { + case 1: + first = vec_sld(second, first, 12); break; + case 2: + first = vec_sld(second, first, 8); break; + case 3: + first = vec_sld(second, first, 4); break; + } +#endif } }; @@ -582,8 +648,25 @@ struct palign_impl { static EIGEN_STRONG_INLINE void run(Packet4i& first, const Packet4i& second) { - if (Offset!=0) - first = vec_sld(first, second, Offset*4); +#ifdef _BIG_ENDIAN + switch (Offset % 4) { + case 1: + first = vec_sld(first, second, 4); break; + case 2: + first = vec_sld(first, second, 8); break; + case 3: + first = vec_sld(first, second, 12); break; + } +#else + switch (Offset % 4) { + case 1: + first = vec_sld(second, first, 12); break; + case 2: + first = vec_sld(second, first, 8); break; + case 3: + first = vec_sld(second, first, 4); break; + } +#endif } }; @@ -625,24 +708,15 @@ 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 Packet2d p2d_COUNTDOWN = { 0.0, 1.0 }; static EIGEN_STRONG_INLINE Packet2d vec_splat_dbl(Packet2d& a, int index) { switch (index) { case 0: - return (Packet2d) vec_perm(a, a, p16uc_SPLATQ1); + return (Packet2d) vec_perm(a, a, p16uc_PSET64_HI); case 1: - return (Packet2d) vec_perm(a, a, p16uc_SPLATQ2); + return (Packet2d) vec_perm(a, a, p16uc_PSET64_LO); } return a; } @@ -665,17 +739,17 @@ template<> struct packet_traits : default_packet_traits template<> struct unpacket_traits { typedef double type; enum {size=2}; typedef Packet2d half; }; -/* + inline std::ostream & operator <<(std::ostream & s, const Packet2d & v) { union { - Packetdf v; + Packet2d 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(const double* from) { EIGEN_DEBUG_ALIGNED_LOAD return (Packet2d) vec_ld(0, (const float *) from); } //FIXME @@ -698,7 +772,7 @@ pbroadcast4(const double *a, a1 = vec_splat_dbl(a1, 1); a3 = pload(a+2); a2 = vec_splat_dbl(a3, 0); - a3 = vec_splat_dbl(a1, 1); + a3 = vec_splat_dbl(a3, 1); } template<> EIGEN_DEVICE_FUNC inline Packet2d pgather(const double* from, DenseIndex stride) { @@ -755,43 +829,30 @@ template<> EIGEN_STRONG_INLINE Packet2d pxor(const Packet2d& a, const template<> EIGEN_STRONG_INLINE Packet2d pandnot(const Packet2d& a, const Packet2d& b) { return vec_and(a, vec_nor(b, b)); } - template<> EIGEN_STRONG_INLINE Packet2d ploadu(const double* from) { EIGEN_DEBUG_ALIGNED_LOAD - return (Packet2d) vec_vsx_ld((int)((size_t)(from) & 15), (const float *)from); // align the data + return (Packet2d) vec_vsx_ld((long)from & 15, (const float*) _EIGEN_ALIGNED_PTR(from)); } template<> EIGEN_STRONG_INLINE Packet2d ploaddup(const double* from) { Packet2d p; if((ptrdiff_t(from) % 16) == 0) p = pload(from); else p = ploadu(from); - return vec_perm(p, p, p16uc_DUPLICATE); + return vec_perm(p, p, p16uc_PSET64_HI); } template<> EIGEN_STRONG_INLINE void pstoreu(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 + EIGEN_DEBUG_ALIGNED_STORE + vec_vsx_st((Packet4f)from, (long)to & 15, (float*) _EIGEN_ALIGNED_PTR(to)); } + template<> EIGEN_STRONG_INLINE void prefetch(const double* addr) { vec_dstt((const float *) addr, DST_CTRL(2,2,32), DST_CHAN); } template<> EIGEN_STRONG_INLINE double pfirst(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 preverse(const Packet2d& a) { return (Packet2d)vec_perm((Packet16uc)a,(Packet16uc)a, p16uc_REVERSE64); } template<> EIGEN_STRONG_INLINE Packet2d pabs(const Packet2d& a) { return vec_abs(a); } @@ -807,9 +868,10 @@ template<> EIGEN_STRONG_INLINE Packet2d preduxp(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]); + v[0] = vec_add(vecs[0], (Packet2d) vec_sld((Packet4ui) vecs[0], (Packet4ui) vecs[0], 8)); + v[1] = vec_add(vecs[1], (Packet2d) vec_sld((Packet4ui) vecs[1], (Packet4ui) vecs[1], 8)); + + sum = (Packet2d) vec_sld((Packet4ui) v[1], (Packet4ui) v[0], 8); return sum; } @@ -817,26 +879,19 @@ template<> EIGEN_STRONG_INLINE Packet2d preduxp(const Packet2d* vecs) // mul template<> EIGEN_STRONG_INLINE double predux_mul(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))); + return pfirst(pmul(a, (Packet2d)vec_sld((Packet4ui) a, (Packet4ui) a, 8))); } // min template<> EIGEN_STRONG_INLINE double predux_min(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); + return pfirst(vec_min(a, (Packet2d) vec_sld((Packet4ui) a, (Packet4ui) a, 8))); } // max template<> EIGEN_STRONG_INLINE double predux_max(const Packet2d& a) { - Packet2d res; - res = vec_max(a, (Packet2d) vec_sld((Packet4ui) a, (Packet4ui) a, 8)); - return pfirst(res); + return pfirst(vec_max(a, (Packet2d) vec_sld((Packet4ui) a, (Packet4ui) a, 8))); } template @@ -844,16 +899,20 @@ struct palign_impl { static EIGEN_STRONG_INLINE void run(Packet2d& first, const Packet2d& second) { - if (Offset!=0) - first = (Packet2d) vec_sld((Packet4ui) first, (Packet4ui) second, Offset*8); + if (Offset == 1) +#ifdef _BIG_ENDIAN + first = (Packet2d) vec_sld((Packet4ui) first, (Packet4ui) second, 8); +#else + first = (Packet2d) vec_sld((Packet4ui) second, (Packet4ui) first, 8); +#endif } }; EIGEN_DEVICE_FUNC inline void ptranspose(PacketBlock& 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); + t0 = vec_perm(kernel.packet[0], kernel.packet[1], p16uc_TRANSPOSE64_HI); + t1 = vec_perm(kernel.packet[0], kernel.packet[1], p16uc_TRANSPOSE64_LO); kernel.packet[0] = t0; kernel.packet[1] = t1; } From 10f8aabb611e7b727ff67141b030689c0a53ca78 Mon Sep 17 00:00:00 2001 From: Konstantinos Margaritis Date: Sat, 20 Sep 2014 22:31:31 +0000 Subject: [PATCH 4/8] VSX port passes packetmath_[1-5] tests! --- Eigen/src/Core/arch/AltiVec/Complex.h | 80 ++++++++---------------- Eigen/src/Core/arch/AltiVec/PacketMath.h | 22 ++++--- 2 files changed, 38 insertions(+), 64 deletions(-) diff --git a/Eigen/src/Core/arch/AltiVec/Complex.h b/Eigen/src/Core/arch/AltiVec/Complex.h index 94ecec2de..58f5fbd27 100644 --- a/Eigen/src/Core/arch/AltiVec/Complex.h +++ b/Eigen/src/Core/arch/AltiVec/Complex.h @@ -14,11 +14,11 @@ namespace Eigen { namespace internal { -#ifdef _BIG_ENDIAN static Packet4ui p4ui_CONJ_XOR = vec_mergeh((Packet4ui)p4i_ZERO, (Packet4ui)p4f_ZERO_);//{ 0x00000000, 0x80000000, 0x00000000, 0x80000000 }; +#ifdef _BIG_ENDIAN #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_XOR1 = (Packet2ul) vec_sld((Packet4ui) p2l_ZERO, (Packet4ui) p2d_ZERO_, 8);//{ 0x8000000000000000, 0x0000000000000000 }; +static Packet2ul p2ul_CONJ_XOR2 = (Packet2ul) vec_sld((Packet4ui) p2d_ZERO_, (Packet4ui) p2l_ZERO, 8);//{ 0x8000000000000000, 0x0000000000000000 }; #endif //---------- float ---------- @@ -90,24 +90,17 @@ template<> EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, con { Packet4f v1, v2; - std::cout << "a: " << a.v << ", b: " << b.v << std::endl; // Permute and multiply the real parts of a and b v1 = vec_perm(a.v, a.v, p16uc_PSET32_WODD); - std::cout << "v1 = a(re): " << v1 << std::endl; // Get the imaginary parts of a v2 = vec_perm(a.v, a.v, p16uc_PSET32_WEVEN); - std::cout << "v2 = a(im): " << v2 << std::endl; // multiply a_re * b v1 = vec_madd(v1, b.v, p4f_ZERO); - std::cout << "v1+ b : " << v1 << std::endl; // multiply a_im * b and get the conjugate result v2 = vec_madd(v2, b.v, p4f_ZERO); - std::cout << "v2+ b : " << v2 << std::endl; v2 = (Packet4f) vec_xor((Packet4ui)v2, p4ui_CONJ_XOR); - std::cout << "~v2 : " << v2 << std::endl; // permute back to a proper order v2 = vec_perm(v2, v2, p16uc_COMPLEX32_REV); - std::cout << "v2 : " << v2 << std::endl; return Packet2cf(vec_add(v1, v2)); } @@ -156,9 +149,13 @@ template<> EIGEN_STRONG_INLINE std::complex predux(const Packe template<> EIGEN_STRONG_INLINE Packet2cf preduxp(const Packet2cf* vecs) { Packet4f b1, b2; - +#ifdef _BIG_ENDIAN b1 = (Packet4f) vec_sld(vecs[0].v, vecs[1].v, 8); b2 = (Packet4f) vec_sld(vecs[1].v, vecs[0].v, 8); +#else + b1 = (Packet4f) vec_sld(vecs[1].v, vecs[0].v, 8); + b2 = (Packet4f) vec_sld(vecs[0].v, vecs[1].v, 8); +#endif b2 = (Packet4f) vec_sld(b2, b2, 8); b2 = padd(b1, b2); @@ -283,16 +280,8 @@ template<> EIGEN_STRONG_INLINE void pstore >(std::complex< template<> EIGEN_STRONG_INLINE void pstoreu >(std::complex * to, const Packet1cd& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, from.v); } template<> EIGEN_STRONG_INLINE Packet1cd pset1(const std::complex& 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((const double *)&from); - else - res.v = ploadu((const double *)&from); - res.v = vec_perm(res.v, res.v, p16uc_PSET64_HI); - return res; -} +{ /* here we really have to use unaligned loads :( */ return ploadu(&from); } + template<> EIGEN_DEVICE_FUNC inline Packet1cd pgather, Packet1cd>(const std::complex* from, DenseIndex stride) { std::complex EIGEN_ALIGN16 af[2]; @@ -311,7 +300,7 @@ template<> EIGEN_DEVICE_FUNC inline void pscatter, Packet1c template<> EIGEN_STRONG_INLINE Packet1cd padd(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_add(a.v,b.v)); } template<> EIGEN_STRONG_INLINE Packet1cd psub(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((Packet2d)a.v, (Packet2d)p2ul_CONJ_XOR)); } +template<> EIGEN_STRONG_INLINE Packet1cd pconj(const Packet1cd& a) { return Packet1cd((Packet2d)vec_xor((Packet2d)a.v, (Packet2d)p2ul_CONJ_XOR2)); } template<> EIGEN_STRONG_INLINE Packet1cd pmul(const Packet1cd& a, const Packet1cd& b) { @@ -320,13 +309,13 @@ template<> EIGEN_STRONG_INLINE Packet1cd pmul(const Packet1cd& a, con // Permute and multiply the real parts of a and b a_re = vec_perm(a.v, a.v, p16uc_PSET64_HI); // Get the imaginary parts of a - a_im = vec_perm(a.v, a.v, p16uc_PSET64_HI); + a_im = vec_perm(a.v, a.v, p16uc_PSET64_LO); // multiply a_re * b - v1 = vec_madd(a_re, b.v, p2d_ZERO_); + v1 = vec_madd(a_re, b.v, p2d_ZERO); // multiply a_im * b and get the conjugate result - v2 = vec_madd(a_im, b.v, p2d_ZERO_); + 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); + v2 = (Packet2d) vec_xor((Packet2d)v2, (Packet2d) p2ul_CONJ_XOR1); return Packet1cd(vec_add(v1, v2)); } @@ -351,48 +340,31 @@ template<> EIGEN_STRONG_INLINE std::complex pfirst(const Pac return res[0]; } -template<> EIGEN_STRONG_INLINE Packet1cd preverse(const Packet1cd& a) -{ - Packet2d rev_a; - rev_a = vec_perm(a.v, a.v, p16uc_COMPLEX32_REV2); - return Packet1cd(rev_a); -} +template<> EIGEN_STRONG_INLINE Packet1cd preverse(const Packet1cd& a) { return a; } template<> EIGEN_STRONG_INLINE std::complex predux(const Packet1cd& a) { - Packet2d b; - b = (Packet2d) vec_sld((Packet4ui) a.v, (Packet4ui) a.v, 8); - b = padd(a.v, b); - return pfirst(Packet1cd(b)); + return pfirst(a); } template<> EIGEN_STRONG_INLINE Packet1cd preduxp(const Packet1cd* vecs) { - Packet2d b1, b2; - - 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); + return vecs[0]; } template<> EIGEN_STRONG_INLINE std::complex predux_mul(const Packet1cd& a) { - Packet2d b; - Packet1cd prod; - b = (Packet2d) vec_sld((Packet4ui) a.v, (Packet4ui) a.v, 8); - prod = pmul(a, Packet1cd(b)); - - return pfirst(prod); + return pfirst(a); } template struct palign_impl { - static EIGEN_STRONG_INLINE void run(Packet1cd& first, const Packet1cd& second) - { } + static EIGEN_STRONG_INLINE void run(Packet1cd& /*first*/, const Packet1cd& /*second*/) + { + // FIXME is it sure we never have to align a Packet1cd? + // Even though a std::complex has 16 bytes, it is not necessarily aligned on a 16 bytes boundary... + } }; template<> struct conj_helper @@ -436,9 +408,9 @@ template<> EIGEN_STRONG_INLINE Packet1cd pdiv(const Packet1cd& a, con return Packet1cd(pdiv(res.v, vec_add(s,vec_perm(s, s, p16uc_COMPLEX32_REV)))); } -template<> EIGEN_STRONG_INLINE Packet1cd pcplxflip(const Packet1cd& x) +EIGEN_STRONG_INLINE Packet1cd pcplxflip/**/(const Packet1cd& x) { - return Packet1cd(vec_perm(x.v, x.v, p16uc_COMPLEX32_REV)); + return Packet1cd(preverse(Packet2d(x.v))); } EIGEN_STRONG_INLINE void ptranspose(PacketBlock& kernel) diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h index aea3e79dc..f319f4266 100755 --- a/Eigen/src/Core/arch/AltiVec/PacketMath.h +++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h @@ -94,18 +94,21 @@ static Packet16uc p16uc_FORWARD = p16uc_REVERSE32; static Packet16uc p16uc_PSET32_WODD = vec_sld((Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 1), (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 3), 8);//{ 0,1,2,3, 0,1,2,3, 8,9,10,11, 8,9,10,11 }; static Packet16uc p16uc_PSET32_WEVEN = vec_sld((Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 0), (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 2), 8);//{ 4,5,6,7, 4,5,6,7, 12,13,14,15, 12,13,14,15 }; static Packet16uc p16uc_HALF64_0_16 = vec_sld(vec_splat((Packet16uc) vec_abs(p4i_MINUS16), 0), (Packet16uc)p4i_ZERO, 8); //{ 0,0,0,0, 0,0,0,0, 16,16,16,16, 16,16,16,16}; - #endif // _BIG_ENDIAN 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 }; 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 }; - -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 }; -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 }; - 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}; 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}; +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 }; + +#ifdef _BIG_ENDIAN +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 }; +#else +static Packet16uc p16uc_COMPLEX32_REV2 = vec_sld(p16uc_PSET64_HI, p16uc_PSET64_LO, 8); //{ 8,9,10,11, 12,13,14,15, 0,1,2,3, 4,5,6,7 }; +#endif // _BIG_ENDIAN + template<> struct packet_traits : default_packet_traits { typedef Packet4f type; @@ -703,12 +706,11 @@ 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 Packet2l p2l_ZERO = (Packet2l) p4i_ZERO; 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 = { 0.0, 1.0 }; +static Packet2d p2d_ZERO = (Packet2d) p4f_ZERO; +static Packet2d p2d_ZERO_ = { -0.0, -0.0 }; +static Packet2d p2d_COUNTDOWN = (Packet2d) vec_sld((Packet16uc) p2d_ONE, (Packet16uc) p2d_ZERO, 8); static EIGEN_STRONG_INLINE Packet2d vec_splat_dbl(Packet2d& a, int index) { From c0205ca4af402eadcf1a21bc0949afa4c9d9712a Mon Sep 17 00:00:00 2001 From: Konstantinos Margaritis Date: Sun, 21 Sep 2014 08:12:22 +0000 Subject: [PATCH 5/8] VSX supports vec_div, implement where appropriate (float/doubles) --- Eigen/src/Core/arch/AltiVec/PacketMath.h | 22 ++++++---------------- 1 file changed, 6 insertions(+), 16 deletions(-) diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h index f319f4266..3555c521d 100755 --- a/Eigen/src/Core/arch/AltiVec/PacketMath.h +++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h @@ -336,6 +336,7 @@ template<> EIGEN_STRONG_INLINE Packet4i pmul(const Packet4i& a, const */ template<> EIGEN_STRONG_INLINE Packet4f pdiv(const Packet4f& a, const Packet4f& b) { +#ifndef __VSX__ // VSX actually provides a div instruction Packet4f t, y_0, y_1, res; // Altivec does not offer a divide instruction, we have to do a reciprocal approximation @@ -345,8 +346,10 @@ template<> EIGEN_STRONG_INLINE Packet4f pdiv(const Packet4f& a, const t = vec_nmsub(y_0, b, p4f_ONE); y_1 = vec_madd(y_0, t, y_0); - res = vec_madd(a, y_1, p4f_ZERO); - return res; + return vec_madd(a, y_1, p4f_ZERO); +#else + return vec_div(a, b); +#endif } template<> EIGEN_STRONG_INLINE Packet4i pdiv(const Packet4i& /*a*/, const Packet4i& /*b*/) @@ -801,20 +804,7 @@ template<> EIGEN_STRONG_INLINE Packet2d pnegate(const Packet2d& a) { return psub template<> EIGEN_STRONG_INLINE Packet2d pconj(const Packet2d& a) { return a; } template<> EIGEN_STRONG_INLINE Packet2d pmul(const Packet2d& a, const Packet2d& b) { return vec_madd(a,b,p2d_ZERO); } -template<> EIGEN_STRONG_INLINE Packet2d pdiv(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; -} +template<> EIGEN_STRONG_INLINE Packet2d pdiv(const Packet2d& a, const Packet2d& b) { return vec_div(a,b); } // 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); } From 974fe38ca38fb812f3ce4be105cf2ff69f40aef1 Mon Sep 17 00:00:00 2001 From: Konstantinos Margaritis Date: Sun, 21 Sep 2014 11:24:30 +0000 Subject: [PATCH 6/8] prefetch are noops on VSX --- Eigen/src/Core/arch/AltiVec/PacketMath.h | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h index 3555c521d..1274f72dd 100755 --- a/Eigen/src/Core/arch/AltiVec/PacketMath.h +++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h @@ -28,7 +28,7 @@ namespace internal { // NOTE Altivec has 32 registers, but Eigen only accepts a value of 8 or 16 #ifndef EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS -#define EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS 16 +#define EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS 32 #endif typedef __vector float Packet4f; @@ -482,8 +482,16 @@ template<> EIGEN_STRONG_INLINE void pstoreu(float* to, const Packet4f& } #endif -template<> EIGEN_STRONG_INLINE void prefetch(const float* addr) { vec_dstt(addr, DST_CTRL(2,2,32), DST_CHAN); } -template<> EIGEN_STRONG_INLINE void prefetch(const int* addr) { vec_dstt(addr, DST_CTRL(2,2,32), DST_CHAN); } +template<> EIGEN_STRONG_INLINE void prefetch(const float* addr) { +#ifndef __VSX__ + vec_dstt(addr, DST_CTRL(2,2,32), DST_CHAN); +#endif +} +template<> EIGEN_STRONG_INLINE void prefetch(const int* addr) { +#ifndef __VSX__ + vec_dstt(addr, DST_CTRL(2,2,32), DST_CHAN); +#endif +} template<> EIGEN_STRONG_INLINE float pfirst(const Packet4f& a) { float EIGEN_ALIGN16 x[4]; vec_st(a, 0, x); return x[0]; } template<> EIGEN_STRONG_INLINE int pfirst(const Packet4i& a) { int EIGEN_ALIGN16 x[4]; vec_st(a, 0, x); return x[0]; } From 56408504e4e3fa5f9c59d9edac14ca1ba1255e5a Mon Sep 17 00:00:00 2001 From: Konstantinos Margaritis Date: Sun, 21 Sep 2014 13:59:30 +0300 Subject: [PATCH 7/8] fix compile error on big endian altivec --- Eigen/src/Core/arch/AltiVec/Complex.h | 3 ++- Eigen/src/Core/arch/AltiVec/PacketMath.h | 5 +++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/Eigen/src/Core/arch/AltiVec/Complex.h b/Eigen/src/Core/arch/AltiVec/Complex.h index 58f5fbd27..30e2088ef 100644 --- a/Eigen/src/Core/arch/AltiVec/Complex.h +++ b/Eigen/src/Core/arch/AltiVec/Complex.h @@ -242,6 +242,7 @@ EIGEN_STRONG_INLINE void ptranspose(PacketBlock& kernel) } //---------- double ---------- +#ifdef __VSX__ struct Packet1cd { EIGEN_STRONG_INLINE Packet1cd() {} @@ -419,7 +420,7 @@ EIGEN_STRONG_INLINE void ptranspose(PacketBlock& kernel) kernel.packet[1].v = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_TRANSPOSE64_LO); kernel.packet[0].v = tmp; } - +#endif // __VSX__ } // end namespace internal } // end namespace Eigen diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h index 1274f72dd..728f8b4d6 100755 --- a/Eigen/src/Core/arch/AltiVec/PacketMath.h +++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h @@ -76,12 +76,12 @@ static Packet4f p4f_COUNTDOWN = { 0.0, 1.0, 2.0, 3.0 }; static Packet4i p4i_COUNTDOWN = { 0, 1, 2, 3 }; static Packet16uc p16uc_REVERSE32 = { 12,13,14,15, 8,9,10,11, 4,5,6,7, 0,1,2,3 }; -static Packet16uc p16uc_REVERSE64 = { 8,9,10,11, 12,13,14,15, 0,1,2,3, 4,5,6,7 }; static Packet16uc p16uc_DUPLICATE32_HI = { 0,1,2,3, 0,1,2,3, 4,5,6,7, 4,5,6,7 }; // Handle endianness properly while loading constants // Define global static constants: #ifdef _BIG_ENDIAN +static Packet16uc p16uc_FORWARD = vec_lvsl(0, (float*)0); 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 }; 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 }; 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}; @@ -91,6 +91,7 @@ static Packet16uc p16uc_HALF64_0_16 = vec_sld((Packet16uc)p4i_ZERO, vec_splat((P #define _EIGEN_ALIGNED_PTR(x) ((ptrdiff_t)(x) & _EIGEN_MASK_ALIGNMENT) static Packet16uc p16uc_FORWARD = p16uc_REVERSE32; +static Packet16uc p16uc_REVERSE64 = { 8,9,10,11, 12,13,14,15, 0,1,2,3, 4,5,6,7 }; static Packet16uc p16uc_PSET32_WODD = vec_sld((Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 1), (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 3), 8);//{ 0,1,2,3, 0,1,2,3, 8,9,10,11, 8,9,10,11 }; static Packet16uc p16uc_PSET32_WEVEN = vec_sld((Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 0), (Packet16uc) vec_splat((Packet4ui)p16uc_FORWARD, 2), 8);//{ 4,5,6,7, 4,5,6,7, 12,13,14,15, 12,13,14,15 }; static Packet16uc p16uc_HALF64_0_16 = vec_sld(vec_splat((Packet16uc) vec_abs(p4i_MINUS16), 0), (Packet16uc)p4i_ZERO, 8); //{ 0,0,0,0, 0,0,0,0, 16,16,16,16, 16,16,16,16}; @@ -337,7 +338,7 @@ template<> EIGEN_STRONG_INLINE Packet4i pmul(const Packet4i& a, const template<> EIGEN_STRONG_INLINE Packet4f pdiv(const Packet4f& a, const Packet4f& b) { #ifndef __VSX__ // VSX actually provides a div instruction - Packet4f t, y_0, y_1, res; + Packet4f t, y_0, y_1; // Altivec does not offer a divide instruction, we have to do a reciprocal approximation y_0 = vec_re(b); From de38ff24993af77f188ab97d458d2884eeffd1ba Mon Sep 17 00:00:00 2001 From: Konstantinos Margaritis Date: Sun, 21 Sep 2014 11:56:07 +0000 Subject: [PATCH 8/8] prefetch are noops on VSX, actually disable the prefetch trait --- Eigen/src/Core/arch/AltiVec/PacketMath.h | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/Eigen/src/Core/arch/AltiVec/PacketMath.h b/Eigen/src/Core/arch/AltiVec/PacketMath.h index 728f8b4d6..1b86e1227 100755 --- a/Eigen/src/Core/arch/AltiVec/PacketMath.h +++ b/Eigen/src/Core/arch/AltiVec/PacketMath.h @@ -483,16 +483,10 @@ template<> EIGEN_STRONG_INLINE void pstoreu(float* to, const Packet4f& } #endif -template<> EIGEN_STRONG_INLINE void prefetch(const float* addr) { #ifndef __VSX__ - vec_dstt(addr, DST_CTRL(2,2,32), DST_CHAN); +template<> EIGEN_STRONG_INLINE void prefetch(const float* addr) { vec_dstt(addr, DST_CTRL(2,2,32), DST_CHAN); } +template<> EIGEN_STRONG_INLINE void prefetch(const int* addr) { vec_dstt(addr, DST_CTRL(2,2,32), DST_CHAN); } #endif -} -template<> EIGEN_STRONG_INLINE void prefetch(const int* addr) { -#ifndef __VSX__ - vec_dstt(addr, DST_CTRL(2,2,32), DST_CHAN); -#endif -} template<> EIGEN_STRONG_INLINE float pfirst(const Packet4f& a) { float EIGEN_ALIGN16 x[4]; vec_st(a, 0, x); return x[0]; } template<> EIGEN_STRONG_INLINE int pfirst(const Packet4i& a) { int EIGEN_ALIGN16 x[4]; vec_st(a, 0, x); return x[0]; }