ground/gcs/src/libs/eigen/bench/bench_norm.cpp

   1 #include <typeinfo>
   2 #include <iostream>
   3 #include <Eigen/Core>
   4 #include "BenchTimer.h"
   5 using namespace Eigen;
   6 using namespace std;
   7
   8 template<typename T>
   9 EIGEN_DONT_INLINE typename T::Scalar sqsumNorm(const T& v)
  10 {
  11   return v.norm();
  12 }
  13
  14 template<typename T>
  15 EIGEN_DONT_INLINE typename T::Scalar hypotNorm(const T& v)
  16 {
  17   return v.hypotNorm();
  18 }
  19
  20 template<typename T>
  21 EIGEN_DONT_INLINE typename T::Scalar blueNorm(const T& v)
  22 {
  23   return v.blueNorm();
  24 }
  25
  26 template<typename T>
  27 EIGEN_DONT_INLINE typename T::Scalar lapackNorm(T& v)
  28 {
  29   typedef typename T::Scalar Scalar;
  30   int n = v.size();
  31   Scalar scale = 0;
  32   Scalar ssq = 1;
  33   for (int i=0;i<n;++i)
  34   {
  35     Scalar ax = internal::abs(v.coeff(i));
  36     if (scale >= ax)
  37     {
  38       ssq += internal::abs2(ax/scale);
  39     }
  40     else
  41     {
  42       ssq = Scalar(1) + ssq * internal::abs2(scale/ax);
  43       scale = ax;
  44     }
  45   }
  46   return scale * internal::sqrt(ssq);
  47 }
  48
  49 template<typename T>
  50 EIGEN_DONT_INLINE typename T::Scalar twopassNorm(T& v)
  51 {
  52   typedef typename T::Scalar Scalar;
  53   Scalar s = v.cwise().abs().maxCoeff();
  54   return s*(v/s).norm();
  55 }
  56
  57 template<typename T>
  58 EIGEN_DONT_INLINE typename T::Scalar bl2passNorm(T& v)
  59 {
  60   return v.stableNorm();
  61 }
  62
  63 template<typename T>
  64 EIGEN_DONT_INLINE typename T::Scalar divacNorm(T& v)
  65 {
  66   int n =v.size() / 2;
  67   for (int i=0;i<n;++i)
  68     v(i) = v(2*i)*v(2*i) + v(2*i+1)*v(2*i+1);
  69   n = n/2;
  70   while (n>0)
  71   {
  72     for (int i=0;i<n;++i)
  73       v(i) = v(2*i) + v(2*i+1);
  74     n = n/2;
  75   }
  76   return internal::sqrt(v(0));
  77 }
  78
  79 #ifdef EIGEN_VECTORIZE
  80 Packet4f internal::plt(const Packet4f& a, Packet4f& b) { return _mm_cmplt_ps(a,b); }
  81 Packet2d internal::plt(const Packet2d& a, Packet2d& b) { return _mm_cmplt_pd(a,b); }
  82
  83 Packet4f internal::pandnot(const Packet4f& a, Packet4f& b) { return _mm_andnot_ps(a,b); }
  84 Packet2d internal::pandnot(const Packet2d& a, Packet2d& b) { return _mm_andnot_pd(a,b); }
  85 #endif
  86
  87 template<typename T>
  88 EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v)
  89 {
  90   #ifndef EIGEN_VECTORIZE
  91   return v.blueNorm();
  92   #else
  93   typedef typename T::Scalar Scalar;
  94
  95   static int nmax = 0;
  96   static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr;
  97   int n;
  98
  99   if(nmax <= 0)
 100   {
 101     int nbig, ibeta, it, iemin, iemax, iexp;
 102     Scalar abig, eps;
 103
 104     nbig  = std::numeric_limits<int>::max();            // largest integer
 105     ibeta = std::numeric_limits<Scalar>::radix; //NumTraits<Scalar>::Base;                    // base for floating-point numbers
 106     it    = std::numeric_limits<Scalar>::digits; //NumTraits<Scalar>::Mantissa;                // number of base-beta digits in mantissa
 107     iemin = std::numeric_limits<Scalar>::min_exponent;  // minimum exponent
 108     iemax = std::numeric_limits<Scalar>::max_exponent;  // maximum exponent
 109     rbig  = std::numeric_limits<Scalar>::max();         // largest floating-point number
 110
 111     // Check the basic machine-dependent constants.
 112     if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5)
 113       || (it<=4 && ibeta <= 3 ) || it<2)
 114     {
 115       eigen_assert(false && "the algorithm cannot be guaranteed on this computer");
 116     }
 117     iexp  = -((1-iemin)/2);
 118     b1    = std::pow(ibeta, iexp);  // lower boundary of midrange
 119     iexp  = (iemax + 1 - it)/2;
 120     b2    = std::pow(ibeta,iexp);   // upper boundary of midrange
 121
 122     iexp  = (2-iemin)/2;
 123     s1m   = std::pow(ibeta,iexp);   // scaling factor for lower range
 124     iexp  = - ((iemax+it)/2);
 125     s2m   = std::pow(ibeta,iexp);   // scaling factor for upper range
 126
 127     overfl  = rbig*s2m;          // overfow boundary for abig
 128     eps     = std::pow(ibeta, 1-it);
 129     relerr  = internal::sqrt(eps);      // tolerance for neglecting asml
 130     abig    = 1.0/eps - 1.0;
 131     if (Scalar(nbig)>abig)  nmax = abig;  // largest safe n
 132     else                    nmax = nbig;
 133   }
 134
 135   typedef typename internal::packet_traits<Scalar>::type Packet;
 136   const int ps = internal::packet_traits<Scalar>::size;
 137   Packet pasml = internal::pset1(Scalar(0));
 138   Packet pamed = internal::pset1(Scalar(0));
 139   Packet pabig = internal::pset1(Scalar(0));
 140   Packet ps2m = internal::pset1(s2m);
 141   Packet ps1m = internal::pset1(s1m);
 142   Packet pb2  = internal::pset1(b2);
 143   Packet pb1  = internal::pset1(b1);
 144   for(int j=0; j<v.size(); j+=ps)
 145   {
 146     Packet ax = internal::pabs(v.template packet<Aligned>(j));
 147     Packet ax_s2m = internal::pmul(ax,ps2m);
 148     Packet ax_s1m = internal::pmul(ax,ps1m);
 149     Packet maskBig = internal::plt(pb2,ax);
 150     Packet maskSml = internal::plt(ax,pb1);
 151
 152 //     Packet maskMed = internal::pand(maskSml,maskBig);
 153 //     Packet scale = internal::pset1(Scalar(0));
 154 //     scale = internal::por(scale, internal::pand(maskBig,ps2m));
 155 //     scale = internal::por(scale, internal::pand(maskSml,ps1m));
 156 //     scale = internal::por(scale, internal::pandnot(internal::pset1(Scalar(1)),maskMed));
 157 //     ax = internal::pmul(ax,scale);
 158 //     ax = internal::pmul(ax,ax);
 159 //     pabig = internal::padd(pabig, internal::pand(maskBig, ax));
 160 //     pasml = internal::padd(pasml, internal::pand(maskSml, ax));
 161 //     pamed = internal::padd(pamed, internal::pandnot(ax,maskMed));
 162
 163
 164     pabig = internal::padd(pabig, internal::pand(maskBig, internal::pmul(ax_s2m,ax_s2m)));
 165     pasml = internal::padd(pasml, internal::pand(maskSml, internal::pmul(ax_s1m,ax_s1m)));
 166     pamed = internal::padd(pamed, internal::pandnot(internal::pmul(ax,ax),internal::pand(maskSml,maskBig)));
 167   }
 168   Scalar abig = internal::predux(pabig);
 169   Scalar asml = internal::predux(pasml);
 170   Scalar amed = internal::predux(pamed);
 171   if(abig > Scalar(0))
 172   {
 173     abig = internal::sqrt(abig);
 174     if(abig > overfl)
 175     {
 176       eigen_assert(false && "overflow");
 177       return rbig;
 178     }
 179     if(amed > Scalar(0))
 180     {
 181       abig = abig/s2m;
 182       amed = internal::sqrt(amed);
 183     }
 184     else
 185     {
 186       return abig/s2m;
 187     }
 188
 189   }
 190   else if(asml > Scalar(0))
 191   {
 192     if (amed > Scalar(0))
 193     {
 194       abig = internal::sqrt(amed);
 195       amed = internal::sqrt(asml) / s1m;
 196     }
 197     else
 198     {
 199       return internal::sqrt(asml)/s1m;
 200     }
 201   }
 202   else
 203   {
 204     return internal::sqrt(amed);
 205   }
 206   asml = std::min(abig, amed);
 207   abig = std::max(abig, amed);
 208   if(asml <= abig*relerr)
 209     return abig;
 210   else
 211     return abig * internal::sqrt(Scalar(1) + internal::abs2(asml/abig));
 212   #endif
 213 }
 214
 215 #define BENCH_PERF(NRM) { \
 216   Eigen::BenchTimer tf, td, tcf; tf.reset(); td.reset(); tcf.reset();\
 217   for (int k=0; k<tries; ++k) { \
 218     tf.start(); \
 219     for (int i=0; i<iters; ++i) NRM(vf); \
 220     tf.stop(); \
 221   } \
 222   for (int k=0; k<tries; ++k) { \
 223     td.start(); \
 224     for (int i=0; i<iters; ++i) NRM(vd); \
 225     td.stop(); \
 226   } \
 227   for (int k=0; k<std::max(1,tries/3); ++k) { \
 228     tcf.start(); \
 229     for (int i=0; i<iters; ++i) NRM(vcf); \
 230     tcf.stop(); \
 231   } \
 232   std::cout << #NRM << "\t" << tf.value() << "   " << td.value() <<  "    " << tcf.value() << "\n"; \
 233 }
 234
 235 void check_accuracy(double basef, double based, int s)
 236 {
 237   double yf = basef * internal::abs(internal::random<double>());
 238   double yd = based * internal::abs(internal::random<double>());
 239   VectorXf vf = VectorXf::Ones(s) * yf;
 240   VectorXd vd = VectorXd::Ones(s) * yd;
 241
 242   std::cout << "reference\t" << internal::sqrt(double(s))*yf << "\t" << internal::sqrt(double(s))*yd << "\n";
 243   std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\n";
 244   std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\n";
 245   std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\n";
 246   std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\n";
 247   std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\n";
 248   std::cout << "twopassNorm\t" << twopassNorm(vf) << "\t" << twopassNorm(vd) << "\n";
 249   std::cout << "bl2passNorm\t" << bl2passNorm(vf) << "\t" << bl2passNorm(vd) << "\n";
 250 }
 251
 252 void check_accuracy_var(int ef0, int ef1, int ed0, int ed1, int s)
 253 {
 254   VectorXf vf(s);
 255   VectorXd vd(s);
 256   for (int i=0; i<s; ++i)
 257   {
 258     vf[i] = internal::abs(internal::random<double>()) * std::pow(double(10), internal::random<int>(ef0,ef1));
 259     vd[i] = internal::abs(internal::random<double>()) * std::pow(double(10), internal::random<int>(ed0,ed1));
 260   }
 261
 262   //std::cout << "reference\t" << internal::sqrt(double(s))*yf << "\t" << internal::sqrt(double(s))*yd << "\n";
 263   std::cout << "sqsumNorm\t"  << sqsumNorm(vf)  << "\t" << sqsumNorm(vd)  << "\t" << sqsumNorm(vf.cast<long double>()) << "\t" << sqsumNorm(vd.cast<long double>()) << "\n";
 264   std::cout << "hypotNorm\t"  << hypotNorm(vf)  << "\t" << hypotNorm(vd)  << "\t" << hypotNorm(vf.cast<long double>()) << "\t" << hypotNorm(vd.cast<long double>()) << "\n";
 265   std::cout << "blueNorm\t"   << blueNorm(vf)   << "\t" << blueNorm(vd)   << "\t" << blueNorm(vf.cast<long double>()) << "\t" << blueNorm(vd.cast<long double>()) << "\n";
 266   std::cout << "pblueNorm\t"  << pblueNorm(vf)  << "\t" << pblueNorm(vd)  << "\t" << blueNorm(vf.cast<long double>()) << "\t" << blueNorm(vd.cast<long double>()) << "\n";
 267   std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\t" << lapackNorm(vf.cast<long double>()) << "\t" << lapackNorm(vd.cast<long double>()) << "\n";
 268   std::cout << "twopassNorm\t" << twopassNorm(vf) << "\t" << twopassNorm(vd) << "\t" << twopassNorm(vf.cast<long double>()) << "\t" << twopassNorm(vd.cast<long double>()) << "\n";
 269 //   std::cout << "bl2passNorm\t" << bl2passNorm(vf) << "\t" << bl2passNorm(vd) << "\t" << bl2passNorm(vf.cast<long double>()) << "\t" << bl2passNorm(vd.cast<long double>()) << "\n";
 270 }
 271
 272 int main(int argc, char** argv)
 273 {
 274   int tries = 10;
 275   int iters = 100000;
 276   double y = 1.1345743233455785456788e12 * internal::random<double>();
 277   VectorXf v = VectorXf::Ones(1024) * y;
 278
 279 // return 0;
 280   int s = 10000;
 281   double basef_ok = 1.1345743233455785456788e15;
 282   double based_ok = 1.1345743233455785456788e95;
 283
 284   double basef_under = 1.1345743233455785456788e-27;
 285   double based_under = 1.1345743233455785456788e-303;
 286
 287   double basef_over = 1.1345743233455785456788e+27;
 288   double based_over = 1.1345743233455785456788e+302;
 289
 290   std::cout.precision(20);
 291
 292   std::cerr << "\nNo under/overflow:\n";
 293   check_accuracy(basef_ok, based_ok, s);
 294
 295   std::cerr << "\nUnderflow:\n";
 296   check_accuracy(basef_under, based_under, s);
 297
 298   std::cerr << "\nOverflow:\n";
 299   check_accuracy(basef_over, based_over, s);
 300
 301   std::cerr << "\nVarying (over):\n";
 302   for (int k=0; k<1; ++k)
 303   {
 304     check_accuracy_var(20,27,190,302,s);
 305     std::cout << "\n";
 306   }
 307
 308   std::cerr << "\nVarying (under):\n";
 309   for (int k=0; k<1; ++k)
 310   {
 311     check_accuracy_var(-27,20,-302,-190,s);
 312     std::cout << "\n";
 313   }
 314
 315   std::cout.precision(4);
 316   std::cerr << "Performance (out of cache):\n";
 317   {
 318     int iters = 1;
 319     VectorXf vf = VectorXf::Random(1024*1024*32) * y;
 320     VectorXd vd = VectorXd::Random(1024*1024*32) * y;
 321     VectorXcf vcf = VectorXcf::Random(1024*1024*32) * y;
 322     BENCH_PERF(sqsumNorm);
 323     BENCH_PERF(blueNorm);
 324 //     BENCH_PERF(pblueNorm);
 325 //     BENCH_PERF(lapackNorm);
 326 //     BENCH_PERF(hypotNorm);
 327 //     BENCH_PERF(twopassNorm);
 328     BENCH_PERF(bl2passNorm);
 329   }
 330
 331   std::cerr << "\nPerformance (in cache):\n";
 332   {
 333     int iters = 100000;
 334     VectorXf vf = VectorXf::Random(512) * y;
 335     VectorXd vd = VectorXd::Random(512) * y;
 336     VectorXcf vcf = VectorXcf::Random(512) * y;
 337     BENCH_PERF(sqsumNorm);
 338     BENCH_PERF(blueNorm);
 339 //     BENCH_PERF(pblueNorm);
 340 //     BENCH_PERF(lapackNorm);
 341 //     BENCH_PERF(hypotNorm);
 342 //     BENCH_PERF(twopassNorm);
 343     BENCH_PERF(bl2passNorm);
 344   }
 345 }