GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-03 10:14:04 +08:00
Code Issues Packages Projects Releases Wiki Activity

Eigenvalues module: Implement setMaxIterations() methods.

Browse Source
This commit is contained in:
Jitse Niesen 2012-07-28 21:30:09 +01:00
parent 6f54269829
commit 696b2f999f
8 changed files with 116 additions and 126 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
Eigen/src/Eigenvalues/ComplexEigenSolver.h
View File

@ -208,27 +208,7 @@ template<typename _MatrixType> class ComplexEigenSolver
* Example: \include ComplexEigenSolver_compute.cpp
* Output: \verbinclude ComplexEigenSolver_compute.out
*/
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);
}
ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
/** \brief Reports whether previous computation was successful.
*
@ -240,6 +220,19 @@ template<typename _MatrixType> class ComplexEigenSolver
return m_schur.info();
}
/** \brief Sets the maximum number of iterations allowed. */
ComplexEigenSolver& setMaxIterations(Index maxIters)
{
m_schur.setMaxIterations(maxIters);
return *this;
}
/** \brief Returns the maximum number of iterations. */
Index getMaxIterations()
{
return m_schur.getMaxIterations();
}
protected:
EigenvectorType m_eivec;
EigenvalueType m_eivalues;
@ -251,26 +244,19 @@ 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>::computeInternal(const MatrixType& matrix, bool computeEigenvectors,
bool maxIterSpecified, Index maxIter)
ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
{
// 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.
if (maxIterSpecified)
m_schur.compute(matrix, computeEigenvectors, maxIter);
else
m_schur.compute(matrix, computeEigenvectors);
m_schur.compute(matrix, computeEigenvectors);
if(m_schur.info() == Success)
{

71
Eigen/src/Eigenvalues/ComplexSchur.h
View File

@ -96,7 +96,8 @@ template<typename _MatrixType> class ComplexSchur
m_matU(size,size),
m_hess(size),
m_isInitialized(false),
m_matUisUptodate(false)
m_matUisUptodate(false),
m_maxIters(-1)
{}
/** \brief Constructor; computes Schur decomposition of given matrix.
@ -109,11 +110,12 @@ template<typename _MatrixType> class ComplexSchur
* \sa matrixT() and matrixU() for examples.
*/
ComplexSchur(const MatrixType& matrix, bool computeU = true)
: m_matT(matrix.rows(),matrix.cols()),
m_matU(matrix.rows(),matrix.cols()),
m_hess(matrix.rows()),
m_isInitialized(false),
m_matUisUptodate(false)
: m_matT(matrix.rows(),matrix.cols()),
m_matU(matrix.rows(),matrix.cols()),
m_hess(matrix.rows()),
m_isInitialized(false),
m_matUisUptodate(false),
m_maxIters(-1)
{
compute(matrix, computeU);
}
@ -184,24 +186,7 @@ template<typename _MatrixType> class ComplexSchur
*
* \sa compute(const MatrixType&, bool, Index)
*/
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);
ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
/** \brief Reports whether previous computation was successful.
*
@ -213,12 +198,29 @@ template<typename _MatrixType> class ComplexSchur
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.
*/
ComplexSchur& 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.
*
* 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;
static const int m_maxIterationsPerRow = 30;
protected:
ComplexMatrixType m_matT, m_matU;
@ -226,11 +228,12 @@ template<typename _MatrixType> class ComplexSchur
ComputationInfo m_info;
bool m_isInitialized;
bool m_matUisUptodate;
Index m_maxIters;
private:
bool subdiagonalEntryIsNeglegible(Index i);
ComplexScalar computeShift(Index iu, Index iter);
void reduceToTriangularForm(bool computeU, Index maxIter);
void reduceToTriangularForm(bool computeU);
friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
};
@ -289,7 +292,7 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
template<typename MatrixType>
ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU, Index maxIter)
ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
{
m_matUisUptodate = false;
eigen_assert(matrix.cols() == matrix.rows());
@ -305,7 +308,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, maxIter);
reduceToTriangularForm(computeU);
return *this;
}
@ -348,8 +351,12 @@ 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, Index maxIter)
void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
{
Index maxIters = m_maxIters;
if (maxIters == -1)
maxIters = m_maxIterationsPerRow * m_matT.rows();
// 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.
// Rows il,...,iu is the part we are working on (the active submatrix).
@ -375,7 +382,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU, Index maxIt
// if we spent too many iterations, we give up
iter++;
totalIter++;
if(totalIter > maxIter) break;
if(totalIter > maxIters) break;
// find il, the top row of the active submatrix
il = iu-1;
@ -405,7 +412,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU, Index maxIt
}
}
if(totalIter <= maxIter)
if(totalIter <= maxIters)
m_info = Success;
else
m_info = NoConvergence;

45
Eigen/src/Eigenvalues/EigenSolver.h
View File

@ -273,27 +273,7 @@ template<typename _MatrixType> class EigenSolver
* Example: \include EigenSolver_compute.cpp
* Output: \verbinclude EigenSolver_compute.out
*/
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);
}
EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
ComputationInfo info() const
{
@ -301,10 +281,21 @@ template<typename _MatrixType> class EigenSolver
return m_realSchur.info();
}
/** \brief Sets the maximum number of iterations allowed. */
EigenSolver& setMaxIterations(Index maxIters)
{
m_realSchur.setMaxIterations(maxIters);
return *this;
}
/** \brief Returns the maximum number of iterations. */
Index getMaxIterations()
{
return m_realSchur.getMaxIterations();
}
private:
void doComputeEigenvectors();
EigenSolver& computeInternal(const MatrixType& matrix, bool computeEigenvectors,
bool maxIterSpecified, Index maxIter);
protected:
MatrixType m_eivec;
@ -371,16 +362,12 @@ typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige
template<typename MatrixType>
EigenSolver<MatrixType>&
EigenSolver<MatrixType>::computeInternal(const MatrixType& matrix, bool computeEigenvectors,
bool maxIterSpecified, Index maxIter)
EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
{
assert(matrix.cols() == matrix.rows());
// Reduce to real Schur form.
if (maxIterSpecified)
m_realSchur.compute(matrix, computeEigenvectors, maxIter);
else
m_realSchur.compute(matrix, computeEigenvectors);
m_realSchur.compute(matrix, computeEigenvectors);
if (m_realSchur.info() == Success)
{

56
Eigen/src/Eigenvalues/RealSchur.h
View File

@ -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);
}
@ -163,24 +165,7 @@ template<typename _MatrixType> class RealSchur
*
* \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, Index maxIter);
RealSchur& compute(const MatrixType& matrix, bool computeU = true);
/** \brief Reports whether previous computation was successful.
*
@ -192,12 +177,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.
*
* 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:
@ -208,6 +210,7 @@ template<typename _MatrixType> class RealSchur
ComputationInfo m_info;
bool m_isInitialized;
bool m_matUisUptodate;
Index m_maxIters;
typedef Matrix<Scalar,3,1> Vector3s;
@ -221,9 +224,12 @@ template<typename _MatrixType> class RealSchur
template<typename MatrixType>
RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU, Index maxIter)
RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
{
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);
@ -273,14 +279,14 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix,
computeShift(iu, iter, exshift, shiftInfo);
iter = iter + 1;
totalIter = totalIter + 1;
if (totalIter > maxIter) break;
if (totalIter > maxIters) break;
Index im;
initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
}
}
}
if(totalIter <= maxIter)
if(totalIter <= maxIters)
m_info = Success;
else
m_info = NoConvergence;

5
test/eigensolver_complex.cpp
View File

@ -60,13 +60,14 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
verify_is_approx_upto_permutation(a.eigenvalues(), ei1.eigenvalues());
ComplexEigenSolver<MatrixType> ei2;
ei2.compute(a, true, ComplexSchur<MatrixType>::m_maxIterations * rows);
ei2.setMaxIterations(ComplexSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a);
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);
ei2.setMaxIterations(1).compute(a);
VERIFY_IS_EQUAL(ei2.info(), NoConvergence);
VERIFY_IS_EQUAL(ei2.getMaxIterations(), 1);
}
ComplexEigenSolver<MatrixType> eiNoEivecs(a, false);

5
test/eigensolver_generic.cpp
View File

@ -46,13 +46,14 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues());
EigenSolver<MatrixType> ei2;
ei2.compute(a, true, RealSchur<MatrixType>::m_maxIterations * rows);
ei2.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * rows).compute(a);
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);
ei2.setMaxIterations(1).compute(a);
VERIFY_IS_EQUAL(ei2.info(), NoConvergence);
VERIFY_IS_EQUAL(ei2.getMaxIterations(), 1);
}
EigenSolver<MatrixType> eiNoEivecs(a, false);

7
test/schur_complex.cpp
View File

@ -49,16 +49,17 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
// Test maximum number of iterations
ComplexSchur<MatrixType> cs3;
cs3.compute(A, true, ComplexSchur<MatrixType>::m_maxIterations * size);
cs3.setMaxIterations(ComplexSchur<MatrixType>::m_maxIterationsPerRow * size).compute(A);
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);
cs3.setMaxIterations(1).compute(A);
VERIFY_IS_EQUAL(cs3.info(), size > 1 ? NoConvergence : Success);
VERIFY_IS_EQUAL(cs3.getMaxIterations(), 1);
MatrixType Atriangular = A;
Atriangular.template triangularView<StrictlyLower>().setZero();
cs3.compute(Atriangular, true, 1); // triangular matrices do not need any iterations
cs3.setMaxIterations(1).compute(Atriangular); // 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));

7
test/schur_real.cpp
View File

@ -68,18 +68,19 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
// Test maximum number of iterations
RealSchur<MatrixType> rs3;
rs3.compute(A, true, RealSchur<MatrixType>::m_maxIterations * size);
rs3.setMaxIterations(RealSchur<MatrixType>::m_maxIterationsPerRow * size).compute(A);
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);
rs3.setMaxIterations(1).compute(A);
VERIFY_IS_EQUAL(rs3.info(), NoConvergence);
VERIFY_IS_EQUAL(rs3.getMaxIterations(), 1);
}
MatrixType Atriangular = A;
Atriangular.template triangularView<StrictlyLower>().setZero();
rs3.compute(Atriangular, true, 1); // triangular matrices do not need any iterations
rs3.setMaxIterations(1).compute(Atriangular); // 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));