Home | History | Annotate | Download | only in Eigenvalues
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      5 // Copyright (C) 2010 Jitse Niesen <jitse (at) 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_TRIDIAGONALIZATION_H
     12 #define EIGEN_TRIDIAGONALIZATION_H
     13 
     14 namespace Eigen {
     15 
     16 namespace internal {
     17 
     18 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType;
     19 template<typename MatrixType>
     20 struct traits<TridiagonalizationMatrixTReturnType<MatrixType> >
     21 {
     22   typedef typename MatrixType::PlainObject ReturnType;
     23 };
     24 
     25 template<typename MatrixType, typename CoeffVectorType>
     26 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
     27 }
     28 
     29 /** \eigenvalues_module \ingroup Eigenvalues_Module
     30   *
     31   *
     32   * \class Tridiagonalization
     33   *
     34   * \brief Tridiagonal decomposition of a selfadjoint matrix
     35   *
     36   * \tparam _MatrixType the type of the matrix of which we are computing the
     37   * tridiagonal decomposition; this is expected to be an instantiation of the
     38   * Matrix class template.
     39   *
     40   * This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that:
     41   * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real symmetric tridiagonal matrix.
     42   *
     43   * A tridiagonal matrix is a matrix which has nonzero elements only on the
     44   * main diagonal and the first diagonal below and above it. The Hessenberg
     45   * decomposition of a selfadjoint matrix is in fact a tridiagonal
     46   * decomposition. This class is used in SelfAdjointEigenSolver to compute the
     47   * eigenvalues and eigenvectors of a selfadjoint matrix.
     48   *
     49   * Call the function compute() to compute the tridiagonal decomposition of a
     50   * given matrix. Alternatively, you can use the Tridiagonalization(const MatrixType&)
     51   * constructor which computes the tridiagonal Schur decomposition at
     52   * construction time. Once the decomposition is computed, you can use the
     53   * matrixQ() and matrixT() functions to retrieve the matrices Q and T in the
     54   * decomposition.
     55   *
     56   * The documentation of Tridiagonalization(const MatrixType&) contains an
     57   * example of the typical use of this class.
     58   *
     59   * \sa class HessenbergDecomposition, class SelfAdjointEigenSolver
     60   */
     61 template<typename _MatrixType> class Tridiagonalization
     62 {
     63   public:
     64 
     65     /** \brief Synonym for the template parameter \p _MatrixType. */
     66     typedef _MatrixType MatrixType;
     67 
     68     typedef typename MatrixType::Scalar Scalar;
     69     typedef typename NumTraits<Scalar>::Real RealScalar;
     70     typedef typename MatrixType::Index Index;
     71 
     72     enum {
     73       Size = MatrixType::RowsAtCompileTime,
     74       SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
     75       Options = MatrixType::Options,
     76       MaxSize = MatrixType::MaxRowsAtCompileTime,
     77       MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
     78     };
     79 
     80     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
     81     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
     82     typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
     83     typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
     84     typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
     85 
     86     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
     87               typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
     88               const Diagonal<const MatrixType>
     89             >::type DiagonalReturnType;
     90 
     91     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
     92               typename internal::add_const_on_value_type<typename Diagonal<
     93                 Block<const MatrixType,SizeMinusOne,SizeMinusOne> >::RealReturnType>::type,
     94               const Diagonal<
     95                 Block<const MatrixType,SizeMinusOne,SizeMinusOne> >
     96             >::type SubDiagonalReturnType;
     97 
     98     /** \brief Return type of matrixQ() */
     99     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType;
    100 
    101     /** \brief Default constructor.
    102       *
    103       * \param [in]  size  Positive integer, size of the matrix whose tridiagonal
    104       * decomposition will be computed.
    105       *
    106       * The default constructor is useful in cases in which the user intends to
    107       * perform decompositions via compute().  The \p size parameter is only
    108       * used as a hint. It is not an error to give a wrong \p size, but it may
    109       * impair performance.
    110       *
    111       * \sa compute() for an example.
    112       */
    113     Tridiagonalization(Index size = Size==Dynamic ? 2 : Size)
    114       : m_matrix(size,size),
    115         m_hCoeffs(size > 1 ? size-1 : 1),
    116         m_isInitialized(false)
    117     {}
    118 
    119     /** \brief Constructor; computes tridiagonal decomposition of given matrix.
    120       *
    121       * \param[in]  matrix  Selfadjoint matrix whose tridiagonal decomposition
    122       * is to be computed.
    123       *
    124       * This constructor calls compute() to compute the tridiagonal decomposition.
    125       *
    126       * Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp
    127       * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out
    128       */
    129     Tridiagonalization(const MatrixType& matrix)
    130       : m_matrix(matrix),
    131         m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
    132         m_isInitialized(false)
    133     {
    134       internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
    135       m_isInitialized = true;
    136     }
    137 
    138     /** \brief Computes tridiagonal decomposition of given matrix.
    139       *
    140       * \param[in]  matrix  Selfadjoint matrix whose tridiagonal decomposition
    141       * is to be computed.
    142       * \returns    Reference to \c *this
    143       *
    144       * The tridiagonal decomposition is computed by bringing the columns of
    145       * the matrix successively in the required form using Householder
    146       * reflections. The cost is \f$ 4n^3/3 \f$ flops, where \f$ n \f$ denotes
    147       * the size of the given matrix.
    148       *
    149       * This method reuses of the allocated data in the Tridiagonalization
    150       * object, if the size of the matrix does not change.
    151       *
    152       * Example: \include Tridiagonalization_compute.cpp
    153       * Output: \verbinclude Tridiagonalization_compute.out
    154       */
    155     Tridiagonalization& compute(const MatrixType& matrix)
    156     {
    157       m_matrix = matrix;
    158       m_hCoeffs.resize(matrix.rows()-1, 1);
    159       internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
    160       m_isInitialized = true;
    161       return *this;
    162     }
    163 
    164     /** \brief Returns the Householder coefficients.
    165       *
    166       * \returns a const reference to the vector of Householder coefficients
    167       *
    168       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
    169       * the member function compute(const MatrixType&) has been called before
    170       * to compute the tridiagonal decomposition of a matrix.
    171       *
    172       * The Householder coefficients allow the reconstruction of the matrix
    173       * \f$ Q \f$ in the tridiagonal decomposition from the packed data.
    174       *
    175       * Example: \include Tridiagonalization_householderCoefficients.cpp
    176       * Output: \verbinclude Tridiagonalization_householderCoefficients.out
    177       *
    178       * \sa packedMatrix(), \ref Householder_Module "Householder module"
    179       */
    180     inline CoeffVectorType householderCoefficients() const
    181     {
    182       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
    183       return m_hCoeffs;
    184     }
    185 
    186     /** \brief Returns the internal representation of the decomposition
    187       *
    188       *	\returns a const reference to a matrix with the internal representation
    189       *	         of the decomposition.
    190       *
    191       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
    192       * the member function compute(const MatrixType&) has been called before
    193       * to compute the tridiagonal decomposition of a matrix.
    194       *
    195       * The returned matrix contains the following information:
    196       *  - the strict upper triangular part is equal to the input matrix A.
    197       *  - the diagonal and lower sub-diagonal represent the real tridiagonal
    198       *    symmetric matrix T.
    199       *  - the rest of the lower part contains the Householder vectors that,
    200       *    combined with Householder coefficients returned by
    201       *    householderCoefficients(), allows to reconstruct the matrix Q as
    202       *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
    203       *    Here, the matrices \f$ H_i \f$ are the Householder transformations
    204       *       \f$ H_i = (I - h_i v_i v_i^T) \f$
    205       *    where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
    206       *    \f$ v_i \f$ is the Householder vector defined by
    207       *       \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
    208       *    with M the matrix returned by this function.
    209       *
    210       * See LAPACK for further details on this packed storage.
    211       *
    212       * Example: \include Tridiagonalization_packedMatrix.cpp
    213       * Output: \verbinclude Tridiagonalization_packedMatrix.out
    214       *
    215       * \sa householderCoefficients()
    216       */
    217     inline const MatrixType& packedMatrix() const
    218     {
    219       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
    220       return m_matrix;
    221     }
    222 
    223     /** \brief Returns the unitary matrix Q in the decomposition
    224       *
    225       * \returns object representing the matrix Q
    226       *
    227       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
    228       * the member function compute(const MatrixType&) has been called before
    229       * to compute the tridiagonal decomposition of a matrix.
    230       *
    231       * This function returns a light-weight object of template class
    232       * HouseholderSequence. You can either apply it directly to a matrix or
    233       * you can convert it to a matrix of type #MatrixType.
    234       *
    235       * \sa Tridiagonalization(const MatrixType&) for an example,
    236       *     matrixT(), class HouseholderSequence
    237       */
    238     HouseholderSequenceType matrixQ() const
    239     {
    240       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
    241       return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
    242              .setLength(m_matrix.rows() - 1)
    243              .setShift(1);
    244     }
    245 
    246     /** \brief Returns an expression of the tridiagonal matrix T in the decomposition
    247       *
    248       * \returns expression object representing the matrix T
    249       *
    250       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
    251       * the member function compute(const MatrixType&) has been called before
    252       * to compute the tridiagonal decomposition of a matrix.
    253       *
    254       * Currently, this function can be used to extract the matrix T from internal
    255       * data and copy it to a dense matrix object. In most cases, it may be
    256       * sufficient to directly use the packed matrix or the vector expressions
    257       * returned by diagonal() and subDiagonal() instead of creating a new
    258       * dense copy matrix with this function.
    259       *
    260       * \sa Tridiagonalization(const MatrixType&) for an example,
    261       * matrixQ(), packedMatrix(), diagonal(), subDiagonal()
    262       */
    263     MatrixTReturnType matrixT() const
    264     {
    265       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
    266       return MatrixTReturnType(m_matrix.real());
    267     }
    268 
    269     /** \brief Returns the diagonal of the tridiagonal matrix T in the decomposition.
    270       *
    271       * \returns expression representing the diagonal of T
    272       *
    273       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
    274       * the member function compute(const MatrixType&) has been called before
    275       * to compute the tridiagonal decomposition of a matrix.
    276       *
    277       * Example: \include Tridiagonalization_diagonal.cpp
    278       * Output: \verbinclude Tridiagonalization_diagonal.out
    279       *
    280       * \sa matrixT(), subDiagonal()
    281       */
    282     DiagonalReturnType diagonal() const;
    283 
    284     /** \brief Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
    285       *
    286       * \returns expression representing the subdiagonal of T
    287       *
    288       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
    289       * the member function compute(const MatrixType&) has been called before
    290       * to compute the tridiagonal decomposition of a matrix.
    291       *
    292       * \sa diagonal() for an example, matrixT()
    293       */
    294     SubDiagonalReturnType subDiagonal() const;
    295 
    296   protected:
    297 
    298     MatrixType m_matrix;
    299     CoeffVectorType m_hCoeffs;
    300     bool m_isInitialized;
    301 };
    302 
    303 template<typename MatrixType>
    304 typename Tridiagonalization<MatrixType>::DiagonalReturnType
    305 Tridiagonalization<MatrixType>::diagonal() const
    306 {
    307   eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
    308   return m_matrix.diagonal();
    309 }
    310 
    311 template<typename MatrixType>
    312 typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
    313 Tridiagonalization<MatrixType>::subDiagonal() const
    314 {
    315   eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
    316   Index n = m_matrix.rows();
    317   return Block<const MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal();
    318 }
    319 
    320 namespace internal {
    321 
    322 /** \internal
    323   * Performs a tridiagonal decomposition of the selfadjoint matrix \a matA in-place.
    324   *
    325   * \param[in,out] matA On input the selfadjoint matrix. Only the \b lower triangular part is referenced.
    326   *                     On output, the strict upper part is left unchanged, and the lower triangular part
    327   *                     represents the T and Q matrices in packed format has detailed below.
    328   * \param[out]    hCoeffs returned Householder coefficients (see below)
    329   *
    330   * On output, the tridiagonal selfadjoint matrix T is stored in the diagonal
    331   * and lower sub-diagonal of the matrix \a matA.
    332   * The unitary matrix Q is represented in a compact way as a product of
    333   * Householder reflectors \f$ H_i \f$ such that:
    334   *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
    335   * The Householder reflectors are defined as
    336   *       \f$ H_i = (I - h_i v_i v_i^T) \f$
    337   * where \f$ h_i = hCoeffs[i]\f$ is the \f$ i \f$th Householder coefficient and
    338   * \f$ v_i \f$ is the Householder vector defined by
    339   *       \f$ v_i = [ 0, \ldots, 0, 1, matA(i+2,i), \ldots, matA(N-1,i) ]^T \f$.
    340   *
    341   * Implemented from Golub's "Matrix Computations", algorithm 8.3.1.
    342   *
    343   * \sa Tridiagonalization::packedMatrix()
    344   */
    345 template<typename MatrixType, typename CoeffVectorType>
    346 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
    347 {
    348   using numext::conj;
    349   typedef typename MatrixType::Index Index;
    350   typedef typename MatrixType::Scalar Scalar;
    351   typedef typename MatrixType::RealScalar RealScalar;
    352   Index n = matA.rows();
    353   eigen_assert(n==matA.cols());
    354   eigen_assert(n==hCoeffs.size()+1 || n==1);
    355 
    356   for (Index i = 0; i<n-1; ++i)
    357   {
    358     Index remainingSize = n-i-1;
    359     RealScalar beta;
    360     Scalar h;
    361     matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
    362 
    363     // Apply similarity transformation to remaining columns,
    364     // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
    365     matA.col(i).coeffRef(i+1) = 1;
    366 
    367     hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
    368                                   * (conj(h) * matA.col(i).tail(remainingSize)));
    369 
    370     hCoeffs.tail(n-i-1) += (conj(h)*Scalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1);
    371 
    372     matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
    373       .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
    374 
    375     matA.col(i).coeffRef(i+1) = beta;
    376     hCoeffs.coeffRef(i) = h;
    377   }
    378 }
    379 
    380 // forward declaration, implementation at the end of this file
    381 template<typename MatrixType,
    382          int Size=MatrixType::ColsAtCompileTime,
    383          bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
    384 struct tridiagonalization_inplace_selector;
    385 
    386 /** \brief Performs a full tridiagonalization in place
    387   *
    388   * \param[in,out]  mat  On input, the selfadjoint matrix whose tridiagonal
    389   *    decomposition is to be computed. Only the lower triangular part referenced.
    390   *    The rest is left unchanged. On output, the orthogonal matrix Q
    391   *    in the decomposition if \p extractQ is true.
    392   * \param[out]  diag  The diagonal of the tridiagonal matrix T in the
    393   *    decomposition.
    394   * \param[out]  subdiag  The subdiagonal of the tridiagonal matrix T in
    395   *    the decomposition.
    396   * \param[in]  extractQ  If true, the orthogonal matrix Q in the
    397   *    decomposition is computed and stored in \p mat.
    398   *
    399   * Computes the tridiagonal decomposition of the selfadjoint matrix \p mat in place
    400   * such that \f$ mat = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real
    401   * symmetric tridiagonal matrix.
    402   *
    403   * The tridiagonal matrix T is passed to the output parameters \p diag and \p subdiag. If
    404   * \p extractQ is true, then the orthogonal matrix Q is passed to \p mat. Otherwise the lower
    405   * part of the matrix \p mat is destroyed.
    406   *
    407   * The vectors \p diag and \p subdiag are not resized. The function
    408   * assumes that they are already of the correct size. The length of the
    409   * vector \p diag should equal the number of rows in \p mat, and the
    410   * length of the vector \p subdiag should be one left.
    411   *
    412   * This implementation contains an optimized path for 3-by-3 matrices
    413   * which is especially useful for plane fitting.
    414   *
    415   * \note Currently, it requires two temporary vectors to hold the intermediate
    416   * Householder coefficients, and to reconstruct the matrix Q from the Householder
    417   * reflectors.
    418   *
    419   * Example (this uses the same matrix as the example in
    420   *    Tridiagonalization::Tridiagonalization(const MatrixType&)):
    421   *    \include Tridiagonalization_decomposeInPlace.cpp
    422   * Output: \verbinclude Tridiagonalization_decomposeInPlace.out
    423   *
    424   * \sa class Tridiagonalization
    425   */
    426 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
    427 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
    428 {
    429   eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
    430   tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
    431 }
    432 
    433 /** \internal
    434   * General full tridiagonalization
    435   */
    436 template<typename MatrixType, int Size, bool IsComplex>
    437 struct tridiagonalization_inplace_selector
    438 {
    439   typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
    440   typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
    441   typedef typename MatrixType::Index Index;
    442   template<typename DiagonalType, typename SubDiagonalType>
    443   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
    444   {
    445     CoeffVectorType hCoeffs(mat.cols()-1);
    446     tridiagonalization_inplace(mat,hCoeffs);
    447     diag = mat.diagonal().real();
    448     subdiag = mat.template diagonal<-1>().real();
    449     if(extractQ)
    450       mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
    451             .setLength(mat.rows() - 1)
    452             .setShift(1);
    453   }
    454 };
    455 
    456 /** \internal
    457   * Specialization for 3x3 real matrices.
    458   * Especially useful for plane fitting.
    459   */
    460 template<typename MatrixType>
    461 struct tridiagonalization_inplace_selector<MatrixType,3,false>
    462 {
    463   typedef typename MatrixType::Scalar Scalar;
    464   typedef typename MatrixType::RealScalar RealScalar;
    465 
    466   template<typename DiagonalType, typename SubDiagonalType>
    467   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
    468   {
    469     using std::sqrt;
    470     diag[0] = mat(0,0);
    471     RealScalar v1norm2 = numext::abs2(mat(2,0));
    472     if(v1norm2 == RealScalar(0))
    473     {
    474       diag[1] = mat(1,1);
    475       diag[2] = mat(2,2);
    476       subdiag[0] = mat(1,0);
    477       subdiag[1] = mat(2,1);
    478       if (extractQ)
    479         mat.setIdentity();
    480     }
    481     else
    482     {
    483       RealScalar beta = sqrt(numext::abs2(mat(1,0)) + v1norm2);
    484       RealScalar invBeta = RealScalar(1)/beta;
    485       Scalar m01 = mat(1,0) * invBeta;
    486       Scalar m02 = mat(2,0) * invBeta;
    487       Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
    488       diag[1] = mat(1,1) + m02*q;
    489       diag[2] = mat(2,2) - m02*q;
    490       subdiag[0] = beta;
    491       subdiag[1] = mat(2,1) - m01 * q;
    492       if (extractQ)
    493       {
    494         mat << 1,   0,    0,
    495                0, m01,  m02,
    496                0, m02, -m01;
    497       }
    498     }
    499   }
    500 };
    501 
    502 /** \internal
    503   * Trivial specialization for 1x1 matrices
    504   */
    505 template<typename MatrixType, bool IsComplex>
    506 struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
    507 {
    508   typedef typename MatrixType::Scalar Scalar;
    509 
    510   template<typename DiagonalType, typename SubDiagonalType>
    511   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ)
    512   {
    513     diag(0,0) = numext::real(mat(0,0));
    514     if(extractQ)
    515       mat(0,0) = Scalar(1);
    516   }
    517 };
    518 
    519 /** \internal
    520   * \eigenvalues_module \ingroup Eigenvalues_Module
    521   *
    522   * \brief Expression type for return value of Tridiagonalization::matrixT()
    523   *
    524   * \tparam MatrixType type of underlying dense matrix
    525   */
    526 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
    527 : public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
    528 {
    529     typedef typename MatrixType::Index Index;
    530   public:
    531     /** \brief Constructor.
    532       *
    533       * \param[in] mat The underlying dense matrix
    534       */
    535     TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { }
    536 
    537     template <typename ResultType>
    538     inline void evalTo(ResultType& result) const
    539     {
    540       result.setZero();
    541       result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
    542       result.diagonal() = m_matrix.diagonal();
    543       result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
    544     }
    545 
    546     Index rows() const { return m_matrix.rows(); }
    547     Index cols() const { return m_matrix.cols(); }
    548 
    549   protected:
    550     typename MatrixType::Nested m_matrix;
    551 };
    552 
    553 } // end namespace internal
    554 
    555 } // end namespace Eigen
    556 
    557 #endif // EIGEN_TRIDIAGONALIZATION_H
    558