$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> 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 using numext::conj; 00349 typedef typename MatrixType::Index Index; 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)*Scalar(-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), -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 typedef typename MatrixType::Index Index; 00442 template<typename DiagonalType, typename SubDiagonalType> 00443 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) 00444 { 00445 CoeffVectorType hCoeffs(mat.cols()-1); 00446 tridiagonalization_inplace(mat,hCoeffs); 00447 diag = mat.diagonal().real(); 00448 subdiag = mat.template diagonal<-1>().real(); 00449 if(extractQ) 00450 mat = HouseholderSequenceType(mat, hCoeffs.conjugate()) 00451 .setLength(mat.rows() - 1) 00452 .setShift(1); 00453 } 00454 }; 00455 00460 template<typename MatrixType> 00461 struct tridiagonalization_inplace_selector<MatrixType,3,false> 00462 { 00463 typedef typename MatrixType::Scalar Scalar; 00464 typedef typename MatrixType::RealScalar RealScalar; 00465 00466 template<typename DiagonalType, typename SubDiagonalType> 00467 static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) 00468 { 00469 using std::sqrt; 00470 diag[0] = mat(0,0); 00471 RealScalar v1norm2 = numext::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(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 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