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 // 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  Scalar max = bl.cwiseAbs().maxCoeff();
21  if (max>scale)
22  {
23  ssq = ssq * numext::abs2(scale/max);
24  scale = max;
25  invScale = Scalar(1)/scale;
26  }
27  // TODO if the max is much much smaller than the current scale,
28  // then we can neglect this sub vector
29  ssq += (bl*invScale).squaredNorm();
30 }
31 
32 template<typename Derived>
33 inline typename NumTraits<typename traits<Derived>::Scalar>::Real
35 {
36  typedef typename Derived::RealScalar RealScalar;
37  typedef typename Derived::Index Index;
38  using std::pow;
39  using std::min;
40  using std::max;
41  using std::sqrt;
42  using std::abs;
43  const Derived& vec(_vec.derived());
44  static bool initialized = false;
45  static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
46  if(!initialized)
47  {
48  int ibeta, it, iemin, iemax, iexp;
49  RealScalar eps;
50  // This program calculates the machine-dependent constants
51  // bl, b2, slm, s2m, relerr overfl
52  // from the "basic" machine-dependent numbers
53  // nbig, ibeta, it, iemin, iemax, rbig.
54  // The following define the basic machine-dependent constants.
55  // For portability, the PORT subprograms "ilmaeh" and "rlmach"
56  // are used. For any specific computer, each of the assignment
57  // statements can be replaced
58  ibeta = std::numeric_limits<RealScalar>::radix; // base for floating-point numbers
59  it = std::numeric_limits<RealScalar>::digits; // number of base-beta digits in mantissa
60  iemin = std::numeric_limits<RealScalar>::min_exponent; // minimum exponent
61  iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent
62  rbig = (std::numeric_limits<RealScalar>::max)(); // largest floating-point number
63 
64  iexp = -((1-iemin)/2);
65  b1 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // lower boundary of midrange
66  iexp = (iemax + 1 - it)/2;
67  b2 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // upper boundary of midrange
68 
69  iexp = (2-iemin)/2;
70  s1m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for lower range
71  iexp = - ((iemax+it)/2);
72  s2m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for upper range
73 
74  overfl = rbig*s2m; // overflow boundary for abig
75  eps = RealScalar(pow(double(ibeta), 1-it));
76  relerr = sqrt(eps); // tolerance for neglecting asml
77  initialized = true;
78  }
79  Index n = vec.size();
80  RealScalar ab2 = b2 / RealScalar(n);
81  RealScalar asml = RealScalar(0);
82  RealScalar amed = RealScalar(0);
83  RealScalar abig = RealScalar(0);
84  for(typename Derived::InnerIterator it(vec, 0); it; ++it)
85  {
86  RealScalar ax = abs(it.value());
87  if(ax > ab2) abig += numext::abs2(ax*s2m);
88  else if(ax < b1) asml += numext::abs2(ax*s1m);
89  else amed += numext::abs2(ax);
90  }
91  if(abig > RealScalar(0))
92  {
93  abig = sqrt(abig);
94  if(abig > overfl)
95  {
96  return rbig;
97  }
98  if(amed > RealScalar(0))
99  {
100  abig = abig/s2m;
101  amed = sqrt(amed);
102  }
103  else
104  return abig/s2m;
105  }
106  else if(asml > RealScalar(0))
107  {
108  if (amed > RealScalar(0))
109  {
110  abig = sqrt(amed);
111  amed = sqrt(asml) / s1m;
112  }
113  else
114  return sqrt(asml)/s1m;
115  }
116  else
117  return sqrt(amed);
118  asml = (min)(abig, amed);
119  abig = (max)(abig, amed);
120  if(asml <= abig*relerr)
121  return abig;
122  else
123  return abig * sqrt(RealScalar(1) + numext::abs2(asml/abig));
124 }
125 
126 } // end namespace internal
127 
138 template<typename Derived>
141 {
142  using std::min;
143  using std::sqrt;
144  const Index blockSize = 4096;
145  RealScalar scale(0);
146  RealScalar invScale(1);
147  RealScalar ssq(0); // sum of square
148  enum {
149  Alignment = (int(Flags)&DirectAccessBit) || (int(Flags)&AlignedBit) ? 1 : 0
150  };
151  Index n = size();
152  Index bi = internal::first_aligned(derived());
153  if (bi>0)
154  internal::stable_norm_kernel(this->head(bi), ssq, scale, invScale);
155  for (; bi<n; bi+=blockSize)
156  internal::stable_norm_kernel(this->segment(bi,(min)(blockSize, n - bi)).template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
157  return scale * sqrt(ssq);
158 }
159 
169 template<typename Derived>
172 {
173  return internal::blueNorm_impl(*this);
174 }
175 
181 template<typename Derived>
184 {
185  return this->cwiseAbs().redux(internal::scalar_hypot_op<RealScalar>());
186 }
187 
188 } // end namespace Eigen
189 
190 #endif // EIGEN_STABLENORM_H
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar >, const Derived > pow(const Eigen::ArrayBase< Derived > &x, const typename Derived::Scalar &exponent)
IntermediateState sqrt(const Expression &arg)
internal::traits< Derived >::Index Index
The type of indices.
Definition: DenseBase.h:61
const unsigned int DirectAccessBit
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
void stable_norm_kernel(const ExpressionType &bl, Scalar &ssq, Scalar &scale, Scalar &invScale)
Definition: StableNorm.h:18
SegmentReturnType segment(Index start, Index vecSize)
Definition: BlockMethods.h:752
IntermediateState pow(const Expression &arg1, const Expression &arg2)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
Derived & derived()
Definition: EigenBase.h:34
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > cwiseAbs() const
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs2_op< Scalar >, const Derived > abs2() const
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
const unsigned int AlignedBit
NumTraits< typename traits< Derived >::Scalar >::Real blueNorm_impl(const EigenBase< Derived > &_vec)
Definition: StableNorm.h:34
RealScalar blueNorm() const
Definition: StableNorm.h:171
RealScalar stableNorm() const
Definition: StableNorm.h:140
SegmentReturnType head(Index vecSize)
Definition: BlockMethods.h:781
NumTraits< Scalar >::Real RealScalar
Definition: DenseBase.h:65
RealScalar hypotNorm() const
Definition: StableNorm.h:183
static Derived::Index first_aligned(const Derived &m)


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:11