← Back to branch summary

~siretart/ubuntu/utopic/blender/libav10

« back to all changes in this revision

Viewing changes to extern/Eigen3/Eigen/src/Eigenvalues/ComplexSchur.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/ComplexSchur.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_SCHUR_H
 
28
#define EIGEN_COMPLEX_SCHUR_H
 
29
 
 
30
#include "./EigenvaluesCommon.h"
 
31
#include "./HessenbergDecomposition.h"
 
32
 
 
33
namespace internal {
 
34
template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
 
35
}
 
36
 
 
37
/** \eigenvalues_module \ingroup Eigenvalues_Module
 
38
  *
 
39
  *
 
40
  * \class ComplexSchur
 
41
  *
 
42
  * \brief Performs a complex Schur decomposition of a real or complex square matrix
 
43
  *
 
44
  * \tparam _MatrixType the type of the matrix of which we are
 
45
  * computing the Schur decomposition; this is expected to be an
 
46
  * instantiation of the Matrix class template.
 
47
  *
 
48
  * Given a real or complex square matrix A, this class computes the
 
49
  * Schur decomposition: \f$ A = U T U^*\f$ where U is a unitary
 
50
  * complex matrix, and T is a complex upper triangular matrix.  The
 
51
  * diagonal of the matrix T corresponds to the eigenvalues of the
 
52
  * matrix A.
 
53
  *
 
54
  * Call the function compute() to compute the Schur decomposition of
 
55
  * a given matrix. Alternatively, you can use the 
 
56
  * ComplexSchur(const MatrixType&, bool) constructor which computes
 
57
  * the Schur decomposition at construction time. Once the
 
58
  * decomposition is computed, you can use the matrixU() and matrixT()
 
59
  * functions to retrieve the matrices U and V in the decomposition.
 
60
  *
 
61
  * \note This code is inspired from Jampack
 
62
  *
 
63
  * \sa class RealSchur, class EigenSolver, class ComplexEigenSolver
 
64
  */
 
65
template<typename _MatrixType> class ComplexSchur
 
66
{
 
67
  public:
 
68
    typedef _MatrixType MatrixType;
 
69
    enum {
 
70
      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
 
71
      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
 
72
      Options = MatrixType::Options,
 
73
      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
 
74
      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
 
75
    };
 
76
 
 
77
    /** \brief Scalar type for matrices of type \p _MatrixType. */
 
78
    typedef typename MatrixType::Scalar Scalar;
 
79
    typedef typename NumTraits<Scalar>::Real RealScalar;
 
80
    typedef typename MatrixType::Index Index;
 
81
 
 
82
    /** \brief Complex scalar type for \p _MatrixType. 
 
83
      *
 
84
      * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
 
85
      * \c float or \c double) and just \c Scalar if #Scalar is
 
86
      * complex.
 
87
      */
 
88
    typedef std::complex<RealScalar> ComplexScalar;
 
89
 
 
90
    /** \brief Type for the matrices in the Schur decomposition.
 
91
      *
 
92
      * This is a square matrix with entries of type #ComplexScalar. 
 
93
      * The size is the same as the size of \p _MatrixType.
 
94
      */
 
95
    typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType;
 
96
 
 
97
    /** \brief Default constructor.
 
98
      *
 
99
      * \param [in] size  Positive integer, size of the matrix whose Schur decomposition will be computed.
 
100
      *
 
101
      * The default constructor is useful in cases in which the user
 
102
      * intends to perform decompositions via compute().  The \p size
 
103
      * parameter is only used as a hint. It is not an error to give a
 
104
      * wrong \p size, but it may impair performance.
 
105
      *
 
106
      * \sa compute() for an example.
 
107
      */
 
108
    ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
 
109
      : m_matT(size,size),
 
110
        m_matU(size,size),
 
111
        m_hess(size),
 
112
        m_isInitialized(false),
 
113
        m_matUisUptodate(false)
 
114
    {}
 
115
 
 
116
    /** \brief Constructor; computes Schur decomposition of given matrix. 
 
117
      * 
 
118
      * \param[in]  matrix    Square matrix whose Schur decomposition is to be computed.
 
119
      * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
 
120
      *
 
121
      * This constructor calls compute() to compute the Schur decomposition.
 
122
      *
 
123
      * \sa matrixT() and matrixU() for examples.
 
124
      */
 
125
    ComplexSchur(const MatrixType& matrix, bool computeU = true)
 
126
            : m_matT(matrix.rows(),matrix.cols()),
 
127
              m_matU(matrix.rows(),matrix.cols()),
 
128
              m_hess(matrix.rows()),
 
129
              m_isInitialized(false),
 
130
              m_matUisUptodate(false)
 
131
    {
 
132
      compute(matrix, computeU);
 
133
    }
 
134
 
 
135
    /** \brief Returns the unitary matrix in the Schur decomposition. 
 
136
      *
 
137
      * \returns A const reference to the matrix U.
 
138
      *
 
139
      * It is assumed that either the constructor
 
140
      * ComplexSchur(const MatrixType& matrix, bool computeU) or the
 
141
      * member function compute(const MatrixType& matrix, bool computeU)
 
142
      * has been called before to compute the Schur decomposition of a
 
143
      * matrix, and that \p computeU was set to true (the default
 
144
      * value).
 
145
      *
 
146
      * Example: \include ComplexSchur_matrixU.cpp
 
147
      * Output: \verbinclude ComplexSchur_matrixU.out
 
148
      */
 
149
    const ComplexMatrixType& matrixU() const
 
150
    {
 
151
      eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
 
152
      eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
 
153
      return m_matU;
 
154
    }
 
155
 
 
156
    /** \brief Returns the triangular matrix in the Schur decomposition. 
 
157
      *
 
158
      * \returns A const reference to the matrix T.
 
159
      *
 
160
      * It is assumed that either the constructor
 
161
      * ComplexSchur(const MatrixType& matrix, bool computeU) or the
 
162
      * member function compute(const MatrixType& matrix, bool computeU)
 
163
      * has been called before to compute the Schur decomposition of a
 
164
      * matrix.
 
165
      *
 
166
      * Note that this function returns a plain square matrix. If you want to reference
 
167
      * only the upper triangular part, use:
 
168
      * \code schur.matrixT().triangularView<Upper>() \endcode 
 
169
      *
 
170
      * Example: \include ComplexSchur_matrixT.cpp
 
171
      * Output: \verbinclude ComplexSchur_matrixT.out
 
172
      */
 
173
    const ComplexMatrixType& matrixT() const
 
174
    {
 
175
      eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
 
176
      return m_matT;
 
177
    }
 
178
 
 
179
    /** \brief Computes Schur decomposition of given matrix. 
 
180
      * 
 
181
      * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
 
182
      * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
 
183
      * \returns    Reference to \c *this
 
184
      *
 
185
      * The Schur decomposition is computed by first reducing the
 
186
      * matrix to Hessenberg form using the class
 
187
      * HessenbergDecomposition. The Hessenberg matrix is then reduced
 
188
      * to triangular form by performing QR iterations with a single
 
189
      * shift. The cost of computing the Schur decomposition depends
 
190
      * on the number of iterations; as a rough guide, it may be taken
 
191
      * on the number of iterations; as a rough guide, it may be taken
 
192
      * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops
 
193
      * if \a computeU is false.
 
194
      *
 
195
      * Example: \include ComplexSchur_compute.cpp
 
196
      * Output: \verbinclude ComplexSchur_compute.out
 
197
      */
 
198
    ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
 
199
 
 
200
    /** \brief Reports whether previous computation was successful.
 
201
      *
 
202
      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
 
203
      */
 
204
    ComputationInfo info() const
 
205
    {
 
206
      eigen_assert(m_isInitialized && "RealSchur is not initialized.");
 
207
      return m_info;
 
208
    }
 
209
 
 
210
    /** \brief Maximum number of iterations.
 
211
      *
 
212
      * Maximum number of iterations allowed for an eigenvalue to converge. 
 
213
      */
 
214
    static const int m_maxIterations = 30;
 
215
 
 
216
  protected:
 
217
    ComplexMatrixType m_matT, m_matU;
 
218
    HessenbergDecomposition<MatrixType> m_hess;
 
219
    ComputationInfo m_info;
 
220
    bool m_isInitialized;
 
221
    bool m_matUisUptodate;
 
222
 
 
223
  private:  
 
224
    bool subdiagonalEntryIsNeglegible(Index i);
 
225
    ComplexScalar computeShift(Index iu, Index iter);
 
226
    void reduceToTriangularForm(bool computeU);
 
227
    friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
 
228
};
 
229
 
 
230
namespace internal {
 
231
 
 
232
/** Computes the principal value of the square root of the complex \a z. */
 
233
template<typename RealScalar>
 
234
std::complex<RealScalar> sqrt(const std::complex<RealScalar> &z)
 
235
{
 
236
  RealScalar t, tre, tim;
 
237
 
 
238
  t = abs(z);
 
239
 
 
240
  if (abs(real(z)) <= abs(imag(z)))
 
241
  {
 
242
    // No cancellation in these formulas
 
243
    tre = sqrt(RealScalar(0.5)*(t + real(z)));
 
244
    tim = sqrt(RealScalar(0.5)*(t - real(z)));
 
245
  }
 
246
  else
 
247
  {
 
248
    // Stable computation of the above formulas
 
249
    if (z.real() > RealScalar(0))
 
250
    {
 
251
      tre = t + z.real();
 
252
      tim = abs(imag(z))*sqrt(RealScalar(0.5)/tre);
 
253
      tre = sqrt(RealScalar(0.5)*tre);
 
254
    }
 
255
    else
 
256
    {
 
257
      tim = t - z.real();
 
258
      tre = abs(imag(z))*sqrt(RealScalar(0.5)/tim);
 
259
      tim = sqrt(RealScalar(0.5)*tim);
 
260
    }
 
261
  }
 
262
  if(z.imag() < RealScalar(0))
 
263
    tim = -tim;
 
264
 
 
265
  return (std::complex<RealScalar>(tre,tim));
 
266
}
 
267
} // end namespace internal
 
268
 
 
269
 
 
270
/** If m_matT(i+1,i) is neglegible in floating point arithmetic
 
271
  * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and
 
272
  * return true, else return false. */
 
273
template<typename MatrixType>
 
274
inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
 
275
{
 
276
  RealScalar d = internal::norm1(m_matT.coeff(i,i)) + internal::norm1(m_matT.coeff(i+1,i+1));
 
277
  RealScalar sd = internal::norm1(m_matT.coeff(i+1,i));
 
278
  if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
 
279
  {
 
280
    m_matT.coeffRef(i+1,i) = ComplexScalar(0);
 
281
    return true;
 
282
  }
 
283
  return false;
 
284
}
 
285
 
 
286
 
 
287
/** Compute the shift in the current QR iteration. */
 
288
template<typename MatrixType>
 
289
typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
 
290
{
 
291
  if (iter == 10 || iter == 20) 
 
292
  {
 
293
    // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
 
294
    return internal::abs(internal::real(m_matT.coeff(iu,iu-1))) + internal::abs(internal::real(m_matT.coeff(iu-1,iu-2)));
 
295
  }
 
296
 
 
297
  // compute the shift as one of the eigenvalues of t, the 2x2
 
298
  // diagonal block on the bottom of the active submatrix
 
299
  Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
 
300
  RealScalar normt = t.cwiseAbs().sum();
 
301
  t /= normt;     // the normalization by sf is to avoid under/overflow
 
302
 
 
303
  ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
 
304
  ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
 
305
  ComplexScalar disc = internal::sqrt(c*c + RealScalar(4)*b);
 
306
  ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
 
307
  ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
 
308
  ComplexScalar eival1 = (trace + disc) / RealScalar(2);
 
309
  ComplexScalar eival2 = (trace - disc) / RealScalar(2);
 
310
 
 
311
  if(internal::norm1(eival1) > internal::norm1(eival2))
 
312
    eival2 = det / eival1;
 
313
  else
 
314
    eival1 = det / eival2;
 
315
 
 
316
  // choose the eigenvalue closest to the bottom entry of the diagonal
 
317
  if(internal::norm1(eival1-t.coeff(1,1)) < internal::norm1(eival2-t.coeff(1,1)))
 
318
    return normt * eival1;
 
319
  else
 
320
    return normt * eival2;
 
321
}
 
322
 
 
323
 
 
324
template<typename MatrixType>
 
325
ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
 
326
{
 
327
  m_matUisUptodate = false;
 
328
  eigen_assert(matrix.cols() == matrix.rows());
 
329
 
 
330
  if(matrix.cols() == 1)
 
331
  {
 
332
    m_matT = matrix.template cast<ComplexScalar>();
 
333
    if(computeU)  m_matU = ComplexMatrixType::Identity(1,1);
 
334
    m_info = Success;
 
335
    m_isInitialized = true;
 
336
    m_matUisUptodate = computeU;
 
337
    return *this;
 
338
  }
 
339
 
 
340
  internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
 
341
  reduceToTriangularForm(computeU);
 
342
  return *this;
 
343
}
 
344
 
 
345
namespace internal {
 
346
 
 
347
/* Reduce given matrix to Hessenberg form */
 
348
template<typename MatrixType, bool IsComplex>
 
349
struct complex_schur_reduce_to_hessenberg
 
350
{
 
351
  // this is the implementation for the case IsComplex = true
 
352
  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
 
353
  {
 
354
    _this.m_hess.compute(matrix);
 
355
    _this.m_matT = _this.m_hess.matrixH();
 
356
    if(computeU)  _this.m_matU = _this.m_hess.matrixQ();
 
357
  }
 
358
};
 
359
 
 
360
template<typename MatrixType>
 
361
struct complex_schur_reduce_to_hessenberg<MatrixType, false>
 
362
{
 
363
  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
 
364
  {
 
365
    typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
 
366
    typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType;
 
367
 
 
368
    // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
 
369
    _this.m_hess.compute(matrix);
 
370
    _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
 
371
    if(computeU)  
 
372
    {
 
373
      // This may cause an allocation which seems to be avoidable
 
374
      MatrixType Q = _this.m_hess.matrixQ(); 
 
375
      _this.m_matU = Q.template cast<ComplexScalar>();
 
376
    }
 
377
  }
 
378
};
 
379
 
 
380
} // end namespace internal
 
381
 
 
382
// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
 
383
template<typename MatrixType>
 
384
void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
 
385
{  
 
386
  // The matrix m_matT is divided in three parts. 
 
387
  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 
 
388
  // Rows il,...,iu is the part we are working on (the active submatrix).
 
389
  // Rows iu+1,...,end are already brought in triangular form.
 
390
  Index iu = m_matT.cols() - 1;
 
391
  Index il;
 
392
  Index iter = 0; // number of iterations we are working on the (iu,iu) element
 
393
 
 
394
  while(true)
 
395
  {
 
396
    // find iu, the bottom row of the active submatrix
 
397
    while(iu > 0)
 
398
    {
 
399
      if(!subdiagonalEntryIsNeglegible(iu-1)) break;
 
400
      iter = 0;
 
401
      --iu;
 
402
    }
 
403
 
 
404
    // if iu is zero then we are done; the whole matrix is triangularized
 
405
    if(iu==0) break;
 
406
 
 
407
    // if we spent too many iterations on the current element, we give up
 
408
    iter++;
 
409
    if(iter > m_maxIterations) break;
 
410
 
 
411
    // find il, the top row of the active submatrix
 
412
    il = iu-1;
 
413
    while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
 
414
    {
 
415
      --il;
 
416
    }
 
417
 
 
418
    /* perform the QR step using Givens rotations. The first rotation
 
419
       creates a bulge; the (il+2,il) element becomes nonzero. This
 
420
       bulge is chased down to the bottom of the active submatrix. */
 
421
 
 
422
    ComplexScalar shift = computeShift(iu, iter);
 
423
    JacobiRotation<ComplexScalar> rot;
 
424
    rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
 
425
    m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
 
426
    m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
 
427
    if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
 
428
 
 
429
    for(Index i=il+1 ; i<iu ; i++)
 
430
    {
 
431
      rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
 
432
      m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
 
433
      m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
 
434
      m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
 
435
      if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
 
436
    }
 
437
  }
 
438
 
 
439
  if(iter <= m_maxIterations) 
 
440
    m_info = Success;
 
441
  else
 
442
    m_info = NoConvergence;
 
443
 
 
444
  m_isInitialized = true;
 
445
  m_matUisUptodate = computeU;
 
446
}
 
447
 
 
448
#endif // EIGEN_COMPLEX_SCHUR_H