Eigen  3.2.91
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 
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 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
29 }
30 
63 template<typename _MatrixType> class Tridiagonalization
64 {
65  public:
66 
68  typedef _MatrixType MatrixType;
69 
70  typedef typename MatrixType::Scalar Scalar;
71  typedef typename NumTraits<Scalar>::Real RealScalar;
72  typedef Eigen::Index Index;
73 
74  enum {
75  Size = MatrixType::RowsAtCompileTime,
76  SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
77  Options = MatrixType::Options,
78  MaxSize = MatrixType::MaxRowsAtCompileTime,
79  MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
80  };
81 
83  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
85  typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
86  typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
87 
88  typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
89  typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
91  >::type DiagonalReturnType;
92 
93  typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
94  typename internal::add_const_on_value_type<typename Diagonal<const MatrixType, -1>::RealReturnType>::type,
95  const Diagonal<const MatrixType, -1>
96  >::type SubDiagonalReturnType;
97 
100 
113  explicit 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  explicit 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 
303 template<typename MatrixType>
304 typename Tridiagonalization<MatrixType>::DiagonalReturnType
306 {
307  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
308  return m_matrix.diagonal().real();
309 }
310 
311 template<typename MatrixType>
312 typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
314 {
315  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
316  return m_matrix.template diagonal<-1>().real();
317 }
318 
319 namespace internal {
320 
344 template<typename MatrixType, typename CoeffVectorType>
345 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
346 {
347  using numext::conj;
348  typedef typename MatrixType::Scalar Scalar;
349  typedef typename MatrixType::RealScalar RealScalar;
350  Index n = matA.rows();
351  eigen_assert(n==matA.cols());
352  eigen_assert(n==hCoeffs.size()+1 || n==1);
353 
354  for (Index i = 0; i<n-1; ++i)
355  {
356  Index remainingSize = n-i-1;
357  RealScalar beta;
358  Scalar h;
359  matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
360 
361  // Apply similarity transformation to remaining columns,
362  // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
363  matA.col(i).coeffRef(i+1) = 1;
364 
365  hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
366  * (conj(h) * matA.col(i).tail(remainingSize)));
367 
368  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);
369 
370  matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
371  .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
372 
373  matA.col(i).coeffRef(i+1) = beta;
374  hCoeffs.coeffRef(i) = h;
375  }
376 }
377 
378 // forward declaration, implementation at the end of this file
379 template<typename MatrixType,
380  int Size=MatrixType::ColsAtCompileTime,
381  bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
382 struct tridiagonalization_inplace_selector;
383 
424 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
425 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
426 {
427  eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
428  tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
429 }
430 
434 template<typename MatrixType, int Size, bool IsComplex>
435 struct tridiagonalization_inplace_selector
436 {
437  typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
438  typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
439  template<typename DiagonalType, typename SubDiagonalType>
440  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
441  {
442  CoeffVectorType hCoeffs(mat.cols()-1);
443  tridiagonalization_inplace(mat,hCoeffs);
444  diag = mat.diagonal().real();
445  subdiag = mat.template diagonal<-1>().real();
446  if(extractQ)
447  mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
448  .setLength(mat.rows() - 1)
449  .setShift(1);
450  }
451 };
452 
457 template<typename MatrixType>
458 struct tridiagonalization_inplace_selector<MatrixType,3,false>
459 {
460  typedef typename MatrixType::Scalar Scalar;
461  typedef typename MatrixType::RealScalar RealScalar;
462 
463  template<typename DiagonalType, typename SubDiagonalType>
464  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
465  {
466  using std::sqrt;
467  const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
468  diag[0] = mat(0,0);
469  RealScalar v1norm2 = numext::abs2(mat(2,0));
470  if(v1norm2 <= tol)
471  {
472  diag[1] = mat(1,1);
473  diag[2] = mat(2,2);
474  subdiag[0] = mat(1,0);
475  subdiag[1] = mat(2,1);
476  if (extractQ)
477  mat.setIdentity();
478  }
479  else
480  {
481  RealScalar beta = sqrt(numext::abs2(mat(1,0)) + v1norm2);
482  RealScalar invBeta = RealScalar(1)/beta;
483  Scalar m01 = mat(1,0) * invBeta;
484  Scalar m02 = mat(2,0) * invBeta;
485  Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
486  diag[1] = mat(1,1) + m02*q;
487  diag[2] = mat(2,2) - m02*q;
488  subdiag[0] = beta;
489  subdiag[1] = mat(2,1) - m01 * q;
490  if (extractQ)
491  {
492  mat << 1, 0, 0,
493  0, m01, m02,
494  0, m02, -m01;
495  }
496  }
497  }
498 };
499 
503 template<typename MatrixType, bool IsComplex>
504 struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
505 {
506  typedef typename MatrixType::Scalar Scalar;
507 
508  template<typename DiagonalType, typename SubDiagonalType>
509  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ)
510  {
511  diag(0,0) = numext::real(mat(0,0));
512  if(extractQ)
513  mat(0,0) = Scalar(1);
514  }
515 };
516 
524 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
525 : public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
526 {
527  public:
532  TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { }
533 
534  template <typename ResultType>
535  inline void evalTo(ResultType& result) const
536  {
537  result.setZero();
538  result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
539  result.diagonal() = m_matrix.diagonal();
540  result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
541  }
542 
543  Index rows() const { return m_matrix.rows(); }
544  Index cols() const { return m_matrix.cols(); }
545 
546  protected:
547  typename MatrixType::Nested m_matrix;
548 };
549 
550 } // end namespace internal
551 
552 } // end namespace Eigen
553 
554 #endif // EIGEN_TRIDIAGONALIZATION_H
HouseholderSequence< MatrixType, typename internal::remove_all< typename CoeffVectorType::ConjugateReturnType >::type > HouseholderSequenceType
Return type of matrixQ()
Definition: Tridiagonalization.h:99
HouseholderSequenceType matrixQ() const
Returns the unitary matrix Q in the decomposition.
Definition: Tridiagonalization.h:238
Eigen::Index Index
Definition: Tridiagonalization.h:72
Definition: LDLT.h:16
MatrixTReturnType matrixT() const
Returns an expression of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:263
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:63
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:252
Tridiagonalization(const MatrixType &matrix)
Constructor; computes tridiagonal decomposition of given matrix.
Definition: Tridiagonalization.h:129
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: ForwardDeclarations.h:258
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
Definition: Tridiagonalization.h:217
Tridiagonalization & compute(const MatrixType &matrix)
Computes tridiagonal decomposition of given matrix.
Definition: Tridiagonalization.h:155
DiagonalReturnType diagonal() const
Returns the diagonal of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:305
SubDiagonalReturnType subDiagonal() const
Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:313
Definition: Eigen_Colamd.h:54
CoeffVectorType householderCoefficients() const
Returns the Householder coefficients.
Definition: Tridiagonalization.h:180
Tridiagonalization(Index size=Size==Dynamic?2:Size)
Default constructor.
Definition: Tridiagonalization.h:113
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:63
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: Tridiagonalization.h:68