Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
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 : public traits<typename MatrixType::PlainObject>
22{
23 typedef typename MatrixType::PlainObject ReturnType; // FIXME shall it be a BandMatrix?
24 enum { Flags = 0 };
25};
26
27template<typename MatrixType, typename CoeffVectorType>
28EIGEN_DEVICE_FUNC
29void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
30}
31
64template<typename _MatrixType> class Tridiagonalization
65{
66 public:
67
69 typedef _MatrixType MatrixType;
70
71 typedef typename MatrixType::Scalar Scalar;
72 typedef typename NumTraits<Scalar>::Real RealScalar;
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
84 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
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,
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
101
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
130 template<typename InputType>
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
157 template<typename InputType>
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
184 {
185 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
186 return m_hCoeffs;
187 }
188
220 inline const MatrixType& packedMatrix() const
221 {
222 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
223 return m_matrix;
224 }
225
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
266 MatrixTReturnType matrixT() const
267 {
268 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
269 return MatrixTReturnType(m_matrix.real());
270 }
271
285 DiagonalReturnType diagonal() const;
286
297 SubDiagonalReturnType subDiagonal() const;
298
299 protected:
300
301 MatrixType m_matrix;
302 CoeffVectorType m_hCoeffs;
303 bool m_isInitialized;
304};
305
306template<typename MatrixType>
307typename Tridiagonalization<MatrixType>::DiagonalReturnType
309{
310 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
311 return m_matrix.diagonal().real();
312}
313
314template<typename MatrixType>
315typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
317{
318 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
319 return m_matrix.template diagonal<-1>().real();
320}
321
322namespace internal {
323
347template<typename MatrixType, typename CoeffVectorType>
348EIGEN_DEVICE_FUNC
349void 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
383template<typename MatrixType,
384 int Size=MatrixType::ColsAtCompileTime,
385 bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
386struct tridiagonalization_inplace_selector;
387
428template<typename MatrixType, typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType>
429EIGEN_DEVICE_FUNC
430void 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
440template<typename MatrixType, int Size, bool IsComplex>
441struct 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
463template<typename MatrixType>
464struct 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
509template<typename MatrixType, bool IsComplex>
510struct 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
531template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
532: public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
533{
534 public:
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
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:65
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:121
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:65
HouseholderSequenceType matrixQ() const
Returns the unitary matrix Q in the decomposition.
Definition: Tridiagonalization.h:241
Tridiagonalization(const EigenBase< InputType > &matrix)
Constructor; computes tridiagonal decomposition of given matrix.
Definition: Tridiagonalization.h:131
DiagonalReturnType diagonal() const
Returns the diagonal of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:308
MatrixTReturnType matrixT() const
Returns an expression of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:266
Eigen::Index Index
Definition: Tridiagonalization.h:73
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
Definition: Tridiagonalization.h:220
Tridiagonalization & compute(const EigenBase< InputType > &matrix)
Computes tridiagonal decomposition of given matrix.
Definition: Tridiagonalization.h:158
Tridiagonalization(Index size=Size==Dynamic ? 2 :Size)
Default constructor.
Definition: Tridiagonalization.h:114
SubDiagonalReturnType subDiagonal() const
Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:316
CoeffVectorType householderCoefficients() const
Returns the Householder coefficients.
Definition: Tridiagonalization.h:183
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: Tridiagonalization.h:69
HouseholderSequence< MatrixType, typename internal::remove_all< typename CoeffVectorType::ConjugateReturnType >::type > HouseholderSequenceType
Return type of matrixQ()
Definition: Tridiagonalization.h:100
Namespace containing all symbols from the Eigen library.
Definition: Core:141
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
Definition: Constants.h:22
Definition: EigenBase.h:30
Derived & derived()
Definition: EigenBase.h:46
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60