Change return type of matrixH() method to HouseholderSequence.

This method is a member of Tridiagonalization and HessenbergDecomposition.

5 files changed, 92 insertions, 81 deletions

diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h
index 52d0fc661..c69e3eafd 100644
--- a/Eigen/src/Eigenvalues/ComplexSchur.h
+++ b/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -3,6 +3,7 @@
 //
 // Copyright (C) 2009 Claire Maurice
 // 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
@@ -26,6 +27,8 @@
 #ifndef EIGEN_COMPLEX_SCHUR_H
 #define EIGEN_COMPLEX_SCHUR_H
 
+template<typename MatrixType, bool IsComplex> struct ei_complex_schur_reduce_to_hessenberg;
+
 /** \eigenvalues_module \ingroup Eigenvalues_Module
   * \nonstableyet
   *
@@ -50,6 +53,8 @@
   * decomposition is computed, you can use the matrixU() and matrixT()
   * functions to retrieve the matrices U and V in the decomposition.
   *
+  * \note This code is inspired from Jampack
+  *
   * \sa class RealSchur, class EigenSolver, class ComplexEigenSolver
   */
 template<typename _MatrixType> class ComplexSchur
@@ -194,6 +199,8 @@ template<typename _MatrixType> class ComplexSchur
   private:  
     bool subdiagonalEntryIsNeglegible(int i);
     ComplexScalar computeShift(int iu, int iter);
+    void reduceToTriangularForm(bool skipU);
+    friend struct ei_complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
 };
 
 /** Computes the principal value of the square root of the complex \a z. */
@@ -290,12 +297,10 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
 template<typename MatrixType>
 void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
 {
-  // this code is inspired from Jampack
   m_matUisUptodate = false;
   ei_assert(matrix.cols() == matrix.rows());
-  int n = matrix.cols();
 
-  if(n==1)
+  if(matrix.cols() == 1)
   {
     m_matU = ComplexMatrixType::Identity(1,1);
     if(!skipU) m_matT = matrix.template cast<ComplexScalar>();
@@ -304,15 +309,49 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
     return;
   }
 
-  // Reduce to Hessenberg form
-  // TODO skip Q if skipU = true
-  m_hess.compute(matrix);
+  ei_complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, skipU);
+  reduceToTriangularForm(skipU);
+}
 
-  m_matT = m_hess.matrixH().template cast<ComplexScalar>();
-  if(!skipU)  m_matU = m_hess.matrixQ().template cast<ComplexScalar>();
+/* Reduce given matrix to Hessenberg form */
+template<typename MatrixType, bool IsComplex>
+struct ei_complex_schur_reduce_to_hessenberg 
+{
+  // this is the implementation for the case IsComplex = true
+  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool skipU)
+  {
+    // TODO skip Q if skipU = true
+    _this.m_hess.compute(matrix);
+    _this.m_matT = _this.m_hess.matrixH();
+    if(!skipU)  _this.m_matU = _this.m_hess.matrixQ();
+  }
+};
 
-  // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
+template<typename MatrixType>
+struct ei_complex_schur_reduce_to_hessenberg<MatrixType, false>
+{
+  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool skipU)
+  {
+    typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
+    typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType;
+
+    // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
+    // TODO skip Q if skipU = true
+    _this.m_hess.compute(matrix);
+    _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
+    if(!skipU)  
+    {
+      // This may cause an allocation which seems to be avoidable
+      MatrixType Q = _this.m_hess.matrixQ(); 
+      _this.m_matU = Q.template cast<ComplexScalar>();
+    }
+  }
+};
 
+// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
+template<typename MatrixType>
+void ComplexSchur<MatrixType>::reduceToTriangularForm(bool skipU)
+{  
   // 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 submatrix).
@@ -352,7 +391,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
     ComplexScalar shift = computeShift(iu, iter);
     PlanarRotation<ComplexScalar> rot;
     rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
-    m_matT.rightCols(n-il).applyOnTheLeft(il, il+1, rot.adjoint());
+    m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
     m_matT.topRows(std::min(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
     if(!skipU) m_matU.applyOnTheRight(il, il+1, rot);
 
@@ -360,7 +399,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
     {
       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.rightCols(n-i).applyOnTheLeft(i, i+1, rot.adjoint());
+      m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
       m_matT.topRows(std::min(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
       if(!skipU) m_matU.applyOnTheRight(i, i+1, rot);
     }
diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
index a1ac31981..30d657dfd 100644
--- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h
+++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
@@ -2,6 +2,7 @@
 // for linear algebra.
 //
 // 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
@@ -59,7 +60,9 @@ template<typename _MatrixType> class HessenbergDecomposition
 {
   public:
 
+    /** \brief Synonym for the template parameter \p _MatrixType. */
     typedef _MatrixType MatrixType;
+
     enum {
       Size = MatrixType::RowsAtCompileTime,
       SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
@@ -68,17 +71,20 @@ template<typename _MatrixType> class HessenbergDecomposition
       MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
     };
 
-    /** \brief Scalar type for matrices of type \p _MatrixType. */
+    /** \brief Scalar type for matrices of type #MatrixType. */
     typedef typename MatrixType::Scalar Scalar;
 
     /** \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.
+      * vector is one less than the size of #MatrixType, if it is a fixed-side
+      * type.
       */
     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
 
+    /** \brief Return type of matrixQ() */
+    typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
+
     /** \brief Default constructor; the decomposition will be computed later.
       *
       * \param [in] size  The size of the matrix whose Hessenberg decomposition will be computed.
@@ -190,19 +196,22 @@ template<typename _MatrixType> class HessenbergDecomposition
 
     /** \brief Reconstructs the orthogonal matrix Q in the decomposition 
       *
-      * \returns the matrix Q
+      * \returns object representing 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.
+      * This function returns a light-weight object of template class
+      * HouseholderSequence. You can either apply it directly to a matrix or
+      * you can convert it to a matrix of type #MatrixType.
       *
-      * \sa matrixH() for an example
+      * \sa matrixH() for an example, class HouseholderSequence
       */
-    MatrixType matrixQ() const;
+    HouseholderSequenceType matrixQ() const
+    {
+      return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1);
+    }
 
     /** \brief Constructs the Hessenberg matrix H in the decomposition
       *
@@ -281,21 +290,6 @@ void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVector
   }
 }
 
-template<typename MatrixType>
-typename HessenbergDecomposition<MatrixType>::MatrixType
-HessenbergDecomposition<MatrixType>::matrixQ() const
-{
-  int n = m_matrix.rows();
-  MatrixType matQ = MatrixType::Identity(n,n);
-  VectorType temp(n);
-  for (int i = n-2; i>=0; i--)
-  {
-    matQ.bottomRightCorner(n-i-1,n-i-1)
-        .applyHouseholderOnTheLeft(m_matrix.col(i).tail(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &temp.coeffRef(0,0));
-  }
-  return matQ;
-}
-
 #endif // EIGEN_HIDE_HEAVY_CODE
 
 template<typename MatrixType>
diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h
index 024590f3c..43509863a 100644
--- a/Eigen/src/Eigenvalues/Tridiagonalization.h
+++ b/Eigen/src/Eigenvalues/Tridiagonalization.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
@@ -61,7 +62,9 @@ template<typename _MatrixType> class Tridiagonalization
 {
   public:
 
+    /** \brief Synonym for the template parameter \p _MatrixType. */
     typedef _MatrixType MatrixType;
+
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<Scalar>::Real RealScalar;
 
@@ -89,6 +92,9 @@ template<typename _MatrixType> class Tridiagonalization
                 Block<MatrixType,SizeMinusOne,SizeMinusOne>,0 >
             >::ret SubDiagonalReturnType;
 
+    /** \brief Return type of matrixQ() */
+    typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
+
     /** \brief Default constructor.
       *
       * \param [in]  size  Positive integer, size of the matrix whose tridiagonal
@@ -195,29 +201,25 @@ template<typename _MatrixType> class Tridiagonalization
       */
     inline const MatrixType& packedMatrix() const { return m_matrix; }
 
-    /** \brief Reconstructs the unitary matrix Q in the decomposition 
+    /** \brief Returns the unitary matrix Q in the decomposition 
       *
-      * \returns the matrix Q
+      * \returns object representing the matrix Q
       *
       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
       * the member function compute(const MatrixType&) has been called before
       * to compute the tridiagonal 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.
+      * This function returns a light-weight object of template class
+      * HouseholderSequence. You can either apply it directly to a matrix or
+      * you can convert it to a matrix of type #MatrixType.
       *
       * \sa Tridiagonalization(const MatrixType&) for an example,
-      * matrixT(), matrixQInPlace()
-      */
-    MatrixType matrixQ() const;
-
-    /** \brief Reconstructs the unitary matrix Q in the decomposition 
-      *
-      * This is an in-place variant of matrixQ() which avoids the copy. 
-      * This function will probably be deleted soon.
+      *     matrixT(), class HouseholderSequence
       */
-    template<typename QDerived> void matrixQInPlace(MatrixBase<QDerived>* q) const;
+    HouseholderSequenceType matrixQ() const
+    {
+      return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1);
+    }
 
     /** \brief Constructs the tridiagonal matrix T in the decomposition
       *
@@ -387,31 +389,6 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
 }
 
 template<typename MatrixType>
-typename Tridiagonalization<MatrixType>::MatrixType
-Tridiagonalization<MatrixType>::matrixQ() const
-{
-  MatrixType matQ;
-  matrixQInPlace(&matQ);
-  return matQ;
-}
-
-template<typename MatrixType>
-template<typename QDerived>
-void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) const
-{
-  QDerived& matQ = q->derived();
-  int n = m_matrix.rows();
-  matQ = MatrixType::Identity(n,n);
-  typedef typename ei_plain_row_type<MatrixType>::type RowVectorType;
-  RowVectorType aux(n);
-  for (int i = n-2; i>=0; i--)
-  {
-    matQ.bottomRightCorner(n-i-1,n-i-1)
-        .applyHouseholderOnTheLeft(m_matrix.col(i).tail(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &aux.coeffRef(0,0));
-  }
-}
-
-template<typename MatrixType>
 void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
 {
   int n = mat.rows();
@@ -426,7 +403,7 @@ void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalT
     diag = tridiag.diagonal();
     subdiag = tridiag.subDiagonal();
     if (extractQ)
-      tridiag.matrixQInPlace(&mat);
+      mat = tridiag.matrixQ();
   }
 }
 
diff --git a/doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp b/doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp
index 109650041..a26012433 100644
--- a/doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp
+++ b/doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp
@@ -1,11 +1,9 @@
 MatrixXd X = MatrixXd::Random(5,5);
 MatrixXd A = X + X.transpose();
 cout << "Here is a random symmetric 5x5 matrix:" << endl << A << endl << endl;
-
 Tridiagonalization<MatrixXd> triOfA(A);
-cout << "The orthogonal matrix Q is:" << endl << triOfA.matrixQ() << endl;
-cout << "The tridiagonal matrix T is:" << endl << triOfA.matrixT() << endl << endl;
-
 MatrixXd Q = triOfA.matrixQ();
+cout << "The orthogonal matrix Q is:" << endl << Q << endl;
 MatrixXd T = triOfA.matrixT();
+cout << "The tridiagonal matrix T is:" << endl << T << endl << endl;
 cout << "Q * T * Q^T = " << endl << Q * T * Q.transpose() << endl;
diff --git a/test/hessenberg.cpp b/test/hessenberg.cpp
index ec1148bfc..61ff98150 100644
--- a/test/hessenberg.cpp
+++ b/test/hessenberg.cpp
@@ -34,8 +34,9 @@ template<typename Scalar,int Size> void hessenberg(int size = Size)
   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 Q = hess.matrixQ();
     MatrixType H = hess.matrixH();
+    VERIFY_IS_APPROX(m, Q * H * Q.adjoint());
     for(int row = 2; row < size; ++row) {
       for(int col = 0; col < row-1; ++col) {
 	VERIFY(H(row,col) == (typename MatrixType::Scalar)0);
@@ -48,8 +49,10 @@ template<typename Scalar,int Size> void hessenberg(int 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());
+  MatrixType cs1Q = cs1.matrixQ();
+  MatrixType cs2Q = cs2.matrixQ();  
+  VERIFY_IS_EQUAL(cs1Q, cs2Q);
 
   // TODO: Add tests for packedMatrix() and householderCoefficients()
 }