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 typename Base::Index Index; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default constructor for fixed-size matrices. 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default constructor is useful in cases in which the user intends to 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * perform decompositions via compute(). This constructor 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * can only be used if \p _MatrixType is a fixed-size matrix; use 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices. 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath GeneralizedSelfAdjointEigenSolver() : Base() {} 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor, pre-allocates memory for dynamic-size matrices. 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param [in] size Positive integer, size of the matrix whose 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues and eigenvectors will be computed. 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor is useful for dynamic-size matrices, when the user 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * intends to perform decompositions via compute(). The \p size 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * parameter is only used as a hint. It is not an error to give a wrong 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p size, but it may impair performance. 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute() for an example 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath GeneralizedSelfAdjointEigenSolver(Index size) 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base(size) 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil. 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matA Selfadjoint matrix in matrix pencil. 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the lower triangular part of the matrix is referenced. 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matB Positive-definite matrix in matrix pencil. 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the lower triangular part of the matrix is referenced. 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Default is #ComputeEigenvectors|#Ax_lBx. 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor calls compute(const MatrixType&, const MatrixType&, int) 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the eigenvalues and (if requested) the eigenvectors of the 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ x^* B x = 1 \f$. The eigenvectors are computed if 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \a options contains ComputeEigenvectors. 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * In addition, the two following variants can be solved via \p options: 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - \c ABx_lx: \f$ ABx = \lambda x \f$ 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - \c BAx_lx: \f$ BAx = \lambda x \f$ 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute(const MatrixType&, const MatrixType&, int) 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int options = ComputeEigenvectors|Ax_lBx) 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : Base(matA.cols()) 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matA, matB, options); 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Computes generalized eigendecomposition of given matrix pencil. 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matA Selfadjoint matrix in matrix pencil. 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the lower triangular part of the matrix is referenced. 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matB Positive-definite matrix in matrix pencil. 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the lower triangular part of the matrix is referenced. 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Default is #ComputeEigenvectors|#Ax_lBx. 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Reference to \c *this 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Accoring to \p options, this function computes eigenvalues and (if requested) 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the eigenvectors of one of the following three generalized eigenproblems: 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - \c Ax_lBx: \f$ Ax = \lambda B x \f$ 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - \c ABx_lx: \f$ ABx = \lambda x \f$ 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - \c BAx_lx: \f$ BAx = \lambda x \f$ 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix \f$ B \f$. 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$. 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The eigenvalues() function can be used to retrieve 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the eigenvalues. If \p options contains ComputeEigenvectors, then the 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvectors are also computed and can be retrieved by calling 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvectors(). 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The implementation uses LLT to compute the Cholesky decomposition 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ B = LL^* \f$ and computes the classical eigendecomposition 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and of \f$ L^{*} A L \f$ otherwise. This solves the 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * generalized eigenproblem, because any solution of the generalized 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * can be made for the two other variants. 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int options = ComputeEigenvectors|Ax_lBx); 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathGeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>:: 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathcompute(const MatrixType& matA, const MatrixType& matB, int options) 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert((options&~(EigVecMask|GenEigMask))==0 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && (options&EigVecMask)!=EigVecMask 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx) 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "invalid option parameter"); 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors); 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the cholesky decomposition of matB = L L' = U'U 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT<MatrixType> cholB(matB); 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int type = (options&GenEigMask); 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(type==0) 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath type = Ax_lBx; 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(type==Ax_lBx) 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // compute C = inv(L) A inv(L') 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matC = matA.template selfadjointView<Lower>(); 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cholB.matrixL().template solveInPlace<OnTheLeft>(matC); 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cholB.matrixU().template solveInPlace<OnTheRight>(matC); 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly ); 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // transform back the eigen vectors: evecs = inv(U) * evecs 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeEigVecs) 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cholB.matrixU().solveInPlace(Base::m_eivec); 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(type==ABx_lx) 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // compute C = L' A L 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matC = matA.template selfadjointView<Lower>(); 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matC = matC * cholB.matrixL(); 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matC = cholB.matrixU() * matC; 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // transform back the eigen vectors: evecs = inv(U) * evecs 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeEigVecs) 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cholB.matrixU().solveInPlace(Base::m_eivec); 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(type==BAx_lx) 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // compute C = L' A L 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matC = matA.template selfadjointView<Lower>(); 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matC = matC * cholB.matrixL(); 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matC = cholB.matrixU() * matC; 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // transform back the eigen vectors: evecs = L * evecs 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeEigVecs) 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Base::m_eivec = cholB.matrixL() * Base::m_eivec; 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 228