Lines 89-95
Link Here
|
89 |
*/ |
89 |
*/ |
90 |
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType; |
90 |
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType; |
91 |
|
91 |
|
92 |
/** \brief Type for vector of complex scalar values eigenvalues as returned by betas(). |
92 |
/** \brief Type for vector of complex scalar values eigenvalues as returned by alphas(). |
93 |
* |
93 |
* |
94 |
* This is a column vector with entries of type #ComplexScalar. |
94 |
* This is a column vector with entries of type #ComplexScalar. |
95 |
* The length of the vector is the size of #MatrixType. |
95 |
* The length of the vector is the size of #MatrixType. |
Lines 114-120
Link Here
|
114 |
* |
114 |
* |
115 |
* \sa compute() for an example. |
115 |
* \sa compute() for an example. |
116 |
*/ |
116 |
*/ |
117 |
GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {} |
117 |
GeneralizedEigenSolver() |
|
|
118 |
: m_eivec(), |
119 |
m_alphas(), |
120 |
m_betas(), |
121 |
m_valuesOkay(false), |
122 |
m_vectorsOkay(false), |
123 |
m_realQZ() |
124 |
{} |
118 |
|
125 |
|
119 |
/** \brief Default constructor with memory preallocation |
126 |
/** \brief Default constructor with memory preallocation |
120 |
* |
127 |
* |
Lines 126-136
Link Here
|
126 |
: m_eivec(size, size), |
133 |
: m_eivec(size, size), |
127 |
m_alphas(size), |
134 |
m_alphas(size), |
128 |
m_betas(size), |
135 |
m_betas(size), |
129 |
m_isInitialized(false), |
136 |
m_valuesOkay(false), |
130 |
m_eigenvectorsOk(false), |
137 |
m_vectorsOkay(false), |
131 |
m_realQZ(size), |
138 |
m_realQZ(size) |
132 |
m_matS(size, size), |
|
|
133 |
m_tmp(size) |
134 |
{} |
139 |
{} |
135 |
|
140 |
|
136 |
/** \brief Constructor; computes the generalized eigendecomposition of given matrix pair. |
141 |
/** \brief Constructor; computes the generalized eigendecomposition of given matrix pair. |
Lines 149-181
Link Here
|
149 |
: m_eivec(A.rows(), A.cols()), |
154 |
: m_eivec(A.rows(), A.cols()), |
150 |
m_alphas(A.cols()), |
155 |
m_alphas(A.cols()), |
151 |
m_betas(A.cols()), |
156 |
m_betas(A.cols()), |
152 |
m_isInitialized(false), |
157 |
m_valuesOkay(false), |
153 |
m_eigenvectorsOk(false), |
158 |
m_vectorsOkay(false), |
154 |
m_realQZ(A.cols()), |
159 |
m_realQZ(A.cols()) |
155 |
m_matS(A.rows(), A.cols()), |
|
|
156 |
m_tmp(A.cols()) |
157 |
{ |
160 |
{ |
158 |
compute(A, B, computeEigenvectors); |
161 |
compute(A, B, computeEigenvectors); |
159 |
} |
162 |
} |
160 |
|
163 |
|
161 |
/* \brief Returns the computed generalized eigenvectors. |
164 |
/* \brief Returns the computed generalized eigenvectors. |
162 |
* |
165 |
* |
163 |
* \returns %Matrix whose columns are the (possibly complex) eigenvectors. |
166 |
* \returns %Matrix whose columns are the (possibly complex) right eigenvectors. |
|
|
167 |
* i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues. |
164 |
* |
168 |
* |
165 |
* \pre Either the constructor |
169 |
* \pre Either the constructor |
166 |
* GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function |
170 |
* GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function |
167 |
* compute(const MatrixType&, const MatrixType& bool) has been called before, and |
171 |
* compute(const MatrixType&, const MatrixType& bool) has been called before, and |
168 |
* \p computeEigenvectors was set to true (the default). |
172 |
* \p computeEigenvectors was set to true (the default). |
169 |
* |
173 |
* |
170 |
* Column \f$ k \f$ of the returned matrix is an eigenvector corresponding |
|
|
171 |
* to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The |
172 |
* eigenvectors are normalized to have (Euclidean) norm equal to one. The |
173 |
* matrix returned by this function is the matrix \f$ V \f$ in the |
174 |
* generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists. |
175 |
* |
176 |
* \sa eigenvalues() |
174 |
* \sa eigenvalues() |
177 |
*/ |
175 |
*/ |
178 |
// EigenvectorsType eigenvectors() const; |
176 |
EigenvectorsType eigenvectors() const { |
|
|
177 |
eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.") |
178 |
return m_eivec; |
179 |
} |
179 |
|
180 |
|
180 |
/** \brief Returns an expression of the computed generalized eigenvalues. |
181 |
/** \brief Returns an expression of the computed generalized eigenvalues. |
181 |
* |
182 |
* |
Lines 197-203
Link Here
|
197 |
*/ |
198 |
*/ |
198 |
EigenvalueType eigenvalues() const |
199 |
EigenvalueType eigenvalues() const |
199 |
{ |
200 |
{ |
200 |
eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); |
201 |
eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); |
201 |
return EigenvalueType(m_alphas,m_betas); |
202 |
return EigenvalueType(m_alphas,m_betas); |
202 |
} |
203 |
} |
203 |
|
204 |
|
Lines 208-214
Link Here
|
208 |
* \sa betas(), eigenvalues() */ |
209 |
* \sa betas(), eigenvalues() */ |
209 |
ComplexVectorType alphas() const |
210 |
ComplexVectorType alphas() const |
210 |
{ |
211 |
{ |
211 |
eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); |
212 |
eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); |
212 |
return m_alphas; |
213 |
return m_alphas; |
213 |
} |
214 |
} |
214 |
|
215 |
|
Lines 219-225
Link Here
|
219 |
* \sa alphas(), eigenvalues() */ |
220 |
* \sa alphas(), eigenvalues() */ |
220 |
VectorType betas() const |
221 |
VectorType betas() const |
221 |
{ |
222 |
{ |
222 |
eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); |
223 |
eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); |
223 |
return m_betas; |
224 |
return m_betas; |
224 |
} |
225 |
} |
225 |
|
226 |
|
Lines 250-256
Link Here
|
250 |
|
251 |
|
251 |
ComputationInfo info() const |
252 |
ComputationInfo info() const |
252 |
{ |
253 |
{ |
253 |
eigen_assert(m_isInitialized && "EigenSolver is not initialized."); |
254 |
eigen_assert(m_valuesOkay && "EigenSolver is not initialized."); |
254 |
return m_realQZ.info(); |
255 |
return m_realQZ.info(); |
255 |
} |
256 |
} |
256 |
|
257 |
|
Lines 270-298
Link Here
|
270 |
EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); |
271 |
EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); |
271 |
} |
272 |
} |
272 |
|
273 |
|
273 |
MatrixType m_eivec; |
274 |
EigenvectorsType m_eivec; |
274 |
ComplexVectorType m_alphas; |
275 |
ComplexVectorType m_alphas; |
275 |
VectorType m_betas; |
276 |
VectorType m_betas; |
276 |
bool m_isInitialized; |
277 |
bool m_valuesOkay, m_vectorsOkay; |
277 |
bool m_eigenvectorsOk; |
|
|
278 |
RealQZ<MatrixType> m_realQZ; |
278 |
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 |
}; |
279 |
}; |
284 |
|
280 |
|
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> |
281 |
template<typename MatrixType> |
297 |
GeneralizedEigenSolver<MatrixType>& |
282 |
GeneralizedEigenSolver<MatrixType>& |
298 |
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) |
283 |
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) |
Lines 302-359
Link Here
|
302 |
using std::sqrt; |
287 |
using std::sqrt; |
303 |
using std::abs; |
288 |
using std::abs; |
304 |
eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows()); |
289 |
eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows()); |
305 |
|
290 |
m_valuesOkay = false; |
|
|
291 |
m_vectorsOkay = false; |
306 |
// Reduce to generalized real Schur form: |
292 |
// Reduce to generalized real Schur form: |
307 |
// A = Q S Z and B = Q T Z |
293 |
// A = Q S Z and B = Q T Z |
308 |
m_realQZ.compute(A, B, computeEigenvectors); |
294 |
m_realQZ.compute(A, B, computeEigenvectors); |
309 |
|
295 |
if (m_realQZ.info() == Success) { |
310 |
if (m_realQZ.info() == Success) |
296 |
// Temp space for the untransformed eigenvectors |
311 |
{ |
297 |
VectorType v; |
312 |
m_matS = m_realQZ.matrixS(); |
298 |
ComplexVectorType cv; |
313 |
const MatrixType &matT = m_realQZ.matrixT(); |
299 |
// Resize storage |
314 |
if (computeEigenvectors) |
|
|
315 |
m_eivec = m_realQZ.matrixZ().transpose(); |
316 |
|
317 |
// Compute eigenvalues from matS |
318 |
m_alphas.resize(A.cols()); |
300 |
m_alphas.resize(A.cols()); |
319 |
m_betas.resize(A.cols()); |
301 |
m_betas.resize(A.cols()); |
|
|
302 |
if (computeEigenvectors) { |
303 |
m_eivec.resize(A.cols(), A.cols()); |
304 |
v.resize(A.cols()); |
305 |
cv.resize(A.cols()); |
306 |
} |
307 |
// Grab some references |
308 |
const MatrixType &mZT = m_realQZ.matrixZ().transpose(); |
309 |
const MatrixType &mS = m_realQZ.matrixS(); |
310 |
const MatrixType &mT = m_realQZ.matrixT(); |
320 |
Index i = 0; |
311 |
Index i = 0; |
321 |
while (i < A.cols()) |
312 |
while (i < A.cols()) { |
322 |
{ |
313 |
if (i == A.cols() - 1 || mS.coeff(i+1, i) == Scalar(0)) { |
323 |
if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0)) |
314 |
// Real eigenvalue |
324 |
{ |
315 |
m_alphas.coeffRef(i) = mS.diagonal().coeff(i); |
325 |
m_alphas.coeffRef(i) = m_matS.coeff(i, i); |
316 |
m_betas.coeffRef(i) = mT.diagonal().coeff(i); |
326 |
m_betas.coeffRef(i) = matT.coeff(i,i); |
317 |
if (computeEigenvectors) { |
|
|
318 |
v.setConstant(Scalar(0.0)); |
319 |
v.coeffRef(i) = Scalar(1.0); |
320 |
// For singular eigenvalues do nothing more |
321 |
if (!internal::isApprox(fabs(m_betas.coeffRef(i)), 0.0)) { |
322 |
//Non-singular eigenvalue |
323 |
const Scalar alpha = real(m_alphas.coeffRef(i)); |
324 |
const Scalar beta = m_betas.coeffRef(i); |
325 |
for (Index j = i-1; j >= 0; j--) { |
326 |
const Index st = j+1; |
327 |
const Index sz = i-j; |
328 |
if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) { // 2x2 block |
329 |
Matrix<Scalar, 2, 1> RHS; |
330 |
RHS[0] = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j-1,st,1,sz) - alpha*mT.block(j-1,st,1,sz)).sum(); |
331 |
RHS[1] = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum(); |
332 |
Matrix<Scalar, 2, 2> LHS; |
333 |
LHS << beta*mS.coeffRef(j-1,j-1) - alpha*mT.coeffRef(j-1,j-1), beta*mS.coeffRef(j-1,j) - alpha*mT.coeffRef(j-1,j), |
334 |
beta*mS.coeffRef(j,j-1) - alpha*mT.coeffRef(j,j-1), beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j); |
335 |
v.segment(j-1,2) = LHS.partialPivLu().solve(RHS); |
336 |
j--; |
337 |
} else { |
338 |
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)); |
339 |
} |
340 |
} |
341 |
} |
342 |
m_eivec.col(i).real() = (mZT * v).normalized(); |
343 |
m_eivec.col(i).imag().setConstant(0); |
344 |
} |
327 |
++i; |
345 |
++i; |
328 |
} |
346 |
} else { |
329 |
else |
|
|
330 |
{ |
331 |
// We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T |
347 |
// We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T |
332 |
// 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): |
348 |
// 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): |
333 |
|
349 |
|
334 |
// T = [a 0] |
350 |
// T = [a 0] |
335 |
// [0 b] |
351 |
// [0 b] |
336 |
RealScalar a = matT.diagonal().coeff(i), |
352 |
RealScalar a = mT.diagonal().coeff(i), |
337 |
b = matT.diagonal().coeff(i+1); |
353 |
b = mT.diagonal().coeff(i+1); |
|
|
354 |
const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b; |
355 |
|
338 |
// ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug. |
356 |
// ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug. |
339 |
Matrix<RealScalar,2,2> S2 = m_matS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal(); |
357 |
Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal(); |
340 |
|
358 |
|
341 |
Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1)); |
359 |
Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1)); |
342 |
Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1))); |
360 |
Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1))); |
343 |
m_alphas.coeffRef(i) = ComplexScalar(S2.coeff(1,1) + p, z); |
361 |
const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z); |
344 |
m_alphas.coeffRef(i+1) = ComplexScalar(S2.coeff(1,1) + p, -z); |
362 |
m_alphas.coeffRef(i) = conj(alpha); |
|
|
363 |
m_alphas.coeffRef(i+1) = alpha; |
345 |
|
364 |
|
346 |
m_betas.coeffRef(i) = |
365 |
if (computeEigenvectors) { |
347 |
m_betas.coeffRef(i+1) = a*b; |
366 |
// Compute eigenvector in position (i+1) and then position (i) is just the conjugate |
348 |
|
367 |
cv.setConstant(ComplexScalar(0.0, 0.0)); |
|
|
368 |
cv.coeffRef(i+1) = ComplexScalar(1.0, 0.0); |
369 |
cv.coeffRef(i) = -(beta*mS.coeffRef(i,i+1) - alpha*mT.coeffRef(i,i+1)) / (beta*mS.coeffRef(i,i) - alpha*mT.coeffRef(i,i)); |
370 |
for (Index j = i-1; j >= 0; j--) { |
371 |
const Index st = j+1; |
372 |
const Index sz = i+1-j; |
373 |
if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) { // 2x2 block |
374 |
Matrix<ComplexScalar, 2, 1> RHS; |
375 |
RHS[0] = -cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j-1,st,1,sz) - alpha*mT.block(j-1,st,1,sz)).sum(); |
376 |
RHS[1] = -cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum(); |
377 |
Matrix<ComplexScalar, 2, 2> LHS; |
378 |
LHS << beta*mS.coeffRef(j-1,j-1) - alpha*mT.coeffRef(j-1,j-1), beta*mS.coeffRef(j-1,j) - alpha*mT.coeffRef(j-1,j), |
379 |
beta*mS.coeffRef(j,j-1) - alpha*mT.coeffRef(j,j-1), beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j); |
380 |
cv.segment(j-1,2) = LHS.partialPivLu().solve(RHS); |
381 |
j--; |
382 |
} else { |
383 |
cv.coeffRef(j) = -cv.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)); |
384 |
} |
385 |
} |
386 |
m_eivec.col(i+1) = (mZT * cv).normalized(); |
387 |
m_eivec.col(i) = m_eivec.col(i+1).conjugate(); |
388 |
} |
349 |
i += 2; |
389 |
i += 2; |
350 |
} |
390 |
} |
351 |
} |
391 |
} |
|
|
392 |
m_valuesOkay = true; |
393 |
m_vectorsOkay = computeEigenvectors; |
352 |
} |
394 |
} |
353 |
|
|
|
354 |
m_isInitialized = true; |
355 |
m_eigenvectorsOk = false;//computeEigenvectors; |
356 |
|
357 |
return *this; |
395 |
return *this; |
358 |
} |
396 |
} |
359 |
|
397 |
|