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

bug #701: workaround (min) and (max) blocking ADL by introducing numext::mini and numext::maxi internal functions and a EIGEN_NOT_A_MACRO macro.

Browse Source
This commit is contained in:
Gael Guennebaud 2014-10-20 15:55:32 +02:00
parent c12b7896d0
commit fe57b2f963
17 changed files with 43 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Eigen/src/Cholesky/LDLT.h
View File

@ -488,11 +488,10 @@ void LDLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) cons
// dst = D^-1 (L^-1 P b)
// more precisely, use pseudo-inverse of D (see bug 241)
using std::abs;
EIGEN_USING_STD_MATH(max);
const typename Diagonal<const MatrixType>::RealReturnType vecD(vectorD());
// In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon
// as motivated by LAPACK's xGELSS:
// RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() *NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest());
// RealScalar tolerance = numext::maxi(vectorD.array().abs().maxCoeff() *NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest());
// However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the highest
// diagonal element is not well justified and to numerical issues in some cases.
// Moreover, Lapack's xSYTRS routines use 0 for the tolerance.

5
Eigen/src/Core/Diagonal.h
View File

@ -77,9 +77,8 @@ template<typename MatrixType, int _DiagIndex> class Diagonal
EIGEN_DEVICE_FUNC
inline Index rows() const
{
EIGEN_USING_STD_MATH(min);
return m_index.value()<0 ? (min)(Index(m_matrix.cols()),Index(m_matrix.rows()+m_index.value()))
: (min)(Index(m_matrix.rows()),Index(m_matrix.cols()-m_index.value()));
return m_index.value()<0 ? numext::mini(Index(m_matrix.cols()),Index(m_matrix.rows()+m_index.value()))
: numext::mini(Index(m_matrix.rows()),Index(m_matrix.cols()-m_index.value()));
}
EIGEN_DEVICE_FUNC

3
Eigen/src/Core/Fuzzy.h
View File

@ -22,10 +22,9 @@ struct isApprox_selector
EIGEN_DEVICE_FUNC
static bool run(const Derived& x, const OtherDerived& y, const typename Derived::RealScalar& prec)
{
EIGEN_USING_STD_MATH(min);
typename internal::nested_eval<Derived,2>::type nested(x);
typename internal::nested_eval<OtherDerived,2>::type otherNested(y);
return (nested - otherNested).cwiseAbs2().sum() <= prec * prec * (min)(nested.cwiseAbs2().sum(), otherNested.cwiseAbs2().sum());
return (nested - otherNested).cwiseAbs2().sum() <= prec * prec * numext::mini(nested.cwiseAbs2().sum(), otherNested.cwiseAbs2().sum());
}
};

4
Eigen/src/Core/GenericPacketMath.h
View File

@ -126,12 +126,12 @@ pdiv(const Packet& a,
/** \internal \returns the min of \a a and \a b (coeff-wise) */
template<typename Packet> EIGEN_DEVICE_FUNC inline Packet
pmin(const Packet& a,
const Packet& b) { EIGEN_USING_STD_MATH(min); return (min)(a, b); }
const Packet& b) { return numext::mini(a, b); }
/** \internal \returns the max of \a a and \a b (coeff-wise) */
template<typename Packet> EIGEN_DEVICE_FUNC inline Packet
pmax(const Packet& a,
const Packet& b) { EIGEN_USING_STD_MATH(max); return (max)(a, b); }
const Packet& b) { return numext::maxi(a, b); }
/** \internal \returns the absolute value of \a a */
template<typename Packet> EIGEN_DEVICE_FUNC inline Packet

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

@ -591,6 +591,22 @@ inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random()
****************************************************************************/
namespace numext {
template<typename T>
EIGEN_DEVICE_FUNC
inline T mini(const T& x, const T& y)
{
using std::min;
return min EIGEN_NOT_A_MACRO (x,y);
}
template<typename T>
EIGEN_DEVICE_FUNC
inline T maxi(const T& x, const T& y)
{
using std::max;
return max EIGEN_NOT_A_MACRO (x,y);
}
template<typename Scalar>
EIGEN_DEVICE_FUNC

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

@ -17,7 +17,6 @@ namespace internal {
template<typename ExpressionType, typename Scalar>
inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale)
{
using std::max;
Scalar maxCoeff = bl.cwiseAbs().maxCoeff();
if(maxCoeff>scale)
@ -58,8 +57,6 @@ blueNorm_impl(const EigenBase<Derived>& _vec)
typedef typename Derived::RealScalar RealScalar;
typedef typename Derived::Index Index;
using std::pow;
EIGEN_USING_STD_MATH(min);
EIGEN_USING_STD_MATH(max);
using std::sqrt;
using std::abs;
const Derived& vec(_vec.derived());
@ -136,8 +133,8 @@ blueNorm_impl(const EigenBase<Derived>& _vec)
}
else
return sqrt(amed);
asml = (min)(abig, amed);
abig = (max)(abig, amed);
asml = numext::mini(abig, amed);
abig = numext::maxi(abig, amed);
if(asml <= abig*relerr)
return abig;
else
@ -160,7 +157,6 @@ template<typename Derived>
inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
MatrixBase<Derived>::stableNorm() const
{
EIGEN_USING_STD_MATH(min);
using std::sqrt;
const Index blockSize = 4096;
RealScalar scale(0);
@ -174,7 +170,7 @@ MatrixBase<Derived>::stableNorm() const
if (bi>0)
internal::stable_norm_kernel(this->head(bi), ssq, scale, invScale);
for (; bi<n; bi+=blockSize)
internal::stable_norm_kernel(this->segment(bi,(min)(blockSize, n - bi)).template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
internal::stable_norm_kernel(this->segment(bi,numext::mini(blockSize, n - bi)).template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
return scale * sqrt(ssq);
}

6
Eigen/src/Core/functors/BinaryFunctors.h
View File

@ -115,7 +115,7 @@ struct functor_traits<scalar_conj_product_op<LhsScalar,RhsScalar> > {
*/
template<typename Scalar> struct scalar_min_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_min_op)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { EIGEN_USING_STD_MATH(min); return (min)(a, b); }
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return numext::mini(a, b); }
template<typename Packet>
EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
{ return internal::pmin(a,b); }
@ -138,7 +138,7 @@ struct functor_traits<scalar_min_op<Scalar> > {
*/
template<typename Scalar> struct scalar_max_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_max_op)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { EIGEN_USING_STD_MATH(max); return (max)(a, b); }
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return numext::maxi(a, b); }
template<typename Packet>
EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
{ return internal::pmax(a,b); }
@ -164,8 +164,6 @@ template<typename Scalar> struct scalar_hypot_op {
// typedef typename NumTraits<Scalar>::Real result_type;
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& _x, const Scalar& _y) const
{
EIGEN_USING_STD_MATH(max);
EIGEN_USING_STD_MATH(min);
using std::sqrt;
Scalar p, qp;
if(_x>_y)

5
Eigen/src/Core/util/Macros.h
View File

@ -86,6 +86,11 @@
#define EIGEN_ALIGN 0
#endif
// This macro can be used to prevent from macro expansion, e.g.:
// std::max EIGNE_NOT_A_MACRO(a,b)
#define EIGEN_NOT_A_MACRO
// EIGEN_ALIGN_STATICALLY is the true test whether we want to align arrays on the stack or not. It takes into account both the user choice to explicitly disable
// alignment (EIGEN_DONT_ALIGN_STATICALLY) and the architecture config (EIGEN_ARCH_WANTS_STACK_ALIGNMENT). Henceforth, only EIGEN_ALIGN_STATICALLY should be used.
#if EIGEN_ARCH_WANTS_STACK_ALIGNMENT && !defined(EIGEN_DONT_ALIGN_STATICALLY)

6
Eigen/src/Eigenvalues/EigenSolver.h
View File

@ -368,7 +368,6 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect
{
using std::sqrt;
using std::abs;
using std::max;
using numext::isfinite;
eigen_assert(matrix.cols() == matrix.rows());
@ -409,7 +408,7 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect
{
Scalar t0 = m_matT.coeff(i+1, i);
Scalar t1 = m_matT.coeff(i, i+1);
Scalar maxval = (max)(abs(p),(max)(abs(t0),abs(t1)));
Scalar maxval = numext::maxi(abs(p),numext::maxi(abs(t0),abs(t1)));
t0 /= maxval;
t1 /= maxval;
Scalar p0 = p/maxval;
@ -600,8 +599,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
}
// Overflow control
EIGEN_USING_STD_MATH(max);
Scalar t = (max)(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
Scalar t = numext::maxi(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
if ((eps * t) * t > Scalar(1))
m_matT.block(i, n-1, size-i, 2) /= t;

3
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -732,7 +732,6 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false>
EIGEN_DEVICE_FUNC
static inline void run(SolverType& solver, const MatrixType& mat, int options)
{
EIGEN_USING_STD_MATH(max)
EIGEN_USING_STD_MATH(sqrt);
eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
@ -746,7 +745,7 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false>
// map the matrix coefficients to [-1:1] to avoid over- and underflow.
Scalar scale = mat.cwiseAbs().maxCoeff();
scale = (max)(scale,Scalar(1));
scale = numext::maxi(scale,Scalar(1));
MatrixType scaledMat = mat / scale;
// Compute the eigenvalues

3
Eigen/src/Geometry/Quaternion.h
View File

@ -571,7 +571,6 @@ template<class Derived>
template<typename Derived1, typename Derived2>
inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b)
{
EIGEN_USING_STD_MATH(max);
using std::sqrt;
Vector3 v0 = a.normalized();
Vector3 v1 = b.normalized();
@ -587,7 +586,7 @@ inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Deri
// which yields a singular value problem
if (c < Scalar(-1)+NumTraits<Scalar>::dummy_precision())
{
c = (max)(c,Scalar(-1));
c = numext::maxi(c,Scalar(-1));
Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV);
Vector3 axis = svd.matrixV().col(2);

3
Eigen/src/SVD/JacobiSVD.h
View File

@ -723,8 +723,7 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
// if this 2x2 sub-matrix is not diagonal already...
// notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
// keep us iterating forever. Similarly, small denormal numbers are considered zero.
EIGEN_USING_STD_MATH(max);
RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)),
RealScalar threshold = numext::maxi(considerAsZero, precision * numext::maxi(abs(m_workMatrix.coeff(p,p)),
abs(m_workMatrix.coeff(q,q))));
// We compare both values to threshold instead of calling max to be robust to NaN (See bug 791)
if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)

3
Eigen/src/SparseCore/SparseFuzzy.h
View File

@ -16,13 +16,12 @@ template<typename Derived>
template<typename OtherDerived>
bool SparseMatrixBase<Derived>::isApprox(const SparseMatrixBase<OtherDerived>& other, const RealScalar &prec) const
{
using std::min;
const typename internal::nested_eval<Derived,2,PlainObject>::type actualA(derived());
typename internal::conditional<bool(IsRowMajor)==bool(OtherDerived::IsRowMajor),
const typename internal::nested_eval<OtherDerived,2,PlainObject>::type,
const PlainObject>::type actualB(other.derived());
return (actualA - actualB).squaredNorm() <= prec * prec * (min)(actualA.squaredNorm(), actualB.squaredNorm());
return (actualA - actualB).squaredNorm() <= prec * prec * numext::mini(actualA.squaredNorm(), actualB.squaredNorm());
}
} // end namespace Eigen

3
Eigen/src/SparseQR/SparseQR.h
View File

@ -325,7 +325,6 @@ template <typename MatrixType, typename OrderingType>
void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
{
using std::abs;
using std::max;
eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
Index m = mat.rows();
@ -377,7 +376,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
if(m_useDefaultThreshold)
{
RealScalar max2Norm = 0.0;
for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm());
for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm());
if(max2Norm==RealScalar(0))
max2Norm = RealScalar(1);
pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();

1
unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h
View File

@ -593,7 +593,6 @@ inline const AutoDiffScalar<Matrix<typename internal::traits<DerTypeA>::Scalar,D
atan2(const AutoDiffScalar<DerTypeA>& a, const AutoDiffScalar<DerTypeB>& b)
{
using std::atan2;
using std::max;
typedef typename internal::traits<DerTypeA>::Scalar Scalar;
typedef AutoDiffScalar<Matrix<Scalar,Dynamic,1> > PlainADS;
PlainADS ret;

8
unsupported/Eigen/src/BDCSVD/BDCSVD.h
View File

@ -649,7 +649,6 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& dia
{
using std::abs;
using std::swap;
using std::max;
Index n = col0.size();
Index actual_n = n;
@ -728,7 +727,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& dia
// rational interpolation: fit a function of the form a / mu + b through the two previous
// iterates and use its zero to compute the next iterate
bool useBisection = fPrev*fCur>0;
while (fCur!=0 && abs(muCur - muPrev) > 8 * NumTraits<RealScalar>::epsilon() * (max)(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection)
while (fCur!=0 && abs(muCur - muPrev) > 8 * NumTraits<RealScalar>::epsilon() * numext::maxi(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection)
{
++m_numIters;
@ -779,7 +778,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayXr& col0, const ArrayXr& dia
#endif
eigen_internal_assert(fLeft * fRight < 0);
while (rightShifted - leftShifted > 2 * NumTraits<RealScalar>::epsilon() * (max)(abs(leftShifted), abs(rightShifted)))
while (rightShifted - leftShifted > 2 * NumTraits<RealScalar>::epsilon() * numext::maxi(abs(leftShifted), abs(rightShifted)))
{
RealScalar midShifted = (leftShifted + rightShifted) / 2;
RealScalar fMid = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
@ -981,7 +980,6 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
{
using std::sqrt;
using std::abs;
using std::max;
const Index length = lastCol + 1 - firstCol;
Block<MatrixXr,Dynamic,1> col0(m_computed, firstCol+shift, firstCol+shift, length, 1);
@ -990,7 +988,7 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
RealScalar maxDiag = diag.tail((std::max)(Index(1),length-1)).cwiseAbs().maxCoeff();
RealScalar epsilon_strict = NumTraits<RealScalar>::epsilon() * maxDiag;
RealScalar epsilon_coarse = 8 * NumTraits<RealScalar>::epsilon() * (max)(col0.cwiseAbs().maxCoeff(), maxDiag);
RealScalar epsilon_coarse = 8 * NumTraits<RealScalar>::epsilon() * numext::maxi(col0.cwiseAbs().maxCoeff(), maxDiag);
#ifdef EIGEN_BDCSVD_SANITY_CHECKS
assert(m_naiveU.allFinite());

3
unsupported/Eigen/src/IterativeSolvers/IncompleteCholesky.h
View File

@ -126,7 +126,6 @@ template<typename _MatrixType>
void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat)
{
using std::sqrt;
using std::min;
eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
// Dropping strategies : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added
@ -160,7 +159,7 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType
for (int j = 0; j < n; j++){
for (int k = colPtr[j]; k < colPtr[j+1]; k++)
vals[k] /= (m_scal(j) * m_scal(rowIdx[k]));
mindiag = (min)(vals[colPtr[j]], mindiag);
mindiag = numext::mini(vals[colPtr[j]], mindiag);
}
if(mindiag < Scalar(0.)) m_shift = m_shift - mindiag;