diff options
author | Gael Guennebaud <g.gael@free.fr> | 2014-02-13 15:44:01 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2014-02-13 15:44:01 +0100 |
commit | 32915806305081d837711305bcf57508714d0068 (patch) | |
tree | fc4843eb8ed812d2858148e388c6ad5e82f7ee88 /Eigen/src | |
parent | 14422decc2583e556352408983f4f8fda3054c70 (diff) |
Fix bug #740: overflow issue in stableNorm
Diffstat (limited to 'Eigen/src')
-rw-r--r-- | Eigen/src/Core/StableNorm.h | 27 |
1 files changed, 20 insertions, 7 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> |