GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-22 01:29:35 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix stableNorm() with respect to NaN and inf, and add respective unit tests. blueNorm() and hypotNorm() are broken wrt to NaN/inf

Browse Source
This commit is contained in:
Gael Guennebaud 2014-09-02 14:49:23 +02:00
parent 3eb5253ca1
commit 18fbe7e7d4
2 changed files with 78 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
Eigen/src/Core/StableNorm.h
View File

@ -20,7 +20,7 @@ inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& sc
using std::max; using std::max;
Scalar maxCoeff = bl.cwiseAbs().maxCoeff(); Scalar maxCoeff = bl.cwiseAbs().maxCoeff();
if (maxCoeff>scale) if(maxCoeff>scale)
{ {
ssq = ssq * numext::abs2(scale/maxCoeff); ssq = ssq * numext::abs2(scale/maxCoeff);
Scalar tmp = Scalar(1)/maxCoeff; Scalar tmp = Scalar(1)/maxCoeff;
@ -29,12 +29,21 @@ inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& sc
invScale = NumTraits<Scalar>::highest(); invScale = NumTraits<Scalar>::highest();
scale = Scalar(1)/invScale; scale = Scalar(1)/invScale;
} }
else if(maxCoeff>NumTraits<Scalar>::highest()) // we got a INF
{
invScale = Scalar(1);
scale = maxCoeff;
}
else else
{ {
scale = maxCoeff; scale = maxCoeff;
invScale = tmp; invScale = tmp;
} }
} }
else if(maxCoeff!=maxCoeff) // we got a NaN
{
scale = maxCoeff;
}
// TODO if the maxCoeff is much much smaller than the current scale, // TODO if the maxCoeff is much much smaller than the current scale,
// then we can neglect this sub vector // then we can neglect this sub vector

69
test/stable_norm.cpp
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library // This file is part of Eigen, a lightweight C++ template library
// for linear algebra. // 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 // 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 // Public License v. 2.0. If a copy of the MPL was not distributed
@ -14,6 +14,21 @@ template<typename T> bool isNotNaN(const T& x)
return x==x; return x==x;
} }
template<typename T> bool isNaN(const T& x)
{
return x!=x;
}
template<typename T> bool isInf(const T& x)
{
return x > NumTraits<T>::highest();
}
template<typename T> bool isMinusInf(const T& x)
{
return x < NumTraits<T>::lowest();
}
// workaround aggressive optimization in ICC // workaround aggressive optimization in ICC
template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; }
@ -106,6 +121,58 @@ template<typename MatrixType> void stable_norm(const MatrixType& m)
VERIFY_IS_APPROX(vrand.rowwise().stableNorm(), vrand.rowwise().norm()); VERIFY_IS_APPROX(vrand.rowwise().stableNorm(), vrand.rowwise().norm());
VERIFY_IS_APPROX(vrand.rowwise().blueNorm(), vrand.rowwise().norm()); VERIFY_IS_APPROX(vrand.rowwise().blueNorm(), vrand.rowwise().norm());
VERIFY_IS_APPROX(vrand.rowwise().hypotNorm(), 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) = RealScalar(0)/RealScalar(0);
VERIFY(!isFinite(v.squaredNorm())); VERIFY(isNaN(v.squaredNorm()));
VERIFY(!isFinite(v.norm())); VERIFY(isNaN(v.norm()));
VERIFY(!isFinite(v.stableNorm())); VERIFY(isNaN(v.stableNorm()));
VERIFY(!isFinite(v.blueNorm())); VERIFY(isNaN(v.blueNorm()));
// VERIFY(!isFinite(v.hypotNorm())); //VERIFY(isNaN(v.hypotNorm()));
}
// +inf
{
v = vrand;
v(i,j) = RealScalar(1)/RealScalar(0);
VERIFY(!isFinite(v.squaredNorm())); VERIFY(isInf(v.squaredNorm()));
VERIFY(!isFinite(v.norm())); VERIFY(isInf(v.norm()));
VERIFY(!isFinite(v.stableNorm())); VERIFY(isInf(v.stableNorm()));
// VERIFY(!isFinite(v.blueNorm())); //VERIFY(isInf(v.blueNorm()));
// VERIFY(!isFinite(v.hypotNorm())); //VERIFY(isInf(v.hypotNorm()));
}
// -inf
{
v = vrand;
v(i,j) = RealScalar(-1)/RealScalar(0);
VERIFY(!isFinite(v.squaredNorm())); VERIFY(isInf(v.squaredNorm()));
VERIFY(!isFinite(v.norm())); VERIFY(isInf(v.norm()));
VERIFY(!isFinite(v.stableNorm())); VERIFY(isInf(v.stableNorm()));
// VERIFY(!isFinite(v.blueNorm())); VERIFY(isInf(v.blueNorm()));
// VERIFY(!isFinite(v.hypotNorm())); VERIFY(isInf(v.hypotNorm()));
}
// mix
{
Index i2 = internal::random<Index>(0,rows-1);
Index j2 = internal::random<Index>(0,cols-1);
v = vrand;
v(i,j) = RealScalar(-1)/RealScalar(0);
v(i2,j2) = RealScalar(0)/RealScalar(0);
VERIFY(!isFinite(v.squaredNorm())); VERIFY(isNaN(v.squaredNorm()));
VERIFY(!isFinite(v.norm())); VERIFY(isNaN(v.norm()));
VERIFY(!isFinite(v.stableNorm())); VERIFY(isNaN(v.stableNorm()));
// VERIFY(!isFinite(v.blueNorm())); //VERIFY(isNaN(v.blueNorm()));
// VERIFY(!isFinite(v.hypotNorm())); VERIFY(isNaN(v.hypotNorm()));
}
} }
void test_stable_norm() void test_stable_norm()