Update eigenvalues() and operatorNorm() methods in MatrixBase.

* use SelfAdjointView instead of Eigen2's SelfAdjoint flag.
* add tests and documentation.
* allow eigenvalues() for non-selfadjoint matrices.
* they no longer depend only on SelfAdjointEigenSolver, so move them to
  a separate file

3 files changed, 29 insertions, 1 deletions

diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp
index b3d9ac24b..5c5d7b38f 100644
--- a/test/eigensolver_complex.cpp
+++ b/test/eigensolver_complex.cpp
@@ -26,6 +26,21 @@
 #include <Eigen/Eigenvalues>
 #include <Eigen/LU>
 
+/* Check that two column vectors are approximately equal upto permutations,
+   by checking that the k-th power sums are equal for k = 1, ..., vec1.rows() */
+template<typename VectorType>
+void verify_is_approx_upto_permutation(const VectorType& vec1, const VectorType& vec2)
+{
+  VERIFY(vec1.cols() == 1);
+  VERIFY(vec2.cols() == 1);
+  VERIFY(vec1.rows() == vec2.rows());
+  for (int k = 1; k <= vec1.rows(); ++k)
+  {
+    VERIFY_IS_APPROX(vec1.array().pow(k).sum(), vec2.array().pow(k).sum());
+  }
+}
+
+
 template<typename MatrixType> void eigensolver(const MatrixType& m)
 {
   /* this test covers the following files:
@@ -48,11 +63,17 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
 
   ComplexEigenSolver<MatrixType> ei1(a);
   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());
+  
   // Regression test for issue #66
   MatrixType z = MatrixType::Zero(rows,cols);
   ComplexEigenSolver<MatrixType> eiz(z);
   VERIFY((eiz.eigenvalues().cwiseEqual(0)).all());
+
+  MatrixType id = MatrixType::Identity(rows, cols);
+  VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
 }
 
 void test_eigensolver_complex()
diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp
index f24c3b4ed..d70f37ea4 100644
--- a/test/eigensolver_generic.cpp
+++ b/test/eigensolver_generic.cpp
@@ -58,7 +58,10 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
   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());
 
+  MatrixType id = MatrixType::Identity(rows, cols);
+  VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1));
 }
 
 template<typename MatrixType> void eigensolver_verify_assert()
diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index 70b3e6791..25ef280a1 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -103,6 +103,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
 
   VERIFY((symmA * eiSymm.eigenvectors()).isApprox(
           eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps));
+  VERIFY_IS_APPROX(symmA.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues());
 
   // generalized eigen problem Ax = lBx
   VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox(
@@ -111,6 +112,9 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   MatrixType sqrtSymmA = eiSymm.operatorSqrt();
   VERIFY_IS_APPROX(symmA, sqrtSymmA*sqrtSymmA);
   VERIFY_IS_APPROX(sqrtSymmA, symmA*eiSymm.operatorInverseSqrt());
+
+  MatrixType id = MatrixType::Identity(rows, cols);
+  VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1));
 }
 
 void test_eigensolver_selfadjoint()