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> >
22 typedef typename MatrixType::PlainObject ReturnType;
25 template<
typename MatrixType,
typename CoeffVectorType>
26 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
68 typedef typename MatrixType::Scalar Scalar;
70 typedef typename MatrixType::Index Index;
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)
81 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
83 typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
84 typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
86 typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
87 typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
89 >::type DiagonalReturnType;
91 typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
92 typename internal::add_const_on_value_type<
typename Diagonal<
96 >::type SubDiagonalReturnType;
114 : m_matrix(size,size),
115 m_hCoeffs(size > 1 ? size-1 : 1),
116 m_isInitialized(false)
131 m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
132 m_isInitialized(false)
134 internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
135 m_isInitialized =
true;
158 m_hCoeffs.
resize(matrix.rows()-1, 1);
159 internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
160 m_isInitialized =
true;
182 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
219 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
240 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
242 .setLength(m_matrix.rows() - 1)
265 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
266 return MatrixTReturnType(m_matrix.real());
282 DiagonalReturnType
diagonal()
const;
299 CoeffVectorType m_hCoeffs;
300 bool m_isInitialized;
303 template<
typename MatrixType>
304 typename Tridiagonalization<MatrixType>::DiagonalReturnType
307 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
308 return m_matrix.diagonal();
311 template<
typename MatrixType>
312 typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
315 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
316 Index n = m_matrix.rows();
345 template<
typename MatrixType,
typename CoeffVectorType>
346 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
348 typedef typename MatrixType::Index Index;
349 typedef typename MatrixType::Scalar Scalar;
350 typedef typename MatrixType::RealScalar RealScalar;
351 Index n = matA.rows();
352 eigen_assert(n==matA.cols());
353 eigen_assert(n==hCoeffs.size()+1 || n==1);
355 for (Index i = 0; i<n-1; ++i)
357 Index remainingSize = n-i-1;
360 matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
364 matA.col(i).coeffRef(i+1) = 1;
366 hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
367 * (
conj(h) * matA.col(i).tail(remainingSize)));
369 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 matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
372 .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
374 matA.col(i).coeffRef(i+1) = beta;
375 hCoeffs.coeffRef(i) = h;
380 template<
typename MatrixType,
381 int Size=MatrixType::ColsAtCompileTime,
382 bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
383 struct tridiagonalization_inplace_selector;
425 template<
typename MatrixType,
typename DiagonalType,
typename SubDiagonalType>
426 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag,
bool extractQ)
428 typedef typename MatrixType::Index Index;
430 eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
431 tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
437 template<
typename MatrixType,
int Size,
bool IsComplex>
438 struct tridiagonalization_inplace_selector
442 typedef typename MatrixType::Index Index;
443 template<
typename DiagonalType,
typename SubDiagonalType>
444 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag,
bool extractQ)
446 CoeffVectorType hCoeffs(mat.cols()-1);
447 tridiagonalization_inplace(mat,hCoeffs);
448 diag = mat.diagonal().real();
449 subdiag = mat.template diagonal<-1>().
real();
451 mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
452 .setLength(mat.rows() - 1)
461 template<
typename MatrixType>
462 struct tridiagonalization_inplace_selector<MatrixType,3,false>
464 typedef typename MatrixType::Scalar Scalar;
465 typedef typename MatrixType::RealScalar RealScalar;
467 template<
typename DiagonalType,
typename SubDiagonalType>
468 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag,
bool extractQ)
471 RealScalar v1norm2 =
abs2(mat(2,0));
472 if(v1norm2 == RealScalar(0))
476 subdiag[0] = mat(1,0);
477 subdiag[1] = mat(2,1);
483 RealScalar beta = sqrt(
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;
491 subdiag[1] = mat(2,1) - m01 * q;
505 template<
typename MatrixType,
bool IsComplex>
506 struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
508 typedef typename MatrixType::Scalar Scalar;
510 template<
typename DiagonalType,
typename SubDiagonalType>
511 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&,
bool extractQ)
513 diag(0,0) =
real(mat(0,0));
515 mat(0,0) = Scalar(1);
526 template<
typename MatrixType>
struct TridiagonalizationMatrixTReturnType
527 :
public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
529 typedef typename MatrixType::Index Index;
535 TridiagonalizationMatrixTReturnType(
const MatrixType& mat) : m_matrix(mat) { }
537 template <
typename ResultType>
538 inline void evalTo(ResultType& result)
const
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>();
546 Index rows()
const {
return m_matrix.rows(); }
547 Index cols()
const {
return m_matrix.cols(); }
550 typename MatrixType::Nested m_matrix;
557 #endif // EIGEN_TRIDIAGONALIZATION_H