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  if (i2 != i || j2 != j) {
174  // hypot propagates inf over NaN.
175  VERIFY(!(numext::isfinite)(v.hypotNorm())); VERIFY((numext::isinf)(v.hypotNorm()));
176  } else {
177  // inf is overwritten by NaN, expect norm to be NaN.
178  VERIFY(!(numext::isfinite)(v.hypotNorm())); VERIFY((numext::isnan)(v.hypotNorm()));
179  }
180  }
181 
182  // stableNormalize[d]
183  {
184  VERIFY_IS_APPROX(vrand.stableNormalized(), vrand.normalized());
185  MatrixType vcopy(vrand);
186  vcopy.stableNormalize();
187  VERIFY_IS_APPROX(vcopy, vrand.normalized());
188  VERIFY_IS_APPROX((vrand.stableNormalized()).norm(), RealScalar(1));
189  VERIFY_IS_APPROX(vcopy.norm(), RealScalar(1));
190  VERIFY_IS_APPROX((vbig.stableNormalized()).norm(), RealScalar(1));
191  VERIFY_IS_APPROX((vsmall.stableNormalized()).norm(), RealScalar(1));
192  RealScalar big_scaling = ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
193  VERIFY_IS_APPROX(vbig/big_scaling, (vbig.stableNorm() * vbig.stableNormalized()).eval()/big_scaling);
194  VERIFY_IS_APPROX(vsmall, vsmall.stableNorm() * vsmall.stableNormalized());
195  }
196 }
197 
198 template<typename Scalar>
199 void test_hypot()
200 {
201  typedef typename NumTraits<Scalar>::Real RealScalar;
202  Scalar factor = internal::random<Scalar>();
203  while(numext::abs2(factor)<RealScalar(1e-4))
204  factor = internal::random<Scalar>();
205  Scalar big = factor * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
206 
207  factor = internal::random<Scalar>();
208  while(numext::abs2(factor)<RealScalar(1e-4))
209  factor = internal::random<Scalar>();
210  Scalar small = factor * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4));
211 
212  Scalar one (1),
213  zero (0),
214  sqrt2 (std::sqrt(2)),
215  nan (std::numeric_limits<RealScalar>::quiet_NaN());
216 
217  Scalar a = internal::random<Scalar>(-1,1);
218  Scalar b = internal::random<Scalar>(-1,1);
219  VERIFY_IS_APPROX(numext::hypot(a,b),std::sqrt(numext::abs2(a)+numext::abs2(b)));
220  VERIFY_IS_EQUAL(numext::hypot(zero,zero), zero);
221  VERIFY_IS_APPROX(numext::hypot(one, one), sqrt2);
222  VERIFY_IS_APPROX(numext::hypot(big,big), sqrt2*numext::abs(big));
223  VERIFY_IS_APPROX(numext::hypot(small,small), sqrt2*numext::abs(small));
224  VERIFY_IS_APPROX(numext::hypot(small,big), numext::abs(big));
225  VERIFY((numext::isnan)(numext::hypot(nan,a)));
226  VERIFY((numext::isnan)(numext::hypot(a,nan)));
227 }
228 
229 EIGEN_DECLARE_TEST(stable_norm)
230 {
231  for(int i = 0; i < g_repeat; i++) {
232  CALL_SUBTEST_3( test_hypot<double>() );
233  CALL_SUBTEST_4( test_hypot<float>() );
234  CALL_SUBTEST_5( test_hypot<std::complex<double> >() );
235  CALL_SUBTEST_6( test_hypot<std::complex<float> >() );
236 
237  CALL_SUBTEST_1( stable_norm(Matrix<float, 1, 1>()) );
238  CALL_SUBTEST_2( stable_norm(Vector4d()) );
239  CALL_SUBTEST_3( stable_norm(VectorXd(internal::random<int>(10,2000))) );
240  CALL_SUBTEST_3( stable_norm(MatrixXd(internal::random<int>(10,200), internal::random<int>(10,200))) );
241  CALL_SUBTEST_4( stable_norm(VectorXf(internal::random<int>(10,2000))) );
242  CALL_SUBTEST_5( stable_norm(VectorXcd(internal::random<int>(10,2000))) );
243  CALL_SUBTEST_6( stable_norm(VectorXcf(internal::random<int>(10,2000))) );
244  }
245 }
VERIFY_IS_MUCH_SMALLER_THAN
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition: main.h:390
stable_norm
void stable_norm(const MatrixType &m)
Definition: stable_norm.cpp:17
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
b
Scalar * b
Definition: benchVecAdd.cpp:17
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
Eigen::isPlusInf
bool isPlusInf(const T &x)
Definition: main.h:713
hypotNorm
EIGEN_DONT_INLINE T::Scalar hypotNorm(T &v)
Definition: bench_norm.cpp:21
abs2
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE Abs2ReturnType abs2() const
Definition: ArrayCwiseUnaryOps.h:80
isnan
#define isnan(X)
Definition: main.h:93
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
test_hypot
void test_hypot()
Definition: stable_norm.cpp:199
VERIFY_IS_NOT_APPROX
#define VERIFY_IS_NOT_APPROX(a, b)
Definition: integer_types.cpp:17
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
copy
EIGEN_DONT_INLINE T copy(const T &x)
Definition: stable_norm.cpp:12
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
zero
EIGEN_DONT_INLINE Scalar zero()
Definition: svd_common.h:296
Eigen::internal::first
EIGEN_CONSTEXPR Index first(const T &x) EIGEN_NOEXCEPT
Definition: IndexedViewHelper.h:81
isfinite
#define isfinite(X)
Definition: main.h:95
stableNorm
EIGEN_DONT_INLINE T::Scalar stableNorm(T &v)
Definition: bench_norm.cpp:15
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(stable_norm)
Definition: stable_norm.cpp:229
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
Eigen::Triplet< double >
CALL_SUBTEST_6
#define CALL_SUBTEST_6(FUNC)
Definition: split_test_helper.h:34
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
main.h
big
static double big
Definition: igam.c:119
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
min
#define min(a, b)
Definition: datatypes.h:19
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
abs
#define abs(x)
Definition: datatypes.h:17
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
inf
static double inf
Definition: testMatrix.cpp:31
max
#define max(a, b)
Definition: datatypes.h:20
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
VERIFY
#define VERIFY(a)
Definition: main.h:380
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_DONT_INLINE
#define EIGEN_DONT_INLINE
Definition: Macros.h:940
isinf
#define isinf(X)
Definition: main.h:94

gtsam
Author(s):
autogenerated on Tue Jun 25 2024 03:03:28