00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef EIGEN_TRIDIAGONALIZATION_H
00027 #define EIGEN_TRIDIAGONALIZATION_H
00028
00029 namespace Eigen {
00030
00031 namespace internal {
00032
00033 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType;
00034 template<typename MatrixType>
00035 struct traits<TridiagonalizationMatrixTReturnType<MatrixType> >
00036 {
00037 typedef typename MatrixType::PlainObject ReturnType;
00038 };
00039
00040 template<typename MatrixType, typename CoeffVectorType>
00041 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
00042 }
00043
00076 template<typename _MatrixType> class Tridiagonalization
00077 {
00078 public:
00079
00081 typedef _MatrixType MatrixType;
00082
00083 typedef typename MatrixType::Scalar Scalar;
00084 typedef typename NumTraits<Scalar>::Real RealScalar;
00085 typedef typename MatrixType::Index Index;
00086
00087 enum {
00088 Size = MatrixType::RowsAtCompileTime,
00089 SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
00090 Options = MatrixType::Options,
00091 MaxSize = MatrixType::MaxRowsAtCompileTime,
00092 MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
00093 };
00094
00095 typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
00096 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
00097 typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
00098 typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
00099 typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
00100
00101 typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
00102 typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
00103 const Diagonal<const MatrixType>
00104 >::type DiagonalReturnType;
00105
00106 typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
00107 typename internal::add_const_on_value_type<typename Diagonal<
00108 Block<const MatrixType,SizeMinusOne,SizeMinusOne> >::RealReturnType>::type,
00109 const Diagonal<
00110 Block<const MatrixType,SizeMinusOne,SizeMinusOne> >
00111 >::type SubDiagonalReturnType;
00112
00114 typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
00115
00128 Tridiagonalization(Index size = Size==Dynamic ? 2 : Size)
00129 : m_matrix(size,size),
00130 m_hCoeffs(size > 1 ? size-1 : 1),
00131 m_isInitialized(false)
00132 {}
00133
00144 Tridiagonalization(const MatrixType& matrix)
00145 : m_matrix(matrix),
00146 m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
00147 m_isInitialized(false)
00148 {
00149 internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
00150 m_isInitialized = true;
00151 }
00152
00170 Tridiagonalization& compute(const MatrixType& matrix)
00171 {
00172 m_matrix = matrix;
00173 m_hCoeffs.resize(matrix.rows()-1, 1);
00174 internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
00175 m_isInitialized = true;
00176 return *this;
00177 }
00178
00195 inline CoeffVectorType householderCoefficients() const
00196 {
00197 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00198 return m_hCoeffs;
00199 }
00200
00232 inline const MatrixType& packedMatrix() const
00233 {
00234 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00235 return m_matrix;
00236 }
00237
00253 HouseholderSequenceType matrixQ() const
00254 {
00255 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00256 return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
00257 .setLength(m_matrix.rows() - 1)
00258 .setShift(1);
00259 }
00260
00278 MatrixTReturnType matrixT() const
00279 {
00280 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00281 return MatrixTReturnType(m_matrix.real());
00282 }
00283
00297 DiagonalReturnType diagonal() const;
00298
00309 SubDiagonalReturnType subDiagonal() const;
00310
00311 protected:
00312
00313 MatrixType m_matrix;
00314 CoeffVectorType m_hCoeffs;
00315 bool m_isInitialized;
00316 };
00317
00318 template<typename MatrixType>
00319 typename Tridiagonalization<MatrixType>::DiagonalReturnType
00320 Tridiagonalization<MatrixType>::diagonal() const
00321 {
00322 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00323 return m_matrix.diagonal();
00324 }
00325
00326 template<typename MatrixType>
00327 typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
00328 Tridiagonalization<MatrixType>::subDiagonal() const
00329 {
00330 eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00331 Index n = m_matrix.rows();
00332 return Block<const MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal();
00333 }
00334
00335 namespace internal {
00336
00360 template<typename MatrixType, typename CoeffVectorType>
00361 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
00362 {
00363 typedef typename MatrixType::Index Index;
00364 typedef typename MatrixType::Scalar Scalar;
00365 typedef typename MatrixType::RealScalar RealScalar;
00366 Index n = matA.rows();
00367 eigen_assert(n==matA.cols());
00368 eigen_assert(n==hCoeffs.size()+1 || n==1);
00369
00370 for (Index i = 0; i<n-1; ++i)
00371 {
00372 Index remainingSize = n-i-1;
00373 RealScalar beta;
00374 Scalar h;
00375 matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
00376
00377
00378
00379 matA.col(i).coeffRef(i+1) = 1;
00380
00381 hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
00382 * (conj(h) * matA.col(i).tail(remainingSize)));
00383
00384 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);
00385
00386 matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
00387 .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
00388
00389 matA.col(i).coeffRef(i+1) = beta;
00390 hCoeffs.coeffRef(i) = h;
00391 }
00392 }
00393
00394
00395 template<typename MatrixType,
00396 int Size=MatrixType::ColsAtCompileTime,
00397 bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
00398 struct tridiagonalization_inplace_selector;
00399
00440 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
00441 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00442 {
00443 typedef typename MatrixType::Index Index;
00444
00445 eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
00446 tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
00447 }
00448
00452 template<typename MatrixType, int Size, bool IsComplex>
00453 struct tridiagonalization_inplace_selector
00454 {
00455 typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
00456 typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
00457 typedef typename MatrixType::Index Index;
00458 template<typename DiagonalType, typename SubDiagonalType>
00459 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00460 {
00461 CoeffVectorType hCoeffs(mat.cols()-1);
00462 tridiagonalization_inplace(mat,hCoeffs);
00463 diag = mat.diagonal().real();
00464 subdiag = mat.template diagonal<-1>().real();
00465 if(extractQ)
00466 mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
00467 .setLength(mat.rows() - 1)
00468 .setShift(1);
00469 }
00470 };
00471
00476 template<typename MatrixType>
00477 struct tridiagonalization_inplace_selector<MatrixType,3,false>
00478 {
00479 typedef typename MatrixType::Scalar Scalar;
00480 typedef typename MatrixType::RealScalar RealScalar;
00481
00482 template<typename DiagonalType, typename SubDiagonalType>
00483 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00484 {
00485 diag[0] = mat(0,0);
00486 RealScalar v1norm2 = abs2(mat(2,0));
00487 if(v1norm2 == RealScalar(0))
00488 {
00489 diag[1] = mat(1,1);
00490 diag[2] = mat(2,2);
00491 subdiag[0] = mat(1,0);
00492 subdiag[1] = mat(2,1);
00493 if (extractQ)
00494 mat.setIdentity();
00495 }
00496 else
00497 {
00498 RealScalar beta = sqrt(abs2(mat(1,0)) + v1norm2);
00499 RealScalar invBeta = RealScalar(1)/beta;
00500 Scalar m01 = mat(1,0) * invBeta;
00501 Scalar m02 = mat(2,0) * invBeta;
00502 Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
00503 diag[1] = mat(1,1) + m02*q;
00504 diag[2] = mat(2,2) - m02*q;
00505 subdiag[0] = beta;
00506 subdiag[1] = mat(2,1) - m01 * q;
00507 if (extractQ)
00508 {
00509 mat << 1, 0, 0,
00510 0, m01, m02,
00511 0, m02, -m01;
00512 }
00513 }
00514 }
00515 };
00516
00520 template<typename MatrixType, bool IsComplex>
00521 struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
00522 {
00523 typedef typename MatrixType::Scalar Scalar;
00524
00525 template<typename DiagonalType, typename SubDiagonalType>
00526 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ)
00527 {
00528 diag(0,0) = real(mat(0,0));
00529 if(extractQ)
00530 mat(0,0) = Scalar(1);
00531 }
00532 };
00533
00541 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
00542 : public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
00543 {
00544 typedef typename MatrixType::Index Index;
00545 public:
00550 TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { }
00551
00552 template <typename ResultType>
00553 inline void evalTo(ResultType& result) const
00554 {
00555 result.setZero();
00556 result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
00557 result.diagonal() = m_matrix.diagonal();
00558 result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
00559 }
00560
00561 Index rows() const { return m_matrix.rows(); }
00562 Index cols() const { return m_matrix.cols(); }
00563
00564 protected:
00565 typename MatrixType::Nested m_matrix;
00566 };
00567
00568 }
00569
00570 }
00571
00572 #endif // EIGEN_TRIDIAGONALIZATION_H