Eigen  3.3.3
Tridiagonalization.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_TRIDIAGONALIZATION_H
00012 #define EIGEN_TRIDIAGONALIZATION_H
00013 
00014 namespace Eigen { 
00015 
00016 namespace internal {
00017   
00018 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType;
00019 template<typename MatrixType>
00020 struct traits<TridiagonalizationMatrixTReturnType<MatrixType> >
00021   : public traits<typename MatrixType::PlainObject>
00022 {
00023   typedef typename MatrixType::PlainObject ReturnType; // FIXME shall it be a BandMatrix?
00024   enum { Flags = 0 };
00025 };
00026 
00027 template<typename MatrixType, typename CoeffVectorType>
00028 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
00029 }
00030 
00063 template<typename _MatrixType> class Tridiagonalization
00064 {
00065   public:
00066 
00068     typedef _MatrixType MatrixType;
00069 
00070     typedef typename MatrixType::Scalar Scalar;
00071     typedef typename NumTraits<Scalar>::Real RealScalar;
00072     typedef Eigen::Index Index; 
00073 
00074     enum {
00075       Size = MatrixType::RowsAtCompileTime,
00076       SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
00077       Options = MatrixType::Options,
00078       MaxSize = MatrixType::MaxRowsAtCompileTime,
00079       MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
00080     };
00081 
00082     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
00083     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
00084     typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
00085     typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
00086     typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
00087 
00088     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
00089               typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
00090               const Diagonal<const MatrixType>
00091             >::type DiagonalReturnType;
00092 
00093     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
00094               typename internal::add_const_on_value_type<typename Diagonal<const MatrixType, -1>::RealReturnType>::type,
00095               const Diagonal<const MatrixType, -1>
00096             >::type SubDiagonalReturnType;
00097 
00099     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType;
00100 
00113     explicit Tridiagonalization(Index size = Size==Dynamic ? 2 : Size)
00114       : m_matrix(size,size),
00115         m_hCoeffs(size > 1 ? size-1 : 1),
00116         m_isInitialized(false)
00117     {}
00118 
00129     template<typename InputType>
00130     explicit Tridiagonalization(const EigenBase<InputType>& matrix)
00131       : m_matrix(matrix.derived()),
00132         m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
00133         m_isInitialized(false)
00134     {
00135       internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
00136       m_isInitialized = true;
00137     }
00138 
00156     template<typename InputType>
00157     Tridiagonalization& compute(const EigenBase<InputType>& matrix)
00158     {
00159       m_matrix = matrix.derived();
00160       m_hCoeffs.resize(matrix.rows()-1, 1);
00161       internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
00162       m_isInitialized = true;
00163       return *this;
00164     }
00165 
00182     inline CoeffVectorType householderCoefficients() const
00183     {
00184       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00185       return m_hCoeffs;
00186     }
00187 
00219     inline const MatrixType& packedMatrix() const
00220     {
00221       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00222       return m_matrix;
00223     }
00224 
00240     HouseholderSequenceType matrixQ() const
00241     {
00242       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00243       return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
00244              .setLength(m_matrix.rows() - 1)
00245              .setShift(1);
00246     }
00247 
00265     MatrixTReturnType matrixT() const
00266     {
00267       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00268       return MatrixTReturnType(m_matrix.real());
00269     }
00270 
00284     DiagonalReturnType diagonal() const;
00285 
00296     SubDiagonalReturnType subDiagonal() const;
00297 
00298   protected:
00299 
00300     MatrixType m_matrix;
00301     CoeffVectorType m_hCoeffs;
00302     bool m_isInitialized;
00303 };
00304 
00305 template<typename MatrixType>
00306 typename Tridiagonalization<MatrixType>::DiagonalReturnType
00307 Tridiagonalization<MatrixType>::diagonal() const
00308 {
00309   eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00310   return m_matrix.diagonal().real();
00311 }
00312 
00313 template<typename MatrixType>
00314 typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
00315 Tridiagonalization<MatrixType>::subDiagonal() const
00316 {
00317   eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00318   return m_matrix.template diagonal<-1>().real();
00319 }
00320 
00321 namespace internal {
00322 
00346 template<typename MatrixType, typename CoeffVectorType>
00347 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
00348 {
00349   using numext::conj;
00350   typedef typename MatrixType::Scalar Scalar;
00351   typedef typename MatrixType::RealScalar RealScalar;
00352   Index n = matA.rows();
00353   eigen_assert(n==matA.cols());
00354   eigen_assert(n==hCoeffs.size()+1 || n==1);
00355   
00356   for (Index i = 0; i<n-1; ++i)
00357   {
00358     Index remainingSize = n-i-1;
00359     RealScalar beta;
00360     Scalar h;
00361     matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
00362 
00363     // Apply similarity transformation to remaining columns,
00364     // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
00365     matA.col(i).coeffRef(i+1) = 1;
00366 
00367     hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
00368                                   * (conj(h) * matA.col(i).tail(remainingSize)));
00369 
00370     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);
00371 
00372     matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
00373       .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), Scalar(-1));
00374 
00375     matA.col(i).coeffRef(i+1) = beta;
00376     hCoeffs.coeffRef(i) = h;
00377   }
00378 }
00379 
00380 // forward declaration, implementation at the end of this file
00381 template<typename MatrixType,
00382          int Size=MatrixType::ColsAtCompileTime,
00383          bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
00384 struct tridiagonalization_inplace_selector;
00385 
00426 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
00427 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00428 {
00429   eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
00430   tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
00431 }
00432 
00436 template<typename MatrixType, int Size, bool IsComplex>
00437 struct tridiagonalization_inplace_selector
00438 {
00439   typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
00440   typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
00441   template<typename DiagonalType, typename SubDiagonalType>
00442   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00443   {
00444     CoeffVectorType hCoeffs(mat.cols()-1);
00445     tridiagonalization_inplace(mat,hCoeffs);
00446     diag = mat.diagonal().real();
00447     subdiag = mat.template diagonal<-1>().real();
00448     if(extractQ)
00449       mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
00450             .setLength(mat.rows() - 1)
00451             .setShift(1);
00452   }
00453 };
00454 
00459 template<typename MatrixType>
00460 struct tridiagonalization_inplace_selector<MatrixType,3,false>
00461 {
00462   typedef typename MatrixType::Scalar Scalar;
00463   typedef typename MatrixType::RealScalar RealScalar;
00464 
00465   template<typename DiagonalType, typename SubDiagonalType>
00466   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00467   {
00468     using std::sqrt;
00469     const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
00470     diag[0] = mat(0,0);
00471     RealScalar v1norm2 = numext::abs2(mat(2,0));
00472     if(v1norm2 <= tol)
00473     {
00474       diag[1] = mat(1,1);
00475       diag[2] = mat(2,2);
00476       subdiag[0] = mat(1,0);
00477       subdiag[1] = mat(2,1);
00478       if (extractQ)
00479         mat.setIdentity();
00480     }
00481     else
00482     {
00483       RealScalar beta = sqrt(numext::abs2(mat(1,0)) + v1norm2);
00484       RealScalar invBeta = RealScalar(1)/beta;
00485       Scalar m01 = mat(1,0) * invBeta;
00486       Scalar m02 = mat(2,0) * invBeta;
00487       Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
00488       diag[1] = mat(1,1) + m02*q;
00489       diag[2] = mat(2,2) - m02*q;
00490       subdiag[0] = beta;
00491       subdiag[1] = mat(2,1) - m01 * q;
00492       if (extractQ)
00493       {
00494         mat << 1,   0,    0,
00495                0, m01,  m02,
00496                0, m02, -m01;
00497       }
00498     }
00499   }
00500 };
00501 
00505 template<typename MatrixType, bool IsComplex>
00506 struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
00507 {
00508   typedef typename MatrixType::Scalar Scalar;
00509 
00510   template<typename DiagonalType, typename SubDiagonalType>
00511   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ)
00512   {
00513     diag(0,0) = numext::real(mat(0,0));
00514     if(extractQ)
00515       mat(0,0) = Scalar(1);
00516   }
00517 };
00518 
00526 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
00527 : public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
00528 {
00529   public:
00534     TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { }
00535 
00536     template <typename ResultType>
00537     inline void evalTo(ResultType& result) const
00538     {
00539       result.setZero();
00540       result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
00541       result.diagonal() = m_matrix.diagonal();
00542       result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
00543     }
00544 
00545     Index rows() const { return m_matrix.rows(); }
00546     Index cols() const { return m_matrix.cols(); }
00547 
00548   protected:
00549     typename MatrixType::Nested m_matrix;
00550 };
00551 
00552 } // end namespace internal
00553 
00554 } // end namespace Eigen
00555 
00556 #endif // EIGEN_TRIDIAGONALIZATION_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends