11 #ifndef EIGEN_TRIDIAGONALIZATION_H 12 #define EIGEN_TRIDIAGONALIZATION_H 18 template<
typename MatrixType>
struct TridiagonalizationMatrixTReturnType;
19 template<
typename MatrixType>
20 struct traits<TridiagonalizationMatrixTReturnType<MatrixType> >
21 :
public traits<typename MatrixType::PlainObject>
23 typedef typename MatrixType::PlainObject ReturnType;
27 template<
typename MatrixType,
typename CoeffVectorType>
29 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
71 typedef typename MatrixType::Scalar Scalar;
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)
84 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
86 typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
87 typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
89 typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
90 typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
92 >::type DiagonalReturnType;
94 typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
95 typename internal::add_const_on_value_type<
typename Diagonal<
const MatrixType, -1>::RealReturnType>::type,
97 >::type SubDiagonalReturnType;
115 : m_matrix(size,size),
116 m_hCoeffs(size > 1 ? size-1 : 1),
117 m_isInitialized(false)
130 template<
typename InputType>
132 : m_matrix(matrix.derived()),
133 m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
134 m_isInitialized(false)
136 internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
137 m_isInitialized =
true;
157 template<
typename InputType>
161 m_hCoeffs.resize(matrix.
rows()-1, 1);
162 internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
163 m_isInitialized =
true;
185 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
222 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
243 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
244 return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
245 .setLength(m_matrix.rows() - 1)
268 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
269 return MatrixTReturnType(m_matrix.real());
285 DiagonalReturnType diagonal()
const;
297 SubDiagonalReturnType subDiagonal()
const;
302 CoeffVectorType m_hCoeffs;
303 bool m_isInitialized;
306 template<
typename MatrixType>
307 typename Tridiagonalization<MatrixType>::DiagonalReturnType
310 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
311 return m_matrix.diagonal().real();
314 template<
typename MatrixType>
315 typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
318 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
319 return m_matrix.template diagonal<-1>().
real();
347 template<
typename MatrixType,
typename CoeffVectorType>
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);
358 for (
Index i = 0; i<n-1; ++i)
360 Index remainingSize = n-i-1;
363 matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
367 matA.col(i).coeffRef(i+1) = 1;
369 hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
370 * (
conj(h) * matA.col(i).tail(remainingSize)));
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);
374 matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
375 .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), Scalar(-1));
384 int Size=MatrixType::ColsAtCompileTime,
386 struct tridiagonalization_inplace_selector;
428 template<
typename MatrixType,
typename DiagonalType,
typename SubDiagonalType,
typename CoeffVectorType>
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);
440 template<
typename MatrixType,
int Size,
bool IsComplex>
441 struct tridiagonalization_inplace_selector
445 template<
typename DiagonalType,
typename SubDiagonalType>
446 static EIGEN_DEVICE_FUNC
449 tridiagonalization_inplace(mat, hCoeffs);
450 diag = mat.diagonal().real();
451 subdiag = mat.template diagonal<-1>().
real();
453 mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
454 .setLength(mat.rows() - 1)
463 template<
typename MatrixType>
464 struct tridiagonalization_inplace_selector<MatrixType,3,false>
466 typedef typename MatrixType::Scalar Scalar;
467 typedef typename MatrixType::RealScalar RealScalar;
469 template<
typename DiagonalType,
typename SubDiagonalType,
typename CoeffVectorType>
473 const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
475 RealScalar v1norm2 = numext::abs2(mat(2,0));
480 subdiag[0] = mat(1,0);
481 subdiag[1] = mat(2,1);
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;
495 subdiag[1] = mat(2,1) - m01 * q;
509 template<
typename MatrixType,
bool IsComplex>
510 struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
512 typedef typename MatrixType::Scalar Scalar;
514 template<
typename DiagonalType,
typename SubDiagonalType,
typename CoeffVectorType>
515 static EIGEN_DEVICE_FUNC
518 diag(0,0) = numext::real(mat(0,0));
520 mat(0,0) = Scalar(1);
531 template<
typename MatrixType>
struct TridiagonalizationMatrixTReturnType
532 :
public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
539 TridiagonalizationMatrixTReturnType(
const MatrixType& mat) : m_matrix(mat) { }
541 template <
typename ResultType>
542 inline void evalTo(ResultType& result)
const 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>();
550 EIGEN_CONSTEXPR
Index rows()
const EIGEN_NOEXCEPT {
return m_matrix.rows(); }
551 EIGEN_CONSTEXPR
Index cols()
const EIGEN_NOEXCEPT {
return m_matrix.cols(); }
554 typename MatrixType::Nested m_matrix;
561 #endif // EIGEN_TRIDIAGONALIZATION_H Tridiagonalization(Index size=Size==Dynamic ? 2 :Size)
Default constructor.
Definition: Tridiagonalization.h:114
DiagonalReturnType diagonal() const
Returns the diagonal of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:308
HouseholderSequence< MatrixType, typename internal::remove_all< typename CoeffVectorType::ConjugateReturnType >::type > HouseholderSequenceType
Return type of matrixQ()
Definition: Tridiagonalization.h:100
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Eigen::Index Index
Definition: Tridiagonalization.h:73
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
HouseholderSequenceType matrixQ() const
Returns the unitary matrix Q in the decomposition.
Definition: Tridiagonalization.h:241
Namespace containing all symbols from the Eigen library.
Definition: Core:141
Tridiagonalization & compute(const EigenBase< InputType > &matrix)
Computes tridiagonal decomposition of given matrix.
Definition: Tridiagonalization.h:158
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:64
Derived & derived()
Definition: EigenBase.h:46
EIGEN_CONSTEXPR Index size() const EIGEN_NOEXCEPT
Definition: EigenBase.h:67
Definition: EigenBase.h:29
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:119
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:175
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
Definition: Tridiagonalization.h:220
CoeffVectorType householderCoefficients() const
Returns the Householder coefficients.
Definition: Tridiagonalization.h:183
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
MatrixTReturnType matrixT() const
Returns an expression of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:266
Tridiagonalization(const EigenBase< InputType > &matrix)
Constructor; computes tridiagonal decomposition of given matrix.
Definition: Tridiagonalization.h:131
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
Definition: Eigen_Colamd.h:50
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:63
const int Dynamic
Definition: Constants.h:22
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: Tridiagonalization.h:69
SubDiagonalReturnType subDiagonal() const
Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
Definition: Tridiagonalization.h:316