Tridiagonalization.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
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     // Apply similarity transformation to remaining columns,
00378     // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
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 // forward declaration, implementation at the end of this file
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   //Index n = mat.rows();
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 } // end namespace internal
00569 
00570 } // end namespace Eigen
00571 
00572 #endif // EIGEN_TRIDIAGONALIZATION_H