From 62670c83a0ba7cb4f45a734a4817a818a7c92bba Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 10 Jun 2013 23:40:56 +0200 Subject: Fix bug #314: move remaining math functions from internal to numext namespace --- Eigen/Core | 2 +- Eigen/src/Cholesky/LDLT.h | 18 +- Eigen/src/Cholesky/LLT.h | 10 +- Eigen/src/Core/Dot.h | 4 +- Eigen/src/Core/Functors.h | 16 +- Eigen/src/Core/Fuzzy.h | 4 +- Eigen/src/Core/GeneralProduct.h | 2 +- Eigen/src/Core/GlobalFunctions.h | 3 +- Eigen/src/Core/MathFunctions.h | 218 +++++++++++---------- Eigen/src/Core/SelfAdjointView.h | 12 +- Eigen/src/Core/StableNorm.h | 10 +- Eigen/src/Core/arch/SSE/Complex.h | 8 +- Eigen/src/Core/products/GeneralMatrixVector.h | 2 +- Eigen/src/Core/products/SelfadjointMatrixMatrix.h | 34 ++-- Eigen/src/Core/products/SelfadjointMatrixVector.h | 12 +- Eigen/src/Core/products/SelfadjointRank2Update.h | 12 +- Eigen/src/Core/products/TriangularMatrixVector.h | 2 +- Eigen/src/Core/util/BlasUtil.h | 10 +- Eigen/src/Eigen2Support/MathFunctions.h | 10 +- Eigen/src/Eigen2Support/SVD.h | 8 +- Eigen/src/Eigenvalues/ComplexEigenSolver.h | 2 +- Eigen/src/Eigenvalues/ComplexSchur.h | 10 +- Eigen/src/Eigenvalues/EigenSolver.h | 26 +-- Eigen/src/Eigenvalues/HessenbergDecomposition.h | 2 +- Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 16 +- Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h | 2 +- Eigen/src/Eigenvalues/Tridiagonalization.h | 8 +- Eigen/src/Geometry/OrthoMethods.h | 26 +-- Eigen/src/Householder/Householder.h | 10 +- .../src/IterativeLinearSolvers/ConjugateGradient.h | 4 +- Eigen/src/IterativeLinearSolvers/IncompleteLUT.h | 2 +- Eigen/src/Jacobi/Jacobi.h | 50 ++--- Eigen/src/MetisSupport/MetisSupport.h | 2 +- Eigen/src/QR/ColPivHouseholderQR.h | 2 +- Eigen/src/QR/FullPivHouseholderQR.h | 2 +- Eigen/src/SVD/JacobiSVD.h | 8 +- Eigen/src/SparseCholesky/SimplicialCholesky.h | 4 +- Eigen/src/SparseCholesky/SimplicialCholesky_impl.h | 8 +- Eigen/src/SparseCore/SparseDot.h | 6 +- Eigen/src/SparseCore/SparseSelfAdjointView.h | 11 +- Eigen/src/SparseCore/SparseView.h | 2 +- Eigen/src/SparseQR/SparseQR.h | 12 +- Eigen/src/UmfPackSupport/UmfPackSupport.h | 14 +- 43 files changed, 323 insertions(+), 303 deletions(-) (limited to 'Eigen') diff --git a/Eigen/Core b/Eigen/Core index 8197c944f..7e068cbbd 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -245,8 +245,8 @@ using std::ptrdiff_t; #include "src/Core/util/Constants.h" #include "src/Core/util/ForwardDeclarations.h" #include "src/Core/util/Meta.h" -#include "src/Core/util/XprHelper.h" #include "src/Core/util/StaticAssert.h" +#include "src/Core/util/XprHelper.h" #include "src/Core/util/Memory.h" #include "src/Core/NumTraits.h" diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index d933b10d4..d19cb3968 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -259,7 +259,7 @@ template<> struct ldlt_inplace { transpositions.setIdentity(); if(sign) - *sign = real(mat.coeff(0,0))>0 ? 1:-1; + *sign = numext::real(mat.coeff(0,0))>0 ? 1:-1; return true; } @@ -300,11 +300,11 @@ template<> struct ldlt_inplace for(int i=k+1;i::IsComplex) - mat.coeffRef(index_of_biggest_in_corner,k) = conj(mat.coeff(index_of_biggest_in_corner,k)); + mat.coeffRef(index_of_biggest_in_corner,k) = numext::conj(mat.coeff(index_of_biggest_in_corner,k)); } // partition the matrix: @@ -329,7 +329,7 @@ template<> struct ldlt_inplace if(sign) { // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right - int newSign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0; + int newSign = numext::real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0; if(k == 0) *sign = newSign; else if(*sign != newSign) @@ -350,7 +350,7 @@ template<> struct ldlt_inplace template static bool updateInPlace(MatrixType& mat, MatrixBase& w, const typename MatrixType::RealScalar& sigma=1) { - using internal::isfinite; + using numext::isfinite; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; @@ -368,9 +368,9 @@ template<> struct ldlt_inplace break; // Update the diagonal terms - RealScalar dj = real(mat.coeff(j,j)); + RealScalar dj = numext::real(mat.coeff(j,j)); Scalar wj = w.coeff(j); - RealScalar swj2 = sigma*abs2(wj); + RealScalar swj2 = sigma*numext::abs2(wj); RealScalar gamma = dj*alpha + swj2; mat.coeffRef(j,j) += swj2/alpha; @@ -381,7 +381,7 @@ template<> struct ldlt_inplace Index rs = size-j-1; w.tail(rs) -= wj * mat.col(j).tail(rs); if(gamma != 0) - mat.col(j).tail(rs) += (sigma*conj(wj)/gamma)*w.tail(rs); + mat.col(j).tail(rs) += (sigma*numext::conj(wj)/gamma)*w.tail(rs); } return true; } diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index db22a2f85..2e6189f7d 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -232,10 +232,10 @@ static typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const V RealScalar beta = 1; for(Index j=0; j struct llt_inplace Block A10(mat,k,0,1,k); Block A20(mat,k+1,0,rs,k); - RealScalar x = real(mat.coeff(k,k)); + RealScalar x = numext::real(mat.coeff(k,k)); if (k>0) x -= A10.squaredNorm(); if (x<=RealScalar(0)) return k; diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index 9d513d699..e4107a1ac 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -112,7 +112,7 @@ MatrixBase::eigen2_dot(const MatrixBase& other) const template EIGEN_STRONG_INLINE typename NumTraits::Scalar>::Real MatrixBase::squaredNorm() const { - return internal::real((*this).cwiseAbs2().sum()); + return numext::real((*this).cwiseAbs2().sum()); } /** \returns, for vectors, the \em l2 norm of \c *this, and for matrices the Frobenius norm. @@ -228,7 +228,7 @@ bool MatrixBase::isOrthogonal { typename internal::nested::type nested(derived()); typename internal::nested::type otherNested(other.derived()); - return internal::abs2(nested.dot(otherNested)) <= prec * prec * nested.squaredNorm() * otherNested.squaredNorm(); + return numext::abs2(nested.dot(otherNested)) <= prec * prec * nested.squaredNorm() * otherNested.squaredNorm(); } /** \returns true if *this is approximately an unitary matrix, diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index 9a84e8f26..04fb21732 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -171,7 +171,7 @@ struct functor_traits > { */ template struct scalar_binary_pow_op { EIGEN_EMPTY_STRUCT_CTOR(scalar_binary_pow_op) - inline Scalar operator() (const Scalar& a, const OtherScalar& b) const { return internal::pow(a, b); } + inline Scalar operator() (const Scalar& a, const OtherScalar& b) const { return numext::pow(a, b); } }; template struct functor_traits > { @@ -310,7 +310,7 @@ struct functor_traits > template struct scalar_abs2_op { EIGEN_EMPTY_STRUCT_CTOR(scalar_abs2_op) typedef typename NumTraits::Real result_type; - EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a) const { return internal::abs2(a); } + EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a) const { return numext::abs2(a); } template EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const { return internal::pmul(a,a); } @@ -326,7 +326,7 @@ struct functor_traits > */ template struct scalar_conjugate_op { EIGEN_EMPTY_STRUCT_CTOR(scalar_conjugate_op) - EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const { using internal::conj; return conj(a); } + EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const { using numext::conj; return conj(a); } template EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const { return internal::pconj(a); } }; @@ -363,7 +363,7 @@ template struct scalar_real_op { EIGEN_EMPTY_STRUCT_CTOR(scalar_real_op) typedef typename NumTraits::Real result_type; - EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return internal::real(a); } + EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return numext::real(a); } }; template struct functor_traits > @@ -378,7 +378,7 @@ template struct scalar_imag_op { EIGEN_EMPTY_STRUCT_CTOR(scalar_imag_op) typedef typename NumTraits::Real result_type; - EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return internal::imag(a); } + EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return numext::imag(a); } }; template struct functor_traits > @@ -393,7 +393,7 @@ template struct scalar_real_ref_op { EIGEN_EMPTY_STRUCT_CTOR(scalar_real_ref_op) typedef typename NumTraits::Real result_type; - EIGEN_STRONG_INLINE result_type& operator() (const Scalar& a) const { return internal::real_ref(*const_cast(&a)); } + EIGEN_STRONG_INLINE result_type& operator() (const Scalar& a) const { return numext::real_ref(*const_cast(&a)); } }; template struct functor_traits > @@ -408,7 +408,7 @@ template struct scalar_imag_ref_op { EIGEN_EMPTY_STRUCT_CTOR(scalar_imag_ref_op) typedef typename NumTraits::Real result_type; - EIGEN_STRONG_INLINE result_type& operator() (const Scalar& a) const { return internal::imag_ref(*const_cast(&a)); } + EIGEN_STRONG_INLINE result_type& operator() (const Scalar& a) const { return numext::imag_ref(*const_cast(&a)); } }; template struct functor_traits > @@ -801,7 +801,7 @@ struct scalar_pow_op { // FIXME default copy constructors seems bugged with std::complex<> inline scalar_pow_op(const scalar_pow_op& other) : m_exponent(other.m_exponent) { } inline scalar_pow_op(const Scalar& exponent) : m_exponent(exponent) {} - inline Scalar operator() (const Scalar& a) const { return internal::pow(a, m_exponent); } + inline Scalar operator() (const Scalar& a) const { return numext::pow(a, m_exponent); } const Scalar m_exponent; }; template diff --git a/Eigen/src/Core/Fuzzy.h b/Eigen/src/Core/Fuzzy.h index 8fb9a01dd..fe63bd298 100644 --- a/Eigen/src/Core/Fuzzy.h +++ b/Eigen/src/Core/Fuzzy.h @@ -42,7 +42,7 @@ struct isMuchSmallerThan_object_selector { static bool run(const Derived& x, const OtherDerived& y, const typename Derived::RealScalar& prec) { - return x.cwiseAbs2().sum() <= abs2(prec) * y.cwiseAbs2().sum(); + return x.cwiseAbs2().sum() <= numext::abs2(prec) * y.cwiseAbs2().sum(); } }; @@ -60,7 +60,7 @@ struct isMuchSmallerThan_scalar_selector { static bool run(const Derived& x, const typename Derived::RealScalar& y, const typename Derived::RealScalar& prec) { - return x.cwiseAbs2().sum() <= abs2(prec * y); + return x.cwiseAbs2().sum() <= numext::abs2(prec * y); } }; diff --git a/Eigen/src/Core/GeneralProduct.h b/Eigen/src/Core/GeneralProduct.h index 4e6448353..2a59d9464 100644 --- a/Eigen/src/Core/GeneralProduct.h +++ b/Eigen/src/Core/GeneralProduct.h @@ -435,7 +435,7 @@ template<> struct gemv_selector gemv_static_vector_if static_dest; - bool alphaIsCompatible = (!ComplexByReal) || (imag(actualAlpha)==RealScalar(0)); + bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0)); bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible; RhsScalar compatibleAlpha = get_factor::run(actualAlpha); diff --git a/Eigen/src/Core/GlobalFunctions.h b/Eigen/src/Core/GlobalFunctions.h index 02cae552e..2acf97723 100644 --- a/Eigen/src/Core/GlobalFunctions.h +++ b/Eigen/src/Core/GlobalFunctions.h @@ -39,6 +39,7 @@ namespace Eigen { EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(real,scalar_real_op) EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(imag,scalar_imag_op) + EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(conj,scalar_conjugate_op) EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(sin,scalar_sin_op) EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(cos,scalar_cos_op) EIGEN_ARRAY_DECLARE_GLOBAL_UNARY(asin,scalar_asin_op) @@ -86,6 +87,6 @@ namespace Eigen } } -// TODO: cleanly disable those functions that are not supported on Array (internal::real_ref, internal::random, internal::isApprox...) +// TODO: cleanly disable those functions that are not supported on Array (numext::real_ref, internal::random, internal::isApprox...) #endif // EIGEN_GLOBAL_FUNCTIONS_H diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index a2c55f255..5df2d8bca 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -51,16 +51,15 @@ struct global_math_functions_filtering_base typedef typename T::Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl type; }; -#define EIGEN_MATHFUNC_IMPL(func, scalar) func##_impl::type> -#define EIGEN_MATHFUNC_RETVAL(func, scalar) typename func##_retval::type>::type - +#define EIGEN_MATHFUNC_IMPL(func, scalar) Eigen::internal::func##_impl::type> +#define EIGEN_MATHFUNC_RETVAL(func, scalar) typename Eigen::internal::func##_retval::type>::type /**************************************************************************** * Implementation of real * ****************************************************************************/ -template -struct real_impl +template::IsComplex> +struct real_default_impl { typedef typename NumTraits::Real RealScalar; static inline RealScalar run(const Scalar& x) @@ -69,34 +68,32 @@ struct real_impl } }; -template -struct real_impl > +template +struct real_default_impl { - static inline RealScalar run(const std::complex& x) + typedef typename NumTraits::Real RealScalar; + static inline RealScalar run(const Scalar& x) { using std::real; return real(x); } }; +template struct real_impl : real_default_impl {}; + template struct real_retval { typedef typename NumTraits::Real type; }; -template -inline EIGEN_MATHFUNC_RETVAL(real, Scalar) real(const Scalar& x) -{ - return EIGEN_MATHFUNC_IMPL(real, Scalar)::run(x); -} /**************************************************************************** * Implementation of imag * ****************************************************************************/ -template -struct imag_impl +template::IsComplex> +struct imag_default_impl { typedef typename NumTraits::Real RealScalar; static inline RealScalar run(const Scalar&) @@ -105,28 +102,25 @@ struct imag_impl } }; -template -struct imag_impl > +template +struct imag_default_impl { - static inline RealScalar run(const std::complex& x) + typedef typename NumTraits::Real RealScalar; + static inline RealScalar run(const Scalar& x) { using std::imag; return imag(x); } }; +template struct imag_impl : imag_default_impl {}; + template struct imag_retval { typedef typename NumTraits::Real type; }; -template -inline EIGEN_MATHFUNC_RETVAL(imag, Scalar) imag(const Scalar& x) -{ - return EIGEN_MATHFUNC_IMPL(imag, Scalar)::run(x); -} - /**************************************************************************** * Implementation of real_ref * ****************************************************************************/ @@ -151,18 +145,6 @@ struct real_ref_retval typedef typename NumTraits::Real & type; }; -template -inline typename add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar& x) -{ - return real_ref_impl::run(x); -} - -template -inline EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) real_ref(Scalar& x) -{ - return EIGEN_MATHFUNC_IMPL(real_ref, Scalar)::run(x); -} - /**************************************************************************** * Implementation of imag_ref * ****************************************************************************/ @@ -203,23 +185,11 @@ struct imag_ref_retval typedef typename NumTraits::Real & type; }; -template -inline typename add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) >::type imag_ref(const Scalar& x) -{ - return imag_ref_impl::run(x); -} - -template -inline EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) imag_ref(Scalar& x) -{ - return EIGEN_MATHFUNC_IMPL(imag_ref, Scalar)::run(x); -} - /**************************************************************************** * Implementation of conj * ****************************************************************************/ -template +template::IsComplex> struct conj_impl { static inline Scalar run(const Scalar& x) @@ -228,10 +198,10 @@ struct conj_impl } }; -template -struct conj_impl > +template +struct conj_impl { - static inline std::complex run(const std::complex& x) + static inline Scalar run(const Scalar& x) { using std::conj; return conj(x); @@ -244,12 +214,6 @@ struct conj_retval typedef Scalar type; }; -template -inline EIGEN_MATHFUNC_RETVAL(conj, Scalar) conj(const Scalar& x) -{ - return EIGEN_MATHFUNC_IMPL(conj, Scalar)::run(x); -} - /**************************************************************************** * Implementation of abs2 * ****************************************************************************/ @@ -279,12 +243,6 @@ struct abs2_retval typedef typename NumTraits::Real type; }; -template -inline EIGEN_MATHFUNC_RETVAL(abs2, Scalar) abs2(const Scalar& x) -{ - return EIGEN_MATHFUNC_IMPL(abs2, Scalar)::run(x); -} - /**************************************************************************** * Implementation of norm1 * ****************************************************************************/ @@ -319,12 +277,6 @@ struct norm1_retval typedef typename NumTraits::Real type; }; -template -inline EIGEN_MATHFUNC_RETVAL(norm1, Scalar) norm1(const Scalar& x) -{ - return EIGEN_MATHFUNC_IMPL(norm1, Scalar)::run(x); -} - /**************************************************************************** * Implementation of hypot * ****************************************************************************/ @@ -342,6 +294,7 @@ struct hypot_impl RealScalar _x = abs(x); RealScalar _y = abs(y); RealScalar p = (max)(_x, _y); + if(p==RealScalar(0)) return 0; RealScalar q = (min)(_x, _y); RealScalar qp = q/p; return p * sqrt(RealScalar(1) + qp*qp); @@ -354,12 +307,6 @@ struct hypot_retval typedef typename NumTraits::Real type; }; -template -inline EIGEN_MATHFUNC_RETVAL(hypot, Scalar) hypot(const Scalar& x, const Scalar& y) -{ - return EIGEN_MATHFUNC_IMPL(hypot, Scalar)::run(x, y); -} - /**************************************************************************** * Implementation of cast * ****************************************************************************/ @@ -422,12 +369,6 @@ struct atanh2_retval typedef Scalar type; }; -template -inline EIGEN_MATHFUNC_RETVAL(atanh2, Scalar) atanh2(const Scalar& x, const Scalar& y) -{ - return EIGEN_MATHFUNC_IMPL(atanh2, Scalar)::run(x, y); -} - /**************************************************************************** * Implementation of pow * ****************************************************************************/ @@ -471,12 +412,6 @@ struct pow_retval typedef Scalar type; }; -template -inline EIGEN_MATHFUNC_RETVAL(pow, Scalar) pow(const Scalar& x, const Scalar& y) -{ - return EIGEN_MATHFUNC_IMPL(pow, Scalar)::run(x, y); -} - /**************************************************************************** * Implementation of random * ****************************************************************************/ @@ -611,6 +546,97 @@ inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random() return EIGEN_MATHFUNC_IMPL(random, Scalar)::run(); } +} // end namespace internal + +/**************************************************************************** +* Generic math function * +****************************************************************************/ + +namespace numext { + +template +inline EIGEN_MATHFUNC_RETVAL(real, Scalar) real(const Scalar& x) +{ + return EIGEN_MATHFUNC_IMPL(real, Scalar)::run(x); +} + +template +inline typename internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar& x) +{ + return internal::real_ref_impl::run(x); +} + +template +inline EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) real_ref(Scalar& x) +{ + return EIGEN_MATHFUNC_IMPL(real_ref, Scalar)::run(x); +} + +template +inline EIGEN_MATHFUNC_RETVAL(imag, Scalar) imag(const Scalar& x) +{ + return EIGEN_MATHFUNC_IMPL(imag, Scalar)::run(x); +} + +template +inline typename internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) >::type imag_ref(const Scalar& x) +{ + return internal::imag_ref_impl::run(x); +} + +template +inline EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) imag_ref(Scalar& x) +{ + return EIGEN_MATHFUNC_IMPL(imag_ref, Scalar)::run(x); +} + +template +inline EIGEN_MATHFUNC_RETVAL(conj, Scalar) conj(const Scalar& x) +{ + return EIGEN_MATHFUNC_IMPL(conj, Scalar)::run(x); +} + +template +inline EIGEN_MATHFUNC_RETVAL(abs2, Scalar) abs2(const Scalar& x) +{ + return EIGEN_MATHFUNC_IMPL(abs2, Scalar)::run(x); +} + +template +inline EIGEN_MATHFUNC_RETVAL(norm1, Scalar) norm1(const Scalar& x) +{ + return EIGEN_MATHFUNC_IMPL(norm1, Scalar)::run(x); +} + +template +inline EIGEN_MATHFUNC_RETVAL(hypot, Scalar) hypot(const Scalar& x, const Scalar& y) +{ + return EIGEN_MATHFUNC_IMPL(hypot, Scalar)::run(x, y); +} + +template +inline EIGEN_MATHFUNC_RETVAL(atanh2, Scalar) atanh2(const Scalar& x, const Scalar& y) +{ + return EIGEN_MATHFUNC_IMPL(atanh2, Scalar)::run(x, y); +} + +template +inline EIGEN_MATHFUNC_RETVAL(pow, Scalar) pow(const Scalar& x, const Scalar& y) +{ + return EIGEN_MATHFUNC_IMPL(pow, Scalar)::run(x, y); +} + +// std::isfinite is non standard, so let's define our own version, +// even though it is not very efficient. +template bool (isfinite)(const T& x) +{ + return x::highest() && x>NumTraits::lowest(); +} + +} // end namespace numext + +namespace internal { + /**************************************************************************** * Implementation of fuzzy comparisons * ****************************************************************************/ @@ -668,12 +694,12 @@ struct scalar_fuzzy_default_impl template static inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y, const RealScalar& prec) { - return abs2(x) <= abs2(y) * prec * prec; + return numext::abs2(x) <= numext::abs2(y) * prec * prec; } static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec) { using std::min; - return abs2(x - y) <= (min)(abs2(x), abs2(y)) * prec * prec; + return numext::abs2(x - y) <= (min)(numext::abs2(x), numext::abs2(y)) * prec * prec; } }; @@ -735,17 +761,7 @@ template<> struct scalar_fuzzy_impl }; -/**************************************************************************** -* Special functions * -****************************************************************************/ - -// std::isfinite is non standard, so let's define our own version, -// even though it is not very efficient. -template bool (isfinite)(const T& x) -{ - return x::highest() && x>NumTraits::lowest(); -} - + } // end namespace internal } // end namespace Eigen diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index d43789123..6fa7cd15e 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -214,9 +214,9 @@ struct triangular_assignment_selector::run(dst, src); if(row == col) - dst.coeffRef(row, col) = real(src.coeff(row, col)); + dst.coeffRef(row, col) = numext::real(src.coeff(row, col)); else if(row < col) - dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col)); + dst.coeffRef(col, row) = numext::conj(dst.coeffRef(row, col) = src.coeff(row, col)); } }; @@ -239,9 +239,9 @@ struct triangular_assignment_selector::run(dst, src); if(row == col) - dst.coeffRef(row, col) = real(src.coeff(row, col)); + dst.coeffRef(row, col) = numext::real(src.coeff(row, col)); else if(row > col) - dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col)); + dst.coeffRef(col, row) = numext::conj(dst.coeffRef(row, col) = src.coeff(row, col)); } }; @@ -262,7 +262,7 @@ struct triangular_assignment_selectorscale) { - ssq = ssq * abs2(scale/max); + ssq = ssq * numext::abs2(scale/max); scale = max; invScale = Scalar(1)/scale; } @@ -84,9 +84,9 @@ blueNorm_impl(const EigenBase& _vec) for(typename Derived::InnerIterator it(vec, 0); it; ++it) { RealScalar ax = abs(it.value()); - if(ax > ab2) abig += internal::abs2(ax*s2m); - else if(ax < b1) asml += internal::abs2(ax*s1m); - else amed += internal::abs2(ax); + if(ax > ab2) abig += numext::abs2(ax*s2m); + else if(ax < b1) asml += numext::abs2(ax*s1m); + else amed += numext::abs2(ax); } if(abig > RealScalar(0)) { @@ -120,7 +120,7 @@ blueNorm_impl(const EigenBase& _vec) if(asml <= abig*relerr) return abig; else - return abig * sqrt(RealScalar(1) + internal::abs2(asml/abig)); + return abig * sqrt(RealScalar(1) + numext::abs2(asml/abig)); } } // end namespace internal diff --git a/Eigen/src/Core/arch/SSE/Complex.h b/Eigen/src/Core/arch/SSE/Complex.h index bd76d75ed..91bba5e38 100644 --- a/Eigen/src/Core/arch/SSE/Complex.h +++ b/Eigen/src/Core/arch/SSE/Complex.h @@ -81,8 +81,8 @@ template<> EIGEN_STRONG_INLINE Packet2cf por (const Packet2cf& a, template<> EIGEN_STRONG_INLINE Packet2cf pxor (const Packet2cf& a, const Packet2cf& b) { return Packet2cf(_mm_xor_ps(a.v,b.v)); } template<> EIGEN_STRONG_INLINE Packet2cf pandnot(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(_mm_andnot_ps(a.v,b.v)); } -template<> EIGEN_STRONG_INLINE Packet2cf pload (const std::complex* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet2cf(pload(&real_ref(*from))); } -template<> EIGEN_STRONG_INLINE Packet2cf ploadu(const std::complex* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet2cf(ploadu(&real_ref(*from))); } +template<> EIGEN_STRONG_INLINE Packet2cf pload (const std::complex* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet2cf(pload(&numext::real_ref(*from))); } +template<> EIGEN_STRONG_INLINE Packet2cf ploadu(const std::complex* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet2cf(ploadu(&numext::real_ref(*from))); } template<> EIGEN_STRONG_INLINE Packet2cf pset1(const std::complex& from) { @@ -104,8 +104,8 @@ template<> EIGEN_STRONG_INLINE Packet2cf pset1(const std::complex EIGEN_STRONG_INLINE Packet2cf ploaddup(const std::complex* from) { return pset1(*from); } -template<> EIGEN_STRONG_INLINE void pstore >(std::complex * to, const Packet2cf& from) { EIGEN_DEBUG_ALIGNED_STORE pstore(&real_ref(*to), from.v); } -template<> EIGEN_STRONG_INLINE void pstoreu >(std::complex * to, const Packet2cf& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu(&real_ref(*to), from.v); } +template<> EIGEN_STRONG_INLINE void pstore >(std::complex * to, const Packet2cf& from) { EIGEN_DEBUG_ALIGNED_STORE pstore(&numext::real_ref(*to), from.v); } +template<> EIGEN_STRONG_INLINE void pstoreu >(std::complex * to, const Packet2cf& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu(&numext::real_ref(*to), from.v); } template<> EIGEN_STRONG_INLINE void prefetch >(const std::complex * addr) { _mm_prefetch((const char*)(addr), _MM_HINT_T0); } diff --git a/Eigen/src/Core/products/GeneralMatrixVector.h b/Eigen/src/Core/products/GeneralMatrixVector.h index 9bdd588df..c1cb78498 100644 --- a/Eigen/src/Core/products/GeneralMatrixVector.h +++ b/Eigen/src/Core/products/GeneralMatrixVector.h @@ -86,7 +86,7 @@ EIGEN_DONT_INLINE void general_matrix_vector_product cj; conj_helper pcj; if(ConjugateRhs) - alpha = conj(alpha); + alpha = numext::conj(alpha); enum { AllAligned = 0, EvenAligned, FirstAligned, NoneAligned }; const Index columnsAtOnce = 4; diff --git a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h index ee619df99..99cf9e0ae 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixMatrix.h +++ b/Eigen/src/Core/products/SelfadjointMatrixMatrix.h @@ -30,9 +30,9 @@ struct symm_pack_lhs for(Index k=i; k::IsComplex && EIGEN_LOGICAL_XOR(ConjugateLhs, IsRowMajor), ConjugateRhs> pcj0; conj_helper::IsComplex && EIGEN_LOGICAL_XOR(ConjugateLhs, !IsRowMajor), ConjugateRhs> pcj1; - Scalar cjAlpha = ConjugateRhs ? conj(alpha) : alpha; + Scalar cjAlpha = ConjugateRhs ? numext::conj(alpha) : alpha; // FIXME this copy is now handled outside product_selfadjoint_vector, so it could probably be removed. // if the rhs is not sequentially stored in memory we copy it to a temporary buffer, @@ -98,8 +98,8 @@ EIGEN_DONT_INLINE void selfadjoint_matrix_vector_product huge speed up) // gcc 4.2 does this optimization automatically. @@ -152,7 +152,7 @@ EIGEN_DONT_INLINE void selfadjoint_matrix_vector_product for (Index i=0; i >(mat+stride*i+i, size-i) += - (conj(alpha) * conj(u.coeff(i))) * v.tail(size-i) - + (alpha * conj(v.coeff(i))) * u.tail(size-i); + (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.tail(size-i) + + (alpha * numext::conj(v.coeff(i))) * u.tail(size-i); } } }; @@ -44,8 +44,8 @@ struct selfadjoint_rank2_update_selector const Index size = u.size(); for (Index i=0; i >(mat+stride*i, i+1) += - (conj(alpha) * conj(u.coeff(i))) * v.head(i+1) - + (alpha * conj(v.coeff(i))) * u.head(i+1); + (numext::conj(alpha) * numext::conj(u.coeff(i))) * v.head(i+1) + + (alpha * numext::conj(v.coeff(i))) * u.head(i+1); } }; @@ -75,9 +75,9 @@ SelfAdjointView& SelfAdjointView enum { IsRowMajor = (internal::traits::Flags&RowMajorBit) ? 1 : 0 }; Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived()) - * internal::conj(VBlasTraits::extractScalarFactor(v.derived())); + * numext::conj(VBlasTraits::extractScalarFactor(v.derived())); if (IsRowMajor) - actualAlpha = internal::conj(actualAlpha); + actualAlpha = numext::conj(actualAlpha); internal::selfadjoint_rank2_update_selector::type>::type, diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h index c8b7d28c4..6117d5a82 100644 --- a/Eigen/src/Core/products/TriangularMatrixVector.h +++ b/Eigen/src/Core/products/TriangularMatrixVector.h @@ -245,7 +245,7 @@ template<> struct trmv_selector gemv_static_vector_if static_dest; - bool alphaIsCompatible = (!ComplexByReal) || (imag(actualAlpha)==RealScalar(0)); + bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0)); bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible; RhsScalar compatibleAlpha = get_factor::run(actualAlpha); diff --git a/Eigen/src/Core/util/BlasUtil.h b/Eigen/src/Core/util/BlasUtil.h index 91496651c..a28f16fa0 100644 --- a/Eigen/src/Core/util/BlasUtil.h +++ b/Eigen/src/Core/util/BlasUtil.h @@ -42,7 +42,7 @@ template struct conj_if; template<> struct conj_if { template - inline T operator()(const T& x) { return conj(x); } + inline T operator()(const T& x) { return numext::conj(x); } template inline T pconj(const T& x) { return internal::pconj(x); } }; @@ -67,7 +67,7 @@ template struct conj_helper, std:: { return c + pmul(x,y); } EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const - { return Scalar(real(x)*real(y) + imag(x)*imag(y), imag(x)*real(y) - real(x)*imag(y)); } + { return Scalar(numext::real(x)*numext::real(y) + numext::imag(x)*numext::imag(y), numext::imag(x)*numext::real(y) - numext::real(x)*numext::imag(y)); } }; template struct conj_helper, std::complex, true,false> @@ -77,7 +77,7 @@ template struct conj_helper, std:: { return c + pmul(x,y); } EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const - { return Scalar(real(x)*real(y) + imag(x)*imag(y), real(x)*imag(y) - imag(x)*real(y)); } + { return Scalar(numext::real(x)*numext::real(y) + numext::imag(x)*numext::imag(y), numext::real(x)*numext::imag(y) - numext::imag(x)*numext::real(y)); } }; template struct conj_helper, std::complex, true,true> @@ -87,7 +87,7 @@ template struct conj_helper, std:: { return c + pmul(x,y); } EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const - { return Scalar(real(x)*real(y) - imag(x)*imag(y), - real(x)*imag(y) - imag(x)*real(y)); } + { return Scalar(numext::real(x)*numext::real(y) - numext::imag(x)*numext::imag(y), - numext::real(x)*numext::imag(y) - numext::imag(x)*numext::real(y)); } }; template struct conj_helper, RealScalar, Conj,false> @@ -113,7 +113,7 @@ template struct get_factor { }; template struct get_factor::Real> { - static EIGEN_STRONG_INLINE typename NumTraits::Real run(const Scalar& x) { return real(x); } + static EIGEN_STRONG_INLINE typename NumTraits::Real run(const Scalar& x) { return numext::real(x); } }; // Lightweight helper class to access matrix coefficients. diff --git a/Eigen/src/Eigen2Support/MathFunctions.h b/Eigen/src/Eigen2Support/MathFunctions.h index bde5dd441..3544af253 100644 --- a/Eigen/src/Eigen2Support/MathFunctions.h +++ b/Eigen/src/Eigen2Support/MathFunctions.h @@ -12,18 +12,18 @@ namespace Eigen { -template inline typename NumTraits::Real ei_real(const T& x) { return internal::real(x); } -template inline typename NumTraits::Real ei_imag(const T& x) { return internal::imag(x); } -template inline T ei_conj(const T& x) { return internal::conj(x); } +template inline typename NumTraits::Real ei_real(const T& x) { return numext::real(x); } +template inline typename NumTraits::Real ei_imag(const T& x) { return numext::imag(x); } +template inline T ei_conj(const T& x) { return numext::conj(x); } template inline typename NumTraits::Real ei_abs (const T& x) { using std::abs; return abs(x); } -template inline typename NumTraits::Real ei_abs2(const T& x) { return internal::abs2(x); } +template inline typename NumTraits::Real ei_abs2(const T& x) { return numext::abs2(x); } template inline T ei_sqrt(const T& x) { using std::sqrt; return sqrt(x); } template inline T ei_exp (const T& x) { using std::exp; return exp(x); } template inline T ei_log (const T& x) { using std::log; return log(x); } template inline T ei_sin (const T& x) { using std::sin; return sin(x); } template inline T ei_cos (const T& x) { using std::cos; return cos(x); } template inline T ei_atan2(const T& x,const T& y) { using std::atan2; return atan2(x,y); } -template inline T ei_pow (const T& x,const T& y) { return internal::pow(x,y); } +template inline T ei_pow (const T& x,const T& y) { return numext::pow(x,y); } template inline T ei_random () { return internal::random(); } template inline T ei_random (const T& x, const T& y) { return internal::random(x, y); } diff --git a/Eigen/src/Eigen2Support/SVD.h b/Eigen/src/Eigen2Support/SVD.h index a08b695a4..077d26d54 100644 --- a/Eigen/src/Eigen2Support/SVD.h +++ b/Eigen/src/Eigen2Support/SVD.h @@ -315,7 +315,7 @@ void SVD::compute(const MatrixType& matrix) e[p-2] = 0.0; for (j = p-2; j >= k; --j) { - Scalar t(internal::hypot(m_sigma[j],f)); + Scalar t(numext::hypot(m_sigma[j],f)); Scalar cs(m_sigma[j]/t); Scalar sn(f/t); m_sigma[j] = t; @@ -344,7 +344,7 @@ void SVD::compute(const MatrixType& matrix) e[k-1] = 0.0; for (j = k; j < p; ++j) { - Scalar t(internal::hypot(m_sigma[j],f)); + Scalar t(numext::hypot(m_sigma[j],f)); Scalar cs( m_sigma[j]/t); Scalar sn(f/t); m_sigma[j] = t; @@ -392,7 +392,7 @@ void SVD::compute(const MatrixType& matrix) for (j = k; j < p-1; ++j) { - Scalar t = internal::hypot(f,g); + Scalar t = numext::hypot(f,g); Scalar cs = f/t; Scalar sn = g/t; if (j != k) @@ -410,7 +410,7 @@ void SVD::compute(const MatrixType& matrix) m_matV(i,j) = t; } } - t = internal::hypot(f,g); + t = numext::hypot(f,g); cs = f/t; sn = g/t; m_sigma[j] = t; diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h index bd41bf7ed..af434bc9b 100644 --- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h +++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h @@ -294,7 +294,7 @@ void ComplexEigenSolver::doComputeEigenvectors(const RealScalar& mat { // If the i-th and k-th eigenvalue are equal, then z equals 0. // Use a small value instead, to prevent division by zero. - internal::real_ref(z) = NumTraits::epsilon() * matrixnorm; + numext::real_ref(z) = NumTraits::epsilon() * matrixnorm; } m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z; } diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h index 62b57ff66..89e6cade3 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/Eigen/src/Eigenvalues/ComplexSchur.h @@ -263,8 +263,8 @@ template class ComplexSchur template inline bool ComplexSchur::subdiagonalEntryIsNeglegible(Index i) { - RealScalar d = internal::norm1(m_matT.coeff(i,i)) + internal::norm1(m_matT.coeff(i+1,i+1)); - RealScalar sd = internal::norm1(m_matT.coeff(i+1,i)); + RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1)); + RealScalar sd = numext::norm1(m_matT.coeff(i+1,i)); if (internal::isMuchSmallerThan(sd, d, NumTraits::epsilon())) { m_matT.coeffRef(i+1,i) = ComplexScalar(0); @@ -282,7 +282,7 @@ typename ComplexSchur::ComplexScalar ComplexSchur::compu if (iter == 10 || iter == 20) { // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f - return abs(internal::real(m_matT.coeff(iu,iu-1))) + abs(internal::real(m_matT.coeff(iu-1,iu-2))); + return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2))); } // compute the shift as one of the eigenvalues of t, the 2x2 @@ -299,13 +299,13 @@ typename ComplexSchur::ComplexScalar ComplexSchur::compu ComplexScalar eival1 = (trace + disc) / RealScalar(2); ComplexScalar eival2 = (trace - disc) / RealScalar(2); - if(internal::norm1(eival1) > internal::norm1(eival2)) + if(numext::norm1(eival1) > numext::norm1(eival2)) eival2 = det / eival1; else eival1 = det / eival2; // choose the eigenvalue closest to the bottom entry of the diagonal - if(internal::norm1(eival1-t.coeff(1,1)) < internal::norm1(eival2-t.coeff(1,1))) + if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1))) return normt * eival1; else return normt * eival2; diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index f0d4e5afa..6e7150685 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -317,12 +317,12 @@ MatrixType EigenSolver::pseudoEigenvalueMatrix() const MatrixType matD = MatrixType::Zero(n,n); for (Index i=0; i(i,i) << internal::real(m_eivalues.coeff(i)), internal::imag(m_eivalues.coeff(i)), - -internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i)); + matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)), + -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)); ++i; } } @@ -338,7 +338,7 @@ typename EigenSolver::EigenvectorsType EigenSolver::eige EigenvectorsType matV(n,n); for (Index j=0; j(); @@ -515,8 +515,8 @@ void EigenSolver::doComputeEigenvectors() else { std::complex cc = cdiv(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q); - m_matT.coeffRef(n-1,n-1) = internal::real(cc); - m_matT.coeffRef(n-1,n) = internal::imag(cc); + m_matT.coeffRef(n-1,n-1) = numext::real(cc); + m_matT.coeffRef(n-1,n) = numext::imag(cc); } m_matT.coeffRef(n,n-1) = 0.0; m_matT.coeffRef(n,n) = 1.0; @@ -538,8 +538,8 @@ void EigenSolver::doComputeEigenvectors() if (m_eivalues.coeff(i).imag() == RealScalar(0)) { std::complex cc = cdiv(-ra,-sa,w,q); - m_matT.coeffRef(i,n-1) = internal::real(cc); - m_matT.coeffRef(i,n) = internal::imag(cc); + m_matT.coeffRef(i,n-1) = numext::real(cc); + m_matT.coeffRef(i,n) = numext::imag(cc); } else { @@ -552,8 +552,8 @@ void EigenSolver::doComputeEigenvectors() vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw)); std::complex cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi); - m_matT.coeffRef(i,n-1) = internal::real(cc); - m_matT.coeffRef(i,n) = internal::imag(cc); + m_matT.coeffRef(i,n-1) = numext::real(cc); + m_matT.coeffRef(i,n) = numext::imag(cc); if (abs(x) > (abs(lastw) + abs(q))) { m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x; @@ -562,8 +562,8 @@ void EigenSolver::doComputeEigenvectors() else { cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q); - m_matT.coeffRef(i+1,n-1) = internal::real(cc); - m_matT.coeffRef(i+1,n) = internal::imag(cc); + m_matT.coeffRef(i+1,n-1) = numext::real(cc); + m_matT.coeffRef(i+1,n) = numext::imag(cc); } } diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h index ebd8ae908..afa636ffa 100644 --- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h +++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -313,7 +313,7 @@ void HessenbergDecomposition::_compute(MatrixType& matA, CoeffVector // A = A H' matA.rightCols(remainingSize) - .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), internal::conj(h), &temp.coeffRef(0)); + .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0)); } } diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 03c024927..3993046a8 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -395,7 +395,7 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver if(n==1) { - m_eivalues.coeffRef(0,0) = internal::real(matrix.coeff(0,0)); + m_eivalues.coeffRef(0,0) = numext::real(matrix.coeff(0,0)); if(computeEigenvectors) m_eivec.setOnes(n,n); m_info = Success; @@ -669,7 +669,7 @@ template struct direct_selfadjoint_eigenvalues struct direct_selfadjoint_eigenvaluesc2) { eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0); @@ -744,7 +744,7 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta RealScalar e = subdiag[end-1]; // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still // underflow thus leading to inf/NaN values when using the following commented code: -// RealScalar e2 = abs2(subdiag[end-1]); +// RealScalar e2 = numext::abs2(subdiag[end-1]); // RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2)); // This explain the following, somewhat more complicated, version: RealScalar mu = diag[end]; @@ -752,8 +752,8 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta mu -= abs(e); else { - RealScalar e2 = abs2(subdiag[end-1]); - RealScalar h = hypot(td,e); + RealScalar e2 = numext::abs2(subdiag[end-1]); + RealScalar h = numext::hypot(td,e); if(e2==0) mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h); else mu -= e2 / (td + (td>0 ? h : -h)); } diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h index 5de5f15d6..17c0dadd2 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h @@ -56,7 +56,7 @@ SelfAdjointEigenSolver >::compute(c \ if(n==1) \ { \ - m_eivalues.coeffRef(0,0) = internal::real(matrix.coeff(0,0)); \ + m_eivalues.coeffRef(0,0) = numext::real(matrix.coeff(0,0)); \ if(computeEigenvectors) m_eivec.setOnes(n,n); \ m_info = Success; \ m_isInitialized = true; \ diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h index e8408761d..b55a01d3a 100644 --- a/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -345,7 +345,7 @@ namespace internal { template void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs) { - using internal::conj; + using numext::conj; typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; @@ -468,7 +468,7 @@ struct tridiagonalization_inplace_selector { using std::sqrt; diag[0] = mat(0,0); - RealScalar v1norm2 = abs2(mat(2,0)); + RealScalar v1norm2 = numext::abs2(mat(2,0)); if(v1norm2 == RealScalar(0)) { diag[1] = mat(1,1); @@ -480,7 +480,7 @@ struct tridiagonalization_inplace_selector } else { - RealScalar beta = sqrt(abs2(mat(1,0)) + v1norm2); + RealScalar beta = sqrt(numext::abs2(mat(1,0)) + v1norm2); RealScalar invBeta = RealScalar(1)/beta; Scalar m01 = mat(1,0) * invBeta; Scalar m02 = mat(2,0) * invBeta; @@ -510,7 +510,7 @@ struct tridiagonalization_inplace_selector template static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ) { - diag(0,0) = real(mat(0,0)); + diag(0,0) = numext::real(mat(0,0)); if(extractQ) mat(0,0) = Scalar(1); } diff --git a/Eigen/src/Geometry/OrthoMethods.h b/Eigen/src/Geometry/OrthoMethods.h index 4c1bf5fcd..556bc8160 100644 --- a/Eigen/src/Geometry/OrthoMethods.h +++ b/Eigen/src/Geometry/OrthoMethods.h @@ -33,9 +33,9 @@ MatrixBase::cross(const MatrixBase& other) const typename internal::nested::type lhs(derived()); typename internal::nested::type rhs(other.derived()); return typename cross_product_return_type::type( - internal::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)), - internal::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)), - internal::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)) + numext::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)), + numext::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)), + numext::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)) ); } @@ -49,9 +49,9 @@ struct cross3_impl { run(const VectorLhs& lhs, const VectorRhs& rhs) { return typename internal::plain_matrix_type::type( - internal::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)), - internal::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)), - internal::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)), + numext::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)), + numext::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)), + numext::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0)), 0 ); } @@ -141,8 +141,8 @@ struct unitOrthogonal_selector if (maxi==0) sndi = 1; RealScalar invnm = RealScalar(1)/(Vector2() << src.coeff(sndi),src.coeff(maxi)).finished().norm(); - perp.coeffRef(maxi) = -conj(src.coeff(sndi)) * invnm; - perp.coeffRef(sndi) = conj(src.coeff(maxi)) * invnm; + perp.coeffRef(maxi) = -numext::conj(src.coeff(sndi)) * invnm; + perp.coeffRef(sndi) = numext::conj(src.coeff(maxi)) * invnm; return perp; } @@ -168,8 +168,8 @@ struct unitOrthogonal_selector || (!isMuchSmallerThan(src.y(), src.z()))) { RealScalar invnm = RealScalar(1)/src.template head<2>().norm(); - perp.coeffRef(0) = -conj(src.y())*invnm; - perp.coeffRef(1) = conj(src.x())*invnm; + perp.coeffRef(0) = -numext::conj(src.y())*invnm; + perp.coeffRef(1) = numext::conj(src.x())*invnm; perp.coeffRef(2) = 0; } /* if both x and y are close to zero, then the vector is close @@ -180,8 +180,8 @@ struct unitOrthogonal_selector { RealScalar invnm = RealScalar(1)/src.template tail<2>().norm(); perp.coeffRef(0) = 0; - perp.coeffRef(1) = -conj(src.z())*invnm; - perp.coeffRef(2) = conj(src.y())*invnm; + perp.coeffRef(1) = -numext::conj(src.z())*invnm; + perp.coeffRef(2) = numext::conj(src.y())*invnm; } return perp; @@ -193,7 +193,7 @@ struct unitOrthogonal_selector { typedef typename plain_matrix_type::type VectorType; static inline VectorType run(const Derived& src) - { return VectorType(-conj(src.y()), conj(src.x())).normalized(); } + { return VectorType(-numext::conj(src.y()), numext::conj(src.x())).normalized(); } }; } // end namespace internal diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index b7cfa9b2b..32112af9b 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -68,7 +68,7 @@ void MatrixBase::makeHouseholder( RealScalar& beta) const { using std::sqrt; - using internal::conj; + using numext::conj; EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) VectorBlock tail(derived(), 1, size()-1); @@ -76,16 +76,16 @@ void MatrixBase::makeHouseholder( RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm(); Scalar c0 = coeff(0); - if(tailSqNorm == RealScalar(0) && internal::imag(c0)==RealScalar(0)) + if(tailSqNorm == RealScalar(0) && numext::imag(c0)==RealScalar(0)) { tau = RealScalar(0); - beta = internal::real(c0); + beta = numext::real(c0); essential.setZero(); } else { - beta = sqrt(internal::abs2(c0) + tailSqNorm); - if (internal::real(c0)>=RealScalar(0)) + beta = sqrt(numext::abs2(c0) + tailSqNorm); + if (numext::real(c0)>=RealScalar(0)) beta = -beta; essential = tail / (c0 - beta); tau = conj((beta - c0) / beta); diff --git a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h index 00b5647c6..a74a8155e 100644 --- a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h +++ b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h @@ -63,7 +63,7 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, p = precond.solve(residual); //initial search direction VectorType z(n), tmp(n); - RealScalar absNew = internal::real(residual.dot(p)); // the square of the absolute value of r scaled by invM + RealScalar absNew = numext::real(residual.dot(p)); // the square of the absolute value of r scaled by invM int i = 0; while(i < maxIters) { @@ -80,7 +80,7 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, z = precond.solve(residual); // approximately solve for "A z = residual" RealScalar absOld = absNew; - absNew = internal::real(residual.dot(z)); // update the absolute value of r + absNew = numext::real(residual.dot(z)); // update the absolute value of r RealScalar beta = absNew / absOld; // calculate the Gram-Schmidt value used to create the new search direction p = z + beta * p; // update search direction i++; diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 17d18ef58..50a870aec 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h @@ -310,7 +310,7 @@ void IncompleteLUT::factorize(const _MatrixType& amat) jr(k) = jpos; ++sizeu; } - rownorm += internal::abs2(j_it.value()); + rownorm += numext::abs2(j_it.value()); } // 2 - detect possible zero row diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index d9d75196c..956f72d57 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -50,16 +50,16 @@ template class JacobiRotation /** Concatenates two planar rotation */ JacobiRotation operator*(const JacobiRotation& other) { - using internal::conj; + using numext::conj; return JacobiRotation(m_c * other.m_c - conj(m_s) * other.m_s, conj(m_c * conj(other.m_s) + conj(m_s) * conj(other.m_c))); } /** Returns the transposed transformation */ - JacobiRotation transpose() const { using internal::conj; return JacobiRotation(m_c, -conj(m_s)); } + JacobiRotation transpose() const { using numext::conj; return JacobiRotation(m_c, -conj(m_s)); } /** Returns the adjoint transformation */ - JacobiRotation adjoint() const { using internal::conj; return JacobiRotation(conj(m_c), -m_s); } + JacobiRotation adjoint() const { using numext::conj; return JacobiRotation(conj(m_c), -m_s); } template bool makeJacobi(const MatrixBase&, typename Derived::Index p, typename Derived::Index q); @@ -94,7 +94,7 @@ bool JacobiRotation::makeJacobi(const RealScalar& x, const Scalar& y, co else { RealScalar tau = (x-z)/(RealScalar(2)*abs(y)); - RealScalar w = sqrt(internal::abs2(tau) + RealScalar(1)); + RealScalar w = sqrt(numext::abs2(tau) + RealScalar(1)); RealScalar t; if(tau>RealScalar(0)) { @@ -105,8 +105,8 @@ bool JacobiRotation::makeJacobi(const RealScalar& x, const Scalar& y, co t = RealScalar(1) / (tau - w); } RealScalar sign_t = t > RealScalar(0) ? RealScalar(1) : RealScalar(-1); - RealScalar n = RealScalar(1) / sqrt(internal::abs2(t)+RealScalar(1)); - m_s = - sign_t * (internal::conj(y) / abs(y)) * abs(t) * n; + RealScalar n = RealScalar(1) / sqrt(numext::abs2(t)+RealScalar(1)); + m_s = - sign_t * (numext::conj(y) / abs(y)) * abs(t) * n; m_c = n; return true; } @@ -125,7 +125,7 @@ template template inline bool JacobiRotation::makeJacobi(const MatrixBase& m, typename Derived::Index p, typename Derived::Index q) { - return makeJacobi(internal::real(m.coeff(p,p)), m.coeff(p,q), internal::real(m.coeff(q,q))); + return makeJacobi(numext::real(m.coeff(p,p)), m.coeff(p,q), numext::real(m.coeff(q,q))); } /** Makes \c *this as a Givens rotation \c G such that applying \f$ G^* \f$ to the left of the vector @@ -157,11 +157,11 @@ void JacobiRotation::makeGivens(const Scalar& p, const Scalar& q, Scalar { using std::sqrt; using std::abs; - using internal::conj; + using numext::conj; if(q==Scalar(0)) { - m_c = internal::real(p)<0 ? Scalar(-1) : Scalar(1); + m_c = numext::real(p)<0 ? Scalar(-1) : Scalar(1); m_s = 0; if(r) *r = m_c * p; } @@ -173,17 +173,17 @@ void JacobiRotation::makeGivens(const Scalar& p, const Scalar& q, Scalar } else { - RealScalar p1 = internal::norm1(p); - RealScalar q1 = internal::norm1(q); + RealScalar p1 = numext::norm1(p); + RealScalar q1 = numext::norm1(q); if(p1>=q1) { Scalar ps = p / p1; - RealScalar p2 = internal::abs2(ps); + RealScalar p2 = numext::abs2(ps); Scalar qs = q / p1; - RealScalar q2 = internal::abs2(qs); + RealScalar q2 = numext::abs2(qs); RealScalar u = sqrt(RealScalar(1) + q2/p2); - if(internal::real(p)::makeGivens(const Scalar& p, const Scalar& q, Scalar else { Scalar ps = p / q1; - RealScalar p2 = internal::abs2(ps); + RealScalar p2 = numext::abs2(ps); Scalar qs = q / q1; - RealScalar q2 = internal::abs2(qs); + RealScalar q2 = numext::abs2(qs); RealScalar u = q1 * sqrt(p2 + q2); - if(internal::real(p)::makeGivens(const Scalar& p, const Scalar& q, Scalar else if(abs(p) > abs(q)) { Scalar t = q/p; - Scalar u = sqrt(Scalar(1) + internal::abs2(t)); + Scalar u = sqrt(Scalar(1) + numext::abs2(t)); if(p::makeGivens(const Scalar& p, const Scalar& q, Scalar else { Scalar t = p/q; - Scalar u = sqrt(Scalar(1) + internal::abs2(t)); + Scalar u = sqrt(Scalar(1) + numext::abs2(t)); if(q& ColPivHouseholderQR::compute(const for(Index k = 0; k < cols; ++k) m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); - RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits::epsilon()) / RealScalar(rows); + RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits::epsilon()) / RealScalar(rows); m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) m_maxpivot = RealScalar(0); diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index f82126494..0dd5ad347 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -572,7 +572,7 @@ public: template void evalTo(ResultType& result, WorkVectorType& workspace) const { - using internal::conj; + using numext::conj; // compute the product H'_0 H'_1 ... H'_n-1, // where H_k is the k-th Householder transformation I - h_k v_k v_k' // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 495d3fabf..c7a7eeda0 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -374,7 +374,7 @@ struct svd_precondition_2x2_block_to_be_real using std::sqrt; Scalar z; JacobiRotation rot; - RealScalar n = sqrt(abs2(work_matrix.coeff(p,p)) + abs2(work_matrix.coeff(q,p))); + RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p))); if(n==0) { z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); @@ -413,8 +413,8 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, { using std::sqrt; Matrix m; - m << real(matrix.coeff(p,p)), real(matrix.coeff(p,q)), - real(matrix.coeff(q,p)), real(matrix.coeff(q,q)); + m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), + numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); JacobiRotation rot1; RealScalar t = m.coeff(0,0) + m.coeff(1,1); RealScalar d = m.coeff(1,0) - m.coeff(0,1); @@ -426,7 +426,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, else { RealScalar u = d / t; - rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + abs2(u)); + rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u)); rot1.s() = rot1.c() * u; } m.applyOnTheLeft(0,1,rot1); diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h index 62747279d..f41d7e010 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h @@ -364,7 +364,7 @@ public: Scalar determinant() const { Scalar detL = Base::m_matrix.diagonal().prod(); - return internal::abs2(detL); + return numext::abs2(detL); } }; @@ -599,7 +599,7 @@ public: else { Scalar detL = Diagonal(Base::m_matrix).prod(); - return internal::abs2(detL); + return numext::abs2(detL); } } diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h b/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h index 4b249868f..7aaf702be 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky_impl.h @@ -131,7 +131,7 @@ void SimplicialCholeskyBase::factorize_preordered(const CholMatrixType& Index i = it.index(); if(i <= k) { - y[i] += internal::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */ + y[i] += numext::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */ Index len; for(len = 0; tags[i] != k; i = m_parent[i]) { @@ -145,7 +145,7 @@ void SimplicialCholeskyBase::factorize_preordered(const CholMatrixType& /* compute numerical values kth row of L (a sparse triangular solve) */ - RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k) + RealScalar d = numext::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k) y[k] = 0.0; for(; top < size; ++top) { @@ -163,8 +163,8 @@ void SimplicialCholeskyBase::factorize_preordered(const CholMatrixType& Index p2 = Lp[i] + m_nonZerosPerCol[i]; Index p; for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p) - y[Li[p]] -= internal::conj(Lx[p]) * yi; - d -= internal::real(l_ki * internal::conj(yi)); + y[Li[p]] -= numext::conj(Lx[p]) * yi; + d -= numext::real(l_ki * numext::conj(yi)); Li[p] = k; /* store L(k,i) in column form of L */ Lx[p] = l_ki; ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */ diff --git a/Eigen/src/SparseCore/SparseDot.h b/Eigen/src/SparseCore/SparseDot.h index dfeb3a8df..db39c9aec 100644 --- a/Eigen/src/SparseCore/SparseDot.h +++ b/Eigen/src/SparseCore/SparseDot.h @@ -30,7 +30,7 @@ SparseMatrixBase::dot(const MatrixBase& other) const Scalar res(0); while (i) { - res += internal::conj(i.value()) * other.coeff(i.index()); + res += numext::conj(i.value()) * other.coeff(i.index()); ++i; } return res; @@ -64,7 +64,7 @@ SparseMatrixBase::dot(const SparseMatrixBase& other) cons { if (i.index()==j.index()) { - res += internal::conj(i.value()) * j.value(); + res += numext::conj(i.value()) * j.value(); ++i; ++j; } else if (i.index() inline typename NumTraits::Scalar>::Real SparseMatrixBase::squaredNorm() const { - return internal::real((*this).cwiseAbs2().sum()); + return numext::real((*this).cwiseAbs2().sum()); } template diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 9630b60f5..d2e170410 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -240,7 +240,7 @@ class SparseSelfAdjointTimeDenseProduct Index b = LhsIsRowMajor ? i.index() : j; typename Lhs::Scalar v = i.value(); dest.row(a) += (v) * m_rhs.row(b); - dest.row(b) += internal::conj(v) * m_rhs.row(a); + dest.row(b) += numext::conj(v) * m_rhs.row(a); } if (ProcessFirstHalf && i && (i.index()==j)) dest.row(j) += i.value() * m_rhs.row(j); @@ -367,7 +367,7 @@ void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrixjp))) - dest.valuePtr()[k] = conj(it.value()); + dest.valuePtr()[k] = numext::conj(it.value()); else dest.valuePtr()[k] = it.value(); } @@ -461,7 +461,10 @@ class SparseSymmetricPermutationProduct template void evalTo(SparseMatrix& _dest) const { - internal::permute_symm_to_fullsymm(m_matrix,_dest,m_perm.indices().data()); +// internal::permute_symm_to_fullsymm(m_matrix,_dest,m_perm.indices().data()); + SparseMatrix tmp; + internal::permute_symm_to_fullsymm(m_matrix,tmp,m_perm.indices().data()); + _dest = tmp; } template void evalTo(SparseSelfAdjointView& dest) const diff --git a/Eigen/src/SparseCore/SparseView.h b/Eigen/src/SparseCore/SparseView.h index 4fd0cb3d8..67eb93245 100644 --- a/Eigen/src/SparseCore/SparseView.h +++ b/Eigen/src/SparseCore/SparseView.h @@ -18,7 +18,7 @@ namespace internal { template struct traits > : traits { - typedef int Index; + typedef typename MatrixType::Index Index; typedef Sparse StorageKind; enum { Flags = int(traits::Flags) & (RowMajorBit) diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index 1bf631507..07c46e4b9 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -435,23 +435,23 @@ void SparseQR::factorize(const MatrixType& mat) // First, the squared norm of Q((col+1):m, col) RealScalar sqrNorm = 0.; - for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += internal::abs2(tval(Qidx(itq))); + for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq))); - if(sqrNorm == RealScalar(0) && internal::imag(c0) == RealScalar(0)) + if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) { tau = RealScalar(0); - beta = internal::real(c0); + beta = numext::real(c0); tval(Qidx(0)) = 1; } else { - beta = std::sqrt(internal::abs2(c0) + sqrNorm); - if(internal::real(c0) >= RealScalar(0)) + beta = std::sqrt(numext::abs2(c0) + sqrNorm); + if(numext::real(c0) >= RealScalar(0)) beta = -beta; tval(Qidx(0)) = 1; for (Index itq = 1; itq < nzcolQ; ++itq) tval(Qidx(itq)) /= (c0 - beta); - tau = internal::conj((beta-c0) / beta); + tau = numext::conj((beta-c0) / beta); } diff --git a/Eigen/src/UmfPackSupport/UmfPackSupport.h b/Eigen/src/UmfPackSupport/UmfPackSupport.h index d85b8be85..3a48cecf7 100644 --- a/Eigen/src/UmfPackSupport/UmfPackSupport.h +++ b/Eigen/src/UmfPackSupport/UmfPackSupport.h @@ -39,7 +39,7 @@ inline int umfpack_symbolic(int n_row,int n_col, const int Ap[], const int Ai[], const std::complex Ax[], void **Symbolic, const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) { - return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&internal::real_ref(Ax[0]),0,Symbolic,Control,Info); + return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info); } inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[], @@ -53,7 +53,7 @@ inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex X[], const std::complex B[], void *Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) { - return umfpack_zi_solve(sys,Ap,Ai,&internal::real_ref(Ax[0]),0,&internal::real_ref(X[0]),0,&internal::real_ref(B[0]),0,Numeric,Control,Info); + return umfpack_zi_solve(sys,Ap,Ai,&numext::real_ref(Ax[0]),0,&numext::real_ref(X[0]),0,&numext::real_ref(B[0]),0,Numeric,Control,Info); } inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double) @@ -89,9 +89,9 @@ inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex Lx[], int Up[], int Ui[], std::complex Ux[], int P[], int Q[], std::complex Dx[], int *do_recip, double Rs[], void *Numeric) { - double& lx0_real = internal::real_ref(Lx[0]); - double& ux0_real = internal::real_ref(Ux[0]); - double& dx0_real = internal::real_ref(Dx[0]); + double& lx0_real = numext::real_ref(Lx[0]); + double& ux0_real = numext::real_ref(Ux[0]); + double& dx0_real = numext::real_ref(Dx[0]); return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q, Dx?&dx0_real:0,0,do_recip,Rs,Numeric); } @@ -103,7 +103,7 @@ inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, inline int umfpack_get_determinant(std::complex *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO]) { - double& mx_real = internal::real_ref(*Mx); + double& mx_real = numext::real_ref(*Mx); return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info); } -- cgit v1.2.3