update eigen

9 files changed, 337 insertions, 289 deletions

diff --git a/eigen/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h b/eigen/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h
index d280886..279fe5c 100644
--- a/eigen/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h
+++ b/eigen/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h
@@ -683,4 +683,11 @@ template<typename DerType> struct NumTraits<AutoDiffScalar<DerType> >
 
 }
 
+namespace std {
+template <typename T>
+class numeric_limits<Eigen::AutoDiffScalar<T> >
+  : public numeric_limits<typename T::Scalar> {};
+
+}  // namespace std
+
 #endif // EIGEN_AUTODIFF_SCALAR_H
diff --git a/eigen/unsupported/Eigen/src/EulerAngles/EulerAngles.h b/eigen/unsupported/Eigen/src/EulerAngles/EulerAngles.h
index a5d034d..13a0da1 100644
--- a/eigen/unsupported/Eigen/src/EulerAngles/EulerAngles.h
+++ b/eigen/unsupported/Eigen/src/EulerAngles/EulerAngles.h
@@ -12,6 +12,11 @@
 
 namespace Eigen
 {
+  /*template<typename Other,
+         int OtherRows=Other::RowsAtCompileTime,
+         int OtherCols=Other::ColsAtCompileTime>
+  struct ei_eulerangles_assign_impl;*/
+
   /** \class EulerAngles
     *
     * \ingroup EulerAngles_Module
@@ -31,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 single representation (e.g. unlike rotation matrices).
+    *  by Euler angles, but there is no singular 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).
     *
@@ -50,27 +55,33 @@ namespace Eigen
     * Additionally, some axes related computation is done in compile time.
     *
     * #### Euler angles ranges in conversions ####
-    * Rotations representation as EulerAngles are not single (unlike matrices),
-    *  and even have infinite EulerAngles representations.<BR>
-    * For example, add or subtract 2*PI from either angle of EulerAngles
-    *  and you'll get the same rotation.
-    * 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:
-    * - 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].<BR>
-    *   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]
+    * When converting some rotation to Euler angles, there are some ways you can guarantee
+    *  the Euler angles ranges.
     *
+    * #### implicit ranges ####
+    * When using implicit ranges, all angles are guarantee to be in the range [-PI, +PI],
+    *  unless you convert from some other Euler angles.
+    * In this case, the range is __undefined__ (might be even less than -PI or greater than +2*PI).
     * \sa EulerAngles(const MatrixBase<Derived>&)
     * \sa EulerAngles(const RotationBase<Derived, 3>&)
     *
+    * #### explicit ranges ####
+    * When using explicit ranges, all angles are guarantee to be in the range you choose.
+    * In the range Boolean parameter, you're been ask whether you prefer the positive range or not:
+    * - _true_ - force the range between [0, +2*PI]
+    * - _false_ - force the range between [-PI, +PI]
+    *
+    * ##### compile time ranges #####
+    * This is when you have compile time ranges and you prefer to
+    *  use template parameter. (e.g. for performance)
+    * \sa FromRotation()
+    *
+    * ##### run-time time ranges #####
+    * Run-time ranges are also supported.
+    * \sa EulerAngles(const MatrixBase<Derived>&, bool, bool, bool)
+    * \sa EulerAngles(const RotationBase<Derived, 3>&, bool, bool, bool)
+    *
     * ### Convenient user typedefs ###
     *
     * Convenient typedefs for EulerAngles exist for float and double scalar,
@@ -92,7 +103,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.
     */
@@ -100,11 +111,8 @@ namespace Eigen
   class EulerAngles : public RotationBase<EulerAngles<_Scalar, _System>, 3>
   {
     public:
-      typedef RotationBase<EulerAngles<_Scalar, _System>, 3> Base;
-      
       /** the scalar type of the angles */
       typedef _Scalar Scalar;
-      typedef typename NumTraits<Scalar>::Real RealScalar;
       
       /** the EulerSystem to use, which represents the axes of rotation. */
       typedef _System System;
@@ -138,56 +146,67 @@ namespace Eigen
     public:
       /** Default constructor without initialization. */
       EulerAngles() {}
-      /** Constructs and initialize an EulerAngles (\p alpha, \p beta, \p gamma). */
+      /** Constructs and initialize Euler angles(\p alpha, \p beta, \p gamma). */
       EulerAngles(const Scalar& alpha, const Scalar& beta, const Scalar& gamma) :
         m_angles(alpha, beta, gamma) {}
       
-      // 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 initialize Euler angles from a 3x3 rotation matrix \p m.
+        *
+        * \note All angles will be in the range [-PI, PI].
+      */
+      template<typename Derived>
+      EulerAngles(const MatrixBase<Derived>& m) { *this = m; }
       
-      /** 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.
+      /** 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].
+        * Otherwise, the specified angle will be in the range [-PI, +PI].
         *
-        * \note If \p other is a 3x3 rotation matrix, the angles range rules will be as follow:<BR>
-        *  Alpha and gamma angles will be in the range [-PI, PI].<BR>
-        *  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]
-       */
+        * \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<typename Derived>
-      explicit EulerAngles(const MatrixBase<Derived>& other) { *this = other; }
+      EulerAngles(
+        const MatrixBase<Derived>& m,
+        bool positiveRangeAlpha,
+        bool positiveRangeBeta,
+        bool positiveRangeGamma) {
+        
+        System::CalcEulerAngles(*this, m, positiveRangeAlpha, positiveRangeBeta, positiveRangeGamma);
+      }
       
       /** 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].<BR>
-        *  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]
+        * \note All angles will be in the range [-PI, PI], unless \p rot is an EulerAngles.
+        *  If rot is an EulerAngles, expected EulerAngles range is __undefined__.
+        *  (Use other functions here for enforcing range if this effect is desired)
       */
       template<typename Derived>
-      EulerAngles(const RotationBase<Derived, 3>& rot) { System::CalcEulerAngles(*this, rot.toRotationMatrix()); }
+      EulerAngles(const RotationBase<Derived, 3>& rot) { *this = rot; }
       
-      /*EulerAngles(const QuaternionType& q)
-      {
-        // TODO: Implement it in a faster way for quaternions
-        // According to http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/
-        //  we can compute only the needed matrix cells and then convert to euler angles. (see ZYX example below)
-        // Currently we compute all matrix cells from quaternion.
-
-        // Special case only for ZYX
-        //Scalar y2 = q.y() * q.y();
-        //m_angles[0] = std::atan2(2*(q.w()*q.z() + q.x()*q.y()), (1 - 2*(y2 + q.z()*q.z())));
-        //m_angles[1] = std::asin( 2*(q.w()*q.y() - q.z()*q.x()));
-        //m_angles[2] = std::atan2(2*(q.w()*q.x() + q.y()*q.z()), (1 - 2*(q.x()*q.x() + y2)));
-      }*/
+      /** Constructs and initialize Euler angles from a rotation \p rot,
+        *  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].
+        * Otherwise, the specified angle will be in the range [-PI, +PI].
+        *
+        * \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<typename Derived>
+      EulerAngles(
+        const RotationBase<Derived, 3>& rot,
+        bool positiveRangeAlpha,
+        bool positiveRangeBeta,
+        bool positiveRangeGamma) {
+        
+        System::CalcEulerAngles(*this, rot.toRotationMatrix(), positiveRangeAlpha, positiveRangeBeta, positiveRangeGamma);
+      }
 
       /** \returns The angle values stored in a vector (alpha, beta, gamma). */
       const Vector3& angles() const { return m_angles; }
@@ -227,48 +246,90 @@ namespace Eigen
         return inverse();
       }
       
-      /** 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.
+      /** Constructs and initialize Euler angles from a 3x3 rotation matrix \p m,
+        *  with options to choose for each angle the requested range (__only in compile time__).
         *
-        * See EulerAngles(const MatrixBase<Derived, 3>&) for more information about
-        *  angles ranges output.
+        * 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].
+        *
+        * \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<Derived>& m)
+      {
+        EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(Derived, 3, 3)
+        
+        EulerAngles e;
+        System::template CalcEulerAngles<
+          PositiveRangeAlpha, PositiveRangeBeta, PositiveRangeGamma, _Scalar>(e, m);
+        return e;
+      }
+      
+      /** Constructs and initialize Euler angles from a rotation \p rot,
+        *  with options to choose for each angle the requested range (__only in compile time__).
+        *
+        * 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].
+        *
+        * \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<class Derived>
-      EulerAngles& operator=(const MatrixBase<Derived>& other)
+      template<
+        bool PositiveRangeAlpha,
+        bool PositiveRangeBeta,
+        bool PositiveRangeGamma,
+        typename Derived>
+      static EulerAngles FromRotation(const RotationBase<Derived, 3>& rot)
+      {
+        return FromRotation<PositiveRangeAlpha, PositiveRangeBeta, PositiveRangeGamma>(rot.toRotationMatrix());
+      }
+      
+      /*EulerAngles& fromQuaternion(const QuaternionType& q)
       {
-        EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename Derived::Scalar>::value),
-         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
+        // TODO: Implement it in a faster way for quaternions
+        // According to http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/
+        //  we can compute only the needed matrix cells and then convert to euler angles. (see ZYX example below)
+        // Currently we compute all matrix cells from quaternion.
+
+        // Special case only for ZYX
+        //Scalar y2 = q.y() * q.y();
+        //m_angles[0] = std::atan2(2*(q.w()*q.z() + q.x()*q.y()), (1 - 2*(y2 + q.z()*q.z())));
+        //m_angles[1] = std::asin( 2*(q.w()*q.y() - q.z()*q.x()));
+        //m_angles[2] = std::atan2(2*(q.w()*q.x() + q.y()*q.z()), (1 - 2*(q.x()*q.x() + y2)));
+      }*/
+      
+      /** Set \c *this from a rotation matrix(i.e. pure orthogonal matrix with determinant of +1). */
+      template<typename Derived>
+      EulerAngles& operator=(const MatrixBase<Derived>& m) {
+        EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(Derived, 3, 3)
         
-        internal::eulerangles_assign_impl<System, Derived>::run(*this, other.derived());
+        System::CalcEulerAngles(*this, m);
         return *this;
       }
 
       // TODO: Assign and construct from another EulerAngles (with different system)
       
-      /** Set \c *this from a rotation.
-        *
-        * See EulerAngles(const RotationBase<Derived, 3>&) for more information about
-        *  angles ranges output.
-      */
+      /** Set \c *this from a rotation. */
       template<typename Derived>
       EulerAngles& operator=(const RotationBase<Derived, 3>& rot) {
         System::CalcEulerAngles(*this, rot.toRotationMatrix());
         return *this;
       }
       
-      /** \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<Scalar>::dummy_precision()) const
-      { return angles().isApprox(other.angles(), prec); }
+      // TODO: Support isApprox function
 
       /** \returns an equivalent 3x3 rotation matrix. */
       Matrix3 toRotationMatrix() const
       {
-        // TODO: Calc it faster
         return static_cast<QuaternionType>(*this).toRotationMatrix();
       }
 
@@ -286,15 +347,6 @@ namespace Eigen
         s << eulerAngles.angles().transpose();
         return s;
       }
-      
-      /** \returns \c *this with scalar type casted to \a NewScalarType */
-      template <typename NewScalarType>
-      EulerAngles<NewScalarType, System> cast() const
-      {
-        EulerAngles<NewScalarType, System> e;
-        e.angles() = angles().template cast<NewScalarType>();
-        return e;
-      }
   };
 
 #define EIGEN_EULER_ANGLES_SINGLE_TYPEDEF(AXES, SCALAR_TYPE, SCALAR_POSTFIX) \
@@ -327,29 +379,8 @@ EIGEN_EULER_ANGLES_TYPEDEFS(double, d)
     {
       typedef _Scalar Scalar;
     };
-    
-    // set from a rotation matrix
-    template<class System, class Other>
-    struct eulerangles_assign_impl<System,Other,3,3>
-    {
-      typedef typename Other::Scalar Scalar;
-      static void run(EulerAngles<Scalar, System>& e, const Other& m)
-      {
-        System::CalcEulerAngles(e, m);
-      }
-    };
-    
-    // set from a vector of Euler angles
-    template<class System, class Other>
-    struct eulerangles_assign_impl<System,Other,4,1>
-    {
-      typedef typename Other::Scalar Scalar;
-      static void run(EulerAngles<Scalar, System>& e, const Other& vec)
-      {
-        e.angles() = vec;
-      }
-    };
   }
+  
 }
 
 #endif // EIGEN_EULERANGLESCLASS_H
diff --git a/eigen/unsupported/Eigen/src/EulerAngles/EulerSystem.h b/eigen/unsupported/Eigen/src/EulerAngles/EulerSystem.h
index 28f52da..98f9f64 100644
--- a/eigen/unsupported/Eigen/src/EulerAngles/EulerSystem.h
+++ b/eigen/unsupported/Eigen/src/EulerAngles/EulerSystem.h
@@ -18,7 +18,7 @@ namespace Eigen
   
   namespace internal
   {
-    // TODO: Add this trait to the Eigen internal API?
+    // TODO: Check if already exists on the rest API
     template <int Num, bool IsPositive = (Num > 0)>
     struct Abs
     {
@@ -36,12 +36,6 @@ namespace Eigen
     {
       enum { value = Axis != 0 && Abs<Axis>::value <= 3 };
     };
-    
-    template<typename System,
-            typename Other,
-            int OtherRows=Other::RowsAtCompileTime,
-            int OtherCols=Other::ColsAtCompileTime>
-    struct eulerangles_assign_impl;
   }
   
   #define EIGEN_EULER_ANGLES_CLASS_STATIC_ASSERT(COND,MSG) typedef char static_assertion_##MSG[(COND)?1:-1]
@@ -75,7 +69,7 @@ namespace Eigen
     *
     * You can use this class to get two things:
     *  - Build an Euler system, and then pass it as a template parameter to EulerAngles.
-    *  - Query some compile time data about an Euler system. (e.g. Whether it's Tait-Bryan)
+    *  - Query some compile time data about an Euler system. (e.g. Whether it's tait bryan)
     *
     * Euler rotation is a set of three rotation on fixed axes. (see \ref EulerAngles)
     * This meta-class store constantly those signed axes. (see \ref EulerAxis)
@@ -86,7 +80,7 @@ namespace Eigen
     *  signed axes{+X,+Y,+Z,-X,-Y,-Z} are supported:
     *  - all axes X, Y, Z in each valid order (see below what order is valid)
     *  - rotation over the axis is supported both over the positive and negative directions.
-    *  - both Tait-Bryan and proper/classic Euler angles (i.e. the opposite).
+    *  - both tait bryan and proper/classic Euler angles (i.e. the opposite).
     *
     * Since EulerSystem support both positive and negative directions,
     *  you may call this rotation distinction in other names:
@@ -96,7 +90,7 @@ namespace Eigen
     * Notice all axed combination are valid, and would trigger a static assertion.
     * Same unsigned axes can't be neighbors, e.g. {X,X,Y} is invalid.
     * This yield two and only two classes:
-    *  - _Tait-Bryan_ - all unsigned axes are distinct, e.g. {X,Y,Z}
+    *  - _tait bryan_ - all unsigned axes are distinct, e.g. {X,Y,Z}
     *  - _proper/classic Euler angles_ - The first and the third unsigned axes is equal,
     *     and the second is different, e.g. {X,Y,X}
     *
@@ -118,9 +112,9 @@ namespace Eigen
     *
     * \tparam _AlphaAxis the first fixed EulerAxis
     *
-    * \tparam _BetaAxis the second fixed EulerAxis
+    * \tparam _AlphaAxis the second fixed EulerAxis
     *
-    * \tparam _GammaAxis the third fixed EulerAxis
+    * \tparam _AlphaAxis the third fixed EulerAxis
     */
   template <int _AlphaAxis, int _BetaAxis, int _GammaAxis>
   class EulerSystem
@@ -144,16 +138,14 @@ namespace Eigen
       BetaAxisAbs = internal::Abs<BetaAxis>::value, /*!< the second rotation axis unsigned */
       GammaAxisAbs = internal::Abs<GammaAxis>::value, /*!< the third rotation axis unsigned */
       
-      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 */
-
-      // 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 */
+      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 */
 
-      IsTaitBryan = ((unsigned)AlphaAxisAbs != (unsigned)GammaAxisAbs) ? 1 : 0 /*!< whether the Euler system is Tait-Bryan */
+      IsTaitBryan = ((unsigned)AlphaAxisAbs != (unsigned)GammaAxisAbs) ? 1 : 0 /*!< weather the Euler system is tait bryan */
     };
     
     private:
@@ -188,70 +180,71 @@ namespace Eigen
     static void CalcEulerAngles_imp(Matrix<typename MatrixBase<Derived>::Scalar, 3, 1>& res, const MatrixBase<Derived>& mat, internal::true_type /*isTaitBryan*/)
     {
       using std::atan2;
-      using std::sqrt;
+      using std::sin;
+      using std::cos;
       
       typedef typename Derived::Scalar Scalar;
-
-      const Scalar plusMinus = IsEven? 1 : -1;
-      const Scalar minusPlus = IsOdd?  1 : -1;
-
-      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<Scalar>::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) {// 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 {// 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);
-        res[0] = alphaMinusPlusBeta;
-        res[2] = 0;
+      typedef Matrix<Scalar,2,1> 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)) || ((!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);
       }
+      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 <typename Derived>
-    static void CalcEulerAngles_imp(Matrix<typename MatrixBase<Derived>::Scalar,3,1>& res,
-                                    const MatrixBase<Derived>& mat, internal::false_type /*isTaitBryan*/)
+    static void CalcEulerAngles_imp(Matrix<typename MatrixBase<Derived>::Scalar,3,1>& res, const MatrixBase<Derived>& mat, internal::false_type /*isTaitBryan*/)
     {
       using std::atan2;
-      using std::sqrt;
+      using std::sin;
+      using std::cos;
 
       typedef typename Derived::Scalar Scalar;
-
-      const Scalar plusMinus = IsEven? 1 : -1;
-      const Scalar minusPlus = IsOdd?  1 : -1;
-
-      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));
-
-      // There is a singularity when sin(beta) == 0
-      if(Rsum > 4 * NumTraits<Scalar>::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) {// 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;
+      typedef Matrix<Scalar,2,1> Vector2;
+      
+      res[0] = atan2(mat(J,I), mat(K,I));
+      if((IsOdd && res[0]<Scalar(0)) || ((!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));
       }
-      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[2] = 0;
+      else
+      {
+        Scalar s2 = Vector2(mat(J,I), mat(K,I)).norm();
+        res[1] = atan2(s2, mat(I,I));
       }
+
+      // 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<typename Scalar>
@@ -259,28 +252,55 @@ namespace Eigen
       EulerAngles<Scalar, EulerSystem>& res,
       const typename EulerAngles<Scalar, EulerSystem>::Matrix3& mat)
     {
+      CalcEulerAngles(res, mat, false, false, false);
+    }
+    
+    template<
+      bool PositiveRangeAlpha,
+      bool PositiveRangeBeta,
+      bool PositiveRangeGamma,
+      typename Scalar>
+    static void CalcEulerAngles(
+      EulerAngles<Scalar, EulerSystem>& res,
+      const typename EulerAngles<Scalar, EulerSystem>::Matrix3& mat)
+    {
+      CalcEulerAngles(res, mat, PositiveRangeAlpha, PositiveRangeBeta, PositiveRangeGamma);
+    }
+    
+    template<typename Scalar>
+    static void CalcEulerAngles(
+      EulerAngles<Scalar, EulerSystem>& res,
+      const typename EulerAngles<Scalar, EulerSystem>::Matrix3& mat,
+      bool PositiveRangeAlpha,
+      bool PositiveRangeBeta,
+      bool PositiveRangeGamma)
+    {
       CalcEulerAngles_imp(
         res.angles(), mat,
         typename internal::conditional<IsTaitBryan, internal::true_type, internal::false_type>::type());
 
-      if (IsAlphaOpposite)
+      if (IsAlphaOpposite == IsOdd)
         res.alpha() = -res.alpha();
         
-      if (IsBetaOpposite)
+      if (IsBetaOpposite == IsOdd)
         res.beta() = -res.beta();
         
-      if (IsGammaOpposite)
+      if (IsGammaOpposite == IsOdd)
         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);
     }
     
     template <typename _Scalar, class _System>
     friend class Eigen::EulerAngles;
-    
-    template<typename System,
-            typename Other,
-            int OtherRows,
-            int OtherCols>
-    friend struct internal::eulerangles_assign_impl;
   };
 
 #define EIGEN_EULER_SYSTEM_TYPEDEF(A, B, C) \
diff --git a/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h
index db2449d..3f7d777 100644
--- a/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h
+++ b/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h
@@ -398,8 +398,8 @@ struct matrix_function_compute
 template <typename MatrixType>
 struct matrix_function_compute<MatrixType, 0>
 {  
-  template <typename AtomicType, typename ResultType> 
-  static void run(const MatrixType& A, AtomicType& atomic, ResultType &result)
+  template <typename MatA, typename AtomicType, typename ResultType>
+  static void run(const MatA& A, AtomicType& atomic, ResultType &result)
   {
     typedef internal::traits<MatrixType> Traits;
     typedef typename Traits::Scalar Scalar;
@@ -422,11 +422,10 @@ struct matrix_function_compute<MatrixType, 0>
 template <typename MatrixType>
 struct matrix_function_compute<MatrixType, 1>
 {
-  template <typename AtomicType, typename ResultType> 
-  static void run(const MatrixType& A, AtomicType& atomic, ResultType &result)
+  template <typename MatA, typename AtomicType, typename ResultType>
+  static void run(const MatA& A, AtomicType& atomic, ResultType &result)
   {
     typedef internal::traits<MatrixType> Traits;
-    typedef typename MatrixType::Index Index;
     
     // compute Schur decomposition of A
     const ComplexSchur<MatrixType> schurOfA(A);  
@@ -514,7 +513,7 @@ template<typename Derived> class MatrixFunctionReturnValue
       typedef internal::MatrixFunctionAtomic<DynMatrixType> AtomicType;
       AtomicType atomic(m_f);
 
-      internal::matrix_function_compute<NestedEvalTypeClean>::run(m_A, atomic, result);
+      internal::matrix_function_compute<typename NestedEvalTypeClean::PlainObject>::run(m_A, atomic, result);
     }
 
     Index rows() const { return m_A.rows(); }
diff --git a/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
index 1acfbed..ff8f6e7 100644
--- a/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
+++ b/eigen/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
@@ -339,7 +339,7 @@ public:
     typedef internal::MatrixLogarithmAtomic<DynMatrixType> AtomicType;
     AtomicType atomic;
     
-    internal::matrix_function_compute<DerivedEvalTypeClean>::run(m_A, atomic, result);
+    internal::matrix_function_compute<typename DerivedEvalTypeClean::PlainObject>::run(m_A, atomic, result);
   }
 
   Index rows() const { return m_A.rows(); }
diff --git a/eigen/unsupported/Eigen/src/Polynomials/Companion.h b/eigen/unsupported/Eigen/src/Polynomials/Companion.h
index e0af6eb..b515c29 100644
--- a/eigen/unsupported/Eigen/src/Polynomials/Companion.h
+++ b/eigen/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 = Scalar(-1)/poly[deg] * poly.head(deg);
+      m_monic = -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( RealScalar colNorm, RealScalar rowNorm,
-        bool& isBalanced, RealScalar& colB, RealScalar& rowB );
+    bool balanced( Scalar colNorm, Scalar rowNorm,
+        bool& isBalanced, Scalar& colB, Scalar& 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( RealScalar colNorm, RealScalar rowNorm,
-        bool& isBalanced, RealScalar& colB, RealScalar& rowB );
+    bool balancedR( Scalar colNorm, Scalar rowNorm,
+        bool& isBalanced, Scalar& colB, Scalar& rowB );
 
   public:
     /**
@@ -139,10 +139,10 @@ class companion
 
 template< typename _Scalar, int _Deg >
 inline
-bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm,
-    bool& isBalanced, RealScalar& colB, RealScalar& rowB )
+bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm,
+    bool& isBalanced, Scalar& colB, Scalar& rowB )
 {
-  if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm ){ return true; }
+  if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; }
   else
   {
     //To find the balancing coefficients, if the radix is 2,
@@ -150,29 +150,29 @@ bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar 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<RealScalar>();
-    colB = RealScalar(1);
-    const RealScalar s = colNorm + rowNorm;
+    rowB = rowNorm / radix<Scalar>();
+    colB = Scalar(1);
+    const Scalar s = colNorm + rowNorm;
 
     while (colNorm < rowB)
     {
-      colB *= radix<RealScalar>();
-      colNorm *= radix2<RealScalar>();
+      colB *= radix<Scalar>();
+      colNorm *= radix2<Scalar>();
     }
 
-    rowB = rowNorm * radix<RealScalar>();
+    rowB = rowNorm * radix<Scalar>();
 
     while (colNorm >= rowB)
     {
-      colB /= radix<RealScalar>();
-      colNorm /= radix2<RealScalar>();
+      colB /= radix<Scalar>();
+      colNorm /= radix2<Scalar>();
     }
 
     //This line is used to avoid insubstantial balancing
-    if ((rowNorm + colNorm) < RealScalar(0.95) * s * colB)
+    if ((rowNorm + colNorm) < Scalar(0.95) * s * colB)
     {
       isBalanced = false;
-      rowB = RealScalar(1) / colB;
+      rowB = Scalar(1) / colB;
       return false;
     }
     else{
@@ -182,21 +182,21 @@ bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm,
 
 template< typename _Scalar, int _Deg >
 inline
-bool companion<_Scalar,_Deg>::balancedR( RealScalar colNorm, RealScalar rowNorm,
-    bool& isBalanced, RealScalar& colB, RealScalar& rowB )
+bool companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm,
+    bool& isBalanced, Scalar& colB, Scalar& rowB )
 {
-  if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm ){ return true; }
+  if( Scalar(0) == colNorm || Scalar(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 RealScalar q = colNorm/rowNorm;
+    const _Scalar q = colNorm/rowNorm;
     if( !isApprox( q, _Scalar(1) ) )
     {
       rowB = sqrt( colNorm/rowNorm );
-      colB = RealScalar(1)/rowB;
+      colB = Scalar(1)/rowB;
 
       isBalanced = false;
       return false;
@@ -219,8 +219,8 @@ void companion<_Scalar,_Deg>::balance()
   while( !hasConverged )
   {
     hasConverged = true;
-    RealScalar colNorm,rowNorm;
-    RealScalar colB,rowB;
+    Scalar colNorm,rowNorm;
+    Scalar colB,rowB;
 
     //First row, first column excluding the diagonal
     //==============================================
diff --git a/eigen/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/eigen/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
index 7885942..03198ec 100644
--- a/eigen/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
+++ b/eigen/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
@@ -99,7 +99,7 @@ class PolynomialSolverBase
      */
     inline const RootType& greatestRoot() const
     {
-      std::greater<RealScalar> greater;
+      std::greater<Scalar> greater;
       return selectComplexRoot_withRespectToNorm( greater );
     }
 
@@ -108,7 +108,7 @@ class PolynomialSolverBase
      */
     inline const RootType& smallestRoot() const
     {
-      std::less<RealScalar> less;
+      std::less<Scalar> less;
       return selectComplexRoot_withRespectToNorm( less );
     }
 
@@ -213,7 +213,7 @@ class PolynomialSolverBase
         bool& hasArealRoot,
         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
     {
-      std::greater<RealScalar> greater;
+      std::greater<Scalar> greater;
       return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
     }
 
@@ -236,7 +236,7 @@ class PolynomialSolverBase
         bool& hasArealRoot,
         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
     {
-      std::less<RealScalar> less;
+      std::less<Scalar> less;
       return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
     }
 
@@ -259,7 +259,7 @@ class PolynomialSolverBase
         bool& hasArealRoot,
         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
     {
-      std::greater<RealScalar> greater;
+      std::greater<Scalar> greater;
       return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
     }
 
@@ -282,7 +282,7 @@ class PolynomialSolverBase
         bool& hasArealRoot,
         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
     {
-      std::less<RealScalar> less;
+      std::less<Scalar> 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>
+template< typename _Scalar, int _Deg >
 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
 {
   public:
@@ -337,9 +337,7 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
     EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
 
     typedef Matrix<Scalar,_Deg,_Deg>                 CompanionMatrixType;
-    typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
-                                          ComplexEigenSolver<CompanionMatrixType>,
-                                          EigenSolver<CompanionMatrixType> >::type EigenSolverType;
+    typedef EigenSolver<CompanionMatrixType>         EigenSolverType;
 
   public:
     /** Computes the complex roots of a new polynomial. */
diff --git a/eigen/unsupported/Eigen/src/SparseExtra/MarketIO.h b/eigen/unsupported/Eigen/src/SparseExtra/MarketIO.h
index fc70a24..cdc14f8 100644
--- a/eigen/unsupported/Eigen/src/SparseExtra/MarketIO.h
+++ b/eigen/unsupported/Eigen/src/SparseExtra/MarketIO.h
@@ -12,38 +12,38 @@
 #define EIGEN_SPARSE_MARKET_IO_H
 
 #include <iostream>
-#include <vector>
 
 namespace Eigen { 
 
 namespace internal 
 {
-  template <typename Scalar, typename StorageIndex>
-  inline void GetMarketLine (const char* line, StorageIndex& i, StorageIndex& j, Scalar& value)
+  template <typename Scalar>
+  inline bool GetMarketLine (std::stringstream& line, Index& M, Index& N, Index& i, Index& j, Scalar& value)
   {
-    std::stringstream sline(line);
-    sline >> i >> j >> value;
+    line >> i >> j >> value;
+    i--;
+    j--;
+    if(i>=0 && j>=0 && i<M && j<N)
+    {
+      return true; 
+    }
+    else
+      return false;
   }
-
-  template<> inline void GetMarketLine (const char* line, int& i, int& j, float& value)
-  { std::sscanf(line, "%d %d %g", &i, &j, &value); }
-
-  template<> inline void GetMarketLine (const char* line, int& i, int& j, double& value)
-  { std::sscanf(line, "%d %d %lg", &i, &j, &value); }
-
-  template<> inline void GetMarketLine (const char* line, int& i, int& j, std::complex<float>& value)
-  { std::sscanf(line, "%d %d %g %g", &i, &j, &numext::real_ref(value), &numext::imag_ref(value)); }
-
-  template<> inline void GetMarketLine (const char* line, int& i, int& j, std::complex<double>& value)
-  { std::sscanf(line, "%d %d %lg %lg", &i, &j, &numext::real_ref(value), &numext::imag_ref(value)); }
-
-  template <typename Scalar, typename StorageIndex>
-  inline void GetMarketLine (const char* line, StorageIndex& i, StorageIndex& j, std::complex<Scalar>& value)
+  template <typename Scalar>
+  inline bool GetMarketLine (std::stringstream& line, Index& M, Index& N, Index& i, Index& j, std::complex<Scalar>& value)
   {
-    std::stringstream sline(line);
     Scalar valR, valI;
-    sline >> i >> j >> valR >> valI;
-    value = std::complex<Scalar>(valR,valI);
+    line >> i >> j >> valR >> valI;
+    i--;
+    j--;
+    if(i>=0 && j>=0 && i<M && j<N)
+    {
+      value = std::complex<Scalar>(valR, valI);
+      return true; 
+    }
+    else
+      return false;
   }
 
   template <typename RealScalar>
@@ -81,13 +81,13 @@ namespace internal
     }
   }
 
-  template<typename Scalar, typename StorageIndex>
-  inline void PutMatrixElt(Scalar value, StorageIndex row, StorageIndex col, std::ofstream& out)
+  template<typename Scalar>
+  inline void PutMatrixElt(Scalar value, int row, int col, std::ofstream& out)
   {
     out << row << " "<< col << " " << value << "\n";
   }
-  template<typename Scalar, typename StorageIndex>
-  inline void PutMatrixElt(std::complex<Scalar> value, StorageIndex row, StorageIndex col, std::ofstream& out)
+  template<typename Scalar>
+  inline void PutMatrixElt(std::complex<Scalar> value, int row, int col, std::ofstream& out)
   {
     out << row << " " << col << " " << value.real() << " " << value.imag() << "\n";
   }
@@ -133,20 +133,17 @@ template<typename SparseMatrixType>
 bool loadMarket(SparseMatrixType& mat, const std::string& filename)
 {
   typedef typename SparseMatrixType::Scalar Scalar;
-  typedef typename SparseMatrixType::StorageIndex StorageIndex;
+  typedef typename SparseMatrixType::Index Index;
   std::ifstream input(filename.c_str(),std::ios::in);
   if(!input)
     return false;
-
-  char rdbuffer[4096];
-  input.rdbuf()->pubsetbuf(rdbuffer, 4096);
   
   const int maxBuffersize = 2048;
   char buffer[maxBuffersize];
   
   bool readsizes = false;
 
-  typedef Triplet<Scalar,StorageIndex> T;
+  typedef Triplet<Scalar,Index> T;
   std::vector<T> elements;
   
   Index M(-1), N(-1), NNZ(-1);
@@ -157,36 +154,33 @@ bool loadMarket(SparseMatrixType& mat, const std::string& filename)
     //NOTE An appropriate test should be done on the header to get the  symmetry
     if(buffer[0]=='%')
       continue;
-
+    
+    std::stringstream line(buffer);
+    
     if(!readsizes)
     {
-      std::stringstream line(buffer);
       line >> M >> N >> NNZ;
       if(M > 0 && N > 0 && NNZ > 0) 
       {
         readsizes = true;
+        //std::cout << "sizes: " << M << "," << N << "," << NNZ << "\n";
         mat.resize(M,N);
         mat.reserve(NNZ);
       }
     }
     else
     { 
-      StorageIndex i(-1), j(-1);
+      Index i(-1), j(-1);
       Scalar value; 
-      internal::GetMarketLine(buffer, i, j, value);
-
-      i--;
-      j--;
-      if(i>=0 && j>=0 && i<M && j<N)
+      if( internal::GetMarketLine(line, M, N, i, j, value) ) 
       {
-        ++count;
+        ++ count;
         elements.push_back(T(i,j,value));
       }
-      else
+      else 
         std::cerr << "Invalid read: " << i << "," << j << "\n";        
     }
   }
-
   mat.setFromTriplets(elements.begin(), elements.end());
   if(count!=NNZ)
     std::cerr << count << "!=" << NNZ << "\n";
@@ -231,13 +225,12 @@ template<typename SparseMatrixType>
 bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0)
 {
   typedef typename SparseMatrixType::Scalar Scalar;
-  typedef typename SparseMatrixType::RealScalar RealScalar;
   std::ofstream out(filename.c_str(),std::ios::out);
   if(!out)
     return false;
   
   out.flags(std::ios_base::scientific);
-  out.precision(std::numeric_limits<RealScalar>::digits10 + 2);
+  out.precision(64);
   std::string header; 
   internal::putMarketHeader<Scalar>(header, sym); 
   out << header << std::endl; 
@@ -248,6 +241,7 @@ bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sy
     {
       ++ count;
       internal::PutMatrixElt(it.value(), it.row()+1, it.col()+1, out);
+      // out << it.row()+1 << " " << it.col()+1 << " " << it.value() << "\n";
     }
   out.close();
   return true;
@@ -256,14 +250,13 @@ bool saveMarket(const SparseMatrixType& mat, const std::string& filename, int sy
 template<typename VectorType>
 bool saveMarketVector (const VectorType& vec, const std::string& filename)
 {
- typedef typename VectorType::Scalar Scalar;
- typedef typename VectorType::RealScalar RealScalar;
+ typedef typename VectorType::Scalar Scalar; 
  std::ofstream out(filename.c_str(),std::ios::out);
   if(!out)
     return false;
   
   out.flags(std::ios_base::scientific);
-  out.precision(std::numeric_limits<RealScalar>::digits10 + 2);
+  out.precision(64);
   if(internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value)
       out << "%%MatrixMarket matrix array complex general\n"; 
   else
diff --git a/eigen/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/eigen/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
index 369ad97..f524d71 100644
--- a/eigen/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
+++ b/eigen/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
@@ -122,8 +122,8 @@ struct lgamma_impl<float> {
   EIGEN_DEVICE_FUNC
   static EIGEN_STRONG_INLINE float run(float x) {
 #if !defined(__CUDA_ARCH__) && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) && !defined(__APPLE__)
-    int dummy;
-    return ::lgammaf_r(x, &dummy);
+    int signgam;
+    return ::lgammaf_r(x, &signgam);
 #else
     return ::lgammaf(x);
 #endif
@@ -135,8 +135,8 @@ struct lgamma_impl<double> {
   EIGEN_DEVICE_FUNC
   static EIGEN_STRONG_INLINE double run(double x) {
 #if !defined(__CUDA_ARCH__) && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) && !defined(__APPLE__)
-    int dummy;
-    return ::lgamma_r(x, &dummy);
+    int signgam;
+    return ::lgamma_r(x, &signgam);
 #else
     return ::lgamma(x);
 #endif