$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) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2010,2012 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_GENERALIZEDEIGENSOLVER_H 00012 #define EIGEN_GENERALIZEDEIGENSOLVER_H 00013 00014 #include "./RealQZ.h" 00015 00016 namespace Eigen { 00017 00057 template<typename _MatrixType> class GeneralizedEigenSolver 00058 { 00059 public: 00060 00062 typedef _MatrixType MatrixType; 00063 00064 enum { 00065 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00066 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00067 Options = MatrixType::Options, 00068 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00069 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00070 }; 00071 00073 typedef typename MatrixType::Scalar Scalar; 00074 typedef typename NumTraits<Scalar>::Real RealScalar; 00075 typedef typename MatrixType::Index Index; 00076 00083 typedef std::complex<RealScalar> ComplexScalar; 00084 00090 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType; 00091 00097 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType; 00098 00101 typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> EigenvalueType; 00102 00108 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType; 00109 00117 GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {} 00118 00125 GeneralizedEigenSolver(Index size) 00126 : m_eivec(size, size), 00127 m_alphas(size), 00128 m_betas(size), 00129 m_isInitialized(false), 00130 m_eigenvectorsOk(false), 00131 m_realQZ(size), 00132 m_matS(size, size), 00133 m_tmp(size) 00134 {} 00135 00148 GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true) 00149 : m_eivec(A.rows(), A.cols()), 00150 m_alphas(A.cols()), 00151 m_betas(A.cols()), 00152 m_isInitialized(false), 00153 m_eigenvectorsOk(false), 00154 m_realQZ(A.cols()), 00155 m_matS(A.rows(), A.cols()), 00156 m_tmp(A.cols()) 00157 { 00158 compute(A, B, computeEigenvectors); 00159 } 00160 00161 /* \brief Returns the computed generalized eigenvectors. 00162 * 00163 * \returns %Matrix whose columns are the (possibly complex) eigenvectors. 00164 * 00165 * \pre Either the constructor 00166 * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function 00167 * compute(const MatrixType&, const MatrixType& bool) has been called before, and 00168 * \p computeEigenvectors was set to true (the default). 00169 * 00170 * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding 00171 * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The 00172 * eigenvectors are normalized to have (Euclidean) norm equal to one. The 00173 * matrix returned by this function is the matrix \f$ V \f$ in the 00174 * generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists. 00175 * 00176 * \sa eigenvalues() 00177 */ 00178 // EigenvectorsType eigenvectors() const; 00179 00198 EigenvalueType eigenvalues() const 00199 { 00200 eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); 00201 return EigenvalueType(m_alphas,m_betas); 00202 } 00203 00209 ComplexVectorType alphas() const 00210 { 00211 eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); 00212 return m_alphas; 00213 } 00214 00220 VectorType betas() const 00221 { 00222 eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); 00223 return m_betas; 00224 } 00225 00249 GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true); 00250 00251 ComputationInfo info() const 00252 { 00253 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00254 return m_realQZ.info(); 00255 } 00256 00259 GeneralizedEigenSolver& setMaxIterations(Index maxIters) 00260 { 00261 m_realQZ.setMaxIterations(maxIters); 00262 return *this; 00263 } 00264 00265 protected: 00266 00267 static void check_template_parameters() 00268 { 00269 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00270 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); 00271 } 00272 00273 MatrixType m_eivec; 00274 ComplexVectorType m_alphas; 00275 VectorType m_betas; 00276 bool m_isInitialized; 00277 bool m_eigenvectorsOk; 00278 RealQZ<MatrixType> m_realQZ; 00279 MatrixType m_matS; 00280 00281 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; 00282 ColumnVectorType m_tmp; 00283 }; 00284 00285 //template<typename MatrixType> 00286 //typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const 00287 //{ 00288 // eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00289 // eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 00290 // Index n = m_eivec.cols(); 00291 // EigenvectorsType matV(n,n); 00292 // // TODO 00293 // return matV; 00294 //} 00295 00296 template<typename MatrixType> 00297 GeneralizedEigenSolver<MatrixType>& 00298 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) 00299 { 00300 check_template_parameters(); 00301 00302 using std::sqrt; 00303 using std::abs; 00304 eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows()); 00305 00306 // Reduce to generalized real Schur form: 00307 // A = Q S Z and B = Q T Z 00308 m_realQZ.compute(A, B, computeEigenvectors); 00309 00310 if (m_realQZ.info() == Success) 00311 { 00312 m_matS = m_realQZ.matrixS(); 00313 if (computeEigenvectors) 00314 m_eivec = m_realQZ.matrixZ().transpose(); 00315 00316 // Compute eigenvalues from matS 00317 m_alphas.resize(A.cols()); 00318 m_betas.resize(A.cols()); 00319 Index i = 0; 00320 while (i < A.cols()) 00321 { 00322 if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0)) 00323 { 00324 m_alphas.coeffRef(i) = m_matS.coeff(i, i); 00325 m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); 00326 ++i; 00327 } 00328 else 00329 { 00330 Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1)); 00331 Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1))); 00332 m_alphas.coeffRef(i) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z); 00333 m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z); 00334 00335 m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); 00336 m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i); 00337 i += 2; 00338 } 00339 } 00340 } 00341 00342 m_isInitialized = true; 00343 m_eigenvectorsOk = false;//computeEigenvectors; 00344 00345 return *this; 00346 } 00347 00348 } // end namespace Eigen 00349 00350 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H