RealSchur.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
5// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11#ifndef EIGEN_REAL_SCHUR_H
12#define EIGEN_REAL_SCHUR_H
13
14#include "./HessenbergDecomposition.h"
15
16namespace Eigen {
17
54template<typename _MatrixType> class RealSchur
55{
56 public:
57 typedef _MatrixType MatrixType;
58 enum {
59 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
60 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
61 Options = MatrixType::Options,
62 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
63 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
64 };
65 typedef typename MatrixType::Scalar Scalar;
66 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
67 typedef typename MatrixType::Index Index;
68
69 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
70 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
71
83 RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
84 : m_matT(size, size),
85 m_matU(size, size),
86 m_workspaceVector(size),
87 m_hess(size),
88 m_isInitialized(false),
89 m_matUisUptodate(false)
90 { }
91
102 RealSchur(const MatrixType& matrix, bool computeU = true)
103 : m_matT(matrix.rows(),matrix.cols()),
104 m_matU(matrix.rows(),matrix.cols()),
105 m_workspaceVector(matrix.rows()),
106 m_hess(matrix.rows()),
107 m_isInitialized(false),
108 m_matUisUptodate(false)
109 {
110 compute(matrix, computeU);
111 }
112
124 const MatrixType& matrixU() const
125 {
126 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
127 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition.");
128 return m_matU;
129 }
130
141 const MatrixType& matrixT() const
142 {
143 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
144 return m_matT;
145 }
146
164 RealSchur& compute(const MatrixType& matrix, bool computeU = true);
165
170 ComputationInfo info() const
171 {
172 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
173 return m_info;
174 }
175
180 static const int m_maxIterations = 40;
181
182 private:
183
184 MatrixType m_matT;
185 MatrixType m_matU;
186 ColumnVectorType m_workspaceVector;
187 HessenbergDecomposition<MatrixType> m_hess;
188 ComputationInfo m_info;
189 bool m_isInitialized;
190 bool m_matUisUptodate;
191
192 typedef Matrix<Scalar,3,1> Vector3s;
193
194 Scalar computeNormOfT();
195 Index findSmallSubdiagEntry(Index iu, Scalar norm);
196 void splitOffTwoRows(Index iu, bool computeU, Scalar exshift);
197 void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
198 void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
199 void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace);
200};
201
202
203template<typename MatrixType>
204RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
205{
206 assert(matrix.cols() == matrix.rows());
207
208 // Step 1. Reduce to Hessenberg form
209 m_hess.compute(matrix);
210 m_matT = m_hess.matrixH();
211 if (computeU)
212 m_matU = m_hess.matrixQ();
213
214 // Step 2. Reduce to real Schur form
215 m_workspaceVector.resize(m_matT.cols());
216 Scalar* workspace = &m_workspaceVector.coeffRef(0);
217
218 // The matrix m_matT is divided in three parts.
219 // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
220 // Rows il,...,iu is the part we are working on (the active window).
221 // Rows iu+1,...,end are already brought in triangular form.
222 Index iu = m_matT.cols() - 1;
223 Index iter = 0; // iteration count for current eigenvalue
224 Index totalIter = 0; // iteration count for whole matrix
225 Scalar exshift(0); // sum of exceptional shifts
226 Scalar norm = computeNormOfT();
227
228 if(norm!=0)
229 {
230 while (iu >= 0)
231 {
232 Index il = findSmallSubdiagEntry(iu, norm);
233
234 // Check for convergence
235 if (il == iu) // One root found
236 {
237 m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
238 if (iu > 0)
239 m_matT.coeffRef(iu, iu-1) = Scalar(0);
240 iu--;
241 iter = 0;
242 }
243 else if (il == iu-1) // Two roots found
244 {
245 splitOffTwoRows(iu, computeU, exshift);
246 iu -= 2;
247 iter = 0;
248 }
249 else // No convergence yet
250 {
251 // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
252 Vector3s firstHouseholderVector(0,0,0), shiftInfo;
253 computeShift(iu, iter, exshift, shiftInfo);
254 iter = iter + 1;
255 totalIter = totalIter + 1;
256 if (totalIter > m_maxIterations * matrix.cols()) break;
257 Index im;
258 initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
259 performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
260 }
261 }
262 }
263 if(totalIter <= m_maxIterations * matrix.cols())
264 m_info = Success;
265 else
266 m_info = NoConvergence;
267
268 m_isInitialized = true;
269 m_matUisUptodate = computeU;
270 return *this;
271}
272
274template<typename MatrixType>
275inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
276{
277 const Index size = m_matT.cols();
278 // FIXME to be efficient the following would requires a triangular reduxion code
279 // Scalar norm = m_matT.upper().cwiseAbs().sum()
280 // + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
281 Scalar norm(0);
282 for (Index j = 0; j < size; ++j)
283 norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
284 return norm;
285}
286
288template<typename MatrixType>
289inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, Scalar norm)
290{
291 Index res = iu;
292 while (res > 0)
293 {
294 Scalar s = internal::abs(m_matT.coeff(res-1,res-1)) + internal::abs(m_matT.coeff(res,res));
295 if (s == 0.0)
296 s = norm;
297 if (internal::abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
298 break;
299 res--;
300 }
301 return res;
302}
303
305template<typename MatrixType>
306inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scalar exshift)
307{
308 const Index size = m_matT.cols();
309
310 // The eigenvalues of the 2x2 matrix [a b; c d] are
311 // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
312 Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
313 Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); // q = tr^2 / 4 - det = discr/4
314 m_matT.coeffRef(iu,iu) += exshift;
315 m_matT.coeffRef(iu-1,iu-1) += exshift;
316
317 if (q >= Scalar(0)) // Two real eigenvalues
318 {
319 Scalar z = internal::sqrt(internal::abs(q));
320 JacobiRotation<Scalar> rot;
321 if (p >= Scalar(0))
322 rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
323 else
324 rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
325
326 m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
327 m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
328 m_matT.coeffRef(iu, iu-1) = Scalar(0);
329 if (computeU)
330 m_matU.applyOnTheRight(iu-1, iu, rot);
331 }
332
333 if (iu > 1)
334 m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
335}
336
338template<typename MatrixType>
339inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
340{
341 shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
342 shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
343 shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
344
345 // Wilkinson's original ad hoc shift
346 if (iter == 10)
347 {
348 exshift += shiftInfo.coeff(0);
349 for (Index i = 0; i <= iu; ++i)
350 m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
351 Scalar s = internal::abs(m_matT.coeff(iu,iu-1)) + internal::abs(m_matT.coeff(iu-1,iu-2));
352 shiftInfo.coeffRef(0) = Scalar(0.75) * s;
353 shiftInfo.coeffRef(1) = Scalar(0.75) * s;
354 shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
355 }
356
357 // MATLAB's new ad hoc shift
358 if (iter == 30)
359 {
360 Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
361 s = s * s + shiftInfo.coeff(2);
362 if (s > Scalar(0))
363 {
364 s = internal::sqrt(s);
365 if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
366 s = -s;
367 s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
368 s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
369 exshift += s;
370 for (Index i = 0; i <= iu; ++i)
371 m_matT.coeffRef(i,i) -= s;
372 shiftInfo.setConstant(Scalar(0.964));
373 }
374 }
375}
376
378template<typename MatrixType>
379inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
380{
381 Vector3s& v = firstHouseholderVector; // alias to save typing
382
383 for (im = iu-2; im >= il; --im)
384 {
385 const Scalar Tmm = m_matT.coeff(im,im);
386 const Scalar r = shiftInfo.coeff(0) - Tmm;
387 const Scalar s = shiftInfo.coeff(1) - Tmm;
388 v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
389 v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
390 v.coeffRef(2) = m_matT.coeff(im+2,im+1);
391 if (im == il) {
392 break;
393 }
394 const Scalar lhs = m_matT.coeff(im,im-1) * (internal::abs(v.coeff(1)) + internal::abs(v.coeff(2)));
395 const Scalar rhs = v.coeff(0) * (internal::abs(m_matT.coeff(im-1,im-1)) + internal::abs(Tmm) + internal::abs(m_matT.coeff(im+1,im+1)));
396 if (internal::abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
397 {
398 break;
399 }
400 }
401}
402
404template<typename MatrixType>
405inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
406{
407 assert(im >= il);
408 assert(im <= iu-2);
409
410 const Index size = m_matT.cols();
411
412 for (Index k = im; k <= iu-2; ++k)
413 {
414 bool firstIteration = (k == im);
415
416 Vector3s v;
417 if (firstIteration)
418 v = firstHouseholderVector;
419 else
420 v = m_matT.template block<3,1>(k,k-1);
421
422 Scalar tau, beta;
423 Matrix<Scalar, 2, 1> ess;
424 v.makeHouseholder(ess, tau, beta);
425
426 if (beta != Scalar(0)) // if v is not zero
427 {
428 if (firstIteration && k > il)
429 m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
430 else if (!firstIteration)
431 m_matT.coeffRef(k,k-1) = beta;
432
433 // These Householder transformations form the O(n^3) part of the algorithm
434 m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
435 m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
436 if (computeU)
437 m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
438 }
439 }
440
441 Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
442 Scalar tau, beta;
443 Matrix<Scalar, 1, 1> ess;
444 v.makeHouseholder(ess, tau, beta);
445
446 if (beta != Scalar(0)) // if v is not zero
447 {
448 m_matT.coeffRef(iu-1, iu-2) = beta;
449 m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
450 m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
451 if (computeU)
452 m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
453 }
454
455 // clean up pollution due to round-off errors
456 for (Index i = im+2; i <= iu; ++i)
457 {
458 m_matT.coeffRef(i,i-2) = Scalar(0);
459 if (i > im+2)
460 m_matT.coeffRef(i,i-3) = Scalar(0);
461 }
462}
463
464} // end namespace Eigen
465
466#endif // EIGEN_REAL_SCHUR_H
Eigen::HessenbergDecomposition
Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation.
Definition HessenbergDecomposition.h:58
Eigen::JacobiRotation
Rotation given by a cosine-sine pair.
Definition Jacobi.h:35
Eigen::JacobiRotation::makeGivens
void makeGivens(const Scalar &p, const Scalar &q, Scalar *z=0)
Definition Jacobi.h:145
Eigen::MatrixBase::makeHouseholder
void makeHouseholder(EssentialPart &essential, Scalar &tau, RealScalar &beta) const
Definition Householder.h:65
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition Matrix.h:129
Eigen::RealSchur::RealSchur
RealSchur(const MatrixType &matrix, bool computeU=true)
Constructor; computes real Schur decomposition of given matrix.
Definition RealSchur.h:102
Eigen::RealSchur::matrixU
const MatrixType & matrixU() const
Returns the orthogonal matrix in the Schur decomposition.
Definition RealSchur.h:124
Eigen::RealSchur::RealSchur
RealSchur(Index size=RowsAtCompileTime==Dynamic ? 1 :RowsAtCompileTime)
Default constructor.
Definition RealSchur.h:83
Eigen::RealSchur::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition RealSchur.h:170
Eigen::RealSchur::matrixT
const MatrixType & matrixT() const
Returns the quasi-triangular matrix in the Schur decomposition.
Definition RealSchur.h:141
Eigen::RealSchur::compute
RealSchur & compute(const MatrixType &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
Definition RealSchur.h:204
Eigen::RealSchur::m_maxIterations
static const int m_maxIterations
Maximum number of iterations.
Definition RealSchur.h:180
Eigen::ComputationInfo
ComputationInfo
Definition Constants.h:367
Eigen::NoConvergence
@ NoConvergence
Definition Constants.h:373
Eigen::Success
@ Success
Definition Constants.h:369