← Back to branch summary

~siretart/ubuntu/utopic/blender/libav10

« back to all changes in this revision

Viewing changes to extern/Eigen3/Eigen/src/Eigenvalues/ComplexEigenSolver.h

  • Committer: Package Import Robot
  • Author(s): Matteo F. Vescovi
  • Date: 2012-07-23 08:54:18 UTC
  • mfrom: (14.2.16 sid)
  • mto: (14.2.19 sid)
  • mto: This revision was merged to the branch mainline in revision 42.
  • Revision ID: package-import@ubuntu.com-20120723085418-9foz30v6afaf5ffs
Tags: 2.63a-2
http://bugs.debian.org/658075
* debian/: Cycles support added (Closes: #658075)
  For now, this top feature has been enabled only
  on [any-amd64 any-i386] architectures because
  of OpenImageIO failing on all others
* debian/: scripts installation path changed
  from /usr/lib to /usr/share:
  + debian/patches/: patchset re-worked for path changing
  + debian/control: "Breaks" field added on yafaray-exporter

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/Eigen3/Eigen/src/Eigenvalues/ComplexEigenSolver.h
 
1
// This file is part of Eigen, a lightweight C++ template library
 
2
// for linear algebra.
 
3
//
 
4
// Copyright (C) 2009 Claire Maurice
 
5
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
 
6
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
 
7
//
 
8
// Eigen is free software; you can redistribute it and/or
 
9
// modify it under the terms of the GNU Lesser General Public
 
10
// License as published by the Free Software Foundation; either
 
11
// version 3 of the License, or (at your option) any later version.
 
12
//
 
13
// Alternatively, you can redistribute it and/or
 
14
// modify it under the terms of the GNU General Public License as
 
15
// published by the Free Software Foundation; either version 2 of
 
16
// the License, or (at your option) any later version.
 
17
//
 
18
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
 
19
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 
20
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
 
21
// GNU General Public License for more details.
 
22
//
 
23
// You should have received a copy of the GNU Lesser General Public
 
24
// License and a copy of the GNU General Public License along with
 
25
// Eigen. If not, see <http://www.gnu.org/licenses/>.
 
26
 
 
27
#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
 
28
#define EIGEN_COMPLEX_EIGEN_SOLVER_H
 
29
 
 
30
#include "./EigenvaluesCommon.h"
 
31
#include "./ComplexSchur.h"
 
32
 
 
33
/** \eigenvalues_module \ingroup Eigenvalues_Module
 
34
  *
 
35
  *
 
36
  * \class ComplexEigenSolver
 
37
  *
 
38
  * \brief Computes eigenvalues and eigenvectors of general complex matrices
 
39
  *
 
40
  * \tparam _MatrixType the type of the matrix of which we are
 
41
  * computing the eigendecomposition; this is expected to be an
 
42
  * instantiation of the Matrix class template.
 
43
  *
 
44
  * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
 
45
  * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v
 
46
  * \f$.  If \f$ D \f$ is a diagonal matrix with the eigenvalues on
 
47
  * the diagonal, and \f$ V \f$ is a matrix with the eigenvectors as
 
48
  * its columns, then \f$ A V = V D \f$. The matrix \f$ V \f$ is
 
49
  * almost always invertible, in which case we have \f$ A = V D V^{-1}
 
50
  * \f$. This is called the eigendecomposition.
 
51
  *
 
52
  * The main function in this class is compute(), which computes the
 
53
  * eigenvalues and eigenvectors of a given function. The
 
54
  * documentation for that function contains an example showing the
 
55
  * main features of the class.
 
56
  *
 
57
  * \sa class EigenSolver, class SelfAdjointEigenSolver
 
58
  */
 
59
template<typename _MatrixType> class ComplexEigenSolver
 
60
{
 
61
  public:
 
62
 
 
63
    /** \brief Synonym for the template parameter \p _MatrixType. */
 
64
    typedef _MatrixType MatrixType;
 
65
 
 
66
    enum {
 
67
      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
 
68
      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
 
69
      Options = MatrixType::Options,
 
70
      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
 
71
      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
 
72
    };
 
73
 
 
74
    /** \brief Scalar type for matrices of type #MatrixType. */
 
75
    typedef typename MatrixType::Scalar Scalar;
 
76
    typedef typename NumTraits<Scalar>::Real RealScalar;
 
77
    typedef typename MatrixType::Index Index;
 
78
 
 
79
    /** \brief Complex scalar type for #MatrixType.
 
80
      *
 
81
      * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
 
82
      * \c float or \c double) and just \c Scalar if #Scalar is
 
83
      * complex.
 
84
      */
 
85
    typedef std::complex<RealScalar> ComplexScalar;
 
86
 
 
87
    /** \brief Type for vector of eigenvalues as returned by eigenvalues().
 
88
      *
 
89
      * This is a column vector with entries of type #ComplexScalar.
 
90
      * The length of the vector is the size of #MatrixType.
 
91
      */
 
92
    typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType;
 
93
 
 
94
    /** \brief Type for matrix of eigenvectors as returned by eigenvectors().
 
95
      *
 
96
      * This is a square matrix with entries of type #ComplexScalar.
 
97
      * The size is the same as the size of #MatrixType.
 
98
      */
 
99
    typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType;
 
100
 
 
101
    /** \brief Default constructor.
 
102
      *
 
103
      * The default constructor is useful in cases in which the user intends to
 
104
      * perform decompositions via compute().
 
105
      */
 
106
    ComplexEigenSolver()
 
107
            : m_eivec(),
 
108
              m_eivalues(),
 
109
              m_schur(),
 
110
              m_isInitialized(false),
 
111
              m_eigenvectorsOk(false),
 
112
              m_matX()
 
113
    {}
 
114
 
 
115
    /** \brief Default Constructor with memory preallocation
 
116
      *
 
117
      * Like the default constructor but with preallocation of the internal data
 
118
      * according to the specified problem \a size.
 
119
      * \sa ComplexEigenSolver()
 
120
      */
 
121
    ComplexEigenSolver(Index size)
 
122
            : m_eivec(size, size),
 
123
              m_eivalues(size),
 
124
              m_schur(size),
 
125
              m_isInitialized(false),
 
126
              m_eigenvectorsOk(false),
 
127
              m_matX(size, size)
 
128
    {}
 
129
 
 
130
    /** \brief Constructor; computes eigendecomposition of given matrix.
 
131
      *
 
132
      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
 
133
      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
 
134
      *    eigenvalues are computed; if false, only the eigenvalues are
 
135
      *    computed.
 
136
      *
 
137
      * This constructor calls compute() to compute the eigendecomposition.
 
138
      */
 
139
      ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
 
140
            : m_eivec(matrix.rows(),matrix.cols()),
 
141
              m_eivalues(matrix.cols()),
 
142
              m_schur(matrix.rows()),
 
143
              m_isInitialized(false),
 
144
              m_eigenvectorsOk(false),
 
145
              m_matX(matrix.rows(),matrix.cols())
 
146
    {
 
147
      compute(matrix, computeEigenvectors);
 
148
    }
 
149
 
 
150
    /** \brief Returns the eigenvectors of given matrix.
 
151
      *
 
152
      * \returns  A const reference to the matrix whose columns are the eigenvectors.
 
153
      *
 
154
      * \pre Either the constructor
 
155
      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
 
156
      * function compute(const MatrixType& matrix, bool) has been called before
 
157
      * to compute the eigendecomposition of a matrix, and
 
158
      * \p computeEigenvectors was set to true (the default).
 
159
      *
 
160
      * This function returns a matrix whose columns are the eigenvectors. Column
 
161
      * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k
 
162
      * \f$ as returned by eigenvalues().  The eigenvectors are normalized to
 
163
      * have (Euclidean) norm equal to one. The matrix returned by this
 
164
      * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D
 
165
      * V^{-1} \f$, if it exists.
 
166
      *
 
167
      * Example: \include ComplexEigenSolver_eigenvectors.cpp
 
168
      * Output: \verbinclude ComplexEigenSolver_eigenvectors.out
 
169
      */
 
170
    const EigenvectorType& eigenvectors() const
 
171
    {
 
172
      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
 
173
      eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
 
174
      return m_eivec;
 
175
    }
 
176
 
 
177
    /** \brief Returns the eigenvalues of given matrix.
 
178
      *
 
179
      * \returns A const reference to the column vector containing the eigenvalues.
 
180
      *
 
181
      * \pre Either the constructor
 
182
      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
 
183
      * function compute(const MatrixType& matrix, bool) has been called before
 
184
      * to compute the eigendecomposition of a matrix.
 
185
      *
 
186
      * This function returns a column vector containing the
 
187
      * eigenvalues. Eigenvalues are repeated according to their
 
188
      * algebraic multiplicity, so there are as many eigenvalues as
 
189
      * rows in the matrix. The eigenvalues are not sorted in any particular
 
190
      * order.
 
191
      *
 
192
      * Example: \include ComplexEigenSolver_eigenvalues.cpp
 
193
      * Output: \verbinclude ComplexEigenSolver_eigenvalues.out
 
194
      */
 
195
    const EigenvalueType& eigenvalues() const
 
196
    {
 
197
      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
 
198
      return m_eivalues;
 
199
    }
 
200
 
 
201
    /** \brief Computes eigendecomposition of given matrix.
 
202
      *
 
203
      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
 
204
      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
 
205
      *    eigenvalues are computed; if false, only the eigenvalues are
 
206
      *    computed.
 
207
      * \returns    Reference to \c *this
 
208
      *
 
209
      * This function computes the eigenvalues of the complex matrix \p matrix.
 
210
      * The eigenvalues() function can be used to retrieve them.  If
 
211
      * \p computeEigenvectors is true, then the eigenvectors are also computed
 
212
      * and can be retrieved by calling eigenvectors().
 
213
      *
 
214
      * The matrix is first reduced to Schur form using the
 
215
      * ComplexSchur class. The Schur decomposition is then used to
 
216
      * compute the eigenvalues and eigenvectors.
 
217
      *
 
218
      * The cost of the computation is dominated by the cost of the
 
219
      * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$
 
220
      * is the size of the matrix.
 
221
      *
 
222
      * Example: \include ComplexEigenSolver_compute.cpp
 
223
      * Output: \verbinclude ComplexEigenSolver_compute.out
 
224
      */
 
225
    ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
 
226
 
 
227
    /** \brief Reports whether previous computation was successful.
 
228
      *
 
229
      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
 
230
      */
 
231
    ComputationInfo info() const
 
232
    {
 
233
      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
 
234
      return m_schur.info();
 
235
    }
 
236
 
 
237
  protected:
 
238
    EigenvectorType m_eivec;
 
239
    EigenvalueType m_eivalues;
 
240
    ComplexSchur<MatrixType> m_schur;
 
241
    bool m_isInitialized;
 
242
    bool m_eigenvectorsOk;
 
243
    EigenvectorType m_matX;
 
244
 
 
245
  private:
 
246
    void doComputeEigenvectors(RealScalar matrixnorm);
 
247
    void sortEigenvalues(bool computeEigenvectors);
 
248
};
 
249
 
 
250
 
 
251
template<typename MatrixType>
 
252
ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
 
253
{
 
254
  // this code is inspired from Jampack
 
255
  assert(matrix.cols() == matrix.rows());
 
256
 
 
257
  // Do a complex Schur decomposition, A = U T U^*
 
258
  // The eigenvalues are on the diagonal of T.
 
259
  m_schur.compute(matrix, computeEigenvectors);
 
260
 
 
261
  if(m_schur.info() == Success)
 
262
  {
 
263
    m_eivalues = m_schur.matrixT().diagonal();
 
264
    if(computeEigenvectors)
 
265
      doComputeEigenvectors(matrix.norm());
 
266
    sortEigenvalues(computeEigenvectors);
 
267
  }
 
268
 
 
269
  m_isInitialized = true;
 
270
  m_eigenvectorsOk = computeEigenvectors;
 
271
  return *this;
 
272
}
 
273
 
 
274
 
 
275
template<typename MatrixType>
 
276
void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
 
277
{
 
278
  const Index n = m_eivalues.size();
 
279
 
 
280
  // Compute X such that T = X D X^(-1), where D is the diagonal of T.
 
281
  // The matrix X is unit triangular.
 
282
  m_matX = EigenvectorType::Zero(n, n);
 
283
  for(Index k=n-1 ; k>=0 ; k--)
 
284
  {
 
285
    m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
 
286
    // Compute X(i,k) using the (i,k) entry of the equation X T = D X
 
287
    for(Index i=k-1 ; i>=0 ; i--)
 
288
    {
 
289
      m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
 
290
      if(k-i-1>0)
 
291
        m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
 
292
      ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
 
293
      if(z==ComplexScalar(0))
 
294
      {
 
295
        // If the i-th and k-th eigenvalue are equal, then z equals 0.
 
296
        // Use a small value instead, to prevent division by zero.
 
297
        internal::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
 
298
      }
 
299
      m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
 
300
    }
 
301
  }
 
302
 
 
303
  // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
 
304
  m_eivec.noalias() = m_schur.matrixU() * m_matX;
 
305
  // .. and normalize the eigenvectors
 
306
  for(Index k=0 ; k<n ; k++)
 
307
  {
 
308
    m_eivec.col(k).normalize();
 
309
  }
 
310
}
 
311
 
 
312
 
 
313
template<typename MatrixType>
 
314
void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
 
315
{
 
316
  const Index n =  m_eivalues.size();
 
317
  for (Index i=0; i<n; i++)
 
318
  {
 
319
    Index k;
 
320
    m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
 
321
    if (k != 0)
 
322
    {
 
323
      k += i;
 
324
      std::swap(m_eivalues[k],m_eivalues[i]);
 
325
      if(computeEigenvectors)
 
326
        m_eivec.col(i).swap(m_eivec.col(k));
 
327
    }
 
328
  }
 
329
}
 
330
 
 
331
 
 
332
#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H