Lines 12-17
Link Here
|
12 |
#define EIGEN_GENERALIZEDEIGENSOLVER_H |
12 |
#define EIGEN_GENERALIZEDEIGENSOLVER_H |
13 |
|
13 |
|
14 |
#include "./RealQZ.h" |
14 |
#include "./RealQZ.h" |
|
|
15 |
#include <iostream> |
15 |
|
16 |
|
16 |
namespace Eigen { |
17 |
namespace Eigen { |
17 |
|
18 |
|
Lines 114-120
Link Here
|
114 |
* |
115 |
* |
115 |
* \sa compute() for an example. |
116 |
* \sa compute() for an example. |
116 |
*/ |
117 |
*/ |
117 |
GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {} |
118 |
GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ() {} |
118 |
|
119 |
|
119 |
/** \brief Default constructor with memory preallocation |
120 |
/** \brief Default constructor with memory preallocation |
120 |
* |
121 |
* |
Lines 127-136
Link Here
|
127 |
m_alphas(size), |
128 |
m_alphas(size), |
128 |
m_betas(size), |
129 |
m_betas(size), |
129 |
m_isInitialized(false), |
130 |
m_isInitialized(false), |
130 |
m_eigenvectorsOk(false), |
131 |
m_realQZ(size) |
131 |
m_realQZ(size), |
|
|
132 |
m_matS(size, size), |
133 |
m_tmp(size) |
134 |
{} |
132 |
{} |
135 |
|
133 |
|
136 |
/** \brief Constructor; computes the generalized eigendecomposition of given matrix pair. |
134 |
/** \brief Constructor; computes the generalized eigendecomposition of given matrix pair. |
Lines 150-159
Link Here
|
150 |
m_alphas(A.cols()), |
148 |
m_alphas(A.cols()), |
151 |
m_betas(A.cols()), |
149 |
m_betas(A.cols()), |
152 |
m_isInitialized(false), |
150 |
m_isInitialized(false), |
153 |
m_eigenvectorsOk(false), |
151 |
m_realQZ(A.cols()) |
154 |
m_realQZ(A.cols()), |
|
|
155 |
m_matS(A.rows(), A.cols()), |
156 |
m_tmp(A.cols()) |
157 |
{ |
152 |
{ |
158 |
compute(A, B, computeEigenvectors); |
153 |
compute(A, B, computeEigenvectors); |
159 |
} |
154 |
} |
Lines 201-206
Link Here
|
201 |
return EigenvalueType(m_alphas,m_betas); |
196 |
return EigenvalueType(m_alphas,m_betas); |
202 |
} |
197 |
} |
203 |
|
198 |
|
|
|
199 |
/** \brief Returns the computed generalized eigenvectors. |
200 |
* |
201 |
* \returns Returns a matrix whose columns are the (possibly complex) right eigenvectors |
202 |
* i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues. |
203 |
* |
204 |
* \pre Either the constructor |
205 |
* GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function |
206 |
* compute(const MatrixType&,const MatrixType&,bool) has been called before. |
207 |
* |
208 |
* \sa eigenvalues() |
209 |
*/ |
210 |
EigenvectorsType eigenvectors() const |
211 |
{ |
212 |
eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); |
213 |
return m_eivec; |
214 |
} |
215 |
|
204 |
/** \returns A const reference to the vectors containing the alpha values |
216 |
/** \returns A const reference to the vectors containing the alpha values |
205 |
* |
217 |
* |
206 |
* This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). |
218 |
* This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). |
Lines 270-298
Link Here
|
270 |
EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); |
282 |
EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); |
271 |
} |
283 |
} |
272 |
|
284 |
|
273 |
MatrixType m_eivec; |
285 |
EigenvectorsType m_eivec; |
274 |
ComplexVectorType m_alphas; |
286 |
ComplexVectorType m_alphas; |
275 |
VectorType m_betas; |
287 |
VectorType m_betas; |
276 |
bool m_isInitialized; |
288 |
bool m_isInitialized; |
277 |
bool m_eigenvectorsOk; |
|
|
278 |
RealQZ<MatrixType> m_realQZ; |
289 |
RealQZ<MatrixType> m_realQZ; |
279 |
MatrixType m_matS; |
|
|
280 |
|
281 |
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; |
282 |
ColumnVectorType m_tmp; |
283 |
}; |
290 |
}; |
284 |
|
291 |
|
285 |
//template<typename MatrixType> |
|
|
286 |
//typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const |
287 |
//{ |
288 |
// eigen_assert(m_isInitialized && "EigenSolver is not initialized."); |
289 |
// eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); |
290 |
// Index n = m_eivec.cols(); |
291 |
// EigenvectorsType matV(n,n); |
292 |
// // TODO |
293 |
// return matV; |
294 |
//} |
295 |
|
296 |
template<typename MatrixType> |
292 |
template<typename MatrixType> |
297 |
GeneralizedEigenSolver<MatrixType>& |
293 |
GeneralizedEigenSolver<MatrixType>& |
298 |
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) |
294 |
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) |
Lines 307-347
Link Here
|
307 |
// A = Q S Z and B = Q T Z |
303 |
// A = Q S Z and B = Q T Z |
308 |
m_realQZ.compute(A, B, computeEigenvectors); |
304 |
m_realQZ.compute(A, B, computeEigenvectors); |
309 |
|
305 |
|
310 |
if (m_realQZ.info() == Success) |
306 |
if (m_realQZ.info() == Success) { |
311 |
{ |
307 |
// Resize storage |
312 |
m_matS = m_realQZ.matrixS(); |
|
|
313 |
if (computeEigenvectors) |
314 |
m_eivec = m_realQZ.matrixZ().transpose(); |
315 |
|
316 |
// Compute eigenvalues from matS |
317 |
m_alphas.resize(A.cols()); |
308 |
m_alphas.resize(A.cols()); |
318 |
m_betas.resize(A.cols()); |
309 |
m_betas.resize(A.cols()); |
|
|
310 |
if (computeEigenvectors) { |
311 |
m_eivec.resize(A.cols(), A.cols()); |
312 |
} |
313 |
|
314 |
// Grab some references |
315 |
const MatrixType &mZT = m_realQZ.matrixZ().transpose(); |
316 |
const MatrixType &mS = m_realQZ.matrixS(); |
317 |
const MatrixType &mT = m_realQZ.matrixT(); |
318 |
|
319 |
Index i = 0; |
319 |
Index i = 0; |
320 |
while (i < A.cols()) |
320 |
while (i < A.cols()) { |
321 |
{ |
321 |
if (i == A.cols() - 1 || mS.coeff(i+1, i) == Scalar(0)) { |
322 |
if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0)) |
322 |
// Real eigenvalue |
323 |
{ |
323 |
m_alphas.coeffRef(i) = mS.coeff(i, i); |
324 |
m_alphas.coeffRef(i) = m_matS.coeff(i, i); |
324 |
m_betas.coeffRef(i) = mT.coeff(i,i); |
325 |
m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); |
325 |
if (computeEigenvectors) { |
|
|
326 |
VectorXd v = VectorXd::Zero(A.cols()); |
327 |
v.coeffRef(i) = 1.0; |
328 |
if (fabs(m_betas.coeffRef(i)) < NumTraits<Scalar>::epsilon()) { |
329 |
// Singular eigenvalue, do nothing more |
330 |
} else { |
331 |
Scalar lambda = real(m_alphas.coeffRef(i)) / m_betas.coeffRef(i); |
332 |
MatrixType mC = mS - lambda*mT; |
333 |
for (Index j = i-1; j >= 0; j--) { |
334 |
for (Index k = j+1; k <= i; k++) { |
335 |
v.coeffRef(j) -= mC.coeffRef(j, k) * v.coeffRef(k); |
336 |
} |
337 |
v.coeffRef(j) /= mC.coeffRef(j, j); |
338 |
} |
339 |
} |
340 |
m_eivec.col(i).real() = (mZT * v).normalized(); |
341 |
m_eivec.col(i).imag().setConstant(0.); |
342 |
} |
326 |
++i; |
343 |
++i; |
327 |
} |
344 |
} else { |
328 |
else |
345 |
Scalar p = Scalar(0.5) * (mS.coeffRef(i, i) - mS.coeffRef(i+1, i+1)); |
329 |
{ |
346 |
Scalar z = sqrt(abs(p * p + mS.coeffRef(i+1, i) * mS.coeffRef(i, i+1))); |
330 |
Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1)); |
347 |
m_alphas.coeffRef(i) = ComplexScalar(mS.coeffRef(i+1, i+1) + p, z); |
331 |
Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1))); |
348 |
m_alphas.coeffRef(i+1) = ComplexScalar(mS.coeffRef(i+1, i+1) + p, -z); |
332 |
m_alphas.coeffRef(i) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z); |
|
|
333 |
m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z); |
334 |
|
349 |
|
335 |
m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); |
350 |
m_betas.coeffRef(i) = mT.coeffRef(i,i); |
336 |
m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i); |
351 |
m_betas.coeffRef(i+1) = mT.coeffRef(i,i); |
337 |
i += 2; |
352 |
i += 2; |
338 |
} |
353 |
} |
339 |
} |
354 |
} |
|
|
355 |
m_isInitialized = true; |
340 |
} |
356 |
} |
341 |
|
|
|
342 |
m_isInitialized = true; |
343 |
m_eigenvectorsOk = false;//computeEigenvectors; |
344 |
|
345 |
return *this; |
357 |
return *this; |
346 |
} |
358 |
} |
347 |
|
359 |
|