Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
GeneralizedSelfAdjointEigenSolver.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
12 #define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
13 
14 #include "./Tridiagonalization.h"
15 
16 #include "./InternalHeaderCheck.h"
17 
18 namespace Eigen {
19 
49 template<typename MatrixType_>
51 {
53  public:
54 
55  typedef MatrixType_ MatrixType;
56 
65 
79  : Base(size)
80  {}
81 
108  GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
109  int options = ComputeEigenvectors|Ax_lBx)
110  : Base(matA.cols())
111  {
112  compute(matA, matB, options);
113  }
114 
155  GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
156  int options = ComputeEigenvectors|Ax_lBx);
157 
158  protected:
159 
160 };
161 
162 
163 template<typename MatrixType>
165 compute(const MatrixType& matA, const MatrixType& matB, int options)
166 {
167  eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
168  eigen_assert((options&~(EigVecMask|GenEigMask))==0
169  && (options&EigVecMask)!=EigVecMask
170  && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
171  || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
172  && "invalid option parameter");
173 
174  bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
175 
176  // Compute the cholesky decomposition of matB = L L' = U'U
177  LLT<MatrixType> cholB(matB);
178 
179  int type = (options&GenEigMask);
180  if(type==0)
181  type = Ax_lBx;
182 
183  if(type==Ax_lBx)
184  {
185  // compute C = inv(L) A inv(L')
186  MatrixType matC = matA.template selfadjointView<Lower>();
187  cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
188  cholB.matrixU().template solveInPlace<OnTheRight>(matC);
189 
190  Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
191 
192  // transform back the eigen vectors: evecs = inv(U) * evecs
193  if(computeEigVecs)
194  cholB.matrixU().solveInPlace(Base::m_eivec);
195  }
196  else if(type==ABx_lx)
197  {
198  // compute C = L' A L
199  MatrixType matC = matA.template selfadjointView<Lower>();
200  matC = matC * cholB.matrixL();
201  matC = cholB.matrixU() * matC;
202 
203  Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
204 
205  // transform back the eigen vectors: evecs = inv(U) * evecs
206  if(computeEigVecs)
207  cholB.matrixU().solveInPlace(Base::m_eivec);
208  }
209  else if(type==BAx_lx)
210  {
211  // compute C = L' A L
212  MatrixType matC = matA.template selfadjointView<Lower>();
213  matC = matC * cholB.matrixL();
214  matC = cholB.matrixU() * matC;
215 
216  Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
217 
218  // transform back the eigen vectors: evecs = L * evecs
219  if(computeEigVecs)
220  Base::m_eivec = cholB.matrixL() * Base::m_eivec;
221  }
222 
223  return *this;
224 }
225 
226 } // end namespace Eigen
227 
228 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem.
Definition: GeneralizedSelfAdjointEigenSolver.h:51
GeneralizedSelfAdjointEigenSolver()
Default constructor for fixed-size matrices.
Definition: GeneralizedSelfAdjointEigenSolver.h:64
GeneralizedSelfAdjointEigenSolver(Index size)
Constructor, pre-allocates memory for dynamic-size matrices.
Definition: GeneralizedSelfAdjointEigenSolver.h:78
GeneralizedSelfAdjointEigenSolver & compute(const MatrixType &matA, const MatrixType &matB, int options=ComputeEigenvectors|Ax_lBx)
Computes generalized eigendecomposition of given matrix pencil.
Definition: GeneralizedSelfAdjointEigenSolver.h:165
GeneralizedSelfAdjointEigenSolver(const MatrixType &matA, const MatrixType &matB, int options=ComputeEigenvectors|Ax_lBx)
Constructor; computes generalized eigendecomposition of given matrix pencil.
Definition: GeneralizedSelfAdjointEigenSolver.h:108
Standard Cholesky decomposition (LL^T) of a matrix and associated features.
Definition: LLT.h:70
Traits::MatrixU matrixU() const
Definition: LLT.h:130
Traits::MatrixL matrixL() const
Definition: LLT.h:137
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Definition: SelfAdjointEigenSolver.h:79
Eigen::Index Index
Definition: SelfAdjointEigenSolver.h:92
@ Ax_lBx
Definition: Constants.h:412
@ ComputeEigenvectors
Definition: Constants.h:407
@ BAx_lx
Definition: Constants.h:418
@ ABx_lx
Definition: Constants.h:415
@ EigenvaluesOnly
Definition: Constants.h:404
Namespace containing all symbols from the Eigen library.
Definition: Core:139