diff options
Diffstat (limited to 'eigen/Eigen/src/Eigenvalues/HessenbergDecomposition.h')
| -rw-r--r-- | eigen/Eigen/src/Eigenvalues/HessenbergDecomposition.h | 373 |
1 files changed, 373 insertions, 0 deletions
diff --git a/eigen/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/eigen/Eigen/src/Eigenvalues/HessenbergDecomposition.h new file mode 100644 index 0000000..3db0c01 --- /dev/null +++ b/eigen/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -0,0 +1,373 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_HESSENBERGDECOMPOSITION_H +#define EIGEN_HESSENBERGDECOMPOSITION_H + +namespace Eigen { + +namespace internal { + +template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType; +template<typename MatrixType> +struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> > +{ + typedef MatrixType ReturnType; +}; + +} + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * + * + * \class HessenbergDecomposition + * + * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation + * + * \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$. 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$). + * + * 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 +{ + public: + + /** \brief Synonym for the template parameter \p _MatrixType. */ + typedef _MatrixType MatrixType; + + enum { + Size = MatrixType::RowsAtCompileTime, + SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1, + Options = MatrixType::Options, + MaxSize = MatrixType::MaxRowsAtCompileTime, + MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1 + }; + + /** \brief Scalar type for matrices of type #MatrixType. */ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::Index Index; + + /** \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 #MatrixType, if it is a fixed-side + * type. + */ + typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType; + + /** \brief Return type of matrixQ() */ + typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType; + + typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType; + + /** \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(Index size = Size==Dynamic ? 2 : Size) + : m_matrix(size,size), + m_temp(size), + m_isInitialized(false) + { + if(size>1) + 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), + m_temp(matrix.rows()), + m_isInitialized(false) + { + if(matrix.rows()<2) + { + m_isInitialized = true; + return; + } + m_hCoeffs.resize(matrix.rows()-1,1); + _compute(m_matrix, m_hCoeffs, m_temp); + m_isInitialized = true; + } + + /** \brief Computes Hessenberg decomposition of given matrix. + * + * \param[in] matrix Square matrix whose Hessenberg decomposition is to be computed. + * \returns Reference to \c *this + * + * 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 reuses of the allocated data in the HessenbergDecomposition + * object. + * + * Example: \include HessenbergDecomposition_compute.cpp + * Output: \verbinclude HessenbergDecomposition_compute.out + */ + HessenbergDecomposition& compute(const MatrixType& matrix) + { + m_matrix = matrix; + if(matrix.rows()<2) + { + m_isInitialized = true; + return *this; + } + m_hCoeffs.resize(matrix.rows()-1,1); + _compute(m_matrix, m_hCoeffs, m_temp); + m_isInitialized = true; + return *this; + } + + /** \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(), \ref Householder_Module "Householder module" + */ + const CoeffVectorType& householderCoefficients() const + { + eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); + return m_hCoeffs; + } + + /** \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 + * \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() const + { + eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); + return m_matrix; + } + + /** \brief Reconstructs the orthogonal matrix Q in the decomposition + * + * \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 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, class HouseholderSequence + */ + HouseholderSequenceType matrixQ() const + { + eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); + return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate()) + .setLength(m_matrix.rows() - 1) + .setShift(1); + } + + /** \brief Constructs the Hessenberg matrix H in the decomposition + * + * \returns expression object representing 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. + * + * The object returned by this function constructs the Hessenberg matrix H + * when it is assigned to a matrix or otherwise evaluated. The matrix H is + * constructed from the packed matrix as returned by packedMatrix(): The + * upper part (including the subdiagonal) of the packed matrix contains + * the matrix H. It may sometimes be better to directly use the packed + * matrix instead of constructing the matrix H. + * + * Example: \include HessenbergDecomposition_matrixH.cpp + * Output: \verbinclude HessenbergDecomposition_matrixH.out + * + * \sa matrixQ(), packedMatrix() + */ + MatrixHReturnType matrixH() const + { + eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); + return MatrixHReturnType(*this); + } + + private: + + typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType; + typedef typename NumTraits<Scalar>::Real RealScalar; + static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp); + + protected: + MatrixType m_matrix; + CoeffVectorType m_hCoeffs; + VectorType m_temp; + bool m_isInitialized; +}; + +/** \internal + * Performs a tridiagonal decomposition of \a matA in place. + * + * \param matA the input selfadjoint matrix + * \param hCoeffs returned Householder coefficients + * + * The result is written in the lower triangular part of \a matA. + * + * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1. + * + * \sa packedMatrix() + */ +template<typename MatrixType> +void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp) +{ + eigen_assert(matA.rows()==matA.cols()); + Index n = matA.rows(); + temp.resize(n); + for (Index i = 0; i<n-1; ++i) + { + // let's consider the vector v = i-th column starting at position i+1 + Index remainingSize = n-i-1; + RealScalar beta; + Scalar h; + matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta); + matA.col(i).coeffRef(i+1) = beta; + hCoeffs.coeffRef(i) = h; + + // Apply similarity transformation to remaining columns, + // i.e., compute A = H A H' + + // A = H A + matA.bottomRightCorner(remainingSize, remainingSize) + .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0)); + + // A = A H' + matA.rightCols(remainingSize) + .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0)); + } +} + +namespace internal { + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * + * + * \brief Expression type for return value of HessenbergDecomposition::matrixH() + * + * \tparam MatrixType type of matrix in the Hessenberg decomposition + * + * Objects of this type represent the Hessenberg matrix in the Hessenberg + * decomposition of some matrix. The object holds a reference to the + * HessenbergDecomposition class until the it is assigned or evaluated for + * some other reason (the reference should remain valid during the life time + * of this object). This class is the return type of + * HessenbergDecomposition::matrixH(); there is probably no other use for this + * class. + */ +template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType +: public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> > +{ + typedef typename MatrixType::Index Index; + public: + /** \brief Constructor. + * + * \param[in] hess Hessenberg decomposition + */ + HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { } + + /** \brief Hessenberg matrix in decomposition. + * + * \param[out] result Hessenberg matrix in decomposition \p hess which + * was passed to the constructor + */ + template <typename ResultType> + inline void evalTo(ResultType& result) const + { + result = m_hess.packedMatrix(); + Index n = result.rows(); + if (n>2) + result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero(); + } + + Index rows() const { return m_hess.packedMatrix().rows(); } + Index cols() const { return m_hess.packedMatrix().cols(); } + + protected: + const HessenbergDecomposition<MatrixType>& m_hess; +}; + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_HESSENBERGDECOMPOSITION_H |