Eigen  3.3.3
GeneralizedEigenSolver.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
00006 // Copyright (C) 2016 Tobias Wood <tobias@spinicist.org.uk>
00007 //
00008 // This Source Code Form is subject to the terms of the Mozilla
00009 // Public License v. 2.0. If a copy of the MPL was not distributed
00010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00011 
00012 #ifndef EIGEN_GENERALIZEDEIGENSOLVER_H
00013 #define EIGEN_GENERALIZEDEIGENSOLVER_H
00014 
00015 #include "./RealQZ.h"
00016 
00017 namespace Eigen { 
00018 
00058 template<typename _MatrixType> class GeneralizedEigenSolver
00059 {
00060   public:
00061 
00063     typedef _MatrixType MatrixType;
00064 
00065     enum {
00066       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00067       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00068       Options = MatrixType::Options,
00069       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00070       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00071     };
00072 
00074     typedef typename MatrixType::Scalar Scalar;
00075     typedef typename NumTraits<Scalar>::Real RealScalar;
00076     typedef Eigen::Index Index; 
00077 
00084     typedef std::complex<RealScalar> ComplexScalar;
00085 
00091     typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType;
00092 
00098     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType;
00099 
00102     typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> EigenvalueType;
00103 
00109     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
00110 
00118     GeneralizedEigenSolver()
00119       : m_eivec(),
00120         m_alphas(),
00121         m_betas(),
00122         m_valuesOkay(false),
00123         m_vectorsOkay(false),
00124         m_realQZ()
00125     {}
00126 
00133     explicit GeneralizedEigenSolver(Index size)
00134       : m_eivec(size, size),
00135         m_alphas(size),
00136         m_betas(size),
00137         m_valuesOkay(false),
00138         m_vectorsOkay(false),
00139         m_realQZ(size),
00140         m_tmp(size)
00141     {}
00142 
00155     GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true)
00156       : m_eivec(A.rows(), A.cols()),
00157         m_alphas(A.cols()),
00158         m_betas(A.cols()),
00159         m_valuesOkay(false),
00160         m_vectorsOkay(false),
00161         m_realQZ(A.cols()),
00162         m_tmp(A.cols())
00163     {
00164       compute(A, B, computeEigenvectors);
00165     }
00166 
00167     /* \brief Returns the computed generalized eigenvectors.
00168       *
00169       * \returns  %Matrix whose columns are the (possibly complex) right eigenvectors.
00170       * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues.
00171       *
00172       * \pre Either the constructor 
00173       * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
00174       * compute(const MatrixType&, const MatrixType& bool) has been called before, and
00175       * \p computeEigenvectors was set to true (the default).
00176       *
00177       * \sa eigenvalues()
00178       */
00179     EigenvectorsType eigenvectors() const {
00180       eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.");
00181       return m_eivec;
00182     }
00183 
00202     EigenvalueType eigenvalues() const
00203     {
00204       eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
00205       return EigenvalueType(m_alphas,m_betas);
00206     }
00207 
00213     ComplexVectorType alphas() const
00214     {
00215       eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
00216       return m_alphas;
00217     }
00218 
00224     VectorType betas() const
00225     {
00226       eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
00227       return m_betas;
00228     }
00229 
00253     GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true);
00254 
00255     ComputationInfo info() const
00256     {
00257       eigen_assert(m_valuesOkay && "EigenSolver is not initialized.");
00258       return m_realQZ.info();
00259     }
00260 
00263     GeneralizedEigenSolver& setMaxIterations(Index maxIters)
00264     {
00265       m_realQZ.setMaxIterations(maxIters);
00266       return *this;
00267     }
00268 
00269   protected:
00270     
00271     static void check_template_parameters()
00272     {
00273       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00274       EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
00275     }
00276     
00277     EigenvectorsType m_eivec;
00278     ComplexVectorType m_alphas;
00279     VectorType m_betas;
00280     bool m_valuesOkay, m_vectorsOkay;
00281     RealQZ<MatrixType> m_realQZ;
00282     ComplexVectorType m_tmp;
00283 };
00284 
00285 template<typename MatrixType>
00286 GeneralizedEigenSolver<MatrixType>&
00287 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
00288 {
00289   check_template_parameters();
00290   
00291   using std::sqrt;
00292   using std::abs;
00293   eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
00294   Index size = A.cols();
00295   m_valuesOkay = false;
00296   m_vectorsOkay = false;
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   if (m_realQZ.info() == Success)
00301   {
00302     // Resize storage
00303     m_alphas.resize(size);
00304     m_betas.resize(size);
00305     if (computeEigenvectors)
00306     {
00307       m_eivec.resize(size,size);
00308       m_tmp.resize(size);
00309     }
00310 
00311     // Aliases:
00312     Map<VectorType> v(reinterpret_cast<Scalar*>(m_tmp.data()), size);
00313     ComplexVectorType &cv = m_tmp;
00314     const MatrixType &mZ = m_realQZ.matrixZ();
00315     const MatrixType &mS = m_realQZ.matrixS();
00316     const MatrixType &mT = m_realQZ.matrixT();
00317 
00318     Index i = 0;
00319     while (i < size)
00320     {
00321       if (i == size - 1 || mS.coeff(i+1, i) == Scalar(0))
00322       {
00323         // Real eigenvalue
00324         m_alphas.coeffRef(i) = mS.diagonal().coeff(i);
00325         m_betas.coeffRef(i)  = mT.diagonal().coeff(i);
00326         if (computeEigenvectors)
00327         {
00328           v.setConstant(Scalar(0.0));
00329           v.coeffRef(i) = Scalar(1.0);
00330           // For singular eigenvalues do nothing more
00331           if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
00332           {
00333             // Non-singular eigenvalue
00334             const Scalar alpha = real(m_alphas.coeffRef(i));
00335             const Scalar beta = m_betas.coeffRef(i);
00336             for (Index j = i-1; j >= 0; j--)
00337             {
00338               const Index st = j+1;
00339               const Index sz = i-j;
00340               if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
00341               {
00342                 // 2x2 block
00343                 Matrix<Scalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( v.segment(st,sz) );
00344                 Matrix<Scalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1);
00345                 v.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
00346                 j--;
00347               }
00348               else
00349               {
00350                 v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j));
00351               }
00352             }
00353           }
00354           m_eivec.col(i).real().noalias() = mZ.transpose() * v;
00355           m_eivec.col(i).real().normalize();
00356           m_eivec.col(i).imag().setConstant(0);
00357         }
00358         ++i;
00359       }
00360       else
00361       {
00362         // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
00363         // Then taking beta=T_00*T_11, we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00):
00364 
00365         // T =  [a 0]
00366         //      [0 b]
00367         RealScalar a = mT.diagonal().coeff(i),
00368                    b = mT.diagonal().coeff(i+1);
00369         const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
00370 
00371         // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
00372         Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
00373 
00374         Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
00375         Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
00376         const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z);
00377         m_alphas.coeffRef(i)   = conj(alpha);
00378         m_alphas.coeffRef(i+1) = alpha;
00379 
00380         if (computeEigenvectors) {
00381           // Compute eigenvector in position (i+1) and then position (i) is just the conjugate
00382           cv.setZero();
00383           cv.coeffRef(i+1) = Scalar(1.0);
00384           // here, the "static_cast" workaound expression template issues.
00385           cv.coeffRef(i) = -(static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1))
00386                           / (static_cast<Scalar>(beta*mS.coeffRef(i,i))   - alpha*mT.coeffRef(i,i));
00387           for (Index j = i-1; j >= 0; j--)
00388           {
00389             const Index st = j+1;
00390             const Index sz = i+1-j;
00391             if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
00392             {
00393               // 2x2 block
00394               Matrix<ComplexScalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( cv.segment(st,sz) );
00395               Matrix<ComplexScalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1);
00396               cv.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
00397               j--;
00398             } else {
00399               cv.coeffRef(j) =  cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum()
00400                               / (alpha*mT.coeffRef(j,j) - static_cast<Scalar>(beta*mS.coeffRef(j,j)));
00401             }
00402           }
00403           m_eivec.col(i+1).noalias() = (mZ.transpose() * cv);
00404           m_eivec.col(i+1).normalize();
00405           m_eivec.col(i) = m_eivec.col(i+1).conjugate();
00406         }
00407         i += 2;
00408       }
00409     }
00410 
00411     m_valuesOkay = true;
00412     m_vectorsOkay = computeEigenvectors;
00413   }
00414   return *this;
00415 }
00416 
00417 } // end namespace Eigen
00418 
00419 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends