add MatrixBase::stableNorm() avoiding over/under-flow

using it in QR reduced the error of Keir test from 1e-12 to 1e-24 but
that's much more expensive !

6 files changed, 46 insertions, 6 deletions

diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h
index fc8b254d7..b877cbee6 100644
--- a/Eigen/src/Core/DiagonalMatrix.h
+++ b/Eigen/src/Core/DiagonalMatrix.h
@@ -1,6 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra. Eigen itself is part of the KDE project.
 //
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h
index e165a841b..98262fca4 100644
--- a/Eigen/src/Core/Dot.h
+++ b/Eigen/src/Core/Dot.h
@@ -292,6 +292,18 @@ inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase<
   return ei_sqrt(squaredNorm());
 }
 
+/** \returns the \em l2 norm of \c *this using a numerically more stable
+  * algorithm.
+  *
+  * \sa norm(), dot(), squaredNorm()
+  */
+template<typename Derived>
+inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
+MatrixBase<Derived>::stableNorm() const
+{
+  return this->cwise().abs().redux(ei_scalar_hypot_op<RealScalar>());
+}
+
 /** \returns an expression of the quotient of *this by its own norm.
   *
   * \only_for_vectors
diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h
index c8ca3dac1..8b9516ae6 100644
--- a/Eigen/src/Core/Functors.h
+++ b/Eigen/src/Core/Functors.h
@@ -103,6 +103,28 @@ struct ei_functor_traits<ei_scalar_max_op<Scalar> > {
   };
 };
 
+/** \internal
+  * \brief Template functor to compute the hypot of two scalars
+  *
+  * \sa MatrixBase::stableNorm(), class Redux
+  */
+template<typename Scalar> struct ei_scalar_hypot_op EIGEN_EMPTY_STRUCT {
+//   typedef typename NumTraits<Scalar>::Real result_type;
+  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& _x, const Scalar& _y) const
+  {
+//     typedef typename NumTraits<T>::Real RealScalar;
+//     RealScalar _x = ei_abs(x);
+//     RealScalar _y = ei_abs(y);
+    Scalar p = std::max(_x, _y);
+    Scalar q = std::min(_x, _y);
+    Scalar qp = q/p;
+    return p * ei_sqrt(Scalar(1) + qp*qp);
+  }
+};
+template<typename Scalar>
+struct ei_functor_traits<ei_scalar_hypot_op<Scalar> > {
+  enum { Cost = 5 * NumTraits<Scalar>::MulCost };
+};
 
 // other binary functors:
 
diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index 8f56ef40e..975f64e18 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -34,10 +34,11 @@ template<typename T> inline T ei_random_amplitude()
   else return static_cast<T>(10);
 }
 
-template<typename T> inline T ei_hypot(T x, T y)
+template<typename T> inline typename NumTraits<T>::Real ei_hypot(T x, T y)
 {
-  T _x = ei_abs(x);
-  T _y = ei_abs(y);
+  typedef typename NumTraits<T>::Real RealScalar;
+  RealScalar _x = ei_abs(x);
+  RealScalar _y = ei_abs(y);
   T p = std::max(_x, _y);
   T q = std::min(_x, _y);
   T qp = q/p;
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 4dd319985..6b7b002f0 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -351,6 +351,7 @@ template<typename Derived> class MatrixBase
     Scalar dot(const MatrixBase<OtherDerived>& other) const;
     RealScalar squaredNorm() const;
     RealScalar norm()  const;
+    RealScalar stableNorm()  const;
     const PlainMatrixType normalized() const;
     void normalize();
 
diff --git a/test/adjoint.cpp b/test/adjoint.cpp
index f553bad02..f0c321057 100644
--- a/test/adjoint.cpp
+++ b/test/adjoint.cpp
@@ -65,15 +65,18 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
 
   // check basic properties of dot, norm, norm2
   typedef typename NumTraits<Scalar>::Real RealScalar;
-  VERIFY(ei_isApprox((s1 * v1 + s2 * v2).dot(v3),   s1 * v1.dot(v3) + s2 * v2.dot(v3), largerEps));
-  VERIFY(ei_isApprox(v3.dot(s1 * v1 + s2 * v2),     ei_conj(s1)*v3.dot(v1)+ei_conj(s2)*v3.dot(v2), largerEps));
+  VERIFY(ei_isApprox((s1 * v1 + s2 * v2).dot(v3),     s1 * v1.dot(v3) + s2 * v2.dot(v3), largerEps));
+  VERIFY(ei_isApprox(v3.dot(s1 * v1 + s2 * v2),       ei_conj(s1)*v3.dot(v1)+ei_conj(s2)*v3.dot(v2), largerEps));
   VERIFY_IS_APPROX(ei_conj(v1.dot(v2)),               v2.dot(v1));
   VERIFY_IS_APPROX(ei_abs(v1.dot(v1)),                v1.squaredNorm());
   if(NumTraits<Scalar>::HasFloatingPoint)
-    VERIFY_IS_APPROX(v1.squaredNorm(),                      v1.norm() * v1.norm());
+    VERIFY_IS_APPROX(v1.squaredNorm(),                v1.norm() * v1.norm());
   VERIFY_IS_MUCH_SMALLER_THAN(ei_abs(vzero.dot(v1)),  static_cast<RealScalar>(1));
   if(NumTraits<Scalar>::HasFloatingPoint)
+  {
     VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(),         static_cast<RealScalar>(1));
+    VERIFY_IS_APPROX(v1.norm(),                       v1.stableNorm());
+  }
 
   // check compatibility of dot and adjoint
   VERIFY(ei_isApprox(v1.dot(square * v2), (square.adjoint() * v1).dot(v2), largerEps));