1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11#ifndef EIGEN_HESSENBERGDECOMPOSITION_H
12#define EIGEN_HESSENBERGDECOMPOSITION_H
13
14namespace Eigen {
15
16namespace internal {
17
18template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
19template<typename MatrixType>
20struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
21{
22  typedef MatrixType ReturnType;
23};
24
25}
26
27/** \eigenvalues_module \ingroup Eigenvalues_Module
28  *
29  *
30  * \class HessenbergDecomposition
31  *
32  * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation
33  *
34  * \tparam _MatrixType the type of the matrix of which we are computing the Hessenberg decomposition
35  *
36  * This class performs an Hessenberg decomposition of a matrix \f$ A \f$. In
37  * the real case, the Hessenberg decomposition consists of an orthogonal
38  * matrix \f$ Q \f$ and a Hessenberg matrix \f$ H \f$ such that \f$ A = Q H
39  * Q^T \f$. An orthogonal matrix is a matrix whose inverse equals its
40  * transpose (\f$ Q^{-1} = Q^T \f$). A Hessenberg matrix has zeros below the
41  * subdiagonal, so it is almost upper triangular. The Hessenberg decomposition
42  * of a complex matrix is \f$ A = Q H Q^* \f$ with \f$ Q \f$ unitary (that is,
43  * \f$ Q^{-1} = Q^* \f$).
44  *
45  * Call the function compute() to compute the Hessenberg decomposition of a
46  * given matrix. Alternatively, you can use the
47  * HessenbergDecomposition(const MatrixType&) constructor which computes the
48  * Hessenberg decomposition at construction time. Once the decomposition is
49  * computed, you can use the matrixH() and matrixQ() functions to construct
50  * the matrices H and Q in the decomposition.
51  *
52  * The documentation for matrixH() contains an example of the typical use of
53  * this class.
54  *
55  * \sa class ComplexSchur, class Tridiagonalization, \ref QR_Module "QR Module"
56  */
57template<typename _MatrixType> class HessenbergDecomposition
58{
59  public:
60
61    /** \brief Synonym for the template parameter \p _MatrixType. */
62    typedef _MatrixType MatrixType;
63
64    enum {
65      Size = MatrixType::RowsAtCompileTime,
66      SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
67      Options = MatrixType::Options,
68      MaxSize = MatrixType::MaxRowsAtCompileTime,
69      MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
70    };
71
72    /** \brief Scalar type for matrices of type #MatrixType. */
73    typedef typename MatrixType::Scalar Scalar;
74    typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
75
76    /** \brief Type for vector of Householder coefficients.
77      *
78      * This is column vector with entries of type #Scalar. The length of the
79      * vector is one less than the size of #MatrixType, if it is a fixed-side
80      * type.
81      */
82    typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
83
84    /** \brief Return type of matrixQ() */
85    typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType;
86
87    typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType;
88
89    /** \brief Default constructor; the decomposition will be computed later.
90      *
91      * \param [in] size  The size of the matrix whose Hessenberg decomposition will be computed.
92      *
93      * The default constructor is useful in cases in which the user intends to
94      * perform decompositions via compute().  The \p size parameter is only
95      * used as a hint. It is not an error to give a wrong \p size, but it may
96      * impair performance.
97      *
98      * \sa compute() for an example.
99      */
100    explicit HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
101      : m_matrix(size,size),
102        m_temp(size),
103        m_isInitialized(false)
104    {
105      if(size>1)
106        m_hCoeffs.resize(size-1);
107    }
108
109    /** \brief Constructor; computes Hessenberg decomposition of given matrix.
110      *
111      * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
112      *
113      * This constructor calls compute() to compute the Hessenberg
114      * decomposition.
115      *
116      * \sa matrixH() for an example.
117      */
118    template<typename InputType>
119    explicit HessenbergDecomposition(const EigenBase<InputType>& matrix)
120      : m_matrix(matrix.derived()),
121        m_temp(matrix.rows()),
122        m_isInitialized(false)
123    {
124      if(matrix.rows()<2)
125      {
126        m_isInitialized = true;
127        return;
128      }
129      m_hCoeffs.resize(matrix.rows()-1,1);
130      _compute(m_matrix, m_hCoeffs, m_temp);
131      m_isInitialized = true;
132    }
133
134    /** \brief Computes Hessenberg decomposition of given matrix.
135      *
136      * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
137      * \returns    Reference to \c *this
138      *
139      * The Hessenberg decomposition is computed by bringing the columns of the
140      * matrix successively in the required form using Householder reflections
141      * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix
142      * Computations</i>). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$
143      * denotes the size of the given matrix.
144      *
145      * This method reuses of the allocated data in the HessenbergDecomposition
146      * object.
147      *
148      * Example: \include HessenbergDecomposition_compute.cpp
149      * Output: \verbinclude HessenbergDecomposition_compute.out
150      */
151    template<typename InputType>
152    HessenbergDecomposition& compute(const EigenBase<InputType>& matrix)
153    {
154      m_matrix = matrix.derived();
155      if(matrix.rows()<2)
156      {
157        m_isInitialized = true;
158        return *this;
159      }
160      m_hCoeffs.resize(matrix.rows()-1,1);
161      _compute(m_matrix, m_hCoeffs, m_temp);
162      m_isInitialized = true;
163      return *this;
164    }
165
166    /** \brief Returns the Householder coefficients.
167      *
168      * \returns a const reference to the vector of Householder coefficients
169      *
170      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
171      * or the member function compute(const MatrixType&) has been called
172      * before to compute the Hessenberg decomposition of a matrix.
173      *
174      * The Householder coefficients allow the reconstruction of the matrix
175      * \f$ Q \f$ in the Hessenberg decomposition from the packed data.
176      *
177      * \sa packedMatrix(), \ref Householder_Module "Householder module"
178      */
179    const CoeffVectorType& householderCoefficients() const
180    {
181      eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
182      return m_hCoeffs;
183    }
184
185    /** \brief Returns the internal representation of the decomposition
186      *
187      *	\returns a const reference to a matrix with the internal representation
188      *	         of the decomposition.
189      *
190      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
191      * or the member function compute(const MatrixType&) has been called
192      * before to compute the Hessenberg decomposition of a matrix.
193      *
194      * The returned matrix contains the following information:
195      *  - the upper part and lower sub-diagonal represent the Hessenberg matrix H
196      *  - the rest of the lower part contains the Householder vectors that, combined with
197      *    Householder coefficients returned by householderCoefficients(),
198      *    allows to reconstruct the matrix Q as
199      *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
200      *    Here, the matrices \f$ H_i \f$ are the Householder transformations
201      *       \f$ H_i = (I - h_i v_i v_i^T) \f$
202      *    where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
203      *    \f$ v_i \f$ is the Householder vector defined by
204      *       \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
205      *    with M the matrix returned by this function.
206      *
207      * See LAPACK for further details on this packed storage.
208      *
209      * Example: \include HessenbergDecomposition_packedMatrix.cpp
210      * Output: \verbinclude HessenbergDecomposition_packedMatrix.out
211      *
212      * \sa householderCoefficients()
213      */
214    const MatrixType& packedMatrix() const
215    {
216      eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
217      return m_matrix;
218    }
219
220    /** \brief Reconstructs the orthogonal matrix Q in the decomposition
221      *
222      * \returns object representing the matrix Q
223      *
224      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
225      * or the member function compute(const MatrixType&) has been called
226      * before to compute the Hessenberg decomposition of a matrix.
227      *
228      * This function returns a light-weight object of template class
229      * HouseholderSequence. You can either apply it directly to a matrix or
230      * you can convert it to a matrix of type #MatrixType.
231      *
232      * \sa matrixH() for an example, class HouseholderSequence
233      */
234    HouseholderSequenceType matrixQ() const
235    {
236      eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
237      return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
238             .setLength(m_matrix.rows() - 1)
239             .setShift(1);
240    }
241
242    /** \brief Constructs the Hessenberg matrix H in the decomposition
243      *
244      * \returns expression object representing the matrix H
245      *
246      * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
247      * or the member function compute(const MatrixType&) has been called
248      * before to compute the Hessenberg decomposition of a matrix.
249      *
250      * The object returned by this function constructs the Hessenberg matrix H
251      * when it is assigned to a matrix or otherwise evaluated. The matrix H is
252      * constructed from the packed matrix as returned by packedMatrix(): The
253      * upper part (including the subdiagonal) of the packed matrix contains
254      * the matrix H. It may sometimes be better to directly use the packed
255      * matrix instead of constructing the matrix H.
256      *
257      * Example: \include HessenbergDecomposition_matrixH.cpp
258      * Output: \verbinclude HessenbergDecomposition_matrixH.out
259      *
260      * \sa matrixQ(), packedMatrix()
261      */
262    MatrixHReturnType matrixH() const
263    {
264      eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
265      return MatrixHReturnType(*this);
266    }
267
268  private:
269
270    typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
271    typedef typename NumTraits<Scalar>::Real RealScalar;
272    static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
273
274  protected:
275    MatrixType m_matrix;
276    CoeffVectorType m_hCoeffs;
277    VectorType m_temp;
278    bool m_isInitialized;
279};
280
281/** \internal
282  * Performs a tridiagonal decomposition of \a matA in place.
283  *
284  * \param matA the input selfadjoint matrix
285  * \param hCoeffs returned Householder coefficients
286  *
287  * The result is written in the lower triangular part of \a matA.
288  *
289  * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1.
290  *
291  * \sa packedMatrix()
292  */
293template<typename MatrixType>
294void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
295{
296  eigen_assert(matA.rows()==matA.cols());
297  Index n = matA.rows();
298  temp.resize(n);
299  for (Index i = 0; i<n-1; ++i)
300  {
301    // let's consider the vector v = i-th column starting at position i+1
302    Index remainingSize = n-i-1;
303    RealScalar beta;
304    Scalar h;
305    matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
306    matA.col(i).coeffRef(i+1) = beta;
307    hCoeffs.coeffRef(i) = h;
308
309    // Apply similarity transformation to remaining columns,
310    // i.e., compute A = H A H'
311
312    // A = H A
313    matA.bottomRightCorner(remainingSize, remainingSize)
314        .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
315
316    // A = A H'
317    matA.rightCols(remainingSize)
318        .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0));
319  }
320}
321
322namespace internal {
323
324/** \eigenvalues_module \ingroup Eigenvalues_Module
325  *
326  *
327  * \brief Expression type for return value of HessenbergDecomposition::matrixH()
328  *
329  * \tparam MatrixType type of matrix in the Hessenberg decomposition
330  *
331  * Objects of this type represent the Hessenberg matrix in the Hessenberg
332  * decomposition of some matrix. The object holds a reference to the
333  * HessenbergDecomposition class until the it is assigned or evaluated for
334  * some other reason (the reference should remain valid during the life time
335  * of this object). This class is the return type of
336  * HessenbergDecomposition::matrixH(); there is probably no other use for this
337  * class.
338  */
339template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
340: public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
341{
342  public:
343    /** \brief Constructor.
344      *
345      * \param[in] hess  Hessenberg decomposition
346      */
347    HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
348
349    /** \brief Hessenberg matrix in decomposition.
350      *
351      * \param[out] result  Hessenberg matrix in decomposition \p hess which
352      *                     was passed to the constructor
353      */
354    template <typename ResultType>
355    inline void evalTo(ResultType& result) const
356    {
357      result = m_hess.packedMatrix();
358      Index n = result.rows();
359      if (n>2)
360        result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
361    }
362
363    Index rows() const { return m_hess.packedMatrix().rows(); }
364    Index cols() const { return m_hess.packedMatrix().cols(); }
365
366  protected:
367    const HessenbergDecomposition<MatrixType>& m_hess;
368};
369
370} // end namespace internal
371
372} // end namespace Eigen
373
374#endif // EIGEN_HESSENBERGDECOMPOSITION_H
375