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.2 for MRPT 0.9.4 SVN: at Mon Jan 10 22:46:17 UTC 2011 |