diff options
author | Gael Guennebaud <g.gael@free.fr> | 2009-07-14 23:06:25 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2009-07-14 23:06:25 +0200 |
commit | 8120a5cecd982ed182d6e18b072fb9a95dc925dc (patch) | |
tree | 8d4f78cd34a78a63cd40ae454afe6cef3725a38c /Eigen/src/Core/Dot.h | |
parent | 279cedc1cec146cdaef6a93c734297797cdac2a8 (diff) | |
parent | f5d2317b12f3d6eea68db5fb48743ca1801d10d3 (diff) |
synch with main devel branch
Diffstat (limited to 'Eigen/src/Core/Dot.h')
-rw-r--r-- | Eigen/src/Core/Dot.h | 138 |
1 files changed, 137 insertions, 1 deletions
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index 3861984da..4f185ea5b 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -295,7 +295,7 @@ inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase< /** \returns the \em l2 norm of \c *this using a numerically more stable * algorithm. * - * \sa norm(), dot(), squaredNorm() + * \sa norm(), dot(), squaredNorm(), blueNorm() */ template<typename Derived> inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real @@ -304,6 +304,142 @@ MatrixBase<Derived>::stableNorm() const return this->cwise().abs().redux(ei_scalar_hypot_op<RealScalar>()); } +/** \internal Computes ibeta^iexp by binary expansion of iexp, + * exact if ibeta is the machine base */ +template<typename T> 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. + * + * \sa norm(), dot(), squaredNorm(), stableNorm() + */ +template<typename Derived> +inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real +MatrixBase<Derived>::blueNorm() const +{ + static int nmax; + 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<int>::max(); // largest integer + ibeta = NumTraits<Scalar>::Base; // base for floating-point numbers + it = 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) + { + ei_assert(false && "the algorithm cannot be guaranteed on this computer"); + } + iexp = -((1-iemin)/2); + b1 = bexp<Scalar>(ibeta, iexp); // lower boundary of midrange + iexp = (iemax + 1 - it)/2; + b2 = bexp<Scalar>(ibeta,iexp); // upper boundary of midrange + + iexp = (2-iemin)/2; + s1m = bexp<Scalar>(ibeta,iexp); // scaling factor for lower range + iexp = - ((iemax+it)/2); + s2m = bexp<Scalar>(ibeta,iexp); // scaling factor for upper range + + overfl = rbig*s2m; // overfow boundary for abig + eps = bexp<Scalar>(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<n; ++j) + { + ax = ei_abs(coeff(j)); + if(ax > b2) abig += ei_abs2(ax*s2m); + else if(ax < b2) 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 |