From 139529e97b76c250c13474a1ebc9f6dbf27b6115 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 14 Nov 2008 09:55:25 +0000 Subject: [PATCH] * add .imag() function * fix a very old bug in EigenSolver that I had completely forgotten (thanks to Timothy to refresh my mind) * fix doc of Matrix::Map --- Eigen/src/Core/CwiseUnaryOp.h | 14 +++++++++----- Eigen/src/Core/Functors.h | 16 +++++++++++++++- Eigen/src/Core/Matrix.h | 6 +++--- Eigen/src/Core/MatrixBase.h | 3 +++ Eigen/src/Core/util/ForwardDeclarations.h | 1 + Eigen/src/QR/EigenSolver.h | 17 ++++++++--------- test/eigensolver.cpp | 21 +++------------------ 7 files changed, 42 insertions(+), 36 deletions(-) diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h index a50a9c30d..4ceb6d98a 100644 --- a/Eigen/src/Core/CwiseUnaryOp.h +++ b/Eigen/src/Core/CwiseUnaryOp.h @@ -174,13 +174,17 @@ MatrixBase::conjugate() const /** \returns an expression of the real part of \c *this. * - * \sa adjoint() */ + * \sa imag() */ template inline const typename MatrixBase::RealReturnType -MatrixBase::real() const -{ - return derived(); -} +MatrixBase::real() const { return derived(); } + +/** \returns an expression of the imaginary part of \c *this. + * + * \sa real() */ +template +inline const typename MatrixBase::ImagReturnType +MatrixBase::imag() const { return derived(); } /** \returns an expression of *this with the \a Scalar type casted to * \a NewScalar. diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index 4295aa4da..d68558e7b 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -240,7 +240,21 @@ struct ei_scalar_real_op EIGEN_EMPTY_STRUCT { }; template struct ei_functor_traits > -{ enum { Cost = 0, PacketAccess = false }; }; +{ enum { Cost = 0, PacketAccess = false }; }; + +/** \internal + * \brief Template functor to extract the imaginary part of a complex + * + * \sa class CwiseUnaryOp, MatrixBase::imag() + */ +template +struct ei_scalar_imag_op EIGEN_EMPTY_STRUCT { + typedef typename NumTraits::Real result_type; + inline result_type operator() (const Scalar& a) const { return ei_imag(a); } +}; +template +struct ei_functor_traits > +{ enum { Cost = 0, PacketAccess = false }; }; /** \internal * \brief Template functor to multiply a scalar by a fixed other one diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 335150318..ba24b4b19 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -437,14 +437,14 @@ class Matrix this->Base::swap(other); } - /** + /** \name Map * These are convenience functions returning Map objects. The Map() static functions return unaligned Map objects, * while the AlignedMap() functions return aligned Map objects and thus should be called only with 16-byte-aligned * \a data pointers. - * + * * \see class Map */ - //@} + //@{ inline static const UnalignedMapType Map(const Scalar* data) { return UnalignedMapType(data); } inline static UnalignedMapType Map(Scalar* data) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index dbf8da544..d231f973c 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -189,6 +189,8 @@ template class MatrixBase >::ret ConjugateReturnType; /** \internal the return type of MatrixBase::real() */ typedef CwiseUnaryOp, Derived> RealReturnType; + /** \internal the return type of MatrixBase::imag() */ + typedef CwiseUnaryOp, Derived> ImagReturnType; /** \internal the return type of MatrixBase::adjoint() */ typedef Eigen::Transpose::type> > AdjointReturnType; @@ -496,6 +498,7 @@ template class MatrixBase ConjugateReturnType conjugate() const; const RealReturnType real() const; + const ImagReturnType imag() const; template const CwiseUnaryOp unaryExpr(const CustomUnaryOp& func = CustomUnaryOp()) const; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 98d15b415..0d2a457b6 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -64,6 +64,7 @@ template struct ei_scalar_quotient_op; template struct ei_scalar_opposite_op; template struct ei_scalar_conjugate_op; template struct ei_scalar_real_op; +template struct ei_scalar_imag_op; template struct ei_scalar_abs_op; template struct ei_scalar_abs2_op; template struct ei_scalar_sqrt_op; diff --git a/Eigen/src/QR/EigenSolver.h b/Eigen/src/QR/EigenSolver.h index d7d891951..38a383f14 100644 --- a/Eigen/src/QR/EigenSolver.h +++ b/Eigen/src/QR/EigenSolver.h @@ -59,7 +59,7 @@ template class EigenSolver compute(matrix); } - + EigenvectorType eigenvectors(void) const; /** \returns a real matrix V of pseudo eigenvectors. @@ -200,18 +200,18 @@ void EigenSolver::orthes(MatrixType& matH, RealVectorType& ort) for (int m = low+1; m <= high-1; m++) { // Scale column. - Scalar scale = matH.block(m, m-1, high-m+1, 1).cwise().abs().sum(); + RealScalar scale = matH.block(m, m-1, high-m+1, 1).cwise().abs().sum(); if (scale != 0.0) { // Compute Householder transformation. - Scalar h = 0.0; + RealScalar h = 0.0; // FIXME could be rewritten, but this one looks better wrt cache for (int i = high; i >= m; i--) { ort.coeffRef(i) = matH.coeff(i,m-1)/scale; h += ort.coeff(i) * ort.coeff(i); } - Scalar g = ei_sqrt(h); + RealScalar g = ei_sqrt(h); if (ort.coeff(m) > 0) g = -g; h = h - ort.coeff(m) * g; @@ -247,7 +247,6 @@ void EigenSolver::orthes(MatrixType& matH, RealVectorType& ort) } } - // Complex scalar division. template std::complex cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi) @@ -298,7 +297,7 @@ void EigenSolver::hqr2(MatrixType& matH) m_eivalues.coeffRef(j).real() = matH.coeff(j,j); m_eivalues.coeffRef(j).imag() = 0.0; } - norm += matH.col(j).start(std::min(j+1,nn)).cwise().abs().sum(); + norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwise().abs().sum(); } // Outer loop over eigenvalue index @@ -571,7 +570,7 @@ void EigenSolver::hqr2(MatrixType& matH) for (int i = n-1; i >= 0; i--) { w = matH.coeff(i,i) - p; - r = (matH.row(i).end(nn-l) * matH.col(n).end(nn-l))(0,0); + r = (matH.row(i).segment(l,n-l+1) * matH.col(n).segment(l, n-l+1))(0,0); if (m_eivalues.coeff(i).imag() < 0.0) { @@ -630,8 +629,8 @@ void EigenSolver::hqr2(MatrixType& matH) for (int i = n-2; i >= 0; i--) { Scalar ra,sa,vr,vi; - ra = (matH.row(i).end(nn-l) * matH.col(n-1).end(nn-l)).lazy()(0,0); - sa = (matH.row(i).end(nn-l) * matH.col(n).end(nn-l)).lazy()(0,0); + ra = (matH.block(i,l, 1, n-l+1) * matH.block(l,n-1, n-l+1, 1)).lazy()(0,0); + sa = (matH.block(i,l, 1, n-l+1) * matH.block(l,n, n-l+1, 1)).lazy()(0,0); w = matH.coeff(i,i) - p; if (m_eivalues.coeff(i).imag() < 0.0) diff --git a/test/eigensolver.cpp b/test/eigensolver.cpp index fed6ba9ba..aeec1dc97 100644 --- a/test/eigensolver.cpp +++ b/test/eigensolver.cpp @@ -130,27 +130,13 @@ template void eigensolver(const MatrixType& m) MatrixType a1 = MatrixType::Random(rows,cols); MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1; -// MatrixType b = MatrixType::Random(rows,cols); -// MatrixType b1 = MatrixType::Random(rows,cols); -// MatrixType symmB = b.adjoint() * b + b1.adjoint() * b1; - EigenSolver ei0(symmA); + VERIFY_IS_APPROX(symmA * ei0.pseudoEigenvectors(), ei0.pseudoEigenvectors() * ei0.pseudoEigenvalueMatrix()); VERIFY_IS_APPROX((symmA.template cast()) * (ei0.pseudoEigenvectors().template cast()), (ei0.pseudoEigenvectors().template cast()) * (ei0.eigenvalues().asDiagonal())); -// a = a /*+ symmA*/; EigenSolver ei1(a); - - IOFormat OctaveFmt(4, AlignCols, ", ", ";\n", "", "", "[", "]"); -// std::cerr << "==============\n" << a.format(OctaveFmt) << "\n\n" << ei1.eigenvalues().transpose() << "\n\n"; -// std::cerr << a * ei1.pseudoEigenvectors() << "\n\n" << ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix() << "\n\n\n"; - -// VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix()); - - -// std::cerr << a.format(OctaveFmt) << "\n\n"; -// std::cerr << ei1.eigenvectors().format(OctaveFmt) << "\n\n"; -// std::cerr << a.template cast() * ei1.eigenvectors() << "\n\n" << ei1.eigenvectors() * ei1.eigenvalues().asDiagonal().eval() << "\n\n"; + VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix()); VERIFY_IS_APPROX(a.template cast() * ei1.eigenvectors(), ei1.eigenvectors() * ei1.eigenvalues().asDiagonal().eval()); @@ -167,8 +153,7 @@ void test_eigensolver() CALL_SUBTEST( selfadjointeigensolver(MatrixXd(19,19)) ); CALL_SUBTEST( eigensolver(Matrix4f()) ); - // FIXME the test fails for larger matrices -// CALL_SUBTEST( eigensolver(MatrixXd(7,7)) ); + CALL_SUBTEST( eigensolver(MatrixXd(17,17)) ); } }