source: pacpussensors/trunk/Vislab/lib3dv/eigen/Eigen/src/Eigenvalues/ComplexSchur.h@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

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