StableNorm.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_STABLENORM_H
26 #define EIGEN_STABLENORM_H
27 
28 namespace Eigen {
29 
30 namespace internal {
31 template<typename ExpressionType, typename Scalar>
32 inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale)
33 {
34  Scalar max = bl.cwiseAbs().maxCoeff();
35  if (max>scale)
36  {
37  ssq = ssq * abs2(scale/max);
38  scale = max;
39  invScale = Scalar(1)/scale;
40  }
41  // TODO if the max is much much smaller than the current scale,
42  // then we can neglect this sub vector
43  ssq += (bl*invScale).squaredNorm();
44 }
45 }
46 
57 template<typename Derived>
58 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
60 {
61  using std::min;
62  const Index blockSize = 4096;
63  RealScalar scale(0);
64  RealScalar invScale(1);
65  RealScalar ssq(0); // sum of square
66  enum {
67  Alignment = (int(Flags)&DirectAccessBit) || (int(Flags)&AlignedBit) ? 1 : 0
68  };
69  Index n = size();
70  Index bi = internal::first_aligned(derived());
71  if (bi>0)
72  internal::stable_norm_kernel(this->head(bi), ssq, scale, invScale);
73  for (; bi<n; bi+=blockSize)
74  internal::stable_norm_kernel(this->segment(bi,(min)(blockSize, n - bi)).template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
75  return scale * internal::sqrt(ssq);
76 }
77 
87 template<typename Derived>
90 {
91  using std::pow;
92  using std::min;
93  using std::max;
94  static Index nmax = -1;
95  static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
96  if(nmax <= 0)
97  {
98  int nbig, ibeta, it, iemin, iemax, iexp;
99  RealScalar abig, eps;
100  // This program calculates the machine-dependent constants
101  // bl, b2, slm, s2m, relerr overfl, nmax
102  // from the "basic" machine-dependent numbers
103  // nbig, ibeta, it, iemin, iemax, rbig.
104  // The following define the basic machine-dependent constants.
105  // For portability, the PORT subprograms "ilmaeh" and "rlmach"
106  // are used. For any specific computer, each of the assignment
107  // statements can be replaced
108  nbig = (std::numeric_limits<Index>::max)(); // largest integer
109  ibeta = std::numeric_limits<RealScalar>::radix; // base for floating-point numbers
110  it = std::numeric_limits<RealScalar>::digits; // number of base-beta digits in mantissa
111  iemin = std::numeric_limits<RealScalar>::min_exponent; // minimum exponent
112  iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent
113  rbig = (std::numeric_limits<RealScalar>::max)(); // largest floating-point number
114 
115  iexp = -((1-iemin)/2);
116  b1 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // lower boundary of midrange
117  iexp = (iemax + 1 - it)/2;
118  b2 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // upper boundary of midrange
119 
120  iexp = (2-iemin)/2;
121  s1m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for lower range
122  iexp = - ((iemax+it)/2);
123  s2m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for upper range
124 
125  overfl = rbig*s2m; // overflow boundary for abig
126  eps = RealScalar(pow(double(ibeta), 1-it));
127  relerr = internal::sqrt(eps); // tolerance for neglecting asml
128  abig = RealScalar(1.0/eps - 1.0);
129  if (RealScalar(nbig)>abig) nmax = int(abig); // largest safe n
130  else nmax = nbig;
131  }
132  Index n = size();
133  RealScalar ab2 = b2 / RealScalar(n);
134  RealScalar asml = RealScalar(0);
135  RealScalar amed = RealScalar(0);
136  RealScalar abig = RealScalar(0);
137  for(Index j=0; j<n; ++j)
138  {
139  RealScalar ax = internal::abs(coeff(j));
140  if(ax > ab2) abig += internal::abs2(ax*s2m);
141  else if(ax < b1) asml += internal::abs2(ax*s1m);
142  else amed += internal::abs2(ax);
143  }
144  if(abig > RealScalar(0))
145  {
146  abig = internal::sqrt(abig);
147  if(abig > overfl)
148  {
149  eigen_assert(false && "overflow");
150  return rbig;
151  }
152  if(amed > RealScalar(0))
153  {
154  abig = abig/s2m;
155  amed = internal::sqrt(amed);
156  }
157  else
158  return abig/s2m;
159  }
160  else if(asml > RealScalar(0))
161  {
162  if (amed > RealScalar(0))
163  {
164  abig = internal::sqrt(amed);
165  amed = internal::sqrt(asml) / s1m;
166  }
167  else
168  return internal::sqrt(asml)/s1m;
169  }
170  else
171  return internal::sqrt(amed);
172  asml = (min)(abig, amed);
173  abig = (max)(abig, amed);
174  if(asml <= abig*relerr)
175  return abig;
176  else
177  return abig * internal::sqrt(RealScalar(1) + internal::abs2(asml/abig));
178 }
179 
185 template<typename Derived>
188 {
189  return this->cwiseAbs().redux(internal::scalar_hypot_op<RealScalar>());
190 }
191 
192 } // end namespace Eigen
193 
194 #endif // EIGEN_STABLENORM_H