From 0d6321225786bc3d95f1dbe5236f07c5e5e96179 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 26 Nov 2010 15:31:47 +0100 Subject: [PATCH] add a TridiagonalizationMatrixTReturnType class to make Tridiagonalization::matrixT() more efficient and future proof. --- Eigen/src/Eigenvalues/Tridiagonalization.h | 83 +++++++++++++++------- test/eigensolver_selfadjoint.cpp | 5 ++ 2 files changed, 62 insertions(+), 26 deletions(-) diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h index d70c4d433..b97710d61 100644 --- a/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -26,6 +26,16 @@ #ifndef EIGEN_TRIDIAGONALIZATION_H #define EIGEN_TRIDIAGONALIZATION_H +template struct TridiagonalizationMatrixTReturnType; + +namespace internal { +template +struct traits > +{ + typedef typename MatrixType::PlainObject ReturnType; +}; +} + namespace internal { template void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs); @@ -85,6 +95,7 @@ template class Tridiagonalization typedef Matrix CoeffVectorType; typedef typename internal::plain_col_type::type DiagonalType; typedef Matrix SubDiagonalType; + typedef typename internal::remove_all::type MatrixTypeRealView; typedef typename internal::conditional::IsComplex, typename Diagonal::RealReturnType, @@ -244,24 +255,28 @@ template class Tridiagonalization return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1); } - /** \brief Constructs the tridiagonal matrix T in the decomposition + /** \brief Returns an expression of the tridiagonal matrix T in the decomposition * - * \returns the matrix T + * \returns expression object representing the matrix T * * \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 copies the matrix T from internal data. The diagonal and - * subdiagonal of the packed matrix as returned by packedMatrix() - * represents the matrix T. It may sometimes be sufficient to directly use - * the packed matrix or the vector expressions returned by diagonal() - * and subDiagonal() instead of creating a new matrix with this function. + * Currently, this function can be used to extract the matrix T from internal + * data and copy it to a dense matrix object. In most cases, it may be + * sufficient to directly use the packed matrix or the vector expressions + * returned by diagonal() and subDiagonal() instead of creating a new + * dense copy matrix with this function. * * \sa Tridiagonalization(const MatrixType&) for an example, * matrixQ(), packedMatrix(), diagonal(), subDiagonal() */ - MatrixType matrixT() const; + TridiagonalizationMatrixTReturnType matrixT() const + { + eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); + return TridiagonalizationMatrixTReturnType(m_matrix.real()); + } /** \brief Returns the diagonal of the tridiagonal matrix T in the decomposition. * @@ -314,24 +329,6 @@ Tridiagonalization::subDiagonal() const return Block(m_matrix, 1, 0, n-1,n-1).diagonal(); } -template -typename Tridiagonalization::MatrixType -Tridiagonalization::matrixT() const -{ - // FIXME should this function (and other similar ones) rather take a matrix as argument - // and fill it ? (to avoid temporaries) - eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); - Index n = m_matrix.rows(); - MatrixType matT = m_matrix; - matT.topRightCorner(n-1, n-1).diagonal() = subDiagonal().template cast().conjugate(); - if (n>2) - { - matT.topRightCorner(n-2, n-2).template triangularView().setZero(); - matT.bottomLeftCorner(n-2, n-2).template triangularView().setZero(); - } - return matT; -} - namespace internal { /** \internal @@ -530,4 +527,38 @@ struct tridiagonalization_inplace_selector } // end namespace internal +/** \eigenvalues_module \ingroup Eigenvalues_Module + * + * + * \brief Expression type for return value of Tridiagonalization::matrixT() + * + * \tparam MatrixType type of underlying dense matrix + */ +template struct TridiagonalizationMatrixTReturnType +: public ReturnByValue > +{ + typedef typename MatrixType::Index Index; + public: + /** \brief Constructor. + * + * \param[in] mat The underlying dense matrix + */ + TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { } + + template + inline void evalTo(ResultType& result) const + { + result.setZero(); + result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate(); + result.template diagonal() = m_matrix.template diagonal(); + result.template diagonal<-1>() = m_matrix.template diagonal<-1>(); + } + + Index rows() const { return m_matrix.rows(); } + Index cols() const { return m_matrix.cols(); } + + protected: + const typename MatrixType::Nested m_matrix; +}; + #endif // EIGEN_TRIDIAGONALIZATION_H diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 806fe5831..b85bcc289 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -155,6 +155,11 @@ template void selfadjointeigensolver(const MatrixType& m) VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt()); VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt()); + // test Tridiagonalization's methods + Tridiagonalization tridiag(symmA); + // FIXME tridiag.matrixQ().adjoint() does not work + VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint()); + if (rows > 1) { // Test matrix with NaN