ComplexSchur.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
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>
7 //
8 // This Source Code Form is subject to the terms of the Mozilla
9 // Public License v. 2.0. If a copy of the MPL was not distributed
10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11 
12 #ifndef EIGEN_COMPLEX_SCHUR_H
13 #define EIGEN_COMPLEX_SCHUR_H
14 
15 #include "./HessenbergDecomposition.h"
16 
17 namespace Eigen {
18 
19 namespace internal {
20 template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
21 }
22 
51 template<typename _MatrixType> class ComplexSchur
52 {
53  public:
54  typedef _MatrixType MatrixType;
55  enum {
56  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
57  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
58  Options = MatrixType::Options,
59  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
60  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
61  };
62 
64  typedef typename MatrixType::Scalar Scalar;
65  typedef typename NumTraits<Scalar>::Real RealScalar;
66  typedef typename MatrixType::Index Index;
67 
74  typedef std::complex<RealScalar> ComplexScalar;
75 
82 
94  ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
95  : m_matT(size,size),
96  m_matU(size,size),
97  m_hess(size),
98  m_isInitialized(false),
99  m_matUisUptodate(false)
100  {}
101 
111  ComplexSchur(const MatrixType& matrix, bool computeU = true)
112  : m_matT(matrix.rows(),matrix.cols()),
113  m_matU(matrix.rows(),matrix.cols()),
114  m_hess(matrix.rows()),
115  m_isInitialized(false),
116  m_matUisUptodate(false)
117  {
118  compute(matrix, computeU);
119  }
120 
135  const ComplexMatrixType& matrixU() const
136  {
137  eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
138  eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
139  return m_matU;
140  }
141 
159  const ComplexMatrixType& matrixT() const
160  {
161  eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
162  return m_matT;
163  }
164 
184  ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
185 
191  {
192  eigen_assert(m_isInitialized && "RealSchur is not initialized.");
193  return m_info;
194  }
195 
200  static const int m_maxIterations = 30;
201 
202  protected:
203  ComplexMatrixType m_matT, m_matU;
205  ComputationInfo m_info;
206  bool m_isInitialized;
207  bool m_matUisUptodate;
208 
209  private:
210  bool subdiagonalEntryIsNeglegible(Index i);
211  ComplexScalar computeShift(Index iu, Index iter);
212  void reduceToTriangularForm(bool computeU);
213  friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
214 };
215 
219 template<typename MatrixType>
220 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
221 {
222  RealScalar d = internal::norm1(m_matT.coeff(i,i)) + internal::norm1(m_matT.coeff(i+1,i+1));
223  RealScalar sd = internal::norm1(m_matT.coeff(i+1,i));
224  if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
225  {
226  m_matT.coeffRef(i+1,i) = ComplexScalar(0);
227  return true;
228  }
229  return false;
230 }
231 
232 
234 template<typename MatrixType>
235 typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
236 {
237  if (iter == 10 || iter == 20)
238  {
239  // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
240  return internal::abs(internal::real(m_matT.coeff(iu,iu-1))) + internal::abs(internal::real(m_matT.coeff(iu-1,iu-2)));
241  }
242 
243  // compute the shift as one of the eigenvalues of t, the 2x2
244  // diagonal block on the bottom of the active submatrix
245  Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
246  RealScalar normt = t.cwiseAbs().sum();
247  t /= normt; // the normalization by sf is to avoid under/overflow
248 
249  ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
250  ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
251  ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
252  ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
253  ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
254  ComplexScalar eival1 = (trace + disc) / RealScalar(2);
255  ComplexScalar eival2 = (trace - disc) / RealScalar(2);
256 
257  if(internal::norm1(eival1) > internal::norm1(eival2))
258  eival2 = det / eival1;
259  else
260  eival1 = det / eival2;
261 
262  // choose the eigenvalue closest to the bottom entry of the diagonal
263  if(internal::norm1(eival1-t.coeff(1,1)) < internal::norm1(eival2-t.coeff(1,1)))
264  return normt * eival1;
265  else
266  return normt * eival2;
267 }
268 
269 
270 template<typename MatrixType>
271 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
272 {
273  m_matUisUptodate = false;
274  eigen_assert(matrix.cols() == matrix.rows());
275 
276  if(matrix.cols() == 1)
277  {
278  m_matT = matrix.template cast<ComplexScalar>();
279  if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
280  m_info = Success;
281  m_isInitialized = true;
282  m_matUisUptodate = computeU;
283  return *this;
284  }
285 
286  internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
287  reduceToTriangularForm(computeU);
288  return *this;
289 }
290 
291 namespace internal {
292 
293 /* Reduce given matrix to Hessenberg form */
294 template<typename MatrixType, bool IsComplex>
295 struct complex_schur_reduce_to_hessenberg
296 {
297  // this is the implementation for the case IsComplex = true
298  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
299  {
300  _this.m_hess.compute(matrix);
301  _this.m_matT = _this.m_hess.matrixH();
302  if(computeU) _this.m_matU = _this.m_hess.matrixQ();
303  }
304 };
305 
306 template<typename MatrixType>
307 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
308 {
309  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
310  {
313 
314  // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
315  _this.m_hess.compute(matrix);
316  _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
317  if(computeU)
318  {
319  // This may cause an allocation which seems to be avoidable
320  MatrixType Q = _this.m_hess.matrixQ();
321  _this.m_matU = Q.template cast<ComplexScalar>();
322  }
323  }
324 };
325 
326 } // end namespace internal
327 
328 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
329 template<typename MatrixType>
330 void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
331 {
332  // The matrix m_matT is divided in three parts.
333  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
334  // Rows il,...,iu is the part we are working on (the active submatrix).
335  // Rows iu+1,...,end are already brought in triangular form.
336  Index iu = m_matT.cols() - 1;
337  Index il;
338  Index iter = 0; // number of iterations we are working on the (iu,iu) element
339  Index totalIter = 0; // number of iterations for whole matrix
340 
341  while(true)
342  {
343  // find iu, the bottom row of the active submatrix
344  while(iu > 0)
345  {
346  if(!subdiagonalEntryIsNeglegible(iu-1)) break;
347  iter = 0;
348  --iu;
349  }
350 
351  // if iu is zero then we are done; the whole matrix is triangularized
352  if(iu==0) break;
353 
354  // if we spent too many iterations, we give up
355  iter++;
356  totalIter++;
357  if(totalIter > m_maxIterations * m_matT.cols()) break;
358 
359  // find il, the top row of the active submatrix
360  il = iu-1;
361  while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
362  {
363  --il;
364  }
365 
366  /* perform the QR step using Givens rotations. The first rotation
367  creates a bulge; the (il+2,il) element becomes nonzero. This
368  bulge is chased down to the bottom of the active submatrix. */
369 
370  ComplexScalar shift = computeShift(iu, iter);
371  JacobiRotation<ComplexScalar> rot;
372  rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
373  m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
374  m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
375  if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
376 
377  for(Index i=il+1 ; i<iu ; i++)
378  {
379  rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
380  m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
381  m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
382  m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
383  if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
384  }
385  }
386 
387  if(totalIter <= m_maxIterations * m_matT.cols())
388  m_info = Success;
389  else
390  m_info = NoConvergence;
391 
392  m_isInitialized = true;
393  m_matUisUptodate = computeU;
394 }
395 
396 } // end namespace Eigen
397 
398 #endif // EIGEN_COMPLEX_SCHUR_H