diff options
Diffstat (limited to 'test/stable_norm.cpp')
| -rw-r--r-- | test/stable_norm.cpp | 105 |
1 files changed, 91 insertions, 14 deletions
diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp index 231dd9189..c3eb5ff31 100644 --- a/test/stable_norm.cpp +++ b/test/stable_norm.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -9,14 +9,6 @@ #include "main.h" -// workaround aggressive optimization in ICC -template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } - -template<typename T> bool isFinite(const T& x) -{ - return isNotNaN(sub(x,x)); -} - template<typename T> EIGEN_DONT_INLINE T copy(const T& x) { return x; @@ -32,6 +24,8 @@ template<typename MatrixType> void stable_norm(const MatrixType& m) typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; + + bool complex_real_product_ok = true; // Check the basic machine-dependent constants. { @@ -44,6 +38,16 @@ template<typename MatrixType> void stable_norm(const MatrixType& m) VERIFY( (!(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5) || (it<=4 && ibeta <= 3 ) || it<2)) && "the stable norm algorithm cannot be guaranteed on this computer"); + + Scalar inf = std::numeric_limits<RealScalar>::infinity(); + if(NumTraits<Scalar>::IsComplex && (numext::isnan)(inf*RealScalar(1)) ) + { + complex_real_product_ok = false; + static bool first = true; + if(first) + std::cerr << "WARNING: compiler mess up complex*real product, " << inf << " * " << 1.0 << " = " << inf*RealScalar(1) << std::endl; + first = false; + } } @@ -76,19 +80,19 @@ template<typename MatrixType> void stable_norm(const MatrixType& m) RealScalar size = static_cast<RealScalar>(m.size()); - // test isFinite - VERIFY(!isFinite( std::numeric_limits<RealScalar>::infinity())); - VERIFY(!isFinite(sqrt(-abs(big)))); + // test numext::isfinite + VERIFY(!(numext::isfinite)( std::numeric_limits<RealScalar>::infinity())); + VERIFY(!(numext::isfinite)(sqrt(-abs(big)))); // test overflow - VERIFY(isFinite(sqrt(size)*abs(big))); + VERIFY((numext::isfinite)(sqrt(size)*abs(big))); VERIFY_IS_NOT_APPROX(sqrt(copy(vbig.squaredNorm())), abs(sqrt(size)*big)); // here the default norm must fail VERIFY_IS_APPROX(vbig.stableNorm(), sqrt(size)*abs(big)); VERIFY_IS_APPROX(vbig.blueNorm(), sqrt(size)*abs(big)); VERIFY_IS_APPROX(vbig.hypotNorm(), sqrt(size)*abs(big)); // test underflow - VERIFY(isFinite(sqrt(size)*abs(small))); + VERIFY((numext::isfinite)(sqrt(size)*abs(small))); VERIFY_IS_NOT_APPROX(sqrt(copy(vsmall.squaredNorm())), abs(sqrt(size)*small)); // here the default norm must fail VERIFY_IS_APPROX(vsmall.stableNorm(), sqrt(size)*abs(small)); VERIFY_IS_APPROX(vsmall.blueNorm(), sqrt(size)*abs(small)); @@ -101,6 +105,79 @@ template<typename MatrixType> void stable_norm(const MatrixType& m) VERIFY_IS_APPROX(vrand.rowwise().stableNorm(), vrand.rowwise().norm()); VERIFY_IS_APPROX(vrand.rowwise().blueNorm(), vrand.rowwise().norm()); VERIFY_IS_APPROX(vrand.rowwise().hypotNorm(), vrand.rowwise().norm()); + + // test NaN, +inf, -inf + MatrixType v; + Index i = internal::random<Index>(0,rows-1); + Index j = internal::random<Index>(0,cols-1); + + // NaN + { + v = vrand; + v(i,j) = std::numeric_limits<RealScalar>::quiet_NaN(); + VERIFY(!(numext::isfinite)(v.squaredNorm())); VERIFY((numext::isnan)(v.squaredNorm())); + VERIFY(!(numext::isfinite)(v.norm())); VERIFY((numext::isnan)(v.norm())); + VERIFY(!(numext::isfinite)(v.stableNorm())); VERIFY((numext::isnan)(v.stableNorm())); + VERIFY(!(numext::isfinite)(v.blueNorm())); VERIFY((numext::isnan)(v.blueNorm())); + VERIFY(!(numext::isfinite)(v.hypotNorm())); VERIFY((numext::isnan)(v.hypotNorm())); + } + + // +inf + { + v = vrand; + v(i,j) = std::numeric_limits<RealScalar>::infinity(); + VERIFY(!(numext::isfinite)(v.squaredNorm())); VERIFY(isPlusInf(v.squaredNorm())); + VERIFY(!(numext::isfinite)(v.norm())); VERIFY(isPlusInf(v.norm())); + VERIFY(!(numext::isfinite)(v.stableNorm())); + if(complex_real_product_ok){ + VERIFY(isPlusInf(v.stableNorm())); + } + VERIFY(!(numext::isfinite)(v.blueNorm())); VERIFY(isPlusInf(v.blueNorm())); + VERIFY(!(numext::isfinite)(v.hypotNorm())); VERIFY(isPlusInf(v.hypotNorm())); + } + + // -inf + { + v = vrand; + v(i,j) = -std::numeric_limits<RealScalar>::infinity(); + VERIFY(!(numext::isfinite)(v.squaredNorm())); VERIFY(isPlusInf(v.squaredNorm())); + VERIFY(!(numext::isfinite)(v.norm())); VERIFY(isPlusInf(v.norm())); + VERIFY(!(numext::isfinite)(v.stableNorm())); + if(complex_real_product_ok) { + VERIFY(isPlusInf(v.stableNorm())); + } + VERIFY(!(numext::isfinite)(v.blueNorm())); VERIFY(isPlusInf(v.blueNorm())); + VERIFY(!(numext::isfinite)(v.hypotNorm())); VERIFY(isPlusInf(v.hypotNorm())); + } + + // mix + { + Index i2 = internal::random<Index>(0,rows-1); + Index j2 = internal::random<Index>(0,cols-1); + v = vrand; + v(i,j) = -std::numeric_limits<RealScalar>::infinity(); + v(i2,j2) = std::numeric_limits<RealScalar>::quiet_NaN(); + VERIFY(!(numext::isfinite)(v.squaredNorm())); VERIFY((numext::isnan)(v.squaredNorm())); + VERIFY(!(numext::isfinite)(v.norm())); VERIFY((numext::isnan)(v.norm())); + VERIFY(!(numext::isfinite)(v.stableNorm())); VERIFY((numext::isnan)(v.stableNorm())); + VERIFY(!(numext::isfinite)(v.blueNorm())); VERIFY((numext::isnan)(v.blueNorm())); + VERIFY(!(numext::isfinite)(v.hypotNorm())); VERIFY((numext::isnan)(v.hypotNorm())); + } + + // stableNormalize[d] + { + VERIFY_IS_APPROX(vrand.stableNormalized(), vrand.normalized()); + MatrixType vcopy(vrand); + vcopy.stableNormalize(); + VERIFY_IS_APPROX(vcopy, vrand.normalized()); + VERIFY_IS_APPROX((vrand.stableNormalized()).norm(), RealScalar(1)); + VERIFY_IS_APPROX(vcopy.norm(), RealScalar(1)); + VERIFY_IS_APPROX((vbig.stableNormalized()).norm(), RealScalar(1)); + VERIFY_IS_APPROX((vsmall.stableNormalized()).norm(), RealScalar(1)); + RealScalar big_scaling = ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4)); + VERIFY_IS_APPROX(vbig/big_scaling, (vbig.stableNorm() * vbig.stableNormalized()).eval()/big_scaling); + VERIFY_IS_APPROX(vsmall, vsmall.stableNorm() * vsmall.stableNormalized()); + } } void test_stable_norm() |