![]() |
Eigen
3.3.3
|
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