Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
00027 #define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
00028
00029 #include "./Tridiagonalization.h"
00030
00031 namespace Eigen {
00032
00062 template<typename _MatrixType>
00063 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
00064 {
00065 typedef SelfAdjointEigenSolver<_MatrixType> Base;
00066 public:
00067
00068 typedef typename Base::Index Index;
00069 typedef _MatrixType MatrixType;
00070
00078 GeneralizedSelfAdjointEigenSolver() : Base() {}
00079
00092 GeneralizedSelfAdjointEigenSolver(Index size)
00093 : Base(size)
00094 {}
00095
00122 GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
00123 int options = ComputeEigenvectors|Ax_lBx)
00124 : Base(matA.cols())
00125 {
00126 compute(matA, matB, options);
00127 }
00128
00169 GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
00170 int options = ComputeEigenvectors|Ax_lBx);
00171
00172 protected:
00173
00174 };
00175
00176
00177 template<typename MatrixType>
00178 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
00179 compute(const MatrixType& matA, const MatrixType& matB, int options)
00180 {
00181 eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
00182 eigen_assert((options&~(EigVecMask|GenEigMask))==0
00183 && (options&EigVecMask)!=EigVecMask
00184 && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
00185 || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
00186 && "invalid option parameter");
00187
00188 bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
00189
00190
00191 LLT<MatrixType> cholB(matB);
00192
00193 int type = (options&GenEigMask);
00194 if(type==0)
00195 type = Ax_lBx;
00196
00197 if(type==Ax_lBx)
00198 {
00199
00200 MatrixType matC = matA.template selfadjointView<Lower>();
00201 cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
00202 cholB.matrixU().template solveInPlace<OnTheRight>(matC);
00203
00204 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
00205
00206
00207 if(computeEigVecs)
00208 cholB.matrixU().solveInPlace(Base::m_eivec);
00209 }
00210 else if(type==ABx_lx)
00211 {
00212
00213 MatrixType matC = matA.template selfadjointView<Lower>();
00214 matC = matC * cholB.matrixL();
00215 matC = cholB.matrixU() * matC;
00216
00217 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
00218
00219
00220 if(computeEigVecs)
00221 cholB.matrixU().solveInPlace(Base::m_eivec);
00222 }
00223 else if(type==BAx_lx)
00224 {
00225
00226 MatrixType matC = matA.template selfadjointView<Lower>();
00227 matC = matC * cholB.matrixL();
00228 matC = cholB.matrixU() * matC;
00229
00230 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
00231
00232
00233 if(computeEigVecs)
00234 Base::m_eivec = cholB.matrixL() * Base::m_eivec;
00235 }
00236
00237 return *this;
00238 }
00239
00240 }
00241
00242 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H