* decouple the generalized selfadjoint eigenvalue problem to the standard one

* uses named values instead of bools

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"

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 {

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

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
+      * \param[in]  options Can be ComputeEigenvectors (default) or EigenvaluesOnly.
-      *    eigenvalues are computed; if false, only the eigenvalues are
-      *    computed.
       *
       * 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);
+      compute(matrix, options);
-    }
-
-    /** \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);
     }
 
     /** \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
+      * \param[in]  options Can be ComputeEigenvectors (default) or EigenvaluesOnly.
-      *    eigenvalues are computed; if false, only the eigenvalues are
-      *    computed.
       * \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
+      * function can be used to retrieve them.  If \p options equals ComputeEigenvectors,
-      * true, then the eigenvectors are also computed and can be retrieved by
+      * 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);
+    SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors);
-
-    /** \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);
 
     /** \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)
 {

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)