Eigen  3.3.3
StableNorm.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_STABLENORM_H
00011 #define EIGEN_STABLENORM_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 
00017 template<typename ExpressionType, typename Scalar>
00018 inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale)
00019 {
00020   Scalar maxCoeff = bl.cwiseAbs().maxCoeff();
00021   
00022   if(maxCoeff>scale)
00023   {
00024     ssq = ssq * numext::abs2(scale/maxCoeff);
00025     Scalar tmp = Scalar(1)/maxCoeff;
00026     if(tmp > NumTraits<Scalar>::highest())
00027     {
00028       invScale = NumTraits<Scalar>::highest();
00029       scale = Scalar(1)/invScale;
00030     }
00031     else if(maxCoeff>NumTraits<Scalar>::highest()) // we got a INF
00032     {
00033       invScale = Scalar(1);
00034       scale = maxCoeff;
00035     }
00036     else
00037     {
00038       scale = maxCoeff;
00039       invScale = tmp;
00040     }
00041   }
00042   else if(maxCoeff!=maxCoeff) // we got a NaN
00043   {
00044     scale = maxCoeff;
00045   }
00046   
00047   // TODO if the maxCoeff is much much smaller than the current scale,
00048   // then we can neglect this sub vector
00049   if(scale>Scalar(0)) // if scale==0, then bl is 0 
00050     ssq += (bl*invScale).squaredNorm();
00051 }
00052 
00053 template<typename Derived>
00054 inline typename NumTraits<typename traits<Derived>::Scalar>::Real
00055 blueNorm_impl(const EigenBase<Derived>& _vec)
00056 {
00057   typedef typename Derived::RealScalar RealScalar;  
00058   using std::pow;
00059   using std::sqrt;
00060   using std::abs;
00061   const Derived& vec(_vec.derived());
00062   static bool initialized = false;
00063   static RealScalar b1, b2, s1m, s2m, rbig, relerr;
00064   if(!initialized)
00065   {
00066     int ibeta, it, iemin, iemax, iexp;
00067     RealScalar eps;
00068     // This program calculates the machine-dependent constants
00069     // bl, b2, slm, s2m, relerr overfl
00070     // from the "basic" machine-dependent numbers
00071     // nbig, ibeta, it, iemin, iemax, rbig.
00072     // The following define the basic machine-dependent constants.
00073     // For portability, the PORT subprograms "ilmaeh" and "rlmach"
00074     // are used. For any specific computer, each of the assignment
00075     // statements can be replaced
00076     ibeta = std::numeric_limits<RealScalar>::radix;                 // base for floating-point numbers
00077     it    = std::numeric_limits<RealScalar>::digits;                // number of base-beta digits in mantissa
00078     iemin = std::numeric_limits<RealScalar>::min_exponent;          // minimum exponent
00079     iemax = std::numeric_limits<RealScalar>::max_exponent;          // maximum exponent
00080     rbig  = (std::numeric_limits<RealScalar>::max)();               // largest floating-point number
00081 
00082     iexp  = -((1-iemin)/2);
00083     b1    = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));    // lower boundary of midrange
00084     iexp  = (iemax + 1 - it)/2;
00085     b2    = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));    // upper boundary of midrange
00086 
00087     iexp  = (2-iemin)/2;
00088     s1m   = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));    // scaling factor for lower range
00089     iexp  = - ((iemax+it)/2);
00090     s2m   = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));    // scaling factor for upper range
00091 
00092     eps     = RealScalar(pow(double(ibeta), 1-it));
00093     relerr  = sqrt(eps);                                            // tolerance for neglecting asml
00094     initialized = true;
00095   }
00096   Index n = vec.size();
00097   RealScalar ab2 = b2 / RealScalar(n);
00098   RealScalar asml = RealScalar(0);
00099   RealScalar amed = RealScalar(0);
00100   RealScalar abig = RealScalar(0);
00101   for(typename Derived::InnerIterator it(vec, 0); it; ++it)
00102   {
00103     RealScalar ax = abs(it.value());
00104     if(ax > ab2)     abig += numext::abs2(ax*s2m);
00105     else if(ax < b1) asml += numext::abs2(ax*s1m);
00106     else             amed += numext::abs2(ax);
00107   }
00108   if(amed!=amed)
00109     return amed;  // we got a NaN
00110   if(abig > RealScalar(0))
00111   {
00112     abig = sqrt(abig);
00113     if(abig > rbig) // overflow, or *this contains INF values
00114       return abig;  // return INF
00115     if(amed > RealScalar(0))
00116     {
00117       abig = abig/s2m;
00118       amed = sqrt(amed);
00119     }
00120     else
00121       return abig/s2m;
00122   }
00123   else if(asml > RealScalar(0))
00124   {
00125     if (amed > RealScalar(0))
00126     {
00127       abig = sqrt(amed);
00128       amed = sqrt(asml) / s1m;
00129     }
00130     else
00131       return sqrt(asml)/s1m;
00132   }
00133   else
00134     return sqrt(amed);
00135   asml = numext::mini(abig, amed);
00136   abig = numext::maxi(abig, amed);
00137   if(asml <= abig*relerr)
00138     return abig;
00139   else
00140     return abig * sqrt(RealScalar(1) + numext::abs2(asml/abig));
00141 }
00142 
00143 } // end namespace internal
00144 
00155 template<typename Derived>
00156 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
00157 MatrixBase<Derived>::stableNorm() const
00158 {
00159   using std::sqrt;
00160   using std::abs;
00161   const Index blockSize = 4096;
00162   RealScalar scale(0);
00163   RealScalar invScale(1);
00164   RealScalar ssq(0); // sum of square
00165   
00166   typedef typename internal::nested_eval<Derived,2>::type DerivedCopy;
00167   typedef typename internal::remove_all<DerivedCopy>::type DerivedCopyClean;
00168   DerivedCopy copy(derived());
00169   
00170   enum {
00171     CanAlign = (   (int(DerivedCopyClean::Flags)&DirectAccessBit)
00172                 || (int(internal::evaluator<DerivedCopyClean>::Alignment)>0) // FIXME Alignment)>0 might not be enough
00173                ) && (blockSize*sizeof(Scalar)*2<EIGEN_STACK_ALLOCATION_LIMIT) // ifwe cannot allocate on the stack, then let's not bother about this optimization
00174   };
00175   typedef typename internal::conditional<CanAlign, Ref<const Matrix<Scalar,Dynamic,1,0,blockSize,1>, internal::evaluator<DerivedCopyClean>::Alignment>,
00176                                                    typename DerivedCopyClean::ConstSegmentReturnType>::type SegmentWrapper;
00177   Index n = size();
00178   
00179   if(n==1)
00180     return abs(this->coeff(0));
00181   
00182   Index bi = internal::first_default_aligned(copy);
00183   if (bi>0)
00184     internal::stable_norm_kernel(copy.head(bi), ssq, scale, invScale);
00185   for (; bi<n; bi+=blockSize)
00186     internal::stable_norm_kernel(SegmentWrapper(copy.segment(bi,numext::mini(blockSize, n - bi))), ssq, scale, invScale);
00187   return scale * sqrt(ssq);
00188 }
00189 
00199 template<typename Derived>
00200 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
00201 MatrixBase<Derived>::blueNorm() const
00202 {
00203   return internal::blueNorm_impl(*this);
00204 }
00205 
00211 template<typename Derived>
00212 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
00213 MatrixBase<Derived>::hypotNorm() const
00214 {
00215   return this->cwiseAbs().redux(internal::scalar_hypot_op<RealScalar>());
00216 }
00217 
00218 } // end namespace Eigen
00219 
00220 #endif // EIGEN_STABLENORM_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends