Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
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,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
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
64 typedef typename MatrixType::Scalar Scalar;
65 typedef typename NumTraits<Scalar>::Real RealScalar;
67
74 typedef std::complex<RealScalar> ComplexScalar;
75
82
94 explicit 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
112 template<typename InputType>
113 explicit ComplexSchur(const EigenBase<InputType>& matrix, bool computeU = true)
114 : m_matT(matrix.rows(),matrix.cols()),
115 m_matU(matrix.rows(),matrix.cols()),
116 m_hess(matrix.rows()),
117 m_isInitialized(false),
118 m_matUisUptodate(false),
119 m_maxIters(-1)
120 {
121 compute(matrix.derived(), computeU);
122 }
123
139 {
140 eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
141 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
142 return m_matU;
143 }
144
163 {
164 eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
165 return m_matT;
166 }
167
190 template<typename InputType>
191 ComplexSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true);
192
210 template<typename HessMatrixType, typename OrthMatrixType>
211 ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU=true);
212
218 {
219 eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
220 return m_info;
221 }
222
229 {
230 m_maxIters = maxIters;
231 return *this;
232 }
233
236 {
237 return m_maxIters;
238 }
239
245 static const int m_maxIterationsPerRow = 30;
246
247 protected:
248 ComplexMatrixType m_matT, m_matU;
250 ComputationInfo m_info;
251 bool m_isInitialized;
252 bool m_matUisUptodate;
253 Index m_maxIters;
254
255 private:
256 bool subdiagonalEntryIsNeglegible(Index i);
257 ComplexScalar computeShift(Index iu, Index iter);
258 void reduceToTriangularForm(bool computeU);
259 friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
260};
261
265template<typename MatrixType>
266inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
267{
268 RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1));
269 RealScalar sd = numext::norm1(m_matT.coeff(i+1,i));
270 if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
271 {
272 m_matT.coeffRef(i+1,i) = ComplexScalar(0);
273 return true;
274 }
275 return false;
276}
277
278
280template<typename MatrixType>
281typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
282{
283 using std::abs;
284 if (iter == 10 || iter == 20)
285 {
286 // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
287 return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2)));
288 }
289
290 // compute the shift as one of the eigenvalues of t, the 2x2
291 // diagonal block on the bottom of the active submatrix
292 Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
293 RealScalar normt = t.cwiseAbs().sum();
294 t /= normt; // the normalization by sf is to avoid under/overflow
295
296 ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
297 ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
298 ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
299 ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
300 ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
301 ComplexScalar eival1 = (trace + disc) / RealScalar(2);
302 ComplexScalar eival2 = (trace - disc) / RealScalar(2);
303 RealScalar eival1_norm = numext::norm1(eival1);
304 RealScalar eival2_norm = numext::norm1(eival2);
305 // A division by zero can only occur if eival1==eival2==0.
306 // In this case, det==0, and all we have to do is checking that eival2_norm!=0
307 if(eival1_norm > eival2_norm)
308 eival2 = det / eival1;
309 else if(eival2_norm!=RealScalar(0))
310 eival1 = det / eival2;
311
312 // choose the eigenvalue closest to the bottom entry of the diagonal
313 if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1)))
314 return normt * eival1;
315 else
316 return normt * eival2;
317}
318
319
320template<typename MatrixType>
321template<typename InputType>
322ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU)
323{
324 m_matUisUptodate = false;
325 eigen_assert(matrix.cols() == matrix.rows());
326
327 if(matrix.cols() == 1)
328 {
329 m_matT = matrix.derived().template cast<ComplexScalar>();
330 if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
331 m_info = Success;
332 m_isInitialized = true;
333 m_matUisUptodate = computeU;
334 return *this;
335 }
336
337 internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix.derived(), computeU);
338 computeFromHessenberg(m_matT, m_matU, computeU);
339 return *this;
340}
341
342template<typename MatrixType>
343template<typename HessMatrixType, typename OrthMatrixType>
344ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
345{
346 m_matT = matrixH;
347 if(computeU)
348 m_matU = matrixQ;
349 reduceToTriangularForm(computeU);
350 return *this;
351}
352namespace internal {
353
354/* Reduce given matrix to Hessenberg form */
355template<typename MatrixType, bool IsComplex>
356struct complex_schur_reduce_to_hessenberg
357{
358 // this is the implementation for the case IsComplex = true
359 static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
360 {
361 _this.m_hess.compute(matrix);
362 _this.m_matT = _this.m_hess.matrixH();
363 if(computeU) _this.m_matU = _this.m_hess.matrixQ();
364 }
365};
366
367template<typename MatrixType>
368struct complex_schur_reduce_to_hessenberg<MatrixType, false>
369{
370 static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
371 {
372 typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
373
374 // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
375 _this.m_hess.compute(matrix);
376 _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
377 if(computeU)
378 {
379 // This may cause an allocation which seems to be avoidable
380 MatrixType Q = _this.m_hess.matrixQ();
381 _this.m_matU = Q.template cast<ComplexScalar>();
382 }
383 }
384};
385
386} // end namespace internal
387
388// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
389template<typename MatrixType>
390void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
391{
392 Index maxIters = m_maxIters;
393 if (maxIters == -1)
394 maxIters = m_maxIterationsPerRow * m_matT.rows();
395
396 // The matrix m_matT is divided in three parts.
397 // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
398 // Rows il,...,iu is the part we are working on (the active submatrix).
399 // Rows iu+1,...,end are already brought in triangular form.
400 Index iu = m_matT.cols() - 1;
401 Index il;
402 Index iter = 0; // number of iterations we are working on the (iu,iu) element
403 Index totalIter = 0; // number of iterations for whole matrix
404
405 while(true)
406 {
407 // find iu, the bottom row of the active submatrix
408 while(iu > 0)
409 {
410 if(!subdiagonalEntryIsNeglegible(iu-1)) break;
411 iter = 0;
412 --iu;
413 }
414
415 // if iu is zero then we are done; the whole matrix is triangularized
416 if(iu==0) break;
417
418 // if we spent too many iterations, we give up
419 iter++;
420 totalIter++;
421 if(totalIter > maxIters) break;
422
423 // find il, the top row of the active submatrix
424 il = iu-1;
425 while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
426 {
427 --il;
428 }
429
430 /* perform the QR step using Givens rotations. The first rotation
431 creates a bulge; the (il+2,il) element becomes nonzero. This
432 bulge is chased down to the bottom of the active submatrix. */
433
434 ComplexScalar shift = computeShift(iu, iter);
435 JacobiRotation<ComplexScalar> rot;
436 rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
437 m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
438 m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
439 if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
440
441 for(Index i=il+1 ; i<iu ; i++)
442 {
443 rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
444 m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
445 m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
446 m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
447 if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
448 }
449 }
450
451 if(totalIter <= maxIters)
452 m_info = Success;
453 else
454 m_info = NoConvergence;
455
456 m_isInitialized = true;
457 m_matUisUptodate = computeU;
458}
459
460} // end namespace Eigen
461
462#endif // EIGEN_COMPLEX_SCHUR_H
Performs a complex Schur decomposition of a real or complex square matrix.
Definition: ComplexSchur.h:52
const ComplexMatrixType & matrixT() const
Returns the triangular matrix in the Schur decomposition.
Definition: ComplexSchur.h:162
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: ComplexSchur.h:235
ComplexSchur & computeFromHessenberg(const HessMatrixType &matrixH, const OrthMatrixType &matrixQ, bool computeU=true)
Compute Schur decomposition from a given Hessenberg matrix.
Eigen::Index Index
Definition: ComplexSchur.h:66
const ComplexMatrixType & matrixU() const
Returns the unitary matrix in the Schur decomposition.
Definition: ComplexSchur.h:138
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: ComplexSchur.h:217
MatrixType::Scalar Scalar
Scalar type for matrices of type _MatrixType.
Definition: ComplexSchur.h:64
ComplexSchur(const EigenBase< InputType > &matrix, bool computeU=true)
Constructor; computes Schur decomposition of given matrix.
Definition: ComplexSchur.h:113
ComplexSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: ComplexSchur.h:228
static const int m_maxIterationsPerRow
Maximum number of iterations per row.
Definition: ComplexSchur.h:245
ComplexSchur(Index size=RowsAtCompileTime==Dynamic ? 1 :RowsAtCompileTime)
Default constructor.
Definition: ComplexSchur.h:94
std::complex< RealScalar > ComplexScalar
Complex scalar type for _MatrixType.
Definition: ComplexSchur.h:74
ComplexSchur & compute(const EigenBase< InputType > &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > ComplexMatrixType
Type for the matrices in the Schur decomposition.
Definition: ComplexSchur.h:81
Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation.
Definition: HessenbergDecomposition.h:58
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:175
const Scalar & coeff(Index rowId, Index colId) const
Definition: PlainObjectBase.h:152
ComputationInfo
Definition: Constants.h:440
@ Success
Definition: Constants.h:442
@ NoConvergence
Definition: Constants.h:446
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const int Dynamic
Definition: Constants.h:22
Definition: EigenBase.h:30
Derived & derived()
Definition: EigenBase.h:46
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:233