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  if(norm!=0)
243  {
244  while (iu >= 0)
245  {
246  Index il = findSmallSubdiagEntry(iu, norm);
247 
248  // Check for convergence
249  if (il == iu) // One root found
250  {
251  m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
252  if (iu > 0)
253  m_matT.coeffRef(iu, iu-1) = Scalar(0);
254  iu--;
255  iter = 0;
256  }
257  else if (il == iu-1) // Two roots found
258  {
259  splitOffTwoRows(iu, computeU, exshift);
260  iu -= 2;
261  iter = 0;
262  }
263  else // No convergence yet
264  {
265  // The firstHouseholderVector vector has to be initialized to something to get rid of a silly GCC warning (-O1 -Wall -DNDEBUG )
266  Vector3s firstHouseholderVector(0,0,0), shiftInfo;
267  computeShift(iu, iter, exshift, shiftInfo);
268  iter = iter + 1;
269  if (iter > m_maxIterations) break;
270  Index im;
271  initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
272  performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
273  }
274  }
275  }
276  if(iter <= m_maxIterations)
277  m_info = Success;
278  else
279  m_info = NoConvergence;
280 
281  m_isInitialized = true;
282  m_matUisUptodate = computeU;
283  return *this;
284 }
285 
287 template<typename MatrixType>
288 inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
289 {
290  const Index size = m_matT.cols();
291  // FIXME to be efficient the following would requires a triangular reduxion code
292  // Scalar norm = m_matT.upper().cwiseAbs().sum()
293  // + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
294  Scalar norm(0);
295  for (Index j = 0; j < size; ++j)
296  norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
297  return norm;
298 }
299 
301 template<typename MatrixType>
302 inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, Scalar norm)
303 {
304  Index res = iu;
305  while (res > 0)
306  {
307  Scalar s = internal::abs(m_matT.coeff(res-1,res-1)) + internal::abs(m_matT.coeff(res,res));
308  if (s == 0.0)
309  s = norm;
310  if (internal::abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
311  break;
312  res--;
313  }
314  return res;
315 }
316 
318 template<typename MatrixType>
319 inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scalar exshift)
320 {
321  const Index size = m_matT.cols();
322 
323  // The eigenvalues of the 2x2 matrix [a b; c d] are
324  // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
325  Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
326  Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu); // q = tr^2 / 4 - det = discr/4
327  m_matT.coeffRef(iu,iu) += exshift;
328  m_matT.coeffRef(iu-1,iu-1) += exshift;
329 
330  if (q >= Scalar(0)) // Two real eigenvalues
331  {
332  Scalar z = internal::sqrt(internal::abs(q));
333  JacobiRotation<Scalar> rot;
334  if (p >= Scalar(0))
335  rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
336  else
337  rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
338 
339  m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
340  m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
341  m_matT.coeffRef(iu, iu-1) = Scalar(0);
342  if (computeU)
343  m_matU.applyOnTheRight(iu-1, iu, rot);
344  }
345 
346  if (iu > 1)
347  m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
348 }
349 
351 template<typename MatrixType>
352 inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
353 {
354  shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
355  shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
356  shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
357 
358  // Wilkinson's original ad hoc shift
359  if (iter == 10)
360  {
361  exshift += shiftInfo.coeff(0);
362  for (Index i = 0; i <= iu; ++i)
363  m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
364  Scalar s = internal::abs(m_matT.coeff(iu,iu-1)) + internal::abs(m_matT.coeff(iu-1,iu-2));
365  shiftInfo.coeffRef(0) = Scalar(0.75) * s;
366  shiftInfo.coeffRef(1) = Scalar(0.75) * s;
367  shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
368  }
369 
370  // MATLAB's new ad hoc shift
371  if (iter == 30)
372  {
373  Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
374  s = s * s + shiftInfo.coeff(2);
375  if (s > Scalar(0))
376  {
377  s = internal::sqrt(s);
378  if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
379  s = -s;
380  s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
381  s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
382  exshift += s;
383  for (Index i = 0; i <= iu; ++i)
384  m_matT.coeffRef(i,i) -= s;
385  shiftInfo.setConstant(Scalar(0.964));
386  }
387  }
388 }
389 
391 template<typename MatrixType>
392 inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
393 {
394  Vector3s& v = firstHouseholderVector; // alias to save typing
395 
396  for (im = iu-2; im >= il; --im)
397  {
398  const Scalar Tmm = m_matT.coeff(im,im);
399  const Scalar r = shiftInfo.coeff(0) - Tmm;
400  const Scalar s = shiftInfo.coeff(1) - Tmm;
401  v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
402  v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
403  v.coeffRef(2) = m_matT.coeff(im+2,im+1);
404  if (im == il) {
405  break;
406  }
407  const Scalar lhs = m_matT.coeff(im,im-1) * (internal::abs(v.coeff(1)) + internal::abs(v.coeff(2)));
408  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)));
409  if (internal::abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
410  {
411  break;
412  }
413  }
414 }
415 
417 template<typename MatrixType>
418 inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
419 {
420  assert(im >= il);
421  assert(im <= iu-2);
422 
423  const Index size = m_matT.cols();
424 
425  for (Index k = im; k <= iu-2; ++k)
426  {
427  bool firstIteration = (k == im);
428 
429  Vector3s v;
430  if (firstIteration)
431  v = firstHouseholderVector;
432  else
433  v = m_matT.template block<3,1>(k,k-1);
434 
435  Scalar tau, beta;
436  Matrix<Scalar, 2, 1> ess;
437  v.makeHouseholder(ess, tau, beta);
438 
439  if (beta != Scalar(0)) // if v is not zero
440  {
441  if (firstIteration && k > il)
442  m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
443  else if (!firstIteration)
444  m_matT.coeffRef(k,k-1) = beta;
445 
446  // These Householder transformations form the O(n^3) part of the algorithm
447  m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
448  m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
449  if (computeU)
450  m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
451  }
452  }
453 
454  Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
455  Scalar tau, beta;
456  Matrix<Scalar, 1, 1> ess;
457  v.makeHouseholder(ess, tau, beta);
458 
459  if (beta != Scalar(0)) // if v is not zero
460  {
461  m_matT.coeffRef(iu-1, iu-2) = beta;
462  m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
463  m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
464  if (computeU)
465  m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
466  }
467 
468  // clean up pollution due to round-off errors
469  for (Index i = im+2; i <= iu; ++i)
470  {
471  m_matT.coeffRef(i,i-2) = Scalar(0);
472  if (i > im+2)
473  m_matT.coeffRef(i,i-3) = Scalar(0);
474  }
475 }
476 
477 } // end namespace Eigen
478 
479 #endif // EIGEN_REAL_SCHUR_H