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