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  : 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;
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 
87 
92 
93  typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
95  const Diagonal<const MatrixType, -1>
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  template<typename InputType>
131  : m_matrix(matrix.derived()),
132  m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
133  m_isInitialized(false)
134  {
135  internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
136  m_isInitialized = true;
137  }
138 
156  template<typename InputType>
158  {
159  m_matrix = matrix.derived();
160  m_hCoeffs.resize(matrix.rows()-1, 1);
161  internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
162  m_isInitialized = true;
163  return *this;
164  }
165 
182  inline CoeffVectorType householderCoefficients() const
183  {
184  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
185  return m_hCoeffs;
186  }
187 
219  inline const MatrixType& packedMatrix() const
220  {
221  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
222  return m_matrix;
223  }
224 
240  HouseholderSequenceType matrixQ() const
241  {
242  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
243  return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
244  .setLength(m_matrix.rows() - 1)
245  .setShift(1);
246  }
247 
265  MatrixTReturnType matrixT() const
266  {
267  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
268  return MatrixTReturnType(m_matrix.real());
269  }
270 
284  DiagonalReturnType diagonal() const;
285 
296  SubDiagonalReturnType subDiagonal() const;
297 
298  protected:
299 
300  MatrixType m_matrix;
301  CoeffVectorType m_hCoeffs;
303 };
304 
305 template<typename MatrixType>
308 {
309  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
310  return m_matrix.diagonal().real();
311 }
312 
313 template<typename MatrixType>
316 {
317  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
318  return m_matrix.template diagonal<-1>().real();
319 }
320 
321 namespace internal {
322 
346 template<typename MatrixType, typename CoeffVectorType>
347 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
348 {
349  using numext::conj;
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)*RealScalar(-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), Scalar(-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  template<typename DiagonalType, typename SubDiagonalType>
442  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
443  {
444  CoeffVectorType hCoeffs(mat.cols()-1);
445  tridiagonalization_inplace(mat,hCoeffs);
446  diag = mat.diagonal().real();
447  subdiag = mat.template diagonal<-1>().real();
448  if(extractQ)
449  mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
450  .setLength(mat.rows() - 1)
451  .setShift(1);
452  }
453 };
454 
459 template<typename MatrixType>
460 struct tridiagonalization_inplace_selector<MatrixType,3,false>
461 {
462  typedef typename MatrixType::Scalar Scalar;
464 
465  template<typename DiagonalType, typename SubDiagonalType>
466  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
467  {
468  using std::sqrt;
469  const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
470  diag[0] = mat(0,0);
471  RealScalar v1norm2 = numext::abs2(mat(2,0));
472  if(v1norm2 <= tol)
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  public:
534  TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { }
535 
536  template <typename ResultType>
537  inline void evalTo(ResultType& result) const
538  {
539  result.setZero();
540  result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
541  result.diagonal() = m_matrix.diagonal();
542  result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
543  }
544 
545  Index rows() const { return m_matrix.rows(); }
546  Index cols() const { return m_matrix.cols(); }
547 
548  protected:
549  typename MatrixType::Nested m_matrix;
550 };
551 
552 } // end namespace internal
553 
554 } // end namespace Eigen
555 
556 #endif // EIGEN_TRIDIAGONALIZATION_H
SCALAR Scalar
Definition: bench_gemm.cpp:33
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
internal::remove_all< typename MatrixType::RealReturnType >::type MatrixTypeRealView
float real
Definition: datatypes.h:10
Matrix diag(const std::vector< Matrix > &Hs)
Definition: Matrix.cpp:206
#define min(a, b)
Definition: datatypes.h:19
static void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &subdiag, bool extractQ)
HouseholderSequenceType matrixQ() const
Returns the unitary matrix Q in the decomposition.
internal::conditional< NumTraits< Scalar >::IsComplex, typename internal::add_const_on_value_type< typename Diagonal< const MatrixType,-1 >::RealReturnType >::type, const Diagonal< const MatrixType,-1 > >::type SubDiagonalReturnType
static void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &, bool extractQ)
int n
void diagonal(const MatrixType &m)
Definition: diagonal.cpp:12
EIGEN_DEVICE_FUNC const SqrtReturnType sqrt() const
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
MatrixXf MatrixType
Tridiagonalization & compute(const EigenBase< InputType > &matrix)
Computes tridiagonal decomposition of given matrix.
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:150
Tridiagonal decomposition of a selfadjoint matrix.
Sequence of Householder reflections acting on subspaces with decreasing size.
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
Tridiagonalization< MatrixType >::CoeffVectorType CoeffVectorType
Values result
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
#define eigen_assert(x)
Definition: Macros.h:579
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.
EIGEN_DEVICE_FUNC const Scalar & q
MatrixXf matA(2, 2)
Tridiagonalization(const EigenBase< InputType > &matrix)
Constructor; computes tridiagonal decomposition of given matrix.
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Abs2ReturnType abs2() const
const double h
void tridiagonalization_inplace(MatrixType &matA, CoeffVectorType &hCoeffs)
internal::TridiagonalizationMatrixTReturnType< MatrixTypeRealView > MatrixTReturnType
internal::plain_col_type< MatrixType, RealScalar >::type DiagonalType
TridiagonalizationMatrixTReturnType(const MatrixType &mat)
Constructor.
CoeffVectorType householderCoefficients() const
Returns the Householder coefficients.
EIGEN_DEVICE_FUNC Index rows() const
Definition: EigenBase.h:59
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:63
const G double tol
Definition: Group.h:83
const int Dynamic
Definition: Constants.h:21
Matrix< RealScalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 > SubDiagonalType
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
const AutoDiffScalar< DerType > & real(const AutoDiffScalar< DerType > &x)
NumTraits< Scalar >::Real RealScalar
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:45
static void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &subdiag, bool extractQ)
ScalarWithExceptions conj(const ScalarWithExceptions &x)
Definition: exceptions.cpp:74


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:51:19