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