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
13template<typename T> EIGEN_DONT_INLINE  T sub(T a, T b) { return a - b; }
14
15template<typename T> bool isFinite(const T& x)
16{
17  return isNotNaN(sub(x,x));
18}
19
20template<typename T> EIGEN_DONT_INLINE T copy(const T& x)
21{
22  return x;
23}
24
25template<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
106void 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