stable_norm.cpp
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-2014 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 #include "main.h"
11 
12 template<typename T> EIGEN_DONT_INLINE T copy(const T& x)
13 {
14  return x;
15 }
16 
17 template<typename MatrixType> void stable_norm(const MatrixType& m)
18 {
19  /* this test covers the following files:
20  StableNorm.h
21  */
22  using std::sqrt;
23  using std::abs;
24  typedef typename MatrixType::Scalar Scalar;
25  typedef typename NumTraits<Scalar>::Real RealScalar;
26 
27  bool complex_real_product_ok = true;
28 
29  // Check the basic machine-dependent constants.
30  {
31  int ibeta, it, iemin, iemax;
32 
33  ibeta = std::numeric_limits<RealScalar>::radix; // base for floating-point numbers
34  it = std::numeric_limits<RealScalar>::digits; // number of base-beta digits in mantissa
35  iemin = std::numeric_limits<RealScalar>::min_exponent; // minimum exponent
36  iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent
37 
38  VERIFY( (!(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5) || (it<=4 && ibeta <= 3 ) || it<2))
39  && "the stable norm algorithm cannot be guaranteed on this computer");
40 
41  Scalar inf = std::numeric_limits<RealScalar>::infinity();
42  if(NumTraits<Scalar>::IsComplex && (numext::isnan)(inf*RealScalar(1)) )
43  {
44  complex_real_product_ok = false;
45  static bool first = true;
46  if(first)
47  std::cerr << "WARNING: compiler mess up complex*real product, " << inf << " * " << 1.0 << " = " << inf*RealScalar(1) << std::endl;
48  first = false;
49  }
50  }
51 
52 
53  Index rows = m.rows();
54  Index cols = m.cols();
55 
56  // get a non-zero random factor
57  Scalar factor = internal::random<Scalar>();
58  while(numext::abs2(factor)<RealScalar(1e-4))
59  factor = internal::random<Scalar>();
60  Scalar big = factor * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
61 
62  factor = internal::random<Scalar>();
63  while(numext::abs2(factor)<RealScalar(1e-4))
64  factor = internal::random<Scalar>();
65  Scalar small = factor * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4));
66 
67  Scalar one(1);
68 
69  MatrixType vzero = MatrixType::Zero(rows, cols),
70  vrand = MatrixType::Random(rows, cols),
71  vbig(rows, cols),
72  vsmall(rows,cols);
73 
74  vbig.fill(big);
75  vsmall.fill(small);
76 
77  VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1));
78  VERIFY_IS_APPROX(vrand.stableNorm(), vrand.norm());
79  VERIFY_IS_APPROX(vrand.blueNorm(), vrand.norm());
80  VERIFY_IS_APPROX(vrand.hypotNorm(), vrand.norm());
81 
82  // test with expressions as input
83  VERIFY_IS_APPROX((one*vrand).stableNorm(), vrand.norm());
84  VERIFY_IS_APPROX((one*vrand).blueNorm(), vrand.norm());
85  VERIFY_IS_APPROX((one*vrand).hypotNorm(), vrand.norm());
86  VERIFY_IS_APPROX((one*vrand+one*vrand-one*vrand).stableNorm(), vrand.norm());
87  VERIFY_IS_APPROX((one*vrand+one*vrand-one*vrand).blueNorm(), vrand.norm());
88  VERIFY_IS_APPROX((one*vrand+one*vrand-one*vrand).hypotNorm(), vrand.norm());
89 
90  RealScalar size = static_cast<RealScalar>(m.size());
91 
92  // test numext::isfinite
93  VERIFY(!(numext::isfinite)( std::numeric_limits<RealScalar>::infinity()));
94  VERIFY(!(numext::isfinite)(sqrt(-abs(big))));
95 
96  // test overflow
97  VERIFY((numext::isfinite)(sqrt(size)*abs(big)));
98  VERIFY_IS_NOT_APPROX(sqrt(copy(vbig.squaredNorm())), abs(sqrt(size)*big)); // here the default norm must fail
99  VERIFY_IS_APPROX(vbig.stableNorm(), sqrt(size)*abs(big));
100  VERIFY_IS_APPROX(vbig.blueNorm(), sqrt(size)*abs(big));
101  VERIFY_IS_APPROX(vbig.hypotNorm(), sqrt(size)*abs(big));
102 
103  // test underflow
104  VERIFY((numext::isfinite)(sqrt(size)*abs(small)));
105  VERIFY_IS_NOT_APPROX(sqrt(copy(vsmall.squaredNorm())), abs(sqrt(size)*small)); // here the default norm must fail
106  VERIFY_IS_APPROX(vsmall.stableNorm(), sqrt(size)*abs(small));
107  VERIFY_IS_APPROX(vsmall.blueNorm(), sqrt(size)*abs(small));
108  VERIFY_IS_APPROX(vsmall.hypotNorm(), sqrt(size)*abs(small));
109 
110  // Test compilation of cwise() version
111  VERIFY_IS_APPROX(vrand.colwise().stableNorm(), vrand.colwise().norm());
112  VERIFY_IS_APPROX(vrand.colwise().blueNorm(), vrand.colwise().norm());
113  VERIFY_IS_APPROX(vrand.colwise().hypotNorm(), vrand.colwise().norm());
114  VERIFY_IS_APPROX(vrand.rowwise().stableNorm(), vrand.rowwise().norm());
115  VERIFY_IS_APPROX(vrand.rowwise().blueNorm(), vrand.rowwise().norm());
116  VERIFY_IS_APPROX(vrand.rowwise().hypotNorm(), vrand.rowwise().norm());
117 
118  // test NaN, +inf, -inf
119  MatrixType v;
120  Index i = internal::random<Index>(0,rows-1);
121  Index j = internal::random<Index>(0,cols-1);
122 
123  // NaN
124  {
125  v = vrand;
126  v(i,j) = std::numeric_limits<RealScalar>::quiet_NaN();
127  VERIFY(!(numext::isfinite)(v.squaredNorm())); VERIFY((numext::isnan)(v.squaredNorm()));
128  VERIFY(!(numext::isfinite)(v.norm())); VERIFY((numext::isnan)(v.norm()));
129  VERIFY(!(numext::isfinite)(v.stableNorm())); VERIFY((numext::isnan)(v.stableNorm()));
130  VERIFY(!(numext::isfinite)(v.blueNorm())); VERIFY((numext::isnan)(v.blueNorm()));
131  VERIFY(!(numext::isfinite)(v.hypotNorm())); VERIFY((numext::isnan)(v.hypotNorm()));
132  }
133 
134  // +inf
135  {
136  v = vrand;
137  v(i,j) = std::numeric_limits<RealScalar>::infinity();
138  VERIFY(!(numext::isfinite)(v.squaredNorm())); VERIFY(isPlusInf(v.squaredNorm()));
139  VERIFY(!(numext::isfinite)(v.norm())); VERIFY(isPlusInf(v.norm()));
140  VERIFY(!(numext::isfinite)(v.stableNorm()));
141  if(complex_real_product_ok){
142  VERIFY(isPlusInf(v.stableNorm()));
143  }
144  VERIFY(!(numext::isfinite)(v.blueNorm())); VERIFY(isPlusInf(v.blueNorm()));
145  VERIFY(!(numext::isfinite)(v.hypotNorm())); VERIFY(isPlusInf(v.hypotNorm()));
146  }
147 
148  // -inf
149  {
150  v = vrand;
151  v(i,j) = -std::numeric_limits<RealScalar>::infinity();
152  VERIFY(!(numext::isfinite)(v.squaredNorm())); VERIFY(isPlusInf(v.squaredNorm()));
153  VERIFY(!(numext::isfinite)(v.norm())); VERIFY(isPlusInf(v.norm()));
154  VERIFY(!(numext::isfinite)(v.stableNorm()));
155  if(complex_real_product_ok) {
156  VERIFY(isPlusInf(v.stableNorm()));
157  }
158  VERIFY(!(numext::isfinite)(v.blueNorm())); VERIFY(isPlusInf(v.blueNorm()));
159  VERIFY(!(numext::isfinite)(v.hypotNorm())); VERIFY(isPlusInf(v.hypotNorm()));
160  }
161 
162  // mix
163  {
164  Index i2 = internal::random<Index>(0,rows-1);
165  Index j2 = internal::random<Index>(0,cols-1);
166  v = vrand;
167  v(i,j) = -std::numeric_limits<RealScalar>::infinity();
168  v(i2,j2) = std::numeric_limits<RealScalar>::quiet_NaN();
169  VERIFY(!(numext::isfinite)(v.squaredNorm())); VERIFY((numext::isnan)(v.squaredNorm()));
170  VERIFY(!(numext::isfinite)(v.norm())); VERIFY((numext::isnan)(v.norm()));
171  VERIFY(!(numext::isfinite)(v.stableNorm())); VERIFY((numext::isnan)(v.stableNorm()));
172  VERIFY(!(numext::isfinite)(v.blueNorm())); VERIFY((numext::isnan)(v.blueNorm()));
173  VERIFY(!(numext::isfinite)(v.hypotNorm())); VERIFY((numext::isnan)(v.hypotNorm()));
174  }
175 
176  // stableNormalize[d]
177  {
178  VERIFY_IS_APPROX(vrand.stableNormalized(), vrand.normalized());
179  MatrixType vcopy(vrand);
180  vcopy.stableNormalize();
181  VERIFY_IS_APPROX(vcopy, vrand.normalized());
182  VERIFY_IS_APPROX((vrand.stableNormalized()).norm(), RealScalar(1));
183  VERIFY_IS_APPROX(vcopy.norm(), RealScalar(1));
184  VERIFY_IS_APPROX((vbig.stableNormalized()).norm(), RealScalar(1));
185  VERIFY_IS_APPROX((vsmall.stableNormalized()).norm(), RealScalar(1));
186  RealScalar big_scaling = ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
187  VERIFY_IS_APPROX(vbig/big_scaling, (vbig.stableNorm() * vbig.stableNormalized()).eval()/big_scaling);
188  VERIFY_IS_APPROX(vsmall, vsmall.stableNorm() * vsmall.stableNormalized());
189  }
190 }
191 
192 void test_stable_norm()
193 {
194  for(int i = 0; i < g_repeat; i++) {
195  CALL_SUBTEST_1( stable_norm(Matrix<float, 1, 1>()) );
196  CALL_SUBTEST_2( stable_norm(Vector4d()) );
197  CALL_SUBTEST_3( stable_norm(VectorXd(internal::random<int>(10,2000))) );
198  CALL_SUBTEST_4( stable_norm(VectorXf(internal::random<int>(10,2000))) );
199  CALL_SUBTEST_5( stable_norm(VectorXcd(internal::random<int>(10,2000))) );
200  }
201 }
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:33
max
#define max(a, b)
Definition: datatypes.h:20
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
min
#define min(a, b)
Definition: datatypes.h:19
v
ArrayXcf v
Definition: Cwise_arg.cpp:1
hypotNorm
EIGEN_DONT_INLINE T::Scalar hypotNorm(T &v)
Definition: bench_norm.cpp:21
sqrt
EIGEN_DEVICE_FUNC const SqrtReturnType sqrt() const
Definition: ArrayCwiseUnaryOps.h:152
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
isfinite
#define isfinite(X)
Definition: main.h:74
EIGEN_DONT_INLINE
#define EIGEN_DONT_INLINE
Definition: Macros.h:517
copy
EIGEN_DONT_INLINE T copy(const T &x)
Definition: stable_norm.cpp:12
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
first
constexpr int first(int i)
Implementation details for constexpr functions.
Definition: wrap/pybind11/include/pybind11/detail/common.h:570
Eigen::g_repeat
static int g_repeat
Definition: main.h:144
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
stableNorm
EIGEN_DONT_INLINE T::Scalar stableNorm(T &v)
Definition: bench_norm.cpp:15
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
abs2
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Abs2ReturnType abs2() const
Definition: ArrayCwiseUnaryOps.h:71
VERIFY_IS_MUCH_SMALLER_THAN
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition: main.h:335
inf
static double inf
Definition: testMatrix.cpp:30
VERIFY
#define VERIFY(a)
Definition: main.h:325
Eigen::isPlusInf
bool isPlusInf(const T &x)
Definition: main.h:663
stable_norm
void stable_norm(const MatrixType &m)
Definition: stable_norm.cpp:17
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
VERIFY_IS_NOT_APPROX
#define VERIFY_IS_NOT_APPROX(a, b)
Definition: integer_types.cpp:17
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:38
test_stable_norm
void test_stable_norm()
Definition: stable_norm.cpp:192
blueNorm
EIGEN_DONT_INLINE T::Scalar blueNorm(T &v)
Definition: bench_norm.cpp:27
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
abs
#define abs(x)
Definition: datatypes.h:17
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
isnan
#define isnan(X)
Definition: main.h:72

gtsam
Author(s):
autogenerated on Sat May 8 2021 02:44:49