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