* Added a generalized eigen solver for the selfadjoint case.

(as new members to SelfAdjointEigenSolver)
  The QR module now depends on Cholesky.
* Fix Transpose to correctly preserve the *TriangularBit.

Eigen/QR

@@ -2,6 +2,7 @@
 #define EIGEN_QR_MODULE_H
 
 #include "Core"
+#include "Cholesky"
 
 // Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module
 #if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2)

Eigen/src/Core/Transpose.h

@@ -48,7 +48,12 @@ struct ei_traits<Transpose<MatrixType> >
     ColsAtCompileTime = MatrixType::RowsAtCompileTime,
     MaxRowsAtCompileTime = MatrixType::MaxColsAtCompileTime,
     MaxColsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
-    Flags = (int(_MatrixTypeNested::Flags) ^ RowMajorBit) & ~Like1DArrayBit,
+    Flags = (int(_MatrixTypeNested::Flags) ^ RowMajorBit)
+          & ~( Like1DArrayBit | LowerTriangularBit | UpperTriangularBit)
+          | (int(_MatrixTypeNested::Flags)&UpperTriangularBit ? LowerTriangularBit : 0)
+          | (int(_MatrixTypeNested::Flags)&LowerTriangularBit ? UpperTriangularBit : 0),
+    // Note that the above test cannot be simplified using ^ because a diagonal matrix
+    // has both LowerTriangularBit and UpperTriangularBit and both must be preserved.
     CoeffReadCost = _MatrixTypeNested::CoeffReadCost
   };
 };

Eigen/src/QR/EigenSolver.h

@@ -72,7 +72,6 @@ template<typename _MatrixType> class EigenSolver
     EigenvalueType m_eivalues;
 };
 
-
 template<typename MatrixType>
 void EigenSolver<MatrixType>::compute(const MatrixType& matrix)
 {

Eigen/src/QR/QrInstantiations.cpp

@@ -26,6 +26,8 @@
 #define EIGEN_EXTERN_INSTANTIATIONS
 #endif
 #include "../../Core"
+// commented because of -pedantic
+// #include "../../Cholesky"
 #undef EIGEN_EXTERN_INSTANTIATIONS
 
 #include "../../QR"

Eigen/src/QR/SelfAdjointEigenSolver.h

@@ -49,6 +49,11 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
     typedef Tridiagonalization<MatrixType> TridiagonalizationType;
 
+    /** Constructors computing the eigenvalues of the selfadjoint matrix \a matrix,
+      * as well as the eigenvectors if \a computeEigenvectors is true.
+      *
+      * \sa compute(MatrixType,bool), SelfAdjointEigenSolver(MatrixType,MatrixType,bool)
+      */
     SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
       : m_eivec(matrix.rows(), matrix.cols()),
         m_eivalues(matrix.cols())
@@ -56,8 +61,24 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       compute(matrix, computeEigenvectors);
     }
 
+    /** Constructors computing the eigenvalues of the generalized eigen problem
+      * \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$ . The eigenvectors
+      * are computed if \a computeEigenvectors is true.
+      *
+      * \sa compute(MatrixType,MatrixType,bool), SelfAdjointEigenSolver(MatrixType,bool)
+      */
+    SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
+      : m_eivec(matA.rows(), matA.cols()),
+        m_eivalues(matA.cols())
+    {
+      compute(matA, matB, computeEigenvectors);
+    }
+
     void compute(const MatrixType& matrix, bool computeEigenvectors = true);
 
+    void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true);
+
     MatrixType eigenvectors(void) const
     {
       #ifndef NDEBUG
@@ -76,6 +97,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     #endif
 };
 
+#ifndef EIGEN_HIDE_HEAVY_CODE
+
 // from Golub's "Matrix Computations", algorithm 5.1.3
 template<typename Scalar>
 static void ei_givens_rotation(Scalar a, Scalar b, Scalar& c, Scalar& s)
@@ -117,6 +140,11 @@ static void ei_givens_rotation(Scalar a, Scalar b, Scalar& c, Scalar& s)
 template<typename RealScalar, typename Scalar>
 static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n);
 
+/** Computes the eigenvalues of the selfadjoint matrix \a matrix,
+  * as well as the eigenvectors if \a computeEigenvectors is true.
+  *
+  * \sa SelfAdjointEigenSolver(MatrixType,bool), compute(MatrixType,MatrixType,bool)
+  */
 template<typename MatrixType>
 void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
 {
@@ -170,6 +198,39 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool
   }
 }
 
+/** Computes the eigenvalues of the generalized eigen problem
+  * \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$ . The eigenvectors
+  * are computed if \a computeEigenvectors is true.
+  *
+  * \sa SelfAdjointEigenSolver(MatrixType,MatrixType,bool), compute(MatrixType,bool)
+  */
+template<typename MatrixType>
+void 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 = U'U
+  Cholesky<MatrixType> cholB(matB);
+
+  // compute C = inv(U') A inv(U)
+  MatrixType matC = cholB.matrixL().inverseProduct(matA);
+  // FIXME since we currently do not support A * inv(U),
+  // let's do (inv(U') A')' :
+  matC = (cholB.matrixL().inverseProduct(matC.adjoint())).adjoint();
+
+  compute(matC, computeEigenvectors);
+
+  if (computeEigenvectors)
+  {
+    // transform back the eigen vectors: evecs = inv(U) * evecs
+    m_eivec = cholB.matrixL().adjoint().template marked<Upper>().inverseProduct(m_eivec);
+  }
+}
+
+#endif // EIGEN_HIDE_HEAVY_CODE
+
 /** \qr_module
   *
   * \returns a vector listing the eigenvalues of this matrix.

test/eigensolver.cpp

@@ -36,10 +36,16 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
   typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex;
 
   MatrixType a = MatrixType::random(rows,cols);
-  MatrixType covMat =  a.adjoint() * a;
+  MatrixType symmA =  a.adjoint() * a;
 
-  SelfAdjointEigenSolver<MatrixType> eiSymm(covMat);
+  SelfAdjointEigenSolver<MatrixType> eiSymm(symmA);
-  VERIFY_IS_APPROX(covMat * eiSymm.eigenvectors(), (eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal().eval()));
+  VERIFY_IS_APPROX(symmA * eiSymm.eigenvectors(), (eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal().eval()));
+
+  // generalized eigen problem Ax = lBx
+  MatrixType b = MatrixType::random(rows,cols);
+  MatrixType symmB =  b.adjoint() * b;
+  eiSymm.compute(symmA,symmB);
+  VERIFY_IS_APPROX(symmA * eiSymm.eigenvectors(), symmB * (eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal().eval()));
 
 //   EigenSolver<MatrixType> eiNotSymmButSymm(covMat);
 //   VERIFY_IS_APPROX((covMat.template cast<Complex>()) * (eiNotSymmButSymm.eigenvectors().template cast<Complex>()),