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 internal {
00030   
00031 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType;
00032 template<typename MatrixType>
00033 struct traits<TridiagonalizationMatrixTReturnType<MatrixType> >
00034 {
00035   typedef typename MatrixType::PlainObject ReturnType;
00036 };
00037 
00038 template<typename MatrixType, typename CoeffVectorType>
00039 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
00040 }
00041 
00074 template<typename _MatrixType> class Tridiagonalization
00075 {
00076   public:
00077 
00079     typedef _MatrixType MatrixType;
00080 
00081     typedef typename MatrixType::Scalar Scalar;
00082     typedef typename NumTraits<Scalar>::Real RealScalar;
00083     typedef typename MatrixType::Index Index;
00084 
00085     enum {
00086       Size = MatrixType::RowsAtCompileTime,
00087       SizeMinusOne = Size == Dynamic ? Dynamic : (Size > 1 ? Size - 1 : 1),
00088       Options = MatrixType::Options,
00089       MaxSize = MatrixType::MaxRowsAtCompileTime,
00090       MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : (MaxSize > 1 ? MaxSize - 1 : 1)
00091     };
00092 
00093     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
00094     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
00095     typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
00096     typedef typename internal::remove_all<typename MatrixType::RealReturnType>::type MatrixTypeRealView;
00097     typedef internal::TridiagonalizationMatrixTReturnType<MatrixTypeRealView> MatrixTReturnType;
00098 
00099     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
00100               const typename Diagonal<const MatrixType>::RealReturnType,
00101               const Diagonal<const MatrixType>
00102             >::type DiagonalReturnType;
00103 
00104     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
00105               const typename Diagonal<
00106                 Block<const MatrixType,SizeMinusOne,SizeMinusOne> >::RealReturnType,
00107               const Diagonal<
00108                 Block<const MatrixType,SizeMinusOne,SizeMinusOne> >
00109             >::type SubDiagonalReturnType;
00110 
00112     typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
00113 
00126     Tridiagonalization(Index size = Size==Dynamic ? 2 : Size)
00127       : m_matrix(size,size),
00128         m_hCoeffs(size > 1 ? size-1 : 1),
00129         m_isInitialized(false)
00130     {}
00131 
00142     Tridiagonalization(const MatrixType& matrix)
00143       : m_matrix(matrix),
00144         m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
00145         m_isInitialized(false)
00146     {
00147       internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
00148       m_isInitialized = true;
00149     }
00150 
00168     Tridiagonalization& compute(const MatrixType& matrix)
00169     {
00170       m_matrix = matrix;
00171       m_hCoeffs.resize(matrix.rows()-1, 1);
00172       internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
00173       m_isInitialized = true;
00174       return *this;
00175     }
00176 
00193     inline CoeffVectorType householderCoefficients() const
00194     {
00195       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00196       return m_hCoeffs;
00197     }
00198 
00230     inline const MatrixType& packedMatrix() const
00231     {
00232       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00233       return m_matrix;
00234     }
00235 
00251     HouseholderSequenceType matrixQ() const
00252     {
00253       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00254       return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
00255              .setLength(m_matrix.rows() - 1)
00256              .setShift(1);
00257     }
00258 
00276     MatrixTReturnType matrixT() const
00277     {
00278       eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00279       return MatrixTReturnType(m_matrix.real());
00280     }
00281 
00295     DiagonalReturnType diagonal() const;
00296 
00307     SubDiagonalReturnType subDiagonal() const;
00308 
00309   protected:
00310 
00311     MatrixType m_matrix;
00312     CoeffVectorType m_hCoeffs;
00313     bool m_isInitialized;
00314 };
00315 
00316 template<typename MatrixType>
00317 typename Tridiagonalization<MatrixType>::DiagonalReturnType
00318 Tridiagonalization<MatrixType>::diagonal() const
00319 {
00320   eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00321   return m_matrix.diagonal();
00322 }
00323 
00324 template<typename MatrixType>
00325 typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
00326 Tridiagonalization<MatrixType>::subDiagonal() const
00327 {
00328   eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
00329   Index n = m_matrix.rows();
00330   return Block<const MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal();
00331 }
00332 
00333 namespace internal {
00334 
00358 template<typename MatrixType, typename CoeffVectorType>
00359 void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
00360 {
00361   typedef typename MatrixType::Index Index;
00362   typedef typename MatrixType::Scalar Scalar;
00363   typedef typename MatrixType::RealScalar RealScalar;
00364   Index n = matA.rows();
00365   eigen_assert(n==matA.cols());
00366   eigen_assert(n==hCoeffs.size()+1 || n==1);
00367   
00368   for (Index i = 0; i<n-1; ++i)
00369   {
00370     Index remainingSize = n-i-1;
00371     RealScalar beta;
00372     Scalar h;
00373     matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
00374 
00375     // Apply similarity transformation to remaining columns,
00376     // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
00377     matA.col(i).coeffRef(i+1) = 1;
00378 
00379     hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
00380                                   * (conj(h) * matA.col(i).tail(remainingSize)));
00381 
00382     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);
00383 
00384     matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
00385       .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
00386 
00387     matA.col(i).coeffRef(i+1) = beta;
00388     hCoeffs.coeffRef(i) = h;
00389   }
00390 }
00391 
00392 // forward declaration, implementation at the end of this file
00393 template<typename MatrixType,
00394          int Size=MatrixType::ColsAtCompileTime,
00395          bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
00396 struct tridiagonalization_inplace_selector;
00397 
00438 template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
00439 void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00440 {
00441   typedef typename MatrixType::Index Index;
00442   //Index n = mat.rows();
00443   eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
00444   tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
00445 }
00446 
00450 template<typename MatrixType, int Size, bool IsComplex>
00451 struct tridiagonalization_inplace_selector
00452 {
00453   typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
00454   typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
00455   typedef typename MatrixType::Index Index;
00456   template<typename DiagonalType, typename SubDiagonalType>
00457   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00458   {
00459     CoeffVectorType hCoeffs(mat.cols()-1);
00460     tridiagonalization_inplace(mat,hCoeffs);
00461     diag = mat.diagonal().real();
00462     subdiag = mat.template diagonal<-1>().real();
00463     if(extractQ)
00464       mat = HouseholderSequenceType(mat, hCoeffs.conjugate())
00465             .setLength(mat.rows() - 1)
00466             .setShift(1);
00467   }
00468 };
00469 
00474 template<typename MatrixType>
00475 struct tridiagonalization_inplace_selector<MatrixType,3,false>
00476 {
00477   typedef typename MatrixType::Scalar Scalar;
00478   typedef typename MatrixType::RealScalar RealScalar;
00479 
00480   template<typename DiagonalType, typename SubDiagonalType>
00481   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
00482   {
00483     diag[0] = mat(0,0);
00484     RealScalar v1norm2 = abs2(mat(2,0));
00485     if(v1norm2 == RealScalar(0))
00486     {
00487       diag[1] = mat(1,1);
00488       diag[2] = mat(2,2);
00489       subdiag[0] = mat(1,0);
00490       subdiag[1] = mat(2,1);
00491       if (extractQ)
00492         mat.setIdentity();
00493     }
00494     else
00495     {
00496       RealScalar beta = sqrt(abs2(mat(1,0)) + v1norm2);
00497       RealScalar invBeta = RealScalar(1)/beta;
00498       Scalar m01 = mat(1,0) * invBeta;
00499       Scalar m02 = mat(2,0) * invBeta;
00500       Scalar q = RealScalar(2)*m01*mat(2,1) + m02*(mat(2,2) - mat(1,1));
00501       diag[1] = mat(1,1) + m02*q;
00502       diag[2] = mat(2,2) - m02*q;
00503       subdiag[0] = beta;
00504       subdiag[1] = mat(2,1) - m01 * q;
00505       if (extractQ)
00506       {
00507         mat << 1,   0,    0,
00508                0, m01,  m02,
00509                0, m02, -m01;
00510       }
00511     }
00512   }
00513 };
00514 
00518 template<typename MatrixType, bool IsComplex>
00519 struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
00520 {
00521   typedef typename MatrixType::Scalar Scalar;
00522 
00523   template<typename DiagonalType, typename SubDiagonalType>
00524   static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ)
00525   {
00526     diag(0,0) = real(mat(0,0));
00527     if(extractQ)
00528       mat(0,0) = Scalar(1);
00529   }
00530 };
00531 
00539 template<typename MatrixType> struct TridiagonalizationMatrixTReturnType
00540 : public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> >
00541 {
00542     typedef typename MatrixType::Index Index;
00543   public:
00548     TridiagonalizationMatrixTReturnType(const MatrixType& mat) : m_matrix(mat) { }
00549 
00550     template <typename ResultType>
00551     inline void evalTo(ResultType& result) const
00552     {
00553       result.setZero();
00554       result.template diagonal<1>() = m_matrix.template diagonal<-1>().conjugate();
00555       result.diagonal() = m_matrix.diagonal();
00556       result.template diagonal<-1>() = m_matrix.template diagonal<-1>();
00557     }
00558 
00559     Index rows() const { return m_matrix.rows(); }
00560     Index cols() const { return m_matrix.cols(); }
00561 
00562   protected:
00563     const typename MatrixType::Nested m_matrix;
00564 };
00565 
00566 } // end namespace internal
00567 
00568 #endif // EIGEN_TRIDIAGONALIZATION_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:55