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
14namespace Eigen {
15
16namespace internal {
17
18template<typename Derived, bool IsComplex>
19struct eigenvalues_selector
20{
21  // this is the implementation for the case IsComplex = true
22  static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
23  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
31template<typename Derived>
32struct 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  */
65template<typename Derived>
66inline typename MatrixBase<Derived>::EigenvaluesReturnType
67MatrixBase<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  */
87template<typename MatrixType, unsigned int UpLo>
88inline typename SelfAdjointView<MatrixType, UpLo>::EigenvaluesReturnType
89SelfAdjointView<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  */
120template<typename Derived>
121inline typename MatrixBase<Derived>::RealScalar
122MatrixBase<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  */
151template<typename MatrixType, unsigned int UpLo>
152inline typename SelfAdjointView<MatrixType, UpLo>::RealScalar
153SelfAdjointView<MatrixType, UpLo>::operatorNorm() const
154{
155  return eigenvalues().cwiseAbs().maxCoeff();
156}
157
158} // end namespace Eigen
159
160#endif
161