11 #ifndef EIGEN_TRIDIAGONALIZATION_H 12 #define EIGEN_TRIDIAGONALIZATION_H 19 template<
typename MatrixType>
25 template<
typename MatrixType,
typename CoeffVectorType>
68 typedef typename MatrixType::Scalar
Scalar;
70 typedef typename MatrixType::Index
Index;
73 Size = MatrixType::RowsAtCompileTime,
76 MaxSize = MatrixType::MaxRowsAtCompileTime,
77 MaxSizeMinusOne = MaxSize ==
Dynamic ?
Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
91 typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
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)
135 m_isInitialized =
true;
158 m_hCoeffs.resize(matrix.rows()-1, 1);
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.");
241 return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
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;
294 SubDiagonalReturnType subDiagonal()
const;
303 template<
typename MatrixType>
307 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
308 return m_matrix.diagonal();
311 template<
typename MatrixType>
315 eigen_assert(m_isInitialized &&
"Tridiagonalization is not initialized.");
316 Index n = m_matrix.rows();
345 template<
typename MatrixType,
typename CoeffVectorType>
349 typedef typename MatrixType::Index Index;
350 typedef typename MatrixType::Scalar Scalar;
351 typedef typename MatrixType::RealScalar RealScalar;
352 Index n = matA.rows();
356 for (Index i = 0; i<n-1; ++i)
358 Index remainingSize = n-i-1;
361 matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
365 matA.col(i).coeffRef(i+1) = 1;
367 hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
368 * (
conj(h) * matA.col(i).tail(remainingSize)));
370 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);
372 matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
373 .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
375 matA.col(i).coeffRef(i+1) = beta;
376 hCoeffs.coeffRef(i) = h;
381 template<
typename MatrixType,
382 int Size=MatrixType::ColsAtCompileTime,
426 template<
typename MatrixType,
typename DiagonalType,
typename SubDiagonalType>
429 eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
436 template<
typename MatrixType,
int Size,
bool IsComplex>
441 typedef typename MatrixType::Index
Index;
442 template<
typename DiagonalType,
typename SubDiagonalType>
443 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag,
bool extractQ)
445 CoeffVectorType hCoeffs(mat.cols()-1);
447 diag = mat.diagonal().real();
448 subdiag = mat.template diagonal<-1>().
real();
450 mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
451 .setLength(mat.rows() - 1)
460 template<
typename MatrixType>
463 typedef typename MatrixType::Scalar
Scalar;
466 template<
typename DiagonalType,
typename SubDiagonalType>
467 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag,
bool extractQ)
472 if(v1norm2 == RealScalar(0))
476 subdiag[0] = mat(1,0);
477 subdiag[1] = mat(2,1);
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>
508 typedef typename MatrixType::Scalar
Scalar;
510 template<
typename DiagonalType,
typename SubDiagonalType>
511 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&,
bool extractQ)
515 mat(0,0) = Scalar(1);
527 :
public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
529 typedef typename MatrixType::Index
Index;
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(); }
557 #endif // EIGEN_TRIDIAGONALIZATION_H HouseholderSequence< MatrixType, typename internal::remove_all< typename CoeffVectorType::ConjugateReturnType >::type > HouseholderSequenceType
Return type of matrixQ()
Matrix< Scalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 > CoeffVectorType
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
Tridiagonalization< MatrixType >::HouseholderSequenceType HouseholderSequenceType
IntermediateState sqrt(const Expression &arg)
internal::remove_all< typename MatrixType::RealReturnType >::type MatrixTypeRealView
static void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &subdiag, bool extractQ)
HouseholderSequenceType matrixQ() const
Returns the unitary matrix Q in the decomposition.
static void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &, bool extractQ)
internal::conditional< NumTraits< Scalar >::IsComplex, const CwiseUnaryOp< internal::scalar_real_op< Scalar >, const Derived >, const Derived & >::type RealReturnType
iterative scaling algorithm to equilibrate rows and column norms in matrices
MatrixTReturnType matrixT() const
Returns an expression of the tridiagonal matrix T in the decomposition.
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Tridiagonal decomposition of a selfadjoint matrix.
internal::conditional< NumTraits< Scalar >::IsComplex, typename internal::add_const_on_value_type< typename Diagonal< Block< const MatrixType, SizeMinusOne, SizeMinusOne > >::RealReturnType >::type, const Diagonal< Block< const MatrixType, SizeMinusOne, SizeMinusOne > > >::type SubDiagonalReturnType
MatrixType::Nested m_matrix
Tridiagonalization(const MatrixType &matrix)
Constructor; computes tridiagonal decomposition of given matrix.
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs2_op< Scalar >, const Derived > abs2() const
Sequence of Householder reflections acting on subspaces with decreasing size.
RealReturnType real() const
const MatrixType & packedMatrix() const
Returns the internal representation of the decomposition.
Tridiagonalization< MatrixType >::CoeffVectorType CoeffVectorType
Tridiagonalization & compute(const MatrixType &matrix)
Computes tridiagonal decomposition of given matrix.
CoeffVectorType m_hCoeffs
DiagonalReturnType diagonal() const
Returns the diagonal of the tridiagonal matrix T in the decomposition.
SubDiagonalReturnType subDiagonal() const
Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
Provides a generic way to set and pass user-specified options.
void evalTo(ResultType &result) const
MatrixType::Scalar Scalar
MatrixType::Scalar Scalar
void tridiagonalization_inplace(MatrixType &matA, CoeffVectorType &hCoeffs)
internal::TridiagonalizationMatrixTReturnType< MatrixTypeRealView > MatrixTReturnType
Expression of a fixed-size or dynamic-size block.
internal::plain_col_type< MatrixType, RealScalar >::type DiagonalType
TridiagonalizationMatrixTReturnType(const MatrixType &mat)
Constructor.
MatrixType::RealScalar RealScalar
CoeffVectorType householderCoefficients() const
Returns the Householder coefficients.
internal::conditional< NumTraits< Scalar >::IsComplex, typename internal::add_const_on_value_type< typename Diagonal< const MatrixType >::RealReturnType >::type, const Diagonal< const MatrixType > >::type DiagonalReturnType
Tridiagonalization(Index size=Size==Dynamic?2:Size)
Default constructor.
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
MatrixType::Scalar Scalar
Matrix< RealScalar, SizeMinusOne, 1, Options &~RowMajor, MaxSizeMinusOne, 1 > SubDiagonalType
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
const AutoDiffScalar< DerType > & real(const AutoDiffScalar< DerType > &x)
NumTraits< Scalar >::Real RealScalar
MatrixType::PlainObject ReturnType
static void run(MatrixType &mat, DiagonalType &diag, SubDiagonalType &subdiag, bool extractQ)