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 #include "main.h"
11
12 // workaround aggressive optimization in ICC
sub(T a,T b)13 template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; }
14
isFinite(const T & x)15 template<typename T> bool isFinite(const T& x)
16 {
17 return isNotNaN(sub(x,x));
18 }
19
copy(const T & x)20 template<typename T> EIGEN_DONT_INLINE T copy(const T& x)
21 {
22 return x;
23 }
24
stable_norm(const MatrixType & m)25 template<typename MatrixType> void stable_norm(const MatrixType& m)
26 {
27 /* this test covers the following files:
28 StableNorm.h
29 */
30 using std::sqrt;
31 using std::abs;
32 typedef typename MatrixType::Index Index;
33 typedef typename MatrixType::Scalar Scalar;
34 typedef typename NumTraits<Scalar>::Real RealScalar;
35
36 // Check the basic machine-dependent constants.
37 {
38 int ibeta, it, iemin, iemax;
39
40 ibeta = std::numeric_limits<RealScalar>::radix; // base for floating-point numbers
41 it = std::numeric_limits<RealScalar>::digits; // number of base-beta digits in mantissa
42 iemin = std::numeric_limits<RealScalar>::min_exponent; // minimum exponent
43 iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent
44
45 VERIFY( (!(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5) || (it<=4 && ibeta <= 3 ) || it<2))
46 && "the stable norm algorithm cannot be guaranteed on this computer");
47 }
48
49
50 Index rows = m.rows();
51 Index cols = m.cols();
52
53 // get a non-zero random factor
54 Scalar factor = internal::random<Scalar>();
55 while(numext::abs2(factor)<RealScalar(1e-4))
56 factor = internal::random<Scalar>();
57 Scalar big = factor * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
58
59 factor = internal::random<Scalar>();
60 while(numext::abs2(factor)<RealScalar(1e-4))
61 factor = internal::random<Scalar>();
62 Scalar small = factor * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4));
63
64 MatrixType vzero = MatrixType::Zero(rows, cols),
65 vrand = MatrixType::Random(rows, cols),
66 vbig(rows, cols),
67 vsmall(rows,cols);
68
69 vbig.fill(big);
70 vsmall.fill(small);
71
72 VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1));
73 VERIFY_IS_APPROX(vrand.stableNorm(), vrand.norm());
74 VERIFY_IS_APPROX(vrand.blueNorm(), vrand.norm());
75 VERIFY_IS_APPROX(vrand.hypotNorm(), vrand.norm());
76
77 RealScalar size = static_cast<RealScalar>(m.size());
78
79 // test isFinite
80 VERIFY(!isFinite( std::numeric_limits<RealScalar>::infinity()));
81 VERIFY(!isFinite(sqrt(-abs(big))));
82
83 // test overflow
84 VERIFY(isFinite(sqrt(size)*abs(big)));
85 VERIFY_IS_NOT_APPROX(sqrt(copy(vbig.squaredNorm())), abs(sqrt(size)*big)); // here the default norm must fail
86 VERIFY_IS_APPROX(vbig.stableNorm(), sqrt(size)*abs(big));
87 VERIFY_IS_APPROX(vbig.blueNorm(), sqrt(size)*abs(big));
88 VERIFY_IS_APPROX(vbig.hypotNorm(), sqrt(size)*abs(big));
89
90 // test underflow
91 VERIFY(isFinite(sqrt(size)*abs(small)));
92 VERIFY_IS_NOT_APPROX(sqrt(copy(vsmall.squaredNorm())), abs(sqrt(size)*small)); // here the default norm must fail
93 VERIFY_IS_APPROX(vsmall.stableNorm(), sqrt(size)*abs(small));
94 VERIFY_IS_APPROX(vsmall.blueNorm(), sqrt(size)*abs(small));
95 VERIFY_IS_APPROX(vsmall.hypotNorm(), sqrt(size)*abs(small));
96
97 // Test compilation of cwise() version
98 VERIFY_IS_APPROX(vrand.colwise().stableNorm(), vrand.colwise().norm());
99 VERIFY_IS_APPROX(vrand.colwise().blueNorm(), vrand.colwise().norm());
100 VERIFY_IS_APPROX(vrand.colwise().hypotNorm(), vrand.colwise().norm());
101 VERIFY_IS_APPROX(vrand.rowwise().stableNorm(), vrand.rowwise().norm());
102 VERIFY_IS_APPROX(vrand.rowwise().blueNorm(), vrand.rowwise().norm());
103 VERIFY_IS_APPROX(vrand.rowwise().hypotNorm(), vrand.rowwise().norm());
104 }
105
test_stable_norm()106 void test_stable_norm()
107 {
108 for(int i = 0; i < g_repeat; i++) {
109 CALL_SUBTEST_1( stable_norm(Matrix<float, 1, 1>()) );
110 CALL_SUBTEST_2( stable_norm(Vector4d()) );
111 CALL_SUBTEST_3( stable_norm(VectorXd(internal::random<int>(10,2000))) );
112 CALL_SUBTEST_4( stable_norm(VectorXf(internal::random<int>(10,2000))) );
113 CALL_SUBTEST_5( stable_norm(VectorXcd(internal::random<int>(10,2000))) );
114 }
115 }
116