GeneralizedEigenSolver.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) 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     MatrixType m_eivec;
00267     ComplexVectorType m_alphas;
00268     VectorType m_betas;
00269     bool m_isInitialized;
00270     bool m_eigenvectorsOk;
00271     RealQZ<MatrixType> m_realQZ;
00272     MatrixType m_matS;
00273 
00274     typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
00275     ColumnVectorType m_tmp;
00276 };
00277 
00278 //template<typename MatrixType>
00279 //typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const
00280 //{
00281 //  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00282 //  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00283 //  Index n = m_eivec.cols();
00284 //  EigenvectorsType matV(n,n);
00285 //  // TODO
00286 //  return matV;
00287 //}
00288 
00289 template<typename MatrixType>
00290 GeneralizedEigenSolver<MatrixType>&
00291 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
00292 {
00293   using std::sqrt;
00294   using std::abs;
00295   eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
00296 
00297   // Reduce to generalized real Schur form:
00298   // A = Q S Z and B = Q T Z
00299   m_realQZ.compute(A, B, computeEigenvectors);
00300 
00301   if (m_realQZ.info() == Success)
00302   {
00303     m_matS = m_realQZ.matrixS();
00304     if (computeEigenvectors)
00305       m_eivec = m_realQZ.matrixZ().transpose();
00306   
00307     // Compute eigenvalues from matS
00308     m_alphas.resize(A.cols());
00309     m_betas.resize(A.cols());
00310     Index i = 0;
00311     while (i < A.cols())
00312     {
00313       if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0))
00314       {
00315         m_alphas.coeffRef(i) = m_matS.coeff(i, i);
00316         m_betas.coeffRef(i)  = m_realQZ.matrixT().coeff(i,i);
00317         ++i;
00318       }
00319       else
00320       {
00321         Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1));
00322         Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1)));
00323         m_alphas.coeffRef(i)   = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z);
00324         m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z);
00325 
00326         m_betas.coeffRef(i)   = m_realQZ.matrixT().coeff(i,i);
00327         m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i);
00328         i += 2;
00329       }
00330     }
00331   }
00332 
00333   m_isInitialized = true;
00334   m_eigenvectorsOk = false;//computeEigenvectors;
00335 
00336   return *this;
00337 }
00338 
00339 } // end namespace Eigen
00340 
00341 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:37:11