aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/Eigenvalues/ComplexEigenSolver.h24
-rw-r--r--Eigen/src/Eigenvalues/ComplexSchur.h27
-rw-r--r--Eigen/src/Eigenvalues/EigenSolver.h56
-rw-r--r--Eigen/src/Eigenvalues/EigenvaluesCommon.h39
-rw-r--r--Eigen/src/Eigenvalues/RealSchur.h27
-rw-r--r--Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h42
-rw-r--r--test/eigensolver_complex.cpp12
-rw-r--r--test/eigensolver_generic.cpp14
-rw-r--r--test/eigensolver_selfadjoint.cpp14
-rw-r--r--test/schur_complex.cpp11
-rw-r--r--test/schur_real.cpp11
11 files changed, 222 insertions, 55 deletions
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index a344c2d3c..a164aaae6 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -27,6 +27,9 @@
#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
#define EIGEN_COMPLEX_EIGEN_SOLVER_H
+#include "./EigenvaluesCommon.h"
+#include "./ComplexSchur.h"
+
/** \eigenvalues_module \ingroup Eigenvalues_Module
* \nonstableyet
*
@@ -220,6 +223,16 @@ template<typename _MatrixType> class ComplexEigenSolver
*/
ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
+ /** \brief Reports whether previous computation was successful.
+ *
+ * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
+ */
+ ComputationInfo info() const
+ {
+ ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
+ return m_schur.info();
+ }
+
protected:
EigenvectorType m_eivec;
EigenvalueType m_eivalues;
@@ -243,11 +256,14 @@ ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const Ma
// Do a complex Schur decomposition, A = U T U^*
// The eigenvalues are on the diagonal of T.
m_schur.compute(matrix, computeEigenvectors);
- m_eivalues = m_schur.matrixT().diagonal();
- if(computeEigenvectors)
- doComputeEigenvectors(matrix.norm());
- sortEigenvalues(computeEigenvectors);
+ if(m_schur.info() == Success)
+ {
+ m_eivalues = m_schur.matrixT().diagonal();
+ if(computeEigenvectors)
+ doComputeEigenvectors(matrix.norm());
+ sortEigenvalues(computeEigenvectors);
+ }
m_isInitialized = true;
m_eigenvectorsOk = computeEigenvectors;
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h
index cf6755bfc..6a8daebb8 100644
--- a/Eigen/src/Eigenvalues/ComplexSchur.h
+++ b/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -27,6 +27,9 @@
#ifndef EIGEN_COMPLEX_SCHUR_H
#define EIGEN_COMPLEX_SCHUR_H
+#include "./EigenvaluesCommon.h"
+#include "./HessenbergDecomposition.h"
+
template<typename MatrixType, bool IsComplex> struct ei_complex_schur_reduce_to_hessenberg;
/** \eigenvalues_module \ingroup Eigenvalues_Module
@@ -192,6 +195,16 @@ template<typename _MatrixType> class ComplexSchur
*/
ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
+ /** \brief Reports whether previous computation was successful.
+ *
+ * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
+ */
+ ComputationInfo info() const
+ {
+ ei_assert(m_isInitialized && "RealSchur is not initialized.");
+ return m_info;
+ }
+
/** \brief Maximum number of iterations.
*
* Maximum number of iterations allowed for an eigenvalue to converge.
@@ -201,6 +214,7 @@ template<typename _MatrixType> class ComplexSchur
protected:
ComplexMatrixType m_matT, m_matU;
HessenbergDecomposition<MatrixType> m_hess;
+ ComputationInfo m_info;
bool m_isInitialized;
bool m_matUisUptodate;
@@ -312,6 +326,7 @@ ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& ma
{
m_matT = matrix.template cast<ComplexScalar>();
if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
+ m_info = Success;
m_isInitialized = true;
m_matUisUptodate = computeU;
return *this;
@@ -382,7 +397,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
// if we spent too many iterations on the current element, we give up
iter++;
- if(iter >= m_maxIterations) break;
+ if(iter > m_maxIterations) break;
// find il, the top row of the active submatrix
il = iu-1;
@@ -412,12 +427,10 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
}
}
- if(iter >= m_maxIterations)
- {
- // FIXME : what to do when iter==MAXITER ??
- // std::cerr << "MAXITER" << std::endl;
- return;
- }
+ if(iter <= m_maxIterations)
+ m_info = Success;
+ else
+ m_info = NoConvergence;
m_isInitialized = true;
m_matUisUptodate = computeU;
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index f745413a8..04caee658 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -26,6 +26,7 @@
#ifndef EIGEN_EIGENSOLVER_H
#define EIGEN_EIGENSOLVER_H
+#include "./EigenvaluesCommon.h"
#include "./RealSchur.h"
/** \eigenvalues_module \ingroup Eigenvalues_Module
@@ -286,6 +287,12 @@ template<typename _MatrixType> class EigenSolver
*/
EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
+ ComputationInfo info() const
+ {
+ ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
+ return m_realSchur.info();
+ }
+
private:
void doComputeEigenvectors();
@@ -358,33 +365,36 @@ EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matr
// Reduce to real Schur form.
m_realSchur.compute(matrix, computeEigenvectors);
- m_matT = m_realSchur.matrixT();
- if (computeEigenvectors)
- m_eivec = m_realSchur.matrixU();
-
- // Compute eigenvalues from matT
- m_eivalues.resize(matrix.cols());
- Index i = 0;
- while (i < matrix.cols())
+ if (m_realSchur.info() == Success)
{
- if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
- {
- m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
- ++i;
- }
- else
+ m_matT = m_realSchur.matrixT();
+ if (computeEigenvectors)
+ m_eivec = m_realSchur.matrixU();
+
+ // Compute eigenvalues from matT
+ m_eivalues.resize(matrix.cols());
+ Index i = 0;
+ while (i < matrix.cols())
{
- Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
- Scalar z = ei_sqrt(ei_abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
- m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
- m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
- i += 2;
+ if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
+ {
+ m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
+ ++i;
+ }
+ else
+ {
+ Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
+ Scalar z = ei_sqrt(ei_abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
+ m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
+ m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
+ i += 2;
+ }
}
+
+ // Compute eigenvectors.
+ if (computeEigenvectors)
+ doComputeEigenvectors();
}
-
- // Compute eigenvectors.
- if (computeEigenvectors)
- doComputeEigenvectors();
m_isInitialized = true;
m_eigenvectorsOk = computeEigenvectors;
diff --git a/Eigen/src/Eigenvalues/EigenvaluesCommon.h b/Eigen/src/Eigenvalues/EigenvaluesCommon.h
new file mode 100644
index 000000000..d5fff9ba1
--- /dev/null
+++ b/Eigen/src/Eigenvalues/EigenvaluesCommon.h
@@ -0,0 +1,39 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_EIGENVALUES_COMMON_H
+#define EIGEN_EIGENVALUES_COMMON_H
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+ * \nonstableyet
+ *
+ * \brief Enum for reporting the status of a computation.
+ */
+enum ComputationInfo {
+ Success = 0, /**< \brief Computation was successful. */
+ NoConvergence = 1 /**< \brief Iterative procedure did not converge. */
+};
+
+#endif // EIGEN_EIGENVALUES_COMMON_H
+
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index 584b1d05e..156706573 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -26,6 +26,7 @@
#ifndef EIGEN_REAL_SCHUR_H
#define EIGEN_REAL_SCHUR_H
+#include "./EigenvaluesCommon.h"
#include "./HessenbergDecomposition.h"
/** \eigenvalues_module \ingroup Eigenvalues_Module
@@ -176,6 +177,16 @@ template<typename _MatrixType> class RealSchur
*/
RealSchur& compute(const MatrixType& matrix, bool computeU = true);
+ /** \brief Reports whether previous computation was successful.
+ *
+ * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
+ */
+ ComputationInfo info() const
+ {
+ ei_assert(m_isInitialized && "RealSchur is not initialized.");
+ return m_info;
+ }
+
/** \brief Maximum number of iterations.
*
* Maximum number of iterations allowed for an eigenvalue to converge.
@@ -188,6 +199,7 @@ template<typename _MatrixType> class RealSchur
MatrixType m_matU;
ColumnVectorType m_workspaceVector;
HessenbergDecomposition<MatrixType> m_hess;
+ ComputationInfo m_info;
bool m_isInitialized;
bool m_matUisUptodate;
@@ -249,20 +261,21 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix,
{
Vector3s firstHouseholderVector, shiftInfo;
computeShift(iu, iter, exshift, shiftInfo);
- iter = iter + 1; // (Could check iteration count here.)
- if (iter >= m_maxIterations) break;
+ iter = iter + 1;
+ if (iter > m_maxIterations) break;
Index im;
initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
}
}
- if(iter < m_maxIterations)
- {
- m_isInitialized = true;
- m_matUisUptodate = computeU;
- }
+ if(iter <= m_maxIterations)
+ m_info = Success;
+ else
+ m_info = NoConvergence;
+ m_isInitialized = true;
+ m_matUisUptodate = computeU;
return *this;
}
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index b08bbd7e2..6a7d46b39 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -26,6 +26,9 @@
#ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
#define EIGEN_SELFADJOINTEIGENSOLVER_H
+#include "./EigenvaluesCommon.h"
+#include "./Tridiagonalization.h"
+
/** \eigenvalues_module \ingroup Eigenvalues_Module
* \nonstableyet
*
@@ -360,6 +363,16 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
}
+ /** \brief Reports whether previous computation was successful.
+ *
+ * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
+ */
+ ComputationInfo info() const
+ {
+ ei_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
+ return m_info;
+ }
+
/** \brief Maximum number of iterations.
*
* Maximum number of iterations allowed for an eigenvalue to converge.
@@ -371,6 +384,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
RealVectorType m_eivalues;
TridiagonalizationType m_tridiag;
typename TridiagonalizationType::SubDiagonalType m_subdiag;
+ ComputationInfo m_info;
bool m_isInitialized;
bool m_eigenvectorsOk;
};
@@ -410,6 +424,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(
m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0));
if(computeEigenvectors)
m_eivec.setOnes();
+ m_info = Success;
m_isInitialized = true;
m_eigenvectorsOk = computeEigenvectors;
return *this;
@@ -443,7 +458,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(
// if we spent too many iterations on the current element, we give up
iter++;
- if(iter >= m_maxIterations) break;
+ if(iter > m_maxIterations) break;
start = end - 1;
while (start>0 && m_subdiag[start-1]!=0)
@@ -452,23 +467,26 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(
ei_tridiagonal_qr_step(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
}
- if(iter >= m_maxIterations)
- {
- return *this;
- }
+ if (iter <= m_maxIterations)
+ m_info = Success;
+ else
+ m_info = NoConvergence;
// Sort eigenvalues and corresponding vectors.
// TODO make the sort optional ?
// TODO use a better sort algorithm !!
- for (Index i = 0; i < n-1; ++i)
+ if (m_info == Success)
{
- Index k;
- m_eivalues.segment(i,n-i).minCoeff(&k);
- if (k > 0)
+ for (Index i = 0; i < n-1; ++i)
{
- std::swap(m_eivalues[i], m_eivalues[k+i]);
- if(computeEigenvectors)
- m_eivec.col(i).swap(m_eivec.col(k+i));
+ Index k;
+ m_eivalues.segment(i,n-i).minCoeff(&k);
+ if (k > 0)
+ {
+ std::swap(m_eivalues[i], m_eivalues[k+i]);
+ if(computeEigenvectors)
+ m_eivec.col(i).swap(m_eivec.col(k+i));
+ }
}
}
diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp
index 3285d26c2..dc3b2cfb0 100644
--- a/test/eigensolver_complex.cpp
+++ b/test/eigensolver_complex.cpp
@@ -24,6 +24,7 @@
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
+#include <limits>
#include <Eigen/Eigenvalues>
#include <Eigen/LU>
@@ -60,15 +61,18 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
MatrixType symmA = a.adjoint() * a;
ComplexEigenSolver<MatrixType> ei0(symmA);
+ VERIFY_IS_EQUAL(ei0.info(), Success);
VERIFY_IS_APPROX(symmA * ei0.eigenvectors(), ei0.eigenvectors() * ei0.eigenvalues().asDiagonal());
ComplexEigenSolver<MatrixType> ei1(a);
+ VERIFY_IS_EQUAL(ei1.info(), Success);
VERIFY_IS_APPROX(a * ei1.eigenvectors(), ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
// Note: If MatrixType is real then a.eigenvalues() uses EigenSolver and thus
// another algorithm so results may differ slightly
verify_is_approx_upto_permutation(a.eigenvalues(), ei1.eigenvalues());
ComplexEigenSolver<MatrixType> eiNoEivecs(a, false);
+ VERIFY_IS_EQUAL(eiNoEivecs.info(), Success);
VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
// Regression test for issue #66
@@ -78,6 +82,14 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
MatrixType id = MatrixType::Identity(rows, cols);
VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
+
+ if (rows > 1)
+ {
+ // Test matrix with NaN
+ a(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
+ ComplexEigenSolver<MatrixType> eiNaN(a);
+ VERIFY_IS_EQUAL(eiNaN.info(), NoConvergence);
+ }
}
template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp
index 79c08ec31..92741a35c 100644
--- a/test/eigensolver_generic.cpp
+++ b/test/eigensolver_generic.cpp
@@ -24,6 +24,7 @@
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
+#include <limits>
#include <Eigen/Eigenvalues>
#ifdef HAS_GSL
@@ -44,29 +45,38 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType;
typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
- // RealScalar largerEps = 10*test_precision<RealScalar>();
-
MatrixType a = MatrixType::Random(rows,cols);
MatrixType a1 = MatrixType::Random(rows,cols);
MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
EigenSolver<MatrixType> ei0(symmA);
+ VERIFY_IS_EQUAL(ei0.info(), Success);
VERIFY_IS_APPROX(symmA * ei0.pseudoEigenvectors(), ei0.pseudoEigenvectors() * ei0.pseudoEigenvalueMatrix());
VERIFY_IS_APPROX((symmA.template cast<Complex>()) * (ei0.pseudoEigenvectors().template cast<Complex>()),
(ei0.pseudoEigenvectors().template cast<Complex>()) * (ei0.eigenvalues().asDiagonal()));
EigenSolver<MatrixType> ei1(a);
+ VERIFY_IS_EQUAL(ei1.info(), Success);
VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues());
EigenSolver<MatrixType> eiNoEivecs(a, false);
+ VERIFY_IS_EQUAL(eiNoEivecs.info(), Success);
VERIFY_IS_APPROX(ei1.eigenvalues(), eiNoEivecs.eigenvalues());
VERIFY_IS_APPROX(ei1.pseudoEigenvalueMatrix(), eiNoEivecs.pseudoEigenvalueMatrix());
MatrixType id = MatrixType::Identity(rows, cols);
VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
+
+ if (rows > 2)
+ {
+ // Test matrix with NaN
+ a(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
+ EigenSolver<MatrixType> eiNaN(a);
+ VERIFY_IS_EQUAL(eiNaN.info(), NoConvergence);
+ }
}
template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)
diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index 3ff84c4e0..9f0c4cf38 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@@ -23,6 +24,7 @@
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
+#include <limits>
#include <Eigen/Eigenvalues>
#ifdef HAS_GSL
@@ -101,14 +103,17 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
}
#endif
+ VERIFY_IS_EQUAL(eiSymm.info(), Success);
VERIFY((symmA * eiSymm.eigenvectors()).isApprox(
eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps));
VERIFY_IS_APPROX(symmA.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues());
SelfAdjointEigenSolver<MatrixType> eiSymmNoEivecs(symmA, false);
+ VERIFY_IS_EQUAL(eiSymmNoEivecs.info(), Success);
VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues());
// generalized eigen problem Ax = lBx
+ VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox(
symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
@@ -120,6 +125,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1));
SelfAdjointEigenSolver<MatrixType> eiSymmUninitialized;
+ VERIFY_RAISES_ASSERT(eiSymmUninitialized.info());
VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvalues());
VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
@@ -129,6 +135,14 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
+
+ if (rows > 1)
+ {
+ // Test matrix with NaN
+ symmA(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
+ SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmA);
+ VERIFY_IS_EQUAL(eiSymmNaN.info(), NoConvergence);
+ }
}
void test_eigensolver_selfadjoint()
diff --git a/test/schur_complex.cpp b/test/schur_complex.cpp
index cc8174d00..67c41d41f 100644
--- a/test/schur_complex.cpp
+++ b/test/schur_complex.cpp
@@ -23,6 +23,7 @@
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
+#include <limits>
#include <Eigen/Eigenvalues>
template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTime)
@@ -34,6 +35,7 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
for(int counter = 0; counter < g_repeat; ++counter) {
MatrixType A = MatrixType::Random(size, size);
ComplexSchur<MatrixType> schurOfA(A);
+ VERIFY_IS_EQUAL(schurOfA.info(), Success);
ComplexMatrixType U = schurOfA.matrixU();
ComplexMatrixType T = schurOfA.matrixT();
for(int row = 1; row < size; ++row) {
@@ -48,19 +50,28 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
ComplexSchur<MatrixType> csUninitialized;
VERIFY_RAISES_ASSERT(csUninitialized.matrixT());
VERIFY_RAISES_ASSERT(csUninitialized.matrixU());
+ VERIFY_RAISES_ASSERT(csUninitialized.info());
// Test whether compute() and constructor returns same result
MatrixType A = MatrixType::Random(size, size);
ComplexSchur<MatrixType> cs1;
cs1.compute(A);
ComplexSchur<MatrixType> cs2(A);
+ VERIFY_IS_EQUAL(cs1.info(), Success);
+ VERIFY_IS_EQUAL(cs2.info(), Success);
VERIFY_IS_EQUAL(cs1.matrixT(), cs2.matrixT());
VERIFY_IS_EQUAL(cs1.matrixU(), cs2.matrixU());
// Test computation of only T, not U
ComplexSchur<MatrixType> csOnlyT(A, false);
+ VERIFY_IS_EQUAL(csOnlyT.info(), Success);
VERIFY_IS_EQUAL(cs1.matrixT(), csOnlyT.matrixT());
VERIFY_RAISES_ASSERT(csOnlyT.matrixU());
+
+ // Test matrix with NaN
+ A(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
+ ComplexSchur<MatrixType> csNaN(A);
+ VERIFY_IS_EQUAL(csNaN.info(), NoConvergence);
}
void test_schur_complex()
diff --git a/test/schur_real.cpp b/test/schur_real.cpp
index 116c8dbce..1e8c4b0ba 100644
--- a/test/schur_real.cpp
+++ b/test/schur_real.cpp
@@ -23,6 +23,7 @@
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#include "main.h"
+#include <limits>
#include <Eigen/Eigenvalues>
template<typename MatrixType> void verifyIsQuasiTriangular(const MatrixType& T)
@@ -55,6 +56,7 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
for(int counter = 0; counter < g_repeat; ++counter) {
MatrixType A = MatrixType::Random(size, size);
RealSchur<MatrixType> schurOfA(A);
+ VERIFY_IS_EQUAL(schurOfA.info(), Success);
MatrixType U = schurOfA.matrixU();
MatrixType T = schurOfA.matrixT();
verifyIsQuasiTriangular(T);
@@ -65,19 +67,28 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim
RealSchur<MatrixType> rsUninitialized;
VERIFY_RAISES_ASSERT(rsUninitialized.matrixT());
VERIFY_RAISES_ASSERT(rsUninitialized.matrixU());
+ VERIFY_RAISES_ASSERT(rsUninitialized.info());
// Test whether compute() and constructor returns same result
MatrixType A = MatrixType::Random(size, size);
RealSchur<MatrixType> rs1;
rs1.compute(A);
RealSchur<MatrixType> rs2(A);
+ VERIFY_IS_EQUAL(rs1.info(), Success);
+ VERIFY_IS_EQUAL(rs2.info(), Success);
VERIFY_IS_EQUAL(rs1.matrixT(), rs2.matrixT());
VERIFY_IS_EQUAL(rs1.matrixU(), rs2.matrixU());
// Test computation of only T, not U
RealSchur<MatrixType> rsOnlyT(A, false);
+ VERIFY_IS_EQUAL(rsOnlyT.info(), Success);
VERIFY_IS_EQUAL(rs1.matrixT(), rsOnlyT.matrixT());
VERIFY_RAISES_ASSERT(rsOnlyT.matrixU());
+
+ // Test matrix with NaN
+ A(0,0) = std::numeric_limits<typename MatrixType::Scalar>::quiet_NaN();
+ RealSchur<MatrixType> rsNaN(A);
+ VERIFY_IS_EQUAL(rsNaN.info(), NoConvergence);
}
void test_schur_real()
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-04 21:09:59 +0000