GeneralizedSelfAdjointEigenSolver.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 #ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
27 #define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
28 
29 #include "./Tridiagonalization.h"
30 
31 namespace Eigen {
32 
62 template<typename _MatrixType>
64 {
66  public:
67 
68  typedef typename Base::Index Index;
69  typedef _MatrixType MatrixType;
70 
79 
93  : Base(size)
94  {}
95 
123  int options = ComputeEigenvectors|Ax_lBx)
124  : Base(matA.cols())
125  {
126  compute(matA, matB, options);
127  }
128 
170  int options = ComputeEigenvectors|Ax_lBx);
171 
172  protected:
173 
174 };
175 
176 
177 template<typename MatrixType>
178 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
179 compute(const MatrixType& matA, const MatrixType& matB, int options)
180 {
181  eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
182  eigen_assert((options&~(EigVecMask|GenEigMask))==0
183  && (options&EigVecMask)!=EigVecMask
184  && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
185  || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
186  && "invalid option parameter");
187 
188  bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
189 
190  // Compute the cholesky decomposition of matB = L L' = U'U
191  LLT<MatrixType> cholB(matB);
192 
193  int type = (options&GenEigMask);
194  if(type==0)
195  type = Ax_lBx;
196 
197  if(type==Ax_lBx)
198  {
199  // compute C = inv(L) A inv(L')
200  MatrixType matC = matA.template selfadjointView<Lower>();
201  cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
202  cholB.matrixU().template solveInPlace<OnTheRight>(matC);
203 
204  Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
205 
206  // transform back the eigen vectors: evecs = inv(U) * evecs
207  if(computeEigVecs)
208  cholB.matrixU().solveInPlace(Base::m_eivec);
209  }
210  else if(type==ABx_lx)
211  {
212  // compute C = L' A L
213  MatrixType matC = matA.template selfadjointView<Lower>();
214  matC = matC * cholB.matrixL();
215  matC = cholB.matrixU() * matC;
216 
217  Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
218 
219  // transform back the eigen vectors: evecs = inv(U) * evecs
220  if(computeEigVecs)
221  cholB.matrixU().solveInPlace(Base::m_eivec);
222  }
223  else if(type==BAx_lx)
224  {
225  // compute C = L' A L
226  MatrixType matC = matA.template selfadjointView<Lower>();
227  matC = matC * cholB.matrixL();
228  matC = cholB.matrixU() * matC;
229 
230  Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
231 
232  // transform back the eigen vectors: evecs = L * evecs
233  if(computeEigVecs)
234  Base::m_eivec = cholB.matrixL() * Base::m_eivec;
235  }
236 
237  return *this;
238 }
239 
240 } // end namespace Eigen
241 
242 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H