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@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>
stable_norm_kernel(const ExpressionType & bl,Scalar & ssq,Scalar & scale,Scalar & invScale)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
stableNorm()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
blueNorm()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
hypotNorm()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