ComplexEigenSolver.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_EIGEN_SOLVER_H
13 #define EIGEN_COMPLEX_EIGEN_SOLVER_H
14 
15 #include "./ComplexSchur.h"
16 
17 namespace Eigen {
18 
45 template<typename _MatrixType> class ComplexEigenSolver
46 {
47  public:
48 
50  typedef _MatrixType MatrixType;
51 
52  enum {
53  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
54  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
55  Options = MatrixType::Options,
56  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
57  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
58  };
59 
61  typedef typename MatrixType::Scalar Scalar;
62  typedef typename NumTraits<Scalar>::Real RealScalar;
63  typedef typename MatrixType::Index Index;
64 
71  typedef std::complex<RealScalar> ComplexScalar;
72 
79 
86 
93  : m_eivec(),
94  m_eivalues(),
95  m_schur(),
96  m_isInitialized(false),
97  m_eigenvectorsOk(false),
98  m_matX()
99  {}
100 
107  ComplexEigenSolver(Index size)
108  : m_eivec(size, size),
109  m_eivalues(size),
110  m_schur(size),
111  m_isInitialized(false),
112  m_eigenvectorsOk(false),
113  m_matX(size, size)
114  {}
115 
125  ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
126  : m_eivec(matrix.rows(),matrix.cols()),
127  m_eivalues(matrix.cols()),
128  m_schur(matrix.rows()),
129  m_isInitialized(false),
130  m_eigenvectorsOk(false),
131  m_matX(matrix.rows(),matrix.cols())
132  {
133  compute(matrix, computeEigenvectors);
134  }
135 
157  {
158  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
159  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
160  return m_eivec;
161  }
162 
182  {
183  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
184  return m_eivalues;
185  }
186 
211  ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
212 
218  {
219  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
220  return m_schur.info();
221  }
222 
223  protected:
224  EigenvectorType m_eivec;
225  EigenvalueType m_eivalues;
226  ComplexSchur<MatrixType> m_schur;
227  bool m_isInitialized;
228  bool m_eigenvectorsOk;
229  EigenvectorType m_matX;
230 
231  private:
232  void doComputeEigenvectors(RealScalar matrixnorm);
233  void sortEigenvalues(bool computeEigenvectors);
234 };
235 
236 
237 template<typename MatrixType>
239 {
240  // this code is inspired from Jampack
241  assert(matrix.cols() == matrix.rows());
242 
243  // Do a complex Schur decomposition, A = U T U^*
244  // The eigenvalues are on the diagonal of T.
245  m_schur.compute(matrix, computeEigenvectors);
246 
247  if(m_schur.info() == Success)
248  {
249  m_eivalues = m_schur.matrixT().diagonal();
250  if(computeEigenvectors)
251  doComputeEigenvectors(matrix.norm());
252  sortEigenvalues(computeEigenvectors);
253  }
254 
255  m_isInitialized = true;
256  m_eigenvectorsOk = computeEigenvectors;
257  return *this;
258 }
259 
260 
261 template<typename MatrixType>
263 {
264  const Index n = m_eivalues.size();
265 
266  // Compute X such that T = X D X^(-1), where D is the diagonal of T.
267  // The matrix X is unit triangular.
268  m_matX = EigenvectorType::Zero(n, n);
269  for(Index k=n-1 ; k>=0 ; k--)
270  {
271  m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
272  // Compute X(i,k) using the (i,k) entry of the equation X T = D X
273  for(Index i=k-1 ; i>=0 ; i--)
274  {
275  m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
276  if(k-i-1>0)
277  m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
278  ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
279  if(z==ComplexScalar(0))
280  {
281  // If the i-th and k-th eigenvalue are equal, then z equals 0.
282  // Use a small value instead, to prevent division by zero.
283  internal::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
284  }
285  m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
286  }
287  }
288 
289  // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
290  m_eivec.noalias() = m_schur.matrixU() * m_matX;
291  // .. and normalize the eigenvectors
292  for(Index k=0 ; k<n ; k++)
293  {
294  m_eivec.col(k).normalize();
295  }
296 }
297 
298 
299 template<typename MatrixType>
300 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
301 {
302  const Index n = m_eivalues.size();
303  for (Index i=0; i<n; i++)
304  {
305  Index k;
306  m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
307  if (k != 0)
308  {
309  k += i;
310  std::swap(m_eivalues[k],m_eivalues[i]);
311  if(computeEigenvectors)
312  m_eivec.col(i).swap(m_eivec.col(k));
313  }
314  }
315 }
316 
317 } // end namespace Eigen
318 
319 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H