This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
View | Details | Raw Unified | Return to bug 645 | Differences between
and this patch

Collapse All | Expand All

(-)a/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h (-53 / +65 lines)
Lines 12-17 Link Here
12
#define EIGEN_GENERALIZEDEIGENSOLVER_H
12
#define EIGEN_GENERALIZEDEIGENSOLVER_H
13
13
14
#include "./RealQZ.h"
14
#include "./RealQZ.h"
15
#include <iostream>
15
16
16
namespace Eigen { 
17
namespace Eigen { 
17
18
Lines 114-120 Link Here
114
      *
115
      *
115
      * \sa compute() for an example.
116
      * \sa compute() for an example.
116
      */
117
      */
117
    GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {}
118
    GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ() {}
118
119
119
    /** \brief Default constructor with memory preallocation
120
    /** \brief Default constructor with memory preallocation
120
      *
121
      *
Lines 127-136 Link Here
127
        m_alphas(size),
128
        m_alphas(size),
128
        m_betas(size),
129
        m_betas(size),
129
        m_isInitialized(false),
130
        m_isInitialized(false),
130
        m_eigenvectorsOk(false),
131
        m_realQZ(size)
131
        m_realQZ(size),
132
        m_matS(size, size),
133
        m_tmp(size)
134
    {}
132
    {}
135
133
136
    /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair.
134
    /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair.
Lines 150-159 Link Here
150
        m_alphas(A.cols()),
148
        m_alphas(A.cols()),
151
        m_betas(A.cols()),
149
        m_betas(A.cols()),
152
        m_isInitialized(false),
150
        m_isInitialized(false),
153
        m_eigenvectorsOk(false),
151
        m_realQZ(A.cols())
154
        m_realQZ(A.cols()),
155
        m_matS(A.rows(), A.cols()),
156
        m_tmp(A.cols())
157
    {
152
    {
158
      compute(A, B, computeEigenvectors);
153
      compute(A, B, computeEigenvectors);
159
    }
154
    }
Lines 201-206 Link Here
201
      return EigenvalueType(m_alphas,m_betas);
196
      return EigenvalueType(m_alphas,m_betas);
202
    }
197
    }
203
198
199
    /** \brief Returns the computed generalized eigenvectors.
200
      *
201
      * \returns Returns a matrix whose columns are the (possibly complex) right eigenvectors
202
      * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues.
203
      *
204
      * \pre Either the constructor 
205
      * GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function
206
      * compute(const MatrixType&,const MatrixType&,bool) has been called before.
207
      *
208
      * \sa eigenvalues()
209
      */
210
    EigenvectorsType eigenvectors() const
211
    {
212
      eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
213
      return m_eivec;
214
    }
215
204
    /** \returns A const reference to the vectors containing the alpha values
216
    /** \returns A const reference to the vectors containing the alpha values
205
      *
217
      *
206
      * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
218
      * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j).
Lines 270-298 Link Here
270
      EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
282
      EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
271
    }
283
    }
272
    
284
    
273
    MatrixType m_eivec;
285
    EigenvectorsType  m_eivec;
274
    ComplexVectorType m_alphas;
286
    ComplexVectorType m_alphas;
275
    VectorType m_betas;
287
    VectorType        m_betas;
276
    bool m_isInitialized;
288
    bool m_isInitialized;
277
    bool m_eigenvectorsOk;
278
    RealQZ<MatrixType> m_realQZ;
289
    RealQZ<MatrixType> m_realQZ;
279
    MatrixType m_matS;
280
281
    typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
282
    ColumnVectorType m_tmp;
283
};
290
};
284
291
285
//template<typename MatrixType>
286
//typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const
287
//{
288
//  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
289
//  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
290
//  Index n = m_eivec.cols();
291
//  EigenvectorsType matV(n,n);
292
//  // TODO
293
//  return matV;
294
//}
295
296
template<typename MatrixType>
292
template<typename MatrixType>
297
GeneralizedEigenSolver<MatrixType>&
293
GeneralizedEigenSolver<MatrixType>&
298
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
294
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
Lines 307-347 Link Here
307
  // A = Q S Z and B = Q T Z
303
  // A = Q S Z and B = Q T Z
308
  m_realQZ.compute(A, B, computeEigenvectors);
304
  m_realQZ.compute(A, B, computeEigenvectors);
309
305
310
  if (m_realQZ.info() == Success)
306
  if (m_realQZ.info() == Success) {
311
  {
307
    // Resize storage
312
    m_matS = m_realQZ.matrixS();
313
    if (computeEigenvectors)
314
      m_eivec = m_realQZ.matrixZ().transpose();
315
  
316
    // Compute eigenvalues from matS
317
    m_alphas.resize(A.cols());
308
    m_alphas.resize(A.cols());
318
    m_betas.resize(A.cols());
309
    m_betas.resize(A.cols());
310
    if (computeEigenvectors) {
311
      m_eivec.resize(A.cols(), A.cols());
312
    }
313
    
314
    // Grab some references
315
    const MatrixType &mZT = m_realQZ.matrixZ().transpose();
316
    const MatrixType &mS = m_realQZ.matrixS();
317
    const MatrixType &mT = m_realQZ.matrixT();
318
319
    Index i = 0;
319
    Index i = 0;
320
    while (i < A.cols())
320
    while (i < A.cols()) {
321
    {
321
      if (i == A.cols() - 1 || mS.coeff(i+1, i) == Scalar(0)) {
322
      if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0))
322
        // Real eigenvalue
323
      {
323
        m_alphas.coeffRef(i) = mS.coeff(i, i);
324
        m_alphas.coeffRef(i) = m_matS.coeff(i, i);
324
        m_betas.coeffRef(i)  = mT.coeff(i,i);
325
        m_betas.coeffRef(i)  = m_realQZ.matrixT().coeff(i,i);
325
        if (computeEigenvectors) {
326
          VectorXd v = VectorXd::Zero(A.cols());
327
          v.coeffRef(i) = 1.0;
328
          if (fabs(m_betas.coeffRef(i)) < NumTraits<Scalar>::epsilon()) {
329
            // Singular eigenvalue, do nothing more
330
          } else {
331
            Scalar lambda = real(m_alphas.coeffRef(i)) / m_betas.coeffRef(i);
332
            MatrixType mC = mS - lambda*mT;
333
            for (Index j = i-1; j >= 0; j--) {
334
              for (Index k = j+1; k <= i; k++) {
335
                v.coeffRef(j) -= mC.coeffRef(j, k) * v.coeffRef(k);
336
              }
337
              v.coeffRef(j) /= mC.coeffRef(j, j);
338
            }
339
          }
340
          m_eivec.col(i).real() = (mZT * v).normalized();
341
          m_eivec.col(i).imag().setConstant(0.);
342
        }
326
        ++i;
343
        ++i;
327
      }
344
      } else {
328
      else
345
        Scalar p = Scalar(0.5) * (mS.coeffRef(i, i) - mS.coeffRef(i+1, i+1));
329
      {
346
        Scalar z = sqrt(abs(p * p + mS.coeffRef(i+1, i) * mS.coeffRef(i, i+1)));
330
        Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1));
347
        m_alphas.coeffRef(i)   = ComplexScalar(mS.coeffRef(i+1, i+1) + p, z);
331
        Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1)));
348
        m_alphas.coeffRef(i+1) = ComplexScalar(mS.coeffRef(i+1, i+1) + p, -z);
332
        m_alphas.coeffRef(i)   = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z);
333
        m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z);
334
349
335
        m_betas.coeffRef(i)   = m_realQZ.matrixT().coeff(i,i);
350
        m_betas.coeffRef(i)   = mT.coeffRef(i,i);
336
        m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i);
351
        m_betas.coeffRef(i+1) = mT.coeffRef(i,i);
337
        i += 2;
352
        i += 2;
338
      }
353
      }
339
    }
354
    }
355
    m_isInitialized = true;
340
  }
356
  }
341
342
  m_isInitialized = true;
343
  m_eigenvectorsOk = false;//computeEigenvectors;
344
345
  return *this;
357
  return *this;
346
}
358
}
347
359

Return to bug 645