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 {
22  typedef typename MatrixType::PlainObject ReturnType;
23 };
24 
25 template<typename MatrixType, typename CoeffVectorType>
26 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
27 }
28 
61 template<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 
81  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
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,
89  >::type DiagonalReturnType;
90 
91  typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
92  typename internal::add_const_on_value_type<typename Diagonal<
94  const Diagonal<
96  >::type SubDiagonalReturnType;
97 
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 
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 
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 
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 
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();
309 }
310 
311 template<typename MatrixType>
312 typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
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 
320 namespace internal {
321 
345 template<typename MatrixType, typename CoeffVectorType>
346 void 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
380 template<typename MatrixType,
381  int Size=MatrixType::ColsAtCompileTime,
382  bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
383 struct tridiagonalization_inplace_selector;
384 
425 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
426 void 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 
437 template<typename MatrixType, int Size, bool IsComplex>
438 struct 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 
461 template<typename MatrixType>
462 struct 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 
505 template<typename MatrixType, bool IsComplex>
506 struct 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 
526 template<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