Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
GeneralizedEigenSolver.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
5// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
6// Copyright (C) 2016 Tobias Wood <tobias@spinicist.org.uk>
7//
8// This Source Code Form is subject to the terms of the Mozilla
9// Public License v. 2.0. If a copy of the MPL was not distributed
10// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11
12#ifndef EIGEN_GENERALIZEDEIGENSOLVER_H
13#define EIGEN_GENERALIZEDEIGENSOLVER_H
14
15#include "./RealQZ.h"
16
17namespace Eigen {
18
58template<typename _MatrixType> class GeneralizedEigenSolver
59{
60 public:
61
63 typedef _MatrixType MatrixType;
64
65 enum {
66 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
67 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
68 Options = MatrixType::Options,
69 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
70 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
71 };
72
74 typedef typename MatrixType::Scalar Scalar;
75 typedef typename NumTraits<Scalar>::Real RealScalar;
77
84 typedef std::complex<RealScalar> ComplexScalar;
85
92
99
103
110
119 : m_eivec(),
120 m_alphas(),
121 m_betas(),
122 m_valuesOkay(false),
123 m_vectorsOkay(false),
124 m_realQZ()
125 {}
126
134 : m_eivec(size, size),
135 m_alphas(size),
136 m_betas(size),
137 m_valuesOkay(false),
138 m_vectorsOkay(false),
139 m_realQZ(size),
140 m_tmp(size)
141 {}
142
155 GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true)
156 : m_eivec(A.rows(), A.cols()),
157 m_alphas(A.cols()),
158 m_betas(A.cols()),
159 m_valuesOkay(false),
160 m_vectorsOkay(false),
161 m_realQZ(A.cols()),
162 m_tmp(A.cols())
163 {
164 compute(A, B, computeEigenvectors);
165 }
166
167 /* \brief Returns the computed generalized eigenvectors.
168 *
169 * \returns %Matrix whose columns are the (possibly complex) right eigenvectors.
170 * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues.
171 *
172 * \pre Either the constructor
173 * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
174 * compute(const MatrixType&, const MatrixType& bool) has been called before, and
175 * \p computeEigenvectors was set to true (the default).
176 *
177 * \sa eigenvalues()
178 */
179 EigenvectorsType eigenvectors() const {
180 eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.");
181 return m_eivec;
182 }
183
203 {
204 eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
205 return EigenvalueType(m_alphas,m_betas);
206 }
207
214 {
215 eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
216 return m_alphas;
217 }
218
225 {
226 eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
227 return m_betas;
228 }
229
253 GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true);
254
255 ComputationInfo info() const
256 {
257 eigen_assert(m_valuesOkay && "EigenSolver is not initialized.");
258 return m_realQZ.info();
259 }
260
264 {
265 m_realQZ.setMaxIterations(maxIters);
266 return *this;
267 }
268
269 protected:
270
271 static void check_template_parameters()
272 {
273 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
274 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
275 }
276
277 EigenvectorsType m_eivec;
278 ComplexVectorType m_alphas;
279 VectorType m_betas;
280 bool m_valuesOkay, m_vectorsOkay;
281 RealQZ<MatrixType> m_realQZ;
282 ComplexVectorType m_tmp;
283};
284
285template<typename MatrixType>
286GeneralizedEigenSolver<MatrixType>&
287GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
288{
289 check_template_parameters();
290
291 using std::sqrt;
292 using std::abs;
293 eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
294 Index size = A.cols();
295 m_valuesOkay = false;
296 m_vectorsOkay = false;
297 // Reduce to generalized real Schur form:
298 // A = Q S Z and B = Q T Z
299 m_realQZ.compute(A, B, computeEigenvectors);
300 if (m_realQZ.info() == Success)
301 {
302 // Resize storage
303 m_alphas.resize(size);
304 m_betas.resize(size);
305 if (computeEigenvectors)
306 {
307 m_eivec.resize(size,size);
308 m_tmp.resize(size);
309 }
310
311 // Aliases:
312 Map<VectorType> v(reinterpret_cast<Scalar*>(m_tmp.data()), size);
313 ComplexVectorType &cv = m_tmp;
314 const MatrixType &mS = m_realQZ.matrixS();
315 const MatrixType &mT = m_realQZ.matrixT();
316
317 Index i = 0;
318 while (i < size)
319 {
320 if (i == size - 1 || mS.coeff(i+1, i) == Scalar(0))
321 {
322 // Real eigenvalue
323 m_alphas.coeffRef(i) = mS.diagonal().coeff(i);
324 m_betas.coeffRef(i) = mT.diagonal().coeff(i);
325 if (computeEigenvectors)
326 {
327 v.setConstant(Scalar(0.0));
328 v.coeffRef(i) = Scalar(1.0);
329 // For singular eigenvalues do nothing more
330 if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
331 {
332 // Non-singular eigenvalue
333 const Scalar alpha = real(m_alphas.coeffRef(i));
334 const Scalar beta = m_betas.coeffRef(i);
335 for (Index j = i-1; j >= 0; j--)
336 {
337 const Index st = j+1;
338 const Index sz = i-j;
339 if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
340 {
341 // 2x2 block
342 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) );
343 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);
344 v.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
345 j--;
346 }
347 else
348 {
349 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));
350 }
351 }
352 }
353 m_eivec.col(i).real().noalias() = m_realQZ.matrixZ().transpose() * v;
354 m_eivec.col(i).real().normalize();
355 m_eivec.col(i).imag().setConstant(0);
356 }
357 ++i;
358 }
359 else
360 {
361 // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
362 // 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):
363
364 // T = [a 0]
365 // [0 b]
366 RealScalar a = mT.diagonal().coeff(i),
367 b = mT.diagonal().coeff(i+1);
368 const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
369
370 // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
371 Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
372
373 Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
374 Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
375 const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z);
376 m_alphas.coeffRef(i) = conj(alpha);
377 m_alphas.coeffRef(i+1) = alpha;
378
379 if (computeEigenvectors) {
380 // Compute eigenvector in position (i+1) and then position (i) is just the conjugate
381 cv.setZero();
382 cv.coeffRef(i+1) = Scalar(1.0);
383 // here, the "static_cast" workaound expression template issues.
384 cv.coeffRef(i) = -(static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1))
385 / (static_cast<Scalar>(beta*mS.coeffRef(i,i)) - alpha*mT.coeffRef(i,i));
386 for (Index j = i-1; j >= 0; j--)
387 {
388 const Index st = j+1;
389 const Index sz = i+1-j;
390 if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
391 {
392 // 2x2 block
393 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) );
394 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);
395 cv.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
396 j--;
397 } else {
398 cv.coeffRef(j) = cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum()
399 / (alpha*mT.coeffRef(j,j) - static_cast<Scalar>(beta*mS.coeffRef(j,j)));
400 }
401 }
402 m_eivec.col(i+1).noalias() = (m_realQZ.matrixZ().transpose() * cv);
403 m_eivec.col(i+1).normalize();
404 m_eivec.col(i) = m_eivec.col(i+1).conjugate();
405 }
406 i += 2;
407 }
408 }
409
410 m_valuesOkay = true;
411 m_vectorsOkay = computeEigenvectors;
412 }
413 return *this;
414}
415
416} // end namespace Eigen
417
418#endif // EIGEN_GENERALIZEDEIGENSOLVER_H
Generic expression where a coefficient-wise binary operator is applied to two expressions.
Definition: CwiseBinaryOp.h:84
Computes the generalized eigenvalues and eigenvectors of a pair of general matrices.
Definition: GeneralizedEigenSolver.h:59
GeneralizedEigenSolver & compute(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Computes generalized eigendecomposition of given matrix.
Definition: GeneralizedEigenSolver.h:287
GeneralizedEigenSolver(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Constructor; computes the generalized eigendecomposition of given matrix pair.
Definition: GeneralizedEigenSolver.h:155
GeneralizedEigenSolver & setMaxIterations(Index maxIters)
Definition: GeneralizedEigenSolver.h:263
Eigen::Index Index
Definition: GeneralizedEigenSolver.h:76
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: GeneralizedEigenSolver.h:63
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > VectorType
Type for vector of real scalar values eigenvalues as returned by betas().
Definition: GeneralizedEigenSolver.h:91
EigenvalueType eigenvalues() const
Returns an expression of the computed generalized eigenvalues.
Definition: GeneralizedEigenSolver.h:202
ComplexVectorType alphas() const
Definition: GeneralizedEigenSolver.h:213
GeneralizedEigenSolver(Index size)
Default constructor with memory preallocation.
Definition: GeneralizedEigenSolver.h:133
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: GeneralizedEigenSolver.h:84
VectorType betas() const
Definition: GeneralizedEigenSolver.h:224
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ComplexVectorType
Type for vector of complex scalar values eigenvalues as returned by alphas().
Definition: GeneralizedEigenSolver.h:98
CwiseBinaryOp< internal::scalar_quotient_op< ComplexScalar, Scalar >, ComplexVectorType, VectorType > EigenvalueType
Expression type for the eigenvalues as returned by eigenvalues().
Definition: GeneralizedEigenSolver.h:102
GeneralizedEigenSolver()
Default constructor.
Definition: GeneralizedEigenSolver.h:118
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: GeneralizedEigenSolver.h:74
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: GeneralizedEigenSolver.h:109
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:96
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:175
const Scalar & coeff(Index rowId, Index colId) const
Definition: PlainObjectBase.h:152
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:562
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealQZ.h:169
RealQZ & setMaxIterations(Index maxIters)
Definition: RealQZ.h:186
ComputationInfo
Definition: Constants.h:440
@ Success
Definition: Constants.h:442
Namespace containing all symbols from the Eigen library.
Definition: Core:141
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:233