GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-15 05:05:58 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix bug #740: overflow issue in stableNorm

Browse Source 
(grafted from 32915806305081d837711305bcf57508714d0068
)
This commit is contained in:
Gael Guennebaud 2014-02-13 15:44:01 +01:00
parent 24e33a0d86
commit 52dc1d7ffd
2 changed files with 31 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -17,15 +17,28 @@ namespace internal {
template<typename ExpressionType, typename Scalar> template<typename ExpressionType, typename Scalar>
inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale) inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale)
{ {
Scalar max = bl.cwiseAbs().maxCoeff(); using std::max;
if (max>scale) Scalar maxCoeff = bl.cwiseAbs().maxCoeff();
if (maxCoeff>scale)
{ {
ssq = ssq * numext::abs2(scale/max); ssq = ssq * numext::abs2(scale/maxCoeff);
scale = max; Scalar tmp = Scalar(1)/maxCoeff;
invScale = Scalar(1)/scale; if(tmp > NumTraits<Scalar>::highest())
{
invScale = NumTraits<Scalar>::highest();
scale = Scalar(1)/invScale;
} }
// TODO if the max is much much smaller than the current scale, else
{
scale = maxCoeff;
invScale = tmp;
}
}
// 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
if(scale>Scalar(0)) // if scale==0, then bl is 0
ssq += (bl*invScale).squaredNorm(); ssq += (bl*invScale).squaredNorm();
} }

12
test/stable_norm.cpp
View File

@ -55,8 +55,16 @@ template<typename MatrixType> void stable_norm(const MatrixType& m)
Index rows = m.rows(); Index rows = m.rows();
Index cols = m.cols(); Index cols = m.cols();
Scalar big = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4)); // get a non-zero random factor
Scalar small = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4)); Scalar factor = internal::random<Scalar>();
while(factor<RealScalar(1e-3))
factor = internal::random<Scalar>();
Scalar big = factor * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
factor = internal::random<Scalar>();
while(factor<RealScalar(1e-3))
factor = internal::random<Scalar>();
Scalar small = factor * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4));
MatrixType vzero = MatrixType::Zero(rows, cols), MatrixType vzero = MatrixType::Zero(rows, cols),
vrand = MatrixType::Random(rows, cols), vrand = MatrixType::Random(rows, cols),