From 65fc70b75039a5cdfc5df67f62d38b317196293b Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 16 Jul 2009 11:33:56 +0200 Subject: [PATCH 01/11] add a benchmark for the different norms --- bench/bench_norm.cpp | 259 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 259 insertions(+) create mode 100644 bench/bench_norm.cpp diff --git a/bench/bench_norm.cpp b/bench/bench_norm.cpp new file mode 100644 index 000000000..76c8c574d --- /dev/null +++ b/bench/bench_norm.cpp @@ -0,0 +1,259 @@ +#include +#include "BenchTimer.h" +using namespace Eigen; +using namespace std; + +template +EIGEN_DONT_INLINE typename T::Scalar sqsumNorm(const T& v) +{ + return v.norm(); +} + +template +EIGEN_DONT_INLINE typename T::Scalar hypotNorm(const T& v) +{ + return v.stableNorm(); +} + +template +EIGEN_DONT_INLINE typename T::Scalar blueNorm(const T& v) +{ + return v.blueNorm(); +} + +template +EIGEN_DONT_INLINE typename T::Scalar lapackNorm(T& v) +{ + typedef typename T::Scalar Scalar; + int n = v.size(); + Scalar scale = 1; + Scalar ssq = 0; + for (int i=0;i +EIGEN_DONT_INLINE typename T::Scalar divacNorm(T& v) +{ + int n =v.size() / 2; + for (int i=0;i0) + { + for (int i=0;i +EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v) +{ + typedef typename T::Scalar Scalar; + + static int nmax; + 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::max(); // largest integer + ibeta = NumTraits::Base; // base for floating-point numbers + it = NumTraits::Mantissa; // number of base-beta digits in mantissa + iemin = std::numeric_limits::min_exponent; // minimum exponent + iemax = std::numeric_limits::max_exponent; // maximum exponent + rbig = std::numeric_limits::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) + { + ei_assert(false && "the algorithm cannot be guaranteed on this computer"); + } + iexp = -((1-iemin)/2); + b1 = bexp(ibeta, iexp); // lower boundary of midrange + iexp = (iemax + 1 - it)/2; + b2 = bexp(ibeta,iexp); // upper boundary of midrange + + iexp = (2-iemin)/2; + s1m = bexp(ibeta,iexp); // scaling factor for lower range + iexp = - ((iemax+it)/2); + s2m = bexp(ibeta,iexp); // scaling factor for upper range + + overfl = rbig*s2m; // overfow boundary for abig + eps = bexp(ibeta, 1-it); + relerr = ei_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 ei_packet_traits::type Packet; + const int ps = ei_packet_traits::size; + Packet pasml = ei_pset1(Scalar(0)); + Packet pamed = ei_pset1(Scalar(0)); + Packet pabig = ei_pset1(Scalar(0)); + Packet ps2m = ei_pset1(s2m); + Packet ps1m = ei_pset1(s1m); + Packet pb2 = ei_pset1(b2); + Packet pb1 = ei_pset1(b1); + for(int j=0; j(j)); + Packet ax_s2m = ei_pmul(ax,ps2m); + Packet ax_s1m = ei_pmul(ax,ps1m); + Packet maskBig = ei_plt(pb2,ax); + Packet maskSml = ei_plt(ax,pb1); + pabig = ei_padd(pabig, ei_pand(maskBig, ei_pmul(ax_s2m,ax_s2m))); + pasml = ei_padd(pasml, ei_pand(maskSml, ei_pmul(ax_s1m,ax_s1m))); + pamed = ei_padd(pamed, ei_pandnot(ei_pmul(ax,ax),ei_pand(maskSml,maskBig))); + } + Scalar abig = ei_predux(pabig); + Scalar asml = ei_predux(pasml); + Scalar amed = ei_predux(pamed); + if(abig > Scalar(0)) + { + abig = ei_sqrt(abig); + if(abig > overfl) + { + ei_assert(false && "overflow"); + return rbig; + } + if(amed > Scalar(0)) + { + abig = abig/s2m; + amed = ei_sqrt(amed); + } + else + { + return abig/s2m; + } + + } + else if(asml > Scalar(0)) + { + if (amed > Scalar(0)) + { + abig = ei_sqrt(amed); + amed = ei_sqrt(asml) / s1m; + } + else + { + return ei_sqrt(asml)/s1m; + } + } + else + { + return ei_sqrt(amed); + } + asml = std::min(abig, amed); + abig = std::max(abig, amed); + if(asml <= abig*relerr) + return abig; + else + return abig * ei_sqrt(Scalar(1) + ei_abs2(asml/abig)); +} + +#define BENCH_PERF(NRM) { \ + Eigen::BenchTimer tf, td; tf.reset(); td.reset();\ + for (int k=0; k()); + double yd = based * ei_abs(ei_random()); + VectorXf vf = VectorXf::Ones(s) * yf; + VectorXd vd = VectorXd::Ones(s) * yd; + + std::cout << "reference\t" << ei_sqrt(double(s))*yf << "\t" << ei_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"; +} + +int main(int argc, char** argv) +{ + int tries = 5; + int iters = 100000; + double y = 1.1345743233455785456788e12 * ei_random(); + VectorXf v = VectorXf::Ones(1024) * y; + +// std::cerr << "Performance (out of cache):\n"; +// { +// int iters = 1; +// VectorXf vf = VectorXf::Ones(1024*1024*32) * y; +// VectorXd vd = VectorXd::Ones(1024*1024*32) * y; +// BENCH_PERF(sqsumNorm); +// BENCH_PERF(blueNorm); +// BENCH_PERF(pblueNorm); +// BENCH_PERF(lapackNorm); +// BENCH_PERF(hypotNorm); +// } +// +// std::cerr << "\nPerformance (in cache):\n"; +// { +// int iters = 100000; +// VectorXf vf = VectorXf::Ones(512) * y; +// VectorXd vd = VectorXd::Ones(512) * y; +// BENCH_PERF(sqsumNorm); +// BENCH_PERF(blueNorm); +// BENCH_PERF(pblueNorm); +// BENCH_PERF(lapackNorm); +// BENCH_PERF(hypotNorm); +// } + + int s = 10000; + double basef_ok = 1.1345743233455785456788e12; + double based_ok = 1.1345743233455785456788e32; + + double basef_under = 1.1345743233455785456788e-23; + double based_under = 1.1345743233455785456788e-180; + + double basef_over = 1.1345743233455785456788e+27; + double based_over = 1.1345743233455785456788e+185; + + 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, 1); + + std::cerr << "\nOverflow:\n"; + check_accuracy(basef_over, based_over, s); +} From 525da6a464c897d0fe4e401a65851bcd7567fc5a Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 16 Jul 2009 14:20:36 +0200 Subject: [PATCH 02/11] bugfix in blueNorm --- Eigen/src/Core/Dot.h | 39 +++--------- Eigen/src/Core/NumTraits.h | 8 +-- bench/bench_norm.cpp | 126 +++++++++++++++++++++++-------------- 3 files changed, 89 insertions(+), 84 deletions(-) diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index c5f2e8505..d97f5837f 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -304,29 +304,6 @@ MatrixBase::stableNorm() const return this->cwise().abs().redux(ei_scalar_hypot_op()); } -/** \internal Computes ibeta^iexp by binary expansion of iexp, - * exact if ibeta is the machine base */ -template inline T bexp(int ibeta, int iexp) -{ - T tbeta = T(ibeta); - T res = 1.0; - int n = iexp; - if (n<0) - { - n = - n; - tbeta = 1.0/tbeta; - } - for(;;) - { - if ((n % 2)==0) - res = res * tbeta; - n = n/2; - if (n==0) return res; - tbeta = tbeta*tbeta; - } - return res; -} - /** \returns the \em l2 norm of \c *this using the Blue's algorithm. * A Portable Fortran Program to Find the Euclidean Norm of a Vector, * ACM TOMS, Vol 4, Issue 1, 1978. @@ -337,7 +314,7 @@ template inline typename NumTraits::Scalar>::Real MatrixBase::blueNorm() const { - static int nmax; + static int nmax = -1; static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr; int n; Scalar ax, abig, amed, asml; @@ -355,8 +332,8 @@ MatrixBase::blueNorm() const // are used. For any specific computer, each of the assignment // statements can be replaced nbig = std::numeric_limits::max(); // largest integer - ibeta = NumTraits::Base; // base for floating-point numbers - it = NumTraits::Mantissa; // number of base-beta digits in mantissa + ibeta = std::numeric_limits::radix; //NumTraits::Base; // base for floating-point numbers + it = std::numeric_limits::digits; //NumTraits::Mantissa; // number of base-beta digits in mantissa iemin = std::numeric_limits::min_exponent; // minimum exponent iemax = std::numeric_limits::max_exponent; // maximum exponent rbig = std::numeric_limits::max(); // largest floating-point number @@ -368,17 +345,17 @@ MatrixBase::blueNorm() const ei_assert(false && "the algorithm cannot be guaranteed on this computer"); } iexp = -((1-iemin)/2); - b1 = bexp(ibeta, iexp); // lower boundary of midrange + b1 = std::pow(ibeta, iexp); // lower boundary of midrange iexp = (iemax + 1 - it)/2; - b2 = bexp(ibeta,iexp); // upper boundary of midrange + b2 = std::pow(ibeta,iexp); // upper boundary of midrange iexp = (2-iemin)/2; - s1m = bexp(ibeta,iexp); // scaling factor for lower range + s1m = std::pow(ibeta,iexp); // scaling factor for lower range iexp = - ((iemax+it)/2); - s2m = bexp(ibeta,iexp); // scaling factor for upper range + s2m = std::pow(ibeta,iexp); // scaling factor for upper range overfl = rbig*s2m; // overfow boundary for abig - eps = bexp(ibeta, 1-it); + eps = std::pow(ibeta, 1-it); relerr = ei_sqrt(eps); // tolerance for neglecting asml abig = 1.0/eps - 1.0; if (Scalar(nbig)>abig) nmax = abig; // largest safe n diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h index dec07a036..24afe54c5 100644 --- a/Eigen/src/Core/NumTraits.h +++ b/Eigen/src/Core/NumTraits.h @@ -70,9 +70,7 @@ template<> struct NumTraits HasFloatingPoint = 1, ReadCost = 1, AddCost = 1, - MulCost = 1, - Base = 2, - Mantissa = 23 + MulCost = 1 }; }; @@ -85,9 +83,7 @@ template<> struct NumTraits HasFloatingPoint = 1, ReadCost = 1, AddCost = 1, - MulCost = 1, - Base = 2, - Mantissa = 52 + MulCost = 1 }; }; diff --git a/bench/bench_norm.cpp b/bench/bench_norm.cpp index 76c8c574d..e06d06417 100644 --- a/bench/bench_norm.cpp +++ b/bench/bench_norm.cpp @@ -1,3 +1,4 @@ +#include #include #include "BenchTimer.h" using namespace Eigen; @@ -58,18 +59,23 @@ EIGEN_DONT_INLINE typename T::Scalar divacNorm(T& v) return ei_sqrt(v(0)); } +#ifdef EIGEN_VECTORIZE Packet4f ei_plt(const Packet4f& a, Packet4f& b) { return _mm_cmplt_ps(a,b); } Packet2d ei_plt(const Packet2d& a, Packet2d& b) { return _mm_cmplt_pd(a,b); } Packet4f ei_pandnot(const Packet4f& a, Packet4f& b) { return _mm_andnot_ps(a,b); } Packet2d ei_pandnot(const Packet2d& a, Packet2d& b) { return _mm_andnot_pd(a,b); } +#endif template 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; + static int nmax = 0; static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr; int n; @@ -79,8 +85,8 @@ EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v) Scalar abig, eps; nbig = std::numeric_limits::max(); // largest integer - ibeta = NumTraits::Base; // base for floating-point numbers - it = NumTraits::Mantissa; // number of base-beta digits in mantissa + ibeta = std::numeric_limits::radix; //NumTraits::Base; // base for floating-point numbers + it = std::numeric_limits::digits; //NumTraits::Mantissa; // number of base-beta digits in mantissa iemin = std::numeric_limits::min_exponent; // minimum exponent iemax = std::numeric_limits::max_exponent; // maximum exponent rbig = std::numeric_limits::max(); // largest floating-point number @@ -92,23 +98,23 @@ EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v) ei_assert(false && "the algorithm cannot be guaranteed on this computer"); } iexp = -((1-iemin)/2); - b1 = bexp(ibeta, iexp); // lower boundary of midrange + b1 = std::pow(ibeta, iexp); // lower boundary of midrange iexp = (iemax + 1 - it)/2; - b2 = bexp(ibeta,iexp); // upper boundary of midrange + b2 = std::pow(ibeta,iexp); // upper boundary of midrange iexp = (2-iemin)/2; - s1m = bexp(ibeta,iexp); // scaling factor for lower range + s1m = std::pow(ibeta,iexp); // scaling factor for lower range iexp = - ((iemax+it)/2); - s2m = bexp(ibeta,iexp); // scaling factor for upper range + s2m = std::pow(ibeta,iexp); // scaling factor for upper range overfl = rbig*s2m; // overfow boundary for abig - eps = bexp(ibeta, 1-it); + eps = std::pow(ibeta, 1-it); relerr = ei_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 ei_packet_traits::type Packet; const int ps = ei_packet_traits::size; Packet pasml = ei_pset1(Scalar(0)); @@ -173,6 +179,7 @@ EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v) return abig; else return abig * ei_sqrt(Scalar(1) + ei_abs2(asml/abig)); + #endif } #define BENCH_PERF(NRM) { \ @@ -196,7 +203,7 @@ void check_accuracy(double basef, double based, int s) double yd = based * ei_abs(ei_random()); VectorXf vf = VectorXf::Ones(s) * yf; VectorXd vd = VectorXd::Ones(s) * yd; - + std::cout << "reference\t" << ei_sqrt(double(s))*yf << "\t" << ei_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"; @@ -205,55 +212,80 @@ void check_accuracy(double basef, double based, int s) std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\n"; } -int main(int argc, char** argv) +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()) * std::pow(double(10), ei_random(ef0,ef1)); + vd[i] = ei_abs(ei_random()) * std::pow(double(10), ei_random(ed0,ed1)); + } + + //std::cout << "reference\t" << ei_sqrt(double(s))*yf << "\t" << ei_sqrt(double(s))*yd << "\n"; + std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\t" << sqsumNorm(vf.cast()) << "\t" << sqsumNorm(vd.cast()) << "\n"; + std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\t" << hypotNorm(vf.cast()) << "\t" << hypotNorm(vd.cast()) << "\n"; + std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\t" << blueNorm(vf.cast()) << "\t" << blueNorm(vd.cast()) << "\n"; + std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\t" << blueNorm(vf.cast()) << "\t" << blueNorm(vd.cast()) << "\n"; + std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\t" << lapackNorm(vf.cast()) << "\t" << lapackNorm(vd.cast()) << "\n"; +} + +int main(int argc, char** argv) { int tries = 5; int iters = 100000; double y = 1.1345743233455785456788e12 * ei_random(); VectorXf v = VectorXf::Ones(1024) * y; - -// std::cerr << "Performance (out of cache):\n"; -// { -// int iters = 1; -// VectorXf vf = VectorXf::Ones(1024*1024*32) * y; -// VectorXd vd = VectorXd::Ones(1024*1024*32) * y; -// BENCH_PERF(sqsumNorm); -// BENCH_PERF(blueNorm); -// BENCH_PERF(pblueNorm); -// BENCH_PERF(lapackNorm); -// BENCH_PERF(hypotNorm); -// } -// -// std::cerr << "\nPerformance (in cache):\n"; -// { -// int iters = 100000; -// VectorXf vf = VectorXf::Ones(512) * y; -// VectorXd vd = VectorXd::Ones(512) * y; -// BENCH_PERF(sqsumNorm); -// BENCH_PERF(blueNorm); -// BENCH_PERF(pblueNorm); -// BENCH_PERF(lapackNorm); -// BENCH_PERF(hypotNorm); -// } - + + std::cerr << "Performance (out of cache):\n"; + { + int iters = 1; + VectorXf vf = VectorXf::Ones(1024*1024*32) * y; + VectorXd vd = VectorXd::Ones(1024*1024*32) * y; + BENCH_PERF(sqsumNorm); + BENCH_PERF(blueNorm); + BENCH_PERF(pblueNorm); + BENCH_PERF(lapackNorm); + BENCH_PERF(hypotNorm); + } + + std::cerr << "\nPerformance (in cache):\n"; + { + int iters = 100000; + VectorXf vf = VectorXf::Ones(512) * y; + VectorXd vd = VectorXd::Ones(512) * y; + BENCH_PERF(sqsumNorm); + BENCH_PERF(blueNorm); + BENCH_PERF(pblueNorm); + BENCH_PERF(lapackNorm); + BENCH_PERF(hypotNorm); + } + int s = 10000; - double basef_ok = 1.1345743233455785456788e12; - double based_ok = 1.1345743233455785456788e32; - - double basef_under = 1.1345743233455785456788e-23; - double based_under = 1.1345743233455785456788e-180; - + double basef_ok = 1.1345743233455785456788e15; + double based_ok = 1.1345743233455785456788e95; + + double basef_under = 1.1345743233455785456788e-27; + double based_under = 1.1345743233455785456788e-315; + double basef_over = 1.1345743233455785456788e+27; - double based_over = 1.1345743233455785456788e+185; - + 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, 1); - + std::cerr << "\nOverflow:\n"; check_accuracy(basef_over, based_over, s); + + std::cerr << "\nVarying (over):\n"; + for (int k=0; k<5; ++k) + { + check_accuracy_var(20,27,190,302,s); + std::cout << "\n"; + } } From 15ed32dd6ef53111e6e29428418a65e1d3547c12 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 16 Jul 2009 16:21:26 +0200 Subject: [PATCH 03/11] add other stable norm impl. in the benchmark --- bench/bench_norm.cpp | 75 +++++++++++++++++++++++++++++++++++--------- 1 file changed, 61 insertions(+), 14 deletions(-) diff --git a/bench/bench_norm.cpp b/bench/bench_norm.cpp index e06d06417..8e4fefdd7 100644 --- a/bench/bench_norm.cpp +++ b/bench/bench_norm.cpp @@ -1,5 +1,5 @@ #include -#include +#include #include "BenchTimer.h" using namespace Eigen; using namespace std; @@ -27,22 +27,54 @@ EIGEN_DONT_INLINE typename T::Scalar lapackNorm(T& v) { typedef typename T::Scalar Scalar; int n = v.size(); - Scalar scale = 1; - Scalar ssq = 0; + Scalar scale = 0; + Scalar ssq = 1; for (int i=0;i= ax) + { + ssq += ei_abs2(ax/scale); + } + else { ssq = Scalar(1) + ssq * ei_abs2(scale/ax); scale = ax; } - else - ssq += ei_abs2(ax/scale); } return scale * ei_sqrt(ssq); } +template +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 +EIGEN_DONT_INLINE typename T::Scalar bl2passNorm(T& v) +{ + const int blockSize = 4096; + typedef typename T::Scalar Scalar; + Scalar s = 0; + Scalar ssq = 0; + for (int bi=0; bi::type,Eigen::Dynamic,1,Eigen::ForceAligned> sv(v,bi,0,r,1); + Scalar m = sv.cwise().abs().maxCoeff(); + if (m>s) + { + ssq = ssq * ei_abs2(s/m); + s = m; + } + ssq += (sv/s).squaredNorm(); + } + return s*ei_sqrt(ssq); +} + template EIGEN_DONT_INLINE typename T::Scalar divacNorm(T& v) { @@ -210,6 +242,8 @@ void check_accuracy(double basef, double based, int s) 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) @@ -228,11 +262,13 @@ void check_accuracy_var(int ef0, int ef1, int ed0, int ed1, int s) std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\t" << blueNorm(vf.cast()) << "\t" << blueNorm(vd.cast()) << "\n"; std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\t" << blueNorm(vf.cast()) << "\t" << blueNorm(vd.cast()) << "\n"; std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\t" << lapackNorm(vf.cast()) << "\t" << lapackNorm(vd.cast()) << "\n"; + std::cout << "twopassNorm\t" << twopassNorm(vf) << "\t" << twopassNorm(vd) << "\t" << twopassNorm(vf.cast()) << "\t" << twopassNorm(vd.cast()) << "\n"; + std::cout << "bl2passNorm\t" << bl2passNorm(vf) << "\t" << bl2passNorm(vd) << "\t" << bl2passNorm(vf.cast()) << "\t" << bl2passNorm(vd.cast()) << "\n"; } int main(int argc, char** argv) { - int tries = 5; + int tries = 10; int iters = 100000; double y = 1.1345743233455785456788e12 * ei_random(); VectorXf v = VectorXf::Ones(1024) * y; @@ -240,33 +276,37 @@ int main(int argc, char** argv) std::cerr << "Performance (out of cache):\n"; { int iters = 1; - VectorXf vf = VectorXf::Ones(1024*1024*32) * y; - VectorXd vd = VectorXd::Ones(1024*1024*32) * y; + VectorXf vf = VectorXf::Random(1024*1024*32) * y; + VectorXd vd = VectorXd::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::Ones(512) * y; - VectorXd vd = VectorXd::Ones(512) * y; + VectorXf vf = VectorXf::Random(512) * y; + VectorXd vd = VectorXd::Random(512) * y; BENCH_PERF(sqsumNorm); BENCH_PERF(blueNorm); BENCH_PERF(pblueNorm); BENCH_PERF(lapackNorm); BENCH_PERF(hypotNorm); + BENCH_PERF(twopassNorm); + BENCH_PERF(bl2passNorm); } - +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-315; + double based_under = 1.1345743233455785456788e-303; double basef_over = 1.1345743233455785456788e+27; double based_over = 1.1345743233455785456788e+302; @@ -283,9 +323,16 @@ int main(int argc, char** argv) check_accuracy(basef_over, based_over, s); std::cerr << "\nVarying (over):\n"; - for (int k=0; k<5; ++k) + 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"; + } } From 32b08ac971b2ab66cf83360a9e2d42a99bfe3b70 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 17 Jul 2009 16:22:39 +0200 Subject: [PATCH 04/11] re-implement stableNorm using a homemade blocky and vectorization friendly algorithm (slow if no vectorization) --- Eigen/Core | 1 + Eigen/src/Core/Block.h | 4 +- Eigen/src/Core/Dot.h | 125 ----------------------- Eigen/src/Core/Functors.h | 3 - Eigen/src/Core/MatrixBase.h | 1 + Eigen/src/Core/StableNorm.h | 194 ++++++++++++++++++++++++++++++++++++ bench/bench_norm.cpp | 106 ++++++++++---------- test/adjoint.cpp | 4 +- 8 files changed, 256 insertions(+), 182 deletions(-) create mode 100644 Eigen/src/Core/StableNorm.h diff --git a/Eigen/Core b/Eigen/Core index a7a2a768a..ad606ec38 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -161,6 +161,7 @@ namespace Eigen { #include "src/Core/CwiseNullaryOp.h" #include "src/Core/CwiseUnaryView.h" #include "src/Core/Dot.h" +#include "src/Core/StableNorm.h" #include "src/Core/DiagonalProduct.h" #include "src/Core/SolveTriangular.h" #include "src/Core/MapBase.h" diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h index 382518696..ffc0707ee 100644 --- a/Eigen/src/Core/Block.h +++ b/Eigen/src/Core/Block.h @@ -33,8 +33,8 @@ * \param MatrixType the type of the object in which we are taking a block * \param BlockRows the number of rows of the block we are taking at compile time (optional) * \param BlockCols the number of columns of the block we are taking at compile time (optional) - * \param _PacketAccess allows to enforce aligned loads and stores if set to ForceAligned. - * The default is AsRequested. This parameter is internaly used by Eigen + * \param _PacketAccess allows to enforce aligned loads and stores if set to \b ForceAligned. + * The default is \b AsRequested. This parameter is internaly used by Eigen * in expressions such as \code mat.block() += other; \endcode and most of * the time this is the only way it is used. * \param _DirectAccessStatus \internal used for partial specialization diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index d97f5837f..9e84d72bb 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -292,131 +292,6 @@ inline typename NumTraits::Scalar>::Real MatrixBase< return ei_sqrt(squaredNorm()); } -/** \returns the \em l2 norm of \c *this using a numerically more stable - * algorithm. - * - * \sa norm(), dot(), squaredNorm(), blueNorm() - */ -template -inline typename NumTraits::Scalar>::Real -MatrixBase::stableNorm() const -{ - return this->cwise().abs().redux(ei_scalar_hypot_op()); -} - -/** \returns the \em l2 norm of \c *this using the Blue's algorithm. - * A Portable Fortran Program to Find the Euclidean Norm of a Vector, - * ACM TOMS, Vol 4, Issue 1, 1978. - * - * \sa norm(), dot(), squaredNorm(), stableNorm() - */ -template -inline typename NumTraits::Scalar>::Real -MatrixBase::blueNorm() const -{ - static int nmax = -1; - static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr; - int n; - Scalar ax, abig, amed, asml; - - if(nmax <= 0) - { - int nbig, ibeta, it, iemin, iemax, iexp; - Scalar abig, eps; - // This program calculates the machine-dependent constants - // bl, b2, slm, s2m, relerr overfl, nmax - // from the "basic" machine-dependent numbers - // nbig, ibeta, it, iemin, iemax, rbig. - // The following define the basic machine-dependent constants. - // For portability, the PORT subprograms "ilmaeh" and "rlmach" - // are used. For any specific computer, each of the assignment - // statements can be replaced - nbig = std::numeric_limits::max(); // largest integer - ibeta = std::numeric_limits::radix; //NumTraits::Base; // base for floating-point numbers - it = std::numeric_limits::digits; //NumTraits::Mantissa; // number of base-beta digits in mantissa - iemin = std::numeric_limits::min_exponent; // minimum exponent - iemax = std::numeric_limits::max_exponent; // maximum exponent - rbig = std::numeric_limits::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) - { - ei_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 = ei_sqrt(eps); // tolerance for neglecting asml - abig = 1.0/eps - 1.0; - if (Scalar(nbig)>abig) nmax = abig; // largest safe n - else nmax = nbig; - } - n = size(); - if(n==0) - return 0; - asml = Scalar(0); - amed = Scalar(0); - abig = Scalar(0); - for(int j=0; j b2) abig += ei_abs2(ax*s2m); - else if(ax < b1) asml += ei_abs2(ax*s1m); - else amed += ei_abs2(ax); - } - if(abig > Scalar(0)) - { - abig = ei_sqrt(abig); - if(abig > overfl) - { - ei_assert(false && "overflow"); - return rbig; - } - if(amed > Scalar(0)) - { - abig = abig/s2m; - amed = ei_sqrt(amed); - } - else - { - return abig/s2m; - } - - } - else if(asml > Scalar(0)) - { - if (amed > Scalar(0)) - { - abig = ei_sqrt(amed); - amed = ei_sqrt(asml) / s1m; - } - else - { - return ei_sqrt(asml)/s1m; - } - } - else - { - return ei_sqrt(amed); - } - asml = std::min(abig, amed); - abig = std::max(abig, amed); - if(asml <= abig*relerr) - return abig; - else - return abig * ei_sqrt(Scalar(1) + ei_abs2(asml/abig)); -} - /** \returns an expression of the quotient of *this by its own norm. * * \only_for_vectors diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index 89badb353..a4c9604df 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -124,9 +124,6 @@ template struct ei_scalar_hypot_op EIGEN_EMPTY_STRUCT { // typedef typename NumTraits::Real result_type; EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& _x, const Scalar& _y) const { -// typedef typename NumTraits::Real RealScalar; -// RealScalar _x = ei_abs(x); -// RealScalar _y = ei_abs(y); Scalar p = std::max(_x, _y); Scalar q = std::min(_x, _y); Scalar qp = q/p; diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index b8273ca22..5298ae271 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -384,6 +384,7 @@ template class MatrixBase RealScalar norm() const; RealScalar stableNorm() const; RealScalar blueNorm() const; + RealScalar hypotNorm() const; const PlainMatrixType normalized() const; void normalize(); diff --git a/Eigen/src/Core/StableNorm.h b/Eigen/src/Core/StableNorm.h new file mode 100644 index 000000000..383c64f01 --- /dev/null +++ b/Eigen/src/Core/StableNorm.h @@ -0,0 +1,194 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_STABLENORM_H +#define EIGEN_STABLENORM_H + +template +inline void ei_stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale) +{ + Scalar max = bl.cwise().abs().maxCoeff(); + if (max>scale) + { + ssq = ssq * ei_abs2(scale/max); + scale = max; + invScale = Scalar(1)/scale; + } + // TODO if the max is much much smaller than the current scale, + // then we can neglect this sub vector + ssq += (bl*invScale).squaredNorm(); +} + +/** \returns the \em l2 norm of \c *this avoiding underflow and overflow. + * This version use a blockwise two passes algorithm: + * 1 - find the absolute largest coefficient \c s + * 2 - compute \f$ s \Vert \frac{*this}{s} \Vert \f$ in a standard way + * + * For architecture/scalar types supporting vectorization, this version + * is faster than blueNorm(). Otherwise the blueNorm() is much faster. + * + * \sa norm(), blueNorm(), hypotNorm() + */ +template +inline typename NumTraits::Scalar>::Real +MatrixBase::stableNorm() const +{ + const int blockSize = 4096; + RealScalar scale = 0; + RealScalar invScale; + RealScalar ssq = 0; // sum of square + enum { + Alignment = (int(Flags)&DirectAccessBit) || (int(Flags)&AlignedBit) ? ForceAligned : AsRequested + }; + int n = size(); + int bi=0; + if ((int(Flags)&DirectAccessBit) && !(int(Flags)&AlignedBit)) + { + bi = ei_alignmentOffset(&const_cast_derived().coeffRef(0), n); + if (bi>0) + ei_stable_norm_kernel(start(bi), ssq, scale, invScale); + } + for (; bi(derived(),bi,std::min(blockSize, n - bi)), ssq, scale, invScale); + return scale * ei_sqrt(ssq); +} + +/** \returns the \em l2 norm of \c *this using the Blue's algorithm. + * A Portable Fortran Program to Find the Euclidean Norm of a Vector, + * ACM TOMS, Vol 4, Issue 1, 1978. + * + * For architecture/scalar types without vectorization, this version + * is much faster than stableNorm(). Otherwise the stableNorm() is faster. + * + * \sa norm(), stableNorm(), hypotNorm() + */ +template +inline typename NumTraits::Scalar>::Real +MatrixBase::blueNorm() const +{ + static int nmax = -1; + static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr; + if(nmax <= 0) + { + int nbig, ibeta, it, iemin, iemax, iexp; + RealScalar abig, eps; + // This program calculates the machine-dependent constants + // bl, b2, slm, s2m, relerr overfl, nmax + // from the "basic" machine-dependent numbers + // nbig, ibeta, it, iemin, iemax, rbig. + // The following define the basic machine-dependent constants. + // For portability, the PORT subprograms "ilmaeh" and "rlmach" + // are used. For any specific computer, each of the assignment + // statements can be replaced + nbig = std::numeric_limits::max(); // largest integer + ibeta = std::numeric_limits::radix; // base for floating-point numbers + it = std::numeric_limits::digits; // number of base-beta digits in mantissa + iemin = std::numeric_limits::min_exponent; // minimum exponent + iemax = std::numeric_limits::max_exponent; // maximum exponent + rbig = std::numeric_limits::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) + { + ei_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 = ei_sqrt(eps); // tolerance for neglecting asml + abig = 1.0/eps - 1.0; + if (RealScalar(nbig)>abig) nmax = abig; // largest safe n + else nmax = nbig; + } + int n = size(); + RealScalar ab2 = b2 / RealScalar(n); + RealScalar asml = RealScalar(0); + RealScalar amed = RealScalar(0); + RealScalar abig = RealScalar(0); + for(int j=0; j ab2) abig += ei_abs2(ax*s2m); + else if(ax < b1) asml += ei_abs2(ax*s1m); + else amed += ei_abs2(ax); + } + if(abig > RealScalar(0)) + { + abig = ei_sqrt(abig); + if(abig > overfl) + { + ei_assert(false && "overflow"); + return rbig; + } + if(amed > RealScalar(0)) + { + abig = abig/s2m; + amed = ei_sqrt(amed); + } + else + return abig/s2m; + } + else if(asml > RealScalar(0)) + { + if (amed > RealScalar(0)) + { + abig = ei_sqrt(amed); + amed = ei_sqrt(asml) / s1m; + } + else + return ei_sqrt(asml)/s1m; + } + else + return ei_sqrt(amed); + asml = std::min(abig, amed); + abig = std::max(abig, amed); + if(asml <= abig*relerr) + return abig; + else + return abig * ei_sqrt(RealScalar(1) + ei_abs2(asml/abig)); +} + +/** \returns the \em l2 norm of \c *this avoiding undeflow and overflow. + * This version use a concatenation of hypot() calls, and it is very slow. + * + * \sa norm(), stableNorm() + */ +template +inline typename NumTraits::Scalar>::Real +MatrixBase::hypotNorm() const +{ + return this->cwise().abs().redux(ei_scalar_hypot_op()); +} + +#endif // EIGEN_STABLENORM_H diff --git a/bench/bench_norm.cpp b/bench/bench_norm.cpp index 8e4fefdd7..7a3dc2e68 100644 --- a/bench/bench_norm.cpp +++ b/bench/bench_norm.cpp @@ -13,7 +13,7 @@ EIGEN_DONT_INLINE typename T::Scalar sqsumNorm(const T& v) template EIGEN_DONT_INLINE typename T::Scalar hypotNorm(const T& v) { - return v.stableNorm(); + return v.hypotNorm(); } template @@ -56,23 +56,7 @@ EIGEN_DONT_INLINE typename T::Scalar twopassNorm(T& v) template EIGEN_DONT_INLINE typename T::Scalar bl2passNorm(T& v) { - const int blockSize = 4096; - typedef typename T::Scalar Scalar; - Scalar s = 0; - Scalar ssq = 0; - for (int bi=0; bi::type,Eigen::Dynamic,1,Eigen::ForceAligned> sv(v,bi,0,r,1); - Scalar m = sv.cwise().abs().maxCoeff(); - if (m>s) - { - ssq = ssq * ei_abs2(s/m); - s = m; - } - ssq += (sv/s).squaredNorm(); - } - return s*ei_sqrt(ssq); + return v.stableNorm(); } template @@ -163,6 +147,19 @@ EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v) Packet ax_s1m = ei_pmul(ax,ps1m); Packet maskBig = ei_plt(pb2,ax); Packet maskSml = ei_plt(ax,pb1); + +// Packet maskMed = ei_pand(maskSml,maskBig); +// Packet scale = ei_pset1(Scalar(0)); +// scale = ei_por(scale, ei_pand(maskBig,ps2m)); +// scale = ei_por(scale, ei_pand(maskSml,ps1m)); +// scale = ei_por(scale, ei_pandnot(ei_pset1(Scalar(1)),maskMed)); +// ax = ei_pmul(ax,scale); +// ax = ei_pmul(ax,ax); +// pabig = ei_padd(pabig, ei_pand(maskBig, ax)); +// pasml = ei_padd(pasml, ei_pand(maskSml, ax)); +// pamed = ei_padd(pamed, ei_pandnot(ax,maskMed)); + + pabig = ei_padd(pabig, ei_pand(maskBig, ei_pmul(ax_s2m,ax_s2m))); pasml = ei_padd(pasml, ei_pand(maskSml, ei_pmul(ax_s1m,ax_s1m))); pamed = ei_padd(pamed, ei_pandnot(ei_pmul(ax,ax),ei_pand(maskSml,maskBig))); @@ -215,7 +212,7 @@ EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v) } #define BENCH_PERF(NRM) { \ - Eigen::BenchTimer tf, td; tf.reset(); td.reset();\ + Eigen::BenchTimer tf, td, tcf; tf.reset(); td.reset(); tcf.reset();\ for (int k=0; k()) << "\t" << blueNorm(vd.cast()) << "\n"; std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\t" << lapackNorm(vf.cast()) << "\t" << lapackNorm(vd.cast()) << "\n"; std::cout << "twopassNorm\t" << twopassNorm(vf) << "\t" << twopassNorm(vd) << "\t" << twopassNorm(vf.cast()) << "\t" << twopassNorm(vd.cast()) << "\n"; - std::cout << "bl2passNorm\t" << bl2passNorm(vf) << "\t" << bl2passNorm(vd) << "\t" << bl2passNorm(vf.cast()) << "\t" << bl2passNorm(vd.cast()) << "\n"; +// std::cout << "bl2passNorm\t" << bl2passNorm(vf) << "\t" << bl2passNorm(vd) << "\t" << bl2passNorm(vf.cast()) << "\t" << bl2passNorm(vd.cast()) << "\n"; } int main(int argc, char** argv) @@ -273,34 +275,7 @@ int main(int argc, char** argv) double y = 1.1345743233455785456788e12 * ei_random(); VectorXf v = VectorXf::Ones(1024) * y; - 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; - 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; - BENCH_PERF(sqsumNorm); - BENCH_PERF(blueNorm); - BENCH_PERF(pblueNorm); - BENCH_PERF(lapackNorm); - BENCH_PERF(hypotNorm); - BENCH_PERF(twopassNorm); - BENCH_PERF(bl2passNorm); - } -return 0; +// return 0; int s = 10000; double basef_ok = 1.1345743233455785456788e15; double based_ok = 1.1345743233455785456788e95; @@ -317,7 +292,7 @@ return 0; check_accuracy(basef_ok, based_ok, s); std::cerr << "\nUnderflow:\n"; - check_accuracy(basef_under, based_under, 1); + check_accuracy(basef_under, based_under, s); std::cerr << "\nOverflow:\n"; check_accuracy(basef_over, based_over, s); @@ -335,4 +310,35 @@ return 0; 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); + } } diff --git a/test/adjoint.cpp b/test/adjoint.cpp index 1f4aa7427..47e1bd740 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -76,8 +76,8 @@ template void adjoint(const MatrixType& m) { VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast(1)); VERIFY_IS_APPROX(v1.norm(), v1.stableNorm()); - // NOTE disabled because it currently compiles for float and double only - // VERIFY_IS_APPROX(v1.blueNorm(), v1.stableNorm()); + VERIFY_IS_APPROX(v1.blueNorm(), v1.stableNorm()); + VERIFY_IS_APPROX(v1.hypotNorm(), v1.stableNorm()); } // check compatibility of dot and adjoint From a551107ccea8fe027d2672cb82f6b70e741bb996 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 20 Jul 2009 10:35:47 +0200 Subject: [PATCH 05/11] bugfix for a = a * b; when a has to be resized --- Eigen/src/Core/Matrix.h | 21 +++++++++++++++++++-- test/product_large.cpp | 6 ++++++ 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 5301f4849..32b526ef3 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -536,6 +536,9 @@ class Matrix resizeLike(other); } + template + struct ei_matrix_set_selector; + /** \internal Copies the value of the expression \a other into \c *this with automatic resizing. * * *this might be resized to match the dimensions of \a other. If *this was a null matrix (not already initialized), @@ -550,8 +553,8 @@ class Matrix template EIGEN_STRONG_INLINE Matrix& _set(const MatrixBase& other) { - _resize_to_match(other); - return Base::operator=(other); + ei_matrix_set_selector::run(*this,other.derived()); + return *this; } /** \internal Like _set() but additionally makes the assumption that no aliasing effect can happen (which @@ -600,6 +603,20 @@ class Matrix } }; +template +template +struct Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::ei_matrix_set_selector +{ + static void run(MatrixType& dst, const OtherDerived& src) { dst._set_noalias(src.eval()); } +}; + +template +template +struct Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::ei_matrix_set_selector +{ + static void run(MatrixType& dst, const OtherDerived& src) { dst._set_noalias(src); } +}; + /** \defgroup matrixtypedefs Global matrix typedefs * * \ingroup Core_Module diff --git a/test/product_large.cpp b/test/product_large.cpp index 77ae7b587..9b53e7b89 100644 --- a/test/product_large.cpp +++ b/test/product_large.cpp @@ -42,4 +42,10 @@ void test_product_large() m = (v+v).asDiagonal() * m; VERIFY_IS_APPROX(m, MatrixXf::Constant(N,3,2)); } + + { + // test deferred resizing in Matrix::operator= + MatrixXf a = MatrixXf::Random(10,4), b = MatrixXf::Random(4,10), c = a; + VERIFY_IS_APPROX((a = a * b), (c * b).eval()); + } } From b3ad796d4087ca768452d1bbe0555a85cd7e1dde Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 20 Jul 2009 10:44:07 +0200 Subject: [PATCH 06/11] bugfix in operator*= (matrix product) --- Eigen/src/Core/Product.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 05a8221ee..52fb0db2c 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -300,7 +300,7 @@ template inline Derived & MatrixBase::operator*=(const MatrixBase &other) { - return *this = *this * other; + return derived() = derived() * other.derived(); } /*************************************************************************** From c10b919edb9a362388a98f620f8cef7292ff496a Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 20 Jul 2009 10:56:03 +0200 Subject: [PATCH 07/11] compilation fix --- Eigen/src/Core/Matrix.h | 25 +++++++------------------ 1 file changed, 7 insertions(+), 18 deletions(-) diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 32b526ef3..f41c9e042 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -536,9 +536,6 @@ class Matrix resizeLike(other); } - template - struct ei_matrix_set_selector; - /** \internal Copies the value of the expression \a other into \c *this with automatic resizing. * * *this might be resized to match the dimensions of \a other. If *this was a null matrix (not already initialized), @@ -553,10 +550,16 @@ class Matrix template EIGEN_STRONG_INLINE Matrix& _set(const MatrixBase& other) { - ei_matrix_set_selector::run(*this,other.derived()); + _set_selector(other.derived(), typename ei_meta_if<(int(OtherDerived::Flags) & EvalBeforeAssigningBit), ei_meta_true, ei_meta_false>::ret()); return *this; } + template + EIGEN_STRONG_INLINE void _set_selector(const OtherDerived& other, const ei_meta_true&) { _set_noalias(other.eval()); } + + template + EIGEN_STRONG_INLINE void _set_selector(const OtherDerived& other, const ei_meta_false&) { _set_noalias(other); } + /** \internal Like _set() but additionally makes the assumption that no aliasing effect can happen (which * is the case when creating a new matrix) so one can enforce lazy evaluation. * @@ -603,20 +606,6 @@ class Matrix } }; -template -template -struct Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::ei_matrix_set_selector -{ - static void run(MatrixType& dst, const OtherDerived& src) { dst._set_noalias(src.eval()); } -}; - -template -template -struct Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::ei_matrix_set_selector -{ - static void run(MatrixType& dst, const OtherDerived& src) { dst._set_noalias(src); } -}; - /** \defgroup matrixtypedefs Global matrix typedefs * * \ingroup Core_Module From 4c85fa8c734c8d103343c163fe95e9ffdeb9f5aa Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 20 Jul 2009 10:57:31 +0200 Subject: [PATCH 08/11] compilation fix (sun CC) --- Eigen/src/Core/StableNorm.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Eigen/src/Core/StableNorm.h b/Eigen/src/Core/StableNorm.h index 383c64f01..22809633d 100644 --- a/Eigen/src/Core/StableNorm.h +++ b/Eigen/src/Core/StableNorm.h @@ -115,17 +115,17 @@ MatrixBase::blueNorm() const ei_assert(false && "the algorithm cannot be guaranteed on this computer"); } iexp = -((1-iemin)/2); - b1 = std::pow(ibeta, iexp); // lower boundary of midrange + b1 = std::pow(double(ibeta),iexp); // lower boundary of midrange iexp = (iemax + 1 - it)/2; - b2 = std::pow(ibeta,iexp); // upper boundary of midrange + b2 = std::pow(double(ibeta),iexp); // upper boundary of midrange iexp = (2-iemin)/2; - s1m = std::pow(ibeta,iexp); // scaling factor for lower range + s1m = std::pow(double(ibeta),iexp); // scaling factor for lower range iexp = - ((iemax+it)/2); - s2m = std::pow(ibeta,iexp); // scaling factor for upper range + s2m = std::pow(double(ibeta),iexp); // scaling factor for upper range overfl = rbig*s2m; // overfow boundary for abig - eps = std::pow(ibeta, 1-it); + eps = std::pow(double(ibeta), 1-it); relerr = ei_sqrt(eps); // tolerance for neglecting asml abig = 1.0/eps - 1.0; if (RealScalar(nbig)>abig) nmax = abig; // largest safe n From 4375c043ace44973d077bbe7b8ea1fb8ebf419a6 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 20 Jul 2009 13:27:41 +0200 Subject: [PATCH 09/11] minor compilation fixes for Sun CC and ICC --- Eigen/QR | 2 +- Eigen/src/Cholesky/LDLT.h | 26 +++++++++++++------------- Eigen/src/Core/MapBase.h | 9 ++++++++- test/mixingtypes.cpp | 4 ++-- 4 files changed, 24 insertions(+), 17 deletions(-) diff --git a/Eigen/QR b/Eigen/QR index 97907d1e5..5f36d0987 100644 --- a/Eigen/QR +++ b/Eigen/QR @@ -41,7 +41,7 @@ namespace Eigen { // declare all classes for a given matrix type #define EIGEN_QR_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \ - PREFIX template class QR; \ + PREFIX template class HouseholderQR; \ PREFIX template class Tridiagonalization; \ PREFIX template class HessenbergDecomposition; \ PREFIX template class SelfAdjointEigenSolver diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 94660245e..a308b22d5 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -62,7 +62,7 @@ template class LDLT typedef Matrix IntColVectorType; typedef Matrix IntRowVectorType; - /** + /** * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to @@ -80,10 +80,10 @@ template class LDLT } /** \returns the lower triangular matrix L */ - inline Part matrixL(void) const - { + inline Part matrixL(void) const + { ei_assert(m_isInitialized && "LDLT is not initialized."); - return m_matrix; + return m_matrix; } /** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed, @@ -97,24 +97,24 @@ template class LDLT } /** \returns the coefficients of the diagonal matrix D */ - inline Diagonal vectorD(void) const - { + inline Diagonal vectorD(void) const + { ei_assert(m_isInitialized && "LDLT is not initialized."); - return m_matrix.diagonal(); + return m_matrix.diagonal(); } /** \returns true if the matrix is positive (semidefinite) */ - inline bool isPositive(void) const - { + inline bool isPositive(void) const + { ei_assert(m_isInitialized && "LDLT is not initialized."); - return m_sign == 1; + return m_sign == 1; } /** \returns true if the matrix is negative (semidefinite) */ - inline bool isNegative(void) const - { + inline bool isNegative(void) const + { ei_assert(m_isInitialized && "LDLT is not initialized."); - return m_sign == -1; + return m_sign == -1; } template diff --git a/Eigen/src/Core/MapBase.h b/Eigen/src/Core/MapBase.h index 59bf69ad6..0b7d815fb 100644 --- a/Eigen/src/Core/MapBase.h +++ b/Eigen/src/Core/MapBase.h @@ -178,7 +178,14 @@ template class MapBase } using Base::operator*=; - using Base::operator+=; + + template + Derived& operator+=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) + { return Base::operator+=(other); } + + template + Derived& operator-=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other) + { return Base::operator-=(other); } template Derived& operator+=(const MatrixBase& other) diff --git a/test/mixingtypes.cpp b/test/mixingtypes.cpp index ec8a0f190..d14232bd4 100644 --- a/test/mixingtypes.cpp +++ b/test/mixingtypes.cpp @@ -82,7 +82,7 @@ template void mixingtypes(int size = SizeAtCompileType) void test_mixingtypes() { // check that our operator new is indeed called: - CALL_SUBTEST(mixingtypes<3>()); - CALL_SUBTEST(mixingtypes<4>()); + CALL_SUBTEST(mixingtypes<3>(3)); + CALL_SUBTEST(mixingtypes<4>(4)); CALL_SUBTEST(mixingtypes(20)); } From a012aecbc43669123cf3abab2ceb1561cd5aaa92 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 20 Jul 2009 13:44:52 +0200 Subject: [PATCH 10/11] bugfix in SVD --- Eigen/src/SVD/SVD.h | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index f9f9feb89..b68d1c834 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -125,7 +125,7 @@ template class SVD { return (b >= Scalar(0.0) ? ei_abs(a) : -ei_abs(a)); } - + protected: /** \internal */ MatrixUType m_matU; @@ -254,11 +254,14 @@ void SVD::compute(const MatrixType& matrix) if (g != Scalar(0.0)) { g = Scalar(1.0)/g; - for (j=l; j Date: Mon, 20 Jul 2009 23:06:04 +0200 Subject: [PATCH 11/11] enable our own ctest dashboard --- CTestConfig.cmake | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/CTestConfig.cmake b/CTestConfig.cmake index 263cda480..0ccd5a1ad 100644 --- a/CTestConfig.cmake +++ b/CTestConfig.cmake @@ -5,9 +5,9 @@ ## ENABLE_TESTING() ## INCLUDE(CTest) set(CTEST_PROJECT_NAME "Eigen") -set(CTEST_NIGHTLY_START_TIME "05:00:00 UTC") +set(CTEST_NIGHTLY_START_TIME "06:00:00 UTC") set(CTEST_DROP_METHOD "http") -set(CTEST_DROP_SITE "my.cdash.org") -set(CTEST_DROP_LOCATION "/submit.php?project=Eigen") +set(CTEST_DROP_SITE "eigen.tuxfamily.org") +set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen") set(CTEST_DROP_SITE_CDASH TRUE)