From 9178e2bd54f64febb43025b9710387d2e98fea34 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Thu, 3 Jun 2010 22:59:57 +0100 Subject: [PATCH] Add info() method which can be queried to check whether iteration converged. --- Eigen/src/Eigenvalues/ComplexEigenSolver.h | 24 ++++++-- Eigen/src/Eigenvalues/ComplexSchur.h | 27 ++++++--- Eigen/src/Eigenvalues/EigenSolver.h | 60 +++++++++++-------- Eigen/src/Eigenvalues/EigenvaluesCommon.h | 39 ++++++++++++ Eigen/src/Eigenvalues/RealSchur.h | 27 ++++++--- .../src/Eigenvalues/SelfAdjointEigenSolver.h | 42 +++++++++---- test/eigensolver_complex.cpp | 12 ++++ test/eigensolver_generic.cpp | 14 ++++- test/eigensolver_selfadjoint.cpp | 14 +++++ test/schur_complex.cpp | 11 ++++ test/schur_real.cpp | 11 ++++ 11 files changed, 224 insertions(+), 57 deletions(-) create mode 100644 Eigen/src/Eigenvalues/EigenvaluesCommon.h 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 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& ComplexEigenSolver::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 struct ei_complex_schur_reduce_to_hessenberg; /** \eigenvalues_module \ingroup Eigenvalues_Module @@ -192,6 +195,16 @@ template 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 class ComplexSchur protected: ComplexMatrixType m_matT, m_matU; HessenbergDecomposition m_hess; + ComputationInfo m_info; bool m_isInitialized; bool m_matUisUptodate; @@ -312,6 +326,7 @@ ComplexSchur& ComplexSchur::compute(const MatrixType& ma { m_matT = matrix.template cast(); 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::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::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 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& EigenSolver::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 - { - 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; - } - } + m_matT = m_realSchur.matrixT(); + if (computeEigenvectors) + m_eivec = m_realSchur.matrixU(); - // Compute eigenvectors. - if (computeEigenvectors) - doComputeEigenvectors(); + // Compute eigenvalues from matT + m_eivalues.resize(matrix.cols()); + Index i = 0; + while (i < matrix.cols()) + { + 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(); + } 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 +// +// 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 . + +#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 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 class RealSchur MatrixType m_matU; ColumnVectorType m_workspaceVector; HessenbergDecomposition m_hess; + ComputationInfo m_info; bool m_isInitialized; bool m_matUisUptodate; @@ -249,20 +261,21 @@ RealSchur& RealSchur::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 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 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& SelfAdjointEigenSolver::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& SelfAdjointEigenSolver::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& SelfAdjointEigenSolver::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 . #include "main.h" +#include #include #include @@ -60,15 +61,18 @@ template void eigensolver(const MatrixType& m) MatrixType symmA = a.adjoint() * a; ComplexEigenSolver ei0(symmA); + VERIFY_IS_EQUAL(ei0.info(), Success); VERIFY_IS_APPROX(symmA * ei0.eigenvectors(), ei0.eigenvectors() * ei0.eigenvalues().asDiagonal()); ComplexEigenSolver 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 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 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::quiet_NaN(); + ComplexEigenSolver eiNaN(a); + VERIFY_IS_EQUAL(eiNaN.info(), NoConvergence); + } } template 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 . #include "main.h" +#include #include #ifdef HAS_GSL @@ -44,29 +45,38 @@ template void eigensolver(const MatrixType& m) typedef Matrix RealVectorType; typedef typename std::complex::Real> Complex; - // RealScalar largerEps = 10*test_precision(); - MatrixType a = MatrixType::Random(rows,cols); MatrixType a1 = MatrixType::Random(rows,cols); MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1; EigenSolver ei0(symmA); + VERIFY_IS_EQUAL(ei0.info(), Success); VERIFY_IS_APPROX(symmA * ei0.pseudoEigenvectors(), ei0.pseudoEigenvectors() * ei0.pseudoEigenvalueMatrix()); VERIFY_IS_APPROX((symmA.template cast()) * (ei0.pseudoEigenvectors().template cast()), (ei0.pseudoEigenvectors().template cast()) * (ei0.eigenvalues().asDiagonal())); EigenSolver ei1(a); + VERIFY_IS_EQUAL(ei1.info(), Success); VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix()); VERIFY_IS_APPROX(a.template cast() * ei1.eigenvectors(), ei1.eigenvectors() * ei1.eigenvalues().asDiagonal()); VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues()); EigenSolver 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::quiet_NaN(); + EigenSolver eiNaN(a); + VERIFY_IS_EQUAL(eiNaN.info(), NoConvergence); + } } template 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 +// Copyright (C) 2010 Jitse Niesen // // 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 . #include "main.h" +#include #include #ifdef HAS_GSL @@ -101,14 +103,17 @@ template 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().eigenvalues(), eiSymm.eigenvalues()); SelfAdjointEigenSolver 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 void selfadjointeigensolver(const MatrixType& m) VERIFY_IS_APPROX(id.template selfadjointView().operatorNorm(), RealScalar(1)); SelfAdjointEigenSolver 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 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::quiet_NaN(); + SelfAdjointEigenSolver 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 . #include "main.h" +#include #include template void schur(int size = MatrixType::ColsAtCompileTime) @@ -34,6 +35,7 @@ template void schur(int size = MatrixType::ColsAtCompileTim for(int counter = 0; counter < g_repeat; ++counter) { MatrixType A = MatrixType::Random(size, size); ComplexSchur 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 void schur(int size = MatrixType::ColsAtCompileTim ComplexSchur 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 cs1; cs1.compute(A); ComplexSchur 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 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::quiet_NaN(); + ComplexSchur 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 . #include "main.h" +#include #include template void verifyIsQuasiTriangular(const MatrixType& T) @@ -55,6 +56,7 @@ template void schur(int size = MatrixType::ColsAtCompileTim for(int counter = 0; counter < g_repeat; ++counter) { MatrixType A = MatrixType::Random(size, size); RealSchur schurOfA(A); + VERIFY_IS_EQUAL(schurOfA.info(), Success); MatrixType U = schurOfA.matrixU(); MatrixType T = schurOfA.matrixT(); verifyIsQuasiTriangular(T); @@ -65,19 +67,28 @@ template void schur(int size = MatrixType::ColsAtCompileTim RealSchur 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 rs1; rs1.compute(A); RealSchur 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 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::quiet_NaN(); + RealSchur rsNaN(A); + VERIFY_IS_EQUAL(rsNaN.info(), NoConvergence); } void test_schur_real()