Tridiagonalization.h
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
14namespace Eigen {
15
16namespace internal {
17
18template<typename MatrixType> struct TridiagonalizationMatrixTReturnType;
19template<typename MatrixType>
20struct traits<TridiagonalizationMatrixTReturnType<MatrixType> >
21{
22 typedef typename MatrixType::PlainObject ReturnType;
23};
24
25template<typename MatrixType, typename CoeffVectorType>
26void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
27}
28
61template<typename _MatrixType> class Tridiagonalization
62{
63 public:
64
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
99 typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
100
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
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
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
180 inline CoeffVectorType householderCoefficients() const
181 {
182 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
183 return m_hCoeffs;
184 }
185
217 inline const MatrixType& packedMatrix() const
218 {
219 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
220 return m_matrix;
221 }
222
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
263 MatrixTReturnType matrixT() const
264 {
265 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
266 return MatrixTReturnType(m_matrix.real());
267 }
268
282 DiagonalReturnType diagonal() const;
283
294 SubDiagonalReturnType subDiagonal() const;
295
296 protected:
297
298 MatrixType m_matrix;
299 CoeffVectorType m_hCoeffs;
300 bool m_isInitialized;
301};
302
303template<typename MatrixType>
304typename Tridiagonalization<MatrixType>::DiagonalReturnType
305Tridiagonalization<MatrixType>::diagonal() const
306{
307 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
308 return m_matrix.diagonal();
309}
310
311template<typename MatrixType>
312typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
313Tridiagonalization<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
320namespace internal {
321
345template<typename MatrixType, typename CoeffVectorType>
346void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
347{
348 typedef typename MatrixType::Index Index;
349 typedef typename MatrixType::Scalar Scalar;
350 typedef typename MatrixType::RealScalar RealScalar;
351 Index n = matA.rows();
352 eigen_assert(n==matA.cols());
353 eigen_assert(n==hCoeffs.size()+1 || n==1);
354
355 for (Index i = 0; i<n-1; ++i)
356 {
357 Index remainingSize = n-i-1;
358 RealScalar beta;
359 Scalar h;
360 matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
361
362 // Apply similarity transformation to remaining columns,
363 // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
364 matA.col(i).coeffRef(i+1) = 1;
365
366 hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
367 * (conj(h) * matA.col(i).tail(remainingSize)));
368
369 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);
370
371 matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
372 .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
373
374 matA.col(i).coeffRef(i+1) = beta;
375 hCoeffs.coeffRef(i) = h;
376 }
377}
378
379// forward declaration, implementation at the end of this file
380template<typename MatrixType,
381 int Size=MatrixType::ColsAtCompileTime,
382 bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
383struct tridiagonalization_inplace_selector;
384
425template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
426void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
427{
428 typedef typename MatrixType::Index Index;
429 //Index n = mat.rows();
430 eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
431 tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
432}
433
437template<typename MatrixType, int Size, bool IsComplex>
438struct tridiagonalization_inplace_selector
439{
440 typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
441 typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
442 typedef typename MatrixType::Index Index;
443 template<typename DiagonalType, typename SubDiagonalType>
444 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
445 {
446 CoeffVectorType hCoeffs(mat.cols()-1);
447 tridiagonalization_inplace(mat,hCoeffs);
448 diag = mat.diagonal().real();
449 subdiag = mat.template diagonal<-1>().real();
450 if(extractQ)
451 mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
452 .setLength(mat.rows() - 1)
453 .setShift(1);
454 }
455};
456
461template<typename MatrixType>
462struct tridiagonalization_inplace_selector<MatrixType,3,false>
463{
464 typedef typename MatrixType::Scalar Scalar;
465 typedef typename MatrixType::RealScalar RealScalar;
466
467 template<typename DiagonalType, typename SubDiagonalType>
468 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
469 {
470 diag[0] = mat(0,0);
471 RealScalar v1norm2 = 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(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
505template<typename MatrixType, bool IsComplex>
506struct 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) = real(mat(0,0));
514 if(extractQ)
515 mat(0,0) = Scalar(1);
516 }
517};
518
526template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
527: public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
528{
529 typedef typename MatrixType::Index Index;
530 public:
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
Eigen::Block
Expression of a fixed-size or dynamic-size block.
Definition Block.h:99
Eigen::Diagonal
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition Diagonal.h:66
Eigen::HouseholderSequence::setShift
HouseholderSequence & setShift(Index shift)
Sets the shift of the Householder sequence.
Definition HouseholderSequence.h:363
Eigen::HouseholderSequence::setLength
HouseholderSequence & setLength(Index length)
Sets the length of the Householder sequence.
Definition HouseholderSequence.h:346
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition Matrix.h:129
Eigen::Tridiagonalization::packedMatrix
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
Definition Tridiagonalization.h:217
Eigen::Tridiagonalization::Tridiagonalization
Tridiagonalization(Index size=Size==Dynamic ? 2 :Size)
Default constructor.
Definition Tridiagonalization.h:113
Eigen::Tridiagonalization::diagonal
DiagonalReturnType diagonal() const
Returns the diagonal of the tridiagonal matrix T in the decomposition.
Definition Tridiagonalization.h:305
Eigen::Tridiagonalization::compute
Tridiagonalization & compute(const MatrixType &matrix)
Computes tridiagonal decomposition of given matrix.
Definition Tridiagonalization.h:155
Eigen::Tridiagonalization::subDiagonal
SubDiagonalReturnType subDiagonal() const
Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
Definition Tridiagonalization.h:313
Eigen::Tridiagonalization::Tridiagonalization
Tridiagonalization(const MatrixType &matrix)
Constructor; computes tridiagonal decomposition of given matrix.
Definition Tridiagonalization.h:129
Eigen::Tridiagonalization< MatrixType >::HouseholderSequenceType
HouseholderSequence< MatrixType, CoeffVectorType >::ConjugateReturnType HouseholderSequenceType
Definition Tridiagonalization.h:99
Eigen::Tridiagonalization::matrixT
MatrixTReturnType matrixT() const
Returns an expression of the tridiagonal matrix T in the decomposition.
Definition Tridiagonalization.h:263
Eigen::Tridiagonalization::householderCoefficients
CoeffVectorType householderCoefficients() const
Returns the Householder coefficients.
Definition Tridiagonalization.h:180
Eigen::Tridiagonalization< MatrixType >::MatrixType
MatrixType MatrixType
Definition Tridiagonalization.h:66
Eigen::Tridiagonalization::matrixQ
HouseholderSequenceType matrixQ() const
Returns the unitary matrix Q in the decomposition.
Definition Tridiagonalization.h:238
Eigen::internal
Definition LDLT.h:18