1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "./Tridiagonalization.h" 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \eigenvalues_module \ingroup Eigenvalues_Module 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class GeneralizedSelfAdjointEigenSolver 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the matrix of which we are computing the 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigendecomposition; this is expected to be an instantiation of the Matrix 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * class template. 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class solves the generalized eigenvalue problem 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * selfadjoint and the matrix \f$ B \f$ should be positive definite. 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the \b lower \b triangular \b part of the input matrix is referenced. 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Call the function compute() to compute the eigenvalues and eigenvectors of 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * a given matrix. Alternatively, you can use the 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * constructor which computes the eigenvalues and eigenvectors at construction time. 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues() 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and eigenvectors() functions. 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * contains an example of the typical use of this class. 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType> 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType> 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef SelfAdjointEigenSolver<_MatrixType> Base; 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default constructor for fixed-size matrices. 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default constructor is useful in cases in which the user intends to 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * perform decompositions via compute(). This constructor 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * can only be used if \p _MatrixType is a fixed-size matrix; use 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices. 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath GeneralizedSelfAdjointEigenSolver() : Base() {} 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor, pre-allocates memory for dynamic-size matrices. 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param [in] size Positive integer, size of the matrix whose 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues and eigenvectors will be computed. 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor is useful for dynamic-size matrices, when the user 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * intends to perform decompositions via compute(). The \p size 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * parameter is only used as a hint. It is not an error to give a wrong 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p size, but it may impair performance. 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute() for an example 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit GeneralizedSelfAdjointEigenSolver(Index size) 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base(size) 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil. 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matA Selfadjoint matrix in matrix pencil. 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the lower triangular part of the matrix is referenced. 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matB Positive-definite matrix in matrix pencil. 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the lower triangular part of the matrix is referenced. 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Default is #ComputeEigenvectors|#Ax_lBx. 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor calls compute(const MatrixType&, const MatrixType&, int) 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the eigenvalues and (if requested) the eigenvectors of the 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ x^* B x = 1 \f$. The eigenvectors are computed if 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \a options contains ComputeEigenvectors. 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * In addition, the two following variants can be solved via \p options: 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - \c ABx_lx: \f$ ABx = \lambda x \f$ 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - \c BAx_lx: \f$ BAx = \lambda x \f$ 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute(const MatrixType&, const MatrixType&, int) 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int options = ComputeEigenvectors|Ax_lBx) 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base(matA.cols()) 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matA, matB, options); 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Computes generalized eigendecomposition of given matrix pencil. 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matA Selfadjoint matrix in matrix pencil. 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the lower triangular part of the matrix is referenced. 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matB Positive-definite matrix in matrix pencil. 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the lower triangular part of the matrix is referenced. 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Default is #ComputeEigenvectors|#Ax_lBx. 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Reference to \c *this 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Accoring to \p options, this function computes eigenvalues and (if requested) 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the eigenvectors of one of the following three generalized eigenproblems: 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - \c Ax_lBx: \f$ Ax = \lambda B x \f$ 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - \c ABx_lx: \f$ ABx = \lambda x \f$ 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - \c BAx_lx: \f$ BAx = \lambda x \f$ 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix \f$ B \f$. 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$. 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The eigenvalues() function can be used to retrieve 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the eigenvalues. If \p options contains ComputeEigenvectors, then the 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvectors are also computed and can be retrieved by calling 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvectors(). 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The implementation uses LLT to compute the Cholesky decomposition 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ B = LL^* \f$ and computes the classical eigendecomposition 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and of \f$ L^{*} A L \f$ otherwise. This solves the 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * generalized eigenproblem, because any solution of the generalized 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * can be made for the two other variants. 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int options = ComputeEigenvectors|Ax_lBx); 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathGeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>:: 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathcompute(const MatrixType& matA, const MatrixType& matB, int options) 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert((options&~(EigVecMask|GenEigMask))==0 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && (options&EigVecMask)!=EigVecMask 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx) 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "invalid option parameter"); 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors); 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the cholesky decomposition of matB = L L' = U'U 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT<MatrixType> cholB(matB); 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int type = (options&GenEigMask); 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(type==0) 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath type = Ax_lBx; 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(type==Ax_lBx) 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // compute C = inv(L) A inv(L') 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matC = matA.template selfadjointView<Lower>(); 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cholB.matrixL().template solveInPlace<OnTheLeft>(matC); 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cholB.matrixU().template solveInPlace<OnTheRight>(matC); 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly ); 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // transform back the eigen vectors: evecs = inv(U) * evecs 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeEigVecs) 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cholB.matrixU().solveInPlace(Base::m_eivec); 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(type==ABx_lx) 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // compute C = L' A L 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matC = matA.template selfadjointView<Lower>(); 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matC = matC * cholB.matrixL(); 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matC = cholB.matrixU() * matC; 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // transform back the eigen vectors: evecs = inv(U) * evecs 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeEigVecs) 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cholB.matrixU().solveInPlace(Base::m_eivec); 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(type==BAx_lx) 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // compute C = L' A L 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matC = matA.template selfadjointView<Lower>(); 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matC = matC * cholB.matrixL(); 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matC = cholB.matrixU() * matC; 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // transform back the eigen vectors: evecs = L * evecs 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeEigVecs) 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::m_eivec = cholB.matrixL() * Base::m_eivec; 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 227