This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
View | Details | Raw Unified | Return to bug 1739
Collapse All | Expand All

(-)a/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h (-19 / +21 lines)
Lines 59-65 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT Link Here
59
      * can only be used if \p _MatrixType is a fixed-size matrix; use
59
      * can only be used if \p _MatrixType is a fixed-size matrix; use
60
      * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices.
60
      * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices.
61
      */
61
      */
62
    GeneralizedSelfAdjointEigenSolver() : Base() {}
62
    GeneralizedSelfAdjointEigenSolver() : Base(), m_cholB(), m_matC() {}
63
63
64
    /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
64
    /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
65
      *
65
      *
Lines 74-80 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT Link Here
74
      * \sa compute() for an example
74
      * \sa compute() for an example
75
      */
75
      */
76
    explicit GeneralizedSelfAdjointEigenSolver(Index size)
76
    explicit GeneralizedSelfAdjointEigenSolver(Index size)
77
        : Base(size)
77
        : Base(size), m_cholB(size), m_matC(size, size)
78
    {}
78
    {}
79
79
80
    /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil.
80
    /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil.
Lines 105-111 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT Link Here
105
      */
105
      */
106
    GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
106
    GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
107
                                      int options = ComputeEigenvectors|Ax_lBx)
107
                                      int options = ComputeEigenvectors|Ax_lBx)
108
      : Base(matA.cols())
108
      : Base(matA.cols()), m_cholB(matA.cols()), m_matC(matA.cols(), matA.cols())
109
    {
109
    {
110
      compute(matA, matB, options);
110
      compute(matA, matB, options);
111
    }
111
    }
Lines 155-160 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT Link Here
155
155
156
  protected:
156
  protected:
157
157
158
    LLT<MatrixType> m_cholB;
159
    MatrixType m_matC;
158
};
160
};
159
161
160
162
Lines 172-178 compute(const MatrixType& matA, const MatrixType& matB, int options) Link Here
172
  bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
174
  bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
173
175
174
  // Compute the cholesky decomposition of matB = L L' = U'U
176
  // Compute the cholesky decomposition of matB = L L' = U'U
175
  LLT<MatrixType> cholB(matB);
177
  m_cholB.compute(matB);
176
178
177
  int type = (options&GenEigMask);
179
  int type = (options&GenEigMask);
178
  if(type==0)
180
  if(type==0)
Lines 181-221 compute(const MatrixType& matA, const MatrixType& matB, int options) Link Here
181
  if(type==Ax_lBx)
183
  if(type==Ax_lBx)
182
  {
184
  {
183
    // compute C = inv(L) A inv(L')
185
    // compute C = inv(L) A inv(L')
184
    MatrixType matC = matA.template selfadjointView<Lower>();
186
    m_matC = matA.template selfadjointView<Lower>();
185
    cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
187
    m_cholB.matrixL().template solveInPlace<OnTheLeft>(m_matC);
186
    cholB.matrixU().template solveInPlace<OnTheRight>(matC);
188
    m_cholB.matrixU().template solveInPlace<OnTheRight>(m_matC);
187
189
188
    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
190
    Base::compute(m_matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
189
191
190
    // transform back the eigen vectors: evecs = inv(U) * evecs
192
    // transform back the eigen vectors: evecs = inv(U) * evecs
191
    if(computeEigVecs)
193
    if(computeEigVecs)
192
      cholB.matrixU().solveInPlace(Base::m_eivec);
194
      m_cholB.matrixU().solveInPlace(Base::m_eivec);
193
  }
195
  }
194
  else if(type==ABx_lx)
196
  else if(type==ABx_lx)
195
  {
197
  {
196
    // compute C = L' A L
198
    // compute C = L' A L
197
    MatrixType matC = matA.template selfadjointView<Lower>();
199
    m_matC = matA.template selfadjointView<Lower>();
198
    matC = matC * cholB.matrixL();
200
    m_matC = m_matC * m_cholB.matrixL();
199
    matC = cholB.matrixU() * matC;
201
    m_matC = m_cholB.matrixU() * m_matC;
200
202
201
    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
203
    Base::compute(m_matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
202
204
203
    // transform back the eigen vectors: evecs = inv(U) * evecs
205
    // transform back the eigen vectors: evecs = inv(U) * evecs
204
    if(computeEigVecs)
206
    if(computeEigVecs)
205
      cholB.matrixU().solveInPlace(Base::m_eivec);
207
      m_cholB.matrixU().solveInPlace(Base::m_eivec);
206
  }
208
  }
207
  else if(type==BAx_lx)
209
  else if(type==BAx_lx)
208
  {
210
  {
209
    // compute C = L' A L
211
    // compute C = L' A L
210
    MatrixType matC = matA.template selfadjointView<Lower>();
212
    m_matC = matA.template selfadjointView<Lower>();
211
    matC = matC * cholB.matrixL();
213
    m_matC = m_matC * m_cholB.matrixL();
212
    matC = cholB.matrixU() * matC;
214
    m_matC = m_cholB.matrixU() * m_matC;
213
215
214
    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
216
    Base::compute(m_matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
215
217
216
    // transform back the eigen vectors: evecs = L * evecs
218
    // transform back the eigen vectors: evecs = L * evecs
217
    if(computeEigVecs)
219
    if(computeEigVecs)
218
      Base::m_eivec = cholB.matrixL() * Base::m_eivec;
220
      Base::m_eivec = m_cholB.matrixL() * Base::m_eivec;
219
  }
221
  }
220
222
221
  return *this;
223
  return *this;

Return to bug 1739