diff options
| -rw-r--r-- | Eigen/Eigenvalues | 1 | ||||
| -rw-r--r-- | Eigen/src/Core/util/Constants.h | 14 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h | 210 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 129 | ||||
| -rw-r--r-- | test/eigensolver_selfadjoint.cpp | 2 |
5 files changed, 244 insertions, 112 deletions
diff --git a/Eigen/Eigenvalues b/Eigen/Eigenvalues index 5a9757ad5..381802cb5 100644 --- a/Eigen/Eigenvalues +++ b/Eigen/Eigenvalues @@ -30,6 +30,7 @@ namespace Eigen { #include "src/Eigenvalues/RealSchur.h" #include "src/Eigenvalues/EigenSolver.h" #include "src/Eigenvalues/SelfAdjointEigenSolver.h" +#include "src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h" #include "src/Eigenvalues/HessenbergDecomposition.h" #include "src/Eigenvalues/ComplexSchur.h" #include "src/Eigenvalues/ComplexEigenSolver.h" diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 831b82890..0553c40f4 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -234,6 +234,20 @@ enum { IsSparse }; +enum DecompositionOptions { + Pivoting = 0x01, // LDLT, + NoPivoting = 0x02, // LDLT, + ComputeU = 0x10, // SVD, + ComputeR = 0x20, // SVD, + EigenvaluesOnly = 0x40, // all eigen solvers + ComputeEigenvectors = 0x80, // all eigen solvers + EigVecMask = EigenvaluesOnly | ComputeEigenvectors, + Ax_lBx = 0x100, + ABx_lx = 0x200, + BAx_lx = 0x400, + GenEigMask = Ax_lBx | ABx_lx | BAx_lx +}; + /** \brief Enum for reporting the status of a computation. */ enum ComputationInfo { diff --git a/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h new file mode 100644 index 000000000..66d32dee6 --- /dev/null +++ b/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h @@ -0,0 +1,210 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2010 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H +#define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H + +#include "./EigenvaluesCommon.h" +#include "./Tridiagonalization.h" + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * \nonstableyet + * + * \class GeneralizedSelfAdjointEigenSolver + * + * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem + * + * \tparam _MatrixType the type of the matrix of which we are computing the + * eigendecomposition; this is expected to be an instantiation of the Matrix + * class template. + * + * This class solves the generalized eigenvalue problem + * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be + * selfadjoint and the matrix \f$ B \f$ should be positive definite. + * + * Only the \b lower \b triangular \b part of the input matrix is referenced. + * + * Call the function compute() to compute the eigenvalues and eigenvectors of + * a given matrix. Alternatively, you can use the + * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) + * constructor which computes the eigenvalues and eigenvectors at construction time. + * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues() + * and eigenvectors() functions. + * + * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) + * contains an example of the typical use of this class. + * + * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver + */ +template<typename _MatrixType> +class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType> +{ + typedef SelfAdjointEigenSolver<_MatrixType> Base; + public: + + typedef typename Base::Index Index; + typedef _MatrixType MatrixType; + + /** \brief Default constructor for fixed-size matrices. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via compute(const MatrixType&, bool) or + * compute(const MatrixType&, const MatrixType&, bool). This constructor + * can only be used if \p _MatrixType is a fixed-size matrix; use + * SelfAdjointEigenSolver(Index) for dynamic-size matrices. + * + * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp + * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out + */ + GeneralizedSelfAdjointEigenSolver() : Base() {} + + /** \brief Constructor, pre-allocates memory for dynamic-size matrices. + * + * \param [in] size Positive integer, size of the matrix whose + * eigenvalues and eigenvectors will be computed. + * + * This constructor is useful for dynamic-size matrices, when the user + * intends to perform decompositions via compute(const MatrixType&, bool) + * or compute(const MatrixType&, const MatrixType&, bool). The \p size + * parameter is only used as a hint. It is not an error to give a wrong + * \p size, but it may impair performance. + * + * \sa compute(const MatrixType&, bool) for an example + */ + GeneralizedSelfAdjointEigenSolver(Index size) + : Base(size) + {} + + /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil. + * + * \param[in] matA Selfadjoint matrix in matrix pencil. + * Only the lower triangular part of the matrix is referenced. + * \param[in] matB Positive-definite matrix in matrix pencil. + * Only the lower triangular part of the matrix is referenced. + * \param[in] options A or-ed set of flags {ComputeEigenvectors,EigenvaluesOnly} | {Ax_lBx,ABx_lx,BAx_lx}. + * Default is ComputeEigenvectors|Ax_lBx. + * + * This constructor calls compute(const MatrixType&, const MatrixType&, int) + * to compute the eigenvalues and (if requested) the eigenvectors of the + * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the + * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix + * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property + * \f$ x^* B x = 1 \f$. The eigenvectors are computed if + * \a options contains ComputeEigenvectors. + * + * In addition, the two following variants can be solved via \p options: + * - \c ABx_lx: \f$ ABx = \lambda x \f$ + * - \c BAx_lx: \f$ BAx = \lambda x \f$ + * + * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp + * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out + * + * \sa compute(const MatrixType&, const MatrixType&, int) + */ + GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, + int options = ComputeEigenvectors|Ax_lBx) + : Base(matA.cols()) + { + compute(matA, matB, options); + } + + /** \brief Computes generalized eigendecomposition of given matrix pencil. + * + * \param[in] matA Selfadjoint matrix in matrix pencil. + * Only the lower triangular part of the matrix is referenced. + * \param[in] matB Positive-definite matrix in matrix pencil. + * Only the lower triangular part of the matrix is referenced. + * \param[in] options A or-ed set of flags {ComputeEigenvectors,EigenvaluesOnly} | {Ax_lBx,ABx_lx,BAx_lx}. + * Default is ComputeEigenvectors|Ax_lBx. + * + * \returns Reference to \c *this + * + * If \p options contains Ax_lBx (the default), this function computes eigenvalues + * and (if requested) the eigenvectors of the generalized eigenproblem + * \f$ Ax = \lambda B x \f$ with \a matA the selfadjoint + * matrix \f$ A \f$ and \a matB the positive definite + * matrix \f$ B \f$. In addition, each eigenvector \f$ x \f$ + * satisfies the property \f$ x^* B x = 1 \f$. + * + * In addition, the two following variants can be solved via \p options: + * - \c ABx_lx: \f$ ABx = \lambda x \f$ + * - \c BAx_lx: \f$ BAx = \lambda x \f$ + * + * The eigenvalues() function can be used to retrieve + * the eigenvalues. If \p options contains ComputeEigenvectors, then the + * eigenvectors are also computed and can be retrieved by calling + * eigenvectors(). + * + * The implementation uses LLT to compute the Cholesky decomposition + * \f$ B = LL^* \f$ and calls compute(const MatrixType&, bool) to compute + * the eigendecomposition \f$ L^{-1} A (L^*)^{-1} \f$. This solves the + * generalized eigenproblem, because any solution of the generalized + * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution + * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the + * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. + * + * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp + * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out + * + * \sa SelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) + */ + GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, + int options = ComputeEigenvectors|Ax_lBx); + + protected: + +}; + + +template<typename MatrixType> +GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>:: +compute(const MatrixType& matA, const MatrixType& matB, int options) +{ + ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); + ei_assert((options&~(EigVecMask|GenEigMask))==0 + && (options&EigVecMask)!=EigVecMask + && ((options&GenEigMask)==Ax_lBx || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx) + && "invalid option parameter"); + + bool computeEigVecs = (options&EigVecMask)==ComputeEigenvectors; + + // Compute the cholesky decomposition of matB = L L' = U'U + LLT<MatrixType> cholB(matB); + + // compute C = inv(L) A inv(L') + MatrixType matC = matA.template selfadjointView<Lower>(); + cholB.matrixL().template solveInPlace<OnTheLeft>(matC); + cholB.matrixU().template solveInPlace<OnTheRight>(matC); + + Base::compute(matC, options&EigVecMask); + + // transform back the eigen vectors: evecs = inv(U) * evecs + if(computeEigVecs) + cholB.matrixU().solveInPlace(Base::m_eivec); + + return *this; +} + +#endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index a77b96186..23ad85403 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -57,10 +57,6 @@ * * Only the \b lower \b triangular \b part of the input matrix is referenced. * - * This class can also be used to solve the generalized eigenvalue problem - * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be - * selfadjoint and the matrix \f$ B \f$ should be positive definite. - * * Call the function compute() to compute the eigenvalues and eigenvectors of * a given matrix. Alternatively, you can use the * SelfAdjointEigenSolver(const MatrixType&, bool) constructor which computes @@ -70,6 +66,9 @@ * * The documentation for SelfAdjointEigenSolver(const MatrixType&, bool) * contains an example of the typical use of this class. + * + * To solve the \em generalized eigenvalue problem \f$ Av = \lambda Bv \f$ and + * the like see the class GeneralizedSelfAdjointEigenSolver. * * \sa MatrixBase::eigenvalues(), class EigenSolver, class ComplexEigenSolver */ @@ -147,13 +146,11 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \param[in] matrix Selfadjoint matrix whose eigendecomposition is to * be computed. Only the lower triangular part of the matrix is referenced. - * \param[in] computeEigenvectors If true, both the eigenvectors and the - * eigenvalues are computed; if false, only the eigenvalues are - * computed. + * \param[in] options Can be ComputeEigenvectors (default) or EigenvaluesOnly. * * This constructor calls compute(const MatrixType&, bool) to compute the * eigenvalues of the matrix \p matrix. The eigenvectors are computed if - * \p computeEigenvectors is true. + * \p options equals ComputeEigenvectors. * * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.cpp * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.out @@ -161,60 +158,25 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * \sa compute(const MatrixType&, bool), * SelfAdjointEigenSolver(const MatrixType&, const MatrixType&, bool) */ - SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) + SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors) : m_eivec(matrix.rows(), matrix.cols()), m_eivalues(matrix.cols()), m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1), m_isInitialized(false) { - compute(matrix, computeEigenvectors); - } - - /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil. - * - * \param[in] matA Selfadjoint matrix in matrix pencil. - * Only the lower triangular part of the matrix is referenced. - * \param[in] matB Positive-definite matrix in matrix pencil. - * Only the lower triangular part of the matrix is referenced. - * \param[in] computeEigenvectors If true, both the eigenvectors and the - * eigenvalues are computed; if false, only the eigenvalues are - * computed. - * - * This constructor calls compute(const MatrixType&, const MatrixType&, bool) - * to compute the eigenvalues and (if requested) the eigenvectors of the - * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the - * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix - * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property - * \f$ x^* B x = 1 \f$. The eigenvectors are computed if - * \a computeEigenvectors is true. - * - * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp - * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out - * - * \sa compute(const MatrixType&, const MatrixType&, bool), - * SelfAdjointEigenSolver(const MatrixType&, bool) - */ - SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) - : m_eivec(matA.rows(), matA.cols()), - m_eivalues(matA.cols()), - m_subdiag(matA.rows() > 1 ? matA.rows() - 1 : 1), - m_isInitialized(false) - { - compute(matA, matB, computeEigenvectors); + compute(matrix, options); } /** \brief Computes eigendecomposition of given matrix. * * \param[in] matrix Selfadjoint matrix whose eigendecomposition is to * be computed. Only the lower triangular part of the matrix is referenced. - * \param[in] computeEigenvectors If true, both the eigenvectors and the - * eigenvalues are computed; if false, only the eigenvalues are - * computed. + * \param[in] options Can be ComputeEigenvectors (default) or EigenvaluesOnly. * \returns Reference to \c *this * * This function computes the eigenvalues of \p matrix. The eigenvalues() - * function can be used to retrieve them. If \p computeEigenvectors is - * true, then the eigenvectors are also computed and can be retrieved by + * function can be used to retrieve them. If \p options equals ComputeEigenvectors, + * then the eigenvectors are also computed and can be retrieved by * calling eigenvectors(). * * This implementation uses a symmetric QR algorithm. The matrix is first @@ -235,44 +197,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa SelfAdjointEigenSolver(const MatrixType&, bool) */ - SelfAdjointEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); - - /** \brief Computes generalized eigendecomposition of given matrix pencil. - * - * \param[in] matA Selfadjoint matrix in matrix pencil. - * Only the lower triangular part of the matrix is referenced. - * \param[in] matB Positive-definite matrix in matrix pencil. - * Only the lower triangular part of the matrix is referenced. - * \param[in] computeEigenvectors If true, both the eigenvectors and the - * eigenvalues are computed; if false, only the eigenvalues are - * computed. - * \returns Reference to \c *this - * - * This function computes eigenvalues and (if requested) the eigenvectors - * of the generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA - * the selfadjoint matrix \f$ A \f$ and \a matB the positive definite - * matrix \f$ B \f$. In addition, each eigenvector \f$ x \f$ - * satisfies the property \f$ x^* B x = 1 \f$. - * - * The eigenvalues() function can be used to retrieve - * the eigenvalues. If \p computeEigenvectors is true, then the - * eigenvectors are also computed and can be retrieved by calling - * eigenvectors(). - * - * The implementation uses LLT to compute the Cholesky decomposition - * \f$ B = LL^* \f$ and calls compute(const MatrixType&, bool) to compute - * the eigendecomposition \f$ L^{-1} A (L^*)^{-1} \f$. This solves the - * generalized eigenproblem, because any solution of the generalized - * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution - * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the - * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. - * - * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp - * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out - * - * \sa SelfAdjointEigenSolver(const MatrixType&, const MatrixType&, bool) - */ - SelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true); + SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors); /** \brief Returns the eigenvectors of given matrix (pencil). * @@ -414,9 +339,14 @@ template<typename RealScalar, typename Scalar, typename Index> static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n); template<typename MatrixType> -SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) +SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> +::compute(const MatrixType& matrix, int options) { - assert(matrix.cols() == matrix.rows()); + ei_assert(matrix.cols() == matrix.rows()); + ei_assert((options&~(EigVecMask|GenEigMask))==0 + && (options&EigVecMask)!=EigVecMask + && "invalid option parameter"); + bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; Index n = matrix.cols(); m_eivalues.resize(n,1); @@ -497,29 +427,6 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute( return *this; } -template<typename MatrixType> -SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>:: -compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors) -{ - ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); - - // Compute the cholesky decomposition of matB = L L' = U'U - LLT<MatrixType> cholB(matB); - - // compute C = inv(L) A inv(L') - MatrixType matC = matA.template selfadjointView<Lower>(); - cholB.matrixL().template solveInPlace<OnTheLeft>(matC); - cholB.matrixU().template solveInPlace<OnTheRight>(matC); - - compute(matC, computeEigenvectors); - - // transform back the eigen vectors: evecs = inv(U) * evecs - if(computeEigenvectors) - cholB.matrixU().solveInPlace(m_eivec); - - return *this; -} - template<typename RealScalar, typename Scalar, typename Index> static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n) { diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 2c2d5c985..d3bbb3a7e 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -59,7 +59,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) SelfAdjointEigenSolver<MatrixType> eiSymm(symmA); // generalized eigen pb - SelfAdjointEigenSolver<MatrixType> eiSymmGen(symmA, symmB); + GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmA, symmB); #ifdef HAS_GSL if (ei_is_same_type<RealScalar,double>::ret) |