Eigen  3.3.3
Complex.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2016 Konstantinos Margaritis <markos@freevec.org>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_COMPLEX32_ALTIVEC_H
00012 #define EIGEN_COMPLEX32_ALTIVEC_H
00013 
00014 namespace Eigen {
00015 
00016 namespace internal {
00017 
00018 static Packet2ul  p2ul_CONJ_XOR1 = (Packet2ul) vec_sld((Packet4ui) p2d_ZERO_, (Packet4ui) p2l_ZERO, 8);//{ 0x8000000000000000, 0x0000000000000000 };
00019 static Packet2ul  p2ul_CONJ_XOR2 = (Packet2ul) vec_sld((Packet4ui) p2l_ZERO,  (Packet4ui) p2d_ZERO_, 8);//{ 0x8000000000000000, 0x0000000000000000 };
00020 
00021 struct Packet1cd
00022 {
00023   EIGEN_STRONG_INLINE Packet1cd() {}
00024   EIGEN_STRONG_INLINE explicit Packet1cd(const Packet2d& a) : v(a) {}
00025   Packet2d v;
00026 };
00027 
00028 struct Packet2cf
00029 {
00030   EIGEN_STRONG_INLINE Packet2cf() {}
00031   EIGEN_STRONG_INLINE explicit Packet2cf(const Packet4f& a) : v(a) {}
00032   union {
00033     Packet4f v;
00034     Packet1cd cd[2];
00035   };
00036 };
00037 
00038 template<> struct packet_traits<std::complex<float> >  : default_packet_traits
00039 {
00040   typedef Packet2cf type;
00041   typedef Packet2cf half;
00042   enum {
00043     Vectorizable = 1,
00044     AlignedOnScalar = 1,
00045     size = 2,
00046     HasHalfPacket = 0,
00047 
00048     HasAdd    = 1,
00049     HasSub    = 1,
00050     HasMul    = 1,
00051     HasDiv    = 1,
00052     HasNegate = 1,
00053     HasAbs    = 0,
00054     HasAbs2   = 0,
00055     HasMin    = 0,
00056     HasMax    = 0,
00057     HasBlend  = 1,
00058     HasSetLinear = 0
00059   };
00060 };
00061 
00062 
00063 template<> struct packet_traits<std::complex<double> >  : default_packet_traits
00064 {
00065   typedef Packet1cd type;
00066   typedef Packet1cd half;
00067   enum {
00068     Vectorizable = 1,
00069     AlignedOnScalar = 1,
00070     size = 1,
00071     HasHalfPacket = 0,
00072 
00073     HasAdd    = 1,
00074     HasSub    = 1,
00075     HasMul    = 1,
00076     HasDiv    = 1,
00077     HasNegate = 1,
00078     HasAbs    = 0,
00079     HasAbs2   = 0,
00080     HasMin    = 0,
00081     HasMax    = 0,
00082     HasSetLinear = 0
00083   };
00084 };
00085 
00086 template<> struct unpacket_traits<Packet2cf> { typedef std::complex<float>  type; enum {size=2, alignment=Aligned16}; typedef Packet2cf half; };
00087 template<> struct unpacket_traits<Packet1cd> { typedef std::complex<double> type; enum {size=1, alignment=Aligned16}; typedef Packet1cd half; };
00088 
00089 /* Forward declaration */
00090 EIGEN_STRONG_INLINE void ptranspose(PacketBlock<Packet2cf,2>& kernel);
00091 
00092 template<> EIGEN_STRONG_INLINE Packet2cf pload <Packet2cf>(const std::complex<float>* from)  { EIGEN_DEBUG_ALIGNED_LOAD return Packet2cf(pload<Packet4f>((const float*)from)); }
00093 template<> EIGEN_STRONG_INLINE Packet1cd pload <Packet1cd>(const std::complex<double>* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet1cd(pload<Packet2d>((const double*)from)); }
00094 template<> EIGEN_STRONG_INLINE Packet2cf ploadu<Packet2cf>(const std::complex<float>* from)  { EIGEN_DEBUG_UNALIGNED_LOAD return Packet2cf(ploadu<Packet4f>((const float*)from)); }
00095 template<> EIGEN_STRONG_INLINE Packet1cd ploadu<Packet1cd>(const std::complex<double>* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet1cd(ploadu<Packet2d>((const double*)from)); }
00096 template<> EIGEN_STRONG_INLINE void pstore <std::complex<float> >(std::complex<float> *     to, const Packet2cf& from) { EIGEN_DEBUG_ALIGNED_STORE pstore((float*)to, from.v); }
00097 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); }
00098 template<> EIGEN_STRONG_INLINE void pstoreu<std::complex<float> >(std::complex<float> *     to, const Packet2cf& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((float*)to, from.v); }
00099 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); }
00100 
00101 template<> EIGEN_STRONG_INLINE Packet1cd pset1<Packet1cd>(const std::complex<double>&  from)
00102 { /* here we really have to use unaligned loads :( */ return ploadu<Packet1cd>(&from); }
00103 
00104 template<> EIGEN_STRONG_INLINE Packet2cf pset1<Packet2cf>(const std::complex<float>&  from)
00105 {
00106   Packet2cf res;
00107   res.cd[0] = Packet1cd(vec_ld2f((const float *)&from));
00108   res.cd[1] = res.cd[0];
00109   return res;
00110 }
00111 template<> EIGEN_DEVICE_FUNC inline Packet2cf pgather<std::complex<float>, Packet2cf>(const std::complex<float>* from, Index stride)
00112 {
00113   std::complex<float> EIGEN_ALIGN16 af[2];
00114   af[0] = from[0*stride];
00115   af[1] = from[1*stride];
00116   return pload<Packet2cf>(af);
00117 }
00118 template<> EIGEN_DEVICE_FUNC inline Packet1cd pgather<std::complex<double>, Packet1cd>(const std::complex<double>* from, Index stride EIGEN_UNUSED)
00119 {
00120   return pload<Packet1cd>(from);
00121 }
00122 template<> EIGEN_DEVICE_FUNC inline void pscatter<std::complex<float>, Packet2cf>(std::complex<float>* to, const Packet2cf& from, Index stride)
00123 {
00124   std::complex<float> EIGEN_ALIGN16 af[2];
00125   pstore<std::complex<float> >((std::complex<float> *) af, from);
00126   to[0*stride] = af[0];
00127   to[1*stride] = af[1];
00128 }
00129 template<> EIGEN_DEVICE_FUNC inline void pscatter<std::complex<double>, Packet1cd>(std::complex<double>* to, const Packet1cd& from, Index stride EIGEN_UNUSED)
00130 {
00131   pstore<std::complex<double> >(to, from);
00132 }
00133 
00134 template<> EIGEN_STRONG_INLINE Packet2cf padd<Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(padd<Packet4f>(a.v, b.v)); }
00135 template<> EIGEN_STRONG_INLINE Packet1cd padd<Packet1cd>(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(a.v + b.v); }
00136 template<> EIGEN_STRONG_INLINE Packet2cf psub<Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(psub<Packet4f>(a.v, b.v)); }
00137 template<> EIGEN_STRONG_INLINE Packet1cd psub<Packet1cd>(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(a.v - b.v); }
00138 template<> EIGEN_STRONG_INLINE Packet1cd pnegate(const Packet1cd& a) { return Packet1cd(pnegate(Packet2d(a.v))); }
00139 template<> EIGEN_STRONG_INLINE Packet2cf pnegate(const Packet2cf& a) { return Packet2cf(pnegate(Packet4f(a.v))); }
00140 template<> EIGEN_STRONG_INLINE Packet1cd pconj(const Packet1cd& a) { return Packet1cd((Packet2d)vec_xor((Packet2d)a.v, (Packet2d)p2ul_CONJ_XOR2)); }
00141 template<> EIGEN_STRONG_INLINE Packet2cf pconj(const Packet2cf& a)
00142 {
00143   Packet2cf res;
00144   res.v.v4f[0] = pconj(Packet1cd(reinterpret_cast<Packet2d>(a.v.v4f[0]))).v;
00145   res.v.v4f[1] = pconj(Packet1cd(reinterpret_cast<Packet2d>(a.v.v4f[1]))).v;
00146   return res;
00147 }
00148 
00149 template<> EIGEN_STRONG_INLINE Packet1cd pmul<Packet1cd>(const Packet1cd& a, const Packet1cd& b)
00150 {
00151   Packet2d a_re, a_im, v1, v2;
00152 
00153   // Permute and multiply the real parts of a and b
00154   a_re = vec_perm(a.v, a.v, p16uc_PSET64_HI);
00155   // Get the imaginary parts of a
00156   a_im = vec_perm(a.v, a.v, p16uc_PSET64_LO);
00157   // multiply a_re * b
00158   v1 = vec_madd(a_re, b.v, p2d_ZERO);
00159   // multiply a_im * b and get the conjugate result
00160   v2 = vec_madd(a_im, b.v, p2d_ZERO);
00161   v2 = (Packet2d) vec_sld((Packet4ui)v2, (Packet4ui)v2, 8);
00162   v2 = (Packet2d) vec_xor((Packet2d)v2, (Packet2d) p2ul_CONJ_XOR1);
00163 
00164   return Packet1cd(v1 + v2);
00165 }
00166 template<> EIGEN_STRONG_INLINE Packet2cf pmul<Packet2cf>(const Packet2cf& a, const Packet2cf& b)
00167 {
00168   Packet2cf res;
00169   res.v.v4f[0] = pmul(Packet1cd(reinterpret_cast<Packet2d>(a.v.v4f[0])), Packet1cd(reinterpret_cast<Packet2d>(b.v.v4f[0]))).v;
00170   res.v.v4f[1] = pmul(Packet1cd(reinterpret_cast<Packet2d>(a.v.v4f[1])), Packet1cd(reinterpret_cast<Packet2d>(b.v.v4f[1]))).v;
00171   return res;
00172 }
00173 
00174 template<> EIGEN_STRONG_INLINE Packet1cd pand   <Packet1cd>(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_and(a.v,b.v)); }
00175 template<> EIGEN_STRONG_INLINE Packet2cf pand   <Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(pand<Packet4f>(a.v,b.v)); }
00176 template<> EIGEN_STRONG_INLINE Packet1cd por    <Packet1cd>(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_or(a.v,b.v)); }
00177 template<> EIGEN_STRONG_INLINE Packet2cf por    <Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(por<Packet4f>(a.v,b.v)); }
00178 template<> EIGEN_STRONG_INLINE Packet1cd pxor   <Packet1cd>(const Packet1cd& a, const Packet1cd& b) { return Packet1cd(vec_xor(a.v,b.v)); }
00179 template<> EIGEN_STRONG_INLINE Packet2cf pxor   <Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(pxor<Packet4f>(a.v,b.v)); }
00180 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))); }
00181 template<> EIGEN_STRONG_INLINE Packet2cf pandnot<Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(pandnot<Packet4f>(a.v,b.v)); }
00182 
00183 template<> EIGEN_STRONG_INLINE Packet1cd ploaddup<Packet1cd>(const std::complex<double>*     from) {  return pset1<Packet1cd>(*from); }
00184 template<> EIGEN_STRONG_INLINE Packet2cf ploaddup<Packet2cf>(const std::complex<float>*      from) {  return pset1<Packet2cf>(*from); }
00185 
00186 template<> EIGEN_STRONG_INLINE void prefetch<std::complex<float> >(const std::complex<float> *     addr) { EIGEN_ZVECTOR_PREFETCH(addr); }
00187 template<> EIGEN_STRONG_INLINE void prefetch<std::complex<double> >(const std::complex<double> *   addr) { EIGEN_ZVECTOR_PREFETCH(addr); }
00188 
00189 template<> EIGEN_STRONG_INLINE std::complex<double>  pfirst<Packet1cd>(const Packet1cd& a)
00190 {
00191   std::complex<double> EIGEN_ALIGN16 res;
00192   pstore<std::complex<double> >(&res, a);
00193 
00194   return res;
00195 }
00196 template<> EIGEN_STRONG_INLINE std::complex<float>  pfirst<Packet2cf>(const Packet2cf& a)
00197 {
00198   std::complex<float> EIGEN_ALIGN16 res[2];
00199   pstore<std::complex<float> >(res, a);
00200 
00201   return res[0];
00202 }
00203 
00204 template<> EIGEN_STRONG_INLINE Packet1cd preverse(const Packet1cd& a) { return a; }
00205 template<> EIGEN_STRONG_INLINE Packet2cf preverse(const Packet2cf& a)
00206 {
00207   Packet2cf res;
00208   res.cd[0] = a.cd[1];
00209   res.cd[1] = a.cd[0];
00210   return res;
00211 }
00212 
00213 template<> EIGEN_STRONG_INLINE std::complex<double> predux<Packet1cd>(const Packet1cd& a)
00214 {
00215   return pfirst(a);
00216 }
00217 template<> EIGEN_STRONG_INLINE std::complex<float> predux<Packet2cf>(const Packet2cf& a)
00218 {
00219   std::complex<float> res;
00220   Packet1cd b = padd<Packet1cd>(a.cd[0], a.cd[1]);
00221   vec_st2f(b.v, (float*)&res);
00222   return res;
00223 }
00224 
00225 template<> EIGEN_STRONG_INLINE Packet1cd preduxp<Packet1cd>(const Packet1cd* vecs)
00226 {
00227   return vecs[0];
00228 }
00229 template<> EIGEN_STRONG_INLINE Packet2cf preduxp<Packet2cf>(const Packet2cf* vecs)
00230 {
00231   PacketBlock<Packet2cf,2> transpose;
00232   transpose.packet[0] = vecs[0];
00233   transpose.packet[1] = vecs[1];
00234   ptranspose(transpose);
00235 
00236   return padd<Packet2cf>(transpose.packet[0], transpose.packet[1]);
00237 } 
00238 
00239 template<> EIGEN_STRONG_INLINE std::complex<double> predux_mul<Packet1cd>(const Packet1cd& a)
00240 {
00241   return pfirst(a);
00242 }
00243 template<> EIGEN_STRONG_INLINE std::complex<float> predux_mul<Packet2cf>(const Packet2cf& a)
00244 {
00245   std::complex<float> res;
00246   Packet1cd b = pmul<Packet1cd>(a.cd[0], a.cd[1]);
00247   vec_st2f(b.v, (float*)&res);
00248   return res;
00249 }
00250 
00251 template<int Offset>
00252 struct palign_impl<Offset,Packet1cd>
00253 {
00254   static EIGEN_STRONG_INLINE void run(Packet1cd& /*first*/, const Packet1cd& /*second*/)
00255   {
00256     // FIXME is it sure we never have to align a Packet1cd?
00257     // Even though a std::complex<double> has 16 bytes, it is not necessarily aligned on a 16 bytes boundary...
00258   }
00259 };
00260 
00261 template<int Offset>
00262 struct palign_impl<Offset,Packet2cf>
00263 {
00264   static EIGEN_STRONG_INLINE void run(Packet2cf& first, const Packet2cf& second)
00265   {
00266     if (Offset == 1) {
00267       first.cd[0] = first.cd[1];
00268       first.cd[1] = second.cd[0];
00269     }
00270   }
00271 };
00272 
00273 template<> struct conj_helper<Packet1cd, Packet1cd, false,true>
00274 {
00275   EIGEN_STRONG_INLINE Packet1cd pmadd(const Packet1cd& x, const Packet1cd& y, const Packet1cd& c) const
00276   { return padd(pmul(x,y),c); }
00277 
00278   EIGEN_STRONG_INLINE Packet1cd pmul(const Packet1cd& a, const Packet1cd& b) const
00279   {
00280     return internal::pmul(a, pconj(b));
00281   }
00282 };
00283 
00284 template<> struct conj_helper<Packet1cd, Packet1cd, true,false>
00285 {
00286   EIGEN_STRONG_INLINE Packet1cd pmadd(const Packet1cd& x, const Packet1cd& y, const Packet1cd& c) const
00287   { return padd(pmul(x,y),c); }
00288 
00289   EIGEN_STRONG_INLINE Packet1cd pmul(const Packet1cd& a, const Packet1cd& b) const
00290   {
00291     return internal::pmul(pconj(a), b);
00292   }
00293 };
00294 
00295 template<> struct conj_helper<Packet1cd, Packet1cd, true,true>
00296 {
00297   EIGEN_STRONG_INLINE Packet1cd pmadd(const Packet1cd& x, const Packet1cd& y, const Packet1cd& c) const
00298   { return padd(pmul(x,y),c); }
00299 
00300   EIGEN_STRONG_INLINE Packet1cd pmul(const Packet1cd& a, const Packet1cd& b) const
00301   {
00302     return pconj(internal::pmul(a, b));
00303   }
00304 };
00305 
00306 template<> struct conj_helper<Packet2cf, Packet2cf, false,true>
00307 {
00308   EIGEN_STRONG_INLINE Packet2cf pmadd(const Packet2cf& x, const Packet2cf& y, const Packet2cf& c) const
00309   { return padd(pmul(x,y),c); }
00310 
00311   EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, const Packet2cf& b) const
00312   {
00313     return internal::pmul(a, pconj(b));
00314   }
00315 };
00316 
00317 template<> struct conj_helper<Packet2cf, Packet2cf, true,false>
00318 {
00319   EIGEN_STRONG_INLINE Packet2cf pmadd(const Packet2cf& x, const Packet2cf& y, const Packet2cf& c) const
00320   { return padd(pmul(x,y),c); }
00321 
00322   EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, const Packet2cf& b) const
00323   {
00324     return internal::pmul(pconj(a), b);
00325   }
00326 };
00327 
00328 template<> struct conj_helper<Packet2cf, Packet2cf, true,true>
00329 {
00330   EIGEN_STRONG_INLINE Packet2cf pmadd(const Packet2cf& x, const Packet2cf& y, const Packet2cf& c) const
00331   { return padd(pmul(x,y),c); }
00332 
00333   EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, const Packet2cf& b) const
00334   {
00335     return pconj(internal::pmul(a, b));
00336   }
00337 };
00338 
00339 template<> EIGEN_STRONG_INLINE Packet1cd pdiv<Packet1cd>(const Packet1cd& a, const Packet1cd& b)
00340 {
00341   // TODO optimize it for AltiVec
00342   Packet1cd res = conj_helper<Packet1cd,Packet1cd,false,true>().pmul(a,b);
00343   Packet2d s = vec_madd(b.v, b.v, p2d_ZERO_);
00344   return Packet1cd(pdiv(res.v, s + vec_perm(s, s, p16uc_REVERSE64)));
00345 }
00346 
00347 template<> EIGEN_STRONG_INLINE Packet2cf pdiv<Packet2cf>(const Packet2cf& a, const Packet2cf& b)
00348 {
00349   // TODO optimize it for AltiVec
00350   Packet2cf res;
00351   res.cd[0] = pdiv<Packet1cd>(a.cd[0], b.cd[0]);
00352   res.cd[1] = pdiv<Packet1cd>(a.cd[1], b.cd[1]);
00353   return res;
00354 }
00355 
00356 EIGEN_STRONG_INLINE Packet1cd pcplxflip/*<Packet1cd>*/(const Packet1cd& x)
00357 {
00358   return Packet1cd(preverse(Packet2d(x.v)));
00359 }
00360 
00361 EIGEN_STRONG_INLINE Packet2cf pcplxflip/*<Packet2cf>*/(const Packet2cf& x)
00362 {
00363   Packet2cf res;
00364   res.cd[0] = pcplxflip(x.cd[0]);
00365   res.cd[1] = pcplxflip(x.cd[1]);
00366   return res;
00367 }
00368 
00369 EIGEN_STRONG_INLINE void ptranspose(PacketBlock<Packet1cd,2>& kernel)
00370 {
00371   Packet2d tmp = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_TRANSPOSE64_HI);
00372   kernel.packet[1].v = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_TRANSPOSE64_LO);
00373   kernel.packet[0].v = tmp;
00374 }
00375 
00376 EIGEN_STRONG_INLINE void ptranspose(PacketBlock<Packet2cf,2>& kernel)
00377 {
00378   Packet1cd tmp = kernel.packet[0].cd[1];
00379   kernel.packet[0].cd[1] = kernel.packet[1].cd[0];
00380   kernel.packet[1].cd[0] = tmp;
00381 }
00382 
00383 template<> EIGEN_STRONG_INLINE Packet2cf pblend(const Selector<2>& ifPacket, const Packet2cf& thenPacket, const Packet2cf& elsePacket) {
00384   Packet2cf result;
00385   const Selector<4> ifPacket4 = { ifPacket.select[0], ifPacket.select[0], ifPacket.select[1], ifPacket.select[1] };
00386   result.v = pblend<Packet4f>(ifPacket4, thenPacket.v, elsePacket.v);
00387   return result;
00388 }
00389 
00390 } // end namespace internal
00391 
00392 } // end namespace Eigen
00393 
00394 #endif // EIGEN_COMPLEX32_ALTIVEC_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends