diff options
| -rw-r--r-- | Eigen/QR | 1 | ||||
| -rw-r--r-- | Eigen/src/Core/Transpose.h | 7 | ||||
| -rw-r--r-- | Eigen/src/QR/EigenSolver.h | 1 | ||||
| -rw-r--r-- | Eigen/src/QR/QrInstantiations.cpp | 2 | ||||
| -rw-r--r-- | Eigen/src/QR/SelfAdjointEigenSolver.h | 61 | ||||
| -rw-r--r-- | test/eigensolver.cpp | 12 |
6 files changed, 79 insertions, 5 deletions
@@ -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) diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 5c9cd54ac..1a4843973 100644 --- a/Eigen/src/Core/Transpose.h +++ b/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 }; }; diff --git a/Eigen/src/QR/EigenSolver.h b/Eigen/src/QR/EigenSolver.h index 2433781cd..a09bb210a 100644 --- a/Eigen/src/QR/EigenSolver.h +++ b/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) { diff --git a/Eigen/src/QR/QrInstantiations.cpp b/Eigen/src/QR/QrInstantiations.cpp index dacb05d3d..3b2594e34 100644 --- a/Eigen/src/QR/QrInstantiations.cpp +++ b/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" diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h index 262eba4bf..62d3a38bb 100644 --- a/Eigen/src/QR/SelfAdjointEigenSolver.h +++ b/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. diff --git a/test/eigensolver.cpp b/test/eigensolver.cpp index cb519657c..89cbfe8fc 100644 --- a/test/eigensolver.cpp +++ b/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); - VERIFY_IS_APPROX(covMat * eiSymm.eigenvectors(), (eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal().eval())); + SelfAdjointEigenSolver<MatrixType> eiSymm(symmA); + 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>()), |