diff options
| -rw-r--r-- | Eigen/src/Core/StableNorm.h | 27 | ||||
| -rw-r--r-- | test/stable_norm.cpp | 14 |
2 files changed, 31 insertions, 10 deletions
diff --git a/Eigen/src/Core/StableNorm.h b/Eigen/src/Core/StableNorm.h index c219e2f53..525620bad 100644 --- a/Eigen/src/Core/StableNorm.h +++ b/Eigen/src/Core/StableNorm.h @@ -17,16 +17,29 @@ namespace internal { template<typename ExpressionType, typename Scalar> inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale) { - Scalar max = bl.cwiseAbs().maxCoeff(); - if (max>scale) + using std::max; + Scalar maxCoeff = bl.cwiseAbs().maxCoeff(); + + if (maxCoeff>scale) { - ssq = ssq * numext::abs2(scale/max); - scale = max; - invScale = Scalar(1)/scale; + ssq = ssq * numext::abs2(scale/maxCoeff); + Scalar tmp = Scalar(1)/maxCoeff; + if(tmp > NumTraits<Scalar>::highest()) + { + invScale = NumTraits<Scalar>::highest(); + scale = Scalar(1)/invScale; + } + else + { + scale = maxCoeff; + invScale = tmp; + } } - // TODO if the max 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 - ssq += (bl*invScale).squaredNorm(); + if(scale>Scalar(0)) // if scale==0, then bl is 0 + ssq += (bl*invScale).squaredNorm(); } template<typename Derived> diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp index c09fc17b7..364170acd 100644 --- a/test/stable_norm.cpp +++ b/test/stable_norm.cpp @@ -55,8 +55,16 @@ template<typename MatrixType> void stable_norm(const MatrixType& m) Index rows = m.rows(); Index cols = m.cols(); - Scalar big = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4)); - Scalar small = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4)); + // get a non-zero random factor + 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), vrand = MatrixType::Random(rows, cols), @@ -91,7 +99,7 @@ template<typename MatrixType> void stable_norm(const MatrixType& m) VERIFY_IS_APPROX(vsmall.blueNorm(), sqrt(size)*abs(small)); VERIFY_IS_APPROX(vsmall.hypotNorm(), sqrt(size)*abs(small)); -// Test compilation of cwise() version + // Test compilation of cwise() version VERIFY_IS_APPROX(vrand.colwise().stableNorm(), vrand.colwise().norm()); VERIFY_IS_APPROX(vrand.colwise().blueNorm(), vrand.colwise().norm()); VERIFY_IS_APPROX(vrand.colwise().hypotNorm(), vrand.colwise().norm()); |