1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009-2014 Gael Guennebaud <[email protected]>
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
copy(const T & x)12 template<typename T> EIGEN_DONT_INLINE T copy(const T& x)
13 {
14 return x;
15 }
16
stable_norm(const MatrixType & m)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>
test_hypot()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
EIGEN_DECLARE_TEST(stable_norm)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 }
246