26 #ifndef EIGEN_TRIDIAGONALIZATION_H
27 #define EIGEN_TRIDIAGONALIZATION_H
33 template<
typename MatrixType>
struct TridiagonalizationMatrixTReturnType;
34 template<
typename MatrixType>
35 struct traits<TridiagonalizationMatrixTReturnType<MatrixType> >
37 typedef typename MatrixType::PlainObject ReturnType;
40 template<
typename MatrixType,
typename CoeffVectorType>
83 typedef typename MatrixType::Scalar
Scalar;
85 typedef typename MatrixType::Index
Index;
88 Size = MatrixType::RowsAtCompileTime,
91 MaxSize = MatrixType::MaxRowsAtCompileTime,
96 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type
DiagonalType;
98 typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type
MatrixTypeRealView;
99 typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView>
MatrixTReturnType;
102 typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
107 typename internal::add_const_on_value_type<
typename Diagonal<
146 m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
318 template<
typename MatrixType>
322 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
323 return m_matrix.diagonal();
326 template<
typename MatrixType>
330 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
331 Index n = m_matrix.rows();
360 template<
typename MatrixType,
typename CoeffVectorType>
363 typedef typename MatrixType::Index Index;
364 typedef typename MatrixType::Scalar Scalar;
365 typedef typename MatrixType::RealScalar RealScalar;
366 Index n = matA.rows();
370 for (Index i = 0; i<n-1; ++i)
372 Index remainingSize = n-i-1;
375 matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
379 matA.col(i).coeffRef(i+1) = 1;
381 hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
382 * (
conj(h) * matA.col(i).tail(remainingSize)));
384 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);
386 matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
387 .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
389 matA.col(i).coeffRef(i+1) = beta;
390 hCoeffs.coeffRef(i) = h;
395 template<
typename MatrixType,
396 int Size=MatrixType::ColsAtCompileTime,
398 struct tridiagonalization_inplace_selector;
440 template<
typename MatrixType,
typename DiagonalType,
typename SubDiagonalType>
443 typedef typename MatrixType::Index Index;
445 eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
446 tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
452 template<
typename MatrixType,
int Size,
bool IsComplex>
453 struct tridiagonalization_inplace_selector
457 typedef typename MatrixType::Index Index;
458 template<
typename DiagonalType,
typename SubDiagonalType>
459 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag,
bool extractQ)
461 CoeffVectorType hCoeffs(mat.cols()-1);
463 diag = mat.diagonal().real();
464 subdiag = mat.template diagonal<-1>().
real();
466 mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
467 .setLength(mat.rows() - 1)
476 template<
typename MatrixType>
477 struct tridiagonalization_inplace_selector<MatrixType,3,false>
479 typedef typename MatrixType::Scalar Scalar;
480 typedef typename MatrixType::RealScalar RealScalar;
482 template<
typename DiagonalType,
typename SubDiagonalType>
483 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag,
bool extractQ)
486 RealScalar v1norm2 =
abs2(mat(2,0));
487 if(v1norm2 == RealScalar(0))
491 subdiag[0] = mat(1,0);
492 subdiag[1] = mat(2,1);
498 RealScalar beta =
sqrt(
abs2(mat(1,0)) + v1norm2);
499 RealScalar invBeta = RealScalar(1)/beta;
500 Scalar m01 = mat(1,0) * invBeta;
501 Scalar m02 = mat(2,0) * invBeta;
502 Scalar
q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
503 diag[1] = mat(1,1) + m02*
q;
504 diag[2] = mat(2,2) - m02*
q;
506 subdiag[1] = mat(2,1) - m01 *
q;
520 template<
typename MatrixType,
bool IsComplex>
521 struct tridiagonalization_inplace_selector<MatrixType,1,
IsComplex>
523 typedef typename MatrixType::Scalar Scalar;
525 template<
typename DiagonalType,
typename SubDiagonalType>
526 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&,
bool extractQ)
528 diag(0,0) =
real(mat(0,0));
530 mat(0,0) = Scalar(1);
541 template<
typename MatrixType>
struct TridiagonalizationMatrixTReturnType
542 :
public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
544 typedef typename MatrixType::Index Index;
550 TridiagonalizationMatrixTReturnType(
const MatrixType& mat) : m_matrix(mat) { }
552 template <
typename ResultType>
553 inline void evalTo(ResultType& result)
const
556 result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
557 result.diagonal() = m_matrix.diagonal();
558 result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
561 Index
rows()
const {
return m_matrix.rows(); }
562 Index
cols()
const {
return m_matrix.cols(); }
565 typename MatrixType::Nested m_matrix;
572 #endif // EIGEN_TRIDIAGONALIZATION_H