Change return type of matrixH() method to HouseholderSequence. This method is a member of Tridiagonalization and HessenbergDecomposition.
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h index 52d0fc6..c69e3ea 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 @@ 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 @@ 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 @@ 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 @@ 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 @@ { 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 a1ac319..30d657d 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 @@ { 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 @@ 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 @@ /** \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 @@ } } -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 024590f..4350986 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 @@ { 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 @@ 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 @@ */ 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() + * matrixT(), class HouseholderSequence */ - 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. - */ - 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 @@ } 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 @@ 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 1096500..a260124 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 ec1148b..61ff981 100644 --- a/test/hessenberg.cpp +++ b/test/hessenberg.cpp
@@ -34,8 +34,9 @@ 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 @@ 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() }