GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-25 07:43:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source 
using it in QR reduced the error of Keir test from 1e-12 to 1e-24 but
that's much more expensive !
This commit is contained in:
Gael Guennebaud 2009-01-28 22:11:56 +00:00
parent 42b237b83a
commit 36c8a64923
6 changed files with 46 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Eigen/src/Core/DiagonalMatrix.h
View File

@ -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

12
Eigen/src/Core/Dot.h
View File

@ -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

22
Eigen/src/Core/Functors.h
View File

@ -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:

7
Eigen/src/Core/MathFunctions.h
View File

@ -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;

1
Eigen/src/Core/MatrixBase.h
View File

@ -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();

9
test/adjoint.cpp
View File

@ -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));