Bugzilla – Attachment 843 Details for
Bug 1535
Allow extraction of eigenvalues from a Schur matrixT()
Home
|
New
|
Browse
|
Search
|
[?]
|
Reports
|
Requests
|
Help
|
Log In
[x]
|
Forgot Password
Login:
[x]
This bugzilla service is closed. All entries have been migrated to
https://gitlab.com/libeigen/eigen
[patch]
Implement RealSchur::eigenvalues()
schur-eigenvalues.1.diff (text/plain), 7.31 KB, created by
Jan van Dijk
on 2018-04-04 09:10:11 UTC
(
hide
)
Description:
Implement RealSchur::eigenvalues()
Filename:
MIME Type:
Creator:
Jan van Dijk
Created:
2018-04-04 09:10:11 UTC
Size:
7.31 KB
patch
obsolete
>diff -r a9198af72a41 Eigen/src/Eigenvalues/RealSchur.h >--- a/Eigen/src/Eigenvalues/RealSchur.h Tue Apr 03 23:16:43 2018 +0200 >+++ b/Eigen/src/Eigenvalues/RealSchur.h Wed Apr 04 10:55:17 2018 +0200 >@@ -83,6 +83,7 @@ > explicit RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) > : m_matT(size, size), > m_matU(size, size), >+ m_eigenvalues(size), > m_workspaceVector(size), > m_hess(size), > m_isInitialized(false), >@@ -93,9 +94,10 @@ > /** \brief Constructor; computes real Schur decomposition of given matrix. > * > * \param[in] matrix Square matrix whose Schur decomposition is to be computed. >- * \param[in] computeU If true, both T and U are computed; if false, only T is computed. >+ * \param[in] computeU If true, both T and U are computed; if false, only T and the eigenvalues are computed. > * >- * This constructor calls compute() to compute the Schur decomposition. >+ * This constructor calls compute() to compute the Schur decomposition and the >+ * eigenvalues of the matrix. > * > * Example: \include RealSchur_RealSchur_MatrixType.cpp > * Output: \verbinclude RealSchur_RealSchur_MatrixType.out >@@ -147,10 +149,30 @@ > return m_matT; > } > >+ /** \brief Returns the eigenvalues of the given matrix. >+ * >+ * \returns A column vector containing the eigenvalues. >+ * >+ * \pre Either the constructor RealSchur(const MatrixType&, bool) or the >+ * member function compute(const MatrixType&, bool) has been called before >+ * to compute the Schur decomposition of a matrix. >+ * >+ * The eigenvalues of the matrix are the same as the eigenvalues of the 1x1 >+ * and 2x2 diagonal blocks of the matrix T and are extracted from that matrix >+ * after a factorization has been computed. >+ * >+ * \sa RealSchur(const MatrixType&, bool) for an example >+ */ >+ const EigenvalueType& eigenvalues() const >+ { >+ eigen_assert(m_isInitialized && "RealSchur is not initialized."); >+ return m_eigenvalues; >+ } >+ > /** \brief Computes Schur decomposition of given matrix. > * > * \param[in] matrix Square matrix whose Schur decomposition is to be computed. >- * \param[in] computeU If true, both T and U are computed; if false, only T is computed. >+ * \param[in] computeU If true, both T and U are computed; if false, only T and the eigenvalues are computed. > * \returns Reference to \c *this > * > * The Schur decomposition is computed by first reducing the matrix to >@@ -178,8 +200,9 @@ > * This routine assumes that the matrix is already reduced in Hessenberg form matrixH > * using either the class HessenbergDecomposition or another mean. > * It computes the upper quasi-triangular matrix T of the Schur decomposition of H >+ * and calculates the eigenvalues of T. > * When computeU is true, this routine computes the matrix U such that >- * A = U T U^T = (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix >+ * A = U T U^T = (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix. > * > * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix > * is not available, the user should give an identity matrix (Q.setIdentity()) >@@ -226,6 +249,7 @@ > > MatrixType m_matT; > MatrixType m_matU; >+ EigenvalueType m_eigenvalues; > ColumnVectorType m_workspaceVector; > HessenbergDecomposition<MatrixType> m_hess; > ComputationInfo m_info; >@@ -241,6 +265,11 @@ > void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo); > void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector); > void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace); >+ /** \brief Resizes vector m_eigenvalues and populates it with the >+ * eigenvalues of quasi-triangular matrix m_T. >+ * \return NumericalIssue if a non-finite eigenvalue was detected, Success otherwise >+ */ >+ ComputationInfo extractEigenvalues(); > }; > > >@@ -259,6 +288,7 @@ > if(scale<considerAsZero) > { > m_matT.setZero(matrix.rows(),matrix.cols()); >+ m_eigenvalues.setZero(matrix.rows()); > if(computeU) > m_matU.setIdentity(matrix.rows(),matrix.cols()); > m_info = Success; >@@ -274,9 +304,11 @@ > computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU); > > m_matT *= scale; >- >+ m_eigenvalues *= scale; >+ > return *this; > } >+ > template<typename MatrixType> > template<typename HessMatrixType, typename OrthMatrixType> > RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) >@@ -343,6 +375,13 @@ > else > m_info = NoConvergence; > >+ ComputationInfo ev_info = extractEigenvalues(); >+ // If the factorization failed, we keep the original message (NoConvergence...) >+ // Otherwise, set m_info to the result of extractEigenvalues(), >+ // which may return NumericalIssue if a NaN/infinity is found. >+ if (m_info == Success) >+ m_info = ev_info; >+ > m_isInitialized = true; > m_matUisUptodate = computeU; > return *this; >@@ -541,6 +580,53 @@ > } > } > >+template<typename MatrixType> >+ComputationInfo RealSchur<MatrixType>::extractEigenvalues() >+{ >+ using std::sqrt; >+ using std::abs; >+ using numext::isfinite; >+ m_eigenvalues.setZero(m_matT.cols()); >+ Index i = 0; >+ while (i < m_matT.cols()) >+ { >+ if (i == m_matT.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) >+ { >+ m_eigenvalues.coeffRef(i) = m_matT.coeff(i, i); >+ if(!(isfinite)(m_eigenvalues.coeffRef(i))) >+ { >+ return NumericalIssue; >+ } >+ ++i; >+ } >+ else >+ { >+ Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1)); >+ Scalar z; >+ // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); >+ // without overflow >+ { >+ Scalar t0 = m_matT.coeff(i+1, i); >+ Scalar t1 = m_matT.coeff(i, i+1); >+ Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1))); >+ t0 /= maxval; >+ t1 /= maxval; >+ Scalar p0 = p/maxval; >+ z = maxval * sqrt(abs(p0 * p0 + t0 * t1)); >+ } >+ m_eigenvalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z); >+ m_eigenvalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z); >+ if(!((isfinite)(m_eigenvalues.coeffRef(i)) && (isfinite)(m_eigenvalues.coeffRef(i+1)))) >+ { >+ return NumericalIssue; >+ } >+ i += 2; >+ } >+ } >+ return Success; >+} >+ >+ > } // end namespace Eigen > > #endif // EIGEN_REAL_SCHUR_H >diff -r a9198af72a41 doc/snippets/RealSchur_RealSchur_MatrixType.cpp >--- a/doc/snippets/RealSchur_RealSchur_MatrixType.cpp Tue Apr 03 23:16:43 2018 +0200 >+++ b/doc/snippets/RealSchur_RealSchur_MatrixType.cpp Wed Apr 04 10:55:17 2018 +0200 >@@ -4,6 +4,7 @@ > RealSchur<MatrixXd> schur(A); > cout << "The orthogonal matrix U is:" << endl << schur.matrixU() << endl; > cout << "The quasi-triangular matrix T is:" << endl << schur.matrixT() << endl << endl; >+cout << "The eigenvalues are:" << endl << schur.eigenvalues() << endl << endl; > > MatrixXd U = schur.matrixU(); > MatrixXd T = schur.matrixT();
You cannot view the attachment while viewing its details because your browser does not support IFRAMEs.
View the attachment on a separate page
.
View Attachment As Diff
View Attachment As Raw
Actions:
View
|
Diff
Attachments on
bug 1535
: 843 |
844