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