StableNorm.h
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>
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 
42 template<typename Derived>
43 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
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 
72 template<typename Derived>
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  return rbig;
135  }
136  if(amed > RealScalar(0))
137  {
138  abig = abig/s2m;
139  amed = internal::sqrt(amed);
140  }
141  else
142  return abig/s2m;
143  }
144  else if(asml > RealScalar(0))
145  {
146  if (amed > RealScalar(0))
147  {
148  abig = internal::sqrt(amed);
149  amed = internal::sqrt(asml) / s1m;
150  }
151  else
152  return internal::sqrt(asml)/s1m;
153  }
154  else
155  return internal::sqrt(amed);
156  asml = (min)(abig, amed);
157  abig = (max)(abig, amed);
158  if(asml <= abig*relerr)
159  return abig;
160  else
161  return abig * internal::sqrt(RealScalar(1) + internal::abs2(asml/abig));
162 }
163 
169 template<typename Derived>
172 {
173  return this->cwiseAbs().redux(internal::scalar_hypot_op<RealScalar>());
174 }
175 
176 } // end namespace Eigen
177 
178 #endif // EIGEN_STABLENORM_H