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 (at) 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 } 116
