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 / +64 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
          VectorType v = VectorType::Zero(A.cols());
327
          v.coeffRef(i) = 1.0;
328
          if (internal::isApprox(fabs(m_betas.coeffRef(i)), 0.0)) {
329
            // Singular eigenvalue, do nothing more
330
          } else {
331
            const Scalar a = real(m_alphas.coeffRef(i));
332
            const Scalar b = m_betas.coeffRef(i);
333
            for (Index j = i-1; j >= 0; j--) {
334
              const Index st = j+1;
335
              const Index sz = i-j; 
336
              v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(b*mS.block(j,st,1,sz) - a*mT.block(j,st,1,sz)).sum() / (b*mS.coeffRef(j, j) - a*mT.coeffRef(j,j));
337
            }
338
          }
339
          m_eivec.col(i).real() = (mZT * v).normalized();
340
          m_eivec.col(i).imag().setConstant(0.);
341
        }
326
        ++i;
342
        ++i;
327
      }
343
      } else {
328
      else
344
        Scalar p = Scalar(0.5) * (mS.coeffRef(i, i) - mS.coeffRef(i+1, i+1));
329
      {
345
        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));
346
        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)));
347
        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
348
335
        m_betas.coeffRef(i)   = m_realQZ.matrixT().coeff(i,i);
349
        m_betas.coeffRef(i)   = mT.coeffRef(i,i);
336
        m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i);
350
        m_betas.coeffRef(i+1) = mT.coeffRef(i,i);
337
        i += 2;
351
        i += 2;
338
      }
352
      }
339
    }
353
    }
354
    m_isInitialized = true;
340
  }
355
  }
341
342
  m_isInitialized = true;
343
  m_eigenvectorsOk = false;//computeEigenvectors;
344
345
  return *this;
356
  return *this;
346
}
357
}
347
358
(-)a/test/eigensolver_generalized_real.cpp (-2 / +10 lines)
Lines 22-28 Link Here
22
22
23
  typedef typename MatrixType::Scalar Scalar;
23
  typedef typename MatrixType::Scalar Scalar;
24
  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
24
  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
25
25
  typedef GeneralizedEigenSolver<MatrixType> SolverType;
26
  
26
  MatrixType a = MatrixType::Random(rows,cols);
27
  MatrixType a = MatrixType::Random(rows,cols);
27
  MatrixType b = MatrixType::Random(rows,cols);
28
  MatrixType b = MatrixType::Random(rows,cols);
28
  MatrixType a1 = MatrixType::Random(rows,cols);
29
  MatrixType a1 = MatrixType::Random(rows,cols);
Lines 32-38 Link Here
32
33
33
  // lets compare to GeneralizedSelfAdjointEigenSolver
34
  // lets compare to GeneralizedSelfAdjointEigenSolver
34
  GeneralizedSelfAdjointEigenSolver<MatrixType> symmEig(spdA, spdB);
35
  GeneralizedSelfAdjointEigenSolver<MatrixType> symmEig(spdA, spdB);
35
  GeneralizedEigenSolver<MatrixType> eig(spdA, spdB);
36
  SolverType eig(spdA, spdB);
36
37
37
  VERIFY_IS_EQUAL(eig.eigenvalues().imag().cwiseAbs().maxCoeff(), 0);
38
  VERIFY_IS_EQUAL(eig.eigenvalues().imag().cwiseAbs().maxCoeff(), 0);
38
39
Lines 47-52 Link Here
47
    GeneralizedEigenSolver<MatrixType> eig2(a.adjoint() * a,b.adjoint() * b);
48
    GeneralizedEigenSolver<MatrixType> eig2(a.adjoint() * a,b.adjoint() * b);
48
    eig2.compute(a.adjoint() * a,b.adjoint() * b);
49
    eig2.compute(a.adjoint() * a,b.adjoint() * b);
49
  }
50
  }
51
  
52
  // Verify that eigenvectors match eigenvalues
53
  realEigenvalues = eig.eigenvalues().real(); // Need to get them again to undo sorting
54
  MatrixType realEigenvectors = eig.eigenvectors().real();
55
  MatrixType LHS = spdA * realEigenvectors;
56
  MatrixType RHS = (spdB * realEigenvectors).array().rowwise() * realEigenvalues.transpose().array();
57
  VERIFY_IS_APPROX(LHS, RHS);
50
}
58
}
51
59
52
void test_eigensolver_generalized_real()
60
void test_eigensolver_generalized_real()

Return to bug 645