1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_MATRIXBASEEIGENVALUES_H
12 #define EIGEN_MATRIXBASEEIGENVALUES_H
13 
14 namespace Eigen {
15 
16 namespace internal {
17 
18 template<typename Derived, bool IsComplex>
19 struct eigenvalues_selector
20 {
21   // this is the implementation for the case IsComplex = true
22   static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
runeigenvalues_selector23   run(const MatrixBase<Derived>& m)
24   {
25     typedef typename Derived::PlainObject PlainObject;
26     PlainObject m_eval(m);
27     return ComplexEigenSolver<PlainObject>(m_eval, false).eigenvalues();
28   }
29 };
30 
31 template<typename Derived>
32 struct eigenvalues_selector<Derived, false>
33 {
34   static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
35   run(const MatrixBase<Derived>& m)
36   {
37     typedef typename Derived::PlainObject PlainObject;
38     PlainObject m_eval(m);
39     return EigenSolver<PlainObject>(m_eval, false).eigenvalues();
40   }
41 };
42 
43 } // end namespace internal
44 
45 /** \brief Computes the eigenvalues of a matrix
46   * \returns Column vector containing the eigenvalues.
47   *
48   * \eigenvalues_module
49   * This function computes the eigenvalues with the help of the EigenSolver
50   * class (for real matrices) or the ComplexEigenSolver class (for complex
51   * matrices).
52   *
53   * The eigenvalues are repeated according to their algebraic multiplicity,
54   * so there are as many eigenvalues as rows in the matrix.
55   *
56   * The SelfAdjointView class provides a better algorithm for selfadjoint
57   * matrices.
58   *
59   * Example: \include MatrixBase_eigenvalues.cpp
60   * Output: \verbinclude MatrixBase_eigenvalues.out
61   *
62   * \sa EigenSolver::eigenvalues(), ComplexEigenSolver::eigenvalues(),
63   *     SelfAdjointView::eigenvalues()
64   */
65 template<typename Derived>
66 inline typename MatrixBase<Derived>::EigenvaluesReturnType
67 MatrixBase<Derived>::eigenvalues() const
68 {
69   typedef typename internal::traits<Derived>::Scalar Scalar;
70   return internal::eigenvalues_selector<Derived, NumTraits<Scalar>::IsComplex>::run(derived());
71 }
72 
73 /** \brief Computes the eigenvalues of a matrix
74   * \returns Column vector containing the eigenvalues.
75   *
76   * \eigenvalues_module
77   * This function computes the eigenvalues with the help of the
78   * SelfAdjointEigenSolver class.  The eigenvalues are repeated according to
79   * their algebraic multiplicity, so there are as many eigenvalues as rows in
80   * the matrix.
81   *
82   * Example: \include SelfAdjointView_eigenvalues.cpp
83   * Output: \verbinclude SelfAdjointView_eigenvalues.out
84   *
85   * \sa SelfAdjointEigenSolver::eigenvalues(), MatrixBase::eigenvalues()
86   */
87 template<typename MatrixType, unsigned int UpLo>
88 inline typename SelfAdjointView<MatrixType, UpLo>::EigenvaluesReturnType
89 SelfAdjointView<MatrixType, UpLo>::eigenvalues() const
90 {
91   typedef typename SelfAdjointView<MatrixType, UpLo>::PlainObject PlainObject;
92   PlainObject thisAsMatrix(*this);
93   return SelfAdjointEigenSolver<PlainObject>(thisAsMatrix, false).eigenvalues();
94 }
95 
96 
97 
98 /** \brief Computes the L2 operator norm
99   * \returns Operator norm of the matrix.
100   *
101   * \eigenvalues_module
102   * This function computes the L2 operator norm of a matrix, which is also
103   * known as the spectral norm. The norm of a matrix \f$ A \f$ is defined to be
104   * \f[ \|A\|_2 = \max_x \frac{\|Ax\|_2}{\|x\|_2} \f]
105   * where the maximum is over all vectors and the norm on the right is the
106   * Euclidean vector norm. The norm equals the largest singular value, which is
107   * the square root of the largest eigenvalue of the positive semi-definite
108   * matrix \f$ A^*A \f$.
109   *
110   * The current implementation uses the eigenvalues of \f$ A^*A \f$, as computed
111   * by SelfAdjointView::eigenvalues(), to compute the operator norm of a
112   * matrix.  The SelfAdjointView class provides a better algorithm for
113   * selfadjoint matrices.
114   *
115   * Example: \include MatrixBase_operatorNorm.cpp
116   * Output: \verbinclude MatrixBase_operatorNorm.out
117   *
118   * \sa SelfAdjointView::eigenvalues(), SelfAdjointView::operatorNorm()
119   */
120 template<typename Derived>
121 inline typename MatrixBase<Derived>::RealScalar
122 MatrixBase<Derived>::operatorNorm() const
123 {
124   using std::sqrt;
125   typename Derived::PlainObject m_eval(derived());
126   // FIXME if it is really guaranteed that the eigenvalues are already sorted,
127   // then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
128   return sqrt((m_eval*m_eval.adjoint())
129                  .eval()
130 		 .template selfadjointView<Lower>()
131 		 .eigenvalues()
132 		 .maxCoeff()
133 		 );
134 }
135 
136 /** \brief Computes the L2 operator norm
137   * \returns Operator norm of the matrix.
138   *
139   * \eigenvalues_module
140   * This function computes the L2 operator norm of a self-adjoint matrix. For a
141   * self-adjoint matrix, the operator norm is the largest eigenvalue.
142   *
143   * The current implementation uses the eigenvalues of the matrix, as computed
144   * by eigenvalues(), to compute the operator norm of the matrix.
145   *
146   * Example: \include SelfAdjointView_operatorNorm.cpp
147   * Output: \verbinclude SelfAdjointView_operatorNorm.out
148   *
149   * \sa eigenvalues(), MatrixBase::operatorNorm()
150   */
151 template<typename MatrixType, unsigned int UpLo>
152 inline typename SelfAdjointView<MatrixType, UpLo>::RealScalar
153 SelfAdjointView<MatrixType, UpLo>::operatorNorm() const
154 {
155   return eigenvalues().cwiseAbs().maxCoeff();
156 }
157 
158 } // end namespace Eigen
159 
160 #endif
161