Tridiagonalization.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 #ifndef EIGEN_TRIDIAGONALIZATION_H
27 #define EIGEN_TRIDIAGONALIZATION_H
28 
29 namespace Eigen {
30 
31 namespace internal {
32 
33 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType;
34 template<typename MatrixType>
35 struct traits<TridiagonalizationMatrixTReturnType<MatrixType> >
36 {
37  typedef typename MatrixType::PlainObject ReturnType;
38 };
39 
40 template<typename MatrixType, typename CoeffVectorType>
41 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
42 }
43 
76 template<typename _MatrixType> class Tridiagonalization
77 {
78  public:
79 
81  typedef _MatrixType MatrixType;
82 
83  typedef typename MatrixType::Scalar Scalar;
85  typedef typename MatrixType::Index Index;
86 
87  enum {
88  Size = MatrixType::RowsAtCompileTime,
89  SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
90  Options = MatrixType::Options,
91  MaxSize = MatrixType::MaxRowsAtCompileTime,
93  };
94 
96  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
98  typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
99  typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
100 
102  typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
105 
107  typename internal::add_const_on_value_type<typename Diagonal<
108  Block<const MatrixType,SizeMinusOne,SizeMinusOne> >::RealReturnType>::type,
109  const Diagonal<
112 
115 
129  : m_matrix(size,size),
130  m_hCoeffs(size > 1 ? size-1 : 1),
131  m_isInitialized(false)
132  {}
133 
145  : m_matrix(matrix),
146  m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
147  m_isInitialized(false)
148  {
150  m_isInitialized = true;
151  }
152 
171  {
172  m_matrix = matrix;
173  m_hCoeffs.resize(matrix.rows()-1, 1);
175  m_isInitialized = true;
176  return *this;
177  }
178 
196  {
197  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
198  return m_hCoeffs;
199  }
200 
232  inline const MatrixType& packedMatrix() const
233  {
234  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
235  return m_matrix;
236  }
237 
254  {
255  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
256  return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
257  .setLength(m_matrix.rows() - 1)
258  .setShift(1);
259  }
260 
279  {
280  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
281  return MatrixTReturnType(m_matrix.real());
282  }
283 
298 
310 
311  protected:
312 
316 };
317 
318 template<typename MatrixType>
321 {
322  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
323  return m_matrix.diagonal();
324 }
325 
326 template<typename MatrixType>
329 {
330  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
331  Index n = m_matrix.rows();
332  return Block<const MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal();
333 }
334 
335 namespace internal {
336 
360 template<typename MatrixType, typename CoeffVectorType>
361 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
362 {
363  typedef typename MatrixType::Index Index;
364  typedef typename MatrixType::Scalar Scalar;
365  typedef typename MatrixType::RealScalar RealScalar;
366  Index n = matA.rows();
367  eigen_assert(n==matA.cols());
368  eigen_assert(n==hCoeffs.size()+1 || n==1);
369 
370  for (Index i = 0; i<n-1; ++i)
371  {
372  Index remainingSize = n-i-1;
373  RealScalar beta;
374  Scalar h;
375  matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
376 
377  // Apply similarity transformation to remaining columns,
378  // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
379  matA.col(i).coeffRef(i+1) = 1;
380 
381  hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
382  * (conj(h) * matA.col(i).tail(remainingSize)));
383 
384  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);
385 
386  matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
387  .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
388 
389  matA.col(i).coeffRef(i+1) = beta;
390  hCoeffs.coeffRef(i) = h;
391  }
392 }
393 
394 // forward declaration, implementation at the end of this file
395 template<typename MatrixType,
396  int Size=MatrixType::ColsAtCompileTime,
398 struct tridiagonalization_inplace_selector;
399 
440 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
441 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
442 {
443  typedef typename MatrixType::Index Index;
444  //Index n = mat.rows();
445  eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
446  tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
447 }
448 
452 template<typename MatrixType, int Size, bool IsComplex>
453 struct tridiagonalization_inplace_selector
454 {
455  typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
456  typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
457  typedef typename MatrixType::Index Index;
458  template<typename DiagonalType, typename SubDiagonalType>
459  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
460  {
461  CoeffVectorType hCoeffs(mat.cols()-1);
462  tridiagonalization_inplace(mat,hCoeffs);
463  diag = mat.diagonal().real();
464  subdiag = mat.template diagonal<-1>().real();
465  if(extractQ)
466  mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
467  .setLength(mat.rows() - 1)
468  .setShift(1);
469  }
470 };
471 
476 template<typename MatrixType>
477 struct tridiagonalization_inplace_selector<MatrixType,3,false>
478 {
479  typedef typename MatrixType::Scalar Scalar;
480  typedef typename MatrixType::RealScalar RealScalar;
481 
482  template<typename DiagonalType, typename SubDiagonalType>
483  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
484  {
485  diag[0] = mat(0,0);
486  RealScalar v1norm2 = abs2(mat(2,0));
487  if(v1norm2 == RealScalar(0))
488  {
489  diag[1] = mat(1,1);
490  diag[2] = mat(2,2);
491  subdiag[0] = mat(1,0);
492  subdiag[1] = mat(2,1);
493  if (extractQ)
494  mat.setIdentity();
495  }
496  else
497  {
498  RealScalar beta = sqrt(abs2(mat(1,0)) + v1norm2);
499  RealScalar invBeta = RealScalar(1)/beta;
500  Scalar m01 = mat(1,0) * invBeta;
501  Scalar m02 = mat(2,0) * invBeta;
502  Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
503  diag[1] = mat(1,1) + m02*q;
504  diag[2] = mat(2,2) - m02*q;
505  subdiag[0] = beta;
506  subdiag[1] = mat(2,1) - m01 * q;
507  if (extractQ)
508  {
509  mat << 1, 0, 0,
510  0, m01, m02,
511  0, m02, -m01;
512  }
513  }
514  }
515 };
516 
520 template<typename MatrixType, bool IsComplex>
521 struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
522 {
523  typedef typename MatrixType::Scalar Scalar;
524 
525  template<typename DiagonalType, typename SubDiagonalType>
526  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ)
527  {
528  diag(0,0) = real(mat(0,0));
529  if(extractQ)
530  mat(0,0) = Scalar(1);
531  }
532 };
533 
541 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
542 : public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
543 {
544  typedef typename MatrixType::Index Index;
545  public:
550  TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { }
551 
552  template <typename ResultType>
553  inline void evalTo(ResultType& result) const
554  {
555  result.setZero();
556  result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
557  result.diagonal() = m_matrix.diagonal();
558  result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
559  }
560 
561  Index rows() const { return m_matrix.rows(); }
562  Index cols() const { return m_matrix.cols(); }
563 
564  protected:
565  typename MatrixType::Nested m_matrix;
566 };
567 
568 } // end namespace internal
569 
570 } // end namespace Eigen
571 
572 #endif // EIGEN_TRIDIAGONALIZATION_H