ComplexEigenSolver.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_EIGEN_SOLVER_H
28 #define EIGEN_COMPLEX_EIGEN_SOLVER_H
29 
30 #include "./ComplexSchur.h"
31 
32 namespace Eigen {
33 
60 template<typename _MatrixType> class ComplexEigenSolver
61 {
62  public:
63 
65  typedef _MatrixType MatrixType;
66 
67  enum {
72  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
73  };
74 
76  typedef typename MatrixType::Scalar Scalar;
78  typedef typename MatrixType::Index Index;
79 
86  typedef std::complex<RealScalar> ComplexScalar;
87 
94 
101 
108  : m_eivec(),
109  m_eivalues(),
110  m_schur(),
111  m_isInitialized(false),
112  m_eigenvectorsOk(false),
113  m_matX()
114  {}
115 
123  : m_eivec(size, size),
124  m_eivalues(size),
125  m_schur(size),
126  m_isInitialized(false),
127  m_eigenvectorsOk(false),
128  m_matX(size, size)
129  {}
130 
140  ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
141  : m_eivec(matrix.rows(),matrix.cols()),
142  m_eivalues(matrix.cols()),
143  m_schur(matrix.rows()),
144  m_isInitialized(false),
145  m_eigenvectorsOk(false),
146  m_matX(matrix.rows(),matrix.cols())
147  {
148  compute(matrix, computeEigenvectors);
149  }
150 
172  {
173  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
174  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
175  return m_eivec;
176  }
177 
197  {
198  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
199  return m_eivalues;
200  }
201 
226  ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
227 
233  {
234  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
235  return m_schur.info();
236  }
237 
238  protected:
245 
246  private:
247  void doComputeEigenvectors(RealScalar matrixnorm);
248  void sortEigenvalues(bool computeEigenvectors);
249 };
250 
251 
252 template<typename MatrixType>
254 {
255  // this code is inspired from Jampack
256  assert(matrix.cols() == matrix.rows());
257 
258  // Do a complex Schur decomposition, A = U T U^*
259  // The eigenvalues are on the diagonal of T.
260  m_schur.compute(matrix, computeEigenvectors);
261 
262  if(m_schur.info() == Success)
263  {
264  m_eivalues = m_schur.matrixT().diagonal();
265  if(computeEigenvectors)
266  doComputeEigenvectors(matrix.norm());
267  sortEigenvalues(computeEigenvectors);
268  }
269 
270  m_isInitialized = true;
271  m_eigenvectorsOk = computeEigenvectors;
272  return *this;
273 }
274 
275 
276 template<typename MatrixType>
278 {
279  const Index n = m_eivalues.size();
280 
281  // Compute X such that T = X D X^(-1), where D is the diagonal of T.
282  // The matrix X is unit triangular.
283  m_matX = EigenvectorType::Zero(n, n);
284  for(Index k=n-1 ; k>=0 ; k--)
285  {
286  m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
287  // Compute X(i,k) using the (i,k) entry of the equation X T = D X
288  for(Index i=k-1 ; i>=0 ; i--)
289  {
290  m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
291  if(k-i-1>0)
292  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();
293  ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
294  if(z==ComplexScalar(0))
295  {
296  // If the i-th and k-th eigenvalue are equal, then z equals 0.
297  // Use a small value instead, to prevent division by zero.
299  }
300  m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
301  }
302  }
303 
304  // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
305  m_eivec.noalias() = m_schur.matrixU() * m_matX;
306  // .. and normalize the eigenvectors
307  for(Index k=0 ; k<n ; k++)
308  {
309  m_eivec.col(k).normalize();
310  }
311 }
312 
313 
314 template<typename MatrixType>
315 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
316 {
317  const Index n = m_eivalues.size();
318  for (Index i=0; i<n; i++)
319  {
320  Index k;
321  m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
322  if (k != 0)
323  {
324  k += i;
325  std::swap(m_eivalues[k],m_eivalues[i]);
326  if(computeEigenvectors)
327  m_eivec.col(i).swap(m_eivec.col(k));
328  }
329  }
330 }
331 
332 } // end namespace Eigen
333 
334 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H