ComplexSchur.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) 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 // Eigen is free software; you can redistribute it and/or
9 // modify it under the terms of the GNU Lesser General Public
10 // License as published by the Free Software Foundation; either
11 // version 3 of the License, or (at your option) any later version.
12 //
13 // Alternatively, you can redistribute it and/or
14 // modify it under the terms of the GNU General Public License as
15 // published by the Free Software Foundation; either version 2 of
16 // the License, or (at your option) any later version.
17 //
18 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
19 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
21 // GNU General Public License for more details.
22 //
23 // You should have received a copy of the GNU Lesser General Public
24 // License and a copy of the GNU General Public License along with
25 // Eigen. If not, see <http://www.gnu.org/licenses/>.
26 
27 #ifndef EIGEN_COMPLEX_SCHUR_H
28 #define EIGEN_COMPLEX_SCHUR_H
29 
31 
32 namespace Eigen {
33 
34 namespace internal {
35 template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
36 }
37 
66 template<typename _MatrixType> class ComplexSchur
67 {
68  public:
69  typedef _MatrixType MatrixType;
70  enum {
75  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
76  };
77 
79  typedef typename MatrixType::Scalar Scalar;
81  typedef typename MatrixType::Index Index;
82 
89  typedef std::complex<RealScalar> ComplexScalar;
90 
97 
110  : m_matT(size,size),
111  m_matU(size,size),
112  m_hess(size),
113  m_isInitialized(false),
114  m_matUisUptodate(false)
115  {}
116 
126  ComplexSchur(const MatrixType& matrix, bool computeU = true)
127  : m_matT(matrix.rows(),matrix.cols()),
128  m_matU(matrix.rows(),matrix.cols()),
129  m_hess(matrix.rows()),
130  m_isInitialized(false),
131  m_matUisUptodate(false)
132  {
133  compute(matrix, computeU);
134  }
135 
150  const ComplexMatrixType& matrixU() const
151  {
152  eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
153  eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
154  return m_matU;
155  }
156 
174  const ComplexMatrixType& matrixT() const
175  {
176  eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
177  return m_matT;
178  }
179 
199  ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
200 
206  {
207  eigen_assert(m_isInitialized && "RealSchur is not initialized.");
208  return m_info;
209  }
210 
215  static const int m_maxIterations = 30;
216 
217  protected:
223 
224  private:
225  bool subdiagonalEntryIsNeglegible(Index i);
226  ComplexScalar computeShift(Index iu, Index iter);
227  void reduceToTriangularForm(bool computeU);
228  friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
229 };
230 
234 template<typename MatrixType>
235 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
236 {
237  RealScalar d = internal::norm1(m_matT.coeff(i,i)) + internal::norm1(m_matT.coeff(i+1,i+1));
238  RealScalar sd = internal::norm1(m_matT.coeff(i+1,i));
240  {
241  m_matT.coeffRef(i+1,i) = ComplexScalar(0);
242  return true;
243  }
244  return false;
245 }
246 
247 
249 template<typename MatrixType>
250 typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
251 {
252  if (iter == 10 || iter == 20)
253  {
254  // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
256  }
257 
258  // compute the shift as one of the eigenvalues of t, the 2x2
259  // diagonal block on the bottom of the active submatrix
260  Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
261  RealScalar normt = t.cwiseAbs().sum();
262  t /= normt; // the normalization by sf is to avoid under/overflow
263 
264  ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
265  ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
266  ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
267  ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
268  ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
269  ComplexScalar eival1 = (trace + disc) / RealScalar(2);
270  ComplexScalar eival2 = (trace - disc) / RealScalar(2);
271 
272  if(internal::norm1(eival1) > internal::norm1(eival2))
273  eival2 = det / eival1;
274  else
275  eival1 = det / eival2;
276 
277  // choose the eigenvalue closest to the bottom entry of the diagonal
278  if(internal::norm1(eival1-t.coeff(1,1)) < internal::norm1(eival2-t.coeff(1,1)))
279  return normt * eival1;
280  else
281  return normt * eival2;
282 }
283 
284 
285 template<typename MatrixType>
287 {
288  m_matUisUptodate = false;
289  eigen_assert(matrix.cols() == matrix.rows());
290 
291  if(matrix.cols() == 1)
292  {
293  m_matT = matrix.template cast<ComplexScalar>();
294  if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
295  m_info = Success;
296  m_isInitialized = true;
297  m_matUisUptodate = computeU;
298  return *this;
299  }
300 
302  reduceToTriangularForm(computeU);
303  return *this;
304 }
305 
306 namespace internal {
307 
308 /* Reduce given matrix to Hessenberg form */
309 template<typename MatrixType, bool IsComplex>
310 struct complex_schur_reduce_to_hessenberg
311 {
312  // this is the implementation for the case IsComplex = true
313  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
314  {
315  _this.m_hess.compute(matrix);
316  _this.m_matT = _this.m_hess.matrixH();
317  if(computeU) _this.m_matU = _this.m_hess.matrixQ();
318  }
319 };
320 
321 template<typename MatrixType>
322 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
323 {
324  static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
325  {
328 
329  // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
330  _this.m_hess.compute(matrix);
331  _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
332  if(computeU)
333  {
334  // This may cause an allocation which seems to be avoidable
335  MatrixType Q = _this.m_hess.matrixQ();
336  _this.m_matU = Q.template cast<ComplexScalar>();
337  }
338  }
339 };
340 
341 } // end namespace internal
342 
343 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
344 template<typename MatrixType>
345 void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
346 {
347  // The matrix m_matT is divided in three parts.
348  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
349  // Rows il,...,iu is the part we are working on (the active submatrix).
350  // Rows iu+1,...,end are already brought in triangular form.
351  Index iu = m_matT.cols() - 1;
352  Index il;
353  Index iter = 0; // number of iterations we are working on the (iu,iu) element
354 
355  while(true)
356  {
357  // find iu, the bottom row of the active submatrix
358  while(iu > 0)
359  {
360  if(!subdiagonalEntryIsNeglegible(iu-1)) break;
361  iter = 0;
362  --iu;
363  }
364 
365  // if iu is zero then we are done; the whole matrix is triangularized
366  if(iu==0) break;
367 
368  // if we spent too many iterations on the current element, we give up
369  iter++;
370  if(iter > m_maxIterations) break;
371 
372  // find il, the top row of the active submatrix
373  il = iu-1;
374  while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
375  {
376  --il;
377  }
378 
379  /* perform the QR step using Givens rotations. The first rotation
380  creates a bulge; the (il+2,il) element becomes nonzero. This
381  bulge is chased down to the bottom of the active submatrix. */
382 
383  ComplexScalar shift = computeShift(iu, iter);
384  JacobiRotation<ComplexScalar> rot;
385  rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
386  m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
387  m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
388  if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
389 
390  for(Index i=il+1 ; i<iu ; i++)
391  {
392  rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
393  m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
394  m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
395  m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
396  if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
397  }
398  }
399 
400  if(iter <= m_maxIterations)
401  m_info = Success;
402  else
404 
405  m_isInitialized = true;
406  m_matUisUptodate = computeU;
407 }
408 
409 } // end namespace Eigen
410 
411 #endif // EIGEN_COMPLEX_SCHUR_H