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>
29 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
30 }
31 
64 template<typename _MatrixType> class Tridiagonalization
65 {
66  public:
67 
69  typedef _MatrixType MatrixType;
70 
71  typedef typename MatrixType::Scalar Scalar;
73  typedef Eigen::Index Index;
74 
75  enum {
76  Size = MatrixType::RowsAtCompileTime,
77  SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
78  Options = MatrixType::Options,
79  MaxSize = MatrixType::MaxRowsAtCompileTime,
80  MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
81  };
82 
88 
93 
94  typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
96  const Diagonal<const MatrixType, -1>
98 
101 
114  explicit Tridiagonalization(Index size = Size==Dynamic ? 2 : Size)
115  : m_matrix(size,size),
116  m_hCoeffs(size > 1 ? size-1 : 1),
117  m_isInitialized(false)
118  {}
119 
130  template<typename InputType>
132  : m_matrix(matrix.derived()),
133  m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
134  m_isInitialized(false)
135  {
136  internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
137  m_isInitialized = true;
138  }
139 
157  template<typename InputType>
159  {
160  m_matrix = matrix.derived();
161  m_hCoeffs.resize(matrix.rows()-1, 1);
162  internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
163  m_isInitialized = true;
164  return *this;
165  }
166 
183  inline CoeffVectorType householderCoefficients() const
184  {
185  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
186  return m_hCoeffs;
187  }
188 
220  inline const MatrixType& packedMatrix() const
221  {
222  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
223  return m_matrix;
224  }
225 
241  HouseholderSequenceType matrixQ() const
242  {
243  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
244  return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
245  .setLength(m_matrix.rows() - 1)
246  .setShift(1);
247  }
248 
266  MatrixTReturnType matrixT() const
267  {
268  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
269  return MatrixTReturnType(m_matrix.real());
270  }
271 
285  DiagonalReturnType diagonal() const;
286 
297  SubDiagonalReturnType subDiagonal() const;
298 
299  protected:
300 
301  MatrixType m_matrix;
302  CoeffVectorType m_hCoeffs;
304 };
305 
306 template<typename MatrixType>
309 {
310  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
311  return m_matrix.diagonal().real();
312 }
313 
314 template<typename MatrixType>
317 {
318  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
319  return m_matrix.template diagonal<-1>().real();
320 }
321 
322 namespace internal {
323 
347 template<typename MatrixType, typename CoeffVectorType>
349 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
350 {
351  using numext::conj;
352  typedef typename MatrixType::Scalar Scalar;
353  typedef typename MatrixType::RealScalar RealScalar;
354  Index n = matA.rows();
355  eigen_assert(n==matA.cols());
356  eigen_assert(n==hCoeffs.size()+1 || n==1);
357 
358  for (Index i = 0; i<n-1; ++i)
359  {
360  Index remainingSize = n-i-1;
361  RealScalar beta;
362  Scalar h;
363  matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
364 
365  // Apply similarity transformation to remaining columns,
366  // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
367  matA.col(i).coeffRef(i+1) = 1;
368 
369  hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
370  * (conj(h) * matA.col(i).tail(remainingSize)));
371 
372  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);
373 
374  matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
375  .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), Scalar(-1));
376 
377  matA.col(i).coeffRef(i+1) = beta;
378  hCoeffs.coeffRef(i) = h;
379  }
380 }
381 
382 // forward declaration, implementation at the end of this file
383 template<typename MatrixType,
384  int Size=MatrixType::ColsAtCompileTime,
387 
428 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType>
430 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag,
431  CoeffVectorType& hcoeffs, bool extractQ)
432 {
433  eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
434  tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, hcoeffs, extractQ);
435 }
436 
440 template<typename MatrixType, int Size, bool IsComplex>
442 {
445  template<typename DiagonalType, typename SubDiagonalType>
446  static EIGEN_DEVICE_FUNC
447  void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType& hCoeffs, bool extractQ)
448  {
449  tridiagonalization_inplace(mat, hCoeffs);
450  diag = mat.diagonal().real();
451  subdiag = mat.template diagonal<-1>().real();
452  if(extractQ)
453  mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
454  .setLength(mat.rows() - 1)
455  .setShift(1);
456  }
457 };
458 
463 template<typename MatrixType>
464 struct tridiagonalization_inplace_selector<MatrixType,3,false>
465 {
466  typedef typename MatrixType::Scalar Scalar;
468 
469  template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType>
470  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, CoeffVectorType&, bool extractQ)
471  {
472  using std::sqrt;
473  const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
474  diag[0] = mat(0,0);
475  RealScalar v1norm2 = numext::abs2(mat(2,0));
476  if(v1norm2 <= tol)
477  {
478  diag[1] = mat(1,1);
479  diag[2] = mat(2,2);
480  subdiag[0] = mat(1,0);
481  subdiag[1] = mat(2,1);
482  if (extractQ)
483  mat.setIdentity();
484  }
485  else
486  {
487  RealScalar beta = sqrt(numext::abs2(mat(1,0)) + v1norm2);
488  RealScalar invBeta = RealScalar(1)/beta;
489  Scalar m01 = mat(1,0) * invBeta;
490  Scalar m02 = mat(2,0) * invBeta;
491  Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
492  diag[1] = mat(1,1) + m02*q;
493  diag[2] = mat(2,2) - m02*q;
494  subdiag[0] = beta;
495  subdiag[1] = mat(2,1) - m01 * q;
496  if (extractQ)
497  {
498  mat << 1, 0, 0,
499  0, m01, m02,
500  0, m02, -m01;
501  }
502  }
503  }
504 };
505 
509 template<typename MatrixType, bool IsComplex>
510 struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
511 {
512  typedef typename MatrixType::Scalar Scalar;
513 
514  template<typename DiagonalType, typename SubDiagonalType, typename CoeffVectorType>
515  static EIGEN_DEVICE_FUNC
516  void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, CoeffVectorType&, bool extractQ)
517  {
518  diag(0,0) = numext::real(mat(0,0));
519  if(extractQ)
520  mat(0,0) = Scalar(1);
521  }
522 };
523 
531 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
532 : public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
533 {
534  public:
539  TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { }
540 
541  template <typename ResultType>
542  inline void evalTo(ResultType& result) const
543  {
544  result.setZero();
545  result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
546  result.diagonal() = m_matrix.diagonal();
547  result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
548  }
549 
550  EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_matrix.rows(); }
551  EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_matrix.cols(); }
552 
553  protected:
554  typename MatrixType::Nested m_matrix;
555 };
556 
557 } // end namespace internal
558 
559 } // end namespace Eigen
560 
561 #endif // EIGEN_TRIDIAGONALIZATION_H
Tridiagonalization(Index size=Size==Dynamic ? 2 :Size)
Default constructor.
SCALAR Scalar
Definition: bench_gemm.cpp:46
DiagonalReturnType diagonal() const
Returns the diagonal of the tridiagonal matrix T in the decomposition.
HouseholderSequence< MatrixType, typename internal::remove_all< typename CoeffVectorType::ConjugateReturnType >::type > HouseholderSequenceType
Return type of matrixQ()
Matrix< Scalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 > CoeffVectorType
static EIGEN_DEVICE_FUNC void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &subdiag, CoeffVectorType &hCoeffs, bool extractQ)
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
static EIGEN_DEVICE_FUNC void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &, CoeffVectorType &, bool extractQ)
Matrix diag(const std::vector< Matrix > &Hs)
Definition: Matrix.cpp:206
#define min(a, b)
Definition: datatypes.h:19
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
int n
void diagonal(const MatrixType &m)
Definition: diagonal.cpp:12
HouseholderSequenceType matrixQ() const
Returns the unitary matrix Q in the decomposition.
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.
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
Tridiagonal decomposition of a selfadjoint matrix.
AnnoyingScalar conj(const AnnoyingScalar &x)
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Sequence of Householder reflections acting on subspaces with decreasing size.
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
CoeffVectorType householderCoefficients() const
Returns the Householder coefficients.
Tridiagonalization< MatrixType >::CoeffVectorType CoeffVectorType
Values result
#define EIGEN_NOEXCEPT
Definition: Macros.h:1418
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
#define eigen_assert(x)
Definition: Macros.h:1037
MatrixTReturnType matrixT() const
Returns an expression of the tridiagonal matrix T in the decomposition.
#define EIGEN_CONSTEXPR
Definition: Macros.h:787
EIGEN_DEVICE_FUNC const Scalar & q
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
MatrixXf matA(2, 2)
Tridiagonalization(const EigenBase< InputType > &matrix)
Constructor; computes tridiagonal decomposition of given matrix.
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
#define EIGEN_DEVICE_FUNC
Definition: Macros.h:976
const double h
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
internal::TridiagonalizationMatrixTReturnType< MatrixTypeRealView > MatrixTReturnType
internal::plain_col_type< MatrixType, RealScalar >::type DiagonalType
TridiagonalizationMatrixTReturnType(const MatrixType &mat)
Constructor.
internal::conditional< NumTraits< Scalar >::IsComplex, typename internal::add_const_on_value_type< typename Diagonal< const MatrixType >::RealReturnType >::type, const Diagonal< const MatrixType > >::type DiagonalReturnType
EIGEN_DEVICE_FUNC void tridiagonalization_inplace(MatrixType &matA, CoeffVectorType &hCoeffs)
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:63
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
const G double tol
Definition: Group.h:86
const int Dynamic
Definition: Constants.h:22
Matrix< RealScalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 > SubDiagonalType
Generic expression where a coefficient-wise unary operator is applied to an expression.
Definition: CwiseUnaryOp.h:55
static void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &subdiag, CoeffVectorType &, bool extractQ)
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 bool abs2(bool x)
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
SubDiagonalReturnType subDiagonal() const
Returns the subdiagonal of the tridiagonal matrix T in the decomposition.


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:40:38