Add info() method which can be queried to check whether iteration converged.

5 files changed, 60 insertions, 2 deletions

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()