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; |