1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_STABLENORM_H 11 #define EIGEN_STABLENORM_H 12 13 namespace Eigen { 14 15 namespace internal { 16 17 template<typename ExpressionType, typename Scalar> 18 inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale) 19 { 20 using std::max; 21 Scalar maxCoeff = bl.cwiseAbs().maxCoeff(); 22 23 if (maxCoeff>scale) 24 { 25 ssq = ssq * numext::abs2(scale/maxCoeff); 26 Scalar tmp = Scalar(1)/maxCoeff; 27 if(tmp > NumTraits<Scalar>::highest()) 28 { 29 invScale = NumTraits<Scalar>::highest(); 30 scale = Scalar(1)/invScale; 31 } 32 else 33 { 34 scale = maxCoeff; 35 invScale = tmp; 36 } 37 } 38 39 // TODO if the maxCoeff is much much smaller than the current scale, 40 // then we can neglect this sub vector 41 if(scale>Scalar(0)) // if scale==0, then bl is 0 42 ssq += (bl*invScale).squaredNorm(); 43 } 44 45 template<typename Derived> 46 inline typename NumTraits<typename traits<Derived>::Scalar>::Real 47 blueNorm_impl(const EigenBase<Derived>& _vec) 48 { 49 typedef typename Derived::RealScalar RealScalar; 50 typedef typename Derived::Index Index; 51 using std::pow; 52 using std::min; 53 using std::max; 54 using std::sqrt; 55 using std::abs; 56 const Derived& vec(_vec.derived()); 57 static bool initialized = false; 58 static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr; 59 if(!initialized) 60 { 61 int ibeta, it, iemin, iemax, iexp; 62 RealScalar eps; 63 // This program calculates the machine-dependent constants 64 // bl, b2, slm, s2m, relerr overfl 65 // from the "basic" machine-dependent numbers 66 // nbig, ibeta, it, iemin, iemax, rbig. 67 // The following define the basic machine-dependent constants. 68 // For portability, the PORT subprograms "ilmaeh" and "rlmach" 69 // are used. For any specific computer, each of the assignment 70 // statements can be replaced 71 ibeta = std::numeric_limits<RealScalar>::radix; // base for floating-point numbers 72 it = std::numeric_limits<RealScalar>::digits; // number of base-beta digits in mantissa 73 iemin = std::numeric_limits<RealScalar>::min_exponent; // minimum exponent 74 iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent 75 rbig = (std::numeric_limits<RealScalar>::max)(); // largest floating-point number 76 77 iexp = -((1-iemin)/2); 78 b1 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // lower boundary of midrange 79 iexp = (iemax + 1 - it)/2; 80 b2 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // upper boundary of midrange 81 82 iexp = (2-iemin)/2; 83 s1m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for lower range 84 iexp = - ((iemax+it)/2); 85 s2m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for upper range 86 87 overfl = rbig*s2m; // overflow boundary for abig 88 eps = RealScalar(pow(double(ibeta), 1-it)); 89 relerr = sqrt(eps); // tolerance for neglecting asml 90 initialized = true; 91 } 92 Index n = vec.size(); 93 RealScalar ab2 = b2 / RealScalar(n); 94 RealScalar asml = RealScalar(0); 95 RealScalar amed = RealScalar(0); 96 RealScalar abig = RealScalar(0); 97 for(typename Derived::InnerIterator it(vec, 0); it; ++it) 98 { 99 RealScalar ax = abs(it.value()); 100 if(ax > ab2) abig += numext::abs2(ax*s2m); 101 else if(ax < b1) asml += numext::abs2(ax*s1m); 102 else amed += numext::abs2(ax); 103 } 104 if(abig > RealScalar(0)) 105 { 106 abig = sqrt(abig); 107 if(abig > overfl) 108 { 109 return rbig; 110 } 111 if(amed > RealScalar(0)) 112 { 113 abig = abig/s2m; 114 amed = sqrt(amed); 115 } 116 else 117 return abig/s2m; 118 } 119 else if(asml > RealScalar(0)) 120 { 121 if (amed > RealScalar(0)) 122 { 123 abig = sqrt(amed); 124 amed = sqrt(asml) / s1m; 125 } 126 else 127 return sqrt(asml)/s1m; 128 } 129 else 130 return sqrt(amed); 131 asml = (min)(abig, amed); 132 abig = (max)(abig, amed); 133 if(asml <= abig*relerr) 134 return abig; 135 else 136 return abig * sqrt(RealScalar(1) + numext::abs2(asml/abig)); 137 } 138 139 } // end namespace internal 140 141 /** \returns the \em l2 norm of \c *this avoiding underflow and overflow. 142 * This version use a blockwise two passes algorithm: 143 * 1 - find the absolute largest coefficient \c s 144 * 2 - compute \f$ s \Vert \frac{*this}{s} \Vert \f$ in a standard way 145 * 146 * For architecture/scalar types supporting vectorization, this version 147 * is faster than blueNorm(). Otherwise the blueNorm() is much faster. 148 * 149 * \sa norm(), blueNorm(), hypotNorm() 150 */ 151 template<typename Derived> 152 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real 153 MatrixBase<Derived>::stableNorm() const 154 { 155 using std::min; 156 using std::sqrt; 157 const Index blockSize = 4096; 158 RealScalar scale(0); 159 RealScalar invScale(1); 160 RealScalar ssq(0); // sum of square 161 enum { 162 Alignment = (int(Flags)&DirectAccessBit) || (int(Flags)&AlignedBit) ? 1 : 0 163 }; 164 Index n = size(); 165 Index bi = internal::first_aligned(derived()); 166 if (bi>0) 167 internal::stable_norm_kernel(this->head(bi), ssq, scale, invScale); 168 for (; bi<n; bi+=blockSize) 169 internal::stable_norm_kernel(this->segment(bi,(min)(blockSize, n - bi)).template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale); 170 return scale * sqrt(ssq); 171 } 172 173 /** \returns the \em l2 norm of \c *this using the Blue's algorithm. 174 * A Portable Fortran Program to Find the Euclidean Norm of a Vector, 175 * ACM TOMS, Vol 4, Issue 1, 1978. 176 * 177 * For architecture/scalar types without vectorization, this version 178 * is much faster than stableNorm(). Otherwise the stableNorm() is faster. 179 * 180 * \sa norm(), stableNorm(), hypotNorm() 181 */ 182 template<typename Derived> 183 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real 184 MatrixBase<Derived>::blueNorm() const 185 { 186 return internal::blueNorm_impl(*this); 187 } 188 189 /** \returns the \em l2 norm of \c *this avoiding undeflow and overflow. 190 * This version use a concatenation of hypot() calls, and it is very slow. 191 * 192 * \sa norm(), stableNorm() 193 */ 194 template<typename Derived> 195 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real 196 MatrixBase<Derived>::hypotNorm() const 197 { 198 return this->cwiseAbs().redux(internal::scalar_hypot_op<RealScalar>()); 199 } 200 201 } // end namespace Eigen 202 203 #endif // EIGEN_STABLENORM_H 204
