From 014d9f1d9b60206deaeb7ac5349816cb556fb35b Mon Sep 17 00:00:00 2001 From: Hongkai Dai Date: Thu, 13 Oct 2016 14:45:51 -0700 Subject: [PATCH 01/22] implement euler angles with the right ranges --- .../Eigen/src/EulerAngles/EulerAngles.h | 26 ++-- .../Eigen/src/EulerAngles/EulerSystem.h | 137 +++++++++--------- unsupported/test/EulerAngles.cpp | 77 +++++----- 3 files changed, 113 insertions(+), 127 deletions(-) diff --git a/unsupported/Eigen/src/EulerAngles/EulerAngles.h b/unsupported/Eigen/src/EulerAngles/EulerAngles.h index 13a0da1ab..a737a221a 100644 --- a/unsupported/Eigen/src/EulerAngles/EulerAngles.h +++ b/unsupported/Eigen/src/EulerAngles/EulerAngles.h @@ -79,8 +79,8 @@ namespace Eigen * * ##### run-time time ranges ##### * Run-time ranges are also supported. - * \sa EulerAngles(const MatrixBase&, bool, bool, bool) - * \sa EulerAngles(const RotationBase&, bool, bool, bool) + * \sa EulerAngles(const MatrixBase&, bool, bool) + * \sa EulerAngles(const RotationBase&, bool, bool) * * ### Convenient user typedefs ### * @@ -160,22 +160,24 @@ namespace Eigen /** Constructs and initialize Euler angles from a 3x3 rotation matrix \p m, * with options to choose for each angle the requested range. * - * If positive range is true, then the specified angle will be in the range [0, +2*PI]. + * For angle alpha and gamma, if positive range is true, then the + * specified angle will be in the range [0, +2*PI]. * Otherwise, the specified angle will be in the range [-PI, +PI]. + * For angle beta, depending on whether AlphaAxis is the same as GammaAxis + * if AlphaAxis is the same as Gamma ais, then the range of beta is [0, PI]; + * otherwise the range of beta is [-PI/2, PI/2] * * \param m The 3x3 rotation matrix to convert * \param positiveRangeAlpha If true, alpha will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. - * \param positiveRangeBeta If true, beta will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. * \param positiveRangeGamma If true, gamma will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. */ template EulerAngles( const MatrixBase& m, bool positiveRangeAlpha, - bool positiveRangeBeta, bool positiveRangeGamma) { - System::CalcEulerAngles(*this, m, positiveRangeAlpha, positiveRangeBeta, positiveRangeGamma); + System::CalcEulerAngles(*this, m, positiveRangeAlpha, positiveRangeGamma); } /** Constructs and initialize Euler angles from a rotation \p rot. @@ -195,17 +197,15 @@ namespace Eigen * * \param rot The 3x3 rotation matrix to convert * \param positiveRangeAlpha If true, alpha will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. - * \param positiveRangeBeta If true, beta will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. * \param positiveRangeGamma If true, gamma will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. */ template EulerAngles( const RotationBase& rot, bool positiveRangeAlpha, - bool positiveRangeBeta, bool positiveRangeGamma) { - System::CalcEulerAngles(*this, rot.toRotationMatrix(), positiveRangeAlpha, positiveRangeBeta, positiveRangeGamma); + System::CalcEulerAngles(*this, rot.toRotationMatrix(), positiveRangeAlpha, positiveRangeGamma); } /** \returns The angle values stored in a vector (alpha, beta, gamma). */ @@ -254,12 +254,10 @@ namespace Eigen * * \param m The 3x3 rotation matrix to convert * \tparam positiveRangeAlpha If true, alpha will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. - * \tparam positiveRangeBeta If true, beta will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. * \tparam positiveRangeGamma If true, gamma will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. */ template< bool PositiveRangeAlpha, - bool PositiveRangeBeta, bool PositiveRangeGamma, typename Derived> static EulerAngles FromRotation(const MatrixBase& m) @@ -268,7 +266,7 @@ namespace Eigen EulerAngles e; System::template CalcEulerAngles< - PositiveRangeAlpha, PositiveRangeBeta, PositiveRangeGamma, _Scalar>(e, m); + PositiveRangeAlpha, PositiveRangeGamma, _Scalar>(e, m); return e; } @@ -280,17 +278,15 @@ namespace Eigen * * \param rot The 3x3 rotation matrix to convert * \tparam positiveRangeAlpha If true, alpha will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. - * \tparam positiveRangeBeta If true, beta will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. * \tparam positiveRangeGamma If true, gamma will be in [0, 2*PI]. Otherwise, in [-PI, +PI]. */ template< bool PositiveRangeAlpha, - bool PositiveRangeBeta, bool PositiveRangeGamma, typename Derived> static EulerAngles FromRotation(const RotationBase& rot) { - return FromRotation(rot.toRotationMatrix()); + return FromRotation(rot.toRotationMatrix()); } /*EulerAngles& fromQuaternion(const QuaternionType& q) diff --git a/unsupported/Eigen/src/EulerAngles/EulerSystem.h b/unsupported/Eigen/src/EulerAngles/EulerSystem.h index 98f9f647d..76d0b7c57 100644 --- a/unsupported/Eigen/src/EulerAngles/EulerSystem.h +++ b/unsupported/Eigen/src/EulerAngles/EulerSystem.h @@ -112,9 +112,9 @@ namespace Eigen * * \tparam _AlphaAxis the first fixed EulerAxis * - * \tparam _AlphaAxis the second fixed EulerAxis + * \tparam _BetaAxis the second fixed EulerAxis * - * \tparam _AlphaAxis the third fixed EulerAxis + * \tparam _GammaAxis the third fixed EulerAxis */ template class EulerSystem @@ -138,14 +138,16 @@ namespace Eigen BetaAxisAbs = internal::Abs::value, /*!< the second rotation axis unsigned */ GammaAxisAbs = internal::Abs::value, /*!< the third rotation axis unsigned */ - IsAlphaOpposite = (AlphaAxis < 0) ? 1 : 0, /*!< weather alpha axis is negative */ - IsBetaOpposite = (BetaAxis < 0) ? 1 : 0, /*!< weather beta axis is negative */ - IsGammaOpposite = (GammaAxis < 0) ? 1 : 0, /*!< weather gamma axis is negative */ - - IsOdd = ((AlphaAxisAbs)%3 == (BetaAxisAbs - 1)%3) ? 0 : 1, /*!< weather the Euler system is odd */ - IsEven = IsOdd ? 0 : 1, /*!< weather the Euler system is even */ + IsAlphaOpposite = (AlphaAxis < 0) ? 1 : 0, /*!< whether alpha axis is negative */ + IsBetaOpposite = (BetaAxis < 0) ? 1 : 0, /*!< whether beta axis is negative */ + IsGammaOpposite = (GammaAxis < 0) ? 1 : 0, /*!< whether gamma axis is negative */ - IsTaitBryan = ((unsigned)AlphaAxisAbs != (unsigned)GammaAxisAbs) ? 1 : 0 /*!< weather the Euler system is tait bryan */ + // Parity is even if alpha axis X is followed by beta axis Y, or Y is followed + // by Z, or Z is followed by X; otherwise it is odd. + IsOdd = ((AlphaAxisAbs)%3 == (BetaAxisAbs - 1)%3) ? 0 : 1, /*!< whether the Euler system is odd */ + IsEven = IsOdd ? 0 : 1, /*!< whether the Euler system is even */ + + IsTaitBryan = ((unsigned)AlphaAxisAbs != (unsigned)GammaAxisAbs) ? 1 : 0 /*!< whether the Euler system is tait bryan */ }; private: @@ -180,71 +182,70 @@ namespace Eigen static void CalcEulerAngles_imp(Matrix::Scalar, 3, 1>& res, const MatrixBase& mat, internal::true_type /*isTaitBryan*/) { using std::atan2; - using std::sin; - using std::cos; + using std::sqrt; typedef typename Derived::Scalar Scalar; - typedef Matrix Vector2; - - res[0] = atan2(mat(J,K), mat(K,K)); - Scalar c2 = Vector2(mat(I,I), mat(I,J)).norm(); - if((IsOdd && res[0]Scalar(0))) { - if(res[0] > Scalar(0)) { - res[0] -= Scalar(EIGEN_PI); - } - else { - res[0] += Scalar(EIGEN_PI); - } - res[1] = atan2(-mat(I,K), -c2); + + Scalar plusMinus = IsEven? 1 : -1; + Scalar minusPlus = IsOdd? 1 : -1; + + Scalar Rsum = sqrt((mat(I,I) * mat(I,I) + mat(I,J) * mat(I,J) + mat(J,K) * mat(J,K) + mat(K,K) * mat(K,K))/2); + res[1] = atan2(plusMinus * mat(I,K), Rsum); + + // There is a singularity when cos(beta) = 0 + if(Rsum > 4 * NumTraits::epsilon()) { + res[0] = atan2(minusPlus * mat(J, K), mat(K, K)); + res[2] = atan2(minusPlus * mat(I, J), mat(I, I)); + } + else if(plusMinus * mat(I, K) > 0) { + Scalar spos = mat(J, I) + plusMinus * mat(K, J); // 2*sin(alpha + plusMinus * gamma) + Scalar cpos = mat(J, J) + minusPlus * mat(K, I); // 2*cos(alpha + plusMinus * gamma); + Scalar alphaPlusMinusGamma = atan2(spos, cpos); + res[0] = alphaPlusMinusGamma; + res[2] = 0; + } + else { + Scalar sneg = plusMinus * (mat(K, J) + minusPlus * mat(J, I)); // 2*sin(alpha + minusPlus*gamma) + Scalar cneg = mat(J, J) + plusMinus * mat(K, I); // 2*cos(alpha + minusPlus*gamma) + Scalar alphaMinusPlusBeta = atan2(sneg, cneg); + res[0] = alphaMinusPlusBeta; + res[2] = 0; } - else - res[1] = atan2(-mat(I,K), c2); - Scalar s1 = sin(res[0]); - Scalar c1 = cos(res[0]); - res[2] = atan2(s1*mat(K,I)-c1*mat(J,I), c1*mat(J,J) - s1 * mat(K,J)); } template - static void CalcEulerAngles_imp(Matrix::Scalar,3,1>& res, const MatrixBase& mat, internal::false_type /*isTaitBryan*/) + static void CalcEulerAngles_imp(Matrix::Scalar,3,1>& res, + const MatrixBase& mat, internal::false_type /*isTaitBryan*/) { using std::atan2; - using std::sin; - using std::cos; + using std::sqrt; typedef typename Derived::Scalar Scalar; - typedef Matrix Vector2; - - res[0] = atan2(mat(J,I), mat(K,I)); - if((IsOdd && res[0]Scalar(0))) - { - if(res[0] > Scalar(0)) { - res[0] -= Scalar(EIGEN_PI); - } - else { - res[0] += Scalar(EIGEN_PI); - } - Scalar s2 = Vector2(mat(J,I), mat(K,I)).norm(); - res[1] = -atan2(s2, mat(I,I)); + + Scalar plusMinus = IsEven? 1 : -1; + Scalar minusPlus = IsOdd? 1 : -1; + + Scalar Rsum = sqrt((mat(I, J) * mat(I, J) + mat(I, K) * mat(I, K) + mat(J, I) * mat(J, I) + mat(K, I) * mat(K, I)) / 2); + + res[1] = atan2(Rsum, mat(I, I)); + + if(Rsum > 4 * NumTraits::epsilon()) { + res[0] = atan2(mat(J, I), minusPlus * mat(K, I)); + res[2] = atan2(mat(I, J), plusMinus * mat(I, K)); } - else - { - Scalar s2 = Vector2(mat(J,I), mat(K,I)).norm(); - res[1] = atan2(s2, mat(I,I)); + else if( mat(I, I) > 0) { + Scalar spos = plusMinus * mat(K, J) + minusPlus * mat(J, K); // 2*sin(alpha + gamma) + Scalar cpos = mat(J, J) + mat(K, K); // 2*cos(alpha + gamma) + res[0] = atan2(spos, cpos); + res[2] = 0; + } + else { + Scalar sneg = plusMinus * mat(K, J) + plusMinus * mat(J, K); // 2*sin(alpha - gamma) + Scalar cneg = mat(J, J) - mat(K, K); // 2*cos(alpha - gamma) + res[0] = atan2(sneg, cneg); + res[1] = 0; } - // With a=(0,1,0), we have i=0; j=1; k=2, and after computing the first two angles, - // we can compute their respective rotation, and apply its inverse to M. Since the result must - // be a rotation around x, we have: - // - // c2 s1.s2 c1.s2 1 0 0 - // 0 c1 -s1 * M = 0 c3 s3 - // -s2 s1.c2 c1.c2 0 -s3 c3 - // - // Thus: m11.c1 - m21.s1 = c3 & m12.c1 - m22.s1 = s3 - - Scalar s1 = sin(res[0]); - Scalar c1 = cos(res[0]); - res[2] = atan2(c1*mat(J,K)-s1*mat(K,K), c1*mat(J,J) - s1 * mat(K,J)); } template @@ -257,14 +258,13 @@ namespace Eigen template< bool PositiveRangeAlpha, - bool PositiveRangeBeta, bool PositiveRangeGamma, typename Scalar> static void CalcEulerAngles( EulerAngles& res, const typename EulerAngles::Matrix3& mat) { - CalcEulerAngles(res, mat, PositiveRangeAlpha, PositiveRangeBeta, PositiveRangeGamma); + CalcEulerAngles(res, mat, PositiveRangeAlpha, PositiveRangeGamma); } template @@ -272,28 +272,25 @@ namespace Eigen EulerAngles& res, const typename EulerAngles::Matrix3& mat, bool PositiveRangeAlpha, - bool PositiveRangeBeta, bool PositiveRangeGamma) { CalcEulerAngles_imp( res.angles(), mat, typename internal::conditional::type()); - if (IsAlphaOpposite == IsOdd) + if (IsAlphaOpposite) res.alpha() = -res.alpha(); - if (IsBetaOpposite == IsOdd) + if (IsBetaOpposite) res.beta() = -res.beta(); - if (IsGammaOpposite == IsOdd) + if (IsGammaOpposite) res.gamma() = -res.gamma(); // Saturate results to the requested range if (PositiveRangeAlpha && (res.alpha() < 0)) res.alpha() += Scalar(2 * EIGEN_PI); - - if (PositiveRangeBeta && (res.beta() < 0)) - res.beta() += Scalar(2 * EIGEN_PI); + if (PositiveRangeGamma && (res.gamma() < 0)) res.gamma() += Scalar(2 * EIGEN_PI); diff --git a/unsupported/test/EulerAngles.cpp b/unsupported/test/EulerAngles.cpp index a8cb52864..4d0831dc2 100644 --- a/unsupported/test/EulerAngles.cpp +++ b/unsupported/test/EulerAngles.cpp @@ -15,7 +15,7 @@ using namespace Eigen; template void verify_euler_ranged(const Matrix& ea, - bool positiveRangeAlpha, bool positiveRangeBeta, bool positiveRangeGamma) + bool positiveRangeAlpha, bool positiveRangeGamma) { typedef EulerAngles EulerAnglesType; typedef Matrix Matrix3; @@ -39,10 +39,10 @@ void verify_euler_ranged(const Matrix& ea, alphaRangeEnd = Scalar(EIGEN_PI); } - if (positiveRangeBeta) + if (EulerSystem::IsTaitBryan) { - betaRangeStart = Scalar(0); - betaRangeEnd = Scalar(2 * EIGEN_PI); + betaRangeStart = -Scalar(EIGEN_PI / 2); + betaRangeEnd = Scalar(EIGEN_PI / 2); } else { @@ -61,77 +61,70 @@ void verify_euler_ranged(const Matrix& ea, gammaRangeEnd = Scalar(EIGEN_PI); } - const int i = EulerSystem::AlphaAxisAbs - 1; + /*const int i = EulerSystem::AlphaAxisAbs - 1; const int j = EulerSystem::BetaAxisAbs - 1; const int k = EulerSystem::GammaAxisAbs - 1; const int iFactor = EulerSystem::IsAlphaOpposite ? -1 : 1; const int jFactor = EulerSystem::IsBetaOpposite ? -1 : 1; - const int kFactor = EulerSystem::IsGammaOpposite ? -1 : 1; + const int kFactor = EulerSystem::IsGammaOpposite ? -1 : 1;*/ const Vector3 I = EulerAnglesType::AlphaAxisVector(); const Vector3 J = EulerAnglesType::BetaAxisVector(); const Vector3 K = EulerAnglesType::GammaAxisVector(); EulerAnglesType e(ea[0], ea[1], ea[2]); - + Matrix3 m(e); - Vector3 eabis = EulerAnglesType(m, positiveRangeAlpha, positiveRangeBeta, positiveRangeGamma).angles(); + + + Vector3 eabis = EulerAnglesType(m, positiveRangeAlpha, positiveRangeGamma).angles(); // Check that eabis in range VERIFY(alphaRangeStart <= eabis[0] && eabis[0] <= alphaRangeEnd); VERIFY(betaRangeStart <= eabis[1] && eabis[1] <= betaRangeEnd); VERIFY(gammaRangeStart <= eabis[2] && eabis[2] <= gammaRangeEnd); - - Vector3 eabis2 = m.eulerAngles(i, j, k); - - // Invert the relevant axes - eabis2[0] *= iFactor; - eabis2[1] *= jFactor; - eabis2[2] *= kFactor; - - // Saturate the angles to the correct range - if (positiveRangeAlpha && (eabis2[0] < 0)) - eabis2[0] += Scalar(2 * EIGEN_PI); - if (positiveRangeBeta && (eabis2[1] < 0)) - eabis2[1] += Scalar(2 * EIGEN_PI); - if (positiveRangeGamma && (eabis2[2] < 0)) - eabis2[2] += Scalar(2 * EIGEN_PI); - - VERIFY_IS_APPROX(eabis, eabis2);// Verify that our estimation is the same as m.eulerAngles() is - + Matrix3 mbis(AngleAxisType(eabis[0], I) * AngleAxisType(eabis[1], J) * AngleAxisType(eabis[2], K)); VERIFY_IS_APPROX(m, mbis); - - // Tests that are only relevant for no possitive range - if (!(positiveRangeAlpha || positiveRangeBeta || positiveRangeGamma)) + + // Test if ea and eabis are the same + // Need to check both singular and non-singular cases + // There are two singular cases. + // 1. When I==K and sin(ea(1)) == 0 + // 2. When I!=K and cos(ea(1)) == 0 + + // Tests that are only relevant for no positive range + /*if (!(positiveRangeAlpha || positiveRangeGamma)) { - /* If I==K, and ea[1]==0, then there no unique solution. */ - /* The remark apply in the case where I!=K, and |ea[1]| is close to pi/2. */ + // If I==K, and ea[1]==0, then there no unique solution. + // The remark apply in the case where I!=K, and |ea[1]| is close to pi/2. if( (i!=k || ea[1]!=0) && (i==k || !internal::isApprox(abs(ea[1]),Scalar(EIGEN_PI/2),test_precision())) ) VERIFY((ea-eabis).norm() <= test_precision()); // approx_or_less_than does not work for 0 VERIFY(0 < eabis[0] || test_isMuchSmallerThan(eabis[0], Scalar(1))); - } + }*/ // Quaternions QuaternionType q(e); - eabis = EulerAnglesType(q, positiveRangeAlpha, positiveRangeBeta, positiveRangeGamma).angles(); - VERIFY_IS_APPROX(eabis, eabis2);// Verify that the euler angles are still the same + eabis = EulerAnglesType(q, positiveRangeAlpha, positiveRangeGamma).angles(); + QuaternionType qbis(AngleAxisType(eabis[0], I) * AngleAxisType(eabis[1], J) * AngleAxisType(eabis[2], K)); + VERIFY_IS_APPROX(std::abs(q.dot(qbis)), static_cast(1)); + //VERIFY_IS_APPROX(eabis, eabis2);// Verify that the euler angles are still the same } template void verify_euler(const Matrix& ea) { - verify_euler_ranged(ea, false, false, false); - verify_euler_ranged(ea, false, false, true); - verify_euler_ranged(ea, false, true, false); - verify_euler_ranged(ea, false, true, true); - verify_euler_ranged(ea, true, false, false); - verify_euler_ranged(ea, true, false, true); - verify_euler_ranged(ea, true, true, false); - verify_euler_ranged(ea, true, true, true); + verify_euler_ranged(ea, false, false); + verify_euler_ranged(ea, false, true); + verify_euler_ranged(ea, false, false); + verify_euler_ranged(ea, false, true); + verify_euler_ranged(ea, true, false); + verify_euler_ranged(ea, true, true); + verify_euler_ranged(ea, true, false); + verify_euler_ranged(ea, true, true); } template void check_all_var(const Matrix& ea) From 58f5d7d058e21bec85d902504efe988d17aa28cf Mon Sep 17 00:00:00 2001 From: Tal Hadad Date: Sun, 16 Oct 2016 14:39:26 +0300 Subject: [PATCH 02/22] Fix calc bug, docs and better testing. Test code changes: * better coded * rand and manual numbers * singularity checking --- .../Eigen/src/EulerAngles/EulerAngles.h | 28 +-- .../Eigen/src/EulerAngles/EulerSystem.h | 36 ++-- unsupported/test/EulerAngles.cpp | 177 +++++++++++++----- 3 files changed, 164 insertions(+), 77 deletions(-) diff --git a/unsupported/Eigen/src/EulerAngles/EulerAngles.h b/unsupported/Eigen/src/EulerAngles/EulerAngles.h index da86cc13b..8a723d9ee 100644 --- a/unsupported/Eigen/src/EulerAngles/EulerAngles.h +++ b/unsupported/Eigen/src/EulerAngles/EulerAngles.h @@ -36,7 +36,7 @@ namespace Eigen * ### Rotation representation and conversions ### * * It has been proved(see Wikipedia link below) that every rotation can be represented - * by Euler angles, but there is no singular representation (e.g. unlike rotation matrices). + * by Euler angles, but there is no single representation (e.g. unlike rotation matrices). * Therefore, you can convert from Eigen rotation and to them * (including rotation matrices, which is not called "rotations" by Eigen design). * @@ -55,10 +55,12 @@ namespace Eigen * Additionally, some axes related computation is done in compile time. * * #### Euler angles ranges in conversions #### - * Rotations representation as EulerAngles are not singular (unlike matrices), and even have infinite EulerAngles representations.
+ * Rotations representation as EulerAngles are not single (unlike matrices), + * and even have infinite EulerAngles representations.
* For example, add or subtract 2*PI from either angle of EulerAngles * and you'll get the same rotation. - * This is the reason for infinite representation, but it's not the only reason for non-singularity. + * This is the general reason for infinite representation, + * but it's not the only general reason for not having a single representation. * * When converting rotation to EulerAngles, this class convert it to specific ranges * When converting some rotation to EulerAngles, the rules for ranges are as follow: @@ -66,10 +68,10 @@ namespace Eigen * (even when it represented as RotationBase explicitly), angles ranges are __undefined__. * - otherwise, Alpha and Gamma angles will be in the range [-PI, PI].
* As for Beta angle: - * - If the system is Tait-Bryan, the beta angle will be in the range [-PI, PI]. + * - If the system is Tait-Bryan, the beta angle will be in the range [-PI/2, PI/2]. * - otherwise: - * - If the beta axis is positive, the beta angle will be in the range [0, 2*PI] - * - If the beta axis is negative, the beta angle will be in the range [-2*PI, 0] + * - If the beta axis is positive, the beta angle will be in the range [0, PI] + * - If the beta axis is negative, the beta angle will be in the range [-PI, 0] * * \sa EulerAngles(const MatrixBase&) * \sa EulerAngles(const RotationBase&) @@ -95,7 +97,7 @@ namespace Eigen * * More information about Euler angles: https://en.wikipedia.org/wiki/Euler_angles * - * \tparam _Scalar the scalar type, i.e., the type of the angles. + * \tparam _Scalar the scalar type, i.e. the type of the angles. * * \tparam _System the EulerSystem to use, which represents the axes of rotation. */ @@ -146,10 +148,10 @@ namespace Eigen * * \note Alpha and Gamma angles will be in the range [-PI, PI].
* As for Beta angle: - * - If the system is Tait-Bryan, the beta angle will be in the range [-PI, PI]. + * - If the system is Tait-Bryan, the beta angle will be in the range [-PI/2, PI/2]. * - otherwise: - * - If the beta axis is positive, the beta angle will be in the range [0, 2*PI] - * - If the beta axis is negative, the beta angle will be in the range [-2*PI, 0] + * - If the beta axis is positive, the beta angle will be in the range [0, PI] + * - If the beta axis is negative, the beta angle will be in the range [-PI, 0] */ template EulerAngles(const MatrixBase& m) { System::CalcEulerAngles(*this, m); } @@ -160,10 +162,10 @@ namespace Eigen * angles ranges are __undefined__. * Otherwise, Alpha and Gamma angles will be in the range [-PI, PI].
* As for Beta angle: - * - If the system is Tait-Bryan, the beta angle will be in the range [-PI, PI]. + * - If the system is Tait-Bryan, the beta angle will be in the range [-PI/2, PI/2]. * - otherwise: - * - If the beta axis is positive, the beta angle will be in the range [0, 2*PI] - * - If the beta axis is negative, the beta angle will be in the range [-2*PI, 0] + * - If the beta axis is positive, the beta angle will be in the range [0, PI] + * - If the beta axis is negative, the beta angle will be in the range [-PI, 0] */ template EulerAngles(const RotationBase& rot) { System::CalcEulerAngles(*this, rot.toRotationMatrix()); } diff --git a/unsupported/Eigen/src/EulerAngles/EulerSystem.h b/unsupported/Eigen/src/EulerAngles/EulerSystem.h index aa96461f9..0790e612f 100644 --- a/unsupported/Eigen/src/EulerAngles/EulerSystem.h +++ b/unsupported/Eigen/src/EulerAngles/EulerSystem.h @@ -18,7 +18,7 @@ namespace Eigen namespace internal { - // TODO: Check if already exists on the rest API + // TODO: Add this trait to the Eigen internal API? template 0)> struct Abs { @@ -186,25 +186,25 @@ namespace Eigen typedef typename Derived::Scalar Scalar; - Scalar plusMinus = IsEven? 1 : -1; - Scalar minusPlus = IsOdd? 1 : -1; + const Scalar plusMinus = IsEven? 1 : -1; + const Scalar minusPlus = IsOdd? 1 : -1; - Scalar Rsum = sqrt((mat(I,I) * mat(I,I) + mat(I,J) * mat(I,J) + mat(J,K) * mat(J,K) + mat(K,K) * mat(K,K))/2); + const Scalar Rsum = sqrt((mat(I,I) * mat(I,I) + mat(I,J) * mat(I,J) + mat(J,K) * mat(J,K) + mat(K,K) * mat(K,K))/2); res[1] = atan2(plusMinus * mat(I,K), Rsum); - // There is a singularity when cos(beta) = 0 - if(Rsum > 4 * NumTraits::epsilon()) { + // There is a singularity when cos(beta) == 0 + if(Rsum > 4 * NumTraits::epsilon()) {// cos(beta) != 0 res[0] = atan2(minusPlus * mat(J, K), mat(K, K)); res[2] = atan2(minusPlus * mat(I, J), mat(I, I)); } - else if(plusMinus * mat(I, K) > 0) { - Scalar spos = mat(J, I) + plusMinus * mat(K, J); // 2*sin(alpha + plusMinus * gamma) - Scalar cpos = mat(J, J) + minusPlus * mat(K, I); // 2*cos(alpha + plusMinus * gamma); + else if(plusMinus * mat(I, K) > 0) {// cos(beta) == 0 and sin(beta) == 1 + Scalar spos = mat(J, I) + plusMinus * mat(K, J); // 2*sin(alpha + plusMinus * gamma + Scalar cpos = mat(J, J) + minusPlus * mat(K, I); // 2*cos(alpha + plusMinus * gamma) Scalar alphaPlusMinusGamma = atan2(spos, cpos); res[0] = alphaPlusMinusGamma; res[2] = 0; } - else { + else {// cos(beta) == 0 and sin(beta) == -1 Scalar sneg = plusMinus * (mat(K, J) + minusPlus * mat(J, I)); // 2*sin(alpha + minusPlus*gamma) Scalar cneg = mat(J, J) + plusMinus * mat(K, I); // 2*cos(alpha + minusPlus*gamma) Scalar alphaMinusPlusBeta = atan2(sneg, cneg); @@ -222,30 +222,30 @@ namespace Eigen typedef typename Derived::Scalar Scalar; - Scalar plusMinus = IsEven? 1 : -1; - Scalar minusPlus = IsOdd? 1 : -1; + const Scalar plusMinus = IsEven? 1 : -1; + const Scalar minusPlus = IsOdd? 1 : -1; - Scalar Rsum = sqrt((mat(I, J) * mat(I, J) + mat(I, K) * mat(I, K) + mat(J, I) * mat(J, I) + mat(K, I) * mat(K, I)) / 2); + const Scalar Rsum = sqrt((mat(I, J) * mat(I, J) + mat(I, K) * mat(I, K) + mat(J, I) * mat(J, I) + mat(K, I) * mat(K, I)) / 2); res[1] = atan2(Rsum, mat(I, I)); - if(Rsum > 4 * NumTraits::epsilon()) { + // There is a singularity when sin(beta) == 0 + if(Rsum > 4 * NumTraits::epsilon()) {// sin(beta) != 0 res[0] = atan2(mat(J, I), minusPlus * mat(K, I)); res[2] = atan2(mat(I, J), plusMinus * mat(I, K)); } - else if( mat(I, I) > 0) { + else if(mat(I, I) > 0) {// sin(beta) == 0 and cos(beta) == 1 Scalar spos = plusMinus * mat(K, J) + minusPlus * mat(J, K); // 2*sin(alpha + gamma) Scalar cpos = mat(J, J) + mat(K, K); // 2*cos(alpha + gamma) res[0] = atan2(spos, cpos); res[2] = 0; } - else { + else {// sin(beta) == 0 and cos(beta) == -1 Scalar sneg = plusMinus * mat(K, J) + plusMinus * mat(J, K); // 2*sin(alpha - gamma) Scalar cneg = mat(J, J) - mat(K, K); // 2*cos(alpha - gamma) res[0] = atan2(sneg, cneg); - res[1] = 0; + res[2] = 0; } - } template diff --git a/unsupported/test/EulerAngles.cpp b/unsupported/test/EulerAngles.cpp index 8b4706686..43291b31d 100644 --- a/unsupported/test/EulerAngles.cpp +++ b/unsupported/test/EulerAngles.cpp @@ -15,13 +15,17 @@ using namespace Eigen; // Verify that x is in the approxed range [a, b] #define VERIFY_APPROXED_RANGE(a, x, b) \ - do { \ - VERIFY_IS_APPROX_OR_LESS_THAN(a, x); \ - VERIFY_IS_APPROX_OR_LESS_THAN(x, b); \ - } while(0) + do { \ + VERIFY_IS_APPROX_OR_LESS_THAN(a, x); \ + VERIFY_IS_APPROX_OR_LESS_THAN(x, b); \ + } while(0) -template -void verify_euler(const Matrix& ea) +const char X = EULER_X; +const char Y = EULER_Y; +const char Z = EULER_Z; + +template +void verify_euler(const EulerAngles& e) { typedef EulerAngles EulerAnglesType; typedef Matrix Matrix3; @@ -41,17 +45,24 @@ void verify_euler(const Matrix& ea) } else { - betaRangeStart = -PI; - betaRangeEnd = PI; + if (!EulerSystem::IsBetaOpposite) + { + betaRangeStart = 0; + betaRangeEnd = PI; + } + else + { + betaRangeStart = -PI; + betaRangeEnd = 0; + } } const Vector3 I = EulerAnglesType::AlphaAxisVector(); const Vector3 J = EulerAnglesType::BetaAxisVector(); const Vector3 K = EulerAnglesType::GammaAxisVector(); - - EulerAnglesType e(ea[0], ea[1], ea[2]); - Matrix3 m(e); + const Matrix3 m(e); + VERIFY_IS_APPROX(Scalar(m.determinant()), ONE); Vector3 eabis = static_cast(m).angles(); @@ -60,8 +71,16 @@ void verify_euler(const Matrix& ea) VERIFY_APPROXED_RANGE(betaRangeStart, eabis[1], betaRangeEnd); VERIFY_APPROXED_RANGE(-PI, eabis[2], PI); - Matrix3 mbis(AngleAxisType(eabis[0], I) * AngleAxisType(eabis[1], J) * AngleAxisType(eabis[2], K)); - VERIFY_IS_APPROX(m, mbis); + const Matrix3 mbis(AngleAxisType(eabis[0], I) * AngleAxisType(eabis[1], J) * AngleAxisType(eabis[2], K)); + VERIFY_IS_APPROX(Scalar(mbis.determinant()), ONE); + /*std::cout << "===================\n" << + "e: " << e << std::endl << + "eabis: " << eabis.transpose() << std::endl << + "m: " << m << std::endl << + "mbis: " << mbis << std::endl << + "X: " << (m * Vector3::UnitX()).transpose() << std::endl << + "X: " << (mbis * Vector3::UnitX()).transpose() << std::endl;*/ + VERIFY_IS_APPROX(m, mbis); // Test if ea and eabis are the same // Need to check both singular and non-singular cases @@ -69,47 +88,107 @@ void verify_euler(const Matrix& ea) // 1. When I==K and sin(ea(1)) == 0 // 2. When I!=K and cos(ea(1)) == 0 - // Tests that are only relevant for no positive range - /*if (!(positiveRangeAlpha || positiveRangeGamma)) - { - // If I==K, and ea[1]==0, then there no unique solution. - // The remark apply in the case where I!=K, and |ea[1]| is close to pi/2. - if( (i!=k || ea[1]!=0) && (i==k || !internal::isApprox(abs(ea[1]),Scalar(EIGEN_PI/2),test_precision())) ) - VERIFY((ea-eabis).norm() <= test_precision()); - - // approx_or_less_than does not work for 0 - VERIFY(0 < eabis[0] || VERIFY_IS_MUCH_SMALLER_THAN(eabis[0], Scalar(1))); - }*/ + // TODO: Make this test work well, and use range saturation function. + /*// If I==K, and ea[1]==0, then there no unique solution. + // The remark apply in the case where I!=K, and |ea[1]| is close to +-pi/2. + if( (i!=k || ea[1]!=0) && (i==k || !internal::isApprox(abs(ea[1]),Scalar(EIGEN_PI/2),test_precision())) ) + VERIFY_IS_APPROX(ea, eabis);*/ // Quaternions - QuaternionType q(e); + const QuaternionType q(e); eabis = static_cast(q).angles(); - QuaternionType qbis(AngleAxisType(eabis[0], I) * AngleAxisType(eabis[1], J) * AngleAxisType(eabis[2], K)); + const QuaternionType qbis(AngleAxisType(eabis[0], I) * AngleAxisType(eabis[1], J) * AngleAxisType(eabis[2], K)); VERIFY_IS_APPROX(std::abs(q.dot(qbis)), ONE); //VERIFY_IS_APPROX(eabis, eabis2);// Verify that the euler angles are still the same } -template void check_all_var(const Matrix& ea) +template +void verify_euler_vec(const Matrix& ea) { - verify_euler(ea); - verify_euler(ea); - verify_euler(ea); - verify_euler(ea); - - verify_euler(ea); - verify_euler(ea); - verify_euler(ea); - verify_euler(ea); - - verify_euler(ea); - verify_euler(ea); - verify_euler(ea); - verify_euler(ea); - - // TODO: Test negative axes as well! (only test if the angles get negative when needed) + verify_euler(EulerAngles >(ea[0], ea[1], ea[2])); } -template void eulerangles() +template +void verify_euler_all_neg(const Matrix& ea) +{ + verify_euler_vec<+A,+B,+C>(ea); + verify_euler_vec<+A,+B,-C>(ea); + verify_euler_vec<+A,-B,+C>(ea); + verify_euler_vec<+A,-B,-C>(ea); + + verify_euler_vec<-A,+B,+C>(ea); + verify_euler_vec<-A,+B,-C>(ea); + verify_euler_vec<-A,-B,+C>(ea); + verify_euler_vec<-A,-B,-C>(ea); +} + +template void check_all_var(const Matrix& ea) +{ + verify_euler_all_neg(ea); + verify_euler_all_neg(ea); + verify_euler_all_neg(ea); + verify_euler_all_neg(ea); + + verify_euler_all_neg(ea); + verify_euler_all_neg(ea); + verify_euler_all_neg(ea); + verify_euler_all_neg(ea); + + verify_euler_all_neg(ea); + verify_euler_all_neg(ea); + verify_euler_all_neg(ea); + verify_euler_all_neg(ea); +} + +template void check_singular_cases(const Scalar& singularBeta) +{ + typedef Matrix Vector3; + const Scalar epsilon = std::numeric_limits::epsilon(); + const Scalar PI = Scalar(EIGEN_PI); + + check_all_var(Vector3(PI/4, singularBeta, PI/3)); + check_all_var(Vector3(PI/4, singularBeta - epsilon, PI/3)); + check_all_var(Vector3(PI/4, singularBeta - Scalar(1.5)*epsilon, PI/3)); + check_all_var(Vector3(PI/4, singularBeta - 2*epsilon, PI/3)); + check_all_var(Vector3(PI*Scalar(0.8), singularBeta - epsilon, Scalar(0.9)*PI)); + check_all_var(Vector3(PI*Scalar(-0.9), singularBeta + epsilon, PI*Scalar(0.3))); + check_all_var(Vector3(PI*Scalar(-0.6), singularBeta + Scalar(1.5)*epsilon, PI*Scalar(0.3))); + check_all_var(Vector3(PI*Scalar(-0.5), singularBeta + 2*epsilon, PI*Scalar(0.4))); + check_all_var(Vector3(PI*Scalar(0.9), singularBeta + epsilon, Scalar(0.8)*PI)); +} + +template void eulerangles_manual() +{ + typedef Matrix Vector3; + const Vector3 Zero = Vector3::Zero(); + const Scalar PI = Scalar(EIGEN_PI); + + check_all_var(Zero); + + // singular cases + check_singular_cases(PI/2); + check_singular_cases(-PI/2); + + check_singular_cases(Scalar(0)); + check_singular_cases(Scalar(-0)); + + check_singular_cases(PI); + check_singular_cases(-PI); + + // non-singular cases + VectorXd alpha = VectorXd::LinSpaced(Eigen::Sequential, 20, Scalar(-0.99) * PI, PI); + VectorXd beta = VectorXd::LinSpaced(Eigen::Sequential, 20, Scalar(-0.49) * PI, Scalar(0.49) * PI); + VectorXd gamma = VectorXd::LinSpaced(Eigen::Sequential, 20, Scalar(-0.99) * PI, PI); + for (int i = 0; i < alpha.size(); ++i) { + for (int j = 0; j < beta.size(); ++j) { + for (int k = 0; k < gamma.size(); ++k) { + check_all_var(Vector3d(alpha(i), beta(j), gamma(k))); + } + } + } +} + +template void eulerangles_rand() { typedef Matrix Matrix3; typedef Matrix Vector3; @@ -158,8 +237,14 @@ template void eulerangles() void test_EulerAngles() { + CALL_SUBTEST_1( eulerangles_manual() ); + CALL_SUBTEST_2( eulerangles_manual() ); + for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_1( eulerangles() ); - CALL_SUBTEST_2( eulerangles() ); + CALL_SUBTEST_3( eulerangles_rand() ); + CALL_SUBTEST_4( eulerangles_rand() ); } + + // TODO: Add tests for auto diff + // TODO: Add tests for complex numbers } From 7402cfd4cc3bc129b0fc906eb51347882307cbe0 Mon Sep 17 00:00:00 2001 From: Tal Hadad Date: Mon, 17 Oct 2016 20:42:08 +0300 Subject: [PATCH 03/22] Add safty for near pole cases and test them better. --- unsupported/test/EulerAngles.cpp | 37 +++++++++++++++++++++----------- 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/unsupported/test/EulerAngles.cpp b/unsupported/test/EulerAngles.cpp index 43291b31d..3f4523ccc 100644 --- a/unsupported/test/EulerAngles.cpp +++ b/unsupported/test/EulerAngles.cpp @@ -37,6 +37,14 @@ void verify_euler(const EulerAngles& e) const Scalar HALF_PI = Scalar(EIGEN_PI / 2); const Scalar PI = Scalar(EIGEN_PI); + // It's very important calc the acceptable precision depending on the distance from the pole. + const Scalar longitudeRadius = std::abs( + EulerSystem::IsTaitBryan ? + std::cos(e.beta()) : + std::sin(e.beta()) + ); + const Scalar precision = test_precision() / longitudeRadius; + Scalar betaRangeStart, betaRangeEnd; if (EulerSystem::IsTaitBryan) { @@ -80,7 +88,7 @@ void verify_euler(const EulerAngles& e) "mbis: " << mbis << std::endl << "X: " << (m * Vector3::UnitX()).transpose() << std::endl << "X: " << (mbis * Vector3::UnitX()).transpose() << std::endl;*/ - VERIFY_IS_APPROX(m, mbis); + VERIFY(m.isApprox(mbis, precision)); // Test if ea and eabis are the same // Need to check both singular and non-singular cases @@ -98,7 +106,7 @@ void verify_euler(const EulerAngles& e) const QuaternionType q(e); eabis = static_cast(q).angles(); const QuaternionType qbis(AngleAxisType(eabis[0], I) * AngleAxisType(eabis[1], J) * AngleAxisType(eabis[2], K)); - VERIFY_IS_APPROX(std::abs(q.dot(qbis)), ONE); + VERIFY(internal::isApprox(std::abs(q.dot(qbis)), ONE, precision)); //VERIFY_IS_APPROX(eabis, eabis2);// Verify that the euler angles are still the same } @@ -143,18 +151,23 @@ template void check_all_var(const Matrix& ea) template void check_singular_cases(const Scalar& singularBeta) { typedef Matrix Vector3; - const Scalar epsilon = std::numeric_limits::epsilon(); const Scalar PI = Scalar(EIGEN_PI); - check_all_var(Vector3(PI/4, singularBeta, PI/3)); - check_all_var(Vector3(PI/4, singularBeta - epsilon, PI/3)); - check_all_var(Vector3(PI/4, singularBeta - Scalar(1.5)*epsilon, PI/3)); - check_all_var(Vector3(PI/4, singularBeta - 2*epsilon, PI/3)); - check_all_var(Vector3(PI*Scalar(0.8), singularBeta - epsilon, Scalar(0.9)*PI)); - check_all_var(Vector3(PI*Scalar(-0.9), singularBeta + epsilon, PI*Scalar(0.3))); - check_all_var(Vector3(PI*Scalar(-0.6), singularBeta + Scalar(1.5)*epsilon, PI*Scalar(0.3))); - check_all_var(Vector3(PI*Scalar(-0.5), singularBeta + 2*epsilon, PI*Scalar(0.4))); - check_all_var(Vector3(PI*Scalar(0.9), singularBeta + epsilon, Scalar(0.8)*PI)); + for (Scalar epsilon = std::numeric_limits::epsilon(); epsilon < 1; epsilon *= Scalar(1.2)) + { + check_all_var(Vector3(PI/4, singularBeta, PI/3)); + check_all_var(Vector3(PI/4, singularBeta - epsilon, PI/3)); + check_all_var(Vector3(PI/4, singularBeta - Scalar(1.5)*epsilon, PI/3)); + check_all_var(Vector3(PI/4, singularBeta - 2*epsilon, PI/3)); + check_all_var(Vector3(PI*Scalar(0.8), singularBeta - epsilon, Scalar(0.9)*PI)); + check_all_var(Vector3(PI*Scalar(-0.9), singularBeta + epsilon, PI*Scalar(0.3))); + check_all_var(Vector3(PI*Scalar(-0.6), singularBeta + Scalar(1.5)*epsilon, PI*Scalar(0.3))); + check_all_var(Vector3(PI*Scalar(-0.5), singularBeta + 2*epsilon, PI*Scalar(0.4))); + check_all_var(Vector3(PI*Scalar(0.9), singularBeta + epsilon, Scalar(0.8)*PI)); + } + + // This one for sanity, it had a problem with near pole cases in float scalar. + check_all_var(Vector3(PI*Scalar(0.8), singularBeta - Scalar(1E-6), Scalar(0.9)*PI)); } template void eulerangles_manual() From 6f4f12d1eda685f03a2e4b54b479813d134248cf Mon Sep 17 00:00:00 2001 From: Tal Hadad Date: Mon, 17 Oct 2016 22:23:47 +0300 Subject: [PATCH 04/22] Add isApprox() and cast() functions. test cases included --- .../Eigen/src/EulerAngles/EulerAngles.h | 20 +++++++++++++++- unsupported/test/EulerAngles.cpp | 24 ++++++++++++++++++- 2 files changed, 42 insertions(+), 2 deletions(-) diff --git a/unsupported/Eigen/src/EulerAngles/EulerAngles.h b/unsupported/Eigen/src/EulerAngles/EulerAngles.h index 8a723d9ee..6594e4d13 100644 --- a/unsupported/Eigen/src/EulerAngles/EulerAngles.h +++ b/unsupported/Eigen/src/EulerAngles/EulerAngles.h @@ -105,8 +105,11 @@ namespace Eigen class EulerAngles : public RotationBase, 3> { public: + typedef RotationBase, 3> Base; + /** the scalar type of the angles */ typedef _Scalar Scalar; + typedef typename NumTraits::Real RealScalar; /** the EulerSystem to use, which represents the axes of rotation. */ typedef _System System; @@ -248,7 +251,13 @@ namespace Eigen return *this; } - // TODO: Support isApprox function + /** \returns \c true if \c *this is approximately equal to \a other, within the precision + * determined by \a prec. + * + * \sa MatrixBase::isApprox() */ + bool isApprox(const EulerAngles& other, + const RealScalar& prec = NumTraits::dummy_precision()) const + { return angles().isApprox(other.angles(), prec); } /** \returns an equivalent 3x3 rotation matrix. */ Matrix3 toRotationMatrix() const @@ -271,6 +280,15 @@ namespace Eigen s << eulerAngles.angles().transpose(); return s; } + + /** \returns \c *this with scalar type casted to \a NewScalarType */ + template + EulerAngles cast() const + { + EulerAngles e; + e.angles() = angles().cast(); + return e; + } }; #define EIGEN_EULER_ANGLES_SINGLE_TYPEDEF(AXES, SCALAR_TYPE, SCALAR_POSTFIX) \ diff --git a/unsupported/test/EulerAngles.cpp b/unsupported/test/EulerAngles.cpp index 3f4523ccc..149cf7f94 100644 --- a/unsupported/test/EulerAngles.cpp +++ b/unsupported/test/EulerAngles.cpp @@ -13,6 +13,13 @@ using namespace Eigen; +// Unfortunately, we need to specialize it in order to work. (We could add it in main.h test framework) +template +bool verifyIsApprox(const Eigen::EulerAngles& a, const Eigen::EulerAngles& b) +{ + return verifyIsApprox(a.angles(), b.angles()); +} + // Verify that x is in the approxed range [a, b] #define VERIFY_APPROXED_RANGE(a, x, b) \ do { \ @@ -24,7 +31,7 @@ const char X = EULER_X; const char Y = EULER_Y; const char Z = EULER_Z; -template +template void verify_euler(const EulerAngles& e) { typedef EulerAngles EulerAnglesType; @@ -68,6 +75,11 @@ void verify_euler(const EulerAngles& e) const Vector3 I = EulerAnglesType::AlphaAxisVector(); const Vector3 J = EulerAnglesType::BetaAxisVector(); const Vector3 K = EulerAnglesType::GammaAxisVector(); + + // Is approx checks + VERIFY(e.isApprox(e)); + VERIFY_IS_APPROX(e, e); + VERIFY_IS_NOT_APPROX(e, EulerAnglesType(e.alpha() + ONE, e.beta() + ONE, e.gamma() + ONE)); const Matrix3 m(e); VERIFY_IS_APPROX(Scalar(m.determinant()), ONE); @@ -108,6 +120,11 @@ void verify_euler(const EulerAngles& e) const QuaternionType qbis(AngleAxisType(eabis[0], I) * AngleAxisType(eabis[1], J) * AngleAxisType(eabis[2], K)); VERIFY(internal::isApprox(std::abs(q.dot(qbis)), ONE, precision)); //VERIFY_IS_APPROX(eabis, eabis2);// Verify that the euler angles are still the same + + // A suggestion for simple product test when will be supported. + /*EulerAnglesType e2(PI/2, PI/2, PI/2); + Matrix3 m2(e2); + VERIFY_IS_APPROX(e*e2, m*m2);*/ } template @@ -250,6 +267,11 @@ template void eulerangles_rand() void test_EulerAngles() { + // Simple cast test + EulerAnglesXYZd onesEd(1, 1, 1); + EulerAnglesXYZf onesEf = onesEd.cast(); + VERIFY_IS_APPROX(onesEd, onesEf.cast()); + CALL_SUBTEST_1( eulerangles_manual() ); CALL_SUBTEST_2( eulerangles_manual() ); From 15eca2432a071815884af6ccc2ceb66be91e9db7 Mon Sep 17 00:00:00 2001 From: Tal Hadad Date: Tue, 18 Oct 2016 23:24:57 +0300 Subject: [PATCH 05/22] Euler tests: Tighter precision when no roll exists and clean code. --- unsupported/test/EulerAngles.cpp | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/unsupported/test/EulerAngles.cpp b/unsupported/test/EulerAngles.cpp index 149cf7f94..79ee72847 100644 --- a/unsupported/test/EulerAngles.cpp +++ b/unsupported/test/EulerAngles.cpp @@ -50,7 +50,7 @@ void verify_euler(const EulerAngles& e) std::cos(e.beta()) : std::sin(e.beta()) ); - const Scalar precision = test_precision() / longitudeRadius; + Scalar precision = test_precision() / longitudeRadius; Scalar betaRangeStart, betaRangeEnd; if (EulerSystem::IsTaitBryan) @@ -84,15 +84,22 @@ void verify_euler(const EulerAngles& e) const Matrix3 m(e); VERIFY_IS_APPROX(Scalar(m.determinant()), ONE); - Vector3 eabis = static_cast(m).angles(); + EulerAnglesType ebis(m); + + // When no roll(acting like polar representation), we have the best precision. + // One of those cases is when the Euler angles are on the pole, and because it's singular case, + // the computation returns no roll. + if (ebis.beta() == 0) + precision = test_precision(); // Check that eabis in range - VERIFY_APPROXED_RANGE(-PI, eabis[0], PI); - VERIFY_APPROXED_RANGE(betaRangeStart, eabis[1], betaRangeEnd); - VERIFY_APPROXED_RANGE(-PI, eabis[2], PI); + VERIFY_APPROXED_RANGE(-PI, ebis.alpha(), PI); + VERIFY_APPROXED_RANGE(betaRangeStart, ebis.beta(), betaRangeEnd); + VERIFY_APPROXED_RANGE(-PI, ebis.gamma(), PI); - const Matrix3 mbis(AngleAxisType(eabis[0], I) * AngleAxisType(eabis[1], J) * AngleAxisType(eabis[2], K)); + const Matrix3 mbis(AngleAxisType(ebis.alpha(), I) * AngleAxisType(ebis.beta(), J) * AngleAxisType(ebis.gamma(), K)); VERIFY_IS_APPROX(Scalar(mbis.determinant()), ONE); + VERIFY_IS_APPROX(mbis, ebis.toRotationMatrix()); /*std::cout << "===================\n" << "e: " << e << std::endl << "eabis: " << eabis.transpose() << std::endl << @@ -116,8 +123,8 @@ void verify_euler(const EulerAngles& e) // Quaternions const QuaternionType q(e); - eabis = static_cast(q).angles(); - const QuaternionType qbis(AngleAxisType(eabis[0], I) * AngleAxisType(eabis[1], J) * AngleAxisType(eabis[2], K)); + ebis = q; + const QuaternionType qbis(ebis); VERIFY(internal::isApprox(std::abs(q.dot(qbis)), ONE, precision)); //VERIFY_IS_APPROX(eabis, eabis2);// Verify that the euler angles are still the same @@ -170,7 +177,7 @@ template void check_singular_cases(const Scalar& singularBeta) typedef Matrix Vector3; const Scalar PI = Scalar(EIGEN_PI); - for (Scalar epsilon = std::numeric_limits::epsilon(); epsilon < 1; epsilon *= Scalar(1.2)) + for (Scalar epsilon = NumTraits::epsilon(); epsilon < 1; epsilon *= Scalar(1.2)) { check_all_var(Vector3(PI/4, singularBeta, PI/3)); check_all_var(Vector3(PI/4, singularBeta - epsilon, PI/3)); From 76b2a3e6e70e4760755d7fc5c90e807718db92e4 Mon Sep 17 00:00:00 2001 From: Tal Hadad Date: Fri, 18 Nov 2016 15:01:06 +0200 Subject: [PATCH 06/22] Allow to construct EulerAngles from 3D vector directly. Using assignment template struct to distinguish between 3D vector and 3D rotation matrix. --- .../Eigen/src/EulerAngles/EulerAngles.h | 63 +++++++++++++------ .../Eigen/src/EulerAngles/EulerSystem.h | 12 ++++ 2 files changed, 57 insertions(+), 18 deletions(-) diff --git a/unsupported/Eigen/src/EulerAngles/EulerAngles.h b/unsupported/Eigen/src/EulerAngles/EulerAngles.h index 6594e4d13..2a12c8da3 100644 --- a/unsupported/Eigen/src/EulerAngles/EulerAngles.h +++ b/unsupported/Eigen/src/EulerAngles/EulerAngles.h @@ -12,11 +12,6 @@ namespace Eigen { - /*template - struct ei_eulerangles_assign_impl;*/ - /** \class EulerAngles * * \ingroup EulerAngles_Module @@ -66,7 +61,7 @@ namespace Eigen * When converting some rotation to EulerAngles, the rules for ranges are as follow: * - If the rotation we converting from is an EulerAngles * (even when it represented as RotationBase explicitly), angles ranges are __undefined__. - * - otherwise, Alpha and Gamma angles will be in the range [-PI, PI].
+ * - otherwise, alpha and gamma angles will be in the range [-PI, PI].
* As for Beta angle: * - If the system is Tait-Bryan, the beta angle will be in the range [-PI/2, PI/2]. * - otherwise: @@ -143,27 +138,34 @@ namespace Eigen public: /** Default constructor without initialization. */ EulerAngles() {} - /** Constructs and initialize Euler angles(\p alpha, \p beta, \p gamma). */ + /** Constructs and initialize an EulerAngles (\p alpha, \p beta, \p gamma). */ EulerAngles(const Scalar& alpha, const Scalar& beta, const Scalar& gamma) : m_angles(alpha, beta, gamma) {} - /** Constructs and initialize Euler angles from a 3x3 rotation matrix \p m. + // TODO: Test this constructor + /** Constructs and initialize an EulerAngles from the array data {alpha, beta, gamma} */ + explicit EulerAngles(const Scalar* data) : m_angles(data) {} + + /** Constructs and initializes an EulerAngles from either: + * - a 3x3 rotation matrix expression(i.e. pure orthogonal matrix with determinant of +1), + * - a 3D vector expression representing Euler angles. * - * \note Alpha and Gamma angles will be in the range [-PI, PI].
+ * \note If \p other is a 3x3 rotation matrix, the angles range rules will be as follow:
+ * Alpha and gamma angles will be in the range [-PI, PI].
* As for Beta angle: * - If the system is Tait-Bryan, the beta angle will be in the range [-PI/2, PI/2]. * - otherwise: * - If the beta axis is positive, the beta angle will be in the range [0, PI] * - If the beta axis is negative, the beta angle will be in the range [-PI, 0] - */ + */ template - EulerAngles(const MatrixBase& m) { System::CalcEulerAngles(*this, m); } + explicit EulerAngles(const MatrixBase& other) { *this = other; } /** Constructs and initialize Euler angles from a rotation \p rot. * * \note If \p rot is an EulerAngles (even when it represented as RotationBase explicitly), * angles ranges are __undefined__. - * Otherwise, Alpha and Gamma angles will be in the range [-PI, PI].
+ * Otherwise, alpha and gamma angles will be in the range [-PI, PI].
* As for Beta angle: * - If the system is Tait-Bryan, the beta angle will be in the range [-PI/2, PI/2]. * - otherwise: @@ -225,16 +227,20 @@ namespace Eigen return inverse(); } - /** Set \c *this from a rotation matrix(i.e. pure orthogonal matrix with determinant of +1). + /** Set \c *this from either: + * - a 3x3 rotation matrix expression(i.e. pure orthogonal matrix with determinant of +1), + * - a 3D vector expression representing Euler angles. * * See EulerAngles(const MatrixBase&) for more information about * angles ranges output. */ - template - EulerAngles& operator=(const MatrixBase& m) { - EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(Derived, 3, 3) + template + EulerAngles& operator=(const MatrixBase& other) + { + EIGEN_STATIC_ASSERT((internal::is_same::value), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) - System::CalcEulerAngles(*this, m); + internal::eulerangles_assign_impl::run(*this, other.derived()); return *this; } @@ -321,8 +327,29 @@ EIGEN_EULER_ANGLES_TYPEDEFS(double, d) { typedef _Scalar Scalar; }; + + // set from a rotation matrix + template + struct eulerangles_assign_impl + { + typedef typename Other::Scalar Scalar; + static void run(EulerAngles& e, const Other& m) + { + System::CalcEulerAngles(e, m); + } + }; + + // set from a vector of Euler angles + template + struct eulerangles_assign_impl + { + typedef typename Other::Scalar Scalar; + static void run(EulerAngles& e, const Other& vec) + { + e.angles() = vec; + } + }; } - } #endif // EIGEN_EULERANGLESCLASS_H diff --git a/unsupported/Eigen/src/EulerAngles/EulerSystem.h b/unsupported/Eigen/src/EulerAngles/EulerSystem.h index 0790e612f..28f52da61 100644 --- a/unsupported/Eigen/src/EulerAngles/EulerSystem.h +++ b/unsupported/Eigen/src/EulerAngles/EulerSystem.h @@ -36,6 +36,12 @@ namespace Eigen { enum { value = Axis != 0 && Abs::value <= 3 }; }; + + template + struct eulerangles_assign_impl; } #define EIGEN_EULER_ANGLES_CLASS_STATIC_ASSERT(COND,MSG) typedef char static_assertion_##MSG[(COND)?1:-1] @@ -269,6 +275,12 @@ namespace Eigen template friend class Eigen::EulerAngles; + + template + friend struct internal::eulerangles_assign_impl; }; #define EIGEN_EULER_SYSTEM_TYPEDEF(A, B, C) \ From 6a84246a6a5dd980e78e5609b3098cdbce93807c Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 21 Nov 2016 21:46:42 +0100 Subject: [PATCH 07/22] Fix regression in assigment of sparse block to spasre block. --- Eigen/src/SparseCore/SparseBlock.h | 12 +++++++----- test/sparse_block.cpp | 27 +++++++++++++++++++++++++++ 2 files changed, 34 insertions(+), 5 deletions(-) diff --git a/Eigen/src/SparseCore/SparseBlock.h b/Eigen/src/SparseCore/SparseBlock.h index 13e8b0bf1..acaf933f4 100644 --- a/Eigen/src/SparseCore/SparseBlock.h +++ b/Eigen/src/SparseCore/SparseBlock.h @@ -130,7 +130,7 @@ public: // 2 - let's check whether there is enough allocated memory Index nnz = tmp.nonZeros(); - Index start = m_outerStart==0 ? 0 : matrix.outerIndexPtr()[m_outerStart]; // starting position of the current block + Index start = m_outerStart==0 ? 0 : m_matrix.outerIndexPtr()[m_outerStart]; // starting position of the current block Index end = m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]; // ending position of the current block Index block_size = end - start; // available room in the current block Index tail_size = m_matrix.outerIndexPtr()[m_matrix.outerSize()] - end; @@ -139,6 +139,8 @@ public: ? Index(matrix.data().allocatedSize()) + block_size : block_size; + Index tmp_start = tmp.outerIndexPtr()[0]; + bool update_trailing_pointers = false; if(nnz>free_size) { @@ -148,8 +150,8 @@ public: internal::smart_copy(m_matrix.valuePtr(), m_matrix.valuePtr() + start, newdata.valuePtr()); internal::smart_copy(m_matrix.innerIndexPtr(), m_matrix.innerIndexPtr() + start, newdata.indexPtr()); - internal::smart_copy(tmp.valuePtr(), tmp.valuePtr() + nnz, newdata.valuePtr() + start); - internal::smart_copy(tmp.innerIndexPtr(), tmp.innerIndexPtr() + nnz, newdata.indexPtr() + start); + internal::smart_copy(tmp.valuePtr() + tmp_start, tmp.valuePtr() + tmp_start + nnz, newdata.valuePtr() + start); + internal::smart_copy(tmp.innerIndexPtr() + tmp_start, tmp.innerIndexPtr() + tmp_start + nnz, newdata.indexPtr() + start); internal::smart_copy(matrix.valuePtr()+end, matrix.valuePtr()+end + tail_size, newdata.valuePtr()+start+nnz); internal::smart_copy(matrix.innerIndexPtr()+end, matrix.innerIndexPtr()+end + tail_size, newdata.indexPtr()+start+nnz); @@ -173,8 +175,8 @@ public: update_trailing_pointers = true; } - internal::smart_copy(tmp.valuePtr(), tmp.valuePtr() + nnz, matrix.valuePtr() + start); - internal::smart_copy(tmp.innerIndexPtr(), tmp.innerIndexPtr() + nnz, matrix.innerIndexPtr() + start); + internal::smart_copy(tmp.valuePtr() + tmp_start, tmp.valuePtr() + tmp_start + nnz, matrix.valuePtr() + start); + internal::smart_copy(tmp.innerIndexPtr() + tmp_start, tmp.innerIndexPtr() + tmp_start + nnz, matrix.innerIndexPtr() + start); } // update outer index pointers and innerNonZeros diff --git a/test/sparse_block.cpp b/test/sparse_block.cpp index 49a5f135e..cb5213ade 100644 --- a/test/sparse_block.cpp +++ b/test/sparse_block.cpp @@ -223,6 +223,33 @@ template void sparse_block(const SparseMatrixType& re VERIFY_IS_APPROX(m2.block(r0,c0,r1,c1), refMat2.block(r0,c0,r1,c1)); VERIFY_IS_APPROX((2*m2).block(r0,c0,r1,c1), (2*refMat2).block(r0,c0,r1,c1)); + + if(m2.nonZeros()>0) + { + VERIFY_IS_APPROX(m2, refMat2); + SparseMatrixType m3(rows, cols); + DenseMatrix refMat3(rows, cols); refMat3.setZero(); + Index n = internal::random(1,10); + for(Index k=0; k(0,outer-1); + Index o2 = internal::random(0,outer-1); + if(SparseMatrixType::IsRowMajor) + { + m3.innerVector(o1) = m2.row(o2); + refMat3.row(o1) = refMat2.row(o2); + } + else + { + m3.innerVector(o1) = m2.col(o2); + refMat3.col(o1) = refMat2.col(o2); + } + if(internal::random()) + m3.makeCompressed(); + } + if(m3.nonZeros()>0) + VERIFY_IS_APPROX(m3, refMat3); + } } } From f3fb0a1940c93c2eea2342b20506d652050ff48b Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 22 Nov 2016 16:58:31 +0100 Subject: [PATCH 08/22] Disable usage of SSE3 haddpd that is extremely slow. --- Eigen/src/Core/arch/SSE/PacketMath.h | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index 6f31cf12b..3646abdb1 100755 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h @@ -504,6 +504,7 @@ template<> EIGEN_STRONG_INLINE Packet4f preduxp(const Packet4f* vecs) { return _mm_hadd_ps(_mm_hadd_ps(vecs[0], vecs[1]),_mm_hadd_ps(vecs[2], vecs[3])); } + template<> EIGEN_STRONG_INLINE Packet2d preduxp(const Packet2d* vecs) { return _mm_hadd_pd(vecs[0], vecs[1]); @@ -515,7 +516,6 @@ template<> EIGEN_STRONG_INLINE float predux(const Packet4f& a) return pfirst(_mm_hadd_ps(tmp0, tmp0)); } -template<> EIGEN_STRONG_INLINE double predux(const Packet2d& a) { return pfirst(_mm_hadd_pd(a, a)); } #else // SSE2 versions template<> EIGEN_STRONG_INLINE float predux(const Packet4f& a) @@ -523,10 +523,6 @@ template<> EIGEN_STRONG_INLINE float predux(const Packet4f& a) Packet4f tmp = _mm_add_ps(a, _mm_movehl_ps(a,a)); return pfirst(_mm_add_ss(tmp, _mm_shuffle_ps(tmp,tmp, 1))); } -template<> EIGEN_STRONG_INLINE double predux(const Packet2d& a) -{ - return pfirst(_mm_add_sd(a, _mm_unpackhi_pd(a,a))); -} template<> EIGEN_STRONG_INLINE Packet4f preduxp(const Packet4f* vecs) { @@ -548,6 +544,16 @@ template<> EIGEN_STRONG_INLINE Packet2d preduxp(const Packet2d* vecs) } #endif // SSE3 +template<> EIGEN_STRONG_INLINE double predux(const Packet2d& a) +{ + // Disable SSE3 _mm_hadd_pd that is extremely slow on all existing Intel's architectures + // (from Nehalem to Haswell) +// #ifdef EIGEN_VECTORIZE_SSE3 +// return pfirst(_mm_hadd_pd(a, a)); +// #else + return pfirst(_mm_add_sd(a, _mm_unpackhi_pd(a,a))); +// #endif +} #ifdef EIGEN_VECTORIZE_SSSE3 template<> EIGEN_STRONG_INLINE Packet4i preduxp(const Packet4i* vecs) From 7dd894e40e439a6d1f4aed659d1375d65589cff3 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 22 Nov 2016 21:41:30 +0100 Subject: [PATCH 09/22] Optimize predux (AVX) --- Eigen/src/Core/arch/AVX/PacketMath.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Eigen/src/Core/arch/AVX/PacketMath.h b/Eigen/src/Core/arch/AVX/PacketMath.h index e60ef307b..190e9a42b 100644 --- a/Eigen/src/Core/arch/AVX/PacketMath.h +++ b/Eigen/src/Core/arch/AVX/PacketMath.h @@ -401,8 +401,7 @@ template<> EIGEN_STRONG_INLINE float predux(const Packet8f& a) } template<> EIGEN_STRONG_INLINE double predux(const Packet4d& a) { - Packet4d tmp0 = _mm256_hadd_pd(a,_mm256_permute2f128_pd(a,a,1)); - return pfirst(_mm256_hadd_pd(tmp0,tmp0)); + return predux(_mm_add_pd(_mm256_castpd256_pd128(a),_mm256_extractf128_pd(a,1))); } template<> EIGEN_STRONG_INLINE Packet4f predux_downto4(const Packet8f& a) From 178c084856003f1cfd3020615ab98230d9520a80 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 22 Nov 2016 21:53:14 +0100 Subject: [PATCH 10/22] Disable usage of SSE3 _mm_hadd_ps that is extremely slow. --- Eigen/src/Core/arch/SSE/PacketMath.h | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index 3646abdb1..80cf8af09 100755 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h @@ -510,20 +510,7 @@ template<> EIGEN_STRONG_INLINE Packet2d preduxp(const Packet2d* vecs) return _mm_hadd_pd(vecs[0], vecs[1]); } -template<> EIGEN_STRONG_INLINE float predux(const Packet4f& a) -{ - Packet4f tmp0 = _mm_hadd_ps(a,a); - return pfirst(_mm_hadd_ps(tmp0, tmp0)); -} - #else -// SSE2 versions -template<> EIGEN_STRONG_INLINE float predux(const Packet4f& a) -{ - Packet4f tmp = _mm_add_ps(a, _mm_movehl_ps(a,a)); - return pfirst(_mm_add_ss(tmp, _mm_shuffle_ps(tmp,tmp, 1))); -} - template<> EIGEN_STRONG_INLINE Packet4f preduxp(const Packet4f* vecs) { Packet4f tmp0, tmp1, tmp2; @@ -544,6 +531,19 @@ template<> EIGEN_STRONG_INLINE Packet2d preduxp(const Packet2d* vecs) } #endif // SSE3 +template<> EIGEN_STRONG_INLINE float predux(const Packet4f& a) +{ + // Disable SSE3 _mm_hadd_pd that is extremely slow on all existing Intel's architectures + // (from Nehalem to Haswell) +// #ifdef EIGEN_VECTORIZE_SSE3 +// Packet4f tmp = _mm_add_ps(a, vec4f_swizzle1(a,2,3,2,3)); +// return pfirst(_mm_hadd_ps(tmp, tmp)); +// #else + Packet4f tmp = _mm_add_ps(a, _mm_movehl_ps(a,a)); + return pfirst(_mm_add_ss(tmp, _mm_shuffle_ps(tmp,tmp, 1))); +// #endif +} + template<> EIGEN_STRONG_INLINE double predux(const Packet2d& a) { // Disable SSE3 _mm_hadd_pd that is extremely slow on all existing Intel's architectures From 74637fa4e39afc0c6784588e9545b454bb5ea407 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 22 Nov 2016 21:57:52 +0100 Subject: [PATCH 11/22] Optimize predux (AVX) --- Eigen/src/Core/arch/AVX/PacketMath.h | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/Eigen/src/Core/arch/AVX/PacketMath.h b/Eigen/src/Core/arch/AVX/PacketMath.h index 190e9a42b..58b5dc5ba 100644 --- a/Eigen/src/Core/arch/AVX/PacketMath.h +++ b/Eigen/src/Core/arch/AVX/PacketMath.h @@ -395,9 +395,7 @@ template<> EIGEN_STRONG_INLINE Packet4d preduxp(const Packet4d* vecs) template<> EIGEN_STRONG_INLINE float predux(const Packet8f& a) { - Packet8f tmp0 = _mm256_hadd_ps(a,_mm256_permute2f128_ps(a,a,1)); - tmp0 = _mm256_hadd_ps(tmp0,tmp0); - return pfirst(_mm256_hadd_ps(tmp0, tmp0)); + return predux(_mm_add_ps(_mm256_castps256_ps128(a),_mm256_extractf128_ps(a,1))); } template<> EIGEN_STRONG_INLINE double predux(const Packet4d& a) { From a91de27e983d9f752eb9745be0a53f145eb23d5b Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 23 Nov 2016 12:24:48 +0100 Subject: [PATCH 12/22] Fix compilation issue with MSVC: MSVC always messes up with shadowed template arguments, for instance in: struct B { typedef float T; } template struct A : B { T g; }; The type of A::g will be float and not double. --- Eigen/src/SparseCore/SparseCwiseBinaryOp.h | 52 +++++++++++----------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h index 04cef66fc..4ba4d631d 100644 --- a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h +++ b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h @@ -394,10 +394,10 @@ struct sparse_conjunction_evaluator { protected: typedef typename XprType::Functor BinaryOp; - typedef typename XprType::Lhs Lhs; - typedef typename XprType::Rhs Rhs; - typedef typename evaluator::InnerIterator LhsIterator; - typedef typename evaluator::InnerIterator RhsIterator; + typedef typename XprType::Lhs LhsArg; + typedef typename XprType::Rhs RhsArg; + typedef typename evaluator::InnerIterator LhsIterator; + typedef typename evaluator::InnerIterator RhsIterator; typedef typename XprType::StorageIndex StorageIndex; typedef typename traits::Scalar Scalar; public: @@ -449,7 +449,7 @@ public: enum { - CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, + CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, Flags = XprType::Flags }; @@ -468,8 +468,8 @@ public: protected: const BinaryOp m_functor; - evaluator m_lhsImpl; - evaluator m_rhsImpl; + evaluator m_lhsImpl; + evaluator m_rhsImpl; }; // "dense ^ sparse" @@ -479,10 +479,10 @@ struct sparse_conjunction_evaluator { protected: typedef typename XprType::Functor BinaryOp; - typedef typename XprType::Lhs Lhs; - typedef typename XprType::Rhs Rhs; - typedef evaluator LhsEvaluator; - typedef typename evaluator::InnerIterator RhsIterator; + typedef typename XprType::Lhs LhsArg; + typedef typename XprType::Rhs RhsArg; + typedef evaluator LhsEvaluator; + typedef typename evaluator::InnerIterator RhsIterator; typedef typename XprType::StorageIndex StorageIndex; typedef typename traits::Scalar Scalar; public: @@ -490,7 +490,7 @@ public: class ReverseInnerIterator; class InnerIterator { - enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit }; + enum { IsRowMajor = (int(RhsArg::Flags)&RowMajorBit)==RowMajorBit }; public: @@ -522,9 +522,9 @@ public: enum { - CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, + CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, // Expose storage order of the sparse expression - Flags = (XprType::Flags & ~RowMajorBit) | (int(Rhs::Flags)&RowMajorBit) + Flags = (XprType::Flags & ~RowMajorBit) | (int(RhsArg::Flags)&RowMajorBit) }; explicit sparse_conjunction_evaluator(const XprType& xpr) @@ -542,8 +542,8 @@ public: protected: const BinaryOp m_functor; - evaluator m_lhsImpl; - evaluator m_rhsImpl; + evaluator m_lhsImpl; + evaluator m_rhsImpl; }; // "sparse ^ dense" @@ -553,10 +553,10 @@ struct sparse_conjunction_evaluator { protected: typedef typename XprType::Functor BinaryOp; - typedef typename XprType::Lhs Lhs; - typedef typename XprType::Rhs Rhs; - typedef typename evaluator::InnerIterator LhsIterator; - typedef evaluator RhsEvaluator; + typedef typename XprType::Lhs LhsArg; + typedef typename XprType::Rhs RhsArg; + typedef typename evaluator::InnerIterator LhsIterator; + typedef evaluator RhsEvaluator; typedef typename XprType::StorageIndex StorageIndex; typedef typename traits::Scalar Scalar; public: @@ -564,7 +564,7 @@ public: class ReverseInnerIterator; class InnerIterator { - enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit }; + enum { IsRowMajor = (int(LhsArg::Flags)&RowMajorBit)==RowMajorBit }; public: @@ -590,16 +590,16 @@ public: protected: LhsIterator m_lhsIter; - const evaluator &m_rhsEval; + const evaluator &m_rhsEval; const BinaryOp& m_functor; const Index m_outer; }; enum { - CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, + CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, // Expose storage order of the sparse expression - Flags = (XprType::Flags & ~RowMajorBit) | (int(Lhs::Flags)&RowMajorBit) + Flags = (XprType::Flags & ~RowMajorBit) | (int(LhsArg::Flags)&RowMajorBit) }; explicit sparse_conjunction_evaluator(const XprType& xpr) @@ -617,8 +617,8 @@ public: protected: const BinaryOp m_functor; - evaluator m_lhsImpl; - evaluator m_rhsImpl; + evaluator m_lhsImpl; + evaluator m_rhsImpl; }; } From e340866c816504bf47780c559b6b3678db8b14e0 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 23 Nov 2016 14:04:57 +0100 Subject: [PATCH 13/22] Fix compilation with gcc and old ABI version --- Eigen/src/Core/arch/AVX/PacketMath.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/arch/AVX/PacketMath.h b/Eigen/src/Core/arch/AVX/PacketMath.h index 58b5dc5ba..195d40fb4 100644 --- a/Eigen/src/Core/arch/AVX/PacketMath.h +++ b/Eigen/src/Core/arch/AVX/PacketMath.h @@ -395,11 +395,11 @@ template<> EIGEN_STRONG_INLINE Packet4d preduxp(const Packet4d* vecs) template<> EIGEN_STRONG_INLINE float predux(const Packet8f& a) { - return predux(_mm_add_ps(_mm256_castps256_ps128(a),_mm256_extractf128_ps(a,1))); + return predux(Packet4f(_mm_add_ps(_mm256_castps256_ps128(a),_mm256_extractf128_ps(a,1)))); } template<> EIGEN_STRONG_INLINE double predux(const Packet4d& a) { - return predux(_mm_add_pd(_mm256_castpd256_pd128(a),_mm256_extractf128_pd(a,1))); + return predux(Packet2d(_mm_add_pd(_mm256_castpd256_pd128(a),_mm256_extractf128_pd(a,1)))); } template<> EIGEN_STRONG_INLINE Packet4f predux_downto4(const Packet8f& a) From 9246587122440969d0585dab8c2cc2e45c0a0ec0 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 23 Nov 2016 15:42:26 +0100 Subject: [PATCH 14/22] Patch from Oleg Shirokobrod to extend polynomial solver to complexes --- unsupported/Eigen/src/Polynomials/Companion.h | 50 +++++++++---------- .../Eigen/src/Polynomials/PolynomialSolver.h | 4 +- 2 files changed, 27 insertions(+), 27 deletions(-) diff --git a/unsupported/Eigen/src/Polynomials/Companion.h b/unsupported/Eigen/src/Polynomials/Companion.h index b515c2920..e0af6ebe4 100644 --- a/unsupported/Eigen/src/Polynomials/Companion.h +++ b/unsupported/Eigen/src/Polynomials/Companion.h @@ -75,7 +75,7 @@ class companion void setPolynomial( const VectorType& poly ) { const Index deg = poly.size()-1; - m_monic = -1/poly[deg] * poly.head(deg); + m_monic = Scalar(-1)/poly[deg] * poly.head(deg); //m_bl_diag.setIdentity( deg-1 ); m_bl_diag.setOnes(deg-1); } @@ -107,8 +107,8 @@ class companion * colB and rowB are repectively the multipliers for * the column and the row in order to balance them. * */ - bool balanced( Scalar colNorm, Scalar rowNorm, - bool& isBalanced, Scalar& colB, Scalar& rowB ); + bool balanced( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& rowB ); /** Helper function for the balancing algorithm. * \returns true if the row and the column, having colNorm and rowNorm @@ -116,8 +116,8 @@ class companion * colB and rowB are repectively the multipliers for * the column and the row in order to balance them. * */ - bool balancedR( Scalar colNorm, Scalar rowNorm, - bool& isBalanced, Scalar& colB, Scalar& rowB ); + bool balancedR( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& rowB ); public: /** @@ -139,10 +139,10 @@ class companion template< typename _Scalar, int _Deg > inline -bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm, - bool& isBalanced, Scalar& colB, Scalar& rowB ) +bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& rowB ) { - if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; } + if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm ){ return true; } else { //To find the balancing coefficients, if the radix is 2, @@ -150,29 +150,29 @@ bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm, // \f$ 2^{2\sigma-1} < rowNorm / colNorm \le 2^{2\sigma+1} \f$ // then the balancing coefficient for the row is \f$ 1/2^{\sigma} \f$ // and the balancing coefficient for the column is \f$ 2^{\sigma} \f$ - rowB = rowNorm / radix(); - colB = Scalar(1); - const Scalar s = colNorm + rowNorm; + rowB = rowNorm / radix(); + colB = RealScalar(1); + const RealScalar s = colNorm + rowNorm; while (colNorm < rowB) { - colB *= radix(); - colNorm *= radix2(); + colB *= radix(); + colNorm *= radix2(); } - rowB = rowNorm * radix(); + rowB = rowNorm * radix(); while (colNorm >= rowB) { - colB /= radix(); - colNorm /= radix2(); + colB /= radix(); + colNorm /= radix2(); } //This line is used to avoid insubstantial balancing - if ((rowNorm + colNorm) < Scalar(0.95) * s * colB) + if ((rowNorm + colNorm) < RealScalar(0.95) * s * colB) { isBalanced = false; - rowB = Scalar(1) / colB; + rowB = RealScalar(1) / colB; return false; } else{ @@ -182,21 +182,21 @@ bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm, template< typename _Scalar, int _Deg > inline -bool companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm, - bool& isBalanced, Scalar& colB, Scalar& rowB ) +bool companion<_Scalar,_Deg>::balancedR( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& rowB ) { - if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; } + if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm ){ return true; } else { /** * Set the norm of the column and the row to the geometric mean * of the row and column norm */ - const _Scalar q = colNorm/rowNorm; + const RealScalar q = colNorm/rowNorm; if( !isApprox( q, _Scalar(1) ) ) { rowB = sqrt( colNorm/rowNorm ); - colB = Scalar(1)/rowB; + colB = RealScalar(1)/rowB; isBalanced = false; return false; @@ -219,8 +219,8 @@ void companion<_Scalar,_Deg>::balance() while( !hasConverged ) { hasConverged = true; - Scalar colNorm,rowNorm; - Scalar colB,rowB; + RealScalar colNorm,rowNorm; + RealScalar colB,rowB; //First row, first column excluding the diagonal //============================================== diff --git a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h index 03198ec8e..3ba992cb4 100644 --- a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h +++ b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h @@ -327,7 +327,7 @@ class PolynomialSolverBase * However, almost always, correct accuracy is reached even in these cases for 64bit * (double) floating types and small polynomial degree (<20). */ -template< typename _Scalar, int _Deg > +template< typename _Scalar, int _Deg , typename EigenSolverType = EigenSolver > > class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> { public: @@ -337,7 +337,7 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) typedef Matrix CompanionMatrixType; - typedef EigenSolver EigenSolverType; + //typedef EigenSolver EigenSolverType; public: /** Computes the complex roots of a new polynomial. */ From 56e5ec07c6846585a8d808cab15a0aa5f6cd0cc8 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 23 Nov 2016 16:05:10 +0100 Subject: [PATCH 15/22] Automatically switch between EigenSolver and ComplexEigenSolver, and fix a few Real versus Scalar issues. --- .../Eigen/src/Polynomials/PolynomialSolver.h | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h index 3ba992cb4..788594247 100644 --- a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h +++ b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h @@ -99,7 +99,7 @@ class PolynomialSolverBase */ inline const RootType& greatestRoot() const { - std::greater greater; + std::greater greater; return selectComplexRoot_withRespectToNorm( greater ); } @@ -108,7 +108,7 @@ class PolynomialSolverBase */ inline const RootType& smallestRoot() const { - std::less less; + std::less less; return selectComplexRoot_withRespectToNorm( less ); } @@ -213,7 +213,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const { - std::greater greater; + std::greater greater; return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold ); } @@ -236,7 +236,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const { - std::less less; + std::less less; return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold ); } @@ -259,7 +259,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const { - std::greater greater; + std::greater greater; return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold ); } @@ -282,7 +282,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const { - std::less less; + std::less less; return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold ); } @@ -327,7 +327,7 @@ class PolynomialSolverBase * However, almost always, correct accuracy is reached even in these cases for 64bit * (double) floating types and small polynomial degree (<20). */ -template< typename _Scalar, int _Deg , typename EigenSolverType = EigenSolver > > +template class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> { public: @@ -337,7 +337,9 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) typedef Matrix CompanionMatrixType; - //typedef EigenSolver EigenSolverType; + typedef typename internal::conditional::IsComplex, + ComplexEigenSolver, + EigenSolver >::type EigenSolverType; public: /** Computes the complex roots of a new polynomial. */ From f12b368417992f0974678646f2fb7fa2db44b633 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 23 Nov 2016 16:05:45 +0100 Subject: [PATCH 16/22] Extend polynomial solver unit tests to complexes --- unsupported/test/polynomialsolver.cpp | 34 +++++++++++++-------------- 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/unsupported/test/polynomialsolver.cpp b/unsupported/test/polynomialsolver.cpp index 0c87478dd..7ad4aa69d 100644 --- a/unsupported/test/polynomialsolver.cpp +++ b/unsupported/test/polynomialsolver.cpp @@ -32,9 +32,10 @@ bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve ) { typedef typename POLYNOMIAL::Index Index; typedef typename POLYNOMIAL::Scalar Scalar; + typedef typename POLYNOMIAL::RealScalar RealScalar; typedef typename SOLVER::RootsType RootsType; - typedef Matrix EvalRootsType; + typedef Matrix EvalRootsType; const Index deg = pols.size()-1; @@ -57,7 +58,7 @@ bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve ) cerr << endl; } - std::vector rootModuli( roots.size() ); + std::vector rootModuli( roots.size() ); Map< EvalRootsType > aux( &rootModuli[0], roots.size() ); aux = roots.array().abs(); std::sort( rootModuli.begin(), rootModuli.end() ); @@ -83,7 +84,7 @@ void evalSolver( const POLYNOMIAL& pols ) { typedef typename POLYNOMIAL::Scalar Scalar; - typedef PolynomialSolver PolynomialSolverType; + typedef PolynomialSolver PolynomialSolverType; PolynomialSolverType psolve; aux_evalSolver( pols, psolve ); @@ -97,6 +98,7 @@ void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const { using std::sqrt; typedef typename POLYNOMIAL::Scalar Scalar; + typedef typename POLYNOMIAL::RealScalar RealScalar; typedef PolynomialSolver PolynomialSolverType; @@ -107,15 +109,12 @@ void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const // 1) the roots found are correct // 2) the roots have distinct moduli - typedef typename POLYNOMIAL::Scalar Scalar; - typedef typename REAL_ROOTS::Scalar Real; - //Test realRoots - std::vector< Real > calc_realRoots; - psolve.realRoots( calc_realRoots ); - VERIFY( calc_realRoots.size() == (size_t)real_roots.size() ); + std::vector< RealScalar > calc_realRoots; + psolve.realRoots( calc_realRoots, test_precision()); + VERIFY_IS_EQUAL( calc_realRoots.size() , (size_t)real_roots.size() ); - const Scalar psPrec = sqrt( test_precision() ); + const RealScalar psPrec = sqrt( test_precision() ); for( size_t i=0; i 0 ) ); if( hasRealRoot ){ VERIFY( internal::isApprox( real_roots.array().abs().maxCoeff(), abs(r), psPrec ) ); } @@ -167,9 +166,11 @@ void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const template void polynomialsolver(int deg) { - typedef internal::increment_if_fixed_size<_Deg> Dim; + typedef typename NumTraits<_Scalar>::Real RealScalar; + typedef internal::increment_if_fixed_size<_Deg> Dim; typedef Matrix<_Scalar,Dim::ret,1> PolynomialType; typedef Matrix<_Scalar,_Deg,1> EvalRootsType; + typedef Matrix RealRootsType; cout << "Standard cases" << endl; PolynomialType pols = PolynomialType::Random(deg+1); @@ -182,15 +183,11 @@ void polynomialsolver(int deg) evalSolver<_Deg,PolynomialType>( pols ); cout << "Test sugar" << endl; - EvalRootsType realRoots = EvalRootsType::Random(deg); + RealRootsType realRoots = RealRootsType::Random(deg); roots_to_monicPolynomial( realRoots, pols ); evalSolverSugarFunction<_Deg>( pols, - realRoots.template cast < - std::complex< - typename NumTraits<_Scalar>::Real - > - >(), + realRoots.template cast >().eval(), realRoots ); } @@ -214,5 +211,6 @@ void test_polynomialsolver() internal::random(9,13) )) ); CALL_SUBTEST_11((polynomialsolver(1)) ); + CALL_SUBTEST_12((polynomialsolver,Dynamic>(internal::random(2,13))) ); } } From 21d0286d81bf98242a3b5211622f328c7a641661 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 23 Nov 2016 22:15:03 +0100 Subject: [PATCH 17/22] bug #1348: Document EIGEN_MAX_ALIGN_BYTES and EIGEN_MAX_STATIC_ALIGN_BYTES, and reflect in the doc that EIGEN_DONT_ALIGN* are deprecated. --- doc/PreprocessorDirectives.dox | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/doc/PreprocessorDirectives.dox b/doc/PreprocessorDirectives.dox index 2f9c4c370..c4b9903c1 100644 --- a/doc/PreprocessorDirectives.dox +++ b/doc/PreprocessorDirectives.dox @@ -97,31 +97,35 @@ run time. However, these assertions do cost time and can thus be turned off. \section TopicPreprocessorDirectivesPerformance Alignment, vectorization and performance tweaking - - \b EIGEN_MALLOC_ALREADY_ALIGNED - Can be set to 0 or 1 to tell whether default system \c malloc already + - \b \c EIGEN_MALLOC_ALREADY_ALIGNED - Can be set to 0 or 1 to tell whether default system \c malloc already returns aligned buffers. In not defined, then this information is automatically deduced from the compiler and system preprocessor tokens. - - \b EIGEN_DONT_ALIGN - disables alignment completely. %Eigen will not try to align its objects and does not - expect that any objects passed to it are aligned. This will turn off vectorization. Not defined by default. - - \b EIGEN_DONT_ALIGN_STATICALLY - disables alignment of arrays on the stack. Not defined by default, unless - \c EIGEN_DONT_ALIGN is defined. - - \b EIGEN_DONT_PARALLELIZE - if defined, this disables multi-threading. This is only relevant if you enabled OpenMP. + - \b \c EIGEN_MAX_ALIGN_BYTES - Must be a power of two, or 0. Defines an upper bound on the memory boundary in bytes on which dynamically and statically allocated data may be aligned by %Eigen. If not defined, a default value is automatically computed based on architecture, compiler, and OS. + This option is typically used to enforce binary compatibility between code/libraries compiled with different SIMD options. For instance, one may compile AVX code and enforce ABI compatibility with existing SSE code by defining \c EIGEN_MAX_ALIGN_BYTES=16. In the other way round, since by default AVX implies 32 bytes alignment for best performance, one can compile SSE code to be ABI compatible with AVX code by defining \c EIGEN_MAX_ALIGN_BYTES=32. + - \b \c EIGEN_MAX_STATIC_ALIGN_BYTES - Same as \c EIGEN_MAX_ALIGN_BYTES but for statically allocated data only. By default, if only \c EIGEN_MAX_ALIGN_BYTES is defined, then \c EIGEN_MAX_STATIC_ALIGN_BYTES == \c EIGEN_MAX_ALIGN_BYTES, otherwise a default value is automatically computed based on architecture, compiler, and OS (can be smaller than the default value of EIGEN_MAX_ALIGN_BYTES on architectures that do not support stack alignment). + Let us emphasize that \c EIGEN_MAX_*_ALIGN_BYTES define only a diserable upper bound. In practice data is aligned to largest power-of-two common divisor of \c EIGEN_MAX_STATIC_ALIGN_BYTES and the size of the data, such that memory is not wasted. + - \b \c EIGEN_DONT_PARALLELIZE - if defined, this disables multi-threading. This is only relevant if you enabled OpenMP. See \ref TopicMultiThreading for details. - \b EIGEN_DONT_VECTORIZE - disables explicit vectorization when defined. Not defined by default, unless alignment is disabled by %Eigen's platform test or the user defining \c EIGEN_DONT_ALIGN. - - \b EIGEN_UNALIGNED_VECTORIZE - disables/enables vectorization with unaligned stores. Default is 1 (enabled). + - \b \c EIGEN_UNALIGNED_VECTORIZE - disables/enables vectorization with unaligned stores. Default is 1 (enabled). If set to 0 (disabled), then expression for which the destination cannot be aligned are not vectorized (e.g., unaligned small fixed size vectors or matrices) - - \b EIGEN_FAST_MATH - enables some optimizations which might affect the accuracy of the result. This currently + - \b \c EIGEN_FAST_MATH - enables some optimizations which might affect the accuracy of the result. This currently enables the SSE vectorization of sin() and cos(), and speedups sqrt() for single precision. Defined to 1 by default. Define it to 0 to disable. - - \b EIGEN_UNROLLING_LIMIT - defines the size of a loop to enable meta unrolling. Set it to zero to disable + - \b \c EIGEN_UNROLLING_LIMIT - defines the size of a loop to enable meta unrolling. Set it to zero to disable unrolling. The size of a loop here is expressed in %Eigen's own notion of "number of FLOPS", it does not correspond to the number of iterations or the number of instructions. The default is value 100. - - \b EIGEN_STACK_ALLOCATION_LIMIT - defines the maximum bytes for a buffer to be allocated on the stack. For internal + - \b \c EIGEN_STACK_ALLOCATION_LIMIT - defines the maximum bytes for a buffer to be allocated on the stack. For internal temporary buffers, dynamic memory allocation is employed as a fall back. For fixed-size matrices or arrays, exceeding this threshold raises a compile time assertion. Use 0 to set no limit. Default is 128 KB. + - \c EIGEN_DONT_ALIGN - Deprecated, it is a synonym for \c EIGEN_MAX_ALIGN_BYTES=0. It disables alignment completely. %Eigen will not try to align its objects and does not expect that any objects passed to it are aligned. This will turn off vectorization if \b EIGEN_UNALIGNED_VECTORIZE=1. Not defined by default. + - \c EIGEN_DONT_ALIGN_STATICALLY - Deprecated, it is a synonym for \c EIGEN_MAX_STATIC_ALIGN_BYTES=0. It disables alignment of arrays on the stack. Not defined by default, unless \c EIGEN_DONT_ALIGN is defined. + + \section TopicPreprocessorDirectivesPlugins Plugins It is possible to add new methods to many fundamental classes in %Eigen by writing a plugin. As explained in From 308961c05eeaeb49efccc1d0b915f39de5ee02a8 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 23 Nov 2016 22:17:52 +0100 Subject: [PATCH 18/22] Fix compilation. --- unsupported/Eigen/src/EulerAngles/EulerAngles.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unsupported/Eigen/src/EulerAngles/EulerAngles.h b/unsupported/Eigen/src/EulerAngles/EulerAngles.h index 2a12c8da3..a5d034d71 100644 --- a/unsupported/Eigen/src/EulerAngles/EulerAngles.h +++ b/unsupported/Eigen/src/EulerAngles/EulerAngles.h @@ -292,7 +292,7 @@ namespace Eigen EulerAngles cast() const { EulerAngles e; - e.angles() = angles().cast(); + e.angles() = angles().template cast(); return e; } }; From 3be1afca114532f69e1d58c2a7232f05ab44b9a0 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Wed, 23 Nov 2016 18:49:51 -0800 Subject: [PATCH 19/22] Disabled the "remove the call to 'std::abs' since unsigned values cannot be negative" warning introduced in clang 3.5 --- Eigen/src/Core/util/DisableStupidWarnings.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Eigen/src/Core/util/DisableStupidWarnings.h b/Eigen/src/Core/util/DisableStupidWarnings.h index 7559e129c..21f80d86b 100755 --- a/Eigen/src/Core/util/DisableStupidWarnings.h +++ b/Eigen/src/Core/util/DisableStupidWarnings.h @@ -42,6 +42,9 @@ #pragma clang diagnostic push #endif #pragma clang diagnostic ignored "-Wconstant-logical-operand" + #if __clang_major__ >= 3 && __clang_minor__ >= 5 + #pragma clang diagnostic ignored "-Wabsolute-value" + #endif #elif defined __GNUC__ && __GNUC__>=6 From 7ad37606dd256f049060aa3cc535a04b924c35c8 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 24 Nov 2016 12:31:43 -0800 Subject: [PATCH 20/22] Fixed the documentation of Scalar Tensors --- unsupported/Eigen/CXX11/src/Tensor/README.md | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/unsupported/Eigen/CXX11/src/Tensor/README.md b/unsupported/Eigen/CXX11/src/Tensor/README.md index 02146527b..fbb7f3bfc 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/README.md +++ b/unsupported/Eigen/CXX11/src/Tensor/README.md @@ -1737,11 +1737,9 @@ TODO ## Representation of scalar values -Scalar values are often represented by tensors of size 1 and rank 1. It would be -more logical and user friendly to use tensors of rank 0 instead. For example -Tensor::maximum() currently returns a Tensor. Similarly, the inner -product of 2 1d tensors (through contractions) returns a 1d tensor. In the -future these operations might be updated to return 0d tensors instead. +Scalar values are often represented by tensors of size 1 and rank 0.For example +Tensor::maximum() currently returns a Tensor. Similarly, the inner +product of 2 1d tensors (through contractions) returns a 0d tensor. ## Limitations From 7fe704596a2cfa71eb7eae7fd91bc105f038b13f Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Sat, 26 Nov 2016 11:26:07 -0800 Subject: [PATCH 21/22] Added missing array_get method for numeric_list --- unsupported/Eigen/CXX11/src/util/CXX11Meta.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/unsupported/Eigen/CXX11/src/util/CXX11Meta.h b/unsupported/Eigen/CXX11/src/util/CXX11Meta.h index ec27eddb8..63c2a1def 100644 --- a/unsupported/Eigen/CXX11/src/util/CXX11Meta.h +++ b/unsupported/Eigen/CXX11/src/util/CXX11Meta.h @@ -123,6 +123,10 @@ template struct get<0, type_lis template struct get> : get> {}; template struct get<0, numeric_list> { constexpr static T value = a; }; +template constexpr inline const T array_get(const numeric_list& l) { + return get<(int)n, numeric_list>::value; +} + /* always get type, regardless of dummy; good for parameter pack expansion */ template struct id_numeric { typedef t type; }; From 67b2c41f30a29debcb720fe85c2581901ff36fd2 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Sat, 26 Nov 2016 11:27:29 -0800 Subject: [PATCH 22/22] Avoided unnecessary type conversion --- unsupported/Eigen/CXX11/src/Tensor/TensorDimensions.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorDimensions.h b/unsupported/Eigen/CXX11/src/Tensor/TensorDimensions.h index b24cdebf1..a0f4c9cb7 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorDimensions.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorDimensions.h @@ -33,7 +33,7 @@ namespace Eigen { namespace internal { template struct dget { - static const std::size_t value = get::value; + static const std::ptrdiff_t value = get::value; };