aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/Eigenvalues/ComplexEigenSolver.h36
-rw-r--r--Eigen/src/Eigenvalues/ComplexSchur.h41
-rw-r--r--Eigen/src/Eigenvalues/EigenSolver.h36
-rw-r--r--Eigen/src/Eigenvalues/RealSchur.h32
-rw-r--r--test/eigensolver_complex.cpp10
-rw-r--r--test/eigensolver_generic.cpp12
-rw-r--r--test/schur_complex.cpp18
-rw-r--r--test/schur_real.cpp20
8 files changed, 178 insertions, 27 deletions
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index c4b8a308c..39dfe2a0a 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -3,7 +3,7 @@
//
// Copyright (C) 2009 Claire Maurice
// Copyright (C) 2009 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
@@ -208,7 +208,27 @@ template<typename _MatrixType> class ComplexEigenSolver
* Example: \include ComplexEigenSolver_compute.cpp
* Output: \verbinclude ComplexEigenSolver_compute.out
*/
- ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
+ ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true)
+ {
+ return computeInternal(matrix, computeEigenvectors, false, 0);
+ }
+
+ /** \brief Computes eigendecomposition with specified maximum number of iterations.
+ *
+ * \param[in] matrix Square matrix whose eigendecomposition is to be computed.
+ * \param[in] computeEigenvectors If true, both the eigenvectors and the
+ * eigenvalues are computed; if false, only the eigenvalues are
+ * computed.
+ * \param[in] maxIter Maximum number of iterations.
+ * \returns Reference to \c *this
+ *
+ * This method provides the same functionality as compute(const MatrixType&, bool) but it also allows the
+ * user to specify the maximum number of iterations to be used when computing the Schur decomposition.
+ */
+ ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors, Index maxIter)
+ {
+ return computeInternal(matrix, computeEigenvectors, true, maxIter);
+ }
/** \brief Reports whether previous computation was successful.
*
@@ -231,18 +251,26 @@ template<typename _MatrixType> class ComplexEigenSolver
private:
void doComputeEigenvectors(RealScalar matrixnorm);
void sortEigenvalues(bool computeEigenvectors);
+ ComplexEigenSolver& computeInternal(const MatrixType& matrix, bool computeEigenvectors,
+ bool maxIterSpecified, Index maxIter);
+
};
template<typename MatrixType>
-ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
+ComplexEigenSolver<MatrixType>&
+ComplexEigenSolver<MatrixType>::computeInternal(const MatrixType& matrix, bool computeEigenvectors,
+ bool maxIterSpecified, Index maxIter)
{
// this code is inspired from Jampack
assert(matrix.cols() == matrix.rows());
// Do a complex Schur decomposition, A = U T U^*
// The eigenvalues are on the diagonal of T.
- m_schur.compute(matrix, computeEigenvectors);
+ if (maxIterSpecified)
+ m_schur.compute(matrix, computeEigenvectors, maxIter);
+ else
+ m_schur.compute(matrix, computeEigenvectors);
if(m_schur.info() == Success)
{
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h
index 55aeedb90..710a8a276 100644
--- a/Eigen/src/Eigenvalues/ComplexSchur.h
+++ b/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -3,7 +3,7 @@
//
// Copyright (C) 2009 Claire Maurice
// Copyright (C) 2009 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
@@ -166,6 +166,7 @@ template<typename _MatrixType> class ComplexSchur
*
* \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.
+
* \returns Reference to \c *this
*
* The Schur decomposition is computed by first reducing the
@@ -180,8 +181,27 @@ template<typename _MatrixType> class ComplexSchur
*
* Example: \include ComplexSchur_compute.cpp
* Output: \verbinclude ComplexSchur_compute.out
+ *
+ * \sa compute(const MatrixType&, bool, Index)
*/
- ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
+ ComplexSchur& compute(const MatrixType& matrix, bool computeU = true)
+ {
+ return compute(matrix, computeU, m_maxIterations * matrix.rows());
+ }
+
+ /** \brief Computes Schur decomposition with specified maximum number of iterations.
+ *
+ * \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] maxIter Maximum number of iterations.
+ *
+ * \returns Reference to \c *this
+ *
+ * This method provides the same functionality as compute(const MatrixType&, bool) but it also allows the
+ * user to specify the maximum number of QR iterations to be used. The maximum number of iterations for
+ * compute(const MatrixType&, bool) is #m_maxIterations times the size of the matrix.
+ */
+ ComplexSchur& compute(const MatrixType& matrix, bool computeU, Index maxIter);
/** \brief Reports whether previous computation was successful.
*
@@ -189,13 +209,14 @@ template<typename _MatrixType> class ComplexSchur
*/
ComputationInfo info() const
{
- eigen_assert(m_isInitialized && "RealSchur is not initialized.");
+ eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
return m_info;
}
/** \brief Maximum number of iterations.
*
- * 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 30.
*/
static const int m_maxIterations = 30;
@@ -209,7 +230,7 @@ template<typename _MatrixType> class ComplexSchur
private:
bool subdiagonalEntryIsNeglegible(Index i);
ComplexScalar computeShift(Index iu, Index iter);
- void reduceToTriangularForm(bool computeU);
+ void reduceToTriangularForm(bool computeU, Index maxIter);
friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
};
@@ -268,7 +289,7 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
template<typename MatrixType>
-ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
+ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU, Index maxIter)
{
m_matUisUptodate = false;
eigen_assert(matrix.cols() == matrix.rows());
@@ -284,7 +305,7 @@ ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& ma
}
internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
- reduceToTriangularForm(computeU);
+ reduceToTriangularForm(computeU, maxIter);
return *this;
}
@@ -327,7 +348,7 @@ struct complex_schur_reduce_to_hessenberg<MatrixType, false>
// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
template<typename MatrixType>
-void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
+void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU, Index maxIter)
{
// The matrix m_matT is divided in three parts.
// Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
@@ -354,7 +375,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
// if we spent too many iterations, we give up
iter++;
totalIter++;
- if(totalIter > m_maxIterations * m_matT.cols()) break;
+ if(totalIter > maxIter) break;
// find il, the top row of the active submatrix
il = iu-1;
@@ -384,7 +405,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
}
}
- if(totalIter <= m_maxIterations * m_matT.cols())
+ if(totalIter <= maxIter)
m_info = Success;
else
m_info = NoConvergence;
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index c16ff2b74..315958421 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.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
@@ -273,7 +273,27 @@ template<typename _MatrixType> class EigenSolver
* Example: \include EigenSolver_compute.cpp
* Output: \verbinclude EigenSolver_compute.out
*/
- EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
+ EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true)
+ {
+ return computeInternal(matrix, computeEigenvectors, false, 0);
+ }
+
+ /** \brief Computes eigendecomposition with specified maximum number of iterations.
+ *
+ * \param[in] matrix Square matrix whose eigendecomposition is to be computed.
+ * \param[in] computeEigenvectors If true, both the eigenvectors and the
+ * eigenvalues are computed; if false, only the eigenvalues are
+ * computed.
+ * \param[in] maxIter Maximum number of iterations.
+ * \returns Reference to \c *this
+ *
+ * This method provides the same functionality as compute(const MatrixType&, bool) but it also allows the
+ * user to specify the maximum number of iterations to be used when computing the Schur decomposition.
+ */
+ EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors, Index maxIter)
+ {
+ return computeInternal(matrix, computeEigenvectors, true, maxIter);
+ }
ComputationInfo info() const
{
@@ -283,6 +303,8 @@ template<typename _MatrixType> class EigenSolver
private:
void doComputeEigenvectors();
+ EigenSolver& computeInternal(const MatrixType& matrix, bool computeEigenvectors,
+ bool maxIterSpecified, Index maxIter);
protected:
MatrixType m_eivec;
@@ -348,12 +370,18 @@ typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige
}
template<typename MatrixType>
-EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
+EigenSolver<MatrixType>&
+EigenSolver<MatrixType>::computeInternal(const MatrixType& matrix, bool computeEigenvectors,
+ bool maxIterSpecified, Index maxIter)
{
assert(matrix.cols() == matrix.rows());
// Reduce to real Schur form.
- m_realSchur.compute(matrix, computeEigenvectors);
+ if (maxIterSpecified)
+ m_realSchur.compute(matrix, computeEigenvectors, maxIter);
+ else
+ m_realSchur.compute(matrix, computeEigenvectors);
+
if (m_realSchur.info() == Success)
{
m_matT = m_realSchur.matrixT();
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index d1949b83c..da877dd04 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
@@ -160,8 +160,27 @@ 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)
+ {
+ return compute(matrix, computeU, m_maxIterations * matrix.rows());
+ }
+
+ /** \brief Computes Schur decomposition with specified maximum number of iterations.
+ *
+ * \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] maxIter Maximum number of iterations.
+ *
+ * \returns Reference to \c *this
+ *
+ * This method provides the same functionality as compute(const MatrixType&, bool) but it also allows the
+ * user to specify the maximum number of QR iterations to be used. The maximum number of iterations for
+ * compute(const MatrixType&, bool) is #m_maxIterations times the size of the matrix.
*/
- RealSchur& compute(const MatrixType& matrix, bool computeU = true);
+ RealSchur& compute(const MatrixType& matrix, bool computeU, Index maxIter);
/** \brief Reports whether previous computation was successful.
*
@@ -175,7 +194,8 @@ template<typename _MatrixType> class RealSchur
/** \brief Maximum number of iterations.
*
- * 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;
@@ -201,7 +221,7 @@ template<typename _MatrixType> class RealSchur
template<typename MatrixType>
-RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
+RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU, Index maxIter)
{
assert(matrix.cols() == matrix.rows());
@@ -253,14 +273,14 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix,
computeShift(iu, iter, exshift, shiftInfo);
iter = iter + 1;
totalIter = totalIter + 1;
- if (totalIter > m_maxIterations * matrix.cols()) break;
+ if (totalIter > maxIter) break;
Index im;
initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
}
}
}
- if(totalIter <= m_maxIterations * matrix.cols())
+ if(totalIter <= maxIter)
m_info = Success;
else
m_info = NoConvergence;
diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp
index 0c2059512..a569d8507 100644
--- a/test/eigensolver_complex.cpp
+++ b/test/eigensolver_complex.cpp
@@ -59,6 +59,16 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
// another algorithm so results may differ slightly
verify_is_approx_upto_permutation(a.eigenvalues(), ei1.eigenvalues());
+ ComplexEigenSolver<MatrixType> ei2;
+ ei2.compute(a, true, ComplexSchur<MatrixType>::m_maxIterations * rows);
+ VERIFY_IS_EQUAL(ei2.info(), Success);
+ VERIFY_IS_EQUAL(ei2.eigenvectors(), ei1.eigenvectors());
+ VERIFY_IS_EQUAL(ei2.eigenvalues(), ei1.eigenvalues());
+ if (rows > 2) {
+ ei2.compute(a, true, 1);
+ VERIFY_IS_EQUAL(ei2.info(), NoConvergence);
+ }
+
ComplexEigenSolver<MatrixType> eiNoEivecs(a, false);
VERIFY_IS_EQUAL(eiNoEivecs.info(), Success);
VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp
index 0b55ccd93..1b83afbc4 100644
--- a/test/eigensolver_generic.cpp
+++ b/test/eigensolver_generic.cpp
@@ -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
@@ -45,6 +45,16 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
VERIFY_IS_APPROX(ei1.eigenvectors().colwise().norm(), RealVectorType::Ones(rows).transpose());
VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues());
+ EigenSolver<MatrixType> ei2;
+ ei2.compute(a, true, RealSchur<MatrixType>::m_maxIterations * rows);
+ VERIFY_IS_EQUAL(ei2.info(), Success);
+ VERIFY_IS_EQUAL(ei2.eigenvectors(), ei1.eigenvectors());
+ VERIFY_IS_EQUAL(ei2.eigenvalues(), ei1.eigenvalues());
+ if (rows > 2) {
+ ei2.compute(a, true, 1);
+ VERIFY_IS_EQUAL(ei2.info(), NoConvergence);
+ }
+
EigenSolver<MatrixType> eiNoEivecs(a, false);
VERIFY_IS_EQUAL(eiNoEivecs.info(), Success);
VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
diff --git a/test/schur_complex.cpp b/test/schur_complex.cpp
index a6f66ab02..4d5cbbfdf 100644
--- a/test/schur_complex.cpp
+++ b/test/schur_complex.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// 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
@@ -47,6 +47,22 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
VERIFY_IS_EQUAL(cs1.matrixT(), cs2.matrixT());
VERIFY_IS_EQUAL(cs1.matrixU(), cs2.matrixU());
+ // Test maximum number of iterations
+ ComplexSchur<MatrixType> cs3;
+ cs3.compute(A, true, ComplexSchur<MatrixType>::m_maxIterations * size);
+ VERIFY_IS_EQUAL(cs3.info(), Success);
+ VERIFY_IS_EQUAL(cs3.matrixT(), cs1.matrixT());
+ VERIFY_IS_EQUAL(cs3.matrixU(), cs1.matrixU());
+ cs3.compute(A, true, 1);
+ VERIFY_IS_EQUAL(cs3.info(), size > 1 ? NoConvergence : Success);
+
+ MatrixType Atriangular = A;
+ Atriangular.template triangularView<StrictlyLower>().setZero();
+ cs3.compute(Atriangular, true, 1); // triangular matrices do not need any iterations
+ VERIFY_IS_EQUAL(cs3.info(), Success);
+ VERIFY_IS_EQUAL(cs3.matrixT(), Atriangular.template cast<ComplexScalar>());
+ VERIFY_IS_EQUAL(cs3.matrixU(), ComplexMatrixType::Identity(size, size));
+
// Test computation of only T, not U
ComplexSchur<MatrixType> csOnlyT(A, false);
VERIFY_IS_EQUAL(csOnlyT.info(), Success);
diff --git a/test/schur_real.cpp b/test/schur_real.cpp
index e6351d94a..e81a28d7d 100644
--- a/test/schur_real.cpp
+++ b/test/schur_real.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// 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
@@ -66,6 +66,24 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
VERIFY_IS_EQUAL(rs1.matrixT(), rs2.matrixT());
VERIFY_IS_EQUAL(rs1.matrixU(), rs2.matrixU());
+ // Test maximum number of iterations
+ RealSchur<MatrixType> rs3;
+ rs3.compute(A, true, RealSchur<MatrixType>::m_maxIterations * size);
+ VERIFY_IS_EQUAL(rs3.info(), Success);
+ VERIFY_IS_EQUAL(rs3.matrixT(), rs1.matrixT());
+ VERIFY_IS_EQUAL(rs3.matrixU(), rs1.matrixU());
+ if (size > 2) {
+ rs3.compute(A, true, 1);
+ VERIFY_IS_EQUAL(rs3.info(), NoConvergence);
+ }
+
+ MatrixType Atriangular = A;
+ Atriangular.template triangularView<StrictlyLower>().setZero();
+ rs3.compute(Atriangular, true, 1); // triangular matrices do not need any iterations
+ VERIFY_IS_EQUAL(rs3.info(), Success);
+ VERIFY_IS_EQUAL(rs3.matrixT(), Atriangular);
+ VERIFY_IS_EQUAL(rs3.matrixU(), MatrixType::Identity(size, size));
+
// Test computation of only T, not U
RealSchur<MatrixType> rsOnlyT(A, false);
VERIFY_IS_EQUAL(rsOnlyT.info(), Success);
generated by cgit v1.2.3 (git 2.39.1) at 2024-12-16 13:27:25 +0000