Path Tracer
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  : 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 EIGEN_DEVICE_FUNC
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;
72  typedef typename NumTraits<Scalar>::Real RealScalar;
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 
83  typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
84  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
85  typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
86  typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
87  typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
88 
89  typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
90  typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
91  const Diagonal<const MatrixType>
92  >::type DiagonalReturnType;
93 
94  typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
95  typename internal::add_const_on_value_type<typename Diagonal<const MatrixType, -1>::RealReturnType>::type,
96  const Diagonal<const MatrixType, -1>
97  >::type SubDiagonalReturnType;
98 
100  typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType;
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>
131  explicit Tridiagonalization(const EigenBase<InputType>& matrix)
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>
158  Tridiagonalization& compute(const EigenBase<InputType>& matrix)
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;
303  bool m_isInitialized;
304 };
305 
306 template<typename MatrixType>
307 typename Tridiagonalization<MatrixType>::DiagonalReturnType
308 Tridiagonalization<MatrixType>::diagonal() const
309 {
310  eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
311  return m_matrix.diagonal().real();
312 }
313 
314 template<typename MatrixType>
315 typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
316 Tridiagonalization<MatrixType>::subDiagonal() const
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>
348 EIGEN_DEVICE_FUNC
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,
385  bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
386 struct tridiagonalization_inplace_selector;
387 
428 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
429 EIGEN_DEVICE_FUNC
430 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
431 {
432  eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
433  tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
434 }
435 
439 template<typename MatrixType, int Size, bool IsComplex>
440 struct tridiagonalization_inplace_selector
441 {
442  typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
443  typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
444  template<typename DiagonalType, typename SubDiagonalType>
445  static EIGEN_DEVICE_FUNC
446  void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
447  {
448  CoeffVectorType hCoeffs(mat.cols()-1);
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;
467  typedef typename MatrixType::RealScalar RealScalar;
468 
469  template<typename DiagonalType, typename SubDiagonalType>
470  static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, 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>
515  static EIGEN_DEVICE_FUNC
516  void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, 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  Index rows() const { return m_matrix.rows(); }
551  Index cols() const { 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
Eigen::Tridiagonalization::packedMatrix
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
Definition: Tridiagonalization.h:220
Eigen
Namespace containing all symbols from the Eigen library.
Definition: LDLT.h:16
Eigen::ReturnByValue
Definition: ReturnByValue.h:52
Eigen::EigenBase::derived
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
Eigen::HouseholderSequence::setLength
EIGEN_DEVICE_FUNC HouseholderSequence & setLength(Index length)
Sets the length of the Householder sequence.
Definition: HouseholderSequence.h:443
Eigen::Tridiagonalization
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:65
Eigen::EigenBase::rows
EIGEN_DEVICE_FUNC Index rows() const
Definition: EigenBase.h:60
Eigen::internal::TridiagonalizationMatrixTReturnType
Definition: Tridiagonalization.h:533
Eigen::Tridiagonalization::subDiagonal
SubDiagonalReturnType subDiagonal() const
Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:316
Eigen::EigenBase
Definition: EigenBase.h:30
Eigen::internal::TridiagonalizationMatrixTReturnType::TridiagonalizationMatrixTReturnType
TridiagonalizationMatrixTReturnType(const MatrixType &mat)
Constructor.
Definition: Tridiagonalization.h:539
Eigen::Tridiagonalization::HouseholderSequenceType
HouseholderSequence< MatrixType, typename internal::remove_all< typename CoeffVectorType::ConjugateReturnType >::type > HouseholderSequenceType
Return type of matrixQ()
Definition: Tridiagonalization.h:100
Eigen::Tridiagonalization::MatrixType
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: Tridiagonalization.h:69
Eigen::PlainObjectBase::resize
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:279
Eigen::Tridiagonalization::matrixQ
HouseholderSequenceType matrixQ() const
Returns the unitary matrix Q in the decomposition.
Definition: Tridiagonalization.h:241
Eigen::Tridiagonalization::matrixT
MatrixTReturnType matrixT() const
Returns an expression of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:266
Eigen::Tridiagonalization::compute
Tridiagonalization & compute(const EigenBase< InputType > &matrix)
Computes tridiagonal decomposition of given matrix.
Definition: Tridiagonalization.h:158
Eigen::Tridiagonalization::Index
Eigen::Index Index
Definition: Tridiagonalization.h:73
Eigen::internal::true_type
Definition: Meta.h:63
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:21
Eigen::Diagonal
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:65
Eigen::Tridiagonalization::Tridiagonalization
Tridiagonalization(const EigenBase< InputType > &matrix)
Constructor; computes tridiagonal decomposition of given matrix.
Definition: Tridiagonalization.h:131
Eigen::Tridiagonalization::diagonal
DiagonalReturnType diagonal() const
Returns the diagonal of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:308
Eigen::internal::traits
Definition: ForwardDeclarations.h:17
Eigen::internal::tridiagonalization_inplace_selector
Definition: Tridiagonalization.h:441
Eigen::internal::conditional
Definition: Meta.h:76
Eigen::Tridiagonalization::Tridiagonalization
Tridiagonalization(Index size=Size==Dynamic ? 2 :Size)
Default constructor.
Definition: Tridiagonalization.h:114
Eigen::Matrix< Scalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 >
Eigen::HouseholderSequence::setShift
EIGEN_DEVICE_FUNC HouseholderSequence & setShift(Index shift)
Sets the shift of the Householder sequence.
Definition: HouseholderSequence.h:461
Eigen::Tridiagonalization::householderCoefficients
CoeffVectorType householderCoefficients() const
Returns the Householder coefficients.
Definition: Tridiagonalization.h:183
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:213
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:42
Eigen::HouseholderSequence
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:121
Eigen::internal::add_const_on_value_type
Definition: Meta.h:175
Generated by   doxygen 1.8.19