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 
16 namespace Eigen {
17 
54 template<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 
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 
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;
188  ComputationInfo m_info;
189  bool m_isInitialized;
190  bool m_matUisUptodate;
191 
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 
203 template<typename MatrixType>
204 RealSchur<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 
274 template<typename MatrixType>
275 inline 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 
288 template<typename MatrixType>
289 inline 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 
305 template<typename MatrixType>
306 inline 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 
338 template<typename MatrixType>
339 inline 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 
378 template<typename MatrixType>
379 inline 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 
404 template<typename MatrixType>
405 inline 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