1
// This file is part of Eigen, a lightweight C++ template library
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>
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.
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.
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.
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/>.
27
#ifndef EIGEN_COMPLEX_SCHUR_H
28
#define EIGEN_COMPLEX_SCHUR_H
30
#include "./EigenvaluesCommon.h"
31
#include "./HessenbergDecomposition.h"
34
template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
37
/** \eigenvalues_module \ingroup Eigenvalues_Module
42
* \brief Performs a complex Schur decomposition of a real or complex square matrix
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.
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
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.
61
* \note This code is inspired from Jampack
63
* \sa class RealSchur, class EigenSolver, class ComplexEigenSolver
65
template<typename _MatrixType> class ComplexSchur
68
typedef _MatrixType MatrixType;
70
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
71
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
72
Options = MatrixType::Options,
73
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
74
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
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;
82
/** \brief Complex scalar type for \p _MatrixType.
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
88
typedef std::complex<RealScalar> ComplexScalar;
90
/** \brief Type for the matrices in the Schur decomposition.
92
* This is a square matrix with entries of type #ComplexScalar.
93
* The size is the same as the size of \p _MatrixType.
95
typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType;
97
/** \brief Default constructor.
99
* \param [in] size Positive integer, size of the matrix whose Schur decomposition will be computed.
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.
106
* \sa compute() for an example.
108
ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
112
m_isInitialized(false),
113
m_matUisUptodate(false)
116
/** \brief Constructor; computes Schur decomposition of given matrix.
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.
121
* This constructor calls compute() to compute the Schur decomposition.
123
* \sa matrixT() and matrixU() for examples.
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)
132
compute(matrix, computeU);
135
/** \brief Returns the unitary matrix in the Schur decomposition.
137
* \returns A const reference to the matrix U.
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
146
* Example: \include ComplexSchur_matrixU.cpp
147
* Output: \verbinclude ComplexSchur_matrixU.out
149
const ComplexMatrixType& matrixU() const
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.");
156
/** \brief Returns the triangular matrix in the Schur decomposition.
158
* \returns A const reference to the matrix T.
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
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
170
* Example: \include ComplexSchur_matrixT.cpp
171
* Output: \verbinclude ComplexSchur_matrixT.out
173
const ComplexMatrixType& matrixT() const
175
eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
179
/** \brief Computes Schur decomposition of given matrix.
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
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.
195
* Example: \include ComplexSchur_compute.cpp
196
* Output: \verbinclude ComplexSchur_compute.out
198
ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
200
/** \brief Reports whether previous computation was successful.
202
* \returns \c Success if computation was succesful, \c NoConvergence otherwise.
204
ComputationInfo info() const
206
eigen_assert(m_isInitialized && "RealSchur is not initialized.");
210
/** \brief Maximum number of iterations.
212
* Maximum number of iterations allowed for an eigenvalue to converge.
214
static const int m_maxIterations = 30;
217
ComplexMatrixType m_matT, m_matU;
218
HessenbergDecomposition<MatrixType> m_hess;
219
ComputationInfo m_info;
220
bool m_isInitialized;
221
bool m_matUisUptodate;
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>;
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)
236
RealScalar t, tre, tim;
240
if (abs(real(z)) <= abs(imag(z)))
242
// No cancellation in these formulas
243
tre = sqrt(RealScalar(0.5)*(t + real(z)));
244
tim = sqrt(RealScalar(0.5)*(t - real(z)));
248
// Stable computation of the above formulas
249
if (z.real() > RealScalar(0))
252
tim = abs(imag(z))*sqrt(RealScalar(0.5)/tre);
253
tre = sqrt(RealScalar(0.5)*tre);
258
tre = abs(imag(z))*sqrt(RealScalar(0.5)/tim);
259
tim = sqrt(RealScalar(0.5)*tim);
262
if(z.imag() < RealScalar(0))
265
return (std::complex<RealScalar>(tre,tim));
267
} // end namespace internal
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)
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()))
280
m_matT.coeffRef(i+1,i) = ComplexScalar(0);
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)
291
if (iter == 10 || iter == 20)
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)));
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
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);
311
if(internal::norm1(eival1) > internal::norm1(eival2))
312
eival2 = det / eival1;
314
eival1 = det / eival2;
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;
320
return normt * eival2;
324
template<typename MatrixType>
325
ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
327
m_matUisUptodate = false;
328
eigen_assert(matrix.cols() == matrix.rows());
330
if(matrix.cols() == 1)
332
m_matT = matrix.template cast<ComplexScalar>();
333
if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
335
m_isInitialized = true;
336
m_matUisUptodate = computeU;
340
internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
341
reduceToTriangularForm(computeU);
347
/* Reduce given matrix to Hessenberg form */
348
template<typename MatrixType, bool IsComplex>
349
struct complex_schur_reduce_to_hessenberg
351
// this is the implementation for the case IsComplex = true
352
static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
354
_this.m_hess.compute(matrix);
355
_this.m_matT = _this.m_hess.matrixH();
356
if(computeU) _this.m_matU = _this.m_hess.matrixQ();
360
template<typename MatrixType>
361
struct complex_schur_reduce_to_hessenberg<MatrixType, false>
363
static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
365
typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
366
typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType;
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>();
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>();
380
} // end namespace internal
382
// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
383
template<typename MatrixType>
384
void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
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;
392
Index iter = 0; // number of iterations we are working on the (iu,iu) element
396
// find iu, the bottom row of the active submatrix
399
if(!subdiagonalEntryIsNeglegible(iu-1)) break;
404
// if iu is zero then we are done; the whole matrix is triangularized
407
// if we spent too many iterations on the current element, we give up
409
if(iter > m_maxIterations) break;
411
// find il, the top row of the active submatrix
413
while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
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. */
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);
429
for(Index i=il+1 ; i<iu ; i++)
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);
439
if(iter <= m_maxIterations)
442
m_info = NoConvergence;
444
m_isInitialized = true;
445
m_matUisUptodate = computeU;
448
#endif // EIGEN_COMPLEX_SCHUR_H