Allow user to compute only the eigenvalues and not the eigenvectors.

11 files changed, 235 insertions, 155 deletions

diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index f56815c15..a3a4a4eba 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -56,7 +56,10 @@
 template<typename _MatrixType> class ComplexEigenSolver
 {
   public:
+
+    /** \brief Synonym for the template parameter \p _MatrixType. */
     typedef _MatrixType MatrixType;
+
     enum {
       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
@@ -65,12 +68,12 @@ template<typename _MatrixType> class ComplexEigenSolver
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };
 
-    /** \brief Scalar type for matrices of type \p _MatrixType. */
+    /** \brief Scalar type for matrices of type #MatrixType. */
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<Scalar>::Real RealScalar;
     typedef typename MatrixType::Index Index;
 
-    /** \brief Complex scalar type for \p _MatrixType. 
+    /** \brief Complex scalar type for #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
@@ -81,14 +84,14 @@ template<typename _MatrixType> class ComplexEigenSolver
     /** \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.
+      * The length of the vector is the size of #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.
+      * The size is the same as the size of #MatrixType.
       */
     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, ColsAtCompileTime> EigenvectorType;
 
@@ -102,6 +105,7 @@ template<typename _MatrixType> class ComplexEigenSolver
               m_eivalues(),
               m_schur(),
               m_isInitialized(false),
+              m_eigenvectorsOk(false),
               m_matX()
     {}
     
@@ -116,40 +120,46 @@ template<typename _MatrixType> class ComplexEigenSolver
               m_eivalues(size),
               m_schur(size),
               m_isInitialized(false),
+              m_eigenvectorsOk(false),
               m_matX(size, size)
     {}
 
     /** \brief Constructor; computes eigendecomposition of given matrix. 
       * 
       * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
+      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
+      *    eigenvalues are computed; if false, only the eigenvalues are
+      *    computed. 
       *
       * This constructor calls compute() to compute the eigendecomposition.
       */
-    ComplexEigenSolver(const MatrixType& matrix)
+      ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
             : m_eivec(matrix.rows(),matrix.cols()),
               m_eivalues(matrix.cols()),
               m_schur(matrix.rows()),
               m_isInitialized(false),
+              m_eigenvectorsOk(false),
               m_matX(matrix.rows(),matrix.cols())
     {
-      compute(matrix);
+      compute(matrix, computeEigenvectors);
     }
 
     /** \brief Returns the eigenvectors of given matrix. 
       *
       * \returns  A const reference to the matrix whose columns are the eigenvectors.
       *
-      * 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.
+      * \pre Either the constructor
+      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
+      * function compute(const MatrixType& matrix, bool) has been called before
+      * to compute the eigendecomposition of a matrix, and 
+      * \p computeEigenvectors was set to true (the default).
+      *
+      * 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
@@ -157,6 +167,7 @@ template<typename _MatrixType> class ComplexEigenSolver
     const EigenvectorType& eigenvectors() const
     {
       ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
+      ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
       return m_eivec;
     }
 
@@ -164,11 +175,12 @@ template<typename _MatrixType> class ComplexEigenSolver
       *
       * \returns A const reference to the column vector containing the eigenvalues.
       *
-      * 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
+      * \pre Either the constructor
+      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
+      * function compute(const MatrixType& matrix, bool) 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.
@@ -185,10 +197,14 @@ template<typename _MatrixType> class ComplexEigenSolver
     /** \brief Computes eigendecomposition of given matrix. 
       * 
       * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
+      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
+      *    eigenvalues are computed; if false, only the eigenvalues are
+      *    computed. 
       *
-      * This function computes the eigenvalues and eigenvectors of \p
-      * matrix.  The eigenvalues() and eigenvectors() functions can be
-      * used to retrieve the computed eigendecomposition.
+      * This function computes the eigenvalues of the complex matrix \p matrix.
+      * The eigenvalues() function can be used to retrieve them.  If 
+      * \p computeEigenvectors is true, then the eigenvectors are also computed
+      * and can be retrieved by calling eigenvectors().
       *
       * The matrix is first reduced to Schur form using the
       * ComplexSchur class. The Schur decomposition is then used to
@@ -201,19 +217,20 @@ template<typename _MatrixType> class ComplexEigenSolver
       * Example: \include ComplexEigenSolver_compute.cpp
       * Output: \verbinclude ComplexEigenSolver_compute.out
       */
-    void compute(const MatrixType& matrix);
+    void compute(const MatrixType& matrix, bool computeEigenvectors = true);
 
   protected:
     EigenvectorType m_eivec;
     EigenvalueType m_eivalues;
     ComplexSchur<MatrixType> m_schur;
     bool m_isInitialized;
+    bool m_eigenvectorsOk;
     EigenvectorType m_matX;
 };
 
 
 template<typename MatrixType>
-void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
+void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
 {
   // this code is inspired from Jampack
   assert(matrix.cols() == matrix.rows());
@@ -222,40 +239,45 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
 
   // Step 1: Do a complex Schur decomposition, A = U T U^*
   // The eigenvalues are on the diagonal of T.
-  m_schur.compute(matrix);
+  m_schur.compute(matrix, computeEigenvectors);
   m_eivalues = m_schur.matrixT().diagonal();
 
-  // Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T.
-  // The matrix X is unit triangular.
-  m_matX = EigenvectorType::Zero(n, n);
-  for(Index k=n-1 ; k>=0 ; k--)
+  if(computeEigenvectors)
   {
-    m_matX.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(Index i=k-1 ; i>=0 ; i--)
+    // Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T.
+    // The matrix X is unit triangular.
+    m_matX = EigenvectorType::Zero(n, n);
+    for(Index k=n-1 ; k>=0 ; k--)
     {
-      m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
-      if(k-i-1>0)
-        m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
-      ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
-      if(z==ComplexScalar(0))
+      m_matX.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(Index i=k-1 ; i>=0 ; i--)
       {
-	// 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;
+        m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
+        if(k-i-1>0)
+          m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
+        ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_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;
+        }
+        m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
       }
-      m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
+    }
+  
+    // 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.noalias() = m_schur.matrixU() * m_matX;
+    // .. and normalize the eigenvectors
+    for(Index k=0 ; k<n ; k++)
+    {
+      m_eivec.col(k).normalize();
     }
   }
 
-  // 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.noalias() = m_schur.matrixU() * m_matX;
-  // .. and normalize the eigenvectors
-  for(Index k=0 ; k<n ; k++)
-  {
-    m_eivec.col(k).normalize();
-  }
   m_isInitialized = true;
+  m_eigenvectorsOk = computeEigenvectors;
 
   // Step 4: Sort the eigenvalues
   for (Index i=0; i<n; i++)
@@ -266,7 +288,8 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
     {
       k += i;
       std::swap(m_eivalues[k],m_eivalues[i]);
-      m_eivec.col(i).swap(m_eivec.col(k));
+      if(computeEigenvectors)
+	m_eivec.col(i).swap(m_eivec.col(k));
     }
   }
 }
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index 5400fdaf2..f745413a8 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -57,16 +57,16 @@
   * this variant of the eigendecomposition the pseudo-eigendecomposition.
   *
   * 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
+  * a given matrix. Alternatively, you can use the 
+  * EigenSolver(const MatrixType&, bool) 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.
   *
-  * The documentation for EigenSolver(const MatrixType&) contains an example of
-  * the typical use of this class.
+  * The documentation for EigenSolver(const MatrixType&, bool) 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).
@@ -78,7 +78,9 @@ template<typename _MatrixType> class EigenSolver
 {
   public:
 
+    /** \brief Synonym for the template parameter \p _MatrixType. */
     typedef _MatrixType MatrixType;
+
     enum {
       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
@@ -87,12 +89,12 @@ template<typename _MatrixType> class EigenSolver
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };
 
-    /** \brief Scalar type for matrices of type \p _MatrixType. */
+    /** \brief Scalar type for matrices of type #MatrixType. */
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<Scalar>::Real RealScalar;
     typedef typename MatrixType::Index Index;
 
-    /** \brief Complex scalar type for \p _MatrixType. 
+    /** \brief Complex scalar type for #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
@@ -103,27 +105,27 @@ template<typename _MatrixType> class EigenSolver
     /** \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.
+      * The length of the vector is the size of #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.
+      * The size is the same as the size of #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&).
+      * perform decompositions via EigenSolver::compute(const MatrixType&, bool).
       *
       * \sa compute() for an example.
       */
  EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
 
-    /** \brief Default Constructor with memory preallocation
+    /** \brief Default constructor with memory preallocation
       *
       * Like the default constructor but with preallocation of the internal data
       * according to the specified problem \a size.
@@ -133,6 +135,7 @@ template<typename _MatrixType> class EigenSolver
       : m_eivec(size, size),
         m_eivalues(size),
         m_isInitialized(false),
+        m_eigenvectorsOk(false),
         m_realSchur(size),
         m_matT(size, size), 
         m_tmp(size)
@@ -141,6 +144,9 @@ template<typename _MatrixType> class EigenSolver
     /** \brief Constructor; computes eigendecomposition of given matrix. 
       * 
       * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
+      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
+      *    eigenvalues are computed; if false, only the eigenvalues are
+      *    computed. 
       *
       * This constructor calls compute() to compute the eigenvalues
       * and eigenvectors.
@@ -150,23 +156,26 @@ template<typename _MatrixType> class EigenSolver
       *
       * \sa compute()
       */
-    EigenSolver(const MatrixType& matrix)
+    EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
       : m_eivec(matrix.rows(), matrix.cols()),
         m_eivalues(matrix.cols()),
         m_isInitialized(false),
+        m_eigenvectorsOk(false),
         m_realSchur(matrix.cols()),
         m_matT(matrix.rows(), matrix.cols()), 
         m_tmp(matrix.cols())
     {
-      compute(matrix);
+      compute(matrix, computeEigenvectors);
     }
 
     /** \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.
+      * \pre Either the constructor 
+      * EigenSolver(const MatrixType&,bool) or the member function
+      * compute(const MatrixType&, bool) has been called before, and
+      * \p computeEigenvectors was set to true (the default).
       *
       * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
       * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
@@ -185,9 +194,10 @@ template<typename _MatrixType> class EigenSolver
       *
       * \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.
+      * \pre Either the constructor 
+      * EigenSolver(const MatrixType&,bool) or the member function
+      * compute(const MatrixType&, bool) has been called before, and
+      * \p computeEigenvectors was set to true (the default).
       *
       * The real matrix \f$ V \f$ returned by this function and the
       * block-diagonal matrix \f$ D \f$ returned by pseudoEigenvalueMatrix()
@@ -201,6 +211,7 @@ template<typename _MatrixType> class EigenSolver
     const MatrixType& pseudoEigenvectors() const
     {
       ei_assert(m_isInitialized && "EigenSolver is not initialized.");
+      ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
       return m_eivec;
     }
 
@@ -208,8 +219,9 @@ template<typename _MatrixType> class EigenSolver
       *
       * \returns  A block-diagonal matrix.
       *
-      * \pre Either the constructor EigenSolver(const MatrixType&) or the
-      * member function compute(const MatrixType&) has been called before.
+      * \pre Either the constructor 
+      * EigenSolver(const MatrixType&,bool) or the member function
+      * compute(const MatrixType&, bool) 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
@@ -226,8 +238,9 @@ template<typename _MatrixType> class EigenSolver
       *
       * \returns A const reference to the column vector containing the eigenvalues.
       *
-      * \pre Either the constructor EigenSolver(const MatrixType&) or the
-      * member function compute(const MatrixType&) has been called before.
+      * \pre Either the constructor 
+      * EigenSolver(const MatrixType&,bool) or the member function
+      * compute(const MatrixType&, bool) has been called before.
       *
       * The eigenvalues are repeated according to their algebraic multiplicity,
       * so there are as many eigenvalues as rows in the matrix.
@@ -247,34 +260,40 @@ template<typename _MatrixType> class EigenSolver
     /** \brief Computes eigendecomposition of given matrix. 
       * 
       * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
+      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
+      *    eigenvalues are computed; if false, only the eigenvalues are
+      *    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.
+      * This function computes the eigenvalues of the real matrix \p matrix.
+      * The eigenvalues() function can be used to retrieve them.  If 
+      * \p computeEigenvectors is true, then the eigenvectors are also computed
+      * and can be retrieved by calling eigenvectors().
       *
       * 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.
+      * 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) if \p computeEigenvectors 
+      * is true, and \f$ 10n^3 \f$ if \p computeEigenvectors is false.
       *
       * 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);
+    EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
 
   private:
-    void computeEigenvectors();
+    void doComputeEigenvectors();
 
   protected:
     MatrixType m_eivec;
     EigenvalueType m_eivalues;
     bool m_isInitialized;
+    bool m_eigenvectorsOk;
     RealSchur<MatrixType> m_realSchur;
     MatrixType m_matT;
 
@@ -286,7 +305,7 @@ template<typename MatrixType>
 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
 {
   ei_assert(m_isInitialized && "EigenSolver is not initialized.");
-  Index n = m_eivec.cols();
+  Index n = m_eivalues.rows();
   MatrixType matD = MatrixType::Zero(n,n);
   for (Index i=0; i<n; ++i)
   {
@@ -306,6 +325,7 @@ template<typename MatrixType>
 typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
 {
   ei_assert(m_isInitialized && "EigenSolver is not initialized.");
+  ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
   Index n = m_eivec.cols();
   EigenvectorsType matV(n,n);
   for (Index j=0; j<n; ++j)
@@ -332,14 +352,15 @@ typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige
 }
 
 template<typename MatrixType>
-EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix)
+EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
 {
   assert(matrix.cols() == matrix.rows());
 
   // Reduce to real Schur form.
-  m_realSchur.compute(matrix);
+  m_realSchur.compute(matrix, computeEigenvectors);
   m_matT = m_realSchur.matrixT();
-  m_eivec = m_realSchur.matrixU();
+  if (computeEigenvectors)
+    m_eivec = m_realSchur.matrixU();
 
   // Compute eigenvalues from matT
   m_eivalues.resize(matrix.cols());
@@ -362,9 +383,12 @@ EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matr
   }
   
   // Compute eigenvectors.
-  computeEigenvectors();
+  if (computeEigenvectors)
+    doComputeEigenvectors();
 
   m_isInitialized = true;
+  m_eigenvectorsOk = computeEigenvectors;
+
   return *this;
 }
 
@@ -389,7 +413,7 @@ std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
 
 
 template<typename MatrixType>
-void EigenSolver<MatrixType>::computeEigenvectors()
+void EigenSolver<MatrixType>::doComputeEigenvectors()
 {
   const Index size = m_eivec.cols();
   const Scalar eps = NumTraits<Scalar>::epsilon();
@@ -404,7 +428,7 @@ void EigenSolver<MatrixType>::computeEigenvectors()
   // Backsubstitute to find vectors of upper triangular form
   if (norm == 0.0)
   {
-      return;
+    return;
   }
 
   for (Index n = size-1; n >= 0; n--)
diff --git a/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h b/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h
index 7b04e6ba7..f27481fe1 100644
--- a/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h
+++ b/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h
@@ -37,7 +37,7 @@ struct ei_eigenvalues_selector
   {
     typedef typename Derived::PlainObject PlainObject;
     PlainObject m_eval(m);
-    return ComplexEigenSolver<PlainObject>(m_eval).eigenvalues();
+    return ComplexEigenSolver<PlainObject>(m_eval, false).eigenvalues();
   }
 };
 
@@ -49,7 +49,7 @@ struct ei_eigenvalues_selector<Derived, false>
   {
     typedef typename Derived::PlainObject PlainObject;
     PlainObject m_eval(m);
-    return EigenSolver<PlainObject>(m_eval).eigenvalues();
+    return EigenSolver<PlainObject>(m_eval, false).eigenvalues();
   }
 };
 
@@ -101,7 +101,7 @@ SelfAdjointView<MatrixType, UpLo>::eigenvalues() const
 {
   typedef typename SelfAdjointView<MatrixType, UpLo>::PlainObject PlainObject;
   PlainObject thisAsMatrix(*this);
-  return SelfAdjointEigenSolver<PlainObject>(thisAsMatrix).eigenvalues();
+  return SelfAdjointEigenSolver<PlainObject>(thisAsMatrix, false).eigenvalues();
 }
 
 
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index 92ff448ed..c92b72a94 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -50,13 +50,13 @@
   * 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&)
+  * given matrix. Alternatively, you can use the RealSchur(const MatrixType&, bool)
   * 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.
+  * The documentation of RealSchur(const MatrixType&, bool) 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).
@@ -98,41 +98,46 @@ template<typename _MatrixType> class RealSchur
               m_matU(size, size),
               m_workspaceVector(size),
               m_hess(size),
-              m_isInitialized(false)
+              m_isInitialized(false),
+              m_matUisUptodate(false)
     { }
 
     /** \brief Constructor; computes real Schur decomposition of given matrix. 
       * 
-      * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
+      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
+      * \param[in]  computeU  If true, both T and U are computed; if false, only T is 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)
+    RealSchur(const MatrixType& matrix, bool computeU = true)
             : m_matT(matrix.rows(),matrix.cols()),
               m_matU(matrix.rows(),matrix.cols()),
               m_workspaceVector(matrix.rows()),
               m_hess(matrix.rows()),
-              m_isInitialized(false)
+              m_isInitialized(false),
+              m_matUisUptodate(false)
     {
-      compute(matrix);
+      compute(matrix, computeU);
     }
 
     /** \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.
+      * \pre Either the constructor RealSchur(const MatrixType&, bool) or the
+      * member function compute(const MatrixType&, bool) has been called before
+      * to compute the Schur decomposition of a matrix, and \p computeU was set
+      * to true (the default value).
       *
-      * \sa RealSchur(const MatrixType&) for an example
+      * \sa RealSchur(const MatrixType&, bool) for an example
       */
     const MatrixType& matrixU() const
     {
       ei_assert(m_isInitialized && "RealSchur is not initialized.");
+      ei_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition.");
       return m_matU;
     }
 
@@ -140,11 +145,11 @@ template<typename _MatrixType> class RealSchur
       *
       * \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.
+      * \pre Either the constructor RealSchur(const MatrixType&, bool) or the
+      * member function compute(const MatrixType&, bool) has been called before
+      * to compute the Schur decomposition of a matrix.
       *
-      * \sa RealSchur(const MatrixType&) for an example
+      * \sa RealSchur(const MatrixType&, bool) for an example
       */
     const MatrixType& matrixT() const
     {
@@ -154,19 +159,21 @@ template<typename _MatrixType> class RealSchur
   
     /** \brief Computes Schur decomposition of given matrix. 
       * 
-      * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
+      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
+      * \param[in]  computeU  If true, both T and U are computed; if false, only T is 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.
+      * may be taken to be \f$25n^3\f$ flops if \a computeU is true and
+      * \f$10n^3\f$ flops if \a computeU is false.
       *
       * Example: \include RealSchur_compute.cpp
       * Output: \verbinclude RealSchur_compute.out
       */
-    void compute(const MatrixType& matrix);
+    void compute(const MatrixType& matrix, bool computeU = true);
 
   private:
     
@@ -175,38 +182,39 @@ template<typename _MatrixType> class RealSchur
     ColumnVectorType m_workspaceVector;
     HessenbergDecomposition<MatrixType> m_hess;
     bool m_isInitialized;
+    bool m_matUisUptodate;
 
     typedef Matrix<Scalar,3,1> Vector3s;
 
     Scalar computeNormOfT();
     Index findSmallSubdiagEntry(Index iu, Scalar norm);
-    void splitOffTwoRows(Index iu, Scalar exshift);
+    void splitOffTwoRows(Index iu, bool computeU, Scalar exshift);
     void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
     void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
-    void performFrancisQRStep(Index il, Index im, Index iu, const Vector3s& firstHouseholderVector, Scalar* workspace);
+    void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace);
 };
 
 
 template<typename MatrixType>
-void RealSchur<MatrixType>::compute(const MatrixType& matrix)
+void RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
 {
   assert(matrix.cols() == matrix.rows());
 
   // Step 1. Reduce to Hessenberg form
-  // TODO skip Q if skipU = true
   m_hess.compute(matrix);
   m_matT = m_hess.matrixH();
-  m_matU = m_hess.matrixQ();
+  if (computeU)
+    m_matU = m_hess.matrixQ();
 
   // Step 2. Reduce to real Schur form  
-  m_workspaceVector.resize(m_matU.cols());
+  m_workspaceVector.resize(m_matT.cols());
   Scalar* workspace = &m_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.
-  Index iu = m_matU.cols() - 1;
+  Index iu = m_matT.cols() - 1;
   Index iter = 0; // iteration count
   Scalar exshift = 0.0; // sum of exceptional shifts
   Scalar norm = computeNormOfT();
@@ -226,7 +234,7 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
     }
     else if (il == iu-1) // Two roots found
     {
-      splitOffTwoRows(iu, exshift);
+      splitOffTwoRows(iu, computeU, exshift);
       iu -= 2;
       iter = 0;
     }
@@ -237,18 +245,19 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
       iter = iter + 1;   // (Could check iteration count here.)
       Index im;
       initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
-      performFrancisQRStep(il, im, iu, firstHouseholderVector, workspace);
+      performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
     }
   } 
 
   m_isInitialized = true;
+  m_matUisUptodate = computeU;
 }
 
 /** \internal Computes and returns vector L1 norm of T */
 template<typename MatrixType>
 inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
 {
-  const Index size = m_matU.cols();
+  const Index size = m_matT.cols();
   // FIXME to be efficient the following would requires a triangular reduxion code
   // Scalar norm = m_matT.upper().cwiseAbs().sum() 
   //               + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
@@ -277,9 +286,9 @@ inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(I
 
 /** \internal Update T given that rows iu-1 and iu decouple from the rest. */
 template<typename MatrixType>
-inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, Scalar exshift)
+inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scalar exshift)
 {
-  const Index size = m_matU.cols();
+  const Index size = m_matT.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
@@ -300,7 +309,8 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, Scalar exshift)
     m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
     m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
     m_matT.coeffRef(iu, iu-1) = Scalar(0); 
-    m_matU.applyOnTheRight(iu-1, iu, rot);
+    if (computeU)
+      m_matU.applyOnTheRight(iu-1, iu, rot);
   }
 
   if (iu > 1) 
@@ -375,12 +385,12 @@ inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const V
 
 /** \internal Perform a Francis QR step involving rows il:iu and columns im:iu. */
 template<typename MatrixType>
-inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, const Vector3s& firstHouseholderVector, Scalar* workspace)
+inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
 {
   assert(im >= il);
   assert(im <= iu-2);
 
-  const Index size = m_matU.cols();
+  const Index size = m_matT.cols();
 
   for (Index k = im; k <= iu-2; ++k)
   {
@@ -406,7 +416,8 @@ inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Inde
       // 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);
+      if (computeU)
+	m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
     }
   }
 
@@ -420,7 +431,8 @@ inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Inde
     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);
+    if (computeU)
+      m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
   }
 
   // clean up pollution due to round-off errors
diff --git a/doc/snippets/ComplexEigenSolver_eigenvalues.cpp b/doc/snippets/ComplexEigenSolver_eigenvalues.cpp
index 1afa8b086..5509bd897 100644
--- a/doc/snippets/ComplexEigenSolver_eigenvalues.cpp
+++ b/doc/snippets/ComplexEigenSolver_eigenvalues.cpp
@@ -1,4 +1,4 @@
 MatrixXcf ones = MatrixXcf::Ones(3,3);
-ComplexEigenSolver<MatrixXcf> ces(ones);
+ComplexEigenSolver<MatrixXcf> ces(ones, /* computeEigenvectors = */ false);
 cout << "The eigenvalues of the 3x3 matrix of ones are:" 
      << endl << ces.eigenvalues() << endl;
diff --git a/doc/snippets/EigenSolver_compute.cpp b/doc/snippets/EigenSolver_compute.cpp
index 06138f608..a5c96e9b4 100644
--- a/doc/snippets/EigenSolver_compute.cpp
+++ b/doc/snippets/EigenSolver_compute.cpp
@@ -1,6 +1,6 @@
 EigenSolver<MatrixXf> es;
 MatrixXf A = MatrixXf::Random(4,4);
-es.compute(A);
+es.compute(A, /* computeEigenvectors = */ false);
 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
+es.compute(A + MatrixXf::Identity(4,4), false); // 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
index 8d83ea982..ed28869a0 100644
--- a/doc/snippets/EigenSolver_eigenvalues.cpp
+++ b/doc/snippets/EigenSolver_eigenvalues.cpp
@@ -1,4 +1,4 @@
 MatrixXd ones = MatrixXd::Ones(3,3);
-EigenSolver<MatrixXd> es(ones);
+EigenSolver<MatrixXd> es(ones, false);
 cout << "The eigenvalues of the 3x3 matrix of ones are:" 
      << endl << es.eigenvalues() << endl;
diff --git a/doc/snippets/RealSchur_compute.cpp b/doc/snippets/RealSchur_compute.cpp
index 4dcfaf0f2..20c2611b8 100644
--- a/doc/snippets/RealSchur_compute.cpp
+++ b/doc/snippets/RealSchur_compute.cpp
@@ -1,6 +1,6 @@
 MatrixXf A = MatrixXf::Random(4,4);
 RealSchur<MatrixXf> schur(4);
-schur.compute(A);
+schur.compute(A, /* computeU = */ false);
 cout << "The matrix T in the decomposition of A is:" << endl << schur.matrixT() << endl;
-schur.compute(A.inverse());
+schur.compute(A.inverse(), /* computeU = */ false);
 cout << "The matrix T in the decomposition of A^(-1) is:" << endl << schur.matrixT() << endl;
diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp
index 1440cd700..3285d26c2 100644
--- a/test/eigensolver_complex.cpp
+++ b/test/eigensolver_complex.cpp
@@ -2,6 +2,7 @@
 // for linear algebra. Eigen itself is part of the KDE project.
 //
 // Copyright (C) 2008-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
@@ -66,7 +67,10 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
   // Note: If MatrixType is real then a.eigenvalues() uses EigenSolver and thus
   // another algorithm so results may differ slightly
   verify_is_approx_upto_permutation(a.eigenvalues(), ei1.eigenvalues());
-  
+
+  ComplexEigenSolver<MatrixType> eiNoEivecs(a, false);
+  VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
+
   // Regression test for issue #66
   MatrixType z = MatrixType::Zero(rows,cols);
   ComplexEigenSolver<MatrixType> eiz(z);
@@ -76,11 +80,15 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
   VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
 }
 
-template<typename MatrixType> void eigensolver_verify_assert()
+template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
 {
   ComplexEigenSolver<MatrixType> eig;
-  VERIFY_RAISES_ASSERT(eig.eigenvectors())
-  VERIFY_RAISES_ASSERT(eig.eigenvalues())
+  VERIFY_RAISES_ASSERT(eig.eigenvectors());
+  VERIFY_RAISES_ASSERT(eig.eigenvalues());
+
+  MatrixType a = MatrixType::Random(m.rows(),m.cols());
+  eig.compute(a, false);
+  VERIFY_RAISES_ASSERT(eig.eigenvectors());
 }
 
 void test_eigensolver_complex()
@@ -92,10 +100,10 @@ void test_eigensolver_complex()
     CALL_SUBTEST_4( eigensolver(Matrix3f()) );
   }
 
-  CALL_SUBTEST_1( eigensolver_verify_assert<Matrix4cf>() );
-  CALL_SUBTEST_2( eigensolver_verify_assert<MatrixXcd>() );
-  CALL_SUBTEST_3(( eigensolver_verify_assert<Matrix<std::complex<float>, 1, 1> >() ));
-  CALL_SUBTEST_4( eigensolver_verify_assert<Matrix3f>() );
+  CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4cf()) );
+  CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXcd(14,14)) );
+  CALL_SUBTEST_3( eigensolver_verify_assert(Matrix<std::complex<float>, 1, 1>()) );
+  CALL_SUBTEST_4( eigensolver_verify_assert(Matrix3f()) );
 
   // Test problem size constructors
   CALL_SUBTEST_5(ComplexEigenSolver<MatrixXf>(10));
diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp
index d70f37ea4..79c08ec31 100644
--- a/test/eigensolver_generic.cpp
+++ b/test/eigensolver_generic.cpp
@@ -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
@@ -60,19 +61,26 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
                    ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
   VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues());
 
+  EigenSolver<MatrixType> eiNoEivecs(a, false);
+  VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
+  VERIFY_IS_APPROX(ei1.pseudoEigenvalueMatrix(), eiNoEivecs.pseudoEigenvalueMatrix());
+
   MatrixType id = MatrixType::Identity(rows, cols);
   VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
 }
 
-template<typename MatrixType> void eigensolver_verify_assert()
+template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
 {
-  MatrixType tmp;
-
   EigenSolver<MatrixType> eig;
-  VERIFY_RAISES_ASSERT(eig.eigenvectors())
-  VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors())
-  VERIFY_RAISES_ASSERT(eig.pseudoEigenvalueMatrix())
-  VERIFY_RAISES_ASSERT(eig.eigenvalues())
+  VERIFY_RAISES_ASSERT(eig.eigenvectors());
+  VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
+  VERIFY_RAISES_ASSERT(eig.pseudoEigenvalueMatrix());
+  VERIFY_RAISES_ASSERT(eig.eigenvalues());
+
+  MatrixType a = MatrixType::Random(m.rows(),m.cols());
+  eig.compute(a, false);
+  VERIFY_RAISES_ASSERT(eig.eigenvectors());
+  VERIFY_RAISES_ASSERT(eig.pseudoEigenvectors());
 }
 
 void test_eigensolver_generic()
@@ -88,11 +96,11 @@ void test_eigensolver_generic()
     CALL_SUBTEST_4( eigensolver(Matrix2d()) );
   }
 
-  CALL_SUBTEST_1( eigensolver_verify_assert<Matrix4f>() );
-  CALL_SUBTEST_2( eigensolver_verify_assert<MatrixXd>() );
-  CALL_SUBTEST_4( eigensolver_verify_assert<Matrix2d>() );
-  CALL_SUBTEST_5( eigensolver_verify_assert<MatrixXf>() );
+  CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4f()) );
+  CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXd(17,17)) );
+  CALL_SUBTEST_3( eigensolver_verify_assert(Matrix<double,1,1>()) );
+  CALL_SUBTEST_4( eigensolver_verify_assert(Matrix2d()) );
 
   // Test problem size constructors
-  CALL_SUBTEST_6(EigenSolver<MatrixXf>(10));
+  CALL_SUBTEST_5(EigenSolver<MatrixXf>(10));
 }
diff --git a/test/schur_real.cpp b/test/schur_real.cpp
index bcb19c936..116c8dbce 100644
--- a/test/schur_real.cpp
+++ b/test/schur_real.cpp
@@ -73,6 +73,11 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
   RealSchur<MatrixType> rs2(A);
   VERIFY_IS_EQUAL(rs1.matrixT(), rs2.matrixT());
   VERIFY_IS_EQUAL(rs1.matrixU(), rs2.matrixU());
+
+  // Test computation of only T, not U
+  RealSchur<MatrixType> rsOnlyT(A, false);
+  VERIFY_IS_EQUAL(rs1.matrixT(), rsOnlyT.matrixT());
+  VERIFY_RAISES_ASSERT(rsOnlyT.matrixU());
 }
 
 void test_schur_real()