Main MRPT website > C++ reference
MRPT logo

Tridiagonalization.h

Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_TRIDIAGONALIZATION_H
00027 #define EIGEN_TRIDIAGONALIZATION_H
00028 
00029 namespace internal {
00030   
00031 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType;
00032 template<typename MatrixType>
00033 struct traits<TridiagonalizationMatrixTReturnType<MatrixType> >
00034 {
00035   typedef typename MatrixType::PlainObject ReturnType;
00036 };
00037 
00038 template<typename MatrixType, typename CoeffVectorType>
00039 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
00040 }
00041 
00042 /** \eigenvalues_module \ingroup Eigenvalues_Module
00043   *
00044   *
00045   * \class Tridiagonalization
00046   *
00047   * \brief Tridiagonal decomposition of a selfadjoint matrix
00048   *
00049   * \tparam _MatrixType the type of the matrix of which we are computing the
00050   * tridiagonal decomposition; this is expected to be an instantiation of the
00051   * Matrix class template.
00052   *
00053   * This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that:
00054   * \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real symmetric tridiagonal matrix.
00055   *
00056   * A tridiagonal matrix is a matrix which has nonzero elements only on the
00057   * main diagonal and the first diagonal below and above it. The Hessenberg
00058   * decomposition of a selfadjoint matrix is in fact a tridiagonal
00059   * decomposition. This class is used in SelfAdjointEigenSolver to compute the
00060   * eigenvalues and eigenvectors of a selfadjoint matrix.
00061   *
00062   * Call the function compute() to compute the tridiagonal decomposition of a
00063   * given matrix. Alternatively, you can use the Tridiagonalization(const MatrixType&)
00064   * constructor which computes the tridiagonal Schur decomposition at
00065   * construction time. Once the decomposition is computed, you can use the
00066   * matrixQ() and matrixT() functions to retrieve the matrices Q and T in the
00067   * decomposition.
00068   *
00069   * The documentation of Tridiagonalization(const MatrixType&) contains an
00070   * example of the typical use of this class.
00071   *
00072   * \sa class HessenbergDecomposition, class SelfAdjointEigenSolver
00073   */
00074 template<typename _MatrixType> class Tridiagonalization
00075 {
00076   public:
00077 
00078     /** \brief Synonym for the template parameter \p _MatrixType. */
00079     typedef _MatrixType MatrixType;
00080 
00081     typedef typename MatrixType::Scalar Scalar;
00082     typedef typename NumTraits<Scalar>::Real RealScalar;
00083     typedef typename MatrixType::Index Index;
00084 
00085     enum {
00086       Size = MatrixType::RowsAtCompileTime,
00087       SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
00088       Options = MatrixType::Options,
00089       MaxSize = MatrixType::MaxRowsAtCompileTime,
00090       MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
00091     };
00092 
00093     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
00094     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
00095     typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
00096     typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
00097     typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
00098 
00099     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
00100               const typename Diagonal<const MatrixType>::RealReturnType,
00101               const Diagonal<const MatrixType>
00102             >::type DiagonalReturnType;
00103 
00104     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
00105               const typename Diagonal<
00106                 Block<const MatrixType,SizeMinusOne,SizeMinusOne> >::RealReturnType,
00107               const Diagonal<
00108                 Block<const MatrixType,SizeMinusOne,SizeMinusOne> >
00109             >::type SubDiagonalReturnType;
00110 
00111     /** \brief Return type of matrixQ() */
00112     typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
00113 
00114     /** \brief Default constructor.
00115       *
00116       * \param [in]  size  Positive integer, size of the matrix whose tridiagonal
00117       * decomposition will be computed.
00118       *
00119       * The default constructor is useful in cases in which the user intends to
00120       * perform decompositions via compute().  The \p size parameter is only
00121       * used as a hint. It is not an error to give a wrong \p size, but it may
00122       * impair performance.
00123       *
00124       * \sa compute() for an example.
00125       */
00126     Tridiagonalization(Index size = Size==Dynamic ? 2 : Size)
00127       : m_matrix(size,size),
00128         m_hCoeffs(size > 1 ? size-1 : 1),
00129         m_isInitialized(false)
00130     {}
00131 
00132     /** \brief Constructor; computes tridiagonal decomposition of given matrix.
00133       *
00134       * \param[in]  matrix  Selfadjoint matrix whose tridiagonal decomposition
00135       * is to be computed.
00136       *
00137       * This constructor calls compute() to compute the tridiagonal decomposition.
00138       *
00139       * Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp
00140       * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out
00141       */
00142     Tridiagonalization(const MatrixType& matrix)
00143       : m_matrix(matrix),
00144         m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
00145         m_isInitialized(false)
00146     {
00147       internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
00148       m_isInitialized = true;
00149     }
00150 
00151     /** \brief Computes tridiagonal decomposition of given matrix.
00152       *
00153       * \param[in]  matrix  Selfadjoint matrix whose tridiagonal decomposition
00154       * is to be computed.
00155       * \returns    Reference to \c *this
00156       *
00157       * The tridiagonal decomposition is computed by bringing the columns of
00158       * the matrix successively in the required form using Householder
00159       * reflections. The cost is \f$ 4n^3/3 \f$ flops, where \f$ n \f$ denotes
00160       * the size of the given matrix.
00161       *
00162       * This method reuses of the allocated data in the Tridiagonalization
00163       * object, if the size of the matrix does not change.
00164       *
00165       * Example: \include Tridiagonalization_compute.cpp
00166       * Output: \verbinclude Tridiagonalization_compute.out
00167       */
00168     Tridiagonalization& compute(const MatrixType& matrix)
00169     {
00170       m_matrix = matrix;
00171       m_hCoeffs.resize(matrix.rows()-1, 1);
00172       internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
00173       m_isInitialized = true;
00174       return *this;
00175     }
00176 
00177     /** \brief Returns the Householder coefficients.
00178       *
00179       * \returns a const reference to the vector of Householder coefficients
00180       *
00181       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
00182       * the member function compute(const MatrixType&) has been called before
00183       * to compute the tridiagonal decomposition of a matrix.
00184       *
00185       * The Householder coefficients allow the reconstruction of the matrix
00186       * \f$ Q \f$ in the tridiagonal decomposition from the packed data.
00187       *
00188       * Example: \include Tridiagonalization_householderCoefficients.cpp
00189       * Output: \verbinclude Tridiagonalization_householderCoefficients.out
00190       *
00191       * \sa packedMatrix(), \ref Householder_Module "Householder module"
00192       */
00193     inline CoeffVectorType householderCoefficients() const
00194     {
00195       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00196       return m_hCoeffs;
00197     }
00198 
00199     /** \brief Returns the internal representation of the decomposition
00200       *
00201       * \returns a const reference to a matrix with the internal representation
00202       *          of the decomposition.
00203       *
00204       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
00205       * the member function compute(const MatrixType&) has been called before
00206       * to compute the tridiagonal decomposition of a matrix.
00207       *
00208       * The returned matrix contains the following information:
00209       *  - the strict upper triangular part is equal to the input matrix A.
00210       *  - the diagonal and lower sub-diagonal represent the real tridiagonal
00211       *    symmetric matrix T.
00212       *  - the rest of the lower part contains the Householder vectors that,
00213       *    combined with Householder coefficients returned by
00214       *    householderCoefficients(), allows to reconstruct the matrix Q as
00215       *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
00216       *    Here, the matrices \f$ H_i \f$ are the Householder transformations
00217       *       \f$ H_i = (I - h_i v_i v_i^T) \f$
00218       *    where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
00219       *    \f$ v_i \f$ is the Householder vector defined by
00220       *       \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
00221       *    with M the matrix returned by this function.
00222       *
00223       * See LAPACK for further details on this packed storage.
00224       *
00225       * Example: \include Tridiagonalization_packedMatrix.cpp
00226       * Output: \verbinclude Tridiagonalization_packedMatrix.out
00227       *
00228       * \sa householderCoefficients()
00229       */
00230     inline const MatrixType& packedMatrix() const
00231     {
00232       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00233       return m_matrix;
00234     }
00235 
00236     /** \brief Returns the unitary matrix Q in the decomposition
00237       *
00238       * \returns object representing the matrix Q
00239       *
00240       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
00241       * the member function compute(const MatrixType&) has been called before
00242       * to compute the tridiagonal decomposition of a matrix.
00243       *
00244       * This function returns a light-weight object of template class
00245       * HouseholderSequence. You can either apply it directly to a matrix or
00246       * you can convert it to a matrix of type #MatrixType.
00247       *
00248       * \sa Tridiagonalization(const MatrixType&) for an example,
00249       *     matrixT(), class HouseholderSequence
00250       */
00251     HouseholderSequenceType matrixQ() const
00252     {
00253       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00254       return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
00255              .setLength(m_matrix.rows() - 1)
00256              .setShift(1);
00257     }
00258 
00259     /** \brief Returns an expression of the tridiagonal matrix T in the decomposition
00260       *
00261       * \returns expression object representing the matrix T
00262       *
00263       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
00264       * the member function compute(const MatrixType&) has been called before
00265       * to compute the tridiagonal decomposition of a matrix.
00266       *
00267       * Currently, this function can be used to extract the matrix T from internal
00268       * data and copy it to a dense matrix object. In most cases, it may be
00269       * sufficient to directly use the packed matrix or the vector expressions
00270       * returned by diagonal() and subDiagonal() instead of creating a new
00271       * dense copy matrix with this function.
00272       *
00273       * \sa Tridiagonalization(const MatrixType&) for an example,
00274       * matrixQ(), packedMatrix(), diagonal(), subDiagonal()
00275       */
00276     MatrixTReturnType matrixT() const
00277     {
00278       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00279       return MatrixTReturnType(m_matrix.real());
00280     }
00281 
00282     /** \brief Returns the diagonal of the tridiagonal matrix T in the decomposition.
00283       *
00284       * \returns expression representing the diagonal of T
00285       *
00286       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
00287       * the member function compute(const MatrixType&) has been called before
00288       * to compute the tridiagonal decomposition of a matrix.
00289       *
00290       * Example: \include Tridiagonalization_diagonal.cpp
00291       * Output: \verbinclude Tridiagonalization_diagonal.out
00292       *
00293       * \sa matrixT(), subDiagonal()
00294       */
00295     DiagonalReturnType diagonal() const;
00296 
00297     /** \brief Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
00298       *
00299       * \returns expression representing the subdiagonal of T
00300       *
00301       * \pre Either the constructor Tridiagonalization(const MatrixType&) or
00302       * the member function compute(const MatrixType&) has been called before
00303       * to compute the tridiagonal decomposition of a matrix.
00304       *
00305       * \sa diagonal() for an example, matrixT()
00306       */
00307     SubDiagonalReturnType subDiagonal() const;
00308 
00309   protected:
00310 
00311     MatrixType m_matrix;
00312     CoeffVectorType m_hCoeffs;
00313     bool m_isInitialized;
00314 };
00315 
00316 template<typename MatrixType>
00317 typename Tridiagonalization<MatrixType>::DiagonalReturnType
00318 Tridiagonalization<MatrixType>::diagonal() const
00319 {
00320   eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00321   return m_matrix.diagonal();
00322 }
00323 
00324 template<typename MatrixType>
00325 typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
00326 Tridiagonalization<MatrixType>::subDiagonal() const
00327 {
00328   eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00329   Index n = m_matrix.rows();
00330   return Block<const MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal();
00331 }
00332 
00333 namespace internal {
00334 
00335 /** \internal
00336   * Performs a tridiagonal decomposition of the selfadjoint matrix \a matA in-place.
00337   *
00338   * \param[in,out] matA On input the selfadjoint matrix. Only the \b lower triangular part is referenced.
00339   *                     On output, the strict upper part is left unchanged, and the lower triangular part
00340   *                     represents the T and Q matrices in packed format has detailed below.
00341   * \param[out]    hCoeffs returned Householder coefficients (see below)
00342   *
00343   * On output, the tridiagonal selfadjoint matrix T is stored in the diagonal
00344   * and lower sub-diagonal of the matrix \a matA.
00345   * The unitary matrix Q is represented in a compact way as a product of
00346   * Householder reflectors \f$ H_i \f$ such that:
00347   *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
00348   * The Householder reflectors are defined as
00349   *       \f$ H_i = (I - h_i v_i v_i^T) \f$
00350   * where \f$ h_i = hCoeffs[i]\f$ is the \f$ i \f$th Householder coefficient and
00351   * \f$ v_i \f$ is the Householder vector defined by
00352   *       \f$ v_i = [ 0, \ldots, 0, 1, matA(i+2,i), \ldots, matA(N-1,i) ]^T \f$.
00353   *
00354   * Implemented from Golub's "Matrix Computations", algorithm 8.3.1.
00355   *
00356   * \sa Tridiagonalization::packedMatrix()
00357   */
00358 template<typename MatrixType, typename CoeffVectorType>
00359 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
00360 {
00361   typedef typename MatrixType::Index Index;
00362   typedef typename MatrixType::Scalar Scalar;
00363   typedef typename MatrixType::RealScalar RealScalar;
00364   Index n = matA.rows();
00365   eigen_assert(n==matA.cols());
00366   eigen_assert(n==hCoeffs.size()+1 || n==1);
00367   
00368   for (Index i = 0; i<n-1; ++i)
00369   {
00370     Index remainingSize = n-i-1;
00371     RealScalar beta;
00372     Scalar h;
00373     matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
00374 
00375     // Apply similarity transformation to remaining columns,
00376     // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
00377     matA.col(i).coeffRef(i+1) = 1;
00378 
00379     hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
00380                                   * (conj(h) * matA.col(i).tail(remainingSize)));
00381 
00382     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);
00383 
00384     matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
00385       .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
00386 
00387     matA.col(i).coeffRef(i+1) = beta;
00388     hCoeffs.coeffRef(i) = h;
00389   }
00390 }
00391 
00392 // forward declaration, implementation at the end of this file
00393 template<typename MatrixType,
00394          int Size=MatrixType::ColsAtCompileTime,
00395          bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
00396 struct tridiagonalization_inplace_selector;
00397 
00398 /** \brief Performs a full tridiagonalization in place
00399   *
00400   * \param[in,out]  mat  On input, the selfadjoint matrix whose tridiagonal
00401   *    decomposition is to be computed. Only the lower triangular part referenced.
00402   *    The rest is left unchanged. On output, the orthogonal matrix Q
00403   *    in the decomposition if \p extractQ is true.
00404   * \param[out]  diag  The diagonal of the tridiagonal matrix T in the
00405   *    decomposition.
00406   * \param[out]  subdiag  The subdiagonal of the tridiagonal matrix T in
00407   *    the decomposition.
00408   * \param[in]  extractQ  If true, the orthogonal matrix Q in the
00409   *    decomposition is computed and stored in \p mat.
00410   *
00411   * Computes the tridiagonal decomposition of the selfadjoint matrix \p mat in place
00412   * such that \f$ mat = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real
00413   * symmetric tridiagonal matrix.
00414   *
00415   * The tridiagonal matrix T is passed to the output parameters \p diag and \p subdiag. If
00416   * \p extractQ is true, then the orthogonal matrix Q is passed to \p mat. Otherwise the lower
00417   * part of the matrix \p mat is destroyed.
00418   *
00419   * The vectors \p diag and \p subdiag are not resized. The function
00420   * assumes that they are already of the correct size. The length of the
00421   * vector \p diag should equal the number of rows in \p mat, and the
00422   * length of the vector \p subdiag should be one left.
00423   *
00424   * This implementation contains an optimized path for 3-by-3 matrices
00425   * which is especially useful for plane fitting.
00426   *
00427   * \note Currently, it requires two temporary vectors to hold the intermediate
00428   * Householder coefficients, and to reconstruct the matrix Q from the Householder
00429   * reflectors.
00430   *
00431   * Example (this uses the same matrix as the example in
00432   *    Tridiagonalization::Tridiagonalization(const MatrixType&)):
00433   *    \include Tridiagonalization_decomposeInPlace.cpp
00434   * Output: \verbinclude Tridiagonalization_decomposeInPlace.out
00435   *
00436   * \sa class Tridiagonalization
00437   */
00438 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
00439 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00440 {
00441   typedef typename MatrixType::Index Index;
00442   //Index n = mat.rows();
00443   eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
00444   tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
00445 }
00446 
00447 /** \internal
00448   * General full tridiagonalization
00449   */
00450 template<typename MatrixType, int Size, bool IsComplex>
00451 struct tridiagonalization_inplace_selector
00452 {
00453   typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
00454   typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
00455   typedef typename MatrixType::Index Index;
00456   template<typename DiagonalType, typename SubDiagonalType>
00457   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00458   {
00459     CoeffVectorType hCoeffs(mat.cols()-1);
00460     tridiagonalization_inplace(mat,hCoeffs);
00461     diag = mat.diagonal().real();
00462     subdiag = mat.template diagonal<-1>().real();
00463     if(extractQ)
00464       mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
00465             .setLength(mat.rows() - 1)
00466             .setShift(1);
00467   }
00468 };
00469 
00470 /** \internal
00471   * Specialization for 3x3 real matrices.
00472   * Especially useful for plane fitting.
00473   */
00474 template<typename MatrixType>
00475 struct tridiagonalization_inplace_selector<MatrixType,3,false>
00476 {
00477   typedef typename MatrixType::Scalar Scalar;
00478   typedef typename MatrixType::RealScalar RealScalar;
00479 
00480   template<typename DiagonalType, typename SubDiagonalType>
00481   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00482   {
00483     diag[0] = mat(0,0);
00484     RealScalar v1norm2 = abs2(mat(2,0));
00485     if(v1norm2 == RealScalar(0))
00486     {
00487       diag[1] = mat(1,1);
00488       diag[2] = mat(2,2);
00489       subdiag[0] = mat(1,0);
00490       subdiag[1] = mat(2,1);
00491       if (extractQ)
00492         mat.setIdentity();
00493     }
00494     else
00495     {
00496       RealScalar beta = sqrt(abs2(mat(1,0)) + v1norm2);
00497       RealScalar invBeta = RealScalar(1)/beta;
00498       Scalar m01 = mat(1,0) * invBeta;
00499       Scalar m02 = mat(2,0) * invBeta;
00500       Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
00501       diag[1] = mat(1,1) + m02*q;
00502       diag[2] = mat(2,2) - m02*q;
00503       subdiag[0] = beta;
00504       subdiag[1] = mat(2,1) - m01 * q;
00505       if (extractQ)
00506       {
00507         mat << 1,   0,    0,
00508                0, m01,  m02,
00509                0, m02, -m01;
00510       }
00511     }
00512   }
00513 };
00514 
00515 /** \internal
00516   * Trivial specialization for 1x1 matrices
00517   */
00518 template<typename MatrixType, bool IsComplex>
00519 struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
00520 {
00521   typedef typename MatrixType::Scalar Scalar;
00522 
00523   template<typename DiagonalType, typename SubDiagonalType>
00524   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ)
00525   {
00526     diag(0,0) = real(mat(0,0));
00527     if(extractQ)
00528       mat(0,0) = Scalar(1);
00529   }
00530 };
00531 
00532 /** \internal
00533   * \eigenvalues_module \ingroup Eigenvalues_Module
00534   *
00535   * \brief Expression type for return value of Tridiagonalization::matrixT()
00536   *
00537   * \tparam MatrixType type of underlying dense matrix
00538   */
00539 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
00540 : public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
00541 {
00542     typedef typename MatrixType::Index Index;
00543   public:
00544     /** \brief Constructor.
00545       *
00546       * \param[in] mat The underlying dense matrix
00547       */
00548     TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { }
00549 
00550     template <typename ResultType>
00551     inline void evalTo(ResultType& result) const
00552     {
00553       result.setZero();
00554       result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
00555       result.diagonal() = m_matrix.diagonal();
00556       result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
00557     }
00558 
00559     Index rows() const { return m_matrix.rows(); }
00560     Index cols() const { return m_matrix.cols(); }
00561 
00562   protected:
00563     const typename MatrixType::Nested m_matrix;
00564 };
00565 
00566 } // end namespace internal
00567 
00568 #endif // EIGEN_TRIDIAGONALIZATION_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:16:28 UTC 2011