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 (-73 / +111 lines)
Lines 89-95 Link Here
89
      */
89
      */
90
    typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType;
90
    typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType;
91
91
92
    /** \brief Type for vector of complex scalar values eigenvalues as returned by betas().
92
    /** \brief Type for vector of complex scalar values eigenvalues as returned by alphas().
93
      *
93
      *
94
      * This is a column vector with entries of type #ComplexScalar.
94
      * This is a column vector with entries of type #ComplexScalar.
95
      * The length of the vector is the size of #MatrixType.
95
      * The length of the vector is the size of #MatrixType.
Lines 114-120 Link Here
114
      *
114
      *
115
      * \sa compute() for an example.
115
      * \sa compute() for an example.
116
      */
116
      */
117
    GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {}
117
    GeneralizedEigenSolver()
118
      : m_eivec(),
119
        m_alphas(),
120
        m_betas(),
121
        m_valuesOkay(false),
122
        m_vectorsOkay(false),
123
        m_realQZ()
124
    {}
118
125
119
    /** \brief Default constructor with memory preallocation
126
    /** \brief Default constructor with memory preallocation
120
      *
127
      *
Lines 126-136 Link Here
126
      : m_eivec(size, size),
133
      : m_eivec(size, size),
127
        m_alphas(size),
134
        m_alphas(size),
128
        m_betas(size),
135
        m_betas(size),
129
        m_isInitialized(false),
136
        m_valuesOkay(false),
130
        m_eigenvectorsOk(false),
137
        m_vectorsOkay(false),
131
        m_realQZ(size),
138
        m_realQZ(size)
132
        m_matS(size, size),
133
        m_tmp(size)
134
    {}
139
    {}
135
140
136
    /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair.
141
    /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair.
Lines 149-181 Link Here
149
      : m_eivec(A.rows(), A.cols()),
154
      : m_eivec(A.rows(), A.cols()),
150
        m_alphas(A.cols()),
155
        m_alphas(A.cols()),
151
        m_betas(A.cols()),
156
        m_betas(A.cols()),
152
        m_isInitialized(false),
157
        m_valuesOkay(false),
153
        m_eigenvectorsOk(false),
158
        m_vectorsOkay(false),
154
        m_realQZ(A.cols()),
159
        m_realQZ(A.cols())
155
        m_matS(A.rows(), A.cols()),
156
        m_tmp(A.cols())
157
    {
160
    {
158
      compute(A, B, computeEigenvectors);
161
      compute(A, B, computeEigenvectors);
159
    }
162
    }
160
163
161
    /* \brief Returns the computed generalized eigenvectors.
164
    /* \brief Returns the computed generalized eigenvectors.
162
      *
165
      *
163
      * \returns  %Matrix whose columns are the (possibly complex) eigenvectors.
166
      * \returns  %Matrix whose columns are the (possibly complex) right eigenvectors.
167
      * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues.
164
      *
168
      *
165
      * \pre Either the constructor 
169
      * \pre Either the constructor 
166
      * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
170
      * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
167
      * compute(const MatrixType&, const MatrixType& bool) has been called before, and
171
      * compute(const MatrixType&, const MatrixType& bool) has been called before, and
168
      * \p computeEigenvectors was set to true (the default).
172
      * \p computeEigenvectors was set to true (the default).
169
      *
173
      *
170
      * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
171
      * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
172
      * eigenvectors are normalized to have (Euclidean) norm equal to one. The
173
      * matrix returned by this function is the matrix \f$ V \f$ in the
174
      * generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists.
175
      *
176
      * \sa eigenvalues()
174
      * \sa eigenvalues()
177
      */
175
      */
178
//    EigenvectorsType eigenvectors() const;
176
    EigenvectorsType eigenvectors() const {
177
      eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.")
178
      return m_eivec;
179
    }
179
180
180
    /** \brief Returns an expression of the computed generalized eigenvalues.
181
    /** \brief Returns an expression of the computed generalized eigenvalues.
181
      *
182
      *
Lines 197-203 Link Here
197
      */
198
      */
198
    EigenvalueType eigenvalues() const
199
    EigenvalueType eigenvalues() const
199
    {
200
    {
200
      eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
201
      eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
201
      return EigenvalueType(m_alphas,m_betas);
202
      return EigenvalueType(m_alphas,m_betas);
202
    }
203
    }
203
204
Lines 208-214 Link Here
208
      * \sa betas(), eigenvalues() */
209
      * \sa betas(), eigenvalues() */
209
    ComplexVectorType alphas() const
210
    ComplexVectorType alphas() const
210
    {
211
    {
211
      eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
212
      eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
212
      return m_alphas;
213
      return m_alphas;
213
    }
214
    }
214
215
Lines 219-225 Link Here
219
      * \sa alphas(), eigenvalues() */
220
      * \sa alphas(), eigenvalues() */
220
    VectorType betas() const
221
    VectorType betas() const
221
    {
222
    {
222
      eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
223
      eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
223
      return m_betas;
224
      return m_betas;
224
    }
225
    }
225
226
Lines 250-256 Link Here
250
251
251
    ComputationInfo info() const
252
    ComputationInfo info() const
252
    {
253
    {
253
      eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
254
      eigen_assert(m_valuesOkay && "EigenSolver is not initialized.");
254
      return m_realQZ.info();
255
      return m_realQZ.info();
255
    }
256
    }
256
257
Lines 270-298 Link Here
270
      EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
271
      EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
271
    }
272
    }
272
    
273
    
273
    MatrixType m_eivec;
274
    EigenvectorsType m_eivec;
274
    ComplexVectorType m_alphas;
275
    ComplexVectorType m_alphas;
275
    VectorType m_betas;
276
    VectorType m_betas;
276
    bool m_isInitialized;
277
    bool m_valuesOkay, m_vectorsOkay;
277
    bool m_eigenvectorsOk;
278
    RealQZ<MatrixType> m_realQZ;
278
    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
};
279
};
284
280
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>
281
template<typename MatrixType>
297
GeneralizedEigenSolver<MatrixType>&
282
GeneralizedEigenSolver<MatrixType>&
298
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
283
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
Lines 302-359 Link Here
302
  using std::sqrt;
287
  using std::sqrt;
303
  using std::abs;
288
  using std::abs;
304
  eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
289
  eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
305
290
  m_valuesOkay = false;
291
  m_vectorsOkay = false;
306
  // Reduce to generalized real Schur form:
292
  // Reduce to generalized real Schur form:
307
  // A = Q S Z and B = Q T Z
293
  // A = Q S Z and B = Q T Z
308
  m_realQZ.compute(A, B, computeEigenvectors);
294
  m_realQZ.compute(A, B, computeEigenvectors);
309
295
  if (m_realQZ.info() == Success) {
310
  if (m_realQZ.info() == Success)
296
    // Temp space for the untransformed eigenvectors
311
  {
297
    VectorType v;
312
    m_matS = m_realQZ.matrixS();
298
    ComplexVectorType cv;
313
    const MatrixType &matT = m_realQZ.matrixT();
299
    // Resize storage
314
    if (computeEigenvectors)
315
      m_eivec = m_realQZ.matrixZ().transpose();
316
  
317
    // Compute eigenvalues from matS
318
    m_alphas.resize(A.cols());
300
    m_alphas.resize(A.cols());
319
    m_betas.resize(A.cols());
301
    m_betas.resize(A.cols());
302
    if (computeEigenvectors) {
303
      m_eivec.resize(A.cols(), A.cols());
304
      v.resize(A.cols());
305
      cv.resize(A.cols());
306
    }
307
    // Grab some references
308
    const MatrixType &mZT = m_realQZ.matrixZ().transpose();
309
    const MatrixType &mS = m_realQZ.matrixS();
310
    const MatrixType &mT = m_realQZ.matrixT();
320
    Index i = 0;
311
    Index i = 0;
321
    while (i < A.cols())
312
    while (i < A.cols()) {
322
    {
313
      if (i == A.cols() - 1 || mS.coeff(i+1, i) == Scalar(0)) {
323
      if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0))
314
        // Real eigenvalue
324
      {
315
        m_alphas.coeffRef(i) = mS.diagonal().coeff(i);
325
        m_alphas.coeffRef(i) = m_matS.coeff(i, i);
316
        m_betas.coeffRef(i)  = mT.diagonal().coeff(i);
326
        m_betas.coeffRef(i)  = matT.coeff(i,i);
317
        if (computeEigenvectors) {
318
          v.setConstant(Scalar(0.0));
319
          v.coeffRef(i) = Scalar(1.0);
320
          // For singular eigenvalues do nothing more
321
          if (!internal::isApprox(fabs(m_betas.coeffRef(i)), 0.0)) {
322
            //Non-singular eigenvalue
323
            const Scalar alpha = real(m_alphas.coeffRef(i));
324
            const Scalar beta = m_betas.coeffRef(i);
325
            for (Index j = i-1; j >= 0; j--) {
326
              const Index st = j+1;
327
              const Index sz = i-j;
328
              if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) { // 2x2 block
329
                Matrix<Scalar, 2, 1> RHS;
330
                RHS[0] = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j-1,st,1,sz) - alpha*mT.block(j-1,st,1,sz)).sum();
331
                RHS[1] = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum();
332
                Matrix<Scalar, 2, 2> LHS;
333
                LHS << beta*mS.coeffRef(j-1,j-1) - alpha*mT.coeffRef(j-1,j-1), beta*mS.coeffRef(j-1,j) - alpha*mT.coeffRef(j-1,j),
334
                      beta*mS.coeffRef(j,j-1) - alpha*mT.coeffRef(j,j-1), beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j);
335
                v.segment(j-1,2) = LHS.partialPivLu().solve(RHS);
336
                j--;
337
              } else {
338
                v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j));
339
              }
340
            }
341
          }
342
          m_eivec.col(i).real() = (mZT * v).normalized();
343
          m_eivec.col(i).imag().setConstant(0);
344
        }
327
        ++i;
345
        ++i;
328
      }
346
      } else {
329
      else
330
      {
331
        // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
347
        // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
332
        // Then taking beta=T_00*T_11, we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00):
348
        // Then taking beta=T_00*T_11, we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00):
333
349
334
        // T =  [a 0]
350
        // T =  [a 0]
335
        //      [0 b]
351
        //      [0 b]
336
        RealScalar a = matT.diagonal().coeff(i),
352
        RealScalar a = mT.diagonal().coeff(i),
337
                   b = matT.diagonal().coeff(i+1);
353
                   b = mT.diagonal().coeff(i+1);
354
        const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
355
338
        // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
356
        // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
339
        Matrix<RealScalar,2,2> S2 = m_matS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
357
        Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
340
358
341
        Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
359
        Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
342
        Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
360
        Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
343
        m_alphas.coeffRef(i)   = ComplexScalar(S2.coeff(1,1) + p,  z);
361
        const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z);
344
        m_alphas.coeffRef(i+1) = ComplexScalar(S2.coeff(1,1) + p, -z);
362
        m_alphas.coeffRef(i)   = conj(alpha);
363
        m_alphas.coeffRef(i+1) = alpha;
345
364
346
        m_betas.coeffRef(i)   =
365
        if (computeEigenvectors) {
347
        m_betas.coeffRef(i+1) = a*b;
366
          // Compute eigenvector in position (i+1) and then position (i) is just the conjugate
348
        
367
          cv.setConstant(ComplexScalar(0.0, 0.0));
368
          cv.coeffRef(i+1) = ComplexScalar(1.0, 0.0);
369
          cv.coeffRef(i) = -(beta*mS.coeffRef(i,i+1) - alpha*mT.coeffRef(i,i+1)) / (beta*mS.coeffRef(i,i) - alpha*mT.coeffRef(i,i));
370
          for (Index j = i-1; j >= 0; j--) {
371
            const Index st = j+1;
372
            const Index sz = i+1-j;
373
            if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) { // 2x2 block
374
              Matrix<ComplexScalar, 2, 1> RHS;
375
              RHS[0] = -cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j-1,st,1,sz) - alpha*mT.block(j-1,st,1,sz)).sum();
376
              RHS[1] = -cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum();
377
              Matrix<ComplexScalar, 2, 2> LHS;
378
              LHS << beta*mS.coeffRef(j-1,j-1) - alpha*mT.coeffRef(j-1,j-1), beta*mS.coeffRef(j-1,j) - alpha*mT.coeffRef(j-1,j),
379
                     beta*mS.coeffRef(j,j-1) - alpha*mT.coeffRef(j,j-1), beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j);
380
              cv.segment(j-1,2) = LHS.partialPivLu().solve(RHS);
381
              j--;
382
            } else {
383
              cv.coeffRef(j) = -cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j));
384
            }
385
          }
386
          m_eivec.col(i+1) = (mZT * cv).normalized();
387
          m_eivec.col(i)   = m_eivec.col(i+1).conjugate();
388
        }
349
        i += 2;
389
        i += 2;
350
      }
390
      }
351
    }
391
    }
392
    m_valuesOkay = true;
393
    m_vectorsOkay = computeEigenvectors;
352
  }
394
  }
353
354
  m_isInitialized = true;
355
  m_eigenvectorsOk = false;//computeEigenvectors;
356
357
  return *this;
395
  return *this;
358
}
396
}
359
397
(-)a/test/CMakeLists.txt (+1 lines)
Lines 211-216 Link Here
211
ei_add_test(eigensolver_complex)
211
ei_add_test(eigensolver_complex)
212
ei_add_test(real_qz)
212
ei_add_test(real_qz)
213
ei_add_test(eigensolver_generalized_real)
213
ei_add_test(eigensolver_generalized_real)
214
ei_add_test(eigensolver_generalized_complex)
214
ei_add_test(jacobi)
215
ei_add_test(jacobi)
215
ei_add_test(jacobisvd)
216
ei_add_test(jacobisvd)
216
ei_add_test(bdcsvd)
217
ei_add_test(bdcsvd)
(-)dfb0ec8e0d86 (+45 lines)
Added Link Here
1
// This file is part of Eigen, a lightweight C++ template library
2
// for linear algebra.
3
//
4
// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
5
//
6
// This Source Code Form is subject to the terms of the Mozilla
7
// Public License v. 2.0. If a copy of the MPL was not distributed
8
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10
#include "main.h"
11
#include <limits>
12
#include <Eigen/Eigenvalues>
13
14
template<typename MatrixType> void generalized_eigensolver_complex(const MatrixType& m)
15
{
16
  typedef typename MatrixType::Index Index;
17
  /* this test covers the following files:
18
     GeneralizedEigenSolver.h
19
  */
20
  Index rows = m.rows();
21
  Index cols = m.cols();
22
  MatrixType a = MatrixType::Random(rows,cols);
23
  MatrixType b = MatrixType::Random(rows,cols);
24
  GeneralizedEigenSolver<MatrixType> eig(a, b);
25
  typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType D = eig.eigenvalues().asDiagonal();
26
  typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType V = eig.eigenvectors();
27
  VERIFY_IS_APPROX(a*V, b*V*D);
28
}
29
30
void test_eigensolver_generalized_complex()
31
{
32
  for(int i = 0; i < g_repeat; i++) {
33
    int s = 0;
34
    CALL_SUBTEST_1( generalized_eigensolver_complex(Matrix4f()) );
35
    s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
36
    CALL_SUBTEST_2( generalized_eigensolver_complex(MatrixXd(s,s)) );
37
38
    // some trivial but implementation-wise tricky cases
39
    CALL_SUBTEST_3( generalized_eigensolver_complex(MatrixXd(1,1)) );
40
    CALL_SUBTEST_4( generalized_eigensolver_complex(MatrixXd(2,2)) );
41
    CALL_SUBTEST_5( generalized_eigensolver_complex(Matrix<double,1,1>()) );
42
    CALL_SUBTEST_6( generalized_eigensolver_complex(Matrix2d()) );
43
    TEST_SET_BUT_UNUSED_VARIABLE(s)
44
  }
45
}

Return to bug 645