stable_norm.cpp (10379B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009-2014 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 template<typename T> EIGEN_DONT_INLINE T copy(const T& x) 13 { 14 return x; 15 } 16 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> 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 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 }