packetmath.cpp
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #include "main.h"
00027 
00028 // using namespace Eigen;
00029 
00030 namespace Eigen {
00031 namespace internal {
00032 template<typename T> T negate(const T& x) { return -x; }
00033 }
00034 }
00035 
00036 template<typename Scalar> bool isApproxAbs(const Scalar& a, const Scalar& b, const typename NumTraits<Scalar>::Real& refvalue)
00037 {
00038   return internal::isMuchSmallerThan(a-b, refvalue);
00039 }
00040 
00041 template<typename Scalar> bool areApproxAbs(const Scalar* a, const Scalar* b, int size, const typename NumTraits<Scalar>::Real& refvalue)
00042 {
00043   for (int i=0; i<size; ++i)
00044   {
00045     if (!isApproxAbs(a[i],b[i],refvalue))
00046     {
00047       std::cout << "[" << Map<const Matrix<Scalar,1,Dynamic> >(a,size) << "]" << " != " << Map<const Matrix<Scalar,1,Dynamic> >(b,size) << "\n";
00048       return false;
00049     }
00050   }
00051   return true;
00052 }
00053 
00054 template<typename Scalar> bool areApprox(const Scalar* a, const Scalar* b, int size)
00055 {
00056   for (int i=0; i<size; ++i)
00057   {
00058     if (!internal::isApprox(a[i],b[i]))
00059     {
00060       std::cout << "[" << Map<const Matrix<Scalar,1,Dynamic> >(a,size) << "]" << " != " << Map<const Matrix<Scalar,1,Dynamic> >(b,size) << "\n";
00061       return false;
00062     }
00063   }
00064   return true;
00065 }
00066 
00067 
00068 #define CHECK_CWISE2(REFOP, POP) { \
00069   for (int i=0; i<PacketSize; ++i) \
00070     ref[i] = REFOP(data1[i], data1[i+PacketSize]); \
00071   internal::pstore(data2, POP(internal::pload<Packet>(data1), internal::pload<Packet>(data1+PacketSize))); \
00072   VERIFY(areApprox(ref, data2, PacketSize) && #POP); \
00073 }
00074 
00075 #define CHECK_CWISE1(REFOP, POP) { \
00076   for (int i=0; i<PacketSize; ++i) \
00077     ref[i] = REFOP(data1[i]); \
00078   internal::pstore(data2, POP(internal::pload<Packet>(data1))); \
00079   VERIFY(areApprox(ref, data2, PacketSize) && #POP); \
00080 }
00081 
00082 template<bool Cond,typename Packet>
00083 struct packet_helper
00084 {
00085   template<typename T>
00086   inline Packet load(const T* from) const { return internal::pload<Packet>(from); }
00087 
00088   template<typename T>
00089   inline void store(T* to, const Packet& x) const { internal::pstore(to,x); }
00090 };
00091 
00092 template<typename Packet>
00093 struct packet_helper<false,Packet>
00094 {
00095   template<typename T>
00096   inline T load(const T* from) const { return *from; }
00097 
00098   template<typename T>
00099   inline void store(T* to, const T& x) const { *to = x; }
00100 };
00101 
00102 #define CHECK_CWISE1_IF(COND, REFOP, POP) if(COND) { \
00103   packet_helper<COND,Packet> h; \
00104   for (int i=0; i<PacketSize; ++i) \
00105     ref[i] = REFOP(data1[i]); \
00106   h.store(data2, POP(h.load(data1))); \
00107   VERIFY(areApprox(ref, data2, PacketSize) && #POP); \
00108 }
00109 
00110 #define REF_ADD(a,b) ((a)+(b))
00111 #define REF_SUB(a,b) ((a)-(b))
00112 #define REF_MUL(a,b) ((a)*(b))
00113 #define REF_DIV(a,b) ((a)/(b))
00114 
00115 template<typename Scalar> void packetmath()
00116 {
00117   typedef typename internal::packet_traits<Scalar>::type Packet;
00118   const int PacketSize = internal::packet_traits<Scalar>::size;
00119   typedef typename NumTraits<Scalar>::Real RealScalar;
00120 
00121   const int size = PacketSize*4;
00122   EIGEN_ALIGN16 Scalar data1[internal::packet_traits<Scalar>::size*4];
00123   EIGEN_ALIGN16 Scalar data2[internal::packet_traits<Scalar>::size*4];
00124   EIGEN_ALIGN16 Packet packets[PacketSize*2];
00125   EIGEN_ALIGN16 Scalar ref[internal::packet_traits<Scalar>::size*4];
00126   RealScalar refvalue = 0;
00127   for (int i=0; i<size; ++i)
00128   {
00129     data1[i] = internal::random<Scalar>()/RealScalar(PacketSize);
00130     data2[i] = internal::random<Scalar>()/RealScalar(PacketSize);
00131     refvalue = std::max(refvalue,internal::abs(data1[i]));
00132   }
00133 
00134   internal::pstore(data2, internal::pload<Packet>(data1));
00135   VERIFY(areApprox(data1, data2, PacketSize) && "aligned load/store");
00136 
00137   for (int offset=0; offset<PacketSize; ++offset)
00138   {
00139     internal::pstore(data2, internal::ploadu<Packet>(data1+offset));
00140     VERIFY(areApprox(data1+offset, data2, PacketSize) && "internal::ploadu");
00141   }
00142 
00143   for (int offset=0; offset<PacketSize; ++offset)
00144   {
00145     internal::pstoreu(data2+offset, internal::pload<Packet>(data1));
00146     VERIFY(areApprox(data1, data2+offset, PacketSize) && "internal::pstoreu");
00147   }
00148 
00149   for (int offset=0; offset<PacketSize; ++offset)
00150   {
00151     packets[0] = internal::pload<Packet>(data1);
00152     packets[1] = internal::pload<Packet>(data1+PacketSize);
00153          if (offset==0) internal::palign<0>(packets[0], packets[1]);
00154     else if (offset==1) internal::palign<1>(packets[0], packets[1]);
00155     else if (offset==2) internal::palign<2>(packets[0], packets[1]);
00156     else if (offset==3) internal::palign<3>(packets[0], packets[1]);
00157     internal::pstore(data2, packets[0]);
00158 
00159     for (int i=0; i<PacketSize; ++i)
00160       ref[i] = data1[i+offset];
00161 
00162     typedef Matrix<Scalar, PacketSize, 1> Vector;
00163     VERIFY(areApprox(ref, data2, PacketSize) && "internal::palign");
00164   }
00165 
00166   CHECK_CWISE2(REF_ADD,  internal::padd);
00167   CHECK_CWISE2(REF_SUB,  internal::psub);
00168   CHECK_CWISE2(REF_MUL,  internal::pmul);
00169   #ifndef EIGEN_VECTORIZE_ALTIVEC
00170   if (!internal::is_same<Scalar,int>::value)
00171     CHECK_CWISE2(REF_DIV,  internal::pdiv);
00172   #endif
00173   CHECK_CWISE1(internal::negate, internal::pnegate);
00174   CHECK_CWISE1(internal::conj, internal::pconj);
00175 
00176   for(int offset=0;offset<3;++offset)
00177   {
00178     for (int i=0; i<PacketSize; ++i)
00179       ref[i] = data1[offset];
00180     internal::pstore(data2, internal::pset1<Packet>(data1[offset]));
00181     VERIFY(areApprox(ref, data2, PacketSize) && "internal::pset1");
00182   }
00183   
00184   VERIFY(internal::isApprox(data1[0], internal::pfirst(internal::pload<Packet>(data1))) && "internal::pfirst");
00185   
00186   if(PacketSize>1)
00187   {
00188     for(int offset=0;offset<4;++offset)
00189     {
00190       for(int i=0;i<PacketSize/2;++i)
00191         ref[2*i+0] = ref[2*i+1] = data1[offset+i];
00192       internal::pstore(data2,internal::ploaddup<Packet>(data1+offset));
00193       VERIFY(areApprox(ref, data2, PacketSize) && "ploaddup");
00194     }
00195   }
00196 
00197   ref[0] = 0;
00198   for (int i=0; i<PacketSize; ++i)
00199     ref[0] += data1[i];
00200   VERIFY(isApproxAbs(ref[0], internal::predux(internal::pload<Packet>(data1)), refvalue) && "internal::predux");
00201 
00202   ref[0] = 1;
00203   for (int i=0; i<PacketSize; ++i)
00204     ref[0] *= data1[i];
00205   VERIFY(internal::isApprox(ref[0], internal::predux_mul(internal::pload<Packet>(data1))) && "internal::predux_mul");
00206 
00207   for (int j=0; j<PacketSize; ++j)
00208   {
00209     ref[j] = 0;
00210     for (int i=0; i<PacketSize; ++i)
00211       ref[j] += data1[i+j*PacketSize];
00212     packets[j] = internal::pload<Packet>(data1+j*PacketSize);
00213   }
00214   internal::pstore(data2, internal::preduxp(packets));
00215   VERIFY(areApproxAbs(ref, data2, PacketSize, refvalue) && "internal::preduxp");
00216 
00217   for (int i=0; i<PacketSize; ++i)
00218     ref[i] = data1[PacketSize-i-1];
00219   internal::pstore(data2, internal::preverse(internal::pload<Packet>(data1)));
00220   VERIFY(areApprox(ref, data2, PacketSize) && "internal::preverse");
00221 }
00222 
00223 template<typename Scalar> void packetmath_real()
00224 {
00225   typedef typename internal::packet_traits<Scalar>::type Packet;
00226   const int PacketSize = internal::packet_traits<Scalar>::size;
00227 
00228   const int size = PacketSize*4;
00229   EIGEN_ALIGN16 Scalar data1[internal::packet_traits<Scalar>::size*4];
00230   EIGEN_ALIGN16 Scalar data2[internal::packet_traits<Scalar>::size*4];
00231   EIGEN_ALIGN16 Scalar ref[internal::packet_traits<Scalar>::size*4];
00232 
00233   for (int i=0; i<size; ++i)
00234   {
00235     data1[i] = internal::random<Scalar>(-1e3,1e3);
00236     data2[i] = internal::random<Scalar>(-1e3,1e3);
00237   }
00238   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasSin, internal::sin, internal::psin);
00239   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasCos, internal::cos, internal::pcos);
00240   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasTan, internal::tan, internal::ptan);
00241   
00242   for (int i=0; i<size; ++i)
00243   {
00244     data1[i] = internal::random<Scalar>(-1,1);
00245     data2[i] = internal::random<Scalar>(-1,1);
00246   }
00247   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasASin, internal::asin, internal::pasin);
00248   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasACos, internal::acos, internal::pacos);
00249 
00250   for (int i=0; i<size; ++i)
00251   {
00252     data1[i] = internal::random<Scalar>(-87,88);
00253     data2[i] = internal::random<Scalar>(-87,88);
00254   }
00255   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasExp, internal::exp, internal::pexp);
00256 
00257   for (int i=0; i<size; ++i)
00258   {
00259     data1[i] = internal::random<Scalar>(0,1e6);
00260     data2[i] = internal::random<Scalar>(0,1e6);
00261   }
00262   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLog, internal::log, internal::plog);
00263   CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasSqrt, internal::sqrt, internal::psqrt);
00264 
00265   ref[0] = data1[0];
00266   for (int i=0; i<PacketSize; ++i)
00267     ref[0] = std::min(ref[0],data1[i]);
00268   VERIFY(internal::isApprox(ref[0], internal::predux_min(internal::pload<Packet>(data1))) && "internal::predux_min");
00269 
00270   CHECK_CWISE2(std::min, internal::pmin);
00271   CHECK_CWISE2(std::max, internal::pmax);
00272   CHECK_CWISE1(internal::abs, internal::pabs);
00273 
00274   ref[0] = data1[0];
00275   for (int i=0; i<PacketSize; ++i)
00276     ref[0] = std::max(ref[0],data1[i]);
00277   VERIFY(internal::isApprox(ref[0], internal::predux_max(internal::pload<Packet>(data1))) && "internal::predux_max");
00278   
00279   for (int i=0; i<PacketSize; ++i)
00280     ref[i] = data1[0]+Scalar(i);
00281   internal::pstore(data2, internal::plset(data1[0]));
00282   VERIFY(areApprox(ref, data2, PacketSize) && "internal::plset");
00283 }
00284 
00285 template<typename Scalar,bool ConjLhs,bool ConjRhs> void test_conj_helper(Scalar* data1, Scalar* data2, Scalar* ref, Scalar* pval)
00286 {
00287   typedef typename internal::packet_traits<Scalar>::type Packet;
00288   const int PacketSize = internal::packet_traits<Scalar>::size;
00289   
00290   internal::conj_if<ConjLhs> cj0;
00291   internal::conj_if<ConjRhs> cj1;
00292   internal::conj_helper<Scalar,Scalar,ConjLhs,ConjRhs> cj;
00293   internal::conj_helper<Packet,Packet,ConjLhs,ConjRhs> pcj;
00294   
00295   for(int i=0;i<PacketSize;++i)
00296   {
00297     ref[i] = cj0(data1[i]) * cj1(data2[i]);
00298     VERIFY(internal::isApprox(ref[i], cj.pmul(data1[i],data2[i])) && "conj_helper pmul");
00299   }
00300   internal::pstore(pval,pcj.pmul(internal::pload<Packet>(data1),internal::pload<Packet>(data2)));
00301   VERIFY(areApprox(ref, pval, PacketSize) && "conj_helper pmul");
00302   
00303   for(int i=0;i<PacketSize;++i)
00304   {
00305     Scalar tmp = ref[i];
00306     ref[i] += cj0(data1[i]) * cj1(data2[i]);
00307     VERIFY(internal::isApprox(ref[i], cj.pmadd(data1[i],data2[i],tmp)) && "conj_helper pmadd");
00308   }
00309   internal::pstore(pval,pcj.pmadd(internal::pload<Packet>(data1),internal::pload<Packet>(data2),internal::pload<Packet>(pval)));
00310   VERIFY(areApprox(ref, pval, PacketSize) && "conj_helper pmadd");
00311 }
00312 
00313 template<typename Scalar> void packetmath_complex()
00314 {
00315   typedef typename internal::packet_traits<Scalar>::type Packet;
00316   const int PacketSize = internal::packet_traits<Scalar>::size;
00317 
00318   const int size = PacketSize*4;
00319   EIGEN_ALIGN16 Scalar data1[PacketSize*4];
00320   EIGEN_ALIGN16 Scalar data2[PacketSize*4];
00321   EIGEN_ALIGN16 Scalar ref[PacketSize*4];
00322   EIGEN_ALIGN16 Scalar pval[PacketSize*4];
00323 
00324   for (int i=0; i<size; ++i)
00325   {
00326     data1[i] = internal::random<Scalar>() * Scalar(1e2);
00327     data2[i] = internal::random<Scalar>() * Scalar(1e2);
00328   }
00329   
00330   test_conj_helper<Scalar,false,false> (data1,data2,ref,pval);
00331   test_conj_helper<Scalar,false,true>  (data1,data2,ref,pval);
00332   test_conj_helper<Scalar,true,false>  (data1,data2,ref,pval);
00333   test_conj_helper<Scalar,true,true>   (data1,data2,ref,pval);
00334   
00335   {
00336     for(int i=0;i<PacketSize;++i)
00337       ref[i] = Scalar(std::imag(data1[i]),std::real(data1[i]));
00338     internal::pstore(pval,internal::pcplxflip(internal::pload<Packet>(data1)));
00339     VERIFY(areApprox(ref, pval, PacketSize) && "pcplxflip");
00340   }
00341   
00342   
00343 }
00344 
00345 void test_packetmath()
00346 {
00347   for(int i = 0; i < g_repeat; i++) {
00348     CALL_SUBTEST_1( packetmath<float>() );
00349     CALL_SUBTEST_2( packetmath<double>() );
00350     CALL_SUBTEST_3( packetmath<int>() );
00351     CALL_SUBTEST_1( packetmath<std::complex<float> >() );
00352     CALL_SUBTEST_2( packetmath<std::complex<double> >() );
00353 
00354     CALL_SUBTEST_1( packetmath_real<float>() );
00355     CALL_SUBTEST_2( packetmath_real<double>() );
00356 
00357     CALL_SUBTEST_1( packetmath_complex<std::complex<float> >() );
00358     CALL_SUBTEST_2( packetmath_complex<std::complex<double> >() );
00359   }
00360 }


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:08