Lines 83-88
Link Here
|
83 |
explicit RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) |
83 |
explicit RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) |
84 |
: m_matT(size, size), |
84 |
: m_matT(size, size), |
85 |
m_matU(size, size), |
85 |
m_matU(size, size), |
|
|
86 |
m_eigenvalues(size), |
86 |
m_workspaceVector(size), |
87 |
m_workspaceVector(size), |
87 |
m_hess(size), |
88 |
m_hess(size), |
88 |
m_isInitialized(false), |
89 |
m_isInitialized(false), |
Lines 93-101
Link Here
|
93 |
/** \brief Constructor; computes real Schur decomposition of given matrix. |
94 |
/** \brief Constructor; computes real Schur decomposition of given matrix. |
94 |
* |
95 |
* |
95 |
* \param[in] matrix Square matrix whose Schur decomposition is to be computed. |
96 |
* \param[in] matrix Square matrix whose Schur decomposition is to be computed. |
96 |
* \param[in] computeU If true, both T and U are computed; if false, only T is computed. |
97 |
* \param[in] computeU If true, both T and U are computed; if false, only T and the eigenvalues are computed. |
97 |
* |
98 |
* |
98 |
* This constructor calls compute() to compute the Schur decomposition. |
99 |
* This constructor calls compute() to compute the Schur decomposition and the |
|
|
100 |
* eigenvalues of the matrix. |
99 |
* |
101 |
* |
100 |
* Example: \include RealSchur_RealSchur_MatrixType.cpp |
102 |
* Example: \include RealSchur_RealSchur_MatrixType.cpp |
101 |
* Output: \verbinclude RealSchur_RealSchur_MatrixType.out |
103 |
* Output: \verbinclude RealSchur_RealSchur_MatrixType.out |
Lines 147-156
Link Here
|
147 |
return m_matT; |
149 |
return m_matT; |
148 |
} |
150 |
} |
149 |
|
151 |
|
|
|
152 |
/** \brief Returns the eigenvalues of the given matrix. |
153 |
* |
154 |
* \returns A column vector containing the eigenvalues. |
155 |
* |
156 |
* \pre Either the constructor RealSchur(const MatrixType&, bool) or the |
157 |
* member function compute(const MatrixType&, bool) has been called before |
158 |
* to compute the Schur decomposition of a matrix. |
159 |
* |
160 |
* The eigenvalues of the matrix are the same as the eigenvalues of the 1x1 |
161 |
* and 2x2 diagonal blocks of the matrix T and are extracted from that matrix |
162 |
* after a factorization has been computed. |
163 |
* |
164 |
* \sa RealSchur(const MatrixType&, bool) for an example |
165 |
*/ |
166 |
const EigenvalueType& eigenvalues() const |
167 |
{ |
168 |
eigen_assert(m_isInitialized && "RealSchur is not initialized."); |
169 |
return m_eigenvalues; |
170 |
} |
171 |
|
150 |
/** \brief Computes Schur decomposition of given matrix. |
172 |
/** \brief Computes Schur decomposition of given matrix. |
151 |
* |
173 |
* |
152 |
* \param[in] matrix Square matrix whose Schur decomposition is to be computed. |
174 |
* \param[in] matrix Square matrix whose Schur decomposition is to be computed. |
153 |
* \param[in] computeU If true, both T and U are computed; if false, only T is computed. |
175 |
* \param[in] computeU If true, both T and U are computed; if false, only T and the eigenvalues are computed. |
154 |
* \returns Reference to \c *this |
176 |
* \returns Reference to \c *this |
155 |
* |
177 |
* |
156 |
* The Schur decomposition is computed by first reducing the matrix to |
178 |
* The Schur decomposition is computed by first reducing the matrix to |
Lines 178-185
Link Here
|
178 |
* This routine assumes that the matrix is already reduced in Hessenberg form matrixH |
200 |
* This routine assumes that the matrix is already reduced in Hessenberg form matrixH |
179 |
* using either the class HessenbergDecomposition or another mean. |
201 |
* using either the class HessenbergDecomposition or another mean. |
180 |
* It computes the upper quasi-triangular matrix T of the Schur decomposition of H |
202 |
* It computes the upper quasi-triangular matrix T of the Schur decomposition of H |
|
|
203 |
* and calculates the eigenvalues of T. |
181 |
* When computeU is true, this routine computes the matrix U such that |
204 |
* When computeU is true, this routine computes the matrix U such that |
182 |
* A = U T U^T = (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix |
205 |
* A = U T U^T = (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix. |
183 |
* |
206 |
* |
184 |
* NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix |
207 |
* NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix |
185 |
* is not available, the user should give an identity matrix (Q.setIdentity()) |
208 |
* is not available, the user should give an identity matrix (Q.setIdentity()) |
Lines 226-231
Link Here
|
226 |
|
249 |
|
227 |
MatrixType m_matT; |
250 |
MatrixType m_matT; |
228 |
MatrixType m_matU; |
251 |
MatrixType m_matU; |
|
|
252 |
EigenvalueType m_eigenvalues; |
229 |
ColumnVectorType m_workspaceVector; |
253 |
ColumnVectorType m_workspaceVector; |
230 |
HessenbergDecomposition<MatrixType> m_hess; |
254 |
HessenbergDecomposition<MatrixType> m_hess; |
231 |
ComputationInfo m_info; |
255 |
ComputationInfo m_info; |
Lines 241-246
Link Here
|
241 |
void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo); |
265 |
void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo); |
242 |
void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector); |
266 |
void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector); |
243 |
void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace); |
267 |
void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace); |
|
|
268 |
/** \brief Resizes vector m_eigenvalues and populates it with the |
269 |
* eigenvalues of quasi-triangular matrix m_T. |
270 |
* \return NumericalIssue if a non-finite eigenvalue was detected, Success otherwise |
271 |
*/ |
272 |
ComputationInfo extractEigenvalues(); |
244 |
}; |
273 |
}; |
245 |
|
274 |
|
246 |
|
275 |
|
Lines 259-264
Link Here
|
259 |
if(scale<considerAsZero) |
288 |
if(scale<considerAsZero) |
260 |
{ |
289 |
{ |
261 |
m_matT.setZero(matrix.rows(),matrix.cols()); |
290 |
m_matT.setZero(matrix.rows(),matrix.cols()); |
|
|
291 |
m_eigenvalues.setZero(matrix.rows()); |
262 |
if(computeU) |
292 |
if(computeU) |
263 |
m_matU.setIdentity(matrix.rows(),matrix.cols()); |
293 |
m_matU.setIdentity(matrix.rows(),matrix.cols()); |
264 |
m_info = Success; |
294 |
m_info = Success; |
Lines 274-282
Link Here
|
274 |
computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU); |
304 |
computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU); |
275 |
|
305 |
|
276 |
m_matT *= scale; |
306 |
m_matT *= scale; |
277 |
|
307 |
m_eigenvalues *= scale; |
|
|
308 |
|
278 |
return *this; |
309 |
return *this; |
279 |
} |
310 |
} |
|
|
311 |
|
280 |
template<typename MatrixType> |
312 |
template<typename MatrixType> |
281 |
template<typename HessMatrixType, typename OrthMatrixType> |
313 |
template<typename HessMatrixType, typename OrthMatrixType> |
282 |
RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) |
314 |
RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) |
Lines 343-348
Link Here
|
343 |
else |
375 |
else |
344 |
m_info = NoConvergence; |
376 |
m_info = NoConvergence; |
345 |
|
377 |
|
|
|
378 |
ComputationInfo ev_info = extractEigenvalues(); |
379 |
// If the factorization failed, we keep the original message (NoConvergence...) |
380 |
// Otherwise, set m_info to the result of extractEigenvalues(), |
381 |
// which may return NumericalIssue if a NaN/infinity is found. |
382 |
if (m_info == Success) |
383 |
m_info = ev_info; |
384 |
|
346 |
m_isInitialized = true; |
385 |
m_isInitialized = true; |
347 |
m_matUisUptodate = computeU; |
386 |
m_matUisUptodate = computeU; |
348 |
return *this; |
387 |
return *this; |
Lines 541-546
Link Here
|
541 |
} |
580 |
} |
542 |
} |
581 |
} |
543 |
|
582 |
|
|
|
583 |
template<typename MatrixType> |
584 |
ComputationInfo RealSchur<MatrixType>::extractEigenvalues() |
585 |
{ |
586 |
using std::sqrt; |
587 |
using std::abs; |
588 |
using numext::isfinite; |
589 |
m_eigenvalues.setZero(m_matT.cols()); |
590 |
Index i = 0; |
591 |
while (i < m_matT.cols()) |
592 |
{ |
593 |
if (i == m_matT.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) |
594 |
{ |
595 |
m_eigenvalues.coeffRef(i) = m_matT.coeff(i, i); |
596 |
if(!(isfinite)(m_eigenvalues.coeffRef(i))) |
597 |
{ |
598 |
return NumericalIssue; |
599 |
} |
600 |
++i; |
601 |
} |
602 |
else |
603 |
{ |
604 |
Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1)); |
605 |
Scalar z; |
606 |
// Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); |
607 |
// without overflow |
608 |
{ |
609 |
Scalar t0 = m_matT.coeff(i+1, i); |
610 |
Scalar t1 = m_matT.coeff(i, i+1); |
611 |
Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1))); |
612 |
t0 /= maxval; |
613 |
t1 /= maxval; |
614 |
Scalar p0 = p/maxval; |
615 |
z = maxval * sqrt(abs(p0 * p0 + t0 * t1)); |
616 |
} |
617 |
m_eigenvalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z); |
618 |
m_eigenvalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z); |
619 |
if(!((isfinite)(m_eigenvalues.coeffRef(i)) && (isfinite)(m_eigenvalues.coeffRef(i+1)))) |
620 |
{ |
621 |
return NumericalIssue; |
622 |
} |
623 |
i += 2; |
624 |
} |
625 |
} |
626 |
return Success; |
627 |
} |
628 |
|
629 |
|
544 |
} // end namespace Eigen |
630 |
} // end namespace Eigen |
545 |
|
631 |
|
546 |
#endif // EIGEN_REAL_SCHUR_H |
632 |
#endif // EIGEN_REAL_SCHUR_H |