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