* merge with mainline

* adapt Eigenvalues module to the new rule that the RowMajorBit must have the proper value for vectors
* Fix RowMajorBit in ei_traits<ProductBase>
* Fix vectorizability logic in CoeffBasedProduct

62 files changed, 3204 insertions, 804 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index a85bbf222..7c68b0f71 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -152,11 +152,24 @@ option(EIGEN_TEST_C++0x "Enables all C++0x features." OFF)
 
 include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR})
 
-set(INCLUDE_INSTALL_DIR
+# the user modifiable install path for header files
+set(EIGEN_INCLUDE_INSTALL_DIR ${EIGEN_INCLUDE_INSTALL_DIR} CACHE PATH "The directory where we install the header files (optional)")
+
+# set the internal install path for header files which depends on wether the user modifiable
+# EIGEN_INCLUDE_INSTALL_DIR has been set by the user or not.
+if(EIGEN_INCLUDE_INSTALL_DIR)
+  set(INCLUDE_INSTALL_DIR
+    ${EIGEN_INCLUDE_INSTALL_DIR}
+    CACHE INTERNAL
+    "The directory where we install the header files (internal)"
+  )
+else()
+  set(INCLUDE_INSTALL_DIR
     "${CMAKE_INSTALL_PREFIX}/include/eigen3"
-    CACHE PATH
-    "The directory where we install the header files"
-)
+    CACHE INTERNAL
+    "The directory where we install the header files (internal)"
+  )
+endif()
 
 install(FILES
   signature_of_eigen3_matrix_library
@@ -213,6 +226,8 @@ if(cmake_generator_tolower MATCHES "makefile")
   message("--------------+----------------------------------------------------------------")
   message("make install  | Install to ${CMAKE_INSTALL_PREFIX}")
   message("              |   To change that: cmake . -DCMAKE_INSTALL_PREFIX=yourpath")
+  message("              |   Header files are installed to ${INCLUDE_INSTALL_DIR}")
+  message("              |   To change that: cmake . -DEIGEN_INCLUDE_INSTALL_DIR=yourpath")
   message("make doc      | Generate the API documentation, requires Doxygen & LaTeX")
   message("make check    | Build and run the unit-tests. Read this page:")
   message("              |   http://eigen.tuxfamily.org/index.php?title=Tests")
diff --git a/Eigen/Eigenvalues b/Eigen/Eigenvalues
index 1ae0cf098..986f31196 100644
--- a/Eigen/Eigenvalues
+++ b/Eigen/Eigenvalues
@@ -36,6 +36,7 @@ namespace Eigen {
   */
 
 #include "src/Eigenvalues/Tridiagonalization.h"
+#include "src/Eigenvalues/RealSchur.h"
 #include "src/Eigenvalues/EigenSolver.h"
 #include "src/Eigenvalues/SelfAdjointEigenSolver.h"
 #include "src/Eigenvalues/HessenbergDecomposition.h"
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 51a0e44ae..8b01e0ccb 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -295,8 +295,7 @@ template<typename Derived>
 bool LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
 {
   ei_assert(m_isInitialized && "LLT is not initialized.");
-  const int size = m_matrix.rows();
-  ei_assert(size==bAndX.rows());
+  ei_assert(m_matrix.rows()==bAndX.rows());
   matrixL().solveInPlace(bAndX);
   matrixU().solveInPlace(bAndX);
   return true;
diff --git a/Eigen/src/Core/DenseBase.h b/Eigen/src/Core/DenseBase.h
index 6459cd1b1..566b4b410 100644
--- a/Eigen/src/Core/DenseBase.h
+++ b/Eigen/src/Core/DenseBase.h
@@ -528,7 +528,7 @@ template<typename Derived> class DenseBase
     #endif
 
     // disable the use of evalTo for dense objects with a nice compilation error
-    template<typename Dest> inline void evalTo(Dest& dst) const
+    template<typename Dest> inline void evalTo(Dest& ) const
     {
       EIGEN_STATIC_ASSERT((ei_is_same_type<Dest,void>::ret),THE_EVAL_EVALTO_FUNCTION_SHOULD_NEVER_BE_CALLED_FOR_DENSE_OBJECTS);
     }
diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h
index c2b317cc0..d02633cb8 100644
--- a/Eigen/src/Core/Functors.h
+++ b/Eigen/src/Core/Functors.h
@@ -274,7 +274,7 @@ template<typename Scalar, typename NewType>
 struct ei_scalar_cast_op {
   EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_cast_op)
   typedef NewType result_type;
-  EIGEN_STRONG_INLINE const NewType operator() (const Scalar& a) const { return static_cast<NewType>(a); }
+  EIGEN_STRONG_INLINE const NewType operator() (const Scalar& a) const { return ei_cast<Scalar, NewType>(a); }
 };
 template<typename Scalar, typename NewType>
 struct ei_functor_traits<ei_scalar_cast_op<Scalar,NewType> >
diff --git a/Eigen/src/Core/IO.h b/Eigen/src/Core/IO.h
index 3e8d2bc66..c98742246 100644
--- a/Eigen/src/Core/IO.h
+++ b/Eigen/src/Core/IO.h
@@ -126,6 +126,16 @@ DenseBase<Derived>::format(const IOFormat& fmt) const
   return WithFormat<Derived>(derived(), fmt);
 }
 
+template<typename Scalar>
+struct ei_significant_decimals_impl
+{
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  static inline int run()
+  {
+    return ei_cast<RealScalar,int>(std::ceil(-ei_log(NumTraits<RealScalar>::epsilon())/ei_log(RealScalar(10))));
+  }
+};
+
 /** \internal
   * print the matrix \a _m to the output stream \a s using the output format \a fmt */
 template<typename Derived>
@@ -145,9 +155,7 @@ std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOForm
   {
     if (NumTraits<Scalar>::HasFloatingPoint)
     {
-      typedef typename NumTraits<Scalar>::Real RealScalar;
-      RealScalar explicit_precision_fp = std::ceil(-ei_log(NumTraits<Scalar>::epsilon())/ei_log(10.0));
-      explicit_precision = static_cast<std::streamsize>(explicit_precision_fp);
+      explicit_precision = ei_significant_decimals_impl<Scalar>::run();
     }
     else
     {
diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index c97a68e50..4a21ec975 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -44,6 +44,23 @@ template<typename T> inline typename NumTraits<T>::Real ei_hypot(T x, T y)
   return p * ei_sqrt(T(1) + qp*qp);
 }
 
+// the point of wrapping these casts in this helper template struct is to allow users to specialize it to custom types
+// that may not have the needed conversion operators (especially as c++98 doesn't have explicit conversion operators).
+
+template<typename OldType, typename NewType> struct ei_cast_impl
+{
+  static inline NewType run(const OldType& x)
+  {
+    return static_cast<NewType>(x);
+  }
+};
+
+template<typename OldType, typename NewType> inline NewType ei_cast(const OldType& x)
+{
+  return ei_cast_impl<OldType, NewType>::run(x);
+}
+
+
 /**************
 ***   int   ***
 **************/
diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h
index 4013f6ab1..b7c4ac11d 100644
--- a/Eigen/src/Core/ProductBase.h
+++ b/Eigen/src/Core/ProductBase.h
@@ -42,7 +42,7 @@ struct ei_traits<ProductBase<Derived,_Lhs,_Rhs> > //: ei_traits<typename ei_clea
     ColsAtCompileTime = ei_traits<Rhs>::ColsAtCompileTime,
     MaxRowsAtCompileTime = ei_traits<Lhs>::MaxRowsAtCompileTime,
     MaxColsAtCompileTime = ei_traits<Rhs>::MaxColsAtCompileTime,
-    Flags = (RowsAtCompileTime==1 ? RowMajorBit : 0)
+    Flags = (MaxRowsAtCompileTime==1 ? RowMajorBit : 0)
           | EvalBeforeNestingBit | EvalBeforeAssigningBit | NestByRefBit,
                   // Note that EvalBeforeNestingBit and NestByRefBit
                   // are not used in practice because ei_nested is overloaded for products
diff --git a/Eigen/src/Core/arch/SSE/MathFunctions.h b/Eigen/src/Core/arch/SSE/MathFunctions.h
index 3c0020248..99662eb6d 100644
--- a/Eigen/src/Core/arch/SSE/MathFunctions.h
+++ b/Eigen/src/Core/arch/SSE/MathFunctions.h
@@ -369,10 +369,14 @@ static EIGEN_DONT_INLINE EIGEN_UNUSED Packet4f ei_pcos(Packet4f x)
 // For detail see here: http://www.beyond3d.com/content/articles/8/
 static EIGEN_UNUSED Packet4f ei_psqrt(Packet4f _x)
 {
-  Packet4f half = ei_pmul(_x, ei_pset1(.5f));
-  Packet4f x = _mm_rsqrt_ps(_x);
-  x = ei_pmul(x, ei_psub(ei_pset1(1.5f), ei_pmul(half, ei_pmul(x,x))));
-  return ei_pmul(_x,x);
+	Packet4f half = ei_pmul(_x, ei_pset1(.5f));
+	
+	/* select only the inverse sqrt of non-zero inputs */
+	Packet4f non_zero_mask = _mm_cmpgt_ps(_x, ei_pset1(std::numeric_limits<float>::epsilon()));
+	Packet4f x = _mm_and_ps(non_zero_mask, _mm_rsqrt_ps(_x));
+
+	x = ei_pmul(x, ei_psub(ei_pset1(1.5f), ei_pmul(half, ei_pmul(x,x))));
+	return ei_pmul(_x,x);
 }
 
 #endif // EIGEN_MATH_FUNCTIONS_SSE_H
diff --git a/Eigen/src/Core/products/CoeffBasedProduct.h b/Eigen/src/Core/products/CoeffBasedProduct.h
index 17fbc9190..2f7b32c65 100644
--- a/Eigen/src/Core/products/CoeffBasedProduct.h
+++ b/Eigen/src/Core/products/CoeffBasedProduct.h
@@ -72,10 +72,18 @@ struct ei_traits<CoeffBasedProduct<LhsNested,RhsNested,NestingFlags> >
       RhsRowMajor = RhsFlags & RowMajorBit,
 
       CanVectorizeRhs = RhsRowMajor && (RhsFlags & PacketAccessBit)
-                      && (ColsAtCompileTime == Dynamic || (ColsAtCompileTime % ei_packet_traits<Scalar>::size) == 0),
+                      && (ColsAtCompileTime == Dynamic
+                          || ( (ColsAtCompileTime % ei_packet_traits<Scalar>::size) == 0
+                              && (RhsFlags&AlignedBit)
+                             )
+                         ),
 
       CanVectorizeLhs = (!LhsRowMajor) && (LhsFlags & PacketAccessBit)
-                      && (RowsAtCompileTime == Dynamic || (RowsAtCompileTime % ei_packet_traits<Scalar>::size) == 0),
+                      && (RowsAtCompileTime == Dynamic
+                          || ( (RowsAtCompileTime % ei_packet_traits<Scalar>::size) == 0
+                              && (LhsFlags&AlignedBit)
+                             )
+                         ),
 
       EvalToRowMajor = (MaxRowsAtCompileTime==1&&MaxColsAtCompileTime!=1) ? 1
                      : (MaxColsAtCompileTime==1&&MaxRowsAtCompileTime!=1) ? 0
@@ -84,8 +92,7 @@ struct ei_traits<CoeffBasedProduct<LhsNested,RhsNested,NestingFlags> >
       Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & ~RowMajorBit)
             | (EvalToRowMajor ? RowMajorBit : 0)
             | NestingFlags
-            | (CanVectorizeLhs || CanVectorizeRhs ? PacketAccessBit : 0)
-            | (LhsFlags & RhsFlags & AlignedBit),
+            | (CanVectorizeLhs || CanVectorizeRhs ? PacketAccessBit : 0),
 
       CoeffReadCost = InnerSize == Dynamic ? Dynamic
                     : InnerSize * (NumTraits<Scalar>::MulCost + LhsCoeffReadCost + RhsCoeffReadCost)
@@ -96,8 +103,11 @@ struct ei_traits<CoeffBasedProduct<LhsNested,RhsNested,NestingFlags> >
       * loop of the product might be vectorized. This is the meaning of CanVectorizeInner. Since it doesn't affect
       * the Flags, it is safe to make this value depend on ActualPacketAccessBit, that doesn't affect the ABI.
       */
-      CanVectorizeInner = LhsRowMajor && (!RhsRowMajor) && (LhsFlags & RhsFlags & ActualPacketAccessBit)
-                        && (InnerSize % ei_packet_traits<Scalar>::size == 0)
+      CanVectorizeInner =    LhsRowMajor
+                          && (!RhsRowMajor)
+                          && (LhsFlags & RhsFlags & ActualPacketAccessBit)
+                          && (LhsFlags & RhsFlags & AlignedBit)
+                          && (InnerSize % ei_packet_traits<Scalar>::size == 0)
     };
 };
 
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index fe320c3a0..0b2e0b0ec 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -31,9 +31,24 @@
   *
   * \class ComplexEigenSolver
   *
-  * \brief Eigen values/vectors solver for general complex matrices
+  * \brief Computes eigenvalues and eigenvectors of general complex matrices
   *
-  * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
+  * \tparam _MatrixType the type of the matrix of which we are
+  * computing the eigendecomposition; this is expected to be an
+  * instantiation of the Matrix class template.
+  *
+  * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
+  * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v
+  * \f$.  If \f$ D \f$ is a diagonal matrix with the eigenvalues on
+  * the diagonal, and \f$ V \f$ is a matrix with the eigenvectors as
+  * its columns, then \f$ A V = V D \f$. The matrix \f$ V \f$ is
+  * almost always invertible, in which case we have \f$ A = V D V^{-1}
+  * \f$. This is called the eigendecomposition.
+  *
+  * The main function in this class is compute(), which computes the
+  * eigenvalues and eigenvectors of a given function. The
+  * documentation for that function contains an example showing the
+  * main features of the class.
   *
   * \sa class EigenSolver, class SelfAdjointEigenSolver
   */
@@ -48,21 +63,47 @@ template<typename _MatrixType> class ComplexEigenSolver
       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };
+
+    /** \brief Scalar type for matrices of type \p _MatrixType. */
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<Scalar>::Real RealScalar;
-    typedef std::complex<RealScalar> Complex;
-    typedef typename ei_plain_col_type<MatrixType, Complex>::type EigenvalueType;
-    typedef Matrix<Complex, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, ColsAtCompileTime> EigenvectorType;
-
-    /**
-    * \brief Default Constructor.
-    *
-    * The default constructor is useful in cases in which the user intends to
-    * perform decompositions via ComplexEigenSolver::compute(const MatrixType&).
-    */
+
+    /** \brief Complex scalar type for \p _MatrixType. 
+      *
+      * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
+      * \c float or \c double) and just \c Scalar if #Scalar is
+      * complex.
+      */
+    typedef std::complex<RealScalar> ComplexScalar;
+
+    /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 
+      *
+      * This is a column vector with entries of type #ComplexScalar.
+      * The length of the vector is the size of \p _MatrixType.
+      */
+    typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options, MaxColsAtCompileTime, 1> EigenvalueType;
+
+    /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 
+      *
+      * This is a square matrix with entries of type #ComplexScalar. 
+      * The size is the same as the size of \p _MatrixType.
+      */
+    typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, ColsAtCompileTime> EigenvectorType;
+
+    /** \brief Default constructor.
+      *
+      * The default constructor is useful in cases in which the user intends to
+      * perform decompositions via compute().
+      */
     ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false)
     {}
 
+    /** \brief Constructor; computes eigendecomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
+      *
+      * This constructor calls compute() to compute the eigendecomposition.
+      */
     ComplexEigenSolver(const MatrixType& matrix)
             : m_eivec(matrix.rows(),matrix.cols()),
               m_eivalues(matrix.cols()),
@@ -71,22 +112,72 @@ template<typename _MatrixType> class ComplexEigenSolver
       compute(matrix);
     }
 
-    EigenvectorType eigenvectors(void) const
+    /** \brief Returns the eigenvectors of given matrix. 
+      *
+      * It is assumed that either the constructor
+      * ComplexEigenSolver(const MatrixType& matrix) or the member
+      * function compute(const MatrixType& matrix) has been called
+      * before to compute the eigendecomposition of a matrix. This
+      * function returns a matrix whose columns are the
+      * eigenvectors. Column \f$ k \f$ is an eigenvector
+      * corresponding to eigenvalue number \f$ k \f$ as returned by
+      * eigenvalues().  The eigenvectors are normalized to have
+      * (Euclidean) norm equal to one. The matrix returned by this
+      * function is the matrix \f$ V \f$ in the eigendecomposition \f$
+      * A = V D V^{-1} \f$, if it exists.
+      *
+      * Example: \include ComplexEigenSolver_eigenvectors.cpp
+      * Output: \verbinclude ComplexEigenSolver_eigenvectors.out
+      */
+    EigenvectorType eigenvectors() const
     {
       ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
       return m_eivec;
     }
 
+    /** \brief Returns the eigenvalues of given matrix. 
+      *
+      * It is assumed that either the constructor
+      * ComplexEigenSolver(const MatrixType& matrix) or the member
+      * function compute(const MatrixType& matrix) has been called
+      * before to compute the eigendecomposition of a matrix. This
+      * function returns a column vector containing the
+      * eigenvalues. Eigenvalues are repeated according to their
+      * algebraic multiplicity, so there are as many eigenvalues as
+      * rows in the matrix.
+      *
+      * Example: \include ComplexEigenSolver_eigenvalues.cpp
+      * Output: \verbinclude ComplexEigenSolver_eigenvalues.out
+      */
     EigenvalueType eigenvalues() const
     {
       ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
       return m_eivalues;
     }
 
+    /** \brief Computes eigendecomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
+      *
+      * This function computes the eigenvalues and eigenvectors of \p
+      * matrix.  The eigenvalues() and eigenvectors() functions can be
+      * used to retrieve the computed eigendecomposition.
+      *
+      * The matrix is first reduced to Schur form using the
+      * ComplexSchur class. The Schur decomposition is then used to
+      * compute the eigenvalues and eigenvectors.
+      *
+      * The cost of the computation is dominated by the cost of the
+      * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$
+      * is the size of the matrix.
+      *
+      * Example: \include ComplexEigenSolver_compute.cpp
+      * Output: \verbinclude ComplexEigenSolver_compute.out
+      */
     void compute(const MatrixType& matrix);
 
   protected:
-    MatrixType m_eivec;
+    EigenvectorType m_eivec;
     EigenvalueType m_eivalues;
     bool m_isInitialized;
 };
@@ -97,56 +188,56 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
 {
   // this code is inspired from Jampack
   assert(matrix.cols() == matrix.rows());
-  int n = matrix.cols();
-  m_eivalues.resize(n,1);
-  m_eivec.resize(n,n);
-
-  RealScalar eps = NumTraits<RealScalar>::epsilon();
+  const int n = matrix.cols();
+  const RealScalar matrixnorm = matrix.norm();
 
-  // Reduce to complex Schur form
+  // Step 1: Do a complex Schur decomposition, A = U T U^*
+  // The eigenvalues are on the diagonal of T.
   ComplexSchur<MatrixType> schur(matrix);
-
   m_eivalues = schur.matrixT().diagonal();
 
-  m_eivec.setZero();
-
-  Scalar d2, z;
-  RealScalar norm = matrix.norm();
-
-  // compute the (normalized) eigenvectors
+  // Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T.
+  // The matrix X is unit triangular.
+  EigenvectorType X = EigenvectorType::Zero(n, n);
   for(int k=n-1 ; k>=0 ; k--)
   {
-    d2 = schur.matrixT().coeff(k,k);
-    m_eivec.coeffRef(k,k) = Scalar(1.0,0.0);
+    X.coeffRef(k,k) = ComplexScalar(1.0,0.0);
+    // Compute X(i,k) using the (i,k) entry of the equation X T = D X
     for(int i=k-1 ; i>=0 ; i--)
     {
-      m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k);
+      X.coeffRef(i,k) = -schur.matrixT().coeff(i,k);
       if(k-i-1>0)
-        m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value();
-      z = schur.matrixT().coeff(i,i) - d2;
-      if(z==Scalar(0))
-        ei_real_ref(z) = eps * norm;
-      m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z;
-
+        X.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * X.col(k).segment(i+1,k-i-1)).value();
+      ComplexScalar z = schur.matrixT().coeff(i,i) - schur.matrixT().coeff(k,k);
+      if(z==ComplexScalar(0))
+      {
+	// If the i-th and k-th eigenvalue are equal, then z equals 0. 
+	// Use a small value instead, to prevent division by zero.
+        ei_real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
+      }
+      X.coeffRef(i,k) = X.coeff(i,k) / z;
     }
-    m_eivec.col(k).normalize();
   }
 
-  m_eivec = schur.matrixU() * m_eivec;
+  // Step 3: Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
+  m_eivec = schur.matrixU() * X;
+  // .. and normalize the eigenvectors
+  for(int k=0 ; k<n ; k++)
+  {
+    m_eivec.col(k).normalize();
+  }
   m_isInitialized = true;
 
-  // sort the eigenvalues
+  // Step 4: Sort the eigenvalues
+  for (int i=0; i<n; i++)
   {
-    for (int i=0; i<n; i++)
+    int k;
+    m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
+    if (k != 0)
     {
-      int k;
-      m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
-      if (k != 0)
-      {
-        k += i;
-        std::swap(m_eivalues[k],m_eivalues[i]);
-        m_eivec.col(i).swap(m_eivec.col(k));
-      }
+      k += i;
+      std::swap(m_eivalues[k],m_eivalues[i]);
+      m_eivec.col(i).swap(m_eivec.col(k));
     }
   }
 }
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h
index a76b2a9df..7c9ab03c8 100644
--- a/Eigen/src/Eigenvalues/ComplexSchur.h
+++ b/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -29,17 +29,28 @@
 /** \eigenvalues_module \ingroup Eigenvalues_Module
   * \nonstableyet
   *
-  * \class ComplexShur
+  * \class ComplexSchur
   *
   * \brief Performs a complex Schur decomposition of a real or complex square matrix
   *
-  * Given a real or complex square matrix A, this class computes the Schur decomposition:
-  * \f$ A = U T U^*\f$ where U is a unitary complex matrix, and T is a complex upper
-  * triangular matrix.
+  * \tparam _MatrixType the type of the matrix of which we are
+  * computing the Schur decomposition; this is expected to be an
+  * instantiation of the Matrix class template.
   *
-  * The diagonal of the matrix T corresponds to the eigenvalues of the matrix A.
+  * Given a real or complex square matrix A, this class computes the
+  * Schur decomposition: \f$ A = U T U^*\f$ where U is a unitary
+  * complex matrix, and T is a complex upper triangular matrix.  The
+  * diagonal of the matrix T corresponds to the eigenvalues of the
+  * matrix A.
   *
-  * \sa class RealSchur, class EigenSolver
+  * Call the function compute() to compute the Schur decomposition of
+  * a given matrix. Alternatively, you can use the 
+  * ComplexSchur(const MatrixType&, bool) constructor which computes
+  * the Schur decomposition at construction time. Once the
+  * decomposition is computed, you can use the matrixU() and matrixT()
+  * functions to retrieve the matrices U and V in the decomposition.
+  *
+  * \sa class RealSchur, class EigenSolver, class ComplexEigenSolver
   */
 template<typename _MatrixType> class ComplexSchur
 {
@@ -52,25 +63,51 @@ template<typename _MatrixType> class ComplexSchur
       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };
+
+    /** \brief Scalar type for matrices of type \p _MatrixType. */
     typedef typename MatrixType::Scalar Scalar;
+
     typedef typename NumTraits<Scalar>::Real RealScalar;
-    typedef std::complex<RealScalar> Complex;
-    typedef Matrix<Complex, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType;
-    enum {
-      Size = MatrixType::RowsAtCompileTime
-    };
 
-    /** \brief Default Constructor.
+    /** \brief Complex scalar type for \p _MatrixType. 
       *
-      * The default constructor is useful in cases in which the user intends to
-      * perform decompositions via ComplexSchur::compute().
+      * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
+      * \c float or \c double) and just \c Scalar if #Scalar is
+      * complex.
       */
-    ComplexSchur(int size = Size==Dynamic ? 0 : Size)
+    typedef std::complex<RealScalar> ComplexScalar;
+
+    /** \brief Type for the matrices in the Schur decomposition.
+      *
+      * This is a square matrix with entries of type #ComplexScalar. 
+      * The size is the same as the size of \p _MatrixType.
+      */
+    typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType;
+
+    /** \brief Default constructor.
+      *
+      * \param [in] size  Positive integer, size of the matrix whose Schur decomposition will be computed.
+      *
+      * The default constructor is useful in cases in which the user
+      * intends to perform decompositions via compute().  The \p size
+      * parameter is only used as a hint. It is not an error to give a
+      * wrong \p size, but it may impair performance.
+      *
+      * \sa compute() for an example.
+      */
+    ComplexSchur(int size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
       : m_matT(size,size), m_matU(size,size), m_isInitialized(false), m_matUisUptodate(false)
     {}
 
-    /** Constructor computing the Schur decomposition of the matrix \a matrix.
-      * If \a skipU is true, then the matrix U is not computed. */
+    /** \brief Constructor; computes Schur decomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
+      * \param[in]  skipU   If true, then the unitary matrix U in the decomposition is not computed.
+      *
+      * This constructor calls compute() to compute the Schur decomposition.
+      *
+      * \sa matrixT() and matrixU() for examples.
+      */
     ComplexSchur(const MatrixType& matrix, bool skipU = false)
             : m_matT(matrix.rows(),matrix.cols()),
               m_matU(matrix.rows(),matrix.cols()),
@@ -80,7 +117,20 @@ template<typename _MatrixType> class ComplexSchur
       compute(matrix, skipU);
     }
 
-    /** \returns a const reference to the matrix U of the respective Schur decomposition. */
+    /** \brief Returns the unitary matrix in the Schur decomposition. 
+      *
+      * \returns A const reference to the matrix U.
+      *
+      * It is assumed that either the constructor
+      * ComplexSchur(const MatrixType& matrix, bool skipU) or the
+      * member function compute(const MatrixType& matrix, bool skipU)
+      * has been called before to compute the Schur decomposition of a
+      * matrix, and that \p skipU was set to false (the default
+      * value).
+      *
+      * Example: \include ComplexSchur_matrixU.cpp
+      * Output: \verbinclude ComplexSchur_matrixU.out
+      */
     const ComplexMatrixType& matrixU() const
     {
       ei_assert(m_isInitialized && "ComplexSchur is not initialized.");
@@ -88,24 +138,56 @@ template<typename _MatrixType> class ComplexSchur
       return m_matU;
     }
 
-    /** \returns a const reference to the matrix T of the respective Schur decomposition.
+    /** \brief Returns the triangular matrix in the Schur decomposition. 
+      *
+      * \returns A const reference to the matrix T.
+      *
+      * It is assumed that either the constructor
+      * ComplexSchur(const MatrixType& matrix, bool skipU) or the
+      * member function compute(const MatrixType& matrix, bool skipU)
+      * has been called before to compute the Schur decomposition of a
+      * matrix.
+      *
       * Note that this function returns a plain square matrix. If you want to reference
       * only the upper triangular part, use:
-      * \code schur.matrixT().triangularView<Upper>() \endcode. */
+      * \code schur.matrixT().triangularView<Upper>() \endcode 
+      *
+      * Example: \include ComplexSchur_matrixT.cpp
+      * Output: \verbinclude ComplexSchur_matrixT.out
+      */
     const ComplexMatrixType& matrixT() const
     {
-      ei_assert(m_isInitialized && "ComplexShur is not initialized.");
+      ei_assert(m_isInitialized && "ComplexSchur is not initialized.");
       return m_matT;
     }
 
-    /** Computes the Schur decomposition of the matrix \a matrix.
-      * If \a skipU is true, then the matrix U is not computed. */
+    /** \brief Computes Schur decomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
+      * \param[in]  skipU   If true, then the unitary matrix U in the decomposition is not computed.
+      *
+      * The Schur decomposition is computed by first reducing the
+      * matrix to Hessenberg form using the class
+      * HessenbergDecomposition. The Hessenberg matrix is then reduced
+      * to triangular form by performing QR iterations with a single
+      * shift. The cost of computing the Schur decomposition depends
+      * on the number of iterations; as a rough guide, it may be taken
+      * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops
+      * if \a skipU is true.
+      *
+      * Example: \include ComplexSchur_compute.cpp
+      * Output: \verbinclude ComplexSchur_compute.out
+      */
     void compute(const MatrixType& matrix, bool skipU = false);
 
   protected:
     ComplexMatrixType m_matT, m_matU;
     bool m_isInitialized;
     bool m_matUisUptodate;
+
+  private:  
+    bool subdiagonalEntryIsNeglegible(int i);
+    ComplexScalar computeShift(int iu, int iter);
 };
 
 /** Computes the principal value of the square root of the complex \a z. */
@@ -142,9 +224,63 @@ std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z)
     tim = -tim;
 
   return (std::complex<RealScalar>(tre,tim));
+}
+
+
+/** If m_matT(i+1,i) is neglegible in floating point arithmetic
+  * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and
+  * return true, else return false. */
+template<typename MatrixType>
+inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(int i)
+{
+  RealScalar d = ei_norm1(m_matT.coeff(i,i)) + ei_norm1(m_matT.coeff(i+1,i+1));
+  RealScalar sd = ei_norm1(m_matT.coeff(i+1,i));
+  if (ei_isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
+  {
+    m_matT.coeffRef(i+1,i) = ComplexScalar(0);
+    return true;
+  }
+  return false;
+}
 
+
+/** Compute the shift in the current QR iteration. */
+template<typename MatrixType>
+typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(int iu, int iter)
+{
+  if (iter == 10 || iter == 20) 
+  {
+    // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
+    return ei_abs(ei_real(m_matT.coeff(iu,iu-1))) + ei_abs(ei_real(m_matT.coeff(iu-1,iu-2)));
+  }
+
+  // compute the shift as one of the eigenvalues of t, the 2x2
+  // diagonal block on the bottom of the active submatrix
+  Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
+  RealScalar normt = t.cwiseAbs().sum();
+  t /= normt;     // the normalization by sf is to avoid under/overflow
+
+  ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
+  ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
+  ComplexScalar disc = ei_sqrt(c*c + RealScalar(4)*b);
+  ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
+  ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
+  ComplexScalar eival1 = (trace + disc) / RealScalar(2);
+  ComplexScalar eival2 = (trace - disc) / RealScalar(2);
+
+  if(ei_norm1(eival1) > ei_norm1(eival2))
+    eival2 = det / eival1;
+  else
+    eival1 = det / eival2;
+
+  // choose the eigenvalue closest to the bottom entry of the diagonal
+  if(ei_norm1(eival1-t.coeff(1,1)) < ei_norm1(eival2-t.coeff(1,1)))
+    return normt * eival1;
+  else
+    return normt * eival2;
 }
 
+
 template<typename MatrixType>
 void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
 {
@@ -156,7 +292,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
   if(n==1)
   {
     m_matU = ComplexMatrixType::Identity(1,1);
-    if(!skipU) m_matT = matrix;
+    if(!skipU) m_matT = matrix.template cast<ComplexScalar>();
     m_isInitialized = true;
     m_matUisUptodate = !skipU;
     return;
@@ -166,9 +302,8 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
   // TODO skip Q if skipU = true
   HessenbergDecomposition<MatrixType> hess(matrix);
 
-  m_matT = hess.matrixH();
-  if(!skipU) m_matU = hess.matrixQ();
-
+  m_matT = hess.matrixH().template cast<ComplexScalar>();
+  if(!skipU)  m_matU = hess.matrixQ().template cast<ComplexScalar>();
 
   // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
 
@@ -176,109 +311,62 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
   // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 
   // Rows il,...,iu is the part we are working on (the active submatrix).
   // Rows iu+1,...,end are already brought in triangular form.
-
   int iu = m_matT.cols() - 1;
   int il;
-  RealScalar d,sd,sf;
-  Complex c,b,disc,r1,r2,kappa;
-
-  RealScalar eps = NumTraits<RealScalar>::epsilon();
+  int iter = 0; // number of iterations we are working on the (iu,iu) element
 
-  int iter = 0;
   while(true)
   {
     // find iu, the bottom row of the active submatrix
     while(iu > 0)
     {
-      d = ei_norm1(m_matT.coeff(iu,iu)) + ei_norm1(m_matT.coeff(iu-1,iu-1));
-      sd = ei_norm1(m_matT.coeff(iu,iu-1));
-
-      if(!ei_isMuchSmallerThan(sd,d,eps))
-        break;
-
-      m_matT.coeffRef(iu,iu-1) = Complex(0);
+      if(!subdiagonalEntryIsNeglegible(iu-1)) break;
       iter = 0;
       --iu;
     }
+
+    // if iu is zero then we are done; the whole matrix is triangularized
     if(iu==0) break;
-    iter++;
 
-    if(iter >= 30)
-    {
-      // FIXME : what to do when iter==MAXITER ??
-      //std::cerr << "MAXITER" << std::endl;
-      return;
-    }
+    // if we spent 30 iterations on the current element, we give up
+    iter++;
+    if(iter >= 30) break;
 
     // find il, the top row of the active submatrix
     il = iu-1;
-    while(il > 0)
+    while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
     {
-      // check if the current 2x2 block on the diagonal is upper triangular
-      d = ei_norm1(m_matT.coeff(il,il)) + ei_norm1(m_matT.coeff(il-1,il-1));
-      sd = ei_norm1(m_matT.coeff(il,il-1));
-
-      if(ei_isMuchSmallerThan(sd,d,eps))
-        break;
-
       --il;
     }
 
-    if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0);
-
-    // compute the shift kappa as one of the eigenvalues of the 2x2
-    // diagonal block on the bottom of the active submatrix
-
-    Matrix<Scalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
-    sf = t.cwiseAbs().sum();
-    t /= sf;     // the normalization by sf is to avoid under/overflow
+    /* perform the QR step using Givens rotations. The first rotation
+       creates a bulge; the (il+2,il) element becomes nonzero. This
+       bulge is chased down to the bottom of the active submatrix. */
 
-    b = t.coeff(0,1) * t.coeff(1,0);
-    c = t.coeff(0,0) - t.coeff(1,1);
-    disc = ei_sqrt(c*c + RealScalar(4)*b);
+    ComplexScalar shift = computeShift(iu, iter);
+    PlanarRotation<ComplexScalar> rot;
+    rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
+    m_matT.block(0,il,n,n-il).applyOnTheLeft(il, il+1, rot.adjoint());
+    m_matT.block(0,0,std::min(il+2,iu)+1,n).applyOnTheRight(il, il+1, rot);
+    if(!skipU) m_matU.applyOnTheRight(il, il+1, rot);
 
-    c = t.coeff(0,0) * t.coeff(1,1) - b;
-    b = t.coeff(0,0) + t.coeff(1,1);
-    r1 = (b+disc)/RealScalar(2);
-    r2 = (b-disc)/RealScalar(2);
-
-    if(ei_norm1(r1) > ei_norm1(r2))
-      r2 = c/r1;
-    else
-      r1 = c/r2;
-
-    if(ei_norm1(r1-t.coeff(1,1)) < ei_norm1(r2-t.coeff(1,1)))
-      kappa = sf * r1;
-    else
-      kappa = sf * r2;
-    
-    if (iter == 10 || iter == 20) 
-    {
-      // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
-      kappa = ei_abs(ei_real(m_matT.coeff(iu,iu-1))) + ei_abs(ei_real(m_matT.coeff(iu-1,iu-2)));
-    }
-
-    // perform the QR step using Givens rotations
-    PlanarRotation<Complex> rot;
-    rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il));
-
-    for(int i=il ; i<iu ; i++)
+    for(int i=il+1 ; i<iu ; i++)
     {
+      rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
+      m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
       m_matT.block(0,i,n,n-i).applyOnTheLeft(i, i+1, rot.adjoint());
       m_matT.block(0,0,std::min(i+2,iu)+1,n).applyOnTheRight(i, i+1, rot);
       if(!skipU) m_matU.applyOnTheRight(i, i+1, rot);
-
-      if(i != iu-1)
-      {
-        int i1 = i+1;
-        int i2 = i+2;
-
-        rot.makeGivens(m_matT.coeffRef(i1,i), m_matT.coeffRef(i2,i), &m_matT.coeffRef(i1,i));
-        m_matT.coeffRef(i2,i) = Complex(0);
-      }
     }
   }
 
+  if(iter >= 30) 
+  {
+    // FIXME : what to do when iter==MAXITER ??
+    // std::cerr << "MAXITER" << std::endl;
+    return;
+  }
+
   m_isInitialized = true;
   m_matUisUptodate = !skipU;
 }
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index 3fcd8295c..32dee92e9 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -2,6 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -25,20 +26,53 @@
 #ifndef EIGEN_EIGENSOLVER_H
 #define EIGEN_EIGENSOLVER_H
 
+#include "./RealSchur.h"
+
 /** \eigenvalues_module \ingroup Eigenvalues_Module
   * \nonstableyet
   *
   * \class EigenSolver
   *
-  * \brief Eigen values/vectors solver for non selfadjoint matrices
+  * \brief Computes eigenvalues and eigenvectors of general matrices
+  *
+  * \tparam _MatrixType the type of the matrix of which we are computing the
+  * eigendecomposition; this is expected to be an instantiation of the Matrix
+  * class template. Currently, only real matrices are supported.
+  *
+  * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
+  * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v \f$.  If
+  * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and
+  * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V =
+  * V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we
+  * have \f$ A = V D V^{-1} \f$. This is called the eigendecomposition.
   *
-  * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
+  * The eigenvalues and eigenvectors of a matrix may be complex, even when the
+  * matrix is real. However, we can choose real matrices \f$ V \f$ and \f$ D
+  * \f$ satisfying \f$ A V = V D \f$, just like the eigendecomposition, if the
+  * matrix \f$ D \f$ is not required to be diagonal, but if it is allowed to
+  * have blocks of the form
+  * \f[ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f]
+  * (where \f$ u \f$ and \f$ v \f$ are real numbers) on the diagonal.  These
+  * blocks correspond to complex eigenvalue pairs \f$ u \pm iv \f$. We call
+  * this variant of the eigendecomposition the pseudo-eigendecomposition.
   *
-  * Currently it only support real matrices.
+  * Call the function compute() to compute the eigenvalues and eigenvectors of
+  * a given matrix. Alternatively, you can use the
+  * EigenSolver(const MatrixType&) constructor which computes the eigenvalues
+  * and eigenvectors at construction time. Once the eigenvalue and eigenvectors
+  * are computed, they can be retrieved with the eigenvalues() and
+  * eigenvectors() functions. The pseudoEigenvalueMatrix() and
+  * pseudoEigenvectors() methods allow the construction of the
+  * pseudo-eigendecomposition.
   *
-  * \note this code was adapted from JAMA (public domain)
+  * The documentation for EigenSolver(const MatrixType&) contains an example of
+  * the typical use of this class.
   *
-  * \sa MatrixBase::eigenvalues(), SelfAdjointEigenSolver
+  * \note The implementation is adapted from
+  * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
+  * Their code is based on EISPACK.
+  *
+  * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver
   */
 template<typename _MatrixType> class EigenSolver
 {
@@ -52,21 +86,54 @@ template<typename _MatrixType> class EigenSolver
       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };
+
+    /** \brief Scalar type for matrices of type \p _MatrixType. */
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<Scalar>::Real RealScalar;
-    typedef std::complex<RealScalar> Complex;
-    typedef typename ei_plain_col_type<MatrixType, Complex>::type EigenvalueType;
-    typedef Matrix<Complex, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType;
-    typedef typename ei_plain_col_type<MatrixType, RealScalar>::type RealVectorType;
-
-    /**
-    * \brief Default Constructor.
-    *
-    * The default constructor is useful in cases in which the user intends to
-    * perform decompositions via EigenSolver::compute(const MatrixType&).
-    */
+
+    /** \brief Complex scalar type for \p _MatrixType. 
+      *
+      * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
+      * \c float or \c double) and just \c Scalar if #Scalar is
+      * complex.
+      */
+    typedef std::complex<RealScalar> ComplexScalar;
+
+    /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 
+      *
+      * This is a column vector with entries of type #ComplexScalar.
+      * The length of the vector is the size of \p _MatrixType.
+      */
+    typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
+
+    /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 
+      *
+      * This is a square matrix with entries of type #ComplexScalar. 
+      * The size is the same as the size of \p _MatrixType.
+      */
+    typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
+
+    /** \brief Default constructor.
+      *
+      * The default constructor is useful in cases in which the user intends to
+      * perform decompositions via EigenSolver::compute(const MatrixType&).
+      *
+      * \sa compute() for an example.
+      */
     EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {}
 
+    /** \brief Constructor; computes eigendecomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
+      *
+      * This constructor calls compute() to compute the eigenvalues
+      * and eigenvectors.
+      *
+      * Example: \include EigenSolver_EigenSolver_MatrixType.cpp
+      * Output: \verbinclude EigenSolver_EigenSolver_MatrixType.out
+      *
+      * \sa compute()
+      */
     EigenSolver(const MatrixType& matrix)
       : m_eivec(matrix.rows(), matrix.cols()),
         m_eivalues(matrix.cols()),
@@ -75,39 +142,42 @@ template<typename _MatrixType> class EigenSolver
       compute(matrix);
     }
 
+    /** \brief Returns the eigenvectors of given matrix. 
+      *
+      * \returns  %Matrix whose columns are the (possibly complex) eigenvectors.
+      *
+      * \pre Either the constructor EigenSolver(const MatrixType&) or the
+      * member function compute(const MatrixType&) has been called before.
+      *
+      * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
+      * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
+      * eigenvectors are normalized to have (Euclidean) norm equal to one. The
+      * matrix returned by this function is the matrix \f$ V \f$ in the
+      * eigendecomposition \f$ A = V D V^{-1} \f$, if it exists.
+      *
+      * Example: \include EigenSolver_eigenvectors.cpp
+      * Output: \verbinclude EigenSolver_eigenvectors.out
+      *
+      * \sa eigenvalues(), pseudoEigenvectors()
+      */
+    EigenvectorsType eigenvectors() const;
 
-    EigenvectorType eigenvectors(void) const;
-
-    /** \returns a real matrix V of pseudo eigenvectors.
-      *
-      * Let D be the block diagonal matrix with the real eigenvalues in 1x1 blocks,
-      * and any complex values u+iv in 2x2 blocks [u v ; -v u]. Then, the matrices D
-      * and V satisfy A*V = V*D.
-      *
-      * More precisely, if the diagonal matrix of the eigen values is:\n
-      * \f$
-      * \left[ \begin{array}{cccccc}
-      * u+iv &      &      &      &   &   \\
-      *      & u-iv &      &      &   &   \\
-      *      &      & a+ib &      &   &   \\
-      *      &      &      & a-ib &   &   \\
-      *      &      &      &      & x &   \\
-      *      &      &      &      &   & y \\
-      * \end{array} \right]
-      * \f$ \n
-      * then, we have:\n
-      * \f$
-      * D =\left[ \begin{array}{cccccc}
-      *  u & v &    &   &   &   \\
-      * -v & u &    &   &   &   \\
-      *    &   &  a & b &   &   \\
-      *    &   & -b & a &   &   \\
-      *    &   &    &   & x &   \\
-      *    &   &    &   &   & y \\
-      * \end{array} \right]
-      * \f$
-      *
-      * \sa pseudoEigenvalueMatrix()
+    /** \brief Returns the pseudo-eigenvectors of given matrix. 
+      *
+      * \returns  Const reference to matrix whose columns are the pseudo-eigenvectors.
+      *
+      * \pre Either the constructor EigenSolver(const MatrixType&) or
+      * the member function compute(const MatrixType&) has been called
+      * before.
+      *
+      * The real matrix \f$ V \f$ returned by this function and the
+      * block-diagonal matrix \f$ D \f$ returned by pseudoEigenvalueMatrix()
+      * satisfy \f$ AV = VD \f$.
+      *
+      * Example: \include EigenSolver_pseudoEigenvectors.cpp
+      * Output: \verbinclude EigenSolver_pseudoEigenvectors.out
+      *
+      * \sa pseudoEigenvalueMatrix(), eigenvectors()
       */
     const MatrixType& pseudoEigenvectors() const
     {
@@ -115,21 +185,72 @@ template<typename _MatrixType> class EigenSolver
       return m_eivec;
     }
 
+    /** \brief Returns the block-diagonal matrix in the pseudo-eigendecomposition.
+      *
+      * \returns  A block-diagonal matrix.
+      *
+      * \pre Either the constructor EigenSolver(const MatrixType&) or the
+      * member function compute(const MatrixType&) has been called before.
+      *
+      * The matrix \f$ D \f$ returned by this function is real and
+      * block-diagonal. The blocks on the diagonal are either 1-by-1 or 2-by-2
+      * blocks of the form
+      * \f$ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f$.
+      * The matrix \f$ D \f$ and the matrix \f$ V \f$ returned by
+      * pseudoEigenvectors() satisfy \f$ AV = VD \f$.
+      *
+      * \sa pseudoEigenvectors() for an example, eigenvalues()
+      */
     MatrixType pseudoEigenvalueMatrix() const;
 
-    /** \returns the eigenvalues as a column vector */
+    /** \brief Returns the eigenvalues of given matrix. 
+      *
+      * \returns Column vector containing the eigenvalues.
+      *
+      * \pre Either the constructor EigenSolver(const MatrixType&) or the
+      * member function compute(const MatrixType&) has been called before.
+      *
+      * The eigenvalues are repeated according to their algebraic multiplicity,
+      * so there are as many eigenvalues as rows in the matrix.
+      *
+      * Example: \include EigenSolver_eigenvalues.cpp
+      * Output: \verbinclude EigenSolver_eigenvalues.out
+      *
+      * \sa eigenvectors(), pseudoEigenvalueMatrix(),
+      *     MatrixBase::eigenvalues()
+      */
     EigenvalueType eigenvalues() const
     {
       ei_assert(m_isInitialized && "EigenSolver is not initialized.");
       return m_eivalues;
     }
 
+    /** \brief Computes eigendecomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
+      * \returns    Reference to \c *this
+      *
+      * This function computes the eigenvalues and eigenvectors of \p matrix.
+      * The eigenvalues() and eigenvectors() functions can be used to retrieve
+      * the computed eigendecomposition.
+      *
+      * The matrix is first reduced to real Schur form using the RealSchur
+      * class. The Schur decomposition is then used to compute the eigenvalues
+      * and eigenvectors.
+      *
+      * The cost of the computation is dominated by the cost of the Schur
+      * decomposition, which is very approximately \f$ 25n^3 \f$ where 
+      * \f$ n \f$ is the size of the matrix.
+      *
+      * This method reuses of the allocated data in the EigenSolver object.
+      *
+      * Example: \include EigenSolver_compute.cpp
+      * Output: \verbinclude EigenSolver_compute.out
+      */
     EigenSolver& compute(const MatrixType& matrix);
 
   private:
-
-    void orthes(MatrixType& matH, RealVectorType& ort);
-    void hqr2(MatrixType& matH);
+    void hqr2_step2(MatrixType& matH);
 
   protected:
     MatrixType m_eivec;
@@ -137,10 +258,6 @@ template<typename _MatrixType> class EigenSolver
     bool m_isInitialized;
 };
 
-/** \returns the real block diagonal matrix D of the eigenvalues.
-  *
-  * See pseudoEigenvectors() for the details.
-  */
 template<typename MatrixType>
 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
 {
@@ -161,30 +278,26 @@ MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
   return matD;
 }
 
-/** \returns the normalized complex eigenvectors as a matrix of column vectors.
-  *
-  * \sa eigenvalues(), pseudoEigenvectors()
-  */
 template<typename MatrixType>
-typename EigenSolver<MatrixType>::EigenvectorType EigenSolver<MatrixType>::eigenvectors(void) const
+typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
 {
   ei_assert(m_isInitialized && "EigenSolver is not initialized.");
   int n = m_eivec.cols();
-  EigenvectorType matV(n,n);
+  EigenvectorsType matV(n,n);
   for (int j=0; j<n; ++j)
   {
     if (ei_isMuchSmallerThan(ei_abs(ei_imag(m_eivalues.coeff(j))), ei_abs(ei_real(m_eivalues.coeff(j)))))
     {
       // we have a real eigen value
-      matV.col(j) = m_eivec.col(j).template cast<Complex>();
+      matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
     }
     else
     {
       // we have a pair of complex eigen values
       for (int i=0; i<n; ++i)
       {
-        matV.coeffRef(i,j)   = Complex(m_eivec.coeff(i,j),  m_eivec.coeff(i,j+1));
-        matV.coeffRef(i,j+1) = Complex(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
+        matV.coeffRef(i,j)   = ComplexScalar(m_eivec.coeff(i,j),  m_eivec.coeff(i,j+1));
+        matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
       }
       matV.col(j).normalize();
       matV.col(j+1).normalize();
@@ -198,84 +311,37 @@ template<typename MatrixType>
 EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix)
 {
   assert(matrix.cols() == matrix.rows());
-  int n = matrix.cols();
-  m_eivalues.resize(n,1);
-  m_eivec.resize(n,n);
-
-  MatrixType matH = matrix;
-  RealVectorType ort(n);
-
-  // Reduce to Hessenberg form.
-  orthes(matH, ort);
 
-  // Reduce Hessenberg to real Schur form.
-  hqr2(matH);
+  // Reduce to real Schur form.
+  RealSchur<MatrixType> rs(matrix);
+  MatrixType matT = rs.matrixT();
+  m_eivec = rs.matrixU();
 
-  m_isInitialized = true;
-  return *this;
-}
-
-// Nonsymmetric reduction to Hessenberg form.
-template<typename MatrixType>
-void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
-{
-  //  This is derived from the Algol procedures orthes and ortran,
-  //  by Martin and Wilkinson, Handbook for Auto. Comp.,
-  //  Vol.ii-Linear Algebra, and the corresponding
-  //  Fortran subroutines in EISPACK.
-
-  int n = m_eivec.cols();
-  int low = 0;
-  int high = n-1;
-
-  for (int m = low+1; m <= high-1; ++m)
+  // Compute eigenvalues from matT
+  m_eivalues.resize(matrix.cols());
+  int i = 0;
+  while (i < matrix.cols()) 
   {
-    // Scale column.
-    RealScalar scale = matH.block(m, m-1, high-m+1, 1).cwiseAbs().sum();
-    if (scale != 0.0)
+    if (i == matrix.cols() - 1 || matT.coeff(i+1, i) == Scalar(0)) 
     {
-      // Compute Householder transformation.
-      RealScalar h = 0.0;
-      // FIXME could be rewritten, but this one looks better wrt cache
-      for (int i = high; i >= m; i--)
-      {
-        ort.coeffRef(i) = matH.coeff(i,m-1)/scale;
-        h += ort.coeff(i) * ort.coeff(i);
-      }
-      RealScalar g = ei_sqrt(h);
-      if (ort.coeff(m) > 0)
-        g = -g;
-      h = h - ort.coeff(m) * g;
-      ort.coeffRef(m) = ort.coeff(m) - g;
-
-      // Apply Householder similarity transformation
-      // H = (I-u*u'/h)*H*(I-u*u')/h)
-      int bSize = high-m+1;
-      matH.block(m, m, bSize, n-m).noalias() -= ((ort.segment(m, bSize)/h)
-        * (ort.segment(m, bSize).transpose() *  matH.block(m, m, bSize, n-m)));
-
-      matH.block(0, m, high+1, bSize).noalias() -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize))
-        * (ort.segment(m, bSize)/h).transpose());
-
-      ort.coeffRef(m) = scale*ort.coeff(m);
-      matH.coeffRef(m,m-1) = scale*g;
+      m_eivalues.coeffRef(i) = matT.coeff(i, i);
+      ++i;
     }
-  }
-
-  // Accumulate transformations (Algol's ortran).
-  m_eivec.setIdentity();
-
-  for (int m = high-1; m >= low+1; m--)
-  {
-    if (matH.coeff(m,m-1) != 0.0)
+    else
     {
-      ort.segment(m+1, high-m) = matH.col(m-1).segment(m+1, high-m);
-
-      int bSize = high-m+1;
-      m_eivec.block(m, m, bSize, bSize).noalias() += ( (ort.segment(m, bSize) /  (matH.coeff(m,m-1) * ort.coeff(m)))
-        * (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) );
+      Scalar p = Scalar(0.5) * (matT.coeff(i, i) - matT.coeff(i+1, i+1));
+      Scalar z = ei_sqrt(ei_abs(p * p + matT.coeff(i+1, i) * matT.coeff(i, i+1)));
+      m_eivalues.coeffRef(i)   = ComplexScalar(matT.coeff(i+1, i+1) + p, z);
+      m_eivalues.coeffRef(i+1) = ComplexScalar(matT.coeff(i+1, i+1) + p, -z);
+      i += 2;
     }
   }
+  
+  // Compute eigenvectors.
+  hqr2_step2(matT);
+
+  m_isInitialized = true;
+  return *this;
 }
 
 // Complex scalar division.
@@ -298,289 +364,29 @@ std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
 }
 
 
-// Nonsymmetric reduction from Hessenberg to real Schur form.
 template<typename MatrixType>
-void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
+void EigenSolver<MatrixType>::hqr2_step2(MatrixType& matH)
 {
-  //  This is derived from the Algol procedure hqr2,
-  //  by Martin and Wilkinson, Handbook for Auto. Comp.,
-  //  Vol.ii-Linear Algebra, and the corresponding
-  //  Fortran subroutine in EISPACK.
-
-  // Initialize
-  int nn = m_eivec.cols();
-  int n = nn-1;
-  int low = 0;
-  int high = nn-1;
-  Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52));
-  Scalar exshift = 0.0;
-  Scalar p=0,q=0,r=0,s=0,z=0,t,w,x,y;
-
-  // Store roots isolated by balanc and compute matrix norm
-  // FIXME to be efficient the following would requires a triangular reduxion code
-  // Scalar norm = matH.upper().cwiseAbs().sum() + matH.corner(BottomLeft,n,n).diagonal().cwiseAbs().sum();
+  const int nn = m_eivec.cols();
+  const int low = 0;
+  const int high = nn-1;
+  const Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52));
+  Scalar p, q, r=0, s=0, t, w, x, y, z=0;
+
+  // inefficient! this is already computed in RealSchur
   Scalar norm = 0.0;
   for (int j = 0; j < nn; ++j)
   {
-    // FIXME what's the purpose of the following since the condition is always false
-    if ((j < low) || (j > high))
-    {
-      m_eivalues.coeffRef(j) = Complex(matH.coeff(j,j), 0.0);
-    }
     norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwiseAbs().sum();
   }
-
-  // Outer loop over eigenvalue index
-  int iter = 0;
-  while (n >= low)
-  {
-    // Look for single small sub-diagonal element
-    int l = n;
-    while (l > low)
-    {
-      s = ei_abs(matH.coeff(l-1,l-1)) + ei_abs(matH.coeff(l,l));
-      if (s == 0.0)
-          s = norm;
-      if (ei_abs(matH.coeff(l,l-1)) < eps * s)
-        break;
-      l--;
-    }
-
-    // Check for convergence
-    // One root found
-    if (l == n)
-    {
-      matH.coeffRef(n,n) = matH.coeff(n,n) + exshift;
-      m_eivalues.coeffRef(n) = Complex(matH.coeff(n,n), 0.0);
-      n--;
-      iter = 0;
-    }
-    else if (l == n-1) // Two roots found
-    {
-      w = matH.coeff(n,n-1) * matH.coeff(n-1,n);
-      p = (matH.coeff(n-1,n-1) - matH.coeff(n,n)) * Scalar(0.5);
-      q = p * p + w;
-      z = ei_sqrt(ei_abs(q));
-      matH.coeffRef(n,n) = matH.coeff(n,n) + exshift;
-      matH.coeffRef(n-1,n-1) = matH.coeff(n-1,n-1) + exshift;
-      x = matH.coeff(n,n);
-
-      // Scalar pair
-      if (q >= 0)
-      {
-        if (p >= 0)
-          z = p + z;
-        else
-          z = p - z;
-
-        m_eivalues.coeffRef(n-1) = Complex(x + z, 0.0);
-        m_eivalues.coeffRef(n) = Complex(z!=0.0 ? x - w / z : m_eivalues.coeff(n-1).real(), 0.0);
-
-        x = matH.coeff(n,n-1);
-        s = ei_abs(x) + ei_abs(z);
-        p = x / s;
-        q = z / s;
-        r = ei_sqrt(p * p+q * q);
-        p = p / r;
-        q = q / r;
-
-        // Row modification
-        for (int j = n-1; j < nn; ++j)
-        {
-          z = matH.coeff(n-1,j);
-          matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j);
-          matH.coeffRef(n,j) = q * matH.coeff(n,j) - p * z;
-        }
-
-        // Column modification
-        for (int i = 0; i <= n; ++i)
-        {
-          z = matH.coeff(i,n-1);
-          matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n);
-          matH.coeffRef(i,n) = q * matH.coeff(i,n) - p * z;
-        }
-
-        // Accumulate transformations
-        for (int i = low; i <= high; ++i)
-        {
-          z = m_eivec.coeff(i,n-1);
-          m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n);
-          m_eivec.coeffRef(i,n) = q * m_eivec.coeff(i,n) - p * z;
-        }
-      }
-      else // Complex pair
-      {
-        m_eivalues.coeffRef(n-1) = Complex(x + p, z);
-        m_eivalues.coeffRef(n)   = Complex(x + p, -z);
-      }
-      n = n - 2;
-      iter = 0;
-    }
-    else // No convergence yet
-    {
-      // Form shift
-      x = matH.coeff(n,n);
-      y = 0.0;
-      w = 0.0;
-      if (l < n)
-      {
-          y = matH.coeff(n-1,n-1);
-          w = matH.coeff(n,n-1) * matH.coeff(n-1,n);
-      }
-
-      // Wilkinson's original ad hoc shift
-      if (iter == 10)
-      {
-        exshift += x;
-        for (int i = low; i <= n; ++i)
-          matH.coeffRef(i,i) -= x;
-        s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2));
-        x = y = Scalar(0.75) * s;
-        w = Scalar(-0.4375) * s * s;
-      }
-
-      // MATLAB's new ad hoc shift
-      if (iter == 30)
-      {
-        s = Scalar((y - x) / 2.0);
-        s = s * s + w;
-        if (s > 0)
-        {
-          s = ei_sqrt(s);
-          if (y < x)
-            s = -s;
-          s = Scalar(x - w / ((y - x) / 2.0 + s));
-          for (int i = low; i <= n; ++i)
-            matH.coeffRef(i,i) -= s;
-          exshift += s;
-          x = y = w = Scalar(0.964);
-        }
-      }
-
-      iter = iter + 1;   // (Could check iteration count here.)
-
-      // Look for two consecutive small sub-diagonal elements
-      int m = n-2;
-      while (m >= l)
-      {
-        z = matH.coeff(m,m);
-        r = x - z;
-        s = y - z;
-        p = (r * s - w) / matH.coeff(m+1,m) + matH.coeff(m,m+1);
-        q = matH.coeff(m+1,m+1) - z - r - s;
-        r = matH.coeff(m+2,m+1);
-        s = ei_abs(p) + ei_abs(q) + ei_abs(r);
-        p = p / s;
-        q = q / s;
-        r = r / s;
-        if (m == l) {
-          break;
-        }
-        if (ei_abs(matH.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) <
-          eps * (ei_abs(p) * (ei_abs(matH.coeff(m-1,m-1)) + ei_abs(z) +
-          ei_abs(matH.coeff(m+1,m+1)))))
-        {
-          break;
-        }
-        m--;
-      }
-
-      for (int i = m+2; i <= n; ++i)
-      {
-        matH.coeffRef(i,i-2) = 0.0;
-        if (i > m+2)
-          matH.coeffRef(i,i-3) = 0.0;
-      }
-
-      // Double QR step involving rows l:n and columns m:n
-      for (int k = m; k <= n-1; ++k)
-      {
-        int notlast = (k != n-1);
-        if (k != m) {
-          p = matH.coeff(k,k-1);
-          q = matH.coeff(k+1,k-1);
-          r = notlast ? matH.coeff(k+2,k-1) : Scalar(0);
-          x = ei_abs(p) + ei_abs(q) + ei_abs(r);
-          if (x != 0.0)
-          {
-            p = p / x;
-            q = q / x;
-            r = r / x;
-          }
-        }
-
-        if (x == 0.0)
-          break;
-
-        s = ei_sqrt(p * p + q * q + r * r);
-
-        if (p < 0)
-          s = -s;
-
-        if (s != 0)
-        {
-          if (k != m)
-            matH.coeffRef(k,k-1) = -s * x;
-          else if (l != m)
-            matH.coeffRef(k,k-1) = -matH.coeff(k,k-1);
-
-          p = p + s;
-          x = p / s;
-          y = q / s;
-          z = r / s;
-          q = q / p;
-          r = r / p;
-
-          // Row modification
-          for (int j = k; j < nn; ++j)
-          {
-            p = matH.coeff(k,j) + q * matH.coeff(k+1,j);
-            if (notlast)
-            {
-              p = p + r * matH.coeff(k+2,j);
-              matH.coeffRef(k+2,j) = matH.coeff(k+2,j) - p * z;
-            }
-            matH.coeffRef(k,j) = matH.coeff(k,j) - p * x;
-            matH.coeffRef(k+1,j) = matH.coeff(k+1,j) - p * y;
-          }
-
-          // Column modification
-          for (int i = 0; i <= std::min(n,k+3); ++i)
-          {
-            p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1);
-            if (notlast)
-            {
-              p = p + z * matH.coeff(i,k+2);
-              matH.coeffRef(i,k+2) = matH.coeff(i,k+2) - p * r;
-            }
-            matH.coeffRef(i,k) = matH.coeff(i,k) - p;
-            matH.coeffRef(i,k+1) = matH.coeff(i,k+1) - p * q;
-          }
-
-          // Accumulate transformations
-          for (int i = low; i <= high; ++i)
-          {
-            p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1);
-            if (notlast)
-            {
-              p = p + z * m_eivec.coeff(i,k+2);
-              m_eivec.coeffRef(i,k+2) = m_eivec.coeff(i,k+2) - p * r;
-            }
-            m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) - p;
-            m_eivec.coeffRef(i,k+1) = m_eivec.coeff(i,k+1) - p * q;
-          }
-        }  // (s != 0)
-      }  // k loop
-    }  // check convergence
-  }  // while (n >= low)
-
+  
   // Backsubstitute to find vectors of upper triangular form
   if (norm == 0.0)
   {
       return;
   }
 
-  for (n = nn-1; n >= 0; n--)
+  for (int n = nn-1; n >= 0; n--)
   {
     p = m_eivalues.coeff(n).real();
     q = m_eivalues.coeff(n).imag();
diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
index de833528d..3e8c4c64d 100644
--- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h
+++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
@@ -30,14 +30,30 @@
   *
   * \class HessenbergDecomposition
   *
-  * \brief Reduces a squared matrix to an Hessemberg form
+  * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation
   *
-  * \param MatrixType the type of the matrix of which we are computing the Hessenberg decomposition
+  * \tparam _MatrixType the type of the matrix of which we are computing the Hessenberg decomposition
   *
-  * This class performs an Hessenberg decomposition of a matrix \f$ A \f$ such that:
-  * \f$ A = Q H Q^* \f$ where \f$ Q \f$ is unitary and \f$ H \f$ a Hessenberg matrix.
+  * This class performs an Hessenberg decomposition of a matrix \f$ A \f$. In
+  * the real case, the Hessenberg decomposition consists of an orthogonal
+  * matrix \f$ Q \f$ and a Hessenberg matrix \f$ H \f$ such that \f$ A = Q H
+  * Q^T \f$. An orthogonal matrix is a matrix whose inverse equals its
+  * transpose (\f$ Q^{-1} = Q^T \f$). A Hessenberg matrix has zeros below the
+  * subdiagonal, so it is almost upper triangular. The Hessenberg decomposition
+  * of a complex matrix is \f$ A = Q H Q^* \f$ with \f$ Q \f$ unitary (that is,
+  * \f$ Q^{-1} = Q^* \f$).
   *
-  * \sa class Tridiagonalization, class Qr
+  * Call the function compute() to compute the Hessenberg decomposition of a
+  * given matrix. Alternatively, you can use the 
+  * HessenbergDecomposition(const MatrixType&) constructor which computes the
+  * Hessenberg decomposition at construction time. Once the decomposition is
+  * computed, you can use the matrixH() and matrixQ() functions to construct
+  * the matrices H and Q in the decomposition. 
+  *
+  * The documentation for matrixH() contains an example of the typical use of
+  * this class.
+  *
+  * \sa class ComplexSchur, class Tridiagonalization, \ref QR_Module "QR Module"
   */
 template<typename _MatrixType> class HessenbergDecomposition
 {
@@ -51,13 +67,28 @@ template<typename _MatrixType> class HessenbergDecomposition
       MaxSize = MatrixType::MaxRowsAtCompileTime,
       MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
     };
+
+    /** \brief Scalar type for matrices of type \p _MatrixType. */
     typedef typename MatrixType::Scalar Scalar;
-    typedef typename NumTraits<Scalar>::Real RealScalar;
+
+    /** \brief Type for vector of Householder coefficients.
+      *
+      * This is column vector with entries of type #Scalar. The length of the
+      * vector is one less than the size of \p _MatrixType, if it is a
+      * fixed-side type.
+      */
     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
-    typedef typename ei_plain_col_type<MatrixType>::type VectorType;
 
-    /** This constructor initializes a HessenbergDecomposition object for
-      * further use with HessenbergDecomposition::compute()
+    /** \brief Default constructor; the decomposition will be computed later.
+      *
+      * \param [in] size  The size of the matrix whose Hessenberg decomposition will be computed.
+      *
+      * The default constructor is useful in cases in which the user intends to
+      * perform decompositions via compute().  The \p size parameter is only
+      * used as a hint. It is not an error to give a wrong \p size, but it may
+      * impair performance.
+      *
+      * \sa compute() for an example.
       */
     HessenbergDecomposition(int size = Size==Dynamic ? 2 : Size)
       : m_matrix(size,size)
@@ -66,6 +97,15 @@ template<typename _MatrixType> class HessenbergDecomposition
         m_hCoeffs.resize(size-1);
     }
 
+    /** \brief Constructor; computes Hessenberg decomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
+      *
+      * This constructor calls compute() to compute the Hessenberg
+      * decomposition.
+      *
+      * \sa matrixH() for an example.
+      */
     HessenbergDecomposition(const MatrixType& matrix)
       : m_matrix(matrix)
     {
@@ -75,9 +115,21 @@ template<typename _MatrixType> class HessenbergDecomposition
       _compute(m_matrix, m_hCoeffs);
     }
 
-    /** Computes or re-compute the Hessenberg decomposition for the matrix \a matrix.
+    /** \brief Computes Hessenberg decomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
+      *
+      * The Hessenberg decomposition is computed by bringing the columns of the
+      * matrix successively in the required form using Householder reflections
+      * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix
+      * Computations</i>). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$
+      * denotes the size of the given matrix.
       *
-      * This method allows to re-use the allocated data.
+      * This method reuses of the allocated data in the HessenbergDecomposition
+      * object.
+      *
+      * Example: \include HessenbergDecomposition_compute.cpp
+      * Output: \verbinclude HessenbergDecomposition_compute.out
       */
     void compute(const MatrixType& matrix)
     {
@@ -88,36 +140,95 @@ template<typename _MatrixType> class HessenbergDecomposition
       _compute(m_matrix, m_hCoeffs);
     }
 
-    /** \returns a const reference to the householder coefficients allowing to
-      * reconstruct the matrix Q from the packed data.
+    /** \brief Returns the Householder coefficients.
+      *
+      * \returns a const reference to the vector of Householder coefficients
+      *
+      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
+      * or the member function compute(const MatrixType&) has been called
+      * before to compute the Hessenberg decomposition of a matrix.
+      *
+      * The Householder coefficients allow the reconstruction of the matrix 
+      * \f$ Q \f$ in the Hessenberg decomposition from the packed data.
       *
-      * \sa packedMatrix()
+      * \sa packedMatrix(), \ref Householder_Module "Householder module"
       */
     const CoeffVectorType& householderCoefficients() const { return m_hCoeffs; }
 
-    /** \returns a const reference to the internal representation of the decomposition.
+    /** \brief Returns the internal representation of the decomposition 
+      *
+      *	\returns a const reference to a matrix with the internal representation
+      *	         of the decomposition.
+      *
+      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
+      * or the member function compute(const MatrixType&) has been called
+      * before to compute the Hessenberg decomposition of a matrix.
       *
       * The returned matrix contains the following information:
       *  - the upper part and lower sub-diagonal represent the Hessenberg matrix H
       *  - the rest of the lower part contains the Householder vectors that, combined with
       *    Householder coefficients returned by householderCoefficients(),
-      *    allows to reconstruct the matrix Q as follow:
-      *       Q = H_{N-1} ... H_1 H_0
-      *    where the matrices H are the Householder transformation:
-      *       H_i = (I - h_i * v_i * v_i')
-      *    where h_i == householderCoefficients()[i] and v_i is a Householder vector:
-      *       v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ]
+      *    allows to reconstruct the matrix Q as 
+      *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
+      *    Here, the matrices \f$ H_i \f$ are the Householder transformations 
+      *       \f$ H_i = (I - h_i v_i v_i^T) \f$
+      *    where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and 
+      *    \f$ v_i \f$ is the Householder vector defined by
+      *       \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
+      *    with M the matrix returned by this function.
       *
       * See LAPACK for further details on this packed storage.
+      *
+      * Example: \include HessenbergDecomposition_packedMatrix.cpp
+      * Output: \verbinclude HessenbergDecomposition_packedMatrix.out
+      *
+      * \sa householderCoefficients()
       */
     const MatrixType& packedMatrix(void) const { return m_matrix; }
 
+    /** \brief Reconstructs the orthogonal matrix Q in the decomposition 
+      *
+      * \returns the matrix Q
+      *
+      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
+      * or the member function compute(const MatrixType&) has been called
+      * before to compute the Hessenberg decomposition of a matrix.
+      *
+      * This function reconstructs the matrix Q from the Householder
+      * coefficients and the packed matrix stored internally. This
+      * reconstruction requires \f$ 4n^3 / 3 \f$ flops.
+      *
+      * \sa matrixH() for an example
+      */
     MatrixType matrixQ() const;
+
+    /** \brief Constructs the Hessenberg matrix H in the decomposition
+      *
+      * \returns the matrix H
+      *
+      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
+      * or the member function compute(const MatrixType&) has been called
+      * before to compute the Hessenberg decomposition of a matrix.
+      *
+      * This function copies the matrix H from internal data. The upper part
+      * (including the subdiagonal) of the packed matrix as returned by
+      * packedMatrix() contains the matrix H. This function copies those
+      * entries in a newly created matrix and sets the remaining entries to
+      * zero. It may sometimes be sufficient to directly use the packed matrix
+      * instead of creating a new one.
+      *
+      * Example: \include HessenbergDecomposition_matrixH.cpp
+      * Output: \verbinclude HessenbergDecomposition_matrixH.out
+      *
+      * \sa matrixQ(), packedMatrix()
+      */
     MatrixType matrixH() const;
 
   private:
 
     static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs);
+    typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
+    typedef typename NumTraits<Scalar>::Real RealScalar;
 
   protected:
     MatrixType m_matrix;
@@ -134,7 +245,7 @@ template<typename _MatrixType> class HessenbergDecomposition
   *
   * The result is written in the lower triangular part of \a matA.
   *
-  * Implemented from Golub's "Matrix Computations", algorithm 8.3.1.
+  * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1.
   *
   * \sa packedMatrix()
   */
@@ -167,7 +278,6 @@ void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVector
   }
 }
 
-/** reconstructs and returns the matrix Q */
 template<typename MatrixType>
 typename HessenbergDecomposition<MatrixType>::MatrixType
 HessenbergDecomposition<MatrixType>::matrixQ() const
@@ -185,11 +295,6 @@ HessenbergDecomposition<MatrixType>::matrixQ() const
 
 #endif // EIGEN_HIDE_HEAVY_CODE
 
-/** constructs and returns the matrix H.
-  * Note that the matrix H is equivalent to the upper part of the packed matrix
-  * (including the lower sub-diagonal). Therefore, it might be often sufficient
-  * to directly use the packed matrix instead of creating a new one.
-  */
 template<typename MatrixType>
 typename HessenbergDecomposition<MatrixType>::MatrixType
 HessenbergDecomposition<MatrixType>::matrixH() const
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
new file mode 100644
index 000000000..a66ab8ace
--- /dev/null
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -0,0 +1,427 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_REAL_SCHUR_H
+#define EIGEN_REAL_SCHUR_H
+
+#include "./HessenbergDecomposition.h"
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+  * \nonstableyet
+  *
+  * \class RealSchur
+  *
+  * \brief Performs a real Schur decomposition of a square matrix
+  *
+  * \tparam _MatrixType the type of the matrix of which we are computing the
+  * real Schur decomposition; this is expected to be an instantiation of the
+  * Matrix class template.
+  *
+  * Given a real square matrix A, this class computes the real Schur
+  * decomposition: \f$ A = U T U^T \f$ where U is a real orthogonal matrix and
+  * T is a real quasi-triangular matrix. An orthogonal matrix is a matrix whose
+  * inverse is equal to its transpose, \f$ U^{-1} = U^T \f$. A quasi-triangular
+  * matrix is a block-triangular matrix whose diagonal consists of 1-by-1
+  * blocks and 2-by-2 blocks with complex eigenvalues. The eigenvalues of the
+  * blocks on the diagonal of T are the same as the eigenvalues of the matrix
+  * A, and thus the real Schur decomposition is used in EigenSolver to compute
+  * the eigendecomposition of a matrix.
+  *
+  * Call the function compute() to compute the real Schur decomposition of a
+  * given matrix. Alternatively, you can use the RealSchur(const MatrixType&)
+  * constructor which computes the real Schur decomposition at construction
+  * time. Once the decomposition is computed, you can use the matrixU() and
+  * matrixT() functions to retrieve the matrices U and T in the decomposition.
+  *
+  * The documentation of RealSchur(const MatrixType&) contains an example of
+  * the typical use of this class.
+  *
+  * \note The implementation is adapted from
+  * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
+  * Their code is based on EISPACK.
+  *
+  * \sa class ComplexSchur, class EigenSolver, class ComplexEigenSolver
+  */
+template<typename _MatrixType> class RealSchur
+{
+  public:
+    typedef _MatrixType MatrixType;
+    enum {
+      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+      Options = MatrixType::Options,
+      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
+    };
+    typedef typename MatrixType::Scalar Scalar;
+    typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
+    typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
+
+    /** \brief Default constructor.
+      *
+      * \param [in] size  Positive integer, size of the matrix whose Schur decomposition will be computed.
+      *
+      * The default constructor is useful in cases in which the user intends to
+      * perform decompositions via compute().  The \p size parameter is only
+      * used as a hint. It is not an error to give a wrong \p size, but it may
+      * impair performance.
+      *
+      * \sa compute() for an example.
+      */
+    RealSchur(int size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
+            : m_matT(size, size),
+              m_matU(size, size), 
+              m_isInitialized(false)
+    { }
+
+    /** \brief Constructor; computes real Schur decomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
+      *
+      * This constructor calls compute() to compute the Schur decomposition.
+      *
+      * Example: \include RealSchur_RealSchur_MatrixType.cpp
+      * Output: \verbinclude RealSchur_RealSchur_MatrixType.out
+      */
+    RealSchur(const MatrixType& matrix)
+            : m_matT(matrix.rows(),matrix.cols()),
+              m_matU(matrix.rows(),matrix.cols()),
+              m_isInitialized(false)
+    {
+      compute(matrix);
+    }
+
+    /** \brief Returns the orthogonal matrix in the Schur decomposition. 
+      *
+      * \returns A const reference to the matrix U.
+      *
+      * \pre Either the constructor RealSchur(const MatrixType&) or the member
+      * function compute(const MatrixType&) has been called before to compute
+      * the Schur decomposition of a matrix.
+      *
+      * \sa RealSchur(const MatrixType&) for an example
+      */
+    const MatrixType& matrixU() const
+    {
+      ei_assert(m_isInitialized && "RealSchur is not initialized.");
+      return m_matU;
+    }
+
+    /** \brief Returns the quasi-triangular matrix in the Schur decomposition. 
+      *
+      * \returns A const reference to the matrix T.
+      *
+      * \pre Either the constructor RealSchur(const MatrixType&) or the member
+      * function compute(const MatrixType&) has been called before to compute
+      * the Schur decomposition of a matrix.
+      *
+      * \sa RealSchur(const MatrixType&) for an example
+      */
+    const MatrixType& matrixT() const
+    {
+      ei_assert(m_isInitialized && "RealSchur is not initialized.");
+      return m_matT;
+    }
+  
+    /** \brief Computes Schur decomposition of given matrix. 
+      * 
+      * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
+      *
+      * The Schur decomposition is computed by first reducing the matrix to
+      * Hessenberg form using the class HessenbergDecomposition. The Hessenberg
+      * matrix is then reduced to triangular form by performing Francis QR
+      * iterations with implicit double shift. The cost of computing the Schur
+      * decomposition depends on the number of iterations; as a rough guide, it
+      * may be taken to be \f$25n^3\f$ flops.
+      *
+      * Example: \include RealSchur_compute.cpp
+      * Output: \verbinclude RealSchur_compute.out
+      */
+    void compute(const MatrixType& matrix);
+
+  private:
+    
+    MatrixType m_matT;
+    MatrixType m_matU;
+    bool m_isInitialized;
+
+    typedef Matrix<Scalar,3,1> Vector3s;
+
+    Scalar computeNormOfT();
+    int findSmallSubdiagEntry(int iu, Scalar norm);
+    void splitOffTwoRows(int iu, Scalar exshift);
+    void computeShift(int iu, int iter, Scalar& exshift, Vector3s& shiftInfo);
+    void initFrancisQRStep(int il, int iu, const Vector3s& shiftInfo, int& im, Vector3s& firstHouseholderVector);
+    void performFrancisQRStep(int il, int im, int iu, const Vector3s& firstHouseholderVector, Scalar* workspace);
+};
+
+
+template<typename MatrixType>
+void RealSchur<MatrixType>::compute(const MatrixType& matrix)
+{
+  assert(matrix.cols() == matrix.rows());
+
+  // Step 1. Reduce to Hessenberg form
+  // TODO skip Q if skipU = true
+  HessenbergDecomposition<MatrixType> hess(matrix);
+  m_matT = hess.matrixH();
+  m_matU = hess.matrixQ();
+
+  // Step 2. Reduce to real Schur form  
+  typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
+  ColumnVectorType workspaceVector(m_matU.cols());
+  Scalar* workspace = &workspaceVector.coeffRef(0);
+
+  // The matrix m_matT is divided in three parts. 
+  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 
+  // Rows il,...,iu is the part we are working on (the active window).
+  // Rows iu+1,...,end are already brought in triangular form.
+  int iu = m_matU.cols() - 1;
+  int iter = 0; // iteration count
+  Scalar exshift = 0.0; // sum of exceptional shifts
+  Scalar norm = computeNormOfT();
+
+  while (iu >= 0)
+  {
+    int il = findSmallSubdiagEntry(iu, norm);
+
+    // Check for convergence
+    if (il == iu) // One root found
+    {
+      m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
+      if (iu > 0) 
+	m_matT.coeffRef(iu, iu-1) = Scalar(0);
+      iu--;
+      iter = 0;
+    }
+    else if (il == iu-1) // Two roots found
+    {
+      splitOffTwoRows(iu, exshift);
+      iu -= 2;
+      iter = 0;
+    }
+    else // No convergence yet
+    {
+      Vector3s firstHouseholderVector, shiftInfo;
+      computeShift(iu, iter, exshift, shiftInfo);
+      iter = iter + 1;   // (Could check iteration count here.)
+      int im;
+      initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
+      performFrancisQRStep(il, im, iu, firstHouseholderVector, workspace);
+    }
+  } 
+
+  m_isInitialized = true;
+}
+
+/** \internal Computes and returns vector L1 norm of T */
+template<typename MatrixType>
+inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
+{
+  const int size = m_matU.cols();
+  // FIXME to be efficient the following would requires a triangular reduxion code
+  // Scalar norm = m_matT.upper().cwiseAbs().sum() 
+  //               + m_matT.corner(BottomLeft,size-1,size-1).diagonal().cwiseAbs().sum();
+  Scalar norm = 0.0;
+  for (int j = 0; j < size; ++j)
+    norm += m_matT.row(j).segment(std::max(j-1,0), size-std::max(j-1,0)).cwiseAbs().sum();
+  return norm;
+}
+
+/** \internal Look for single small sub-diagonal element and returns its index */
+template<typename MatrixType>
+inline int RealSchur<MatrixType>::findSmallSubdiagEntry(int iu, Scalar norm)
+{
+  int res = iu;
+  while (res > 0)
+  {
+    Scalar s = ei_abs(m_matT.coeff(res-1,res-1)) + ei_abs(m_matT.coeff(res,res));
+    if (s == 0.0)
+      s = norm;
+    if (ei_abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
+      break;
+    res--;
+  }
+  return res;
+}
+
+/** \internal Update T given that rows iu-1 and iu decouple from the rest. */
+template<typename MatrixType>
+inline void RealSchur<MatrixType>::splitOffTwoRows(int iu, Scalar exshift)
+{
+  const int size = m_matU.cols();
+
+  // The eigenvalues of the 2x2 matrix [a b; c d] are 
+  // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
+  Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
+  Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);   // q = tr^2 / 4 - det = discr/4
+  m_matT.coeffRef(iu,iu) += exshift;
+  m_matT.coeffRef(iu-1,iu-1) += exshift;
+
+  if (q >= 0) // Two real eigenvalues
+  {
+    Scalar z = ei_sqrt(ei_abs(q));
+    PlanarRotation<Scalar> rot;
+    if (p >= 0)
+      rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
+    else
+      rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
+
+    m_matT.block(0, iu-1, size, size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
+    m_matT.block(0, 0, iu+1, size).applyOnTheRight(iu-1, iu, rot);
+    m_matT.coeffRef(iu, iu-1) = Scalar(0); 
+    m_matU.applyOnTheRight(iu-1, iu, rot);
+  }
+
+  if (iu > 1) 
+    m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
+}
+
+/** \internal Form shift in shiftInfo, and update exshift if an exceptional shift is performed. */
+template<typename MatrixType>
+inline void RealSchur<MatrixType>::computeShift(int iu, int iter, Scalar& exshift, Vector3s& shiftInfo)
+{
+  shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
+  shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
+  shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
+
+  // Wilkinson's original ad hoc shift
+  if (iter == 10)
+  {
+    exshift += shiftInfo.coeff(0);
+    for (int i = 0; i <= iu; ++i)
+      m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
+    Scalar s = ei_abs(m_matT.coeff(iu,iu-1)) + ei_abs(m_matT.coeff(iu-1,iu-2));
+    shiftInfo.coeffRef(0) = Scalar(0.75) * s;
+    shiftInfo.coeffRef(1) = Scalar(0.75) * s;
+    shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
+  }
+
+  // MATLAB's new ad hoc shift
+  if (iter == 30)
+  {
+    Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
+    s = s * s + shiftInfo.coeff(2);
+    if (s > 0)
+    {
+      s = ei_sqrt(s);
+      if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
+        s = -s;
+      s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
+      s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
+      exshift += s;
+      for (int i = 0; i <= iu; ++i)
+        m_matT.coeffRef(i,i) -= s;
+      shiftInfo.setConstant(Scalar(0.964));
+    }
+  }
+}
+
+/** \internal Compute index im at which Francis QR step starts and the first Householder vector. */
+template<typename MatrixType>
+inline void RealSchur<MatrixType>::initFrancisQRStep(int il, int iu, const Vector3s& shiftInfo, int& im, Vector3s& firstHouseholderVector)
+{
+  Vector3s& v = firstHouseholderVector; // alias to save typing
+
+  for (im = iu-2; im >= il; --im)
+  {
+    const Scalar Tmm = m_matT.coeff(im,im);
+    const Scalar r = shiftInfo.coeff(0) - Tmm;
+    const Scalar s = shiftInfo.coeff(1) - Tmm;
+    v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
+    v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
+    v.coeffRef(2) = m_matT.coeff(im+2,im+1);
+    if (im == il) {
+      break;
+    }
+    const Scalar lhs = m_matT.coeff(im,im-1) * (ei_abs(v.coeff(1)) + ei_abs(v.coeff(2)));
+    const Scalar rhs = v.coeff(0) * (ei_abs(m_matT.coeff(im-1,im-1)) + ei_abs(Tmm) + ei_abs(m_matT.coeff(im+1,im+1)));
+    if (ei_abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
+    {
+      break;
+    }
+  }
+}
+
+/** \internal Perform a Francis QR step involving rows il:iu and columns im:iu. */
+template<typename MatrixType>
+inline void RealSchur<MatrixType>::performFrancisQRStep(int il, int im, int iu, const Vector3s& firstHouseholderVector, Scalar* workspace)
+{
+  assert(im >= il);
+  assert(im <= iu-2);
+
+  const int size = m_matU.cols();
+
+  for (int k = im; k <= iu-2; ++k)
+  {
+    bool firstIteration = (k == im);
+
+    Vector3s v;
+    if (firstIteration)
+      v = firstHouseholderVector;
+    else
+      v = m_matT.template block<3,1>(k,k-1);
+
+    Scalar tau, beta;
+    Matrix<Scalar, 2, 1> ess;
+    v.makeHouseholder(ess, tau, beta);
+    
+    if (beta != Scalar(0)) // if v is not zero
+    {
+      if (firstIteration && k > il)
+        m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
+      else if (!firstIteration)
+        m_matT.coeffRef(k,k-1) = beta;
+
+      // These Householder transformations form the O(n^3) part of the algorithm
+      m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
+      m_matT.block(0, k, std::min(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
+      m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
+    }
+  }
+
+  Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
+  Scalar tau, beta;
+  Matrix<Scalar, 1, 1> ess;
+  v.makeHouseholder(ess, tau, beta);
+
+  if (beta != Scalar(0)) // if v is not zero
+  {
+    m_matT.coeffRef(iu-1, iu-2) = beta;
+    m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
+    m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
+    m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
+  }
+
+  // clean up pollution due to round-off errors
+  for (int i = im+2; i <= iu; ++i)
+  {
+    m_matT.coeffRef(i,i-2) = Scalar(0);
+    if (i > im+2)
+      m_matT.coeffRef(i,i-3) = Scalar(0);
+  }
+}
+
+#endif // EIGEN_REAL_SCHUR_H
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index 63b924926..5062fd0ce 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -439,8 +439,7 @@ struct ei_solve_retval<PartialPivLU<_MatrixType>, Rhs>
     * Step 3: replace c by the solution x to Ux = c.
     */
 
-    const int size = dec().matrixLU().rows();
-    ei_assert(rhs().rows() == size);
+    ei_assert(rhs().rows() == dec().matrixLU().rows());
 
     // Step 1
     dst = dec().permutationP() * rhs();
diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h
index fbd035ce4..a9ef2d3a4 100644
--- a/Eigen/src/Sparse/CholmodSupport.h
+++ b/Eigen/src/Sparse/CholmodSupport.h
@@ -65,12 +65,12 @@ cholmod_sparse SparseMatrixBase<Derived>::asCholmodMatrix()
   res.p       = derived()._outerIndexPtr();
   res.i       = derived()._innerIndexPtr();
   res.x       = derived()._valuePtr();
-  res.xtype = CHOLMOD_REAL;
-  res.itype = CHOLMOD_INT;
-  res.sorted = 1;
-  res.packed = 1;
-  res.dtype = 0;
-  res.stype = -1;
+  res.xtype   = CHOLMOD_REAL;
+  res.itype   = CHOLMOD_INT;
+  res.sorted  = 1;
+  res.packed  = 1;
+  res.dtype   = 0;
+  res.stype   = -1;
 
   ei_cholmod_configure_matrix<Scalar>(res);
 
@@ -84,7 +84,7 @@ cholmod_sparse SparseMatrixBase<Derived>::asCholmodMatrix()
       res.stype = 0;
   }
   else
-    res.stype = 0;
+    res.stype = -1; // by default we consider the lower part
 
   return res;
 }
@@ -177,21 +177,21 @@ void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
   }
 
   cholmod_sparse A = const_cast<MatrixType&>(a).asCholmodMatrix();
-  m_cholmod.supernodal = CHOLMOD_AUTO;
+//   m_cholmod.supernodal = CHOLMOD_AUTO;
   // TODO
-  if (m_flags&IncompleteFactorization)
-  {
-    m_cholmod.nmethods = 1;
-    m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
-    m_cholmod.postorder = 0;
-  }
-  else
-  {
-    m_cholmod.nmethods = 1;
-    m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
-    m_cholmod.postorder = 0;
-  }
-  m_cholmod.final_ll = 1;
+//   if (m_flags&IncompleteFactorization)
+//   {
+//     m_cholmod.nmethods = 1;
+//     m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
+//     m_cholmod.postorder = 0;
+//   }
+//   else
+//   {
+//     m_cholmod.nmethods = 1;
+//     m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
+//     m_cholmod.postorder = 0;
+//   }
+//   m_cholmod.final_ll = 1;
   m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
   cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
 
diff --git a/Eigen/src/Sparse/SparseDot.h b/Eigen/src/Sparse/SparseDot.h
index 2baf09b0e..769bf8832 100644
--- a/Eigen/src/Sparse/SparseDot.h
+++ b/Eigen/src/Sparse/SparseDot.h
@@ -38,7 +38,7 @@ SparseMatrixBase<Derived>::dot(const MatrixBase<OtherDerived>& other) const
 
   ei_assert(size() == other.size());
   ei_assert(other.size()>0 && "you are using a non initialized vector");
-  
+
   typename Derived::InnerIterator i(derived(),0);
   Scalar res = 0;
   while (i)
@@ -59,9 +59,9 @@ SparseMatrixBase<Derived>::dot(const SparseMatrixBase<OtherDerived>& other) cons
   EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived)
   EIGEN_STATIC_ASSERT((ei_is_same_type<Scalar, typename OtherDerived::Scalar>::ret),
     YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
-  
+
   ei_assert(size() == other.size());
-  
+
   typename Derived::InnerIterator i(derived(),0);
   typename OtherDerived::InnerIterator j(other.derived(),0);
   Scalar res = 0;
@@ -84,7 +84,7 @@ template<typename Derived>
 inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
 SparseMatrixBase<Derived>::squaredNorm() const
 {
-  return ei_real((*this).cwise().abs2().sum());
+  return ei_real((*this).cwiseAbs2().sum());
 }
 
 template<typename Derived>
diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h
index f2a9d19dc..0caa8cbed 100644
--- a/Eigen/src/Sparse/TaucsSupport.h
+++ b/Eigen/src/Sparse/TaucsSupport.h
@@ -80,6 +80,7 @@ class SparseLLT<MatrixType,Taucs> : public SparseLLT<MatrixType>
     typedef SparseLLT<MatrixType> Base;
     typedef typename Base::Scalar Scalar;
     typedef typename Base::RealScalar RealScalar;
+    typedef typename Base::CholMatrixType CholMatrixType;
     using Base::MatrixLIsDirty;
     using Base::SupernodalFactorIsDirty;
     using Base::m_flags;
@@ -88,12 +89,12 @@ class SparseLLT<MatrixType,Taucs> : public SparseLLT<MatrixType>
 
   public:
 
-    SparseLLT(int flags = 0)
+    SparseLLT(int flags = SupernodalMultifrontal)
       : Base(flags), m_taucsSupernodalFactor(0)
     {
     }
 
-    SparseLLT(const MatrixType& matrix, int flags = 0)
+    SparseLLT(const MatrixType& matrix, int flags = SupernodalMultifrontal)
       : Base(flags), m_taucsSupernodalFactor(0)
     {
       compute(matrix);
@@ -105,7 +106,7 @@ class SparseLLT<MatrixType,Taucs> : public SparseLLT<MatrixType>
         taucs_supernodal_factor_free(m_taucsSupernodalFactor);
     }
 
-    inline const typename Base::CholMatrixType& matrixL(void) const;
+    inline const CholMatrixType& matrixL() const;
 
     template<typename Derived>
     void solveInPlace(MatrixBase<Derived> &b) const;
@@ -156,7 +157,7 @@ void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a)
 }
 
 template<typename MatrixType>
-inline const typename SparseLLT<MatrixType>::CholMatrixType&
+inline const typename SparseLLT<MatrixType,Taucs>::CholMatrixType&
 SparseLLT<MatrixType,Taucs>::matrixL() const
 {
   if (m_status & MatrixLIsDirty)
diff --git a/blas/CMakeLists.txt b/blas/CMakeLists.txt
index e71076f9d..d1f5b5676 100644
--- a/blas/CMakeLists.txt
+++ b/blas/CMakeLists.txt
@@ -4,9 +4,9 @@ add_custom_target(blas)
 
 set(EigenBlas_SRCS single.cpp double.cpp complex_single.cpp complex_double.cpp xerbla.cpp)
 
-add_library(eigen_blas ${EigenBlas_SRCS})
-# add_library(eigen_blas SHARED ${EigenBlas_SRCS})
-add_dependencies(blas eigen_blas)
+add_library(eigen_blas_static ${EigenBlas_SRCS})
+add_library(eigen_blas SHARED ${EigenBlas_SRCS})
+add_dependencies(blas eigen_blas eigen_blas_static)
 
 install(TARGETS eigen_blas
         RUNTIME DESTINATION bin
diff --git a/blas/level3_impl.h b/blas/level3_impl.h
index 6a0e64392..a4b04ce1c 100644
--- a/blas/level3_impl.h
+++ b/blas/level3_impl.h
@@ -27,7 +27,7 @@
 int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
 {
 //   std::cerr << "in gemm " << *opa << " " << *opb << " " << *m << " " << *n << " " << *k << " " << *lda << " " << *ldb << " " << *ldc << " " << *palpha << " " << *pbeta << "\n";
-  typedef void (*functype)(int, int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar);
+  typedef void (*functype)(int, int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar, Eigen::GemmParallelInfo<Scalar>*);
   static functype func[12];
 
   static bool init = false;
@@ -67,7 +67,7 @@ int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScal
     else
       matrix(c, *m, *n, *ldc) *= beta;
 
-  func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha);
+  func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha, 0);
   return 0;
 }
 
diff --git a/doc/eigendoxy.css b/doc/eigendoxy.css
index 77d062837..c1973289b 100644
--- a/doc/eigendoxy.css
+++ b/doc/eigendoxy.css
@@ -477,3 +477,8 @@ DIV.eimainmenu {
 /*   border-top:     solid; */
 /*   border-bottom:  solid; */
 }
+
+/* center version number on main page */
+H3.version { 
+  text-align:     center;
+}
diff --git a/doc/snippets/ComplexEigenSolver_compute.cpp b/doc/snippets/ComplexEigenSolver_compute.cpp
new file mode 100644
index 000000000..11d6bd399
--- /dev/null
+++ b/doc/snippets/ComplexEigenSolver_compute.cpp
@@ -0,0 +1,16 @@
+MatrixXcf A = MatrixXcf::Random(4,4);
+cout << "Here is a random 4x4 matrix, A:" << endl << A << endl << endl;
+
+ComplexEigenSolver<MatrixXcf> ces;
+ces.compute(A);
+cout << "The eigenvalues of A are:" << endl << ces.eigenvalues() << endl;
+cout << "The matrix of eigenvectors, V, is:" << endl << ces.eigenvectors() << endl << endl;
+
+complex<float> lambda = ces.eigenvalues()[0];
+cout << "Consider the first eigenvalue, lambda = " << lambda << endl;
+VectorXcf v = ces.eigenvectors().col(0);
+cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl;
+cout << "... and A * v = " << endl << A * v << endl << endl;
+
+cout << "Finally, V * D * V^(-1) = " << endl
+     << ces.eigenvectors() * ces.eigenvalues().asDiagonal() * ces.eigenvectors().inverse() << endl;
diff --git a/doc/snippets/ComplexEigenSolver_eigenvalues.cpp b/doc/snippets/ComplexEigenSolver_eigenvalues.cpp
new file mode 100644
index 000000000..1afa8b086
--- /dev/null
+++ b/doc/snippets/ComplexEigenSolver_eigenvalues.cpp
@@ -0,0 +1,4 @@
+MatrixXcf ones = MatrixXcf::Ones(3,3);
+ComplexEigenSolver<MatrixXcf> ces(ones);
+cout << "The eigenvalues of the 3x3 matrix of ones are:" 
+     << endl << ces.eigenvalues() << endl;
diff --git a/doc/snippets/ComplexEigenSolver_eigenvectors.cpp b/doc/snippets/ComplexEigenSolver_eigenvectors.cpp
new file mode 100644
index 000000000..bb1c2ccf1
--- /dev/null
+++ b/doc/snippets/ComplexEigenSolver_eigenvectors.cpp
@@ -0,0 +1,4 @@
+MatrixXcf ones = MatrixXcf::Ones(3,3);
+ComplexEigenSolver<MatrixXcf> ces(ones);
+cout << "The first eigenvector of the 3x3 matrix of ones is:" 
+     << endl << ces.eigenvectors().col(1) << endl;
diff --git a/doc/snippets/ComplexSchur_compute.cpp b/doc/snippets/ComplexSchur_compute.cpp
new file mode 100644
index 000000000..3a5170101
--- /dev/null
+++ b/doc/snippets/ComplexSchur_compute.cpp
@@ -0,0 +1,6 @@
+MatrixXcf A = MatrixXcf::Random(4,4);
+ComplexSchur<MatrixXcf> schur(4);
+schur.compute(A);
+cout << "The matrix T in the decomposition of A is:" << endl << schur.matrixT() << endl;
+schur.compute(A.inverse());
+cout << "The matrix T in the decomposition of A^(-1) is:" << endl << schur.matrixT() << endl;
diff --git a/doc/snippets/ComplexSchur_matrixT.cpp b/doc/snippets/ComplexSchur_matrixT.cpp
new file mode 100644
index 000000000..8380571ac
--- /dev/null
+++ b/doc/snippets/ComplexSchur_matrixT.cpp
@@ -0,0 +1,4 @@
+MatrixXcf A = MatrixXcf::Random(4,4);
+cout << "Here is a random 4x4 matrix, A:" << endl << A << endl << endl;
+ComplexSchur<MatrixXcf> schurOfA(A, false); // false means do not compute U
+cout << "The triangular matrix T is:" << endl << schurOfA.matrixT() << endl;
diff --git a/doc/snippets/ComplexSchur_matrixU.cpp b/doc/snippets/ComplexSchur_matrixU.cpp
new file mode 100644
index 000000000..ba3d9c22e
--- /dev/null
+++ b/doc/snippets/ComplexSchur_matrixU.cpp
@@ -0,0 +1,4 @@
+MatrixXcf A = MatrixXcf::Random(4,4);
+cout << "Here is a random 4x4 matrix, A:" << endl << A << endl << endl;
+ComplexSchur<MatrixXcf> schurOfA(A);
+cout << "The unitary matrix U is:" << endl << schurOfA.matrixU() << endl;
diff --git a/doc/snippets/EigenSolver_EigenSolver_MatrixType.cpp b/doc/snippets/EigenSolver_EigenSolver_MatrixType.cpp
new file mode 100644
index 000000000..c1d9fa879
--- /dev/null
+++ b/doc/snippets/EigenSolver_EigenSolver_MatrixType.cpp
@@ -0,0 +1,16 @@
+MatrixXd A = MatrixXd::Random(6,6);
+cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl;
+
+EigenSolver<MatrixXd> es(A);
+cout << "The eigenvalues of A are:" << endl << es.eigenvalues() << endl;
+cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl;
+
+complex<double> lambda = es.eigenvalues()[0];
+cout << "Consider the first eigenvalue, lambda = " << lambda << endl;
+VectorXcd v = es.eigenvectors().col(0);
+cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl;
+cout << "... and A * v = " << endl << A.cast<complex<double> >() * v << endl << endl;
+
+MatrixXcd D = es.eigenvalues().asDiagonal();
+MatrixXcd V = es.eigenvectors();
+cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl;
diff --git a/doc/snippets/EigenSolver_compute.cpp b/doc/snippets/EigenSolver_compute.cpp
new file mode 100644
index 000000000..06138f608
--- /dev/null
+++ b/doc/snippets/EigenSolver_compute.cpp
@@ -0,0 +1,6 @@
+EigenSolver<MatrixXf> es;
+MatrixXf A = MatrixXf::Random(4,4);
+es.compute(A);
+cout << "The eigenvalues of A are: " << es.eigenvalues().transpose() << endl;
+es.compute(A + MatrixXf::Identity(4,4)); // re-use es to compute eigenvalues of A+I
+cout << "The eigenvalues of A+I are: " << es.eigenvalues().transpose() << endl;
diff --git a/doc/snippets/EigenSolver_eigenvalues.cpp b/doc/snippets/EigenSolver_eigenvalues.cpp
new file mode 100644
index 000000000..8d83ea982
--- /dev/null
+++ b/doc/snippets/EigenSolver_eigenvalues.cpp
@@ -0,0 +1,4 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+EigenSolver<MatrixXd> es(ones);
+cout << "The eigenvalues of the 3x3 matrix of ones are:" 
+     << endl << es.eigenvalues() << endl;
diff --git a/doc/snippets/EigenSolver_eigenvectors.cpp b/doc/snippets/EigenSolver_eigenvectors.cpp
new file mode 100644
index 000000000..0fad4dadb
--- /dev/null
+++ b/doc/snippets/EigenSolver_eigenvectors.cpp
@@ -0,0 +1,4 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+EigenSolver<MatrixXd> es(ones);
+cout << "The first eigenvector of the 3x3 matrix of ones is:" 
+     << endl << es.eigenvectors().col(1) << endl;
diff --git a/doc/snippets/EigenSolver_pseudoEigenvectors.cpp b/doc/snippets/EigenSolver_pseudoEigenvectors.cpp
new file mode 100644
index 000000000..85e2569df
--- /dev/null
+++ b/doc/snippets/EigenSolver_pseudoEigenvectors.cpp
@@ -0,0 +1,9 @@
+MatrixXd A = MatrixXd::Random(6,6);
+cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl;
+
+EigenSolver<MatrixXd> es(A);
+MatrixXd D = es.pseudoEigenvalueMatrix();
+MatrixXd V = es.pseudoEigenvectors();
+cout << "The pseudo-eigenvalue matrix D is:" << endl << D << endl;
+cout << "The pseudo-eigenvector matrix V is:" << endl << V << endl;
+cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl;
diff --git a/doc/snippets/HessenbergDecomposition_compute.cpp b/doc/snippets/HessenbergDecomposition_compute.cpp
new file mode 100644
index 000000000..50e37833a
--- /dev/null
+++ b/doc/snippets/HessenbergDecomposition_compute.cpp
@@ -0,0 +1,6 @@
+MatrixXcf A = MatrixXcf::Random(4,4);
+HessenbergDecomposition<MatrixXcf> hd(4);
+hd.compute(A);
+cout << "The matrix H in the decomposition of A is:" << endl << hd.matrixH() << endl;
+hd.compute(2*A); // re-use hd to compute and store decomposition of 2A
+cout << "The matrix H in the decomposition of 2A is:" << endl << hd.matrixH() << endl;
diff --git a/doc/snippets/HessenbergDecomposition_matrixH.cpp b/doc/snippets/HessenbergDecomposition_matrixH.cpp
new file mode 100644
index 000000000..af0136668
--- /dev/null
+++ b/doc/snippets/HessenbergDecomposition_matrixH.cpp
@@ -0,0 +1,8 @@
+Matrix4f A = MatrixXf::Random(4,4);
+cout << "Here is a random 4x4 matrix:" << endl << A << endl;
+HessenbergDecomposition<MatrixXf> hessOfA(A);
+MatrixXf H = hessOfA.matrixH();
+cout << "The Hessenberg matrix H is:" << endl << H << endl;
+MatrixXf Q = hessOfA.matrixQ();
+cout << "The orthogonal matrix Q is:" << endl << Q << endl;
+cout << "Q H Q^T is:" << endl << Q * H * Q.transpose() << endl;
diff --git a/doc/snippets/HessenbergDecomposition_packedMatrix.cpp b/doc/snippets/HessenbergDecomposition_packedMatrix.cpp
new file mode 100644
index 000000000..4fa5957e8
--- /dev/null
+++ b/doc/snippets/HessenbergDecomposition_packedMatrix.cpp
@@ -0,0 +1,9 @@
+Matrix4d A = Matrix4d::Random(4,4);
+cout << "Here is a random 4x4 matrix:" << endl << A << endl;
+HessenbergDecomposition<Matrix4d> hessOfA(A);
+Matrix4d pm = hessOfA.packedMatrix();
+cout << "The packed matrix M is:" << endl << pm << endl;
+cout << "The upper Hessenberg part corresponds to the matrix H, which is:" 
+     << endl << hessOfA.matrixH() << endl;
+Vector3d hc = hessOfA.householderCoefficients();
+cout << "The vector of Householder coefficients is:" << endl << hc << endl;
diff --git a/doc/snippets/RealSchur_RealSchur_MatrixType.cpp b/doc/snippets/RealSchur_RealSchur_MatrixType.cpp
new file mode 100644
index 000000000..a5530dcc8
--- /dev/null
+++ b/doc/snippets/RealSchur_RealSchur_MatrixType.cpp
@@ -0,0 +1,10 @@
+MatrixXd A = MatrixXd::Random(6,6);
+cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl;
+
+RealSchur<MatrixXd> schur(A);
+cout << "The orthogonal matrix U is:" << endl << schur.matrixU() << endl;
+cout << "The quasi-triangular matrix T is:" << endl << schur.matrixT() << endl << endl;
+
+MatrixXd U = schur.matrixU();
+MatrixXd T = schur.matrixT();
+cout << "U * T * U^T = " << endl << U * T * U.transpose() << endl;
diff --git a/doc/snippets/RealSchur_compute.cpp b/doc/snippets/RealSchur_compute.cpp
new file mode 100644
index 000000000..4dcfaf0f2
--- /dev/null
+++ b/doc/snippets/RealSchur_compute.cpp
@@ -0,0 +1,6 @@
+MatrixXf A = MatrixXf::Random(4,4);
+RealSchur<MatrixXf> schur(4);
+schur.compute(A);
+cout << "The matrix T in the decomposition of A is:" << endl << schur.matrixT() << endl;
+schur.compute(A.inverse());
+cout << "The matrix T in the decomposition of A^(-1) is:" << endl << schur.matrixT() << endl;
diff --git a/doc/unsupported_modules.dox b/doc/unsupported_modules.dox
index 503ecb8b3..2f2b2f6a8 100644
--- a/doc/unsupported_modules.dox
+++ b/doc/unsupported_modules.dox
@@ -44,6 +44,9 @@ namespace Eigen {
   * \defgroup NumericalDiff_Module */
 
 /** \ingroup Unsupported_modules
+  * \defgroup Polynomials_Module */
+
+/** \ingroup Unsupported_modules
   * \defgroup Skyline_Module */
 
 } // namespace Eigen
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 072f63e1d..f7d5ed8e9 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -138,7 +138,9 @@ ei_add_test(qr)
 ei_add_test(qr_colpivoting)
 ei_add_test(qr_fullpivoting)
 ei_add_test(upperbidiagonalization)
-ei_add_test(hessenberg " " "${GSL_LIBRARIES}")
+ei_add_test(hessenberg)
+ei_add_test(schur_real)
+ei_add_test(schur_complex)
 ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}")
 ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}")
 ei_add_test(eigensolver_complex)
diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp
index fa052a9ae..4e973c877 100644
--- a/test/eigensolver_complex.cpp
+++ b/test/eigensolver_complex.cpp
@@ -61,5 +61,6 @@ void test_eigensolver_complex()
     CALL_SUBTEST_1( eigensolver(Matrix4cf()) );
     CALL_SUBTEST_2( eigensolver(MatrixXcd(14,14)) );
     CALL_SUBTEST_3( eigensolver(Matrix<std::complex<float>, 1, 1>()) );
+    CALL_SUBTEST_4( eigensolver(Matrix3f()) );
   }
 }
diff --git a/test/gsl_helper.h b/test/gsl_helper.h
index f8f7a5d79..eefffe792 100644
--- a/test/gsl_helper.h
+++ b/test/gsl_helper.h
@@ -33,6 +33,7 @@
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_complex.h>
 #include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_poly.h>
 
 namespace Eigen {
 
@@ -69,6 +70,27 @@ template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex> struct
     gsl_eigen_gensymmv_free(w);
     free(a);
   }
+
+  template<class EIGEN_VECTOR, class EIGEN_ROOTS>
+  static void eigen_poly_solve(const EIGEN_VECTOR& poly, EIGEN_ROOTS& roots )
+  {
+    const int deg = poly.size()-1;
+    double *z = new double[2*deg];
+    double *a = new double[poly.size()];
+    for( int i=0; i<poly.size(); ++i ){ a[i] = poly[i]; }
+
+    gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (poly.size());
+    gsl_poly_complex_solve(a, poly.size(), w, z);
+    gsl_poly_complex_workspace_free (w);
+    for( int i=0; i<deg; ++i )
+    {
+      roots[i].real() = z[2*i];
+      roots[i].imag() = z[2*i+1];
+    }
+
+    delete[] a;
+    delete[] z;
+  }
 };
 
 template<typename Scalar> struct GslTraits<Scalar,true>
diff --git a/test/hessenberg.cpp b/test/hessenberg.cpp
index d917be357..cba9c8fda 100644
--- a/test/hessenberg.cpp
+++ b/test/hessenberg.cpp
@@ -2,6 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -28,19 +29,36 @@
 template<typename Scalar,int Size> void hessenberg(int size = Size)
 {
   typedef Matrix<Scalar,Size,Size> MatrixType;
-  MatrixType m = MatrixType::Random(size,size);
-  HessenbergDecomposition<MatrixType> hess(m);
 
-  VERIFY_IS_APPROX(m, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
+  // Test basic functionality: A = U H U* and H is Hessenberg
+  for(int counter = 0; counter < g_repeat; ++counter) {
+    MatrixType m = MatrixType::Random(size,size);
+    HessenbergDecomposition<MatrixType> hess(m);
+    VERIFY_IS_APPROX(m, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
+    MatrixType H = hess.matrixH();
+    for(int row = 2; row < size; ++row) {
+      for(int col = 0; col < row-1; ++col) {
+	VERIFY(H(row,col) == (typename MatrixType::Scalar)0);
+      }
+    }
+  }
+
+  // Test whether compute() and constructor returns same result
+  MatrixType A = MatrixType::Random(size, size);
+  HessenbergDecomposition<MatrixType> cs1;
+  cs1.compute(A);
+  HessenbergDecomposition<MatrixType> cs2(A);
+  VERIFY_IS_EQUAL(cs1.matrixQ(), cs2.matrixQ());
+  VERIFY_IS_EQUAL(cs1.matrixH(), cs2.matrixH());
+
+  // TODO: Add tests for packedMatrix() and householderCoefficients()
 }
 
 void test_hessenberg()
 {
-  for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST_1(( hessenberg<std::complex<double>,1>() ));
-    CALL_SUBTEST_2(( hessenberg<std::complex<double>,2>() ));
-    CALL_SUBTEST_3(( hessenberg<std::complex<float>,4>() ));
-    CALL_SUBTEST_4(( hessenberg<float,Dynamic>(ei_random<int>(1,320)) ));
-    CALL_SUBTEST_5(( hessenberg<std::complex<double>,Dynamic>(ei_random<int>(1,320)) ));
-  }
+  CALL_SUBTEST_1(( hessenberg<std::complex<double>,1>() ));
+  CALL_SUBTEST_2(( hessenberg<std::complex<double>,2>() ));
+  CALL_SUBTEST_3(( hessenberg<std::complex<float>,4>() ));
+  CALL_SUBTEST_4(( hessenberg<float,Dynamic>(ei_random<int>(1,320)) ));
+  CALL_SUBTEST_5(( hessenberg<std::complex<double>,Dynamic>(ei_random<int>(1,320)) ));
 }
diff --git a/test/schur_complex.cpp b/test/schur_complex.cpp
new file mode 100644
index 000000000..3659a074c
--- /dev/null
+++ b/test/schur_complex.cpp
@@ -0,0 +1,67 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+#include <Eigen/Eigenvalues>
+
+template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTime)
+{
+  typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
+  typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType;
+
+  // Test basic functionality: T is triangular and A = U T U*
+  for(int counter = 0; counter < g_repeat; ++counter) {
+    MatrixType A = MatrixType::Random(size, size);
+    ComplexSchur<MatrixType> schurOfA(A);
+    ComplexMatrixType U = schurOfA.matrixU();
+    ComplexMatrixType T = schurOfA.matrixT();
+    for(int row = 1; row < size; ++row) {
+      for(int col = 0; col < row; ++col) {
+	VERIFY(T(row,col) == (typename MatrixType::Scalar)0);
+      }
+    }
+    VERIFY_IS_APPROX(A.template cast<ComplexScalar>(), U * T * U.adjoint());
+  }
+
+  // Test asserts when not initialized
+  ComplexSchur<MatrixType> csUninitialized;
+  VERIFY_RAISES_ASSERT(csUninitialized.matrixT());
+  VERIFY_RAISES_ASSERT(csUninitialized.matrixU());
+  
+  // Test whether compute() and constructor returns same result
+  MatrixType A = MatrixType::Random(size, size);
+  ComplexSchur<MatrixType> cs1;
+  cs1.compute(A);
+  ComplexSchur<MatrixType> cs2(A);
+  VERIFY_IS_EQUAL(cs1.matrixT(), cs2.matrixT());
+  VERIFY_IS_EQUAL(cs1.matrixU(), cs2.matrixU());
+}
+
+void test_schur_complex()
+{
+  CALL_SUBTEST_1(( schur<Matrix4cd>() ));
+  CALL_SUBTEST_2(( schur<MatrixXcf>(ei_random<int>(1,50)) ));
+  CALL_SUBTEST_3(( schur<Matrix<std::complex<float>, 1, 1> >() ));
+  CALL_SUBTEST_4(( schur<Matrix<float, 3, 3, Eigen::RowMajor> >() ));
+}
diff --git a/test/schur_real.cpp b/test/schur_real.cpp
new file mode 100644
index 000000000..c1747e2f5
--- /dev/null
+++ b/test/schur_real.cpp
@@ -0,0 +1,85 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+#include <Eigen/Eigenvalues>
+
+template<typename MatrixType> void verifyIsQuasiTriangular(const MatrixType& T)
+{
+  const int size = T.cols();
+  typedef typename MatrixType::Scalar Scalar;
+
+  // Check T is lower Hessenberg
+  for(int row = 2; row < size; ++row) {
+    for(int col = 0; col < row - 1; ++col) {
+      VERIFY(T(row,col) == Scalar(0));
+    }
+  }
+
+  // Check that any non-zero on the subdiagonal is followed by a zero and is
+  // part of a 2x2 diagonal block with imaginary eigenvalues.
+  for(int row = 1; row < size; ++row) {
+    if (T(row,row-1) != Scalar(0)) {
+      VERIFY(row == size-1 || T(row+1,row) == 0);
+      Scalar tr = T(row-1,row-1) + T(row,row);
+      Scalar det = T(row-1,row-1) * T(row,row) - T(row-1,row) * T(row,row-1);
+      VERIFY(4 * det > tr * tr);
+    }
+  }
+}
+
+template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTime)
+{
+  // Test basic functionality: T is quasi-triangular and A = U T U*
+  for(int counter = 0; counter < g_repeat; ++counter) {
+    MatrixType A = MatrixType::Random(size, size);
+    RealSchur<MatrixType> schurOfA(A);
+    MatrixType U = schurOfA.matrixU();
+    MatrixType T = schurOfA.matrixT();
+    std::cout << "T = \n" << T << "\n\n";
+    verifyIsQuasiTriangular(T);
+    VERIFY_IS_APPROX(A, U * T * U.transpose());
+  }
+
+  // Test asserts when not initialized
+  RealSchur<MatrixType> rsUninitialized;
+  VERIFY_RAISES_ASSERT(rsUninitialized.matrixT());
+  VERIFY_RAISES_ASSERT(rsUninitialized.matrixU());
+  
+  // Test whether compute() and constructor returns same result
+  MatrixType A = MatrixType::Random(size, size);
+  RealSchur<MatrixType> rs1;
+  rs1.compute(A);
+  RealSchur<MatrixType> rs2(A);
+  VERIFY_IS_EQUAL(rs1.matrixT(), rs2.matrixT());
+  VERIFY_IS_EQUAL(rs1.matrixU(), rs2.matrixU());
+}
+
+void test_schur_real()
+{
+  CALL_SUBTEST_1(( schur<Matrix4f>() ));
+  CALL_SUBTEST_2(( schur<MatrixXd>(ei_random<int>(1,50)) ));
+  CALL_SUBTEST_3(( schur<Matrix<float, 1, 1> >() ));
+  CALL_SUBTEST_4(( schur<Matrix<double, 3, 3, Eigen::RowMajor> >() ));
+}
diff --git a/test/sparse_vector.cpp b/test/sparse_vector.cpp
index 5c6dadc00..a5101398c 100644
--- a/test/sparse_vector.cpp
+++ b/test/sparse_vector.cpp
@@ -79,13 +79,15 @@ template<typename Scalar> void sparse_vector(int rows, int cols)
 
   VERIFY_IS_APPROX(v1*=s1, refV1*=s1);
   VERIFY_IS_APPROX(v1/=s1, refV1/=s1);
-  
+
   VERIFY_IS_APPROX(v1+=v2, refV1+=refV2);
   VERIFY_IS_APPROX(v1-=v2, refV1-=refV2);
 
   VERIFY_IS_APPROX(v1.dot(v2), refV1.dot(refV2));
   VERIFY_IS_APPROX(v1.dot(refV2), refV1.dot(refV2));
 
+  VERIFY_IS_APPROX(v1.squaredNorm(), refV1.squaredNorm());
+
 }
 
 void test_sparse_vector()
diff --git a/unsupported/Eigen/CMakeLists.txt b/unsupported/Eigen/CMakeLists.txt
index f9d79c51a..87cc4be1e 100644
--- a/unsupported/Eigen/CMakeLists.txt
+++ b/unsupported/Eigen/CMakeLists.txt
@@ -1,4 +1,4 @@
-set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3)
+set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3 Polynomials)
 
 install(FILES
   ${Eigen_HEADERS}
diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT
index 045c44611..2a1e21007 100644
--- a/unsupported/Eigen/FFT
+++ b/unsupported/Eigen/FFT
@@ -320,7 +320,7 @@ class FFT
         // if the vector is strided, then we need to copy it to a packed temporary
         Matrix<src_type,1,Dynamic> tmp;
         if ( resize_input ) {
-          size_t ncopy = min(src.size(),src.size() + resize_input);
+          size_t ncopy = std::min(src.size(),src.size() + resize_input);
           tmp.setZero(src.size() + resize_input);
           if ( realfft && HasFlag(HalfSpectrum) ) {
             // pad at the Nyquist bin
diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions
index 292357fda..e0bc4732c 100644
--- a/unsupported/Eigen/MatrixFunctions
+++ b/unsupported/Eigen/MatrixFunctions
@@ -40,6 +40,22 @@ namespace Eigen {
   * \brief This module aims to provide various methods for the computation of
   * matrix functions. 
   *
+  * To use this module, add 
+  * \code
+  * #include <unsupported/Eigen/MatrixFunctions>
+  * \endcode
+  * at the start of your source file.
+  *
+  * This module defines the following MatrixBase methods.
+  *  - \ref matrixbase_cos "MatrixBase::cos()", for computing the matrix cosine
+  *  - \ref matrixbase_cosh "MatrixBase::cosh()", for computing the matrix hyperbolic cosine
+  *  - \ref matrixbase_exp "MatrixBase::exp()", for computing the matrix exponential
+  *  - \ref matrixbase_matrixfunction "MatrixBase::matrixFunction()", for computing general matrix functions
+  *  - \ref matrixbase_sin "MatrixBase::sin()", for computing the matrix sine
+  *  - \ref matrixbase_sinh "MatrixBase::sinh()", for computing the matrix hyperbolic sine
+  *
+  * These methods are the main entry points to this module. 
+  *
   * %Matrix functions are defined as follows.  Suppose that \f$ f \f$
   * is an entire function (that is, a function on the complex plane
   * that is everywhere complex differentiable).  Then its Taylor
@@ -49,16 +65,205 @@ namespace Eigen {
   * function by the same series:
   * \f[ f(M) = f(0) + f'(0) M + \frac{f''(0)}{2} M^2 + \frac{f'''(0)}{3!} M^3 + \cdots \f]
   *
-  * \code
-  * #include <unsupported/Eigen/MatrixFunctions>
-  * \endcode
   */
 
 #include "src/MatrixFunctions/MatrixExponential.h"
 #include "src/MatrixFunctions/MatrixFunction.h"
 
-}
 
 
+/** 
+\page matrixbaseextra MatrixBase methods defined in the MatrixFunctions module
+\ingroup MatrixFunctions_Module
+
+The remainder of the page documents the following MatrixBase methods
+which are defined in the MatrixFunctions module.
+
+
+
+\section matrixbase_cos MatrixBase::cos()
+
+Compute the matrix cosine.
+
+\code
+const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
+\endcode
+
+\param[in]  M  a square matrix.
+\returns  expression representing \f$ \cos(M) \f$.
+
+This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos().
+
+\sa \ref matrixbase_sin "sin()" for an example.
+
+
+
+\section matrixbase_cosh MatrixBase::cosh()
+
+Compute the matrix hyberbolic cosine.
+
+\code
+const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
+\endcode
+
+\param[in]  M  a square matrix.
+\returns  expression representing \f$ \cosh(M) \f$
+
+This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cosh().
+
+\sa \ref matrixbase_sinh "sinh()" for an example.
+
+
+
+\section matrixbase_exp MatrixBase::exp()
+
+Compute the matrix exponential.
+
+\code
+const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
+\endcode
+
+\param[in]  M  matrix whose exponential is to be computed.
+\returns    expression representing the matrix exponential of \p M.
+
+The matrix exponential of \f$ M \f$ is defined by
+\f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f]
+The matrix exponential can be used to solve linear ordinary
+differential equations: the solution of \f$ y' = My \f$ with the
+initial condition \f$ y(0) = y_0 \f$ is given by
+\f$ y(t) = \exp(M) y_0 \f$.
+
+The cost of the computation is approximately \f$ 20 n^3 \f$ for
+matrices of size \f$ n \f$. The number 20 depends weakly on the
+norm of the matrix.
+
+The matrix exponential is computed using the scaling-and-squaring
+method combined with Pad&eacute; approximation. The matrix is first
+rescaled, then the exponential of the reduced matrix is computed
+approximant, and then the rescaling is undone by repeated
+squaring. The degree of the Pad&eacute; approximant is chosen such
+that the approximation error is less than the round-off
+error. However, errors may accumulate during the squaring phase.
+
+Details of the algorithm can be found in: Nicholas J. Higham, "The
+scaling and squaring method for the matrix exponential revisited,"
+<em>SIAM J. %Matrix Anal. Applic.</em>, <b>26</b>:1179&ndash;1193,
+2005.
+
+Example: The following program checks that
+\f[ \exp \left[ \begin{array}{ccc}
+      0 & \frac14\pi & 0 \\
+      -\frac14\pi & 0 & 0 \\
+      0 & 0 & 0
+    \end{array} \right] = \left[ \begin{array}{ccc}
+      \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
+      \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
+      0 & 0 & 1
+    \end{array} \right]. \f]
+This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
+the z-axis.
+
+\include MatrixExponential.cpp
+Output: \verbinclude MatrixExponential.out
+
+\note \p M has to be a matrix of \c float, \c double,
+\c complex<float> or \c complex<double> .
+
+
+
+\section matrixbase_matrixfunction MatrixBase::matrixFunction()
+
+Compute a matrix function.
+
+\code
+const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f) const
+\endcode
+
+\param[in]  M  argument of matrix function, should be a square matrix.
+\param[in]  f  an entire function; \c f(x,n) should compute the n-th
+derivative of f at x.
+\returns  expression representing \p f applied to \p M.
+
+Suppose that \p M is a matrix whose entries have type \c Scalar. 
+Then, the second argument, \p f, should be a function with prototype
+\code 
+ComplexScalar f(ComplexScalar, int) 
+\endcode
+where \c ComplexScalar = \c std::complex<Scalar> if \c Scalar is
+real (e.g., \c float or \c double) and \c ComplexScalar =
+\c Scalar if \c Scalar is complex. The return value of \c f(x,n)
+should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x.
+
+This routine uses the algorithm described in:
+Philip Davies and Nicholas J. Higham, 
+"A Schur-Parlett algorithm for computing matrix functions", 
+<em>SIAM J. %Matrix Anal. Applic.</em>, <b>25</b>:464&ndash;485, 2003.
+
+The actual work is done by the MatrixFunction class.
+
+Example: The following program checks that
+\f[ \exp \left[ \begin{array}{ccc} 
+      0 & \frac14\pi & 0 \\ 
+      -\frac14\pi & 0 & 0 \\
+      0 & 0 & 0 
+    \end{array} \right] = \left[ \begin{array}{ccc}
+      \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
+      \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
+      0 & 0 & 1
+    \end{array} \right]. \f]
+This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
+the z-axis. This is the same example as used in the documentation
+of \ref matrixbase_exp "exp()".
+
+\include MatrixFunction.cpp
+Output: \verbinclude MatrixFunction.out
+
+Note that the function \c expfn is defined for complex numbers 
+\c x, even though the matrix \c A is over the reals. Instead of
+\c expfn, we could also have used StdStemFunctions::exp:
+\code
+A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B);
+\endcode
+
+
+
+\section matrixbase_sin MatrixBase::sin()
+
+Compute the matrix sine.
+
+\code
+const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
+\endcode
+
+\param[in]  M  a square matrix.
+\returns  expression representing \f$ \sin(M) \f$.
+
+This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin().
+
+Example: \include MatrixSine.cpp
+Output: \verbinclude MatrixSine.out
+
+
+
+\section matrixbase_sinh const MatrixBase::sinh()
+
+Compute the matrix hyperbolic sine.
+
+\code
+MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
+\endcode
+
+\param[in]  M  a square matrix.
+\returns  expression representing \f$ \sinh(M) \f$
+
+This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sinh().
+
+Example: \include MatrixSinh.cpp
+Output: \verbinclude MatrixSinh.out
+
+*/
+
+}
+
 #endif // EIGEN_MATRIX_FUNCTIONS
 
diff --git a/unsupported/Eigen/Polynomials b/unsupported/Eigen/Polynomials
new file mode 100644
index 000000000..d69abb1be
--- /dev/null
+++ b/unsupported/Eigen/Polynomials
@@ -0,0 +1,137 @@
+#ifndef EIGEN_POLYNOMIALS_MODULE_H
+#define EIGEN_POLYNOMIALS_MODULE_H
+
+#include <Eigen/Core>
+
+#include <Eigen/src/Core/util/DisableMSVCWarnings.h>
+
+#include <Eigen/QR>
+
+// Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module
+#if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2)
+  #ifndef EIGEN_HIDE_HEAVY_CODE
+  #define EIGEN_HIDE_HEAVY_CODE
+  #endif
+#elif defined EIGEN_HIDE_HEAVY_CODE
+  #undef EIGEN_HIDE_HEAVY_CODE
+#endif
+
+namespace Eigen {
+
+/** \ingroup Unsupported_modules
+  * \defgroup Polynomials_Module Polynomials module
+  *
+  * \nonstableyet
+  *
+  * \brief This module provides a QR based polynomial solver.
+	*
+  * To use this module, add
+  * \code
+  * #include <unsupported/Eigen/Polynomials>
+  * \endcode
+	* at the start of your source file.
+  */
+
+#include "src/Polynomials/PolynomialUtils.h"
+#include "src/Polynomials/Companion.h"
+#include "src/Polynomials/PolynomialSolver.h"
+
+/**
+	\page polynomials Polynomials defines functions for dealing with polynomials
+	and a QR based polynomial solver.
+	\ingroup Polynomials_Module
+
+	The remainder of the page documents first the functions for evaluating, computing
+	polynomials, computing estimates about polynomials and next the QR based polynomial
+	solver.
+
+	\section polynomialUtils convenient functions to deal with polynomials
+	\subsection roots_to_monicPolynomial
+	The function
+	\code
+	void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
+	\endcode
+	computes the coefficients \f$ a_i \f$ of
+
+	\f$ p(x) = a_0 + a_{1}x + ... + a_{n-1}x^{n-1} + x^n \f$
+
+	where \f$ p \f$ is known through its roots i.e. \f$ p(x) = (x-r_1)(x-r_2)...(x-r_n) \f$.
+
+	\subsection poly_eval
+	The function
+	\code
+	T poly_eval( const Polynomials& poly, const T& x )
+	\endcode
+	evaluates a polynomial at a given point using stabilized H&ouml;rner method.
+
+	The following code: first computes the coefficients in the monomial basis of the monic polynomial that has the provided roots;
+	then, it evaluates the computed polynomial, using a stabilized H&ouml;rner method.
+
+	\include PolynomialUtils1.cpp
+  Output: \verbinclude PolynomialUtils1.out
+
+	\subsection Cauchy bounds
+	The function
+	\code
+	Real cauchy_max_bound( const Polynomial& poly )
+	\endcode
+	provides a maximum bound (the Cauchy one: \f$C(p)\f$) for the absolute value of a root of the given polynomial i.e.
+	\f$ \forall r_i \f$ root of \f$ p(x) = \sum_{k=0}^d a_k x^k \f$,
+	\f$ |r_i| \le C(p) = \sum_{k=0}^{d} \left | \frac{a_k}{a_d} \right | \f$
+	The leading coefficient \f$ p \f$: should be non zero \f$a_d \neq 0\f$.
+
+
+	The function
+	\code
+	Real cauchy_min_bound( const Polynomial& poly )
+	\endcode
+	provides a minimum bound (the Cauchy one: \f$c(p)\f$) for the absolute value of a non zero root of the given polynomial i.e.
+	\f$ \forall r_i \neq 0 \f$ root of \f$ p(x) = \sum_{k=0}^d a_k x^k \f$,
+	\f$ |r_i| \ge c(p) = \left( \sum_{k=0}^{d} \left | \frac{a_k}{a_0} \right | \right)^{-1} \f$
+
+
+
+
+	\section QR polynomial solver class
+	Computes the complex roots of a polynomial by computing the eigenvalues of the associated companion matrix with the QR algorithm.
+	
+	The roots of \f$ p(x) = a_0 + a_1 x + a_2 x^2 + a_{3} x^3 + x^4 \f$ are the eigenvalues of
+	\f$
+	\left [
+	\begin{array}{cccc}
+	0 & 0 &  0 & a_0 \\
+	1 & 0 &  0 & a_1 \\
+	0 & 1 &  0 & a_2 \\
+	0 & 0 &  1 & a_3
+	\end{array} \right ]
+	\f$
+
+	However, the QR algorithm is not guaranteed to converge when there are several eigenvalues with same modulus.
+
+	Therefore the current polynomial solver is guaranteed to provide a correct result only when the complex roots \f$r_1,r_2,...,r_d\f$ have distinct moduli i.e.
+	
+	\f$ \forall i,j \in [1;d],~ \| r_i \| \neq \| r_j \| \f$.
+
+	With 32bit (float) floating types this problem shows up frequently.
+  However, almost always, correct accuracy is reached even in these cases for 64bit
+  (double) floating types and small polynomial degree (<20).
+
+	\include PolynomialSolver1.cpp
+	
+	In the above example:
+	
+	-# a simple use of the polynomial solver is shown;
+	-# the accuracy problem with the QR algorithm is presented: a polynomial with almost conjugate roots is provided to the solver.
+	Those roots have almost same module therefore the QR algorithm failed to converge: the accuracy
+	of the last root is bad;
+	-# a simple way to circumvent the problem is shown: use doubles instead of floats.
+
+  Output: \verbinclude PolynomialSolver1.out
+*/
+
+} // namespace Eigen
+
+#include <Eigen/src/Core/util/EnableMSVCWarnings.h>
+
+#endif // EIGEN_POLYNOMIALS_MODULE_H
+/* vim: set filetype=cpp et sw=2 ts=2 ai: */
diff --git a/unsupported/Eigen/src/CMakeLists.txt b/unsupported/Eigen/src/CMakeLists.txt
index a0870d3af..035a84ca8 100644
--- a/unsupported/Eigen/src/CMakeLists.txt
+++ b/unsupported/Eigen/src/CMakeLists.txt
@@ -5,3 +5,4 @@ ADD_SUBDIRECTORY(MoreVectorization)
 # ADD_SUBDIRECTORY(FFT)
 # ADD_SUBDIRECTORY(Skyline)
 ADD_SUBDIRECTORY(MatrixFunctions)
+ADD_SUBDIRECTORY(Polynomials)
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
index 3ef91a7b2..d7eb1b2fe 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
@@ -330,56 +330,6 @@ struct ei_traits<MatrixExponentialReturnValue<Derived> >
   typedef typename Derived::PlainObject ReturnType;
 };
 
-/** \ingroup MatrixFunctions_Module
- *
- * \brief Compute the matrix exponential.
- *
- * \param[in]  M  matrix whose exponential is to be computed.
- * \returns    expression representing the matrix exponential of \p M.
- *
- * The matrix exponential of \f$ M \f$ is defined by
- * \f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f]
- * The matrix exponential can be used to solve linear ordinary
- * differential equations: the solution of \f$ y' = My \f$ with the
- * initial condition \f$ y(0) = y_0 \f$ is given by
- * \f$ y(t) = \exp(M) y_0 \f$.
- *
- * The cost of the computation is approximately \f$ 20 n^3 \f$ for
- * matrices of size \f$ n \f$. The number 20 depends weakly on the
- * norm of the matrix.
- *
- * The matrix exponential is computed using the scaling-and-squaring
- * method combined with Pad&eacute; approximation. The matrix is first
- * rescaled, then the exponential of the reduced matrix is computed
- * approximant, and then the rescaling is undone by repeated
- * squaring. The degree of the Pad&eacute; approximant is chosen such
- * that the approximation error is less than the round-off
- * error. However, errors may accumulate during the squaring phase.
- *
- * Details of the algorithm can be found in: Nicholas J. Higham, "The
- * scaling and squaring method for the matrix exponential revisited,"
- * <em>SIAM J. %Matrix Anal. Applic.</em>, <b>26</b>:1179&ndash;1193,
- * 2005.
- *
- * Example: The following program checks that
- * \f[ \exp \left[ \begin{array}{ccc}
- *       0 & \frac14\pi & 0 \\
- *       -\frac14\pi & 0 & 0 \\
- *       0 & 0 & 0
- *     \end{array} \right] = \left[ \begin{array}{ccc}
- *       \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
- *       \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
- *       0 & 0 & 1
- *     \end{array} \right]. \f]
- * This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
- * the z-axis.
- *
- * \include MatrixExponential.cpp
- * Output: \verbinclude MatrixExponential.out
- *
- * \note \p M has to be a matrix of \c float, \c double,
- * \c complex<float> or \c complex<double> .
- */
 template <typename Derived>
 const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const
 {
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h
index 72e42aa7c..270d23b8b 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h
@@ -536,56 +536,6 @@ struct ei_traits<MatrixFunctionReturnValue<Derived> >
 /********** MatrixBase methods **********/
 
 
-/** \ingroup MatrixFunctions_Module
-  *
-  * \brief Compute a matrix function.
-  *
-  * \param[in]  M  argument of matrix function, should be a square matrix.
-  * \param[in]  f  an entire function; \c f(x,n) should compute the n-th
-  * derivative of f at x.
-  * \returns  expression representing \p f applied to \p M.
-  *
-  * Suppose that \p M is a matrix whose entries have type \c Scalar. 
-  * Then, the second argument, \p f, should be a function with prototype
-  * \code 
-  * ComplexScalar f(ComplexScalar, int) 
-  * \endcode
-  * where \c ComplexScalar = \c std::complex<Scalar> if \c Scalar is
-  * real (e.g., \c float or \c double) and \c ComplexScalar =
-  * \c Scalar if \c Scalar is complex. The return value of \c f(x,n)
-  * should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x.
-  *
-  * This routine uses the algorithm described in:
-  * Philip Davies and Nicholas J. Higham, 
-  * "A Schur-Parlett algorithm for computing matrix functions", 
-  * <em>SIAM J. %Matrix Anal. Applic.</em>, <b>25</b>:464&ndash;485, 2003.
-  *
-  * The actual work is done by the MatrixFunction class.
-  *
-  * Example: The following program checks that
-  * \f[ \exp \left[ \begin{array}{ccc} 
-  *       0 & \frac14\pi & 0 \\ 
-  *       -\frac14\pi & 0 & 0 \\
-  *       0 & 0 & 0 
-  *     \end{array} \right] = \left[ \begin{array}{ccc}
-  *       \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\
-  *       \frac12\sqrt2 & \frac12\sqrt2 & 0 \\
-  *       0 & 0 & 1
-  *     \end{array} \right]. \f]
-  * This corresponds to a rotation of \f$ \frac14\pi \f$ radians around
-  * the z-axis. This is the same example as used in the documentation
-  * of MatrixBase::exp().
-  *
-  * \include MatrixFunction.cpp
-  * Output: \verbinclude MatrixFunction.out
-  *
-  * Note that the function \c expfn is defined for complex numbers 
-  * \c x, even though the matrix \c A is over the reals. Instead of
-  * \c expfn, we could also have used StdStemFunctions::exp:
-  * \code
-  * A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B);
-  * \endcode
-  */
 template <typename Derived>
 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f) const
 {
@@ -593,18 +543,6 @@ const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typ
   return MatrixFunctionReturnValue<Derived>(derived(), f);
 }
 
-/** \ingroup MatrixFunctions_Module
-  *
-  * \brief Compute the matrix sine.
-  *
-  * \param[in]  M  a square matrix.
-  * \returns  expression representing \f$ \sin(M) \f$.
-  * 
-  * This function calls matrixFunction() with StdStemFunctions::sin().
-  *
-  * \include MatrixSine.cpp
-  * Output: \verbinclude MatrixSine.out
-  */
 template <typename Derived>
 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
 {
@@ -613,17 +551,6 @@ const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const
   return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sin);
 }
 
-/** \ingroup MatrixFunctions_Module
-  *
-  * \brief Compute the matrix cosine.
-  *
-  * \param[in]  M  a square matrix.
-  * \returns  expression representing \f$ \cos(M) \f$.
-  * 
-  * This function calls matrixFunction() with StdStemFunctions::cos().
-  *
-  * \sa ei_matrix_sin() for an example.
-  */
 template <typename Derived>
 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
 {
@@ -632,18 +559,6 @@ const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const
   return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cos);
 }
 
-/** \ingroup MatrixFunctions_Module
-  *
-  * \brief Compute the matrix hyperbolic sine.
-  *
-  * \param[in]  M  a square matrix.
-  * \returns  expression representing \f$ \sinh(M) \f$
-  * 
-  * This function calls matrixFunction() with StdStemFunctions::sinh().
-  *
-  * \include MatrixSinh.cpp
-  * Output: \verbinclude MatrixSinh.out
-  */
 template <typename Derived>
 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
 {
@@ -652,17 +567,6 @@ const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const
   return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sinh);
 }
 
-/** \ingroup MatrixFunctions_Module
-  *
-  * \brief Compute the matrix hyberbolic cosine.
-  *
-  * \param[in]  M  a square matrix.
-  * \returns  expression representing \f$ \cosh(M) \f$
-  * 
-  * This function calls matrixFunction() with StdStemFunctions::cosh().
-  *
-  * \sa ei_matrix_sinh() for an example.
-  */
 template <typename Derived>
 const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const
 {
diff --git a/unsupported/Eigen/src/Polynomials/CMakeLists.txt b/unsupported/Eigen/src/Polynomials/CMakeLists.txt
new file mode 100644
index 000000000..51f13f3cb
--- /dev/null
+++ b/unsupported/Eigen/src/Polynomials/CMakeLists.txt
@@ -0,0 +1,6 @@
+FILE(GLOB Eigen_Polynomials_SRCS "*.h")
+
+INSTALL(FILES
+  ${Eigen_Polynomials_SRCS}
+  DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/Polynomials COMPONENT Devel
+  )
diff --git a/unsupported/Eigen/src/Polynomials/Companion.h b/unsupported/Eigen/src/Polynomials/Companion.h
new file mode 100644
index 000000000..7c9e9c518
--- /dev/null
+++ b/unsupported/Eigen/src/Polynomials/Companion.h
@@ -0,0 +1,281 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_COMPANION_H
+#define EIGEN_COMPANION_H
+
+// This file requires the user to include
+// * Eigen/Core
+// * Eigen/src/PolynomialSolver.h
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+
+template <typename T>
+T ei_radix(){ return 2; }
+
+template <typename T>
+T ei_radix2(){ return ei_radix<T>()*ei_radix<T>(); }
+
+template<int Size>
+struct ei_decrement_if_fixed_size
+{
+  enum {
+    ret = (Size == Dynamic) ? Dynamic : Size-1 };
+};
+
+#endif
+
+template< typename _Scalar, int _Deg >
+class ei_companion
+{
+  public:
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
+
+    enum {
+      Deg = _Deg,
+      Deg_1=ei_decrement_if_fixed_size<Deg>::ret
+    };
+
+    typedef _Scalar                                Scalar;
+    typedef typename NumTraits<Scalar>::Real       RealScalar;
+    typedef Matrix<Scalar, Deg, 1>                 RightColumn;
+    //typedef DiagonalMatrix< Scalar, Deg_1, Deg_1 > BottomLeftDiagonal;
+    typedef Matrix<Scalar, Deg_1, 1>               BottomLeftDiagonal;
+
+    typedef Matrix<Scalar, Deg, Deg>               DenseCompanionMatrixType;
+    typedef Matrix< Scalar, _Deg, Deg_1 >          LeftBlock;
+    typedef Matrix< Scalar, Deg_1, Deg_1 >         BottomLeftBlock;
+    typedef Matrix< Scalar, 1, Deg_1 >             LeftBlockFirstRow;
+
+  public:
+    EIGEN_STRONG_INLINE const _Scalar operator()( int row, int col ) const
+    {
+      if( m_bl_diag.rows() > col )
+      {
+        if( 0 < row ){ return m_bl_diag[col]; }
+        else{ return 0; }
+      }
+      else{ return m_monic[row]; }
+    }
+
+  public:
+    template<typename VectorType>
+    void setPolynomial( const VectorType& poly )
+    {
+      const int deg = poly.size()-1;
+      m_monic = -1/poly[deg] * poly.head(deg);
+      //m_bl_diag.setIdentity( deg-1 );
+      m_bl_diag.setOnes(deg-1);
+    }
+
+    template<typename VectorType>
+    ei_companion( const VectorType& poly ){
+      setPolynomial( poly ); }
+
+  public:
+    DenseCompanionMatrixType denseMatrix() const
+    {
+      const int deg   = m_monic.size();
+      const int deg_1 = deg-1;
+      DenseCompanionMatrixType companion(deg,deg);
+      companion <<
+        ( LeftBlock(deg,deg_1)
+          << LeftBlockFirstRow::Zero(1,deg_1),
+          BottomLeftBlock::Identity(deg-1,deg-1)*m_bl_diag.asDiagonal() ).finished()
+        , m_monic;
+      return companion;
+    }
+
+
+
+  protected:
+    /** Helper function for the balancing algorithm.
+     * \returns true if the row and the column, having colNorm and rowNorm
+     * as norms, are balanced, false otherwise.
+     * 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 );
+
+    /** Helper function for the balancing algorithm.
+     * \returns true if the row and the column, having colNorm and rowNorm
+     * as norms, are balanced, false otherwise.
+     * 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 );
+
+  public:
+    /**
+     * Balancing algorithm from B. N. PARLETT and C. REINSCH (1969)
+     * "Balancing a matrix for calculation of eigenvalues and eigenvectors"
+     * adapted to the case of companion matrices.
+     * A matrix with non zero row and non zero column is balanced
+     * for a certain norm if the i-th row and the i-th column
+     * have same norm for all i.
+     */
+    void balance();
+
+  protected:
+      RightColumn                m_monic;
+      BottomLeftDiagonal         m_bl_diag;
+};
+
+
+
+template< typename _Scalar, int _Deg >
+inline
+bool ei_companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm,
+    bool& isBalanced, Scalar& colB, Scalar& rowB )
+{
+  if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; }
+  else
+  {
+    //To find the balancing coefficients, if the radix is 2,
+    //one finds \f$ \sigma \f$ such that
+    // \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 / ei_radix<Scalar>();
+    colB = Scalar(1);
+    const Scalar s = colNorm + rowNorm;
+
+    while (colNorm < rowB)
+    {
+      colB *= ei_radix<Scalar>();
+      colNorm *= ei_radix2<Scalar>();
+    }
+
+    rowB = rowNorm * ei_radix<Scalar>();
+
+    while (colNorm >= rowB)
+    {
+      colB /= ei_radix<Scalar>();
+      colNorm /= ei_radix2<Scalar>();
+    }
+
+    //This line is used to avoid insubstantial balancing
+    if ((rowNorm + colNorm) < Scalar(0.95) * s * colB)
+    {
+      isBalanced = false;
+      rowB = Scalar(1) / colB;
+      return false;
+    }
+    else{
+      return true; }
+  }
+}
+
+template< typename _Scalar, int _Deg >
+inline
+bool ei_companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm,
+    bool& isBalanced, Scalar& colB, Scalar& rowB )
+{
+  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 _Scalar q = colNorm/rowNorm;
+    if( !ei_isApprox( q, _Scalar(1) ) )
+    {
+      rowB = ei_sqrt( colNorm/rowNorm );
+      colB = Scalar(1)/rowB;
+
+      isBalanced = false;
+      return false;
+    }
+    else{
+      return true; }
+  }
+}
+
+
+template< typename _Scalar, int _Deg >
+void ei_companion<_Scalar,_Deg>::balance()
+{
+  EIGEN_STATIC_ASSERT( 1 < Deg, YOU_MADE_A_PROGRAMMING_MISTAKE );
+  const int deg   = m_monic.size();
+  const int deg_1 = deg-1;
+
+  bool hasConverged=false;
+  while( !hasConverged )
+  {
+    hasConverged = true;
+    Scalar colNorm,rowNorm;
+    Scalar colB,rowB;
+
+    //First row, first column excluding the diagonal
+    //==============================================
+    colNorm = ei_abs(m_bl_diag[0]);
+    rowNorm = ei_abs(m_monic[0]);
+
+    //Compute balancing of the row and the column
+    if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
+    {
+      m_bl_diag[0] *= colB;
+      m_monic[0] *= rowB;
+    }
+
+    //Middle rows and columns excluding the diagonal
+    //==============================================
+    for( int i=1; i<deg_1; ++i )
+    {
+      // column norm, excluding the diagonal
+      colNorm = ei_abs(m_bl_diag[i]);
+
+      // row norm, excluding the diagonal
+      rowNorm = ei_abs(m_bl_diag[i-1]) + ei_abs(m_monic[i]);
+
+      //Compute balancing of the row and the column
+      if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
+      {
+        m_bl_diag[i]   *= colB;
+        m_bl_diag[i-1] *= rowB;
+        m_monic[i]     *= rowB;
+      }
+    }
+
+    //Last row, last column excluding the diagonal
+    //============================================
+    const int ebl = m_bl_diag.size()-1;
+    VectorBlock<RightColumn,Deg_1> headMonic( m_monic, 0, deg_1 );
+    colNorm = headMonic.array().abs().sum();
+    rowNorm = ei_abs( m_bl_diag[ebl] );
+
+    //Compute balancing of the row and the column
+    if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
+    {
+      headMonic      *= colB;
+      m_bl_diag[ebl] *= rowB;
+    }
+  }
+}
+
+
+#endif // EIGEN_COMPANION_H
diff --git a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
new file mode 100644
index 000000000..49a5ffffa
--- /dev/null
+++ b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
@@ -0,0 +1,395 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_POLYNOMIAL_SOLVER_H
+#define EIGEN_POLYNOMIAL_SOLVER_H
+
+/** \ingroup Polynomials_Module
+ *  \class PolynomialSolverBase.
+ *
+ * \brief Defined to be inherited by polynomial solvers: it provides
+ * convenient methods such as
+ *  - real roots,
+ *  - greatest, smallest complex roots,
+ *  - real roots with greatest, smallest absolute real value,
+ *  - greatest, smallest real roots.
+ *
+ * It stores the set of roots as a vector of complexes.
+ *
+ */
+template< typename _Scalar, int _Deg >
+class PolynomialSolverBase
+{
+  public:
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
+
+    typedef _Scalar                             Scalar;
+    typedef typename NumTraits<Scalar>::Real    RealScalar;
+    typedef std::complex<RealScalar>            RootType;
+    typedef Matrix<RootType,_Deg,1>             RootsType;
+
+  protected:
+    template< typename OtherPolynomial >
+    inline void setPolynomial( const OtherPolynomial& poly ){
+      m_roots.resize(poly.size()); }
+
+  public:
+    template< typename OtherPolynomial >
+    inline PolynomialSolverBase( const OtherPolynomial& poly ){
+      setPolynomial( poly() ); }
+
+    inline PolynomialSolverBase(){}
+
+  public:
+    /** \returns the complex roots of the polynomial */
+    inline const RootsType& roots() const { return m_roots; }
+
+  public:
+    /** Clear and fills the back insertion sequence with the real roots of the polynomial
+     * i.e. the real part of the complex roots that have an imaginary part which
+     * absolute value is smaller than absImaginaryThreshold.
+     * absImaginaryThreshold takes the dummy_precision associated
+     * with the _Scalar template parameter of the PolynomialSolver class as the default value.
+     *
+     * \param[out] bi_seq : the back insertion sequence (stl concept)
+     * \param[in]  absImaginaryThreshold : the maximum bound of the imaginary part of a complex
+     *  number that is considered as real.
+     * */
+    template<typename Stl_back_insertion_sequence>
+    inline void realRoots( Stl_back_insertion_sequence& bi_seq,
+        const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+    {
+      bi_seq.clear();
+      for( int i=0; i<m_roots.size(); ++i )
+      {
+        if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold ){
+          bi_seq.push_back( m_roots[i].real() ); }
+      }
+    }
+
+  protected:
+    template<typename squaredNormBinaryPredicate>
+    inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
+    {
+      int res=0;
+      RealScalar norm2 = ei_abs2( m_roots[0] );
+      for( int i=1; i<m_roots.size(); ++i )
+      {
+        const RealScalar currNorm2 = ei_abs2( m_roots[i] );
+        if( pred( currNorm2, norm2 ) ){
+          res=i; norm2=currNorm2; }
+      }
+      return m_roots[res];
+    }
+
+  public:
+    /**
+     * \returns the complex root with greatest norm.
+     */
+    inline const RootType& greatestRoot() const
+    {
+      std::greater<Scalar> greater;
+      return selectComplexRoot_withRespectToNorm( greater );
+    }
+
+    /**
+     * \returns the complex root with smallest norm.
+     */
+    inline const RootType& smallestRoot() const
+    {
+      std::less<Scalar> less;
+      return selectComplexRoot_withRespectToNorm( less );
+    }
+
+  protected:
+    template<typename squaredRealPartBinaryPredicate>
+    inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
+        squaredRealPartBinaryPredicate& pred,
+        bool& hasArealRoot,
+        const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+    {
+      hasArealRoot = false;
+      int res=0;
+      RealScalar abs2;
+
+      for( int i=0; i<m_roots.size(); ++i )
+      {
+        if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold )
+        {
+          if( !hasArealRoot )
+          {
+            hasArealRoot = true;
+            res = i;
+            abs2 = m_roots[i].real() * m_roots[i].real();
+          }
+          else
+          {
+            const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
+            if( pred( currAbs2, abs2 ) )
+            {
+              abs2 = currAbs2;
+              res = i;
+            }
+          }
+        }
+        else
+        {
+          if( ei_abs( m_roots[i].imag() ) < ei_abs( m_roots[res].imag() ) ){
+            res = i; }
+        }
+      }
+      return m_roots[res].real();
+    }
+
+
+    template<typename RealPartBinaryPredicate>
+    inline const RealScalar& selectRealRoot_withRespectToRealPart(
+        RealPartBinaryPredicate& pred,
+        bool& hasArealRoot,
+        const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+    {
+      hasArealRoot = false;
+      int res=0;
+      RealScalar val;
+
+      for( int i=0; i<m_roots.size(); ++i )
+      {
+        if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold )
+        {
+          if( !hasArealRoot )
+          {
+            hasArealRoot = true;
+            res = i;
+            val = m_roots[i].real();
+          }
+          else
+          {
+            const RealScalar curr = m_roots[i].real();
+            if( pred( curr, val ) )
+            {
+              val = curr;
+              res = i;
+            }
+          }
+        }
+        else
+        {
+          if( ei_abs( m_roots[i].imag() ) < ei_abs( m_roots[res].imag() ) ){
+            res = i; }
+        }
+      }
+      return m_roots[res].real();
+    }
+
+  public:
+    /**
+     * \returns a real root with greatest absolute magnitude.
+     * A real root is defined as the real part of a complex root with absolute imaginary
+     * part smallest than absImaginaryThreshold.
+     * absImaginaryThreshold takes the dummy_precision associated
+     * with the _Scalar template parameter of the PolynomialSolver class as the default value.
+     * If no real root is found the boolean hasArealRoot is set to false and the real part of
+     * the root with smallest absolute imaginary part is returned instead.
+     *
+     * \param[out] hasArealRoot : boolean true if a real root is found according to the
+     *  absImaginaryThreshold criterion, false otherwise.
+     * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
+     *  whether or not a root is real.
+     */
+    inline const RealScalar& absGreatestRealRoot(
+        bool& hasArealRoot,
+        const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+    {
+      std::greater<Scalar> greater;
+      return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
+    }
+
+
+    /**
+     * \returns a real root with smallest absolute magnitude.
+     * A real root is defined as the real part of a complex root with absolute imaginary
+     * part smallest than absImaginaryThreshold.
+     * absImaginaryThreshold takes the dummy_precision associated
+     * with the _Scalar template parameter of the PolynomialSolver class as the default value.
+     * If no real root is found the boolean hasArealRoot is set to false and the real part of
+     * the root with smallest absolute imaginary part is returned instead.
+     *
+     * \param[out] hasArealRoot : boolean true if a real root is found according to the
+     *  absImaginaryThreshold criterion, false otherwise.
+     * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
+     *  whether or not a root is real.
+     */
+    inline const RealScalar& absSmallestRealRoot(
+        bool& hasArealRoot,
+        const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+    {
+      std::less<Scalar> less;
+      return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
+    }
+
+
+    /**
+     * \returns the real root with greatest value.
+     * A real root is defined as the real part of a complex root with absolute imaginary
+     * part smallest than absImaginaryThreshold.
+     * absImaginaryThreshold takes the dummy_precision associated
+     * with the _Scalar template parameter of the PolynomialSolver class as the default value.
+     * If no real root is found the boolean hasArealRoot is set to false and the real part of
+     * the root with smallest absolute imaginary part is returned instead.
+     *
+     * \param[out] hasArealRoot : boolean true if a real root is found according to the
+     *  absImaginaryThreshold criterion, false otherwise.
+     * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
+     *  whether or not a root is real.
+     */
+    inline const RealScalar& greatestRealRoot(
+        bool& hasArealRoot,
+        const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+    {
+      std::greater<Scalar> greater;
+      return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
+    }
+
+
+    /**
+     * \returns the real root with smallest value.
+     * A real root is defined as the real part of a complex root with absolute imaginary
+     * part smallest than absImaginaryThreshold.
+     * absImaginaryThreshold takes the dummy_precision associated
+     * with the _Scalar template parameter of the PolynomialSolver class as the default value.
+     * If no real root is found the boolean hasArealRoot is set to false and the real part of
+     * the root with smallest absolute imaginary part is returned instead.
+     *
+     * \param[out] hasArealRoot : boolean true if a real root is found according to the
+     *  absImaginaryThreshold criterion, false otherwise.
+     * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
+     *  whether or not a root is real.
+     */
+    inline const RealScalar& smallestRealRoot(
+        bool& hasArealRoot,
+        const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+    {
+      std::less<Scalar> less;
+      return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
+    }
+
+  protected:
+    RootsType               m_roots;
+};
+
+#define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE )  \
+  typedef typename BASE::Scalar                 Scalar;       \
+  typedef typename BASE::RealScalar             RealScalar;   \
+  typedef typename BASE::RootType               RootType;     \
+  typedef typename BASE::RootsType              RootsType;
+
+
+
+/** \ingroup Polynomials_Module
+  *
+  * \class PolynomialSolver
+  *
+  * \brief A polynomial solver
+  *
+  * Computes the complex roots of a real polynomial.
+  *
+  * \param _Scalar the scalar type, i.e., the type of the polynomial coefficients
+  * \param _Deg the degree of the polynomial, can be a compile time value or Dynamic.
+  *             Notice that the number of polynomial coefficients is _Deg+1.
+  *
+  * This class implements a polynomial solver and provides convenient methods such as
+  * - real roots,
+  * - greatest, smallest complex roots,
+  * - real roots with greatest, smallest absolute real value.
+  * - greatest, smallest real roots.
+  *
+  * WARNING: this polynomial solver is experimental, part of the unsuported Eigen modules.
+  *
+  *
+  * Currently a QR algorithm is used to compute the eigenvalues of the companion matrix of
+  * the polynomial to compute its roots.
+  * This supposes that the complex moduli of the roots are all distinct: e.g. there should
+  * be no multiple roots or conjugate roots for instance.
+  * With 32bit (float) floating types this problem shows up frequently.
+  * 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 >
+class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
+{
+  public:
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
+
+    typedef PolynomialSolverBase<_Scalar,_Deg>    PS_Base;
+    EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
+
+    typedef Matrix<Scalar,_Deg,_Deg>                 CompanionMatrixType;
+    typedef EigenSolver<CompanionMatrixType>         EigenSolverType;
+
+  public:
+    /** Computes the complex roots of a new polynomial. */
+    template< typename OtherPolynomial >
+    void compute( const OtherPolynomial& poly )
+    {
+      assert( Scalar(0) != poly[poly.size()-1] );
+      ei_companion<Scalar,_Deg> companion( poly );
+      companion.balance();
+      m_eigenSolver.compute( companion.denseMatrix() );
+      m_roots = m_eigenSolver.eigenvalues();
+    }
+
+  public:
+    template< typename OtherPolynomial >
+    inline PolynomialSolver( const OtherPolynomial& poly ){
+      compute( poly ); }
+
+    inline PolynomialSolver(){}
+
+  protected:
+    using                   PS_Base::m_roots;
+    EigenSolverType         m_eigenSolver;
+};
+
+
+template< typename _Scalar >
+class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
+{
+  public:
+    typedef PolynomialSolverBase<_Scalar,1>    PS_Base;
+    EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
+
+  public:
+    /** Computes the complex roots of a new polynomial. */
+    template< typename OtherPolynomial >
+    void compute( const OtherPolynomial& poly )
+    {
+      assert( Scalar(0) != poly[poly.size()-1] );
+      m_roots[0] = -poly[0]/poly[poly.size()-1];
+    }
+
+  protected:
+    using                   PS_Base::m_roots;
+};
+
+#endif // EIGEN_POLYNOMIAL_SOLVER_H
diff --git a/unsupported/Eigen/src/Polynomials/PolynomialUtils.h b/unsupported/Eigen/src/Polynomials/PolynomialUtils.h
new file mode 100644
index 000000000..d78821f90
--- /dev/null
+++ b/unsupported/Eigen/src/Polynomials/PolynomialUtils.h
@@ -0,0 +1,153 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_POLYNOMIAL_UTILS_H
+#define EIGEN_POLYNOMIAL_UTILS_H
+
+/** \ingroup Polynomials_Module
+ * \returns the evaluation of the polynomial at x using Horner algorithm.
+ *
+ * \param[in] poly : the vector of coefficients of the polynomial ordered
+ *  by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
+ *  e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
+ * \param[in] x : the value to evaluate the polynomial at.
+ *
+ * <i><b>Note for stability:</b></i>
+ *  <dd> \f$ |x| \le 1 \f$ </dd>
+ */
+template <typename Polynomials, typename T>
+inline
+T poly_eval_horner( const Polynomials& poly, const T& x )
+{
+  T val=poly[poly.size()-1];
+  for( int i=poly.size()-2; i>=0; --i ){
+    val = val*x + poly[i]; }
+  return val;
+}
+
+/** \ingroup Polynomials_Module
+ * \returns the evaluation of the polynomial at x using stabilized Horner algorithm.
+ *
+ * \param[in] poly : the vector of coefficients of the polynomial ordered
+ *  by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
+ *  e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
+ * \param[in] x : the value to evaluate the polynomial at.
+ */
+template <typename Polynomials, typename T>
+inline
+T poly_eval( const Polynomials& poly, const T& x )
+{
+  typedef typename NumTraits<T>::Real Real;
+
+  if( ei_abs2( x ) <= Real(1) ){
+    return poly_eval_horner( poly, x ); }
+  else
+  {
+    T val=poly[0];
+    T inv_x = T(1)/x;
+    for( int i=1; i<poly.size(); ++i ){
+      val = val*inv_x + poly[i]; }
+
+    return std::pow(x,(T)(poly.size()-1)) * val;
+  }
+}
+
+/** \ingroup Polynomials_Module
+ * \returns a maximum bound for the absolute value of any root of the polynomial.
+ *
+ * \param[in] poly : the vector of coefficients of the polynomial ordered
+ *  by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
+ *  e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
+ *
+ *  <i><b>Precondition:</b></i>
+ *  <dd> the leading coefficient of the input polynomial poly must be non zero </dd>
+ */
+template <typename Polynomial>
+inline
+typename NumTraits<typename Polynomial::Scalar>::Real cauchy_max_bound( const Polynomial& poly )
+{
+  typedef typename Polynomial::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real Real;
+
+  assert( Scalar(0) != poly[poly.size()-1] );
+  const Scalar inv_leading_coeff = Scalar(1)/poly[poly.size()-1];
+  Real cb(0);
+
+  for( int i=0; i<poly.size()-1; ++i ){
+    cb += ei_abs(poly[i]*inv_leading_coeff); }
+  return cb + Real(1);
+}
+
+/** \ingroup Polynomials_Module
+ * \returns a minimum bound for the absolute value of any non zero root of the polynomial.
+ * \param[in] poly : the vector of coefficients of the polynomial ordered
+ *  by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
+ *  e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
+ */
+template <typename Polynomial>
+inline
+typename NumTraits<typename Polynomial::Scalar>::Real cauchy_min_bound( const Polynomial& poly )
+{
+  typedef typename Polynomial::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real Real;
+
+  int i=0;
+  while( i<poly.size()-1 && Scalar(0) == poly(i) ){ ++i; }
+  if( poly.size()-1 == i ){
+    return Real(1); }
+
+  const Scalar inv_min_coeff = Scalar(1)/poly[i];
+  Real cb(1);
+  for( int j=i+1; j<poly.size(); ++j ){
+    cb += ei_abs(poly[j]*inv_min_coeff); }
+  return Real(1)/cb;
+}
+
+/** \ingroup Polynomials_Module
+ * Given the roots of a polynomial compute the coefficients in the
+ * monomial basis of the monic polynomial with same roots and minimal degree.
+ * If RootVector is a vector of complexes, Polynomial should also be a vector
+ * of complexes.
+ * \param[in] rv : a vector containing the roots of a polynomial.
+ * \param[out] poly : the vector of coefficients of the polynomial ordered
+ *  by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
+ *  e.g. \f$ 3 + x^2 \f$ is stored as a vector \f$ [ 3, 0, 1 ] \f$.
+ */
+template <typename RootVector, typename Polynomial>
+void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
+{
+
+  typedef typename Polynomial::Scalar Scalar;
+
+  poly.setZero( rv.size()+1 );
+  poly[0] = -rv[0]; poly[1] = Scalar(1);
+  for( int i=1; i<(int)rv.size(); ++i )
+  {
+    for( int j=i+1; j>0; --j ){ poly[j] = poly[j-1] - rv[i]*poly[j]; }
+    poly[0] = -rv[i]*poly[0];
+  }
+}
+
+
+#endif // EIGEN_POLYNOMIAL_UTILS_H
diff --git a/unsupported/doc/examples/PolynomialSolver1.cpp b/unsupported/doc/examples/PolynomialSolver1.cpp
new file mode 100644
index 000000000..c875c9361
--- /dev/null
+++ b/unsupported/doc/examples/PolynomialSolver1.cpp
@@ -0,0 +1,53 @@
+#include <unsupported/Eigen/Polynomials>
+#include <vector>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main()
+{
+  typedef Matrix<double,5,1> Vector5d;
+
+  Vector5d roots = Vector5d::Random();
+  cout << "Roots: " << roots.transpose() << endl;
+  Eigen::Matrix<double,6,1> polynomial;
+  roots_to_monicPolynomial( roots, polynomial );
+
+  PolynomialSolver<double,5> psolve( polynomial );
+  cout << "Complex roots: " << psolve.roots().transpose() << endl;
+
+  std::vector<double> realRoots;
+  psolve.realRoots( realRoots );
+  Map<Vector5d> mapRR( &realRoots[0] );
+  cout << "Real roots: " << mapRR.transpose() << endl;
+
+  cout << endl;
+  cout << "Illustration of the convergence problem with the QR algorithm: " << endl;
+  cout << "---------------------------------------------------------------" << endl;
+  Eigen::Matrix<float,7,1> hardCase_polynomial;
+  hardCase_polynomial <<
+  -0.957, 0.9219, 0.3516, 0.9453, -0.4023, -0.5508, -0.03125;
+  cout << "Hard case polynomial defined by floats: " << hardCase_polynomial.transpose() << endl;
+  PolynomialSolver<float,6> psolvef( hardCase_polynomial );
+  cout << "Complex roots: " << psolvef.roots().transpose() << endl;
+  Eigen::Matrix<float,6,1> evals;
+  for( int i=0; i<6; ++i ){ evals[i] = std::abs( poly_eval( hardCase_polynomial, psolvef.roots()[i] ) ); }
+  cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl;
+
+  cout << "Using double's almost always solves the problem for small degrees: " << endl;
+  cout << "-------------------------------------------------------------------" << endl;
+  PolynomialSolver<double,6> psolve6d( hardCase_polynomial.cast<double>() );
+  cout << "Complex roots: " << psolve6d.roots().transpose() << endl;
+  for( int i=0; i<6; ++i )
+  {
+    std::complex<float> castedRoot( psolve6d.roots()[i].real(), psolve6d.roots()[i].imag() );
+    evals[i] = std::abs( poly_eval( hardCase_polynomial, castedRoot ) );
+  }
+  cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl;
+
+  cout.precision(10);
+  cout << "The last root in float then in double: " << psolvef.roots()[5] << "\t" << psolve6d.roots()[5] << endl;
+  std::complex<float> castedRoot( psolve6d.roots()[5].real(), psolve6d.roots()[5].imag() );
+  cout << "Norm of the difference: " << ei_abs( psolvef.roots()[5] - castedRoot ) << endl;
+}
diff --git a/unsupported/doc/examples/PolynomialUtils1.cpp b/unsupported/doc/examples/PolynomialUtils1.cpp
new file mode 100644
index 000000000..dbfe520b5
--- /dev/null
+++ b/unsupported/doc/examples/PolynomialUtils1.cpp
@@ -0,0 +1,20 @@
+#include <unsupported/Eigen/Polynomials>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main()
+{
+  Vector4d roots = Vector4d::Random();
+  cout << "Roots: " << roots.transpose() << endl;
+  Eigen::Matrix<double,5,1> polynomial;
+  roots_to_monicPolynomial( roots, polynomial );
+  cout << "Polynomial: ";
+  for( int i=0; i<4; ++i ){ cout << polynomial[i] << ".x^" << i << "+ "; }
+  cout << polynomial[4] << ".x^4" << endl;
+  Vector4d evaluation;
+  for( int i=0; i<4; ++i ){
+    evaluation[i] = poly_eval( polynomial, roots[i] ); }
+  cout << "Evaluation of the polynomial at the roots: " << evaluation.transpose();
+}
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index 9618a5a56..f01d99c36 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -24,3 +24,19 @@ if(FFTW_FOUND)
   ei_add_test(FFTW  "-DEIGEN_FFTW_DEFAULT " "-lfftw3 -lfftw3f -lfftw3l" )
 endif(FFTW_FOUND)
 
+find_package(GSL)
+if(GSL_FOUND AND GSL_VERSION_MINOR LESS 9)
+  set(GSL_FOUND "")
+endif(GSL_FOUND AND GSL_VERSION_MINOR LESS 9)
+if(GSL_FOUND)
+  add_definitions("-DHAS_GSL" ${GSL_DEFINITIONS})
+  include_directories(${GSL_INCLUDE_DIR})
+  ei_add_property(EIGEN_TESTED_BACKENDS "GSL, ")
+else(GSL_FOUND)
+  ei_add_property(EIGEN_MISSING_BACKENDS "GSL, ")
+  set(GSL_LIBRARIES " ")
+endif(GSL_FOUND)
+
+ei_add_test(polynomialutils)
+ei_add_test(polynomialsolver " " "${GSL_LIBRARIES}" )
+
diff --git a/unsupported/test/polynomialsolver.cpp b/unsupported/test/polynomialsolver.cpp
new file mode 100644
index 000000000..1f2d7e1f3
--- /dev/null
+++ b/unsupported/test/polynomialsolver.cpp
@@ -0,0 +1,263 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+#include <unsupported/Eigen/Polynomials>
+#include <iostream>
+#include <algorithm>
+
+#ifdef HAS_GSL
+#include "gsl_helper.h"
+#endif
+
+using namespace std;
+
+template<int Size>
+struct ei_increment_if_fixed_size
+{
+  enum {
+    ret = (Size == Dynamic) ? Dynamic : Size+1
+  };
+};
+
+
+
+
+template<int Deg, typename POLYNOMIAL, typename SOLVER>
+bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve )
+{
+  typedef typename POLYNOMIAL::Scalar Scalar;
+
+  typedef typename SOLVER::RootsType    RootsType;
+  typedef Matrix<Scalar,Deg,1>          EvalRootsType;
+
+  const int deg = pols.size()-1;
+
+  psolve.compute( pols );
+  const RootsType& roots( psolve.roots() );
+  EvalRootsType evr( deg );
+  for( int i=0; i<roots.size(); ++i ){
+    evr[i] = std::abs( poly_eval( pols, roots[i] ) ); }
+
+  bool evalToZero = evr.isZero( test_precision<Scalar>() );
+  if( !evalToZero )
+  {
+    cerr << "WRONG root: " << endl;
+    cerr << "Polynomial: " << pols.transpose() << endl;
+    cerr << "Roots found: " << roots.transpose() << endl;
+    cerr << "Abs value of the polynomial at the roots: " << evr.transpose() << endl;
+    cerr << endl;
+  }
+
+  #ifdef HAS_GSL
+  if (ei_is_same_type< Scalar, double>::ret)
+  {
+    typedef GslTraits<Scalar> Gsl;
+    RootsType gslRoots(deg);
+    Gsl::eigen_poly_solve( pols, gslRoots );
+    EvalRootsType gslEvr( deg );
+    for( int i=0; i<gslRoots.size(); ++i )
+    {
+      gslEvr[i] = std::abs( poly_eval( pols, gslRoots[i] ) );
+    }
+    bool gslEvalToZero = gslEvr.isZero( test_precision<Scalar>() );
+    if( !evalToZero )
+    {
+      if( !gslEvalToZero ){
+        cerr << "GSL also failed" << endl; }
+      else{
+        cerr << "GSL did NOT failed" << endl; }
+      cerr << "GSL roots found: " << gslRoots.transpose() << endl;
+      cerr << "Abs value of the polynomial at the GSL roots: " << gslEvr.transpose() << endl;
+      cerr << endl;
+    }
+  }
+  #endif //< HAS_GSL
+
+
+  std::vector<Scalar> rootModuli( roots.size() );
+  Map< EvalRootsType > aux( &rootModuli[0], roots.size() );
+  aux = roots.array().abs();
+  std::sort( rootModuli.begin(), rootModuli.end() );
+  bool distinctModuli=true;
+  for( size_t i=1; i<rootModuli.size() && distinctModuli; ++i )
+  {
+    if( ei_isApprox( rootModuli[i], rootModuli[i-1] ) ){
+      distinctModuli = false; }
+  }
+  VERIFY( evalToZero || !distinctModuli );
+
+  return distinctModuli;
+}
+
+
+
+
+
+
+
+template<int Deg, typename POLYNOMIAL>
+void evalSolver( const POLYNOMIAL& pols )
+{
+  typedef typename POLYNOMIAL::Scalar Scalar;
+
+  typedef PolynomialSolver<Scalar, Deg >              PolynomialSolverType;
+
+  PolynomialSolverType psolve;
+  aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve );
+}
+
+
+
+
+template< int Deg, typename POLYNOMIAL, typename ROOTS, typename REAL_ROOTS >
+void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const REAL_ROOTS& real_roots )
+{
+  typedef typename POLYNOMIAL::Scalar Scalar;
+
+  typedef PolynomialSolver<Scalar, Deg >              PolynomialSolverType;
+
+  PolynomialSolverType psolve;
+  if( aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve ) )
+  {
+    //It is supposed that
+    // 1) the roots found are correct
+    // 2) the roots have distinct moduli
+
+    typedef typename POLYNOMIAL::Scalar                 Scalar;
+    typedef typename REAL_ROOTS::Scalar                 Real;
+
+    typedef PolynomialSolver<Scalar, Deg >              PolynomialSolverType;
+    typedef typename PolynomialSolverType::RootsType    RootsType;
+    typedef Matrix<Scalar,Deg,1>                        EvalRootsType;
+
+    //Test realRoots
+    std::vector< Real > calc_realRoots;
+    psolve.realRoots( calc_realRoots );
+    VERIFY( calc_realRoots.size() == (size_t)real_roots.size() );
+
+    const Scalar psPrec = ei_sqrt( test_precision<Scalar>() );
+
+    for( size_t i=0; i<calc_realRoots.size(); ++i )
+    {
+      bool found = false;
+      for( size_t j=0; j<calc_realRoots.size()&& !found; ++j )
+      {
+        if( ei_isApprox( calc_realRoots[i], real_roots[j] ), psPrec ){
+          found = true; }
+      }
+      VERIFY( found );
+    }
+
+    //Test greatestRoot
+    VERIFY( ei_isApprox( roots.array().abs().maxCoeff(),
+          ei_abs( psolve.greatestRoot() ), psPrec ) );
+
+    //Test smallestRoot
+    VERIFY( ei_isApprox( roots.array().abs().minCoeff(),
+          ei_abs( psolve.smallestRoot() ), psPrec ) );
+
+    bool hasRealRoot;
+    //Test absGreatestRealRoot
+    Real r = psolve.absGreatestRealRoot( hasRealRoot );
+    VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
+    if( hasRealRoot ){
+      VERIFY( ei_isApprox( real_roots.array().abs().maxCoeff(), ei_abs(r), psPrec ) );  }
+
+    //Test absSmallestRealRoot
+    r = psolve.absSmallestRealRoot( hasRealRoot );
+    VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
+    if( hasRealRoot ){
+      VERIFY( ei_isApprox( real_roots.array().abs().minCoeff(), ei_abs( r ), psPrec ) ); }
+
+    //Test greatestRealRoot
+    r = psolve.greatestRealRoot( hasRealRoot );
+    VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
+    if( hasRealRoot ){
+      VERIFY( ei_isApprox( real_roots.array().maxCoeff(), r, psPrec ) ); }
+
+    //Test smallestRealRoot
+    r = psolve.smallestRealRoot( hasRealRoot );
+    VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
+    if( hasRealRoot ){
+    VERIFY( ei_isApprox( real_roots.array().minCoeff(), r, psPrec ) ); }
+  }
+}
+
+
+template<typename _Scalar, int _Deg>
+void polynomialsolver(int deg)
+{
+  typedef ei_increment_if_fixed_size<_Deg>            Dim;
+  typedef Matrix<_Scalar,Dim::ret,1>                  PolynomialType;
+  typedef Matrix<_Scalar,_Deg,1>                      EvalRootsType;
+
+  cout << "Standard cases" << endl;
+  PolynomialType pols = PolynomialType::Random(deg+1);
+  evalSolver<_Deg,PolynomialType>( pols );
+
+  cout << "Hard cases" << endl;
+  _Scalar multipleRoot = ei_random<_Scalar>();
+  EvalRootsType allRoots = EvalRootsType::Constant(deg,multipleRoot);
+  roots_to_monicPolynomial( allRoots, pols );
+  evalSolver<_Deg,PolynomialType>( pols );
+
+  cout << "Test sugar" << endl;
+  EvalRootsType realRoots = EvalRootsType::Random(deg);
+  roots_to_monicPolynomial( realRoots, pols );
+  evalSolverSugarFunction<_Deg>(
+      pols,
+      realRoots.template cast <
+                    std::complex<
+                         typename NumTraits<_Scalar>::Real
+                         >
+                    >(),
+      realRoots );
+}
+
+
+template<typename _Scalar> void polynomialsolver_scalar()
+{
+  CALL_SUBTEST_1( (polynomialsolver<_Scalar,1>(1)) );
+  CALL_SUBTEST_2( (polynomialsolver<_Scalar,2>(2)) );
+  CALL_SUBTEST_3( (polynomialsolver<_Scalar,3>(3)) );
+  CALL_SUBTEST_4( (polynomialsolver<_Scalar,4>(4)) );
+  CALL_SUBTEST_5( (polynomialsolver<_Scalar,5>(5)) );
+  CALL_SUBTEST_6( (polynomialsolver<_Scalar,6>(6)) );
+  CALL_SUBTEST_7( (polynomialsolver<_Scalar,7>(7)) );
+  CALL_SUBTEST_8( (polynomialsolver<_Scalar,8>(8)) );
+
+  CALL_SUBTEST_9( (polynomialsolver<_Scalar,Dynamic>(
+          ei_random<int>(9,45)
+          )) );
+}
+
+void test_polynomialsolver()
+{
+  for(int i = 0; i < g_repeat; i++)
+  {
+    polynomialsolver_scalar<double>();
+    polynomialsolver_scalar<float>();
+  }
+}
diff --git a/unsupported/test/polynomialutils.cpp b/unsupported/test/polynomialutils.cpp
new file mode 100644
index 000000000..7f93c2f0d
--- /dev/null
+++ b/unsupported/test/polynomialutils.cpp
@@ -0,0 +1,124 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+#include <unsupported/Eigen/Polynomials>
+#include <iostream>
+
+using namespace std;
+
+template<int Size>
+struct ei_increment_if_fixed_size
+{
+  enum {
+    ret = (Size == Dynamic) ? Dynamic : Size+1
+  };
+};
+
+template<typename _Scalar, int _Deg>
+void realRoots_to_monicPolynomial_test(int deg)
+{
+  typedef ei_increment_if_fixed_size<_Deg>            Dim;
+  typedef Matrix<_Scalar,Dim::ret,1>                  PolynomialType;
+  typedef Matrix<_Scalar,_Deg,1>                      EvalRootsType;
+
+  PolynomialType pols(deg+1);
+  EvalRootsType roots = EvalRootsType::Random(deg);
+  roots_to_monicPolynomial( roots, pols );
+
+  EvalRootsType evr( deg );
+  for( int i=0; i<roots.size(); ++i ){
+    evr[i] = std::abs( poly_eval( pols, roots[i] ) ); }
+
+  bool evalToZero = evr.isZero( test_precision<_Scalar>() );
+  if( !evalToZero ){
+    cerr << evr.transpose() << endl; }
+  VERIFY( evalToZero );
+}
+
+template<typename _Scalar> void realRoots_to_monicPolynomial_scalar()
+{
+  CALL_SUBTEST_2( (realRoots_to_monicPolynomial_test<_Scalar,2>(2)) );
+  CALL_SUBTEST_3( (realRoots_to_monicPolynomial_test<_Scalar,3>(3)) );
+  CALL_SUBTEST_4( (realRoots_to_monicPolynomial_test<_Scalar,4>(4)) );
+  CALL_SUBTEST_5( (realRoots_to_monicPolynomial_test<_Scalar,5>(5)) );
+  CALL_SUBTEST_6( (realRoots_to_monicPolynomial_test<_Scalar,6>(6)) );
+  CALL_SUBTEST_7( (realRoots_to_monicPolynomial_test<_Scalar,7>(7)) );
+  CALL_SUBTEST_8( (realRoots_to_monicPolynomial_test<_Scalar,17>(17)) );
+
+  CALL_SUBTEST_9( (realRoots_to_monicPolynomial_test<_Scalar,Dynamic>(
+          ei_random<int>(18,26) )) );
+}
+
+
+
+
+template<typename _Scalar, int _Deg>
+void CauchyBounds(int deg)
+{
+  typedef ei_increment_if_fixed_size<_Deg>            Dim;
+  typedef Matrix<_Scalar,Dim::ret,1>                  PolynomialType;
+  typedef Matrix<_Scalar,_Deg,1>                      EvalRootsType;
+
+  PolynomialType pols(deg+1);
+  EvalRootsType roots = EvalRootsType::Random(deg);
+  roots_to_monicPolynomial( roots, pols );
+  _Scalar M = cauchy_max_bound( pols );
+  _Scalar m = cauchy_min_bound( pols );
+  _Scalar Max = roots.array().abs().maxCoeff();
+  _Scalar min = roots.array().abs().minCoeff();
+  bool eval = (M >= Max) && (m <= min);
+  if( !eval )
+  {
+    cerr << "Roots: " << roots << endl;
+    cerr << "Bounds: (" << m << ", " << M << ")" << endl;
+    cerr << "Min,Max: (" << min << ", " << Max << ")" << endl;
+  }
+  VERIFY( eval );
+}
+
+template<typename _Scalar> void CauchyBounds_scalar()
+{
+  CALL_SUBTEST_2( (CauchyBounds<_Scalar,2>(2)) );
+  CALL_SUBTEST_3( (CauchyBounds<_Scalar,3>(3)) );
+  CALL_SUBTEST_4( (CauchyBounds<_Scalar,4>(4)) );
+  CALL_SUBTEST_5( (CauchyBounds<_Scalar,5>(5)) );
+  CALL_SUBTEST_6( (CauchyBounds<_Scalar,6>(6)) );
+  CALL_SUBTEST_7( (CauchyBounds<_Scalar,7>(7)) );
+  CALL_SUBTEST_8( (CauchyBounds<_Scalar,17>(17)) );
+
+  CALL_SUBTEST_9( (CauchyBounds<_Scalar,Dynamic>(
+          ei_random<int>(18,26) )) );
+}
+
+void test_polynomialutils()
+{
+  for(int i = 0; i < g_repeat; i++)
+  {
+    realRoots_to_monicPolynomial_scalar<double>();
+    realRoots_to_monicPolynomial_scalar<float>();
+    CauchyBounds_scalar<double>();
+    CauchyBounds_scalar<float>();
+  }
+}