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 // 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>
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;
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 
85 
90 
91  typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
94  const Diagonal<
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 
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;
301 };
302 
303 template<typename MatrixType>
306 {
307  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
308  return m_matrix.diagonal();
309 }
310 
311 template<typename MatrixType>
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  using numext::conj;
349  typedef typename MatrixType::Index Index;
350  typedef typename MatrixType::Scalar Scalar;
351  typedef typename MatrixType::RealScalar RealScalar;
352  Index n = matA.rows();
353  eigen_assert(n==matA.cols());
354  eigen_assert(n==hCoeffs.size()+1 || n==1);
355 
356  for (Index i = 0; i<n-1; ++i)
357  {
358  Index remainingSize = n-i-1;
359  RealScalar beta;
360  Scalar h;
361  matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
362 
363  // Apply similarity transformation to remaining columns,
364  // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
365  matA.col(i).coeffRef(i+1) = 1;
366 
367  hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
368  * (conj(h) * matA.col(i).tail(remainingSize)));
369 
370  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);
371 
372  matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
373  .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
374 
375  matA.col(i).coeffRef(i+1) = beta;
376  hCoeffs.coeffRef(i) = h;
377  }
378 }
379 
380 // forward declaration, implementation at the end of this file
381 template<typename MatrixType,
382  int Size=MatrixType::ColsAtCompileTime,
385 
426 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
427 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
428 {
429  eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
430  tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
431 }
432 
436 template<typename MatrixType, int Size, bool IsComplex>
438 {
441  typedef typename MatrixType::Index Index;
442  template<typename DiagonalType, typename SubDiagonalType>
443  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
444  {
445  CoeffVectorType hCoeffs(mat.cols()-1);
446  tridiagonalization_inplace(mat,hCoeffs);
447  diag = mat.diagonal().real();
448  subdiag = mat.template diagonal<-1>().real();
449  if(extractQ)
450  mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
451  .setLength(mat.rows() - 1)
452  .setShift(1);
453  }
454 };
455 
460 template<typename MatrixType>
461 struct tridiagonalization_inplace_selector<MatrixType,3,false>
462 {
463  typedef typename MatrixType::Scalar Scalar;
464  typedef typename MatrixType::RealScalar RealScalar;
465 
466  template<typename DiagonalType, typename SubDiagonalType>
467  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
468  {
469  using std::sqrt;
470  diag[0] = mat(0,0);
471  RealScalar v1norm2 = numext::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(numext::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) = numext::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
HouseholderSequence< MatrixType, typename internal::remove_all< typename CoeffVectorType::ConjugateReturnType >::type > HouseholderSequenceType
Return type of matrixQ()
Matrix< Scalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 > CoeffVectorType
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
Tridiagonalization< MatrixType >::HouseholderSequenceType HouseholderSequenceType
IntermediateState sqrt(const Expression &arg)
internal::remove_all< typename MatrixType::RealReturnType >::type MatrixTypeRealView
static void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &subdiag, bool extractQ)
HouseholderSequenceType matrixQ() const
Returns the unitary matrix Q in the decomposition.
static void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &, bool extractQ)
internal::conditional< NumTraits< Scalar >::IsComplex, const CwiseUnaryOp< internal::scalar_real_op< Scalar >, const Derived >, const Derived & >::type RealReturnType
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
MatrixTReturnType matrixT() const
Returns an expression of the tridiagonal matrix T in the decomposition.
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
Tridiagonal decomposition of a selfadjoint matrix.
internal::conditional< NumTraits< Scalar >::IsComplex, typename internal::add_const_on_value_type< typename Diagonal< Block< const MatrixType, SizeMinusOne, SizeMinusOne > >::RealReturnType >::type, const Diagonal< Block< const MatrixType, SizeMinusOne, SizeMinusOne > > >::type SubDiagonalReturnType
Tridiagonalization(const MatrixType &matrix)
Constructor; computes tridiagonal decomposition of given matrix.
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs2_op< Scalar >, const Derived > abs2() const
Sequence of Householder reflections acting on subspaces with decreasing size.
RealReturnType real() const
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
Tridiagonalization< MatrixType >::CoeffVectorType CoeffVectorType
Tridiagonalization & compute(const MatrixType &matrix)
Computes tridiagonal decomposition of given matrix.
DiagonalReturnType diagonal() const
Returns the diagonal of the tridiagonal matrix T in the decomposition.
SubDiagonalReturnType subDiagonal() const
Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
Provides a generic way to set and pass user-specified options.
Definition: options.hpp:65
void tridiagonalization_inplace(MatrixType &matA, CoeffVectorType &hCoeffs)
internal::TridiagonalizationMatrixTReturnType< MatrixTypeRealView > MatrixTReturnType
Expression of a fixed-size or dynamic-size block.
Definition: Core/Block.h:102
internal::plain_col_type< MatrixType, RealScalar >::type DiagonalType
TridiagonalizationMatrixTReturnType(const MatrixType &mat)
Constructor.
CoeffVectorType householderCoefficients() const
Returns the Householder coefficients.
internal::conditional< NumTraits< Scalar >::IsComplex, typename internal::add_const_on_value_type< typename Diagonal< const MatrixType >::RealReturnType >::type, const Diagonal< const MatrixType > >::type DiagonalReturnType
Tridiagonalization(Index size=Size==Dynamic?2:Size)
Default constructor.
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:64
Matrix< RealScalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 > SubDiagonalType
#define eigen_assert(x)
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
const AutoDiffScalar< DerType > & real(const AutoDiffScalar< DerType > &x)
NumTraits< Scalar >::Real RealScalar
static void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &subdiag, bool extractQ)


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:15