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 // 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 {
00022   typedef typename MatrixType::PlainObject ReturnType;
00023 };
00024 
00025 template<typename MatrixType, typename CoeffVectorType>
00026 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
00027 }
00028 
00061 template<typename _MatrixType> class Tridiagonalization
00062 {
00063   public:
00064 
00066     typedef _MatrixType MatrixType;
00067 
00068     typedef typename MatrixType::Scalar Scalar;
00069     typedef typename NumTraits<Scalar>::Real RealScalar;
00070     typedef typename MatrixType::Index Index;
00071 
00072     enum {
00073       Size = MatrixType::RowsAtCompileTime,
00074       SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
00075       Options = MatrixType::Options,
00076       MaxSize = MatrixType::MaxRowsAtCompileTime,
00077       MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
00078     };
00079 
00080     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
00081     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
00082     typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
00083     typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
00084     typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
00085 
00086     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
00087               typename internal::add_const_on_value_type<typename Diagonal<const MatrixType>::RealReturnType>::type,
00088               const Diagonal<const MatrixType>
00089             >::type DiagonalReturnType;
00090 
00091     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
00092               typename internal::add_const_on_value_type<typename Diagonal<
00093                 Block<const MatrixType,SizeMinusOne,SizeMinusOne> >::RealReturnType>::type,
00094               const Diagonal<
00095                 Block<const MatrixType,SizeMinusOne,SizeMinusOne> >
00096             >::type SubDiagonalReturnType;
00097 
00099     typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
00100 
00113     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     Tridiagonalization(const MatrixType& matrix)
00130       : m_matrix(matrix),
00131         m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
00132         m_isInitialized(false)
00133     {
00134       internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
00135       m_isInitialized = true;
00136     }
00137 
00155     Tridiagonalization& compute(const MatrixType& matrix)
00156     {
00157       m_matrix = matrix;
00158       m_hCoeffs.resize(matrix.rows()-1, 1);
00159       internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
00160       m_isInitialized = true;
00161       return *this;
00162     }
00163 
00180     inline CoeffVectorType householderCoefficients() const
00181     {
00182       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00183       return m_hCoeffs;
00184     }
00185 
00217     inline const MatrixType& packedMatrix() const
00218     {
00219       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00220       return m_matrix;
00221     }
00222 
00238     HouseholderSequenceType matrixQ() const
00239     {
00240       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00241       return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
00242              .setLength(m_matrix.rows() - 1)
00243              .setShift(1);
00244     }
00245 
00263     MatrixTReturnType matrixT() const
00264     {
00265       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00266       return MatrixTReturnType(m_matrix.real());
00267     }
00268 
00282     DiagonalReturnType diagonal() const;
00283 
00294     SubDiagonalReturnType subDiagonal() const;
00295 
00296   protected:
00297 
00298     MatrixType m_matrix;
00299     CoeffVectorType m_hCoeffs;
00300     bool m_isInitialized;
00301 };
00302 
00303 template<typename MatrixType>
00304 typename Tridiagonalization<MatrixType>::DiagonalReturnType
00305 Tridiagonalization<MatrixType>::diagonal() const
00306 {
00307   eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00308   return m_matrix.diagonal();
00309 }
00310 
00311 template<typename MatrixType>
00312 typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
00313 Tridiagonalization<MatrixType>::subDiagonal() const
00314 {
00315   eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00316   Index n = m_matrix.rows();
00317   return Block<const MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal();
00318 }
00319 
00320 namespace internal {
00321 
00345 template<typename MatrixType, typename CoeffVectorType>
00346 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
00347 {
00348   typedef typename MatrixType::Index Index;
00349   typedef typename MatrixType::Scalar Scalar;
00350   typedef typename MatrixType::RealScalar RealScalar;
00351   Index n = matA.rows();
00352   eigen_assert(n==matA.cols());
00353   eigen_assert(n==hCoeffs.size()+1 || n==1);
00354   
00355   for (Index i = 0; i<n-1; ++i)
00356   {
00357     Index remainingSize = n-i-1;
00358     RealScalar beta;
00359     Scalar h;
00360     matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
00361 
00362     // Apply similarity transformation to remaining columns,
00363     // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
00364     matA.col(i).coeffRef(i+1) = 1;
00365 
00366     hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
00367                                   * (conj(h) * matA.col(i).tail(remainingSize)));
00368 
00369     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);
00370 
00371     matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
00372       .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
00373 
00374     matA.col(i).coeffRef(i+1) = beta;
00375     hCoeffs.coeffRef(i) = h;
00376   }
00377 }
00378 
00379 // forward declaration, implementation at the end of this file
00380 template<typename MatrixType,
00381          int Size=MatrixType::ColsAtCompileTime,
00382          bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
00383 struct tridiagonalization_inplace_selector;
00384 
00425 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
00426 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00427 {
00428   typedef typename MatrixType::Index Index;
00429   //Index n = mat.rows();
00430   eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
00431   tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
00432 }
00433 
00437 template<typename MatrixType, int Size, bool IsComplex>
00438 struct tridiagonalization_inplace_selector
00439 {
00440   typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
00441   typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
00442   typedef typename MatrixType::Index Index;
00443   template<typename DiagonalType, typename SubDiagonalType>
00444   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00445   {
00446     CoeffVectorType hCoeffs(mat.cols()-1);
00447     tridiagonalization_inplace(mat,hCoeffs);
00448     diag = mat.diagonal().real();
00449     subdiag = mat.template diagonal<-1>().real();
00450     if(extractQ)
00451       mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
00452             .setLength(mat.rows() - 1)
00453             .setShift(1);
00454   }
00455 };
00456 
00461 template<typename MatrixType>
00462 struct tridiagonalization_inplace_selector<MatrixType,3,false>
00463 {
00464   typedef typename MatrixType::Scalar Scalar;
00465   typedef typename MatrixType::RealScalar RealScalar;
00466 
00467   template<typename DiagonalType, typename SubDiagonalType>
00468   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00469   {
00470     diag[0] = mat(0,0);
00471     RealScalar v1norm2 = abs2(mat(2,0));
00472     if(v1norm2 == RealScalar(0))
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(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) = 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     typedef typename MatrixType::Index Index;
00530   public:
00535     TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { }
00536 
00537     template <typename ResultType>
00538     inline void evalTo(ResultType& result) const
00539     {
00540       result.setZero();
00541       result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
00542       result.diagonal() = m_matrix.diagonal();
00543       result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
00544     }
00545 
00546     Index rows() const { return m_matrix.rows(); }
00547     Index cols() const { return m_matrix.cols(); }
00548 
00549   protected:
00550     typename MatrixType::Nested m_matrix;
00551 };
00552 
00553 } // end namespace internal
00554 
00555 } // end namespace Eigen
00556 
00557 #endif // EIGEN_TRIDIAGONALIZATION_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Wed Sep 16 2015 07:12:35