aboutsummaryrefslogtreecommitdiff
path: root/Eigen/src/Eigenvalues/RealSchur.h
diff options
context:
space:
mode:
Diffstat (limited to 'Eigen/src/Eigenvalues/RealSchur.h')
-rw-r--r--Eigen/src/Eigenvalues/RealSchur.h123
1 files changed, 94 insertions, 29 deletions
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index 781692ecc..64d136341 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -2,7 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -86,7 +86,8 @@ template<typename _MatrixType> class RealSchur
m_workspaceVector(size),
m_hess(size),
m_isInitialized(false),
- m_matUisUptodate(false)
+ m_matUisUptodate(false),
+ m_maxIters(-1)
{ }
/** \brief Constructor; computes real Schur decomposition of given matrix.
@@ -105,7 +106,8 @@ template<typename _MatrixType> class RealSchur
m_workspaceVector(matrix.rows()),
m_hess(matrix.rows()),
m_isInitialized(false),
- m_matUisUptodate(false)
+ m_matUisUptodate(false),
+ m_maxIters(-1)
{
compute(matrix, computeU);
}
@@ -160,9 +162,30 @@ template<typename _MatrixType> class RealSchur
*
* Example: \include RealSchur_compute.cpp
* Output: \verbinclude RealSchur_compute.out
+ *
+ * \sa compute(const MatrixType&, bool, Index)
*/
RealSchur& compute(const MatrixType& matrix, bool computeU = true);
+ /** \brief Computes Schur decomposition of a Hessenberg matrix H = Z T Z^T
+ * \param[in] matrixH Matrix in Hessenberg form H
+ * \param[in] matrixQ orthogonal matrix Q that transform a matrix A to H : A = Q H Q^T
+ * \param computeU Computes the matriX U of the Schur vectors
+ * \return Reference to \c *this
+ *
+ * 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
+ * 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
+ *
+ * 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())
+ *
+ * \sa compute(const MatrixType&, bool)
+ */
+ template<typename HessMatrixType, typename OrthMatrixType>
+ RealSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU);
/** \brief Reports whether previous computation was successful.
*
* \returns \c Success if computation was succesful, \c NoConvergence otherwise.
@@ -173,11 +196,29 @@ template<typename _MatrixType> class RealSchur
return m_info;
}
- /** \brief Maximum number of iterations.
+ /** \brief Sets the maximum number of iterations allowed.
+ *
+ * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size
+ * of the matrix.
+ */
+ RealSchur& setMaxIterations(Index maxIters)
+ {
+ m_maxIters = maxIters;
+ return *this;
+ }
+
+ /** \brief Returns the maximum number of iterations. */
+ Index getMaxIterations()
+ {
+ return m_maxIters;
+ }
+
+ /** \brief Maximum number of iterations per row.
*
- * Maximum number of iterations allowed for an eigenvalue to converge.
+ * If not otherwise specified, the maximum number of iterations is this number times the size of the
+ * matrix. It is currently set to 40.
*/
- static const int m_maxIterations = 40;
+ static const int m_maxIterationsPerRow = 40;
private:
@@ -188,12 +229,13 @@ template<typename _MatrixType> class RealSchur
ComputationInfo m_info;
bool m_isInitialized;
bool m_matUisUptodate;
+ Index m_maxIters;
typedef Matrix<Scalar,3,1> Vector3s;
Scalar computeNormOfT();
- Index findSmallSubdiagEntry(Index iu, Scalar norm);
- void splitOffTwoRows(Index iu, bool computeU, Scalar exshift);
+ Index findSmallSubdiagEntry(Index iu, const Scalar& norm);
+ void splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift);
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);
@@ -203,15 +245,30 @@ template<typename _MatrixType> class RealSchur
template<typename MatrixType>
RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
{
- assert(matrix.cols() == matrix.rows());
+ eigen_assert(matrix.cols() == matrix.rows());
+ Index maxIters = m_maxIters;
+ if (maxIters == -1)
+ maxIters = m_maxIterationsPerRow * matrix.rows();
// Step 1. Reduce to Hessenberg form
m_hess.compute(matrix);
- m_matT = m_hess.matrixH();
- if (computeU)
- m_matU = m_hess.matrixQ();
// Step 2. Reduce to real Schur form
+ computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU);
+
+ return *this;
+}
+template<typename MatrixType>
+template<typename HessMatrixType, typename OrthMatrixType>
+RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
+{
+ m_matT = matrixH;
+ if(computeU)
+ m_matU = matrixQ;
+
+ Index maxIters = m_maxIters;
+ if (maxIters == -1)
+ maxIters = m_maxIterationsPerRow * matrixH.rows();
m_workspaceVector.resize(m_matT.cols());
Scalar* workspace = &m_workspaceVector.coeffRef(0);
@@ -220,8 +277,9 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix,
// Rows il,...,iu is the part we are working on (the active window).
// Rows iu+1,...,end are already brought in triangular form.
Index iu = m_matT.cols() - 1;
- Index iter = 0; // iteration count
- Scalar exshift(0); // sum of exceptional shifts
+ Index iter = 0; // iteration count for current eigenvalue
+ Index totalIter = 0; // iteration count for whole matrix
+ Scalar exshift(0); // sum of exceptional shifts
Scalar norm = computeNormOfT();
if(norm!=0)
@@ -251,14 +309,15 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix,
Vector3s firstHouseholderVector(0,0,0), shiftInfo;
computeShift(iu, iter, exshift, shiftInfo);
iter = iter + 1;
- if (iter > m_maxIterations * m_matT.cols()) break;
+ totalIter = totalIter + 1;
+ if (totalIter > maxIters) break;
Index im;
initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
}
}
}
- if(iter <= m_maxIterations * m_matT.cols())
+ if(totalIter <= maxIters)
m_info = Success;
else
m_info = NoConvergence;
@@ -278,21 +337,22 @@ inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
// + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
Scalar norm(0);
for (Index j = 0; j < size; ++j)
- norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
+ norm += m_matT.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum();
return norm;
}
/** \internal Look for single small sub-diagonal element and returns its index */
template<typename MatrixType>
-inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, Scalar norm)
+inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, const Scalar& norm)
{
+ using std::abs;
Index res = iu;
while (res > 0)
{
- Scalar s = internal::abs(m_matT.coeff(res-1,res-1)) + internal::abs(m_matT.coeff(res,res));
+ Scalar s = abs(m_matT.coeff(res-1,res-1)) + abs(m_matT.coeff(res,res));
if (s == 0.0)
s = norm;
- if (internal::abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
+ if (abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
break;
res--;
}
@@ -301,8 +361,10 @@ inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(I
/** \internal Update T given that rows iu-1 and iu decouple from the rest. */
template<typename MatrixType>
-inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scalar exshift)
+inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, const Scalar& exshift)
{
+ using std::sqrt;
+ using std::abs;
const Index size = m_matT.cols();
// The eigenvalues of the 2x2 matrix [a b; c d] are
@@ -314,7 +376,7 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scal
if (q >= Scalar(0)) // Two real eigenvalues
{
- Scalar z = internal::sqrt(internal::abs(q));
+ Scalar z = sqrt(abs(q));
JacobiRotation<Scalar> rot;
if (p >= Scalar(0))
rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
@@ -336,6 +398,8 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scal
template<typename MatrixType>
inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
{
+ using std::sqrt;
+ using std::abs;
shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
@@ -346,7 +410,7 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
exshift += shiftInfo.coeff(0);
for (Index i = 0; i <= iu; ++i)
m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
- Scalar s = internal::abs(m_matT.coeff(iu,iu-1)) + internal::abs(m_matT.coeff(iu-1,iu-2));
+ Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2));
shiftInfo.coeffRef(0) = Scalar(0.75) * s;
shiftInfo.coeffRef(1) = Scalar(0.75) * s;
shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
@@ -359,7 +423,7 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
s = s * s + shiftInfo.coeff(2);
if (s > Scalar(0))
{
- s = internal::sqrt(s);
+ s = sqrt(s);
if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
s = -s;
s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
@@ -376,6 +440,7 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
template<typename MatrixType>
inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
{
+ using std::abs;
Vector3s& v = firstHouseholderVector; // alias to save typing
for (im = iu-2; im >= il; --im)
@@ -389,9 +454,9 @@ inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const V
if (im == il) {
break;
}
- const Scalar lhs = m_matT.coeff(im,im-1) * (internal::abs(v.coeff(1)) + internal::abs(v.coeff(2)));
- const Scalar rhs = v.coeff(0) * (internal::abs(m_matT.coeff(im-1,im-1)) + internal::abs(Tmm) + internal::abs(m_matT.coeff(im+1,im+1)));
- if (internal::abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
+ const Scalar lhs = m_matT.coeff(im,im-1) * (abs(v.coeff(1)) + abs(v.coeff(2)));
+ const Scalar rhs = v.coeff(0) * (abs(m_matT.coeff(im-1,im-1)) + abs(Tmm) + abs(m_matT.coeff(im+1,im+1)));
+ if (abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
{
break;
}
@@ -402,8 +467,8 @@ inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const V
template<typename MatrixType>
inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
{
- assert(im >= il);
- assert(im <= iu-2);
+ eigen_assert(im >= il);
+ eigen_assert(im <= iu-2);
const Index size = m_matT.cols();
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-06 14:52:28 +0000