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

Eigen/Eigenvalues

@@ -42,6 +42,7 @@ namespace Eigen {
 #include "src/Eigenvalues/HessenbergDecomposition.h"
 #include "src/Eigenvalues/ComplexSchur.h"
 #include "src/Eigenvalues/ComplexEigenSolver.h"
+#include "src/Eigenvalues/MatrixBaseEigenvalues.h"
 
 // declare all classes for a given matrix type
 #define EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \

Eigen/src/Core/MatrixBase.h

@@ -136,8 +136,8 @@ template<typename Derived> class MatrixBase
                         CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, Eigen::Transpose<Derived> >,
                         Transpose<Derived>
                      >::ret AdjointReturnType;
-    /** \internal the return type of MatrixBase::eigenvalues() */
-    typedef Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1> EigenvaluesReturnType;
+    /** \internal Return type of eigenvalues() */
+    typedef Matrix<std::complex<RealScalar>, ei_traits<Derived>::ColsAtCompileTime, 1> EigenvaluesReturnType;
     /** \internal the return type of identity */
     typedef CwiseNullaryOp<ei_scalar_identity_op<Scalar>,Derived> IdentityReturnType;
     /** \internal the return type of unit vectors */

Eigen/src/Core/SelfAdjointView.h

@@ -153,6 +153,16 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
     const LLT<PlainObject, UpLo> llt() const;
     const LDLT<PlainObject> ldlt() const;
 
+/////////// Eigenvalue module ///////////
+
+    /** Real part of #Scalar */
+    typedef typename NumTraits<Scalar>::Real RealScalar;
+    /** Return type of eigenvalues() */
+    typedef Matrix<RealScalar, ei_traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType;
+
+    EigenvaluesReturnType eigenvalues() const;
+    RealScalar operatorNorm() const;
+
   protected:
     const typename MatrixType::Nested m_matrix;
 };

Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h

@@ -0,0 +1,168 @@
+// This file is part of Eigen, a lightweight C++ template library
+// 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
+// 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_MATRIXBASEEIGENVALUES_H
+#define EIGEN_MATRIXBASEEIGENVALUES_H
+
+
+
+template<typename Derived, bool IsComplex>
+struct ei_eigenvalues_selector
+{
+  // this is the implementation for the case IsComplex = true
+  static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
+  run(const MatrixBase<Derived>& m)
+  {
+    typedef typename Derived::PlainObject PlainObject;
+    PlainObject m_eval(m);
+    return ComplexEigenSolver<PlainObject>(m_eval).eigenvalues();
+  }
+};
+
+template<typename Derived>
+struct ei_eigenvalues_selector<Derived, false>
+{
+  static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
+  run(const MatrixBase<Derived>& m)
+  {
+    typedef typename Derived::PlainObject PlainObject;
+    PlainObject m_eval(m);
+    return EigenSolver<PlainObject>(m_eval).eigenvalues();
+  }
+};
+
+/** \brief Computes the eigenvalues of a matrix 
+  * \returns Column vector containing the eigenvalues.
+  *
+  * \eigenvalues_module
+  * This function computes the eigenvalues with the help of the EigenSolver
+  * class (for real matrices) or the ComplexEigenSolver class (for complex
+  * matrices). 
+  *
+  * The eigenvalues are repeated according to their algebraic multiplicity,
+  * so there are as many eigenvalues as rows in the matrix.
+  *
+  * The SelfAdjointView class provides a better algorithm for selfadjoint
+  * matrices.
+  *
+  * Example: \include MatrixBase_eigenvalues.cpp
+  * Output: \verbinclude MatrixBase_eigenvalues.out
+  *
+  * \sa EigenSolver::eigenvalues(), ComplexEigenSolver::eigenvalues(),
+  *     SelfAdjointView::eigenvalues()
+  */
+template<typename Derived>
+inline typename MatrixBase<Derived>::EigenvaluesReturnType
+MatrixBase<Derived>::eigenvalues() const
+{
+  typedef typename ei_traits<Derived>::Scalar Scalar;
+  return ei_eigenvalues_selector<Derived, NumTraits<Scalar>::IsComplex>::run(derived());
+}
+
+/** \brief Computes the eigenvalues of a matrix
+  * \returns Column vector containing the eigenvalues.
+  *
+  * \eigenvalues_module
+  * This function computes the eigenvalues with the help of the
+  * SelfAdjointEigenSolver class.  The eigenvalues are repeated according to
+  * their algebraic multiplicity, so there are as many eigenvalues as rows in
+  * the matrix.
+  *
+  * Example: \include SelfAdjointView_eigenvalues.cpp
+  * Output: \verbinclude SelfAdjointView_eigenvalues.out
+  *
+  * \sa SelfAdjointEigenSolver::eigenvalues(), MatrixBase::eigenvalues()
+  */
+template<typename MatrixType, unsigned int UpLo> 
+inline typename SelfAdjointView<MatrixType, UpLo>::EigenvaluesReturnType
+SelfAdjointView<MatrixType, UpLo>::eigenvalues() const
+{
+  typedef typename SelfAdjointView<MatrixType, UpLo>::PlainObject PlainObject;
+  PlainObject thisAsMatrix(*this);
+  return SelfAdjointEigenSolver<PlainObject>(thisAsMatrix).eigenvalues();
+}
+
+
+
+/** \brief Computes the L2 operator norm
+  * \returns Operator norm of the matrix.
+  *
+  * \eigenvalues_module
+  * This function computes the L2 operator norm of a matrix, which is also
+  * known as the spectral norm. The norm of a matrix \f$ A \f$ is defined to be
+  * \f[ \|A\|_2 = \max_x \frac{\|Ax\|_2}{\|x\|_2} \f]
+  * where the maximum is over all vectors and the norm on the right is the
+  * Euclidean vector norm. The norm equals the largest singular value, which is
+  * the square root of the largest eigenvalue of the positive semi-definite
+  * matrix \f$ A^*A \f$.
+  *
+  * The current implementation uses the eigenvalues of \f$ A^*A \f$, as computed
+  * by SelfAdjointView::eigenvalues(), to compute the operator norm of a
+  * matrix.  The SelfAdjointView class provides a better algorithm for
+  * selfadjoint matrices.
+  *
+  * Example: \include MatrixBase_operatorNorm.cpp
+  * Output: \verbinclude MatrixBase_operatorNorm.out
+  *
+  * \sa SelfAdjointView::eigenvalues(), SelfAdjointView::operatorNorm()
+  */
+template<typename Derived>
+inline typename MatrixBase<Derived>::RealScalar
+MatrixBase<Derived>::operatorNorm() const
+{
+  typename Derived::PlainObject m_eval(derived());
+  // FIXME if it is really guaranteed that the eigenvalues are already sorted,
+  // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
+  return ei_sqrt((m_eval*m_eval.adjoint())
+                 .eval()
+		 .template selfadjointView<Lower>()
+		 .eigenvalues()
+		 .maxCoeff()
+		 );
+}
+
+/** \brief Computes the L2 operator norm
+  * \returns Operator norm of the matrix.
+  *
+  * \eigenvalues_module
+  * This function computes the L2 operator norm of a self-adjoint matrix. For a
+  * self-adjoint matrix, the operator norm is the largest eigenvalue.
+  *
+  * The current implementation uses the eigenvalues of the matrix, as computed
+  * by eigenvalues(), to compute the operator norm of the matrix.
+  *
+  * Example: \include SelfAdjointView_operatorNorm.cpp
+  * Output: \verbinclude SelfAdjointView_operatorNorm.out
+  *
+  * \sa eigenvalues(), MatrixBase::operatorNorm()
+  */
+template<typename MatrixType, unsigned int UpLo>
+inline typename SelfAdjointView<MatrixType, UpLo>::RealScalar
+SelfAdjointView<MatrixType, UpLo>::operatorNorm() const
+{
+  return eigenvalues().cwiseAbs().maxCoeff();
+}
+
+#endif

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

@@ -481,59 +481,6 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors
 
 #endif // EIGEN_HIDE_HEAVY_CODE
 
-/** \eigenvalues_module
-  *
-  * \returns a vector listing the eigenvalues of this matrix.
-  */
-template<typename Derived>
-inline Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1>
-MatrixBase<Derived>::eigenvalues() const
-{
-  ei_assert(Flags&SelfAdjoint);
-  return SelfAdjointEigenSolver<typename Derived::PlainObject>(eval(),false).eigenvalues();
-}
-
-template<typename Derived, bool IsSelfAdjoint>
-struct ei_operatorNorm_selector
-{
-  static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
-  operatorNorm(const MatrixBase<Derived>& m)
-  {
-    // FIXME if it is really guaranteed that the eigenvalues are already sorted,
-    // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
-    return m.eigenvalues().cwiseAbs().maxCoeff();
-  }
-};
-
-template<typename Derived> struct ei_operatorNorm_selector<Derived, false>
-{
-  static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
-  operatorNorm(const MatrixBase<Derived>& m)
-  {
-    typename Derived::PlainObject m_eval(m);
-    // FIXME if it is really guaranteed that the eigenvalues are already sorted,
-    // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
-    return ei_sqrt(
-             (m_eval*m_eval.adjoint())
-             .template marked<SelfAdjoint>()
-             .eigenvalues()
-             .maxCoeff()
-           );
-  }
-};
-
-/** \eigenvalues_module
-  *
-  * \returns the matrix norm of this matrix.
-  */
-template<typename Derived>
-inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
-MatrixBase<Derived>::operatorNorm() const
-{
-  return ei_operatorNorm_selector<Derived, Flags&SelfAdjoint>
-       ::operatorNorm(derived());
-}
-
 #ifndef EIGEN_EXTERN_INSTANTIATIONS
 template<typename RealScalar, typename Scalar>
 static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n)

doc/snippets/MatrixBase_eigenvalues.cpp

@@ -0,0 +1,3 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+VectorXcd eivals = ones.eigenvalues();
+cout << "The eigenvalues of the 3x3 matrix of ones are:" << endl << eivals << endl;

doc/snippets/MatrixBase_operatorNorm.cpp

@@ -0,0 +1,3 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+cout << "The operator norm of the 3x3 matrix of ones is "
+     << ones.operatorNorm() << endl;

doc/snippets/SelfAdjointView_eigenvalues.cpp

@@ -0,0 +1,3 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+VectorXd eivals = ones.selfadjointView<Lower>().eigenvalues();
+cout << "The eigenvalues of the 3x3 matrix of ones are:" << endl << eivals << endl;

doc/snippets/SelfAdjointView_operatorNorm.cpp

@@ -0,0 +1,3 @@
+MatrixXd ones = MatrixXd::Ones(3,3);
+cout << "The operator norm of the 3x3 matrix of ones is "
+     << ones.selfadjointView<Lower>().operatorNorm() << endl;

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

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

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