1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_MATRIXBASEEIGENVALUES_H
#define EIGEN_MATRIXBASEEIGENVALUES_H
namespace Eigen {
namespace internal {
template<typename Derived, bool IsComplex>
struct 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, false).eigenvalues();
}
};
template<typename Derived>
struct 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, false).eigenvalues();
}
};
} // end namespace internal
/** \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 internal::traits<Derived>::Scalar Scalar;
return internal::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, false).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
{
using std::sqrt;
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 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();
}
} // end namespace Eigen
#endif
|
