RealSchur.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 #ifndef EIGEN_REAL_SCHUR_H
27 #define EIGEN_REAL_SCHUR_H
28 
30 
31 namespace Eigen {
32 
69 template<typename _MatrixType> class RealSchur
70 {
71  public:
72  typedef _MatrixType MatrixType;
73  enum {
78  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
79  };
80  typedef typename MatrixType::Scalar Scalar;
81  typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
82  typedef typename MatrixType::Index Index;
83 
86 
99  : m_matT(size, size),
100  m_matU(size, size),
101  m_workspaceVector(size),
102  m_hess(size),
103  m_isInitialized(false),
104  m_matUisUptodate(false)
105  { }
106 
117  RealSchur(const MatrixType& matrix, bool computeU = true)
118  : m_matT(matrix.rows(),matrix.cols()),
119  m_matU(matrix.rows(),matrix.cols()),
120  m_workspaceVector(matrix.rows()),
121  m_hess(matrix.rows()),
122  m_isInitialized(false),
123  m_matUisUptodate(false)
124  {
125  compute(matrix, computeU);
126  }
127 
139  const MatrixType& matrixU() const
140  {
141  eigen_assert(m_isInitialized && "RealSchur is not initialized.");
142  eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition.");
143  return m_matU;
144  }
145 
156  const MatrixType& matrixT() const
157  {
158  eigen_assert(m_isInitialized && "RealSchur is not initialized.");
159  return m_matT;
160  }
161 
179  RealSchur& compute(const MatrixType& matrix, bool computeU = true);
180 
186  {
187  eigen_assert(m_isInitialized && "RealSchur is not initialized.");
188  return m_info;
189  }
190 
195  static const int m_maxIterations = 40;
196 
197  private:
198 
199  MatrixType m_matT;
200  MatrixType m_matU;
201  ColumnVectorType m_workspaceVector;
203  ComputationInfo m_info;
204  bool m_isInitialized;
205  bool m_matUisUptodate;
206 
208 
209  Scalar computeNormOfT();
210  Index findSmallSubdiagEntry(Index iu, Scalar norm);
211  void splitOffTwoRows(Index iu, bool computeU, Scalar exshift);
212  void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
213  void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
214  void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace);
215 };
216 
217 
218 template<typename MatrixType>
220 {
221  assert(matrix.cols() == matrix.rows());
222 
223  // Step 1. Reduce to Hessenberg form
224  m_hess.compute(matrix);
225  m_matT = m_hess.matrixH();
226  if (computeU)
227  m_matU = m_hess.matrixQ();
228 
229  // Step 2. Reduce to real Schur form
230  m_workspaceVector.resize(m_matT.cols());
231  Scalar* workspace = &m_workspaceVector.coeffRef(0);
232 
233  // The matrix m_matT is divided in three parts.
234  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
235  // Rows il,...,iu is the part we are working on (the active window).
236  // Rows iu+1,...,end are already brought in triangular form.
237  Index iu = m_matT.cols() - 1;
238  Index iter = 0; // iteration count
239  Scalar exshift(0); // sum of exceptional shifts
240  Scalar norm = computeNormOfT();
241 
242  while (iu >= 0)
243  {
244  Index il = findSmallSubdiagEntry(iu, norm);
245 
246  // Check for convergence
247  if (il == iu) // One root found
248  {
249  m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
250  if (iu > 0)
251  m_matT.coeffRef(iu, iu-1) = Scalar(0);
252  iu--;
253  iter = 0;
254  }
255  else if (il == iu-1) // Two roots found
256  {
257  splitOffTwoRows(iu, computeU, exshift);
258  iu -= 2;
259  iter = 0;
260  }
261  else // No convergence yet
262  {
263  // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
264  Vector3s firstHouseholderVector(0,0,0), shiftInfo;
265  computeShift(iu, iter, exshift, shiftInfo);
266  iter = iter + 1;
267  if (iter > m_maxIterations) break;
268  Index im;
269  initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
270  performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
271  }
272  }
273 
274  if(iter <= m_maxIterations)
275  m_info = Success;
276  else
277  m_info = NoConvergence;
278 
279  m_isInitialized = true;
280  m_matUisUptodate = computeU;
281  return *this;
282 }
283 
285 template<typename MatrixType>
286 inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
287 {
288  const Index size = m_matT.cols();
289  // FIXME to be efficient the following would requires a triangular reduxion code
290  // Scalar norm = m_matT.upper().cwiseAbs().sum()
291  // + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
292  Scalar norm(0);
293  for (Index j = 0; j < size; ++j)
294  norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
295  return norm;
296 }
297 
299 template<typename MatrixType>
300 inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, Scalar norm)
301 {
302  Index res = iu;
303  while (res > 0)
304  {
305  Scalar s = internal::abs(m_matT.coeff(res-1,res-1)) + internal::abs(m_matT.coeff(res,res));
306  if (s == 0.0)
307  s = norm;
308  if (internal::abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
309  break;
310  res--;
311  }
312  return res;
313 }
314 
316 template<typename MatrixType>
317 inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scalar exshift)
318 {
319  const Index size = m_matT.cols();
320 
321  // The eigenvalues of the 2x2 matrix [a b; c d] are
322  // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
323  Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
324  Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); // q = tr^2 / 4 - det = discr/4
325  m_matT.coeffRef(iu,iu) += exshift;
326  m_matT.coeffRef(iu-1,iu-1) += exshift;
327 
328  if (q >= Scalar(0)) // Two real eigenvalues
329  {
330  Scalar z = internal::sqrt(internal::abs(q));
331  JacobiRotation<Scalar> rot;
332  if (p >= Scalar(0))
333  rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
334  else
335  rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
336 
337  m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
338  m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
339  m_matT.coeffRef(iu, iu-1) = Scalar(0);
340  if (computeU)
341  m_matU.applyOnTheRight(iu-1, iu, rot);
342  }
343 
344  if (iu > 1)
345  m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
346 }
347 
349 template<typename MatrixType>
350 inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
351 {
352  shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
353  shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
354  shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
355 
356  // Wilkinson's original ad hoc shift
357  if (iter == 10)
358  {
359  exshift += shiftInfo.coeff(0);
360  for (Index i = 0; i <= iu; ++i)
361  m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
362  Scalar s = internal::abs(m_matT.coeff(iu,iu-1)) + internal::abs(m_matT.coeff(iu-1,iu-2));
363  shiftInfo.coeffRef(0) = Scalar(0.75) * s;
364  shiftInfo.coeffRef(1) = Scalar(0.75) * s;
365  shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
366  }
367 
368  // MATLAB's new ad hoc shift
369  if (iter == 30)
370  {
371  Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
372  s = s * s + shiftInfo.coeff(2);
373  if (s > Scalar(0))
374  {
375  s = internal::sqrt(s);
376  if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
377  s = -s;
378  s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
379  s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
380  exshift += s;
381  for (Index i = 0; i <= iu; ++i)
382  m_matT.coeffRef(i,i) -= s;
383  shiftInfo.setConstant(Scalar(0.964));
384  }
385  }
386 }
387 
389 template<typename MatrixType>
390 inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
391 {
392  Vector3s& v = firstHouseholderVector; // alias to save typing
393 
394  for (im = iu-2; im >= il; --im)
395  {
396  const Scalar Tmm = m_matT.coeff(im,im);
397  const Scalar r = shiftInfo.coeff(0) - Tmm;
398  const Scalar s = shiftInfo.coeff(1) - Tmm;
399  v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
400  v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
401  v.coeffRef(2) = m_matT.coeff(im+2,im+1);
402  if (im == il) {
403  break;
404  }
405  const Scalar lhs = m_matT.coeff(im,im-1) * (internal::abs(v.coeff(1)) + internal::abs(v.coeff(2)));
406  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)));
407  if (internal::abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
408  {
409  break;
410  }
411  }
412 }
413 
415 template<typename MatrixType>
416 inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
417 {
418  assert(im >= il);
419  assert(im <= iu-2);
420 
421  const Index size = m_matT.cols();
422 
423  for (Index k = im; k <= iu-2; ++k)
424  {
425  bool firstIteration = (k == im);
426 
427  Vector3s v;
428  if (firstIteration)
429  v = firstHouseholderVector;
430  else
431  v = m_matT.template block<3,1>(k,k-1);
432 
433  Scalar tau, beta;
434  Matrix<Scalar, 2, 1> ess;
435  v.makeHouseholder(ess, tau, beta);
436 
437  if (beta != Scalar(0)) // if v is not zero
438  {
439  if (firstIteration && k > il)
440  m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
441  else if (!firstIteration)
442  m_matT.coeffRef(k,k-1) = beta;
443 
444  // These Householder transformations form the O(n^3) part of the algorithm
445  m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
446  m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
447  if (computeU)
448  m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
449  }
450  }
451 
452  Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
453  Scalar tau, beta;
454  Matrix<Scalar, 1, 1> ess;
455  v.makeHouseholder(ess, tau, beta);
456 
457  if (beta != Scalar(0)) // if v is not zero
458  {
459  m_matT.coeffRef(iu-1, iu-2) = beta;
460  m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
461  m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
462  if (computeU)
463  m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
464  }
465 
466  // clean up pollution due to round-off errors
467  for (Index i = im+2; i <= iu; ++i)
468  {
469  m_matT.coeffRef(i,i-2) = Scalar(0);
470  if (i > im+2)
471  m_matT.coeffRef(i,i-3) = Scalar(0);
472  }
473 }
474 
475 } // end namespace Eigen
476 
477 #endif // EIGEN_REAL_SCHUR_H