Squashed 'third_party/eigen/' content from commit 61d72f6
Change-Id: Iccc90fa0b55ab44037f018046d2fcffd90d9d025
git-subtree-dir: third_party/eigen
git-subtree-split: 61d72f6383cfa842868c53e30e087b0258177257
diff --git a/bench/bench_norm.cpp b/bench/bench_norm.cpp
new file mode 100644
index 0000000..806db29
--- /dev/null
+++ b/bench/bench_norm.cpp
@@ -0,0 +1,345 @@
+#include <typeinfo>
+#include <iostream>
+#include <Eigen/Core>
+#include "BenchTimer.h"
+using namespace Eigen;
+using namespace std;
+
+template<typename T>
+EIGEN_DONT_INLINE typename T::Scalar sqsumNorm(const T& v)
+{
+ return v.norm();
+}
+
+template<typename T>
+EIGEN_DONT_INLINE typename T::Scalar hypotNorm(const T& v)
+{
+ return v.hypotNorm();
+}
+
+template<typename T>
+EIGEN_DONT_INLINE typename T::Scalar blueNorm(const T& v)
+{
+ return v.blueNorm();
+}
+
+template<typename T>
+EIGEN_DONT_INLINE typename T::Scalar lapackNorm(T& v)
+{
+ typedef typename T::Scalar Scalar;
+ int n = v.size();
+ Scalar scale = 0;
+ Scalar ssq = 1;
+ for (int i=0;i<n;++i)
+ {
+ Scalar ax = internal::abs(v.coeff(i));
+ if (scale >= ax)
+ {
+ ssq += internal::abs2(ax/scale);
+ }
+ else
+ {
+ ssq = Scalar(1) + ssq * internal::abs2(scale/ax);
+ scale = ax;
+ }
+ }
+ return scale * internal::sqrt(ssq);
+}
+
+template<typename T>
+EIGEN_DONT_INLINE typename T::Scalar twopassNorm(T& v)
+{
+ typedef typename T::Scalar Scalar;
+ Scalar s = v.cwise().abs().maxCoeff();
+ return s*(v/s).norm();
+}
+
+template<typename T>
+EIGEN_DONT_INLINE typename T::Scalar bl2passNorm(T& v)
+{
+ return v.stableNorm();
+}
+
+template<typename T>
+EIGEN_DONT_INLINE typename T::Scalar divacNorm(T& v)
+{
+ int n =v.size() / 2;
+ for (int i=0;i<n;++i)
+ v(i) = v(2*i)*v(2*i) + v(2*i+1)*v(2*i+1);
+ n = n/2;
+ while (n>0)
+ {
+ for (int i=0;i<n;++i)
+ v(i) = v(2*i) + v(2*i+1);
+ n = n/2;
+ }
+ return internal::sqrt(v(0));
+}
+
+#ifdef EIGEN_VECTORIZE
+Packet4f internal::plt(const Packet4f& a, Packet4f& b) { return _mm_cmplt_ps(a,b); }
+Packet2d internal::plt(const Packet2d& a, Packet2d& b) { return _mm_cmplt_pd(a,b); }
+
+Packet4f internal::pandnot(const Packet4f& a, Packet4f& b) { return _mm_andnot_ps(a,b); }
+Packet2d internal::pandnot(const Packet2d& a, Packet2d& b) { return _mm_andnot_pd(a,b); }
+#endif
+
+template<typename T>
+EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v)
+{
+ #ifndef EIGEN_VECTORIZE
+ return v.blueNorm();
+ #else
+ typedef typename T::Scalar Scalar;
+
+ static int nmax = 0;
+ static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr;
+ int n;
+
+ if(nmax <= 0)
+ {
+ int nbig, ibeta, it, iemin, iemax, iexp;
+ Scalar abig, eps;
+
+ nbig = std::numeric_limits<int>::max(); // largest integer
+ ibeta = std::numeric_limits<Scalar>::radix; //NumTraits<Scalar>::Base; // base for floating-point numbers
+ it = std::numeric_limits<Scalar>::digits; //NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa
+ iemin = std::numeric_limits<Scalar>::min_exponent; // minimum exponent
+ iemax = std::numeric_limits<Scalar>::max_exponent; // maximum exponent
+ rbig = std::numeric_limits<Scalar>::max(); // largest floating-point number
+
+ // Check the basic machine-dependent constants.
+ if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5)
+ || (it<=4 && ibeta <= 3 ) || it<2)
+ {
+ eigen_assert(false && "the algorithm cannot be guaranteed on this computer");
+ }
+ iexp = -((1-iemin)/2);
+ b1 = std::pow(ibeta, iexp); // lower boundary of midrange
+ iexp = (iemax + 1 - it)/2;
+ b2 = std::pow(ibeta,iexp); // upper boundary of midrange
+
+ iexp = (2-iemin)/2;
+ s1m = std::pow(ibeta,iexp); // scaling factor for lower range
+ iexp = - ((iemax+it)/2);
+ s2m = std::pow(ibeta,iexp); // scaling factor for upper range
+
+ overfl = rbig*s2m; // overfow boundary for abig
+ eps = std::pow(ibeta, 1-it);
+ relerr = internal::sqrt(eps); // tolerance for neglecting asml
+ abig = 1.0/eps - 1.0;
+ if (Scalar(nbig)>abig) nmax = abig; // largest safe n
+ else nmax = nbig;
+ }
+
+ typedef typename internal::packet_traits<Scalar>::type Packet;
+ const int ps = internal::packet_traits<Scalar>::size;
+ Packet pasml = internal::pset1(Scalar(0));
+ Packet pamed = internal::pset1(Scalar(0));
+ Packet pabig = internal::pset1(Scalar(0));
+ Packet ps2m = internal::pset1(s2m);
+ Packet ps1m = internal::pset1(s1m);
+ Packet pb2 = internal::pset1(b2);
+ Packet pb1 = internal::pset1(b1);
+ for(int j=0; j<v.size(); j+=ps)
+ {
+ Packet ax = internal::pabs(v.template packet<Aligned>(j));
+ Packet ax_s2m = internal::pmul(ax,ps2m);
+ Packet ax_s1m = internal::pmul(ax,ps1m);
+ Packet maskBig = internal::plt(pb2,ax);
+ Packet maskSml = internal::plt(ax,pb1);
+
+// Packet maskMed = internal::pand(maskSml,maskBig);
+// Packet scale = internal::pset1(Scalar(0));
+// scale = internal::por(scale, internal::pand(maskBig,ps2m));
+// scale = internal::por(scale, internal::pand(maskSml,ps1m));
+// scale = internal::por(scale, internal::pandnot(internal::pset1(Scalar(1)),maskMed));
+// ax = internal::pmul(ax,scale);
+// ax = internal::pmul(ax,ax);
+// pabig = internal::padd(pabig, internal::pand(maskBig, ax));
+// pasml = internal::padd(pasml, internal::pand(maskSml, ax));
+// pamed = internal::padd(pamed, internal::pandnot(ax,maskMed));
+
+
+ pabig = internal::padd(pabig, internal::pand(maskBig, internal::pmul(ax_s2m,ax_s2m)));
+ pasml = internal::padd(pasml, internal::pand(maskSml, internal::pmul(ax_s1m,ax_s1m)));
+ pamed = internal::padd(pamed, internal::pandnot(internal::pmul(ax,ax),internal::pand(maskSml,maskBig)));
+ }
+ Scalar abig = internal::predux(pabig);
+ Scalar asml = internal::predux(pasml);
+ Scalar amed = internal::predux(pamed);
+ if(abig > Scalar(0))
+ {
+ abig = internal::sqrt(abig);
+ if(abig > overfl)
+ {
+ eigen_assert(false && "overflow");
+ return rbig;
+ }
+ if(amed > Scalar(0))
+ {
+ abig = abig/s2m;
+ amed = internal::sqrt(amed);
+ }
+ else
+ {
+ return abig/s2m;
+ }
+
+ }
+ else if(asml > Scalar(0))
+ {
+ if (amed > Scalar(0))
+ {
+ abig = internal::sqrt(amed);
+ amed = internal::sqrt(asml) / s1m;
+ }
+ else
+ {
+ return internal::sqrt(asml)/s1m;
+ }
+ }
+ else
+ {
+ return internal::sqrt(amed);
+ }
+ asml = std::min(abig, amed);
+ abig = std::max(abig, amed);
+ if(asml <= abig*relerr)
+ return abig;
+ else
+ return abig * internal::sqrt(Scalar(1) + internal::abs2(asml/abig));
+ #endif
+}
+
+#define BENCH_PERF(NRM) { \
+ Eigen::BenchTimer tf, td, tcf; tf.reset(); td.reset(); tcf.reset();\
+ for (int k=0; k<tries; ++k) { \
+ tf.start(); \
+ for (int i=0; i<iters; ++i) NRM(vf); \
+ tf.stop(); \
+ } \
+ for (int k=0; k<tries; ++k) { \
+ td.start(); \
+ for (int i=0; i<iters; ++i) NRM(vd); \
+ td.stop(); \
+ } \
+ for (int k=0; k<std::max(1,tries/3); ++k) { \
+ tcf.start(); \
+ for (int i=0; i<iters; ++i) NRM(vcf); \
+ tcf.stop(); \
+ } \
+ std::cout << #NRM << "\t" << tf.value() << " " << td.value() << " " << tcf.value() << "\n"; \
+}
+
+void check_accuracy(double basef, double based, int s)
+{
+ double yf = basef * internal::abs(internal::random<double>());
+ double yd = based * internal::abs(internal::random<double>());
+ VectorXf vf = VectorXf::Ones(s) * yf;
+ VectorXd vd = VectorXd::Ones(s) * yd;
+
+ std::cout << "reference\t" << internal::sqrt(double(s))*yf << "\t" << internal::sqrt(double(s))*yd << "\n";
+ std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\n";
+ std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\n";
+ std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\n";
+ std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\n";
+ std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\n";
+ std::cout << "twopassNorm\t" << twopassNorm(vf) << "\t" << twopassNorm(vd) << "\n";
+ std::cout << "bl2passNorm\t" << bl2passNorm(vf) << "\t" << bl2passNorm(vd) << "\n";
+}
+
+void check_accuracy_var(int ef0, int ef1, int ed0, int ed1, int s)
+{
+ VectorXf vf(s);
+ VectorXd vd(s);
+ for (int i=0; i<s; ++i)
+ {
+ vf[i] = internal::abs(internal::random<double>()) * std::pow(double(10), internal::random<int>(ef0,ef1));
+ vd[i] = internal::abs(internal::random<double>()) * std::pow(double(10), internal::random<int>(ed0,ed1));
+ }
+
+ //std::cout << "reference\t" << internal::sqrt(double(s))*yf << "\t" << internal::sqrt(double(s))*yd << "\n";
+ std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\t" << sqsumNorm(vf.cast<long double>()) << "\t" << sqsumNorm(vd.cast<long double>()) << "\n";
+ std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\t" << hypotNorm(vf.cast<long double>()) << "\t" << hypotNorm(vd.cast<long double>()) << "\n";
+ std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\t" << blueNorm(vf.cast<long double>()) << "\t" << blueNorm(vd.cast<long double>()) << "\n";
+ std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\t" << blueNorm(vf.cast<long double>()) << "\t" << blueNorm(vd.cast<long double>()) << "\n";
+ std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\t" << lapackNorm(vf.cast<long double>()) << "\t" << lapackNorm(vd.cast<long double>()) << "\n";
+ std::cout << "twopassNorm\t" << twopassNorm(vf) << "\t" << twopassNorm(vd) << "\t" << twopassNorm(vf.cast<long double>()) << "\t" << twopassNorm(vd.cast<long double>()) << "\n";
+// std::cout << "bl2passNorm\t" << bl2passNorm(vf) << "\t" << bl2passNorm(vd) << "\t" << bl2passNorm(vf.cast<long double>()) << "\t" << bl2passNorm(vd.cast<long double>()) << "\n";
+}
+
+int main(int argc, char** argv)
+{
+ int tries = 10;
+ int iters = 100000;
+ double y = 1.1345743233455785456788e12 * internal::random<double>();
+ VectorXf v = VectorXf::Ones(1024) * y;
+
+// return 0;
+ int s = 10000;
+ double basef_ok = 1.1345743233455785456788e15;
+ double based_ok = 1.1345743233455785456788e95;
+
+ double basef_under = 1.1345743233455785456788e-27;
+ double based_under = 1.1345743233455785456788e-303;
+
+ double basef_over = 1.1345743233455785456788e+27;
+ double based_over = 1.1345743233455785456788e+302;
+
+ std::cout.precision(20);
+
+ std::cerr << "\nNo under/overflow:\n";
+ check_accuracy(basef_ok, based_ok, s);
+
+ std::cerr << "\nUnderflow:\n";
+ check_accuracy(basef_under, based_under, s);
+
+ std::cerr << "\nOverflow:\n";
+ check_accuracy(basef_over, based_over, s);
+
+ std::cerr << "\nVarying (over):\n";
+ for (int k=0; k<1; ++k)
+ {
+ check_accuracy_var(20,27,190,302,s);
+ std::cout << "\n";
+ }
+
+ std::cerr << "\nVarying (under):\n";
+ for (int k=0; k<1; ++k)
+ {
+ check_accuracy_var(-27,20,-302,-190,s);
+ std::cout << "\n";
+ }
+
+ std::cout.precision(4);
+ std::cerr << "Performance (out of cache):\n";
+ {
+ int iters = 1;
+ VectorXf vf = VectorXf::Random(1024*1024*32) * y;
+ VectorXd vd = VectorXd::Random(1024*1024*32) * y;
+ VectorXcf vcf = VectorXcf::Random(1024*1024*32) * y;
+ BENCH_PERF(sqsumNorm);
+ BENCH_PERF(blueNorm);
+// BENCH_PERF(pblueNorm);
+// BENCH_PERF(lapackNorm);
+// BENCH_PERF(hypotNorm);
+// BENCH_PERF(twopassNorm);
+ BENCH_PERF(bl2passNorm);
+ }
+
+ std::cerr << "\nPerformance (in cache):\n";
+ {
+ int iters = 100000;
+ VectorXf vf = VectorXf::Random(512) * y;
+ VectorXd vd = VectorXd::Random(512) * y;
+ VectorXcf vcf = VectorXcf::Random(512) * y;
+ BENCH_PERF(sqsumNorm);
+ BENCH_PERF(blueNorm);
+// BENCH_PERF(pblueNorm);
+// BENCH_PERF(lapackNorm);
+// BENCH_PERF(hypotNorm);
+// BENCH_PERF(twopassNorm);
+ BENCH_PERF(bl2passNorm);
+ }
+}