From 0444e3601a598b97bb20e45636e2251ae76ee51b Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Mon, 2 Jun 2008 04:42:45 +0000
Subject: [PATCH] - add MatrixBase::eigenvalues() convenience method - add
 MatrixBase::matrixNorm(); in the non-selfadjoint case, we reduce to the  
 selfadjoint case by using the "C*-identity" a.k.a. 	norm of x = sqrt(norm
 of x * x.adjoint())

---
 Eigen/src/Core/MatrixBase.h           |  4 +++
 Eigen/src/LU/Inverse.h                |  4 +--
 Eigen/src/QR/SelfAdjointEigenSolver.h | 48 ++++++++++++++++++++++++++-
 3 files changed, 53 insertions(+), 3 deletions(-)

diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index a29504b33..e0220fd19 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -202,6 +202,7 @@ template<typename Derived> class MatrixBase : public ArrayBase<Derived>
     /** the return type of MatrixBase::adjoint() */
     typedef Transpose<NestByValue<typename ei_unref<ConjugateReturnType>::type> >
             AdjointReturnType;
+    typedef Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1> EigenvaluesReturnType;
     //@}
 
     /// \name Copying and initialization
@@ -605,6 +606,9 @@ template<typename Derived> class MatrixBase : public ArrayBase<Derived>
       */
     //@{
     const QR<typename ei_eval<Derived>::type> qr() const;
+    
+    EigenvaluesReturnType eigenvalues() const;
+    RealScalar matrixNorm() const;
     //@}
 
 };
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h
index 89140e8ff..8ff723e05 100644
--- a/Eigen/src/LU/Inverse.h
+++ b/Eigen/src/LU/Inverse.h
@@ -269,7 +269,7 @@ template<typename Derived>
 const Inverse<typename ei_eval<Derived>::type, true>
 MatrixBase<Derived>::inverse() const
 {
-  return Inverse<typename ei_eval<Derived>::type, true>(eval());
+  return Inverse<typename Derived::Eval, true>(eval());
 }
 
 /** \return the matrix inverse of \c *this, which is assumed to exist.
@@ -283,7 +283,7 @@ template<typename Derived>
 const Inverse<typename ei_eval<Derived>::type, false>
 MatrixBase<Derived>::quickInverse() const
 {
-  return Inverse<typename ei_eval<Derived>::type, false>(eval());
+  return Inverse<typename Derived::Eval, false>(eval());
 }
 
 #endif // EIGEN_INVERSE_H
diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h
index 0df0db0b7..e62a95356 100644
--- a/Eigen/src/QR/SelfAdjointEigenSolver.h
+++ b/Eigen/src/QR/SelfAdjointEigenSolver.h
@@ -31,6 +31,8 @@
   *
   * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
   *
+  * \note MatrixType must be an actual Matrix type, it can't be an expression type.
+  *
   * \sa MatrixBase::eigenvalues(), class EigenSolver
   */
 template<typename _MatrixType> class SelfAdjointEigenSolver
@@ -58,7 +60,6 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
 
     RealVectorType eigenvalues(void) const { return m_eivalues; }
 
-
   protected:
     MatrixType m_eivec;
     RealVectorType m_eivalues;
@@ -169,4 +170,49 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix)
   std::cout << "ei values = " << m_eivalues.transpose() << "\n\n";
 }
 
+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&SelfAdjointBit);
+  return SelfAdjointEigenSolver<typename Derived::Eval>(eval()).eigenvalues();
+}
+
+template<typename Derived, bool IsSelfAdjoint>
+struct ei_matrixNorm_selector
+{
+  static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
+  matrixNorm(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_matrixNorm_selector<Derived, false>
+{
+  static inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
+  matrixNorm(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 ei_sqrt(
+             (m*m.adjoint())
+             .template marked<SelfAdjoint>()
+             .eigenvalues()
+             .cwiseAbs()
+             .maxCoeff()
+           );
+  }
+};
+
+template<typename Derived>
+inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
+MatrixBase<Derived>::matrixNorm() const
+{
+  return ei_matrixNorm_selector<Derived, Flags&SelfAdjointBit>
+       ::matrixNorm(derived());
+}
+
 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H