From 9a4a08da46fb88a64550c65c43231df7f5193276 Mon Sep 17 00:00:00 2001 From: Manuel Yguel Date: Thu, 25 Mar 2010 03:15:05 +0100 Subject: [PATCH 01/19] Creation of the Polynomials module with the following features: * convenient functions: - Horner and stabilized Horner evaluation - polynomial coefficients from a set of given roots - Cauchy bounds * a QR based polynomial solver --- unsupported/Eigen/CMakeLists.txt | 2 +- unsupported/Eigen/Polynomials | 132 ++++++ unsupported/Eigen/src/CMakeLists.txt | 1 + .../Eigen/src/Polynomials/CMakeLists.txt | 6 + unsupported/Eigen/src/Polynomials/Companion.h | 281 +++++++++++++ .../Eigen/src/Polynomials/PolynomialSolver.h | 395 ++++++++++++++++++ .../Eigen/src/Polynomials/PolynomialUtils.h | 153 +++++++ .../doc/examples/PolynomialSolver1.cpp | 53 +++ unsupported/doc/examples/PolynomialUtils1.cpp | 20 + unsupported/test/CMakeLists.txt | 16 + 10 files changed, 1058 insertions(+), 1 deletion(-) create mode 100644 unsupported/Eigen/Polynomials create mode 100644 unsupported/Eigen/src/Polynomials/CMakeLists.txt create mode 100644 unsupported/Eigen/src/Polynomials/Companion.h create mode 100644 unsupported/Eigen/src/Polynomials/PolynomialSolver.h create mode 100644 unsupported/Eigen/src/Polynomials/PolynomialUtils.h create mode 100644 unsupported/doc/examples/PolynomialSolver1.cpp create mode 100644 unsupported/doc/examples/PolynomialUtils1.cpp diff --git a/unsupported/Eigen/CMakeLists.txt b/unsupported/Eigen/CMakeLists.txt index f9d79c51a..87cc4be1e 100644 --- a/unsupported/Eigen/CMakeLists.txt +++ b/unsupported/Eigen/CMakeLists.txt @@ -1,4 +1,4 @@ -set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3) +set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3 Polynomials) install(FILES ${Eigen_HEADERS} diff --git a/unsupported/Eigen/Polynomials b/unsupported/Eigen/Polynomials new file mode 100644 index 000000000..9e1f6b759 --- /dev/null +++ b/unsupported/Eigen/Polynomials @@ -0,0 +1,132 @@ +#ifndef EIGEN_POLYNOMIALS_MODULE_H +#define EIGEN_POLYNOMIALS_MODULE_H + +#include + +#include + +#include + +// Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module +#if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2) + #ifndef EIGEN_HIDE_HEAVY_CODE + #define EIGEN_HIDE_HEAVY_CODE + #endif +#elif defined EIGEN_HIDE_HEAVY_CODE + #undef EIGEN_HIDE_HEAVY_CODE +#endif + +namespace Eigen { + +/** \ingroup Unsupported_modules + * \defgroup Polynomials_Module Polynomials module + * + * \nonstableyet + * + * \brief This module provides a QR based polynomial solver. + * + * To use this module, add + * \code + * #include + * \endcode + * at the start of your source file. + */ + +#include "src/Polynomials/PolynomialUtils.h" +#include "src/Polynomials/Companion.h" +#include "src/Polynomials/PolynomialSolver.h" + +/** + \page polynomials Polynomials defines functions for dealing with polynomials + and a QR based polynomial solver. + \ingroup Polynomials_Module + + The remainder of the page documents first the functions for evaluating, computing + polynomials, computing estimates about polynomials and next the QR based polynomial + solver. + + \section polynomialUtils convenient functions to deal with polynomials + \subsection roots_to_monicPolynomial + The function + \code + void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly ) + \endcode + computes the coefficients \f$ a_i \f$ of + + \f$ p(x) = a_0 + a_{1}x + ... + a_{n-1}x^{n-1} + x^n \f$ + + where \f$ p \f$ is known through its roots i.e. \f$ p(x) = (x-r_1)(x-r_2)...(x-r_n) \f$. + + \subsection poly_eval + The function + \code + T poly_eval( const Polynomials& poly, const T& x ) + \endcode + evaluates a polynomial at a given point using stabilized Hörner method. + + The following code computes the coefficients in the monomial basis of the monic polynomial given through its roots then evaluate it at those roots. + + \include PolynomialUtils1.cpp + Output: \verbinclude PolynomialUtils1.out + + \subsection Cauchy bounds + The function + \code + Real cauchy_max_bound( const Polynomial& poly ) + \endcode + provides a maximum bound (the Cauchy one: \f$C(p)\f$) for the absolute value of a root of the given polynomial i.e. + \f$ \forall r_i \f$ root of \f$ p(x) = \sum_{k=0}^d a_k x^k \f$, + \f$ |r_i| \le C(p) = \sum_{k=0}^{d} \left | \frac{a_k}{a_d} \right | \f$ + The leading coefficient \f$ p \f$: should be non zero \f$a_d \neq 0\f$. + + + The function + \code + Real cauchy_min_bound( const Polynomial& poly ) + \endcode + provides a minimum bound (the Cauchy one: \f$c(p)\f$) for the absolute value of a non zero root of the given polynomial i.e. + \f$ \forall r_i \neq 0 \f$ root of \f$ p(x) = \sum_{k=0}^d a_k x^k \f$, + \f$ |r_i| \ge c(p) = \left( \sum_{k=0}^{d} \left | \frac{a_k}{a_0} \right | \right)^{-1} \f$ + + + + + \section QR polynomial solver class + Computes the complex roots of a polynomial by computing the eigenvalues of the associated companion matrix with the QR algorithm. + + The roots of \f$ p(x) = a_0 + a_1 x + a_2 x^2 + a_{3} x^3 + x^4 \f$ are the eigenvalues of + \f$ + \left [ + \begin{array}{cccc} + 0 & 0 & 0 & a_0 \\ + 1 & 0 & 0 & a_1 \\ + 0 & 1 & 0 & a_2 \\ + 0 & 0 & 1 & a_3 + \end{array} \right ] + \f$ + + However, the QR algorithm is not guaranteed to converge when there are several eigenvalues with same modulus. + + Therefore the current polynomial solver is guaranteed to provide a correct result only when the complex roots \f$r_1,r_2,...,r_d\f$ have distinct moduli i.e. + + \f$ \forall i,j \in [1;d],~ \| r_i \| \neq \| r_j \| \f$. + + With 32bit (float) floating types this problem shows up frequently. + However, almost always, correct accuracy is reached even in these cases for 64bit + (double) floating types and small polynomial degree (<20). + + \include PolynomialSolver1.cpp + In the example a polynomial with almost conjugate roots is provided to the solver. + Those roots have almost same module therefore the QR algorithm failed to converge: the accuracy + of the last root is bad. + + This problem is less visible with double. + Output: \verbinclude PolynomialSolver1.out +*/ + +} // namespace Eigen + +#include + +#endif // EIGEN_POLYNOMIALS_MODULE_H +/* vim: set filetype=cpp et sw=2 ts=2 ai: */ diff --git a/unsupported/Eigen/src/CMakeLists.txt b/unsupported/Eigen/src/CMakeLists.txt index a0870d3af..035a84ca8 100644 --- a/unsupported/Eigen/src/CMakeLists.txt +++ b/unsupported/Eigen/src/CMakeLists.txt @@ -5,3 +5,4 @@ ADD_SUBDIRECTORY(MoreVectorization) # ADD_SUBDIRECTORY(FFT) # ADD_SUBDIRECTORY(Skyline) ADD_SUBDIRECTORY(MatrixFunctions) +ADD_SUBDIRECTORY(Polynomials) diff --git a/unsupported/Eigen/src/Polynomials/CMakeLists.txt b/unsupported/Eigen/src/Polynomials/CMakeLists.txt new file mode 100644 index 000000000..51f13f3cb --- /dev/null +++ b/unsupported/Eigen/src/Polynomials/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_Polynomials_SRCS "*.h") + +INSTALL(FILES + ${Eigen_Polynomials_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/Polynomials COMPONENT Devel + ) diff --git a/unsupported/Eigen/src/Polynomials/Companion.h b/unsupported/Eigen/src/Polynomials/Companion.h new file mode 100644 index 000000000..7c9e9c518 --- /dev/null +++ b/unsupported/Eigen/src/Polynomials/Companion.h @@ -0,0 +1,281 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Manuel Yguel +// +// 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 . + +#ifndef EIGEN_COMPANION_H +#define EIGEN_COMPANION_H + +// This file requires the user to include +// * Eigen/Core +// * Eigen/src/PolynomialSolver.h + +#ifndef EIGEN_PARSED_BY_DOXYGEN + +template +T ei_radix(){ return 2; } + +template +T ei_radix2(){ return ei_radix()*ei_radix(); } + +template +struct ei_decrement_if_fixed_size +{ + enum { + ret = (Size == Dynamic) ? Dynamic : Size-1 }; +}; + +#endif + +template< typename _Scalar, int _Deg > +class ei_companion +{ + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) + + enum { + Deg = _Deg, + Deg_1=ei_decrement_if_fixed_size::ret + }; + + typedef _Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix RightColumn; + //typedef DiagonalMatrix< Scalar, Deg_1, Deg_1 > BottomLeftDiagonal; + typedef Matrix BottomLeftDiagonal; + + typedef Matrix DenseCompanionMatrixType; + typedef Matrix< Scalar, _Deg, Deg_1 > LeftBlock; + typedef Matrix< Scalar, Deg_1, Deg_1 > BottomLeftBlock; + typedef Matrix< Scalar, 1, Deg_1 > LeftBlockFirstRow; + + public: + EIGEN_STRONG_INLINE const _Scalar operator()( int row, int col ) const + { + if( m_bl_diag.rows() > col ) + { + if( 0 < row ){ return m_bl_diag[col]; } + else{ return 0; } + } + else{ return m_monic[row]; } + } + + public: + template + void setPolynomial( const VectorType& poly ) + { + const int deg = poly.size()-1; + m_monic = -1/poly[deg] * poly.head(deg); + //m_bl_diag.setIdentity( deg-1 ); + m_bl_diag.setOnes(deg-1); + } + + template + ei_companion( const VectorType& poly ){ + setPolynomial( poly ); } + + public: + DenseCompanionMatrixType denseMatrix() const + { + const int deg = m_monic.size(); + const int deg_1 = deg-1; + DenseCompanionMatrixType companion(deg,deg); + companion << + ( LeftBlock(deg,deg_1) + << LeftBlockFirstRow::Zero(1,deg_1), + BottomLeftBlock::Identity(deg-1,deg-1)*m_bl_diag.asDiagonal() ).finished() + , m_monic; + return companion; + } + + + + protected: + /** Helper function for the balancing algorithm. + * \returns true if the row and the column, having colNorm and rowNorm + * as norms, are balanced, false otherwise. + * colB and rowB are repectively the multipliers for + * the column and the row in order to balance them. + * */ + bool balanced( Scalar colNorm, Scalar rowNorm, + bool& isBalanced, Scalar& colB, Scalar& rowB ); + + /** Helper function for the balancing algorithm. + * \returns true if the row and the column, having colNorm and rowNorm + * as norms, are balanced, false otherwise. + * colB and rowB are repectively the multipliers for + * the column and the row in order to balance them. + * */ + bool balancedR( Scalar colNorm, Scalar rowNorm, + bool& isBalanced, Scalar& colB, Scalar& rowB ); + + public: + /** + * Balancing algorithm from B. N. PARLETT and C. REINSCH (1969) + * "Balancing a matrix for calculation of eigenvalues and eigenvectors" + * adapted to the case of companion matrices. + * A matrix with non zero row and non zero column is balanced + * for a certain norm if the i-th row and the i-th column + * have same norm for all i. + */ + void balance(); + + protected: + RightColumn m_monic; + BottomLeftDiagonal m_bl_diag; +}; + + + +template< typename _Scalar, int _Deg > +inline +bool ei_companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm, + bool& isBalanced, Scalar& colB, Scalar& rowB ) +{ + if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; } + else + { + //To find the balancing coefficients, if the radix is 2, + //one finds \f$ \sigma \f$ such that + // \f$ 2^{2\sigma-1} < rowNorm / colNorm \le 2^{2\sigma+1} \f$ + // then the balancing coefficient for the row is \f$ 1/2^{\sigma} \f$ + // and the balancing coefficient for the column is \f$ 2^{\sigma} \f$ + rowB = rowNorm / ei_radix(); + colB = Scalar(1); + const Scalar s = colNorm + rowNorm; + + while (colNorm < rowB) + { + colB *= ei_radix(); + colNorm *= ei_radix2(); + } + + rowB = rowNorm * ei_radix(); + + while (colNorm >= rowB) + { + colB /= ei_radix(); + colNorm /= ei_radix2(); + } + + //This line is used to avoid insubstantial balancing + if ((rowNorm + colNorm) < Scalar(0.95) * s * colB) + { + isBalanced = false; + rowB = Scalar(1) / colB; + return false; + } + else{ + return true; } + } +} + +template< typename _Scalar, int _Deg > +inline +bool ei_companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm, + bool& isBalanced, Scalar& colB, Scalar& rowB ) +{ + if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; } + else + { + /** + * Set the norm of the column and the row to the geometric mean + * of the row and column norm + */ + const _Scalar q = colNorm/rowNorm; + if( !ei_isApprox( q, _Scalar(1) ) ) + { + rowB = ei_sqrt( colNorm/rowNorm ); + colB = Scalar(1)/rowB; + + isBalanced = false; + return false; + } + else{ + return true; } + } +} + + +template< typename _Scalar, int _Deg > +void ei_companion<_Scalar,_Deg>::balance() +{ + EIGEN_STATIC_ASSERT( 1 < Deg, YOU_MADE_A_PROGRAMMING_MISTAKE ); + const int deg = m_monic.size(); + const int deg_1 = deg-1; + + bool hasConverged=false; + while( !hasConverged ) + { + hasConverged = true; + Scalar colNorm,rowNorm; + Scalar colB,rowB; + + //First row, first column excluding the diagonal + //============================================== + colNorm = ei_abs(m_bl_diag[0]); + rowNorm = ei_abs(m_monic[0]); + + //Compute balancing of the row and the column + if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) ) + { + m_bl_diag[0] *= colB; + m_monic[0] *= rowB; + } + + //Middle rows and columns excluding the diagonal + //============================================== + for( int i=1; i headMonic( m_monic, 0, deg_1 ); + colNorm = headMonic.array().abs().sum(); + rowNorm = ei_abs( m_bl_diag[ebl] ); + + //Compute balancing of the row and the column + if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) ) + { + headMonic *= colB; + m_bl_diag[ebl] *= rowB; + } + } +} + + +#endif // EIGEN_COMPANION_H diff --git a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h new file mode 100644 index 000000000..49a5ffffa --- /dev/null +++ b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h @@ -0,0 +1,395 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Manuel Yguel +// +// 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 . + +#ifndef EIGEN_POLYNOMIAL_SOLVER_H +#define EIGEN_POLYNOMIAL_SOLVER_H + +/** \ingroup Polynomials_Module + * \class PolynomialSolverBase. + * + * \brief Defined to be inherited by polynomial solvers: it provides + * convenient methods such as + * - real roots, + * - greatest, smallest complex roots, + * - real roots with greatest, smallest absolute real value, + * - greatest, smallest real roots. + * + * It stores the set of roots as a vector of complexes. + * + */ +template< typename _Scalar, int _Deg > +class PolynomialSolverBase +{ + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) + + typedef _Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef std::complex RootType; + typedef Matrix RootsType; + + protected: + template< typename OtherPolynomial > + inline void setPolynomial( const OtherPolynomial& poly ){ + m_roots.resize(poly.size()); } + + public: + template< typename OtherPolynomial > + inline PolynomialSolverBase( const OtherPolynomial& poly ){ + setPolynomial( poly() ); } + + inline PolynomialSolverBase(){} + + public: + /** \returns the complex roots of the polynomial */ + inline const RootsType& roots() const { return m_roots; } + + public: + /** Clear and fills the back insertion sequence with the real roots of the polynomial + * i.e. the real part of the complex roots that have an imaginary part which + * absolute value is smaller than absImaginaryThreshold. + * absImaginaryThreshold takes the dummy_precision associated + * with the _Scalar template parameter of the PolynomialSolver class as the default value. + * + * \param[out] bi_seq : the back insertion sequence (stl concept) + * \param[in] absImaginaryThreshold : the maximum bound of the imaginary part of a complex + * number that is considered as real. + * */ + template + inline void realRoots( Stl_back_insertion_sequence& bi_seq, + const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const + { + bi_seq.clear(); + for( int i=0; i + inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const + { + int res=0; + RealScalar norm2 = ei_abs2( m_roots[0] ); + for( int i=1; i greater; + return selectComplexRoot_withRespectToNorm( greater ); + } + + /** + * \returns the complex root with smallest norm. + */ + inline const RootType& smallestRoot() const + { + std::less less; + return selectComplexRoot_withRespectToNorm( less ); + } + + protected: + template + inline const RealScalar& selectRealRoot_withRespectToAbsRealPart( + squaredRealPartBinaryPredicate& pred, + bool& hasArealRoot, + const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const + { + hasArealRoot = false; + int res=0; + RealScalar abs2; + + for( int i=0; i + inline const RealScalar& selectRealRoot_withRespectToRealPart( + RealPartBinaryPredicate& pred, + bool& hasArealRoot, + const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const + { + hasArealRoot = false; + int res=0; + RealScalar val; + + for( int i=0; i::dummy_precision() ) const + { + std::greater greater; + return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold ); + } + + + /** + * \returns a real root with smallest absolute magnitude. + * A real root is defined as the real part of a complex root with absolute imaginary + * part smallest than absImaginaryThreshold. + * absImaginaryThreshold takes the dummy_precision associated + * with the _Scalar template parameter of the PolynomialSolver class as the default value. + * If no real root is found the boolean hasArealRoot is set to false and the real part of + * the root with smallest absolute imaginary part is returned instead. + * + * \param[out] hasArealRoot : boolean true if a real root is found according to the + * absImaginaryThreshold criterion, false otherwise. + * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide + * whether or not a root is real. + */ + inline const RealScalar& absSmallestRealRoot( + bool& hasArealRoot, + const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const + { + std::less less; + return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold ); + } + + + /** + * \returns the real root with greatest value. + * A real root is defined as the real part of a complex root with absolute imaginary + * part smallest than absImaginaryThreshold. + * absImaginaryThreshold takes the dummy_precision associated + * with the _Scalar template parameter of the PolynomialSolver class as the default value. + * If no real root is found the boolean hasArealRoot is set to false and the real part of + * the root with smallest absolute imaginary part is returned instead. + * + * \param[out] hasArealRoot : boolean true if a real root is found according to the + * absImaginaryThreshold criterion, false otherwise. + * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide + * whether or not a root is real. + */ + inline const RealScalar& greatestRealRoot( + bool& hasArealRoot, + const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const + { + std::greater greater; + return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold ); + } + + + /** + * \returns the real root with smallest value. + * A real root is defined as the real part of a complex root with absolute imaginary + * part smallest than absImaginaryThreshold. + * absImaginaryThreshold takes the dummy_precision associated + * with the _Scalar template parameter of the PolynomialSolver class as the default value. + * If no real root is found the boolean hasArealRoot is set to false and the real part of + * the root with smallest absolute imaginary part is returned instead. + * + * \param[out] hasArealRoot : boolean true if a real root is found according to the + * absImaginaryThreshold criterion, false otherwise. + * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide + * whether or not a root is real. + */ + inline const RealScalar& smallestRealRoot( + bool& hasArealRoot, + const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const + { + std::less less; + return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold ); + } + + protected: + RootsType m_roots; +}; + +#define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \ + typedef typename BASE::Scalar Scalar; \ + typedef typename BASE::RealScalar RealScalar; \ + typedef typename BASE::RootType RootType; \ + typedef typename BASE::RootsType RootsType; + + + +/** \ingroup Polynomials_Module + * + * \class PolynomialSolver + * + * \brief A polynomial solver + * + * Computes the complex roots of a real polynomial. + * + * \param _Scalar the scalar type, i.e., the type of the polynomial coefficients + * \param _Deg the degree of the polynomial, can be a compile time value or Dynamic. + * Notice that the number of polynomial coefficients is _Deg+1. + * + * This class implements a polynomial solver and provides convenient methods such as + * - real roots, + * - greatest, smallest complex roots, + * - real roots with greatest, smallest absolute real value. + * - greatest, smallest real roots. + * + * WARNING: this polynomial solver is experimental, part of the unsuported Eigen modules. + * + * + * Currently a QR algorithm is used to compute the eigenvalues of the companion matrix of + * the polynomial to compute its roots. + * This supposes that the complex moduli of the roots are all distinct: e.g. there should + * be no multiple roots or conjugate roots for instance. + * With 32bit (float) floating types this problem shows up frequently. + * However, almost always, correct accuracy is reached even in these cases for 64bit + * (double) floating types and small polynomial degree (<20). + */ +template< typename _Scalar, int _Deg > +class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> +{ + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) + + typedef PolynomialSolverBase<_Scalar,_Deg> PS_Base; + EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) + + typedef Matrix CompanionMatrixType; + typedef EigenSolver EigenSolverType; + + public: + /** Computes the complex roots of a new polynomial. */ + template< typename OtherPolynomial > + void compute( const OtherPolynomial& poly ) + { + assert( Scalar(0) != poly[poly.size()-1] ); + ei_companion companion( poly ); + companion.balance(); + m_eigenSolver.compute( companion.denseMatrix() ); + m_roots = m_eigenSolver.eigenvalues(); + } + + public: + template< typename OtherPolynomial > + inline PolynomialSolver( const OtherPolynomial& poly ){ + compute( poly ); } + + inline PolynomialSolver(){} + + protected: + using PS_Base::m_roots; + EigenSolverType m_eigenSolver; +}; + + +template< typename _Scalar > +class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1> +{ + public: + typedef PolynomialSolverBase<_Scalar,1> PS_Base; + EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) + + public: + /** Computes the complex roots of a new polynomial. */ + template< typename OtherPolynomial > + void compute( const OtherPolynomial& poly ) + { + assert( Scalar(0) != poly[poly.size()-1] ); + m_roots[0] = -poly[0]/poly[poly.size()-1]; + } + + protected: + using PS_Base::m_roots; +}; + +#endif // EIGEN_POLYNOMIAL_SOLVER_H diff --git a/unsupported/Eigen/src/Polynomials/PolynomialUtils.h b/unsupported/Eigen/src/Polynomials/PolynomialUtils.h new file mode 100644 index 000000000..d78821f90 --- /dev/null +++ b/unsupported/Eigen/src/Polynomials/PolynomialUtils.h @@ -0,0 +1,153 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Manuel Yguel +// +// 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 . + +#ifndef EIGEN_POLYNOMIAL_UTILS_H +#define EIGEN_POLYNOMIAL_UTILS_H + +/** \ingroup Polynomials_Module + * \returns the evaluation of the polynomial at x using Horner algorithm. + * + * \param[in] poly : the vector of coefficients of the polynomial ordered + * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial + * e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$. + * \param[in] x : the value to evaluate the polynomial at. + * + * Note for stability: + *
\f$ |x| \le 1 \f$
+ */ +template +inline +T poly_eval_horner( const Polynomials& poly, const T& x ) +{ + T val=poly[poly.size()-1]; + for( int i=poly.size()-2; i>=0; --i ){ + val = val*x + poly[i]; } + return val; +} + +/** \ingroup Polynomials_Module + * \returns the evaluation of the polynomial at x using stabilized Horner algorithm. + * + * \param[in] poly : the vector of coefficients of the polynomial ordered + * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial + * e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$. + * \param[in] x : the value to evaluate the polynomial at. + */ +template +inline +T poly_eval( const Polynomials& poly, const T& x ) +{ + typedef typename NumTraits::Real Real; + + if( ei_abs2( x ) <= Real(1) ){ + return poly_eval_horner( poly, x ); } + else + { + T val=poly[0]; + T inv_x = T(1)/x; + for( int i=1; iPrecondition: + *
the leading coefficient of the input polynomial poly must be non zero
+ */ +template +inline +typename NumTraits::Real cauchy_max_bound( const Polynomial& poly ) +{ + typedef typename Polynomial::Scalar Scalar; + typedef typename NumTraits::Real Real; + + assert( Scalar(0) != poly[poly.size()-1] ); + const Scalar inv_leading_coeff = Scalar(1)/poly[poly.size()-1]; + Real cb(0); + + for( int i=0; i +inline +typename NumTraits::Real cauchy_min_bound( const Polynomial& poly ) +{ + typedef typename Polynomial::Scalar Scalar; + typedef typename NumTraits::Real Real; + + int i=0; + while( i +void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly ) +{ + + typedef typename Polynomial::Scalar Scalar; + + poly.setZero( rv.size()+1 ); + poly[0] = -rv[0]; poly[1] = Scalar(1); + for( int i=1; i<(int)rv.size(); ++i ) + { + for( int j=i+1; j>0; --j ){ poly[j] = poly[j-1] - rv[i]*poly[j]; } + poly[0] = -rv[i]*poly[0]; + } +} + + +#endif // EIGEN_POLYNOMIAL_UTILS_H diff --git a/unsupported/doc/examples/PolynomialSolver1.cpp b/unsupported/doc/examples/PolynomialSolver1.cpp new file mode 100644 index 000000000..c875c9361 --- /dev/null +++ b/unsupported/doc/examples/PolynomialSolver1.cpp @@ -0,0 +1,53 @@ +#include +#include +#include + +using namespace Eigen; +using namespace std; + +int main() +{ + typedef Matrix Vector5d; + + Vector5d roots = Vector5d::Random(); + cout << "Roots: " << roots.transpose() << endl; + Eigen::Matrix polynomial; + roots_to_monicPolynomial( roots, polynomial ); + + PolynomialSolver psolve( polynomial ); + cout << "Complex roots: " << psolve.roots().transpose() << endl; + + std::vector realRoots; + psolve.realRoots( realRoots ); + Map mapRR( &realRoots[0] ); + cout << "Real roots: " << mapRR.transpose() << endl; + + cout << endl; + cout << "Illustration of the convergence problem with the QR algorithm: " << endl; + cout << "---------------------------------------------------------------" << endl; + Eigen::Matrix hardCase_polynomial; + hardCase_polynomial << + -0.957, 0.9219, 0.3516, 0.9453, -0.4023, -0.5508, -0.03125; + cout << "Hard case polynomial defined by floats: " << hardCase_polynomial.transpose() << endl; + PolynomialSolver psolvef( hardCase_polynomial ); + cout << "Complex roots: " << psolvef.roots().transpose() << endl; + Eigen::Matrix evals; + for( int i=0; i<6; ++i ){ evals[i] = std::abs( poly_eval( hardCase_polynomial, psolvef.roots()[i] ) ); } + cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl; + + cout << "Using double's almost always solves the problem for small degrees: " << endl; + cout << "-------------------------------------------------------------------" << endl; + PolynomialSolver psolve6d( hardCase_polynomial.cast() ); + cout << "Complex roots: " << psolve6d.roots().transpose() << endl; + for( int i=0; i<6; ++i ) + { + std::complex castedRoot( psolve6d.roots()[i].real(), psolve6d.roots()[i].imag() ); + evals[i] = std::abs( poly_eval( hardCase_polynomial, castedRoot ) ); + } + cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl; + + cout.precision(10); + cout << "The last root in float then in double: " << psolvef.roots()[5] << "\t" << psolve6d.roots()[5] << endl; + std::complex castedRoot( psolve6d.roots()[5].real(), psolve6d.roots()[5].imag() ); + cout << "Norm of the difference: " << ei_abs( psolvef.roots()[5] - castedRoot ) << endl; +} diff --git a/unsupported/doc/examples/PolynomialUtils1.cpp b/unsupported/doc/examples/PolynomialUtils1.cpp new file mode 100644 index 000000000..dbfe520b5 --- /dev/null +++ b/unsupported/doc/examples/PolynomialUtils1.cpp @@ -0,0 +1,20 @@ +#include +#include + +using namespace Eigen; +using namespace std; + +int main() +{ + Vector4d roots = Vector4d::Random(); + cout << "Roots: " << roots.transpose() << endl; + Eigen::Matrix polynomial; + roots_to_monicPolynomial( roots, polynomial ); + cout << "Polynomial: "; + for( int i=0; i<4; ++i ){ cout << polynomial[i] << ".x^" << i << "+ "; } + cout << polynomial[4] << ".x^4" << endl; + Vector4d evaluation; + for( int i=0; i<4; ++i ){ + evaluation[i] = poly_eval( polynomial, roots[i] ); } + cout << "Evaluation of the polynomial at the roots: " << evaluation.transpose(); +} diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index 9618a5a56..f01d99c36 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -24,3 +24,19 @@ if(FFTW_FOUND) ei_add_test(FFTW "-DEIGEN_FFTW_DEFAULT " "-lfftw3 -lfftw3f -lfftw3l" ) endif(FFTW_FOUND) +find_package(GSL) +if(GSL_FOUND AND GSL_VERSION_MINOR LESS 9) + set(GSL_FOUND "") +endif(GSL_FOUND AND GSL_VERSION_MINOR LESS 9) +if(GSL_FOUND) + add_definitions("-DHAS_GSL" ${GSL_DEFINITIONS}) + include_directories(${GSL_INCLUDE_DIR}) + ei_add_property(EIGEN_TESTED_BACKENDS "GSL, ") +else(GSL_FOUND) + ei_add_property(EIGEN_MISSING_BACKENDS "GSL, ") + set(GSL_LIBRARIES " ") +endif(GSL_FOUND) + +ei_add_test(polynomialutils) +ei_add_test(polynomialsolver " " "${GSL_LIBRARIES}" ) + From 5ef83d6a6cd9cc4a471cfc3cc45518a813f47b92 Mon Sep 17 00:00:00 2001 From: Manuel Yguel Date: Thu, 25 Mar 2010 03:21:52 +0100 Subject: [PATCH 02/19] Add missing test files for Polynomials module. --- unsupported/test/polynomialsolver.cpp | 263 ++++++++++++++++++++++++++ unsupported/test/polynomialutils.cpp | 124 ++++++++++++ 2 files changed, 387 insertions(+) create mode 100644 unsupported/test/polynomialsolver.cpp create mode 100644 unsupported/test/polynomialutils.cpp diff --git a/unsupported/test/polynomialsolver.cpp b/unsupported/test/polynomialsolver.cpp new file mode 100644 index 000000000..1f2d7e1f3 --- /dev/null +++ b/unsupported/test/polynomialsolver.cpp @@ -0,0 +1,263 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Manuel Yguel +// +// 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 . + +#include "main.h" +#include +#include +#include + +#ifdef HAS_GSL +#include "gsl_helper.h" +#endif + +using namespace std; + +template +struct ei_increment_if_fixed_size +{ + enum { + ret = (Size == Dynamic) ? Dynamic : Size+1 + }; +}; + + + + +template +bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve ) +{ + typedef typename POLYNOMIAL::Scalar Scalar; + + typedef typename SOLVER::RootsType RootsType; + typedef Matrix EvalRootsType; + + const int deg = pols.size()-1; + + psolve.compute( pols ); + const RootsType& roots( psolve.roots() ); + EvalRootsType evr( deg ); + for( int i=0; i() ); + if( !evalToZero ) + { + cerr << "WRONG root: " << endl; + cerr << "Polynomial: " << pols.transpose() << endl; + cerr << "Roots found: " << roots.transpose() << endl; + cerr << "Abs value of the polynomial at the roots: " << evr.transpose() << endl; + cerr << endl; + } + + #ifdef HAS_GSL + if (ei_is_same_type< Scalar, double>::ret) + { + typedef GslTraits Gsl; + RootsType gslRoots(deg); + Gsl::eigen_poly_solve( pols, gslRoots ); + EvalRootsType gslEvr( deg ); + for( int i=0; i() ); + if( !evalToZero ) + { + if( !gslEvalToZero ){ + cerr << "GSL also failed" << endl; } + else{ + cerr << "GSL did NOT failed" << endl; } + cerr << "GSL roots found: " << gslRoots.transpose() << endl; + cerr << "Abs value of the polynomial at the GSL roots: " << gslEvr.transpose() << endl; + cerr << endl; + } + } + #endif //< HAS_GSL + + + std::vector rootModuli( roots.size() ); + Map< EvalRootsType > aux( &rootModuli[0], roots.size() ); + aux = roots.array().abs(); + std::sort( rootModuli.begin(), rootModuli.end() ); + bool distinctModuli=true; + for( size_t i=1; i +void evalSolver( const POLYNOMIAL& pols ) +{ + typedef typename POLYNOMIAL::Scalar Scalar; + + typedef PolynomialSolver PolynomialSolverType; + + PolynomialSolverType psolve; + aux_evalSolver( pols, psolve ); +} + + + + +template< int Deg, typename POLYNOMIAL, typename ROOTS, typename REAL_ROOTS > +void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const REAL_ROOTS& real_roots ) +{ + typedef typename POLYNOMIAL::Scalar Scalar; + + typedef PolynomialSolver PolynomialSolverType; + + PolynomialSolverType psolve; + if( aux_evalSolver( pols, psolve ) ) + { + //It is supposed that + // 1) the roots found are correct + // 2) the roots have distinct moduli + + typedef typename POLYNOMIAL::Scalar Scalar; + typedef typename REAL_ROOTS::Scalar Real; + + typedef PolynomialSolver PolynomialSolverType; + typedef typename PolynomialSolverType::RootsType RootsType; + typedef Matrix EvalRootsType; + + //Test realRoots + std::vector< Real > calc_realRoots; + psolve.realRoots( calc_realRoots ); + VERIFY( calc_realRoots.size() == (size_t)real_roots.size() ); + + const Scalar psPrec = ei_sqrt( test_precision() ); + + for( size_t i=0; i 0 ) ); + if( hasRealRoot ){ + VERIFY( ei_isApprox( real_roots.array().abs().maxCoeff(), ei_abs(r), psPrec ) ); } + + //Test absSmallestRealRoot + r = psolve.absSmallestRealRoot( hasRealRoot ); + VERIFY( hasRealRoot == (real_roots.size() > 0 ) ); + if( hasRealRoot ){ + VERIFY( ei_isApprox( real_roots.array().abs().minCoeff(), ei_abs( r ), psPrec ) ); } + + //Test greatestRealRoot + r = psolve.greatestRealRoot( hasRealRoot ); + VERIFY( hasRealRoot == (real_roots.size() > 0 ) ); + if( hasRealRoot ){ + VERIFY( ei_isApprox( real_roots.array().maxCoeff(), r, psPrec ) ); } + + //Test smallestRealRoot + r = psolve.smallestRealRoot( hasRealRoot ); + VERIFY( hasRealRoot == (real_roots.size() > 0 ) ); + if( hasRealRoot ){ + VERIFY( ei_isApprox( real_roots.array().minCoeff(), r, psPrec ) ); } + } +} + + +template +void polynomialsolver(int deg) +{ + typedef ei_increment_if_fixed_size<_Deg> Dim; + typedef Matrix<_Scalar,Dim::ret,1> PolynomialType; + typedef Matrix<_Scalar,_Deg,1> EvalRootsType; + + cout << "Standard cases" << endl; + PolynomialType pols = PolynomialType::Random(deg+1); + evalSolver<_Deg,PolynomialType>( pols ); + + cout << "Hard cases" << endl; + _Scalar multipleRoot = ei_random<_Scalar>(); + EvalRootsType allRoots = EvalRootsType::Constant(deg,multipleRoot); + roots_to_monicPolynomial( allRoots, pols ); + evalSolver<_Deg,PolynomialType>( pols ); + + cout << "Test sugar" << endl; + EvalRootsType realRoots = EvalRootsType::Random(deg); + roots_to_monicPolynomial( realRoots, pols ); + evalSolverSugarFunction<_Deg>( + pols, + realRoots.template cast < + std::complex< + typename NumTraits<_Scalar>::Real + > + >(), + realRoots ); +} + + +template void polynomialsolver_scalar() +{ + CALL_SUBTEST_1( (polynomialsolver<_Scalar,1>(1)) ); + CALL_SUBTEST_2( (polynomialsolver<_Scalar,2>(2)) ); + CALL_SUBTEST_3( (polynomialsolver<_Scalar,3>(3)) ); + CALL_SUBTEST_4( (polynomialsolver<_Scalar,4>(4)) ); + CALL_SUBTEST_5( (polynomialsolver<_Scalar,5>(5)) ); + CALL_SUBTEST_6( (polynomialsolver<_Scalar,6>(6)) ); + CALL_SUBTEST_7( (polynomialsolver<_Scalar,7>(7)) ); + CALL_SUBTEST_8( (polynomialsolver<_Scalar,8>(8)) ); + + CALL_SUBTEST_9( (polynomialsolver<_Scalar,Dynamic>( + ei_random(9,45) + )) ); +} + +void test_polynomialsolver() +{ + for(int i = 0; i < g_repeat; i++) + { + polynomialsolver_scalar(); + polynomialsolver_scalar(); + } +} diff --git a/unsupported/test/polynomialutils.cpp b/unsupported/test/polynomialutils.cpp new file mode 100644 index 000000000..b1afaccc6 --- /dev/null +++ b/unsupported/test/polynomialutils.cpp @@ -0,0 +1,124 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Manuel Yguel +// +// 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 . + +#include "main.h" +#include +#include + +using namespace std; + +template +struct ei_increment_if_fixed_size +{ + enum { + ret = (Size == Dynamic) ? Dynamic : Size+1 + }; +}; + +template +void realRoots_to_monicPolynomial_test(int deg) +{ + typedef ei_increment_if_fixed_size<_Deg> Dim; + typedef Matrix<_Scalar,Dim::ret,1> PolynomialType; + typedef Matrix<_Scalar,_Deg,1> EvalRootsType; + + PolynomialType pols(deg+1); + EvalRootsType roots = EvalRootsType::Random(deg); + realRoots_to_monicPolynomial( roots, pols ); + + EvalRootsType evr( deg ); + for( int i=0; i() ); + if( !evalToZero ){ + cerr << evr.transpose() << endl; } + VERIFY( evalToZero ); +} + +template void realRoots_to_monicPolynomial_scalar() +{ + CALL_SUBTEST_2( (realRoots_to_monicPolynomial_test<_Scalar,2>(2)) ); + CALL_SUBTEST_3( (realRoots_to_monicPolynomial_test<_Scalar,3>(3)) ); + CALL_SUBTEST_4( (realRoots_to_monicPolynomial_test<_Scalar,4>(4)) ); + CALL_SUBTEST_5( (realRoots_to_monicPolynomial_test<_Scalar,5>(5)) ); + CALL_SUBTEST_6( (realRoots_to_monicPolynomial_test<_Scalar,6>(6)) ); + CALL_SUBTEST_7( (realRoots_to_monicPolynomial_test<_Scalar,7>(7)) ); + CALL_SUBTEST_8( (realRoots_to_monicPolynomial_test<_Scalar,17>(17)) ); + + int deg = ei_random(18,26); + CALL_SUBTEST_9( (realRoots_to_monicPolynomial_test<_Scalar,Dynamic>(deg)) ); +} + + + + +template +void CauchyBounds(int deg) +{ + typedef ei_increment_if_fixed_size<_Deg> Dim; + typedef Matrix<_Scalar,Dim::ret,1> PolynomialType; + typedef Matrix<_Scalar,_Deg,1> EvalRootsType; + + PolynomialType pols(deg+1); + EvalRootsType roots = EvalRootsType::Random(deg); + realRoots_to_monicPolynomial( roots, pols ); + _Scalar M = cauchy_max_bound( pols ); + _Scalar m = cauchy_min_bound( pols ); + _Scalar Max = roots.array().abs().maxCoeff(); + _Scalar min = roots.array().abs().minCoeff(); + bool eval = (M >= Max) && (m <= min); + if( !eval ) + { + cerr << "Roots: " << roots << endl; + cerr << "Bounds: (" << m << ", " << M << ")" << endl; + cerr << "Min,Max: (" << min << ", " << Max << ")" << endl; + } + VERIFY( eval ); +} + +template void CauchyBounds_scalar() +{ + CALL_SUBTEST_2( (CauchyBounds<_Scalar,2>(2)) ); + CALL_SUBTEST_3( (CauchyBounds<_Scalar,3>(3)) ); + CALL_SUBTEST_4( (CauchyBounds<_Scalar,4>(4)) ); + CALL_SUBTEST_5( (CauchyBounds<_Scalar,5>(5)) ); + CALL_SUBTEST_6( (CauchyBounds<_Scalar,6>(6)) ); + CALL_SUBTEST_7( (CauchyBounds<_Scalar,7>(7)) ); + CALL_SUBTEST_8( (CauchyBounds<_Scalar,17>(17)) ); + + int deg = ei_random(18,26); + CALL_SUBTEST_9( (CauchyBounds<_Scalar,Dynamic>(deg)) ); +} + +void test_polynomialutils() +{ + for(int i = 0; i < g_repeat; i++) + { + realRoots_to_monicPolynomial_scalar(); + realRoots_to_monicPolynomial_scalar(); + CauchyBounds_scalar(); + CauchyBounds_scalar(); + } +} From 89f2d5667fd1ace51ab3e05ad614adbbcff8dab8 Mon Sep 17 00:00:00 2001 From: Manuel Yguel Date: Thu, 25 Mar 2010 03:25:47 +0100 Subject: [PATCH 03/19] Add the possibility to use the polynomial solver of the gsl. --- test/gsl_helper.h | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/test/gsl_helper.h b/test/gsl_helper.h index f8f7a5d79..eefffe792 100644 --- a/test/gsl_helper.h +++ b/test/gsl_helper.h @@ -33,6 +33,7 @@ #include #include #include +#include namespace Eigen { @@ -69,6 +70,27 @@ template::IsComplex> struct gsl_eigen_gensymmv_free(w); free(a); } + + template + static void eigen_poly_solve(const EIGEN_VECTOR& poly, EIGEN_ROOTS& roots ) + { + const int deg = poly.size()-1; + double *z = new double[2*deg]; + double *a = new double[poly.size()]; + for( int i=0; i struct GslTraits From 671cfb7ad0776ae42cedcb3361b9f2238d7c1ffe Mon Sep 17 00:00:00 2001 From: Manuel Yguel Date: Thu, 25 Mar 2010 03:31:33 +0100 Subject: [PATCH 04/19] Fix wrong header and warnings in polynomialutils.cpp --- unsupported/test/polynomialutils.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/unsupported/test/polynomialutils.cpp b/unsupported/test/polynomialutils.cpp index b1afaccc6..7f93c2f0d 100644 --- a/unsupported/test/polynomialutils.cpp +++ b/unsupported/test/polynomialutils.cpp @@ -23,7 +23,7 @@ // Eigen. If not, see . #include "main.h" -#include +#include #include using namespace std; @@ -45,7 +45,7 @@ void realRoots_to_monicPolynomial_test(int deg) PolynomialType pols(deg+1); EvalRootsType roots = EvalRootsType::Random(deg); - realRoots_to_monicPolynomial( roots, pols ); + roots_to_monicPolynomial( roots, pols ); EvalRootsType evr( deg ); for( int i=0; i void realRoots_to_monicPolynomial_scalar() CALL_SUBTEST_7( (realRoots_to_monicPolynomial_test<_Scalar,7>(7)) ); CALL_SUBTEST_8( (realRoots_to_monicPolynomial_test<_Scalar,17>(17)) ); - int deg = ei_random(18,26); - CALL_SUBTEST_9( (realRoots_to_monicPolynomial_test<_Scalar,Dynamic>(deg)) ); + CALL_SUBTEST_9( (realRoots_to_monicPolynomial_test<_Scalar,Dynamic>( + ei_random(18,26) )) ); } @@ -83,7 +83,7 @@ void CauchyBounds(int deg) PolynomialType pols(deg+1); EvalRootsType roots = EvalRootsType::Random(deg); - realRoots_to_monicPolynomial( roots, pols ); + roots_to_monicPolynomial( roots, pols ); _Scalar M = cauchy_max_bound( pols ); _Scalar m = cauchy_min_bound( pols ); _Scalar Max = roots.array().abs().maxCoeff(); @@ -108,8 +108,8 @@ template void CauchyBounds_scalar() CALL_SUBTEST_7( (CauchyBounds<_Scalar,7>(7)) ); CALL_SUBTEST_8( (CauchyBounds<_Scalar,17>(17)) ); - int deg = ei_random(18,26); - CALL_SUBTEST_9( (CauchyBounds<_Scalar,Dynamic>(deg)) ); + CALL_SUBTEST_9( (CauchyBounds<_Scalar,Dynamic>( + ei_random(18,26) )) ); } void test_polynomialutils() From df1cbe8c3daeafa413507cc9885b321afb217d71 Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Thu, 25 Mar 2010 08:50:23 +0100 Subject: [PATCH 05/19] fix display of modules list in documentation --- doc/unsupported_modules.dox | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/unsupported_modules.dox b/doc/unsupported_modules.dox index 503ecb8b3..2f2b2f6a8 100644 --- a/doc/unsupported_modules.dox +++ b/doc/unsupported_modules.dox @@ -43,6 +43,9 @@ namespace Eigen { /** \ingroup Unsupported_modules * \defgroup NumericalDiff_Module */ +/** \ingroup Unsupported_modules + * \defgroup Polynomials_Module */ + /** \ingroup Unsupported_modules * \defgroup Skyline_Module */ From eb0c921a588203323390da4ba08b0c4f454efc8d Mon Sep 17 00:00:00 2001 From: Manuel Yguel Date: Thu, 25 Mar 2010 16:54:00 +0100 Subject: [PATCH 06/19] Fix some doc typos. --- unsupported/Eigen/Polynomials | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/unsupported/Eigen/Polynomials b/unsupported/Eigen/Polynomials index 9e1f6b759..d69abb1be 100644 --- a/unsupported/Eigen/Polynomials +++ b/unsupported/Eigen/Polynomials @@ -64,7 +64,8 @@ namespace Eigen { \endcode evaluates a polynomial at a given point using stabilized Hörner method. - The following code computes the coefficients in the monomial basis of the monic polynomial given through its roots then evaluate it at those roots. + The following code: first computes the coefficients in the monomial basis of the monic polynomial that has the provided roots; + then, it evaluates the computed polynomial, using a stabilized Hörner method. \include PolynomialUtils1.cpp Output: \verbinclude PolynomialUtils1.out @@ -116,11 +117,15 @@ namespace Eigen { (double) floating types and small polynomial degree (<20). \include PolynomialSolver1.cpp - In the example a polynomial with almost conjugate roots is provided to the solver. + + In the above example: + + -# a simple use of the polynomial solver is shown; + -# the accuracy problem with the QR algorithm is presented: a polynomial with almost conjugate roots is provided to the solver. Those roots have almost same module therefore the QR algorithm failed to converge: the accuracy - of the last root is bad. + of the last root is bad; + -# a simple way to circumvent the problem is shown: use doubles instead of floats. - This problem is less visible with double. Output: \verbinclude PolynomialSolver1.out */ From af08770890acdc6f78e9a4730c543b6bd1aadd3a Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Fri, 26 Mar 2010 12:36:42 +0000 Subject: [PATCH 07/19] Center version number on main page of API documentation. --- doc/eigendoxy.css | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/doc/eigendoxy.css b/doc/eigendoxy.css index 77d062837..c1973289b 100644 --- a/doc/eigendoxy.css +++ b/doc/eigendoxy.css @@ -477,3 +477,8 @@ DIV.eimainmenu { /* border-top: solid; */ /* border-bottom: solid; */ } + +/* center version number on main page */ +H3.version { + text-align: center; +} From 0a5c2d8a54bf0bdab7a7c68e824002ba163bbdca Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Sat, 27 Mar 2010 18:58:29 +0100 Subject: [PATCH 08/19] fix misc warnings, more importantly when NDEBUG is defined, assert() is a nop. --- Eigen/src/Cholesky/LLT.h | 3 +-- Eigen/src/Core/DenseBase.h | 2 +- Eigen/src/LU/PartialPivLU.h | 3 +-- unsupported/Eigen/FFT | 2 +- 4 files changed, 4 insertions(+), 6 deletions(-) diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index d552e4e8a..e95b97016 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -296,8 +296,7 @@ template bool LLT::solveInPlace(MatrixBase &bAndX) const { ei_assert(m_isInitialized && "LLT is not initialized."); - const int size = m_matrix.rows(); - ei_assert(size==bAndX.rows()); + ei_assert(m_matrix.rows()==bAndX.rows()); matrixL().solveInPlace(bAndX); matrixU().solveInPlace(bAndX); return true; diff --git a/Eigen/src/Core/DenseBase.h b/Eigen/src/Core/DenseBase.h index 29a95b8db..199d4aa98 100644 --- a/Eigen/src/Core/DenseBase.h +++ b/Eigen/src/Core/DenseBase.h @@ -584,7 +584,7 @@ template class DenseBase #endif // disable the use of evalTo for dense objects with a nice compilation error - template inline void evalTo(Dest& dst) const + template inline void evalTo(Dest& ) const { EIGEN_STATIC_ASSERT((ei_is_same_type::ret),THE_EVAL_EVALTO_FUNCTION_SHOULD_NEVER_BE_CALLED_FOR_DENSE_OBJECTS); } diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 8ff9c5b03..61d59b8cb 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -439,8 +439,7 @@ struct ei_solve_retval, Rhs> * Step 3: replace c by the solution x to Ux = c. */ - const int size = dec().matrixLU().rows(); - ei_assert(rhs().rows() == size); + ei_assert(rhs().rows() == dec().matrixLU().rows()); // Step 1 dst = dec().permutationP() * rhs(); diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT index 83f761028..23e0275b2 100644 --- a/unsupported/Eigen/FFT +++ b/unsupported/Eigen/FFT @@ -320,7 +320,7 @@ class FFT // if the vector is strided, then we need to copy it to a packed temporary Matrix tmp; if ( resize_input ) { - size_t ncopy = min(src.size(),src.size() + resize_input); + size_t ncopy = std::min(src.size(),src.size() + resize_input); tmp.setZero(src.size() + resize_input); if ( realfft && HasFlag(HalfSpectrum) ) { // pad at the Nyquist bin From e6300efb5c97cbd66b58b944441f66147bb375ad Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Sun, 28 Mar 2010 17:33:56 +0100 Subject: [PATCH 09/19] Extend documentation for HessenbergDecomposition. --- .../src/Eigenvalues/HessenbergDecomposition.h | 165 ++++++++++++++---- .../HessenbergDecomposition_compute.cpp | 6 + .../HessenbergDecomposition_matrixH.cpp | 8 + .../HessenbergDecomposition_packedMatrix.cpp | 9 + 4 files changed, 158 insertions(+), 30 deletions(-) create mode 100644 doc/snippets/HessenbergDecomposition_compute.cpp create mode 100644 doc/snippets/HessenbergDecomposition_matrixH.cpp create mode 100644 doc/snippets/HessenbergDecomposition_packedMatrix.cpp diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h index f87c0b842..1b2dcbe58 100644 --- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h +++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -30,14 +30,30 @@ * * \class HessenbergDecomposition * - * \brief Reduces a squared matrix to an Hessemberg form + * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation * - * \param MatrixType the type of the matrix of which we are computing the Hessenberg decomposition + * \tparam _MatrixType the type of the matrix of which we are computing the Hessenberg decomposition * - * This class performs an Hessenberg decomposition of a matrix \f$ A \f$ such that: - * \f$ A = Q H Q^* \f$ where \f$ Q \f$ is unitary and \f$ H \f$ a Hessenberg matrix. + * This class performs an Hessenberg decomposition of a matrix \f$ A \f$. In + * the real case, the Hessenberg decomposition consists of an orthogonal + * matrix \f$ Q \f$ and a Hessenberg matrix \f$ H \f$ such that \f$ A = Q H + * Q^T \f$. An orthogonal matrix is a matrix whose inverse equals its + * transpose (\f$ Q^{-1} = Q^T \f$). A Hessenberg matrix has zeros below the + * subdiagonal, so it is almost upper triangular. The Hessenberg decomposition + * of a complex matrix is \f$ A = Q H Q^* \f$ with \f$ Q \f$ unitary (that is, + * \f$ Q^{-1} = Q^* \f$). * - * \sa class Tridiagonalization, class Qr + * Call the function compute() to compute the Hessenberg decomposition of a + * given matrix. Alternatively, you can use the + * HessenbergDecomposition(const MatrixType&) constructor which computes the + * Hessenberg decomposition at construction time. Once the decomposition is + * computed, you can use the matrixH() and matrixQ() functions to construct + * the matrices H and Q in the decomposition. + * + * The documentation for matrixH() contains an example of the typical use of + * this class. + * + * \sa class ComplexSchur, class Tridiagonalization, \ref QR_Module "QR Module" */ template class HessenbergDecomposition { @@ -51,13 +67,28 @@ template class HessenbergDecomposition MaxSize = MatrixType::MaxRowsAtCompileTime, MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1 }; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef Matrix CoeffVectorType; - typedef Matrix VectorType; - /** This constructor initializes a HessenbergDecomposition object for - * further use with HessenbergDecomposition::compute() + /** \brief Scalar type for matrices of type \p _MatrixType. */ + typedef typename MatrixType::Scalar Scalar; + + /** \brief Type for vector of Householder coefficients. + * + * This is column vector with entries of type #Scalar. The length of the + * vector is one less than the size of \p _MatrixType, if it is a + * fixed-side type. + */ + typedef Matrix CoeffVectorType; + + /** \brief Default constructor; the decomposition will be computed later. + * + * \param [in] size The size of the matrix whose Hessenberg decomposition will be computed. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via compute(). 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() for an example. */ HessenbergDecomposition(int size = Size==Dynamic ? 2 : Size) : m_matrix(size,size) @@ -66,6 +97,15 @@ template class HessenbergDecomposition m_hCoeffs.resize(size-1); } + /** \brief Constructor; computes Hessenberg decomposition of given matrix. + * + * \param[in] matrix Square matrix whose Hessenberg decomposition is to be computed. + * + * This constructor calls compute() to compute the Hessenberg + * decomposition. + * + * \sa matrixH() for an example. + */ HessenbergDecomposition(const MatrixType& matrix) : m_matrix(matrix) { @@ -75,9 +115,21 @@ template class HessenbergDecomposition _compute(m_matrix, m_hCoeffs); } - /** Computes or re-compute the Hessenberg decomposition for the matrix \a matrix. + /** \brief Computes Hessenberg decomposition of given matrix. + * + * \param[in] matrix Square matrix whose Hessenberg decomposition is to be computed. * - * This method allows to re-use the allocated data. + * The Hessenberg decomposition is computed by bringing the columns of the + * matrix successively in the required form using Householder reflections + * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, %Matrix + * Computations). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$ + * denotes the size of the given matrix. + * + * This method reuses of the allocated data in the HessenbergDecomposition + * object. + * + * Example: \include HessenbergDecomposition_compute.cpp + * Output: \verbinclude HessenbergDecomposition_compute.out */ void compute(const MatrixType& matrix) { @@ -88,36 +140,95 @@ template class HessenbergDecomposition _compute(m_matrix, m_hCoeffs); } - /** \returns a const reference to the householder coefficients allowing to - * reconstruct the matrix Q from the packed data. + /** \brief Returns the Householder coefficients. * - * \sa packedMatrix() + * \returns a const reference to the vector of Householder coefficients + * + * \pre Either the constructor HessenbergDecomposition(const MatrixType&) + * or the member function compute(const MatrixType&) has been called + * before to compute the Hessenberg decomposition of a matrix. + * + * The Householder coefficients allow the reconstruction of the matrix + * \f$ Q \f$ in the Hessenberg decomposition from the packed data. + * + * \sa packedMatrix(), \ref Householder_Module "Householder module" */ const CoeffVectorType& householderCoefficients() const { return m_hCoeffs; } - /** \returns a const reference to the internal representation of the decomposition. + /** \brief Returns the internal representation of the decomposition + * + * \returns a const reference to a matrix with the internal representation + * of the decomposition. + * + * \pre Either the constructor HessenbergDecomposition(const MatrixType&) + * or the member function compute(const MatrixType&) has been called + * before to compute the Hessenberg decomposition of a matrix. * * The returned matrix contains the following information: * - the upper part and lower sub-diagonal represent the Hessenberg matrix H * - the rest of the lower part contains the Householder vectors that, combined with * Householder coefficients returned by householderCoefficients(), - * allows to reconstruct the matrix Q as follow: - * Q = H_{N-1} ... H_1 H_0 - * where the matrices H are the Householder transformation: - * H_i = (I - h_i * v_i * v_i') - * where h_i == householderCoefficients()[i] and v_i is a Householder vector: - * v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ] + * allows to reconstruct the matrix Q as + * \f$ Q = H_{N-1} \ldots H_1 H_0 \f$. + * Here, the matrices \f$ H_i \f$ are the Householder transformations + * \f$ H_i = (I - h_i v_i v_i^T) \f$ + * where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and + * \f$ v_i \f$ is the Householder vector defined by + * \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$ + * with M the matrix returned by this function. * * See LAPACK for further details on this packed storage. + * + * Example: \include HessenbergDecomposition_packedMatrix.cpp + * Output: \verbinclude HessenbergDecomposition_packedMatrix.out + * + * \sa householderCoefficients() */ const MatrixType& packedMatrix(void) const { return m_matrix; } + /** \brief Reconstructs the orthogonal matrix Q in the decomposition + * + * \returns the matrix Q + * + * \pre Either the constructor HessenbergDecomposition(const MatrixType&) + * or the member function compute(const MatrixType&) has been called + * before to compute the Hessenberg decomposition of a matrix. + * + * This function reconstructs the matrix Q from the Householder + * coefficients and the packed matrix stored internally. This + * reconstruction requires \f$ 4n^3 / 3 \f$ flops. + * + * \sa matrixH() for an example + */ MatrixType matrixQ() const; + + /** \brief Constructs the Hessenberg matrix H in the decomposition + * + * \returns the matrix H + * + * \pre Either the constructor HessenbergDecomposition(const MatrixType&) + * or the member function compute(const MatrixType&) has been called + * before to compute the Hessenberg decomposition of a matrix. + * + * This function copies the matrix H from internal data. The upper part + * (including the subdiagonal) of the packed matrix as returned by + * packedMatrix() contains the matrix H. This function copies those + * entries in a newly created matrix and sets the remaining entries to + * zero. It may sometimes be sufficient to directly use the packed matrix + * instead of creating a new one. + * + * Example: \include HessenbergDecomposition_matrixH.cpp + * Output: \verbinclude HessenbergDecomposition_matrixH.out + * + * \sa matrixQ(), packedMatrix() + */ MatrixType matrixH() const; private: static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs); + typedef Matrix VectorType; + typedef typename NumTraits::Real RealScalar; protected: MatrixType m_matrix; @@ -134,7 +245,7 @@ template class HessenbergDecomposition * * The result is written in the lower triangular part of \a matA. * - * Implemented from Golub's "Matrix Computations", algorithm 8.3.1. + * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1. * * \sa packedMatrix() */ @@ -167,7 +278,6 @@ void HessenbergDecomposition::_compute(MatrixType& matA, CoeffVector } } -/** reconstructs and returns the matrix Q */ template typename HessenbergDecomposition::MatrixType HessenbergDecomposition::matrixQ() const @@ -185,11 +295,6 @@ HessenbergDecomposition::matrixQ() const #endif // EIGEN_HIDE_HEAVY_CODE -/** constructs and returns the matrix H. - * Note that the matrix H is equivalent to the upper part of the packed matrix - * (including the lower sub-diagonal). Therefore, it might be often sufficient - * to directly use the packed matrix instead of creating a new one. - */ template typename HessenbergDecomposition::MatrixType HessenbergDecomposition::matrixH() const diff --git a/doc/snippets/HessenbergDecomposition_compute.cpp b/doc/snippets/HessenbergDecomposition_compute.cpp new file mode 100644 index 000000000..50e37833a --- /dev/null +++ b/doc/snippets/HessenbergDecomposition_compute.cpp @@ -0,0 +1,6 @@ +MatrixXcf A = MatrixXcf::Random(4,4); +HessenbergDecomposition hd(4); +hd.compute(A); +cout << "The matrix H in the decomposition of A is:" << endl << hd.matrixH() << endl; +hd.compute(2*A); // re-use hd to compute and store decomposition of 2A +cout << "The matrix H in the decomposition of 2A is:" << endl << hd.matrixH() << endl; diff --git a/doc/snippets/HessenbergDecomposition_matrixH.cpp b/doc/snippets/HessenbergDecomposition_matrixH.cpp new file mode 100644 index 000000000..af0136668 --- /dev/null +++ b/doc/snippets/HessenbergDecomposition_matrixH.cpp @@ -0,0 +1,8 @@ +Matrix4f A = MatrixXf::Random(4,4); +cout << "Here is a random 4x4 matrix:" << endl << A << endl; +HessenbergDecomposition hessOfA(A); +MatrixXf H = hessOfA.matrixH(); +cout << "The Hessenberg matrix H is:" << endl << H << endl; +MatrixXf Q = hessOfA.matrixQ(); +cout << "The orthogonal matrix Q is:" << endl << Q << endl; +cout << "Q H Q^T is:" << endl << Q * H * Q.transpose() << endl; diff --git a/doc/snippets/HessenbergDecomposition_packedMatrix.cpp b/doc/snippets/HessenbergDecomposition_packedMatrix.cpp new file mode 100644 index 000000000..4fa5957e8 --- /dev/null +++ b/doc/snippets/HessenbergDecomposition_packedMatrix.cpp @@ -0,0 +1,9 @@ +Matrix4d A = Matrix4d::Random(4,4); +cout << "Here is a random 4x4 matrix:" << endl << A << endl; +HessenbergDecomposition hessOfA(A); +Matrix4d pm = hessOfA.packedMatrix(); +cout << "The packed matrix M is:" << endl << pm << endl; +cout << "The upper Hessenberg part corresponds to the matrix H, which is:" + << endl << hessOfA.matrixH() << endl; +Vector3d hc = hessOfA.householderCoefficients(); +cout << "The vector of Householder coefficients is:" << endl << hc << endl; From 8f99ae5ea44965396a6174ab97a41f0a2101b984 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Tue, 30 Mar 2010 11:38:09 -0400 Subject: [PATCH 10/19] move the computation of the number of significant digits to a templated helper struct, that can be specialized to custom types if needed. Should address this request: http://forum.kde.org/viewtopic.php?f=74&t=86914 --- Eigen/src/Core/IO.h | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/Eigen/src/Core/IO.h b/Eigen/src/Core/IO.h index 3e8d2bc66..ffb214894 100644 --- a/Eigen/src/Core/IO.h +++ b/Eigen/src/Core/IO.h @@ -126,6 +126,16 @@ DenseBase::format(const IOFormat& fmt) const return WithFormat(derived(), fmt); } +template +struct ei_significant_decimals_impl +{ + typedef typename NumTraits::Real RealScalar; + static inline int run() + { + return (int) std::ceil(-ei_log(NumTraits::epsilon())/ei_log(RealScalar(10))); + } +}; + /** \internal * print the matrix \a _m to the output stream \a s using the output format \a fmt */ template @@ -145,9 +155,7 @@ std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOForm { if (NumTraits::HasFloatingPoint) { - typedef typename NumTraits::Real RealScalar; - RealScalar explicit_precision_fp = std::ceil(-ei_log(NumTraits::epsilon())/ei_log(10.0)); - explicit_precision = static_cast(explicit_precision_fp); + explicit_precision = ei_significant_decimals_impl::run(); } else { From 9e0d8697c78fb4668246c59186eca342ba3eff88 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Tue, 30 Mar 2010 14:16:54 -0400 Subject: [PATCH 11/19] add ei_cast_to_int, we are indeed somethings (e.g. in IO.h) casting scalars to int and the only way to allow users to extend that to their own scalar types that don't have int cast operators, was to allow them specialize ei_cast_to_int_impl. --- Eigen/src/Core/IO.h | 2 +- Eigen/src/Core/MathFunctions.h | 14 ++++++++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Core/IO.h b/Eigen/src/Core/IO.h index ffb214894..22db103ed 100644 --- a/Eigen/src/Core/IO.h +++ b/Eigen/src/Core/IO.h @@ -132,7 +132,7 @@ struct ei_significant_decimals_impl typedef typename NumTraits::Real RealScalar; static inline int run() { - return (int) std::ceil(-ei_log(NumTraits::epsilon())/ei_log(RealScalar(10))); + return ei_cast_to_int(std::ceil(-ei_log(NumTraits::epsilon())/ei_log(RealScalar(10)))); } }; diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index c97a68e50..943b44cfa 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -44,6 +44,20 @@ template inline typename NumTraits::Real ei_hypot(T x, T y) return p * ei_sqrt(T(1) + qp*qp); } +template struct ei_cast_to_int_impl +{ + static int run(const T& x) + { + return int(x); + } +}; + +template inline int ei_cast_to_int(const T& x) +{ + return ei_cast_to_int_impl::run(x); +} + + /************** *** int *** **************/ From 16e416b8d745b1fc64111581fec431d843591274 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Tue, 30 Mar 2010 14:47:45 -0400 Subject: [PATCH 12/19] generalize the idea of the previous commit to all kinds of casts, see this forum thread: http://forum.kde.org/viewtopic.php?f=74&t=86914 this is important to allow users to support custom types that don't have the needed conversion operators. --- Eigen/src/Core/IO.h | 2 +- Eigen/src/Core/MathFunctions.h | 13 ++++++++----- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/Eigen/src/Core/IO.h b/Eigen/src/Core/IO.h index 22db103ed..c98742246 100644 --- a/Eigen/src/Core/IO.h +++ b/Eigen/src/Core/IO.h @@ -132,7 +132,7 @@ struct ei_significant_decimals_impl typedef typename NumTraits::Real RealScalar; static inline int run() { - return ei_cast_to_int(std::ceil(-ei_log(NumTraits::epsilon())/ei_log(RealScalar(10)))); + return ei_cast(std::ceil(-ei_log(NumTraits::epsilon())/ei_log(RealScalar(10)))); } }; diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 943b44cfa..4a21ec975 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -44,17 +44,20 @@ template inline typename NumTraits::Real ei_hypot(T x, T y) return p * ei_sqrt(T(1) + qp*qp); } -template struct ei_cast_to_int_impl +// the point of wrapping these casts in this helper template struct is to allow users to specialize it to custom types +// that may not have the needed conversion operators (especially as c++98 doesn't have explicit conversion operators). + +template struct ei_cast_impl { - static int run(const T& x) + static inline NewType run(const OldType& x) { - return int(x); + return static_cast(x); } }; -template inline int ei_cast_to_int(const T& x) +template inline NewType ei_cast(const OldType& x) { - return ei_cast_to_int_impl::run(x); + return ei_cast_impl::run(x); } From 338ec0390fc89d9f202d2eb2a41725b7fc738a0a Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Tue, 30 Mar 2010 14:51:47 -0400 Subject: [PATCH 13/19] let the cast functor use the new ei_cast() --- Eigen/src/Core/Functors.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index c2b317cc0..d02633cb8 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -274,7 +274,7 @@ template struct ei_scalar_cast_op { EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_cast_op) typedef NewType result_type; - EIGEN_STRONG_INLINE const NewType operator() (const Scalar& a) const { return static_cast(a); } + EIGEN_STRONG_INLINE const NewType operator() (const Scalar& a) const { return ei_cast(a); } }; template struct ei_functor_traits > From 1b3f7f2feef34d0506e66e7514130a3994589be8 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Wed, 31 Mar 2010 11:59:11 +0100 Subject: [PATCH 14/19] Extend documentation and add examples for EigenSolver class. --- Eigen/src/Eigenvalues/EigenSolver.h | 216 +++++++++++++----- .../EigenSolver_EigenSolver_MatrixType.cpp | 16 ++ doc/snippets/EigenSolver_compute.cpp | 6 + doc/snippets/EigenSolver_eigenvalues.cpp | 4 + doc/snippets/EigenSolver_eigenvectors.cpp | 4 + .../EigenSolver_pseudoEigenvectors.cpp | 9 + 6 files changed, 202 insertions(+), 53 deletions(-) create mode 100644 doc/snippets/EigenSolver_EigenSolver_MatrixType.cpp create mode 100644 doc/snippets/EigenSolver_compute.cpp create mode 100644 doc/snippets/EigenSolver_eigenvalues.cpp create mode 100644 doc/snippets/EigenSolver_eigenvectors.cpp create mode 100644 doc/snippets/EigenSolver_pseudoEigenvectors.cpp diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index 579585618..84f5c408f 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -30,15 +30,44 @@ * * \class EigenSolver * - * \brief Eigen values/vectors solver for non selfadjoint matrices + * \brief Computes eigenvalues and eigenvectors of general matrices * - * \param MatrixType the type of the matrix of which we are computing the eigen decomposition + * \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. Currently, only real matrices are supported. * - * Currently it only support real matrices. + * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars + * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v \f$. If + * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and + * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V = + * V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we + * have \f$ A = V D V^{-1} \f$. This is called the eigendecomposition. + * + * The eigenvalues and eigenvectors of a matrix may be complex, even when the + * matrix is real. However, we can choose real matrices \f$ V \f$ and \f$ D + * \f$ satisfying \f$ A V = V D \f$, just like the eigendecomposition, if the + * matrix \f$ D \f$ is not required to be diagonal, but if it is allowed to + * have blocks of the form + * \f[ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f] + * (where \f$ u \f$ and \f$ v \f$ are real numbers) on the diagonal. These + * blocks correspond to complex eigenvalue pairs \f$ u \pm iv \f$. We call + * this variant of the eigendecomposition the pseudo-eigendecomposition. + * + * Call the function compute() to compute the eigenvalues and eigenvectors of + * a given matrix. Alternatively, you can use the + * EigenSolver(const MatrixType&) 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 pseudoEigenvalueMatrix() and + * pseudoEigenvectors() methods allow the construction of the + * pseudo-eigendecomposition. + * + * The documentation for EigenSolver(const MatrixType&) contains an example of + * the typical use of this class. * * \note this code was adapted from JAMA (public domain) * - * \sa MatrixBase::eigenvalues(), SelfAdjointEigenSolver + * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver */ template class EigenSolver { @@ -52,21 +81,54 @@ template class EigenSolver MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; + + /** \brief Scalar type for matrices of type \p _MatrixType. */ typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef std::complex Complex; - typedef Matrix EigenvalueType; - typedef Matrix EigenvectorType; - typedef Matrix RealVectorType; - /** - * \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via EigenSolver::compute(const MatrixType&). - */ + /** \brief Complex scalar type for \p _MatrixType. + * + * This is \c std::complex if #Scalar is real (e.g., + * \c float or \c double) and just \c Scalar if #Scalar is + * complex. + */ + typedef std::complex Complex; + + /** \brief Type for vector of eigenvalues as returned by eigenvalues(). + * + * This is a column vector with entries of type #Complex. + * The length of the vector is the size of \p _MatrixType. + */ + typedef Matrix EigenvalueType; + + /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). + * + * This is a square matrix with entries of type #Complex. + * The size is the same as the size of \p _MatrixType. + */ + typedef Matrix EigenvectorType; + + /** \brief Default constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via EigenSolver::compute(const MatrixType&). + * + * \sa compute() for an example. + */ EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {} + /** \brief Constructor; computes eigendecomposition of given matrix. + * + * \param[in] matrix Square matrix whose eigendecomposition is to be computed. + * + * This constructor calls compute() to compute the eigenvalues + * and eigenvectors. + * + * Example: \include EigenSolver_EigenSolver_MatrixType.cpp + * Output: \verbinclude EigenSolver_EigenSolver_MatrixType.out + * + * \sa compute() + */ EigenSolver(const MatrixType& matrix) : m_eivec(matrix.rows(), matrix.cols()), m_eivalues(matrix.cols()), @@ -75,39 +137,42 @@ template class EigenSolver compute(matrix); } + /** \brief Returns the eigenvectors of given matrix. + * + * \returns %Matrix whose columns are the (possibly complex) eigenvectors. + * + * \pre Either the constructor EigenSolver(const MatrixType&) or the + * member function compute(const MatrixType&) has been called before. + * + * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding + * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The + * eigenvectors are normalized to have (Euclidean) norm equal to one. The + * matrix returned by this function is the matrix \f$ V \f$ in the + * eigendecomposition \f$ A = V D V^{-1} \f$, if it exists. + * + * Example: \include EigenSolver_eigenvectors.cpp + * Output: \verbinclude EigenSolver_eigenvectors.out + * + * \sa eigenvalues(), pseudoEigenvectors() + */ + EigenvectorType eigenvectors() const; - EigenvectorType eigenvectors(void) const; - - /** \returns a real matrix V of pseudo eigenvectors. + /** \brief Returns the pseudo-eigenvectors of given matrix. * - * Let D be the block diagonal matrix with the real eigenvalues in 1x1 blocks, - * and any complex values u+iv in 2x2 blocks [u v ; -v u]. Then, the matrices D - * and V satisfy A*V = V*D. + * \returns Const reference to matrix whose columns are the pseudo-eigenvectors. * - * More precisely, if the diagonal matrix of the eigen values is:\n - * \f$ - * \left[ \begin{array}{cccccc} - * u+iv & & & & & \\ - * & u-iv & & & & \\ - * & & a+ib & & & \\ - * & & & a-ib & & \\ - * & & & & x & \\ - * & & & & & y \\ - * \end{array} \right] - * \f$ \n - * then, we have:\n - * \f$ - * D =\left[ \begin{array}{cccccc} - * u & v & & & & \\ - * -v & u & & & & \\ - * & & a & b & & \\ - * & & -b & a & & \\ - * & & & & x & \\ - * & & & & & y \\ - * \end{array} \right] - * \f$ + * \pre Either the constructor EigenSolver(const MatrixType&) or + * the member function compute(const MatrixType&) has been called + * before. * - * \sa pseudoEigenvalueMatrix() + * The real matrix \f$ V \f$ returned by this function and the + * block-diagonal matrix \f$ D \f$ returned by pseudoEigenvalueMatrix() + * satisfy \f$ AV = VD \f$. + * + * Example: \include EigenSolver_pseudoEigenvectors.cpp + * Output: \verbinclude EigenSolver_pseudoEigenvectors.out + * + * \sa pseudoEigenvalueMatrix(), eigenvectors() */ const MatrixType& pseudoEigenvectors() const { @@ -115,19 +180,72 @@ template class EigenSolver return m_eivec; } + /** \brief Returns the block-diagonal matrix in the pseudo-eigendecomposition. + * + * \returns A block-diagonal matrix. + * + * \pre Either the constructor EigenSolver(const MatrixType&) or the + * member function compute(const MatrixType&) has been called before. + * + * The matrix \f$ D \f$ returned by this function is real and + * block-diagonal. The blocks on the diagonal are either 1-by-1 or 2-by-2 + * blocks of the form + * \f$ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f$. + * The matrix \f$ D \f$ and the matrix \f$ V \f$ returned by + * pseudoEigenvectors() satisfy \f$ AV = VD \f$. + * + * \sa pseudoEigenvectors() for an example, eigenvalues() + */ MatrixType pseudoEigenvalueMatrix() const; - /** \returns the eigenvalues as a column vector */ + /** \brief Returns the eigenvalues of given matrix. + * + * \returns Column vector containing the eigenvalues. + * + * \pre Either the constructor EigenSolver(const MatrixType&) or the + * member function compute(const MatrixType&) has been called before. + * + * The eigenvalues are repeated according to their algebraic multiplicity, + * so there are as many eigenvalues as rows in the matrix. + * + * Example: \include EigenSolver_eigenvalues.cpp + * Output: \verbinclude EigenSolver_eigenvalues.out + * + * \sa eigenvectors(), pseudoEigenvalueMatrix(), + * MatrixBase::eigenvalues() + */ EigenvalueType eigenvalues() const { ei_assert(m_isInitialized && "EigenSolver is not initialized."); return m_eivalues; } + /** \brief Computes eigendecomposition of given matrix. + * + * \param[in] matrix Square matrix whose eigendecomposition is to be computed. + * \returns Reference to \c *this + * + * This function computes the eigenvalues and eigenvectors of \p matrix. + * The eigenvalues() and eigenvectors() functions can be used to retrieve + * the computed eigendecomposition. + * + * The matrix is first reduced to Schur form. The Schur decomposition is + * then used to compute the eigenvalues and eigenvectors. + * + * The cost of the computation is dominated by the cost of the Schur + * decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$ is the size of + * the matrix. + * + * This method reuses of the allocated data in the EigenSolver object. + * + * Example: \include EigenSolver_compute.cpp + * Output: \verbinclude EigenSolver_compute.out + */ EigenSolver& compute(const MatrixType& matrix); private: + typedef Matrix RealVectorType; void orthes(MatrixType& matH, RealVectorType& ort); void hqr2(MatrixType& matH); @@ -137,10 +255,6 @@ template class EigenSolver bool m_isInitialized; }; -/** \returns the real block diagonal matrix D of the eigenvalues. - * - * See pseudoEigenvectors() for the details. - */ template MatrixType EigenSolver::pseudoEigenvalueMatrix() const { @@ -161,12 +275,8 @@ MatrixType EigenSolver::pseudoEigenvalueMatrix() const return matD; } -/** \returns the normalized complex eigenvectors as a matrix of column vectors. - * - * \sa eigenvalues(), pseudoEigenvectors() - */ template -typename EigenSolver::EigenvectorType EigenSolver::eigenvectors(void) const +typename EigenSolver::EigenvectorType EigenSolver::eigenvectors() const { ei_assert(m_isInitialized && "EigenSolver is not initialized."); int n = m_eivec.cols(); diff --git a/doc/snippets/EigenSolver_EigenSolver_MatrixType.cpp b/doc/snippets/EigenSolver_EigenSolver_MatrixType.cpp new file mode 100644 index 000000000..c1d9fa879 --- /dev/null +++ b/doc/snippets/EigenSolver_EigenSolver_MatrixType.cpp @@ -0,0 +1,16 @@ +MatrixXd A = MatrixXd::Random(6,6); +cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl; + +EigenSolver es(A); +cout << "The eigenvalues of A are:" << endl << es.eigenvalues() << endl; +cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl; + +complex lambda = es.eigenvalues()[0]; +cout << "Consider the first eigenvalue, lambda = " << lambda << endl; +VectorXcd v = es.eigenvectors().col(0); +cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl; +cout << "... and A * v = " << endl << A.cast >() * v << endl << endl; + +MatrixXcd D = es.eigenvalues().asDiagonal(); +MatrixXcd V = es.eigenvectors(); +cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl; diff --git a/doc/snippets/EigenSolver_compute.cpp b/doc/snippets/EigenSolver_compute.cpp new file mode 100644 index 000000000..06138f608 --- /dev/null +++ b/doc/snippets/EigenSolver_compute.cpp @@ -0,0 +1,6 @@ +EigenSolver es; +MatrixXf A = MatrixXf::Random(4,4); +es.compute(A); +cout << "The eigenvalues of A are: " << es.eigenvalues().transpose() << endl; +es.compute(A + MatrixXf::Identity(4,4)); // re-use es to compute eigenvalues of A+I +cout << "The eigenvalues of A+I are: " << es.eigenvalues().transpose() << endl; diff --git a/doc/snippets/EigenSolver_eigenvalues.cpp b/doc/snippets/EigenSolver_eigenvalues.cpp new file mode 100644 index 000000000..8d83ea982 --- /dev/null +++ b/doc/snippets/EigenSolver_eigenvalues.cpp @@ -0,0 +1,4 @@ +MatrixXd ones = MatrixXd::Ones(3,3); +EigenSolver es(ones); +cout << "The eigenvalues of the 3x3 matrix of ones are:" + << endl << es.eigenvalues() << endl; diff --git a/doc/snippets/EigenSolver_eigenvectors.cpp b/doc/snippets/EigenSolver_eigenvectors.cpp new file mode 100644 index 000000000..0fad4dadb --- /dev/null +++ b/doc/snippets/EigenSolver_eigenvectors.cpp @@ -0,0 +1,4 @@ +MatrixXd ones = MatrixXd::Ones(3,3); +EigenSolver es(ones); +cout << "The first eigenvector of the 3x3 matrix of ones is:" + << endl << es.eigenvectors().col(1) << endl; diff --git a/doc/snippets/EigenSolver_pseudoEigenvectors.cpp b/doc/snippets/EigenSolver_pseudoEigenvectors.cpp new file mode 100644 index 000000000..85e2569df --- /dev/null +++ b/doc/snippets/EigenSolver_pseudoEigenvectors.cpp @@ -0,0 +1,9 @@ +MatrixXd A = MatrixXd::Random(6,6); +cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl; + +EigenSolver es(A); +MatrixXd D = es.pseudoEigenvalueMatrix(); +MatrixXd V = es.pseudoEigenvectors(); +cout << "The pseudo-eigenvalue matrix D is:" << endl << D << endl; +cout << "The pseudo-eigenvector matrix V is:" << endl << V << endl; +cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl; From 8cfa672fe0bb70a3262b0287b298c09c8d982d21 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Wed, 31 Mar 2010 21:50:18 +0100 Subject: [PATCH 15/19] Use HessenbergDecomposition in EigenSolver. --- Eigen/src/Eigenvalues/EigenSolver.h | 73 ++--------------------------- 1 file changed, 5 insertions(+), 68 deletions(-) diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index 84f5c408f..36a9acb4d 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -25,6 +25,8 @@ #ifndef EIGEN_EIGENSOLVER_H #define EIGEN_EIGENSOLVER_H +#include "./HessenbergDecomposition.h" + /** \eigenvalues_module \ingroup Eigenvalues_Module * \nonstableyet * @@ -310,13 +312,11 @@ EigenSolver& EigenSolver::compute(const MatrixType& matr assert(matrix.cols() == matrix.rows()); int n = matrix.cols(); m_eivalues.resize(n,1); - m_eivec.resize(n,n); - - MatrixType matH = matrix; - RealVectorType ort(n); // Reduce to Hessenberg form. - orthes(matH, ort); + HessenbergDecomposition hd(matrix); + MatrixType matH = hd.matrixH(); + m_eivec = hd.matrixQ(); // Reduce Hessenberg to real Schur form. hqr2(matH); @@ -325,69 +325,6 @@ EigenSolver& EigenSolver::compute(const MatrixType& matr return *this; } -// Nonsymmetric reduction to Hessenberg form. -template -void EigenSolver::orthes(MatrixType& matH, RealVectorType& ort) -{ - // This is derived from the Algol procedures orthes and ortran, - // by Martin and Wilkinson, Handbook for Auto. Comp., - // Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutines in EISPACK. - - int n = m_eivec.cols(); - int low = 0; - int high = n-1; - - for (int m = low+1; m <= high-1; ++m) - { - // Scale column. - RealScalar scale = matH.block(m, m-1, high-m+1, 1).cwiseAbs().sum(); - if (scale != 0.0) - { - // Compute Householder transformation. - 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); - } - RealScalar g = ei_sqrt(h); - if (ort.coeff(m) > 0) - g = -g; - h = h - ort.coeff(m) * g; - ort.coeffRef(m) = ort.coeff(m) - g; - - // Apply Householder similarity transformation - // H = (I-u*u'/h)*H*(I-u*u')/h) - int bSize = high-m+1; - matH.block(m, m, bSize, n-m).noalias() -= ((ort.segment(m, bSize)/h) - * (ort.segment(m, bSize).transpose() * matH.block(m, m, bSize, n-m))); - - matH.block(0, m, high+1, bSize).noalias() -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize)) - * (ort.segment(m, bSize)/h).transpose()); - - ort.coeffRef(m) = scale*ort.coeff(m); - matH.coeffRef(m,m-1) = scale*g; - } - } - - // Accumulate transformations (Algol's ortran). - m_eivec.setIdentity(); - - for (int m = high-1; m >= low+1; m--) - { - if (matH.coeff(m,m-1) != 0.0) - { - ort.segment(m+1, high-m) = matH.col(m-1).segment(m+1, high-m); - - int bSize = high-m+1; - m_eivec.block(m, m, bSize, bSize).noalias() += ( (ort.segment(m, bSize) / (matH.coeff(m,m-1) * ort.coeff(m))) - * (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) ); - } - } -} - // Complex scalar division. template std::complex cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi) From 3a14a135331519f85b6490e803403d054598329e Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Thu, 1 Apr 2010 12:32:56 +0100 Subject: [PATCH 16/19] Split computation of real Schur form in EigenSolver to its own class. This is done with the minimal amount of work, so the result is very rough. --- Eigen/Eigenvalues | 1 + Eigen/src/Eigenvalues/EigenSolver.h | 301 ++------------------- Eigen/src/Eigenvalues/RealSchur.h | 394 ++++++++++++++++++++++++++++ 3 files changed, 414 insertions(+), 282 deletions(-) create mode 100644 Eigen/src/Eigenvalues/RealSchur.h diff --git a/Eigen/Eigenvalues b/Eigen/Eigenvalues index 1ae0cf098..986f31196 100644 --- a/Eigen/Eigenvalues +++ b/Eigen/Eigenvalues @@ -36,6 +36,7 @@ namespace Eigen { */ #include "src/Eigenvalues/Tridiagonalization.h" +#include "src/Eigenvalues/RealSchur.h" #include "src/Eigenvalues/EigenSolver.h" #include "src/Eigenvalues/SelfAdjointEigenSolver.h" #include "src/Eigenvalues/HessenbergDecomposition.h" diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index 36a9acb4d..efd31f18c 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -2,6 +2,7 @@ // for linear algebra. // // Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2010 Jitse Niesen // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -25,7 +26,7 @@ #ifndef EIGEN_EIGENSOLVER_H #define EIGEN_EIGENSOLVER_H -#include "./HessenbergDecomposition.h" +#include "./RealSchur.h" /** \eigenvalues_module \ingroup Eigenvalues_Module * \nonstableyet @@ -246,10 +247,7 @@ template class EigenSolver EigenSolver& compute(const MatrixType& matrix); private: - - typedef Matrix RealVectorType; - void orthes(MatrixType& matH, RealVectorType& ort); - void hqr2(MatrixType& matH); + void hqr2_step2(MatrixType& matH); protected: MatrixType m_eivec; @@ -310,16 +308,15 @@ template EigenSolver& EigenSolver::compute(const MatrixType& matrix) { assert(matrix.cols() == matrix.rows()); - int n = matrix.cols(); - m_eivalues.resize(n,1); - // Reduce to Hessenberg form. - HessenbergDecomposition hd(matrix); - MatrixType matH = hd.matrixH(); - m_eivec = hd.matrixQ(); + // Reduce to real Schur form. + RealSchur rs(matrix); + MatrixType matH = rs.matrixT(); + m_eivec = rs.matrixU(); + m_eivalues = rs.eigenvalues(); - // Reduce Hessenberg to real Schur form. - hqr2(matH); + // Compute eigenvectors. + hqr2_step2(matH); m_isInitialized = true; return *this; @@ -345,289 +342,29 @@ std::complex cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi) } -// Nonsymmetric reduction from Hessenberg to real Schur form. template -void EigenSolver::hqr2(MatrixType& matH) +void EigenSolver::hqr2_step2(MatrixType& matH) { - // This is derived from the Algol procedure hqr2, - // by Martin and Wilkinson, Handbook for Auto. Comp., - // Vol.ii-Linear Algebra, and the corresponding - // Fortran subroutine in EISPACK. + const int nn = m_eivec.cols(); + const int low = 0; + const int high = nn-1; + const Scalar eps = ei_pow(Scalar(2),ei_is_same_type::ret ? Scalar(-23) : Scalar(-52)); + Scalar p, q, r=0, s=0, t, w, x, y, z=0; - // Initialize - int nn = m_eivec.cols(); - int n = nn-1; - int low = 0; - int high = nn-1; - Scalar eps = ei_pow(Scalar(2),ei_is_same_type::ret ? Scalar(-23) : Scalar(-52)); - Scalar exshift = 0.0; - Scalar p=0,q=0,r=0,s=0,z=0,t,w,x,y; - - // Store roots isolated by balanc and compute matrix norm - // FIXME to be efficient the following would requires a triangular reduxion code - // Scalar norm = matH.upper().cwiseAbs().sum() + matH.corner(BottomLeft,n,n).diagonal().cwiseAbs().sum(); + // inefficient! this is already computed in RealSchur Scalar norm = 0.0; for (int j = 0; j < nn; ++j) { - // FIXME what's the purpose of the following since the condition is always false - if ((j < low) || (j > high)) - { - m_eivalues.coeffRef(j) = Complex(matH.coeff(j,j), 0.0); - } norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwiseAbs().sum(); } - - // Outer loop over eigenvalue index - int iter = 0; - while (n >= low) - { - // Look for single small sub-diagonal element - int l = n; - while (l > low) - { - s = ei_abs(matH.coeff(l-1,l-1)) + ei_abs(matH.coeff(l,l)); - if (s == 0.0) - s = norm; - if (ei_abs(matH.coeff(l,l-1)) < eps * s) - break; - l--; - } - - // Check for convergence - // One root found - if (l == n) - { - matH.coeffRef(n,n) = matH.coeff(n,n) + exshift; - m_eivalues.coeffRef(n) = Complex(matH.coeff(n,n), 0.0); - n--; - iter = 0; - } - else if (l == n-1) // Two roots found - { - w = matH.coeff(n,n-1) * matH.coeff(n-1,n); - p = (matH.coeff(n-1,n-1) - matH.coeff(n,n)) * Scalar(0.5); - q = p * p + w; - z = ei_sqrt(ei_abs(q)); - matH.coeffRef(n,n) = matH.coeff(n,n) + exshift; - matH.coeffRef(n-1,n-1) = matH.coeff(n-1,n-1) + exshift; - x = matH.coeff(n,n); - - // Scalar pair - if (q >= 0) - { - if (p >= 0) - z = p + z; - else - z = p - z; - - m_eivalues.coeffRef(n-1) = Complex(x + z, 0.0); - m_eivalues.coeffRef(n) = Complex(z!=0.0 ? x - w / z : m_eivalues.coeff(n-1).real(), 0.0); - - x = matH.coeff(n,n-1); - s = ei_abs(x) + ei_abs(z); - p = x / s; - q = z / s; - r = ei_sqrt(p * p+q * q); - p = p / r; - q = q / r; - - // Row modification - for (int j = n-1; j < nn; ++j) - { - z = matH.coeff(n-1,j); - matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j); - matH.coeffRef(n,j) = q * matH.coeff(n,j) - p * z; - } - - // Column modification - for (int i = 0; i <= n; ++i) - { - z = matH.coeff(i,n-1); - matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n); - matH.coeffRef(i,n) = q * matH.coeff(i,n) - p * z; - } - - // Accumulate transformations - for (int i = low; i <= high; ++i) - { - z = m_eivec.coeff(i,n-1); - m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n); - m_eivec.coeffRef(i,n) = q * m_eivec.coeff(i,n) - p * z; - } - } - else // Complex pair - { - m_eivalues.coeffRef(n-1) = Complex(x + p, z); - m_eivalues.coeffRef(n) = Complex(x + p, -z); - } - n = n - 2; - iter = 0; - } - else // No convergence yet - { - // Form shift - x = matH.coeff(n,n); - y = 0.0; - w = 0.0; - if (l < n) - { - y = matH.coeff(n-1,n-1); - w = matH.coeff(n,n-1) * matH.coeff(n-1,n); - } - - // Wilkinson's original ad hoc shift - if (iter == 10) - { - exshift += x; - for (int i = low; i <= n; ++i) - matH.coeffRef(i,i) -= x; - s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2)); - x = y = Scalar(0.75) * s; - w = Scalar(-0.4375) * s * s; - } - - // MATLAB's new ad hoc shift - if (iter == 30) - { - s = Scalar((y - x) / 2.0); - s = s * s + w; - if (s > 0) - { - s = ei_sqrt(s); - if (y < x) - s = -s; - s = Scalar(x - w / ((y - x) / 2.0 + s)); - for (int i = low; i <= n; ++i) - matH.coeffRef(i,i) -= s; - exshift += s; - x = y = w = Scalar(0.964); - } - } - - iter = iter + 1; // (Could check iteration count here.) - - // Look for two consecutive small sub-diagonal elements - int m = n-2; - while (m >= l) - { - z = matH.coeff(m,m); - r = x - z; - s = y - z; - p = (r * s - w) / matH.coeff(m+1,m) + matH.coeff(m,m+1); - q = matH.coeff(m+1,m+1) - z - r - s; - r = matH.coeff(m+2,m+1); - s = ei_abs(p) + ei_abs(q) + ei_abs(r); - p = p / s; - q = q / s; - r = r / s; - if (m == l) { - break; - } - if (ei_abs(matH.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) < - eps * (ei_abs(p) * (ei_abs(matH.coeff(m-1,m-1)) + ei_abs(z) + - ei_abs(matH.coeff(m+1,m+1))))) - { - break; - } - m--; - } - - for (int i = m+2; i <= n; ++i) - { - matH.coeffRef(i,i-2) = 0.0; - if (i > m+2) - matH.coeffRef(i,i-3) = 0.0; - } - - // Double QR step involving rows l:n and columns m:n - for (int k = m; k <= n-1; ++k) - { - int notlast = (k != n-1); - if (k != m) { - p = matH.coeff(k,k-1); - q = matH.coeff(k+1,k-1); - r = notlast ? matH.coeff(k+2,k-1) : Scalar(0); - x = ei_abs(p) + ei_abs(q) + ei_abs(r); - if (x != 0.0) - { - p = p / x; - q = q / x; - r = r / x; - } - } - - if (x == 0.0) - break; - - s = ei_sqrt(p * p + q * q + r * r); - - if (p < 0) - s = -s; - - if (s != 0) - { - if (k != m) - matH.coeffRef(k,k-1) = -s * x; - else if (l != m) - matH.coeffRef(k,k-1) = -matH.coeff(k,k-1); - - p = p + s; - x = p / s; - y = q / s; - z = r / s; - q = q / p; - r = r / p; - - // Row modification - for (int j = k; j < nn; ++j) - { - p = matH.coeff(k,j) + q * matH.coeff(k+1,j); - if (notlast) - { - p = p + r * matH.coeff(k+2,j); - matH.coeffRef(k+2,j) = matH.coeff(k+2,j) - p * z; - } - matH.coeffRef(k,j) = matH.coeff(k,j) - p * x; - matH.coeffRef(k+1,j) = matH.coeff(k+1,j) - p * y; - } - - // Column modification - for (int i = 0; i <= std::min(n,k+3); ++i) - { - p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1); - if (notlast) - { - p = p + z * matH.coeff(i,k+2); - matH.coeffRef(i,k+2) = matH.coeff(i,k+2) - p * r; - } - matH.coeffRef(i,k) = matH.coeff(i,k) - p; - matH.coeffRef(i,k+1) = matH.coeff(i,k+1) - p * q; - } - - // Accumulate transformations - for (int i = low; i <= high; ++i) - { - p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1); - if (notlast) - { - p = p + z * m_eivec.coeff(i,k+2); - m_eivec.coeffRef(i,k+2) = m_eivec.coeff(i,k+2) - p * r; - } - m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) - p; - m_eivec.coeffRef(i,k+1) = m_eivec.coeff(i,k+1) - p * q; - } - } // (s != 0) - } // k loop - } // check convergence - } // while (n >= low) - + // Backsubstitute to find vectors of upper triangular form if (norm == 0.0) { return; } - for (n = nn-1; n >= 0; n--) + for (int n = nn-1; n >= 0; n--) { p = m_eivalues.coeff(n).real(); q = m_eivalues.coeff(n).imag(); diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h new file mode 100644 index 000000000..3aaeaa753 --- /dev/null +++ b/Eigen/src/Eigenvalues/RealSchur.h @@ -0,0 +1,394 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2010 Jitse Niesen +// +// 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 . + +#ifndef EIGEN_REAL_SCHUR_H +#define EIGEN_REAL_SCHUR_H + +#include "./HessenbergDecomposition.h" + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * \nonstableyet + * + * \class RealSchur + * + * \brief Performs a real Schur decomposition of a square matrix + */ +template class RealSchur +{ + public: + typedef _MatrixType MatrixType; + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + typedef typename MatrixType::Scalar Scalar; + typedef std::complex::Real> Complex; + typedef Matrix EigenvalueType; + + /** \brief Constructor; computes Schur decomposition of given matrix. */ + RealSchur(const MatrixType& matrix) + : matH(matrix.rows(),matrix.cols()), + m_eivec(matrix.rows(),matrix.cols()), + m_eivalues(matrix.rows()), + m_isInitialized(false) + { + compute(matrix); + } + + /** \brief Returns the orthogonal matrix in the Schur decomposition. */ + const MatrixType& matrixU() const + { + ei_assert(m_isInitialized && "RealSchur is not initialized."); + return m_eivec; + } + + /** \brief Returns the quasi-triangular matrix in the Schur decomposition. */ + const MatrixType& matrixT() const + { + ei_assert(m_isInitialized && "RealSchur is not initialized."); + return matH; + } + + /** \brief Returns vector of eigenvalues. + * + * This function will likely be removed. */ + const EigenvalueType& eigenvalues() const + { + ei_assert(m_isInitialized && "RealSchur is not initialized."); + return m_eivalues; + } + + /** \brief Computes Schur decomposition of given matrix. */ + void compute(const MatrixType& matrix); + + private: + + MatrixType matH; + MatrixType m_eivec; + EigenvalueType m_eivalues; + bool m_isInitialized; + + void hqr2(); +}; + + +template +void RealSchur::compute(const MatrixType& matrix) +{ + assert(matrix.cols() == matrix.rows()); + + // Reduce to Hessenberg form + // TODO skip Q if skipU = true + HessenbergDecomposition hess(matrix); + matH = hess.matrixH(); + m_eivec = hess.matrixQ(); + + // Reduce to Real Schur form + hqr2(); + + m_isInitialized = true; +} + + +template +void RealSchur::hqr2() +{ + // This is derived from the Algol procedure hqr2, + // by Martin and Wilkinson, Handbook for Auto. Comp., + // Vol.ii-Linear Algebra, and the corresponding + // Fortran subroutine in EISPACK. + + // Initialize + int nn = m_eivec.cols(); + int n = nn-1; + const int low = 0; + const int high = nn-1; + const Scalar eps = ei_pow(Scalar(2),ei_is_same_type::ret ? Scalar(-23) : Scalar(-52)); + Scalar exshift = 0.0; + Scalar p=0,q=0,r=0,s=0,z=0,w,x,y; + + // Store roots isolated by balanc and compute matrix norm + // FIXME to be efficient the following would requires a triangular reduxion code + // Scalar norm = matH.upper().cwiseAbs().sum() + matH.corner(BottomLeft,n,n).diagonal().cwiseAbs().sum(); + Scalar norm = 0.0; + for (int j = 0; j < nn; ++j) + { + // FIXME what's the purpose of the following since the condition is always false + if ((j < low) || (j > high)) + { + m_eivalues.coeffRef(j) = Complex(matH.coeff(j,j), 0.0); + } + norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwiseAbs().sum(); + } + + // Outer loop over eigenvalue index + int iter = 0; + while (n >= low) + { + // Look for single small sub-diagonal element + int l = n; + while (l > low) + { + s = ei_abs(matH.coeff(l-1,l-1)) + ei_abs(matH.coeff(l,l)); + if (s == 0.0) + s = norm; + if (ei_abs(matH.coeff(l,l-1)) < eps * s) + break; + l--; + } + + // Check for convergence + // One root found + if (l == n) + { + matH.coeffRef(n,n) = matH.coeff(n,n) + exshift; + m_eivalues.coeffRef(n) = Complex(matH.coeff(n,n), 0.0); + n--; + iter = 0; + } + else if (l == n-1) // Two roots found + { + w = matH.coeff(n,n-1) * matH.coeff(n-1,n); + p = (matH.coeff(n-1,n-1) - matH.coeff(n,n)) * Scalar(0.5); + q = p * p + w; + z = ei_sqrt(ei_abs(q)); + matH.coeffRef(n,n) = matH.coeff(n,n) + exshift; + matH.coeffRef(n-1,n-1) = matH.coeff(n-1,n-1) + exshift; + x = matH.coeff(n,n); + + // Scalar pair + if (q >= 0) + { + if (p >= 0) + z = p + z; + else + z = p - z; + + m_eivalues.coeffRef(n-1) = Complex(x + z, 0.0); + m_eivalues.coeffRef(n) = Complex(z!=0.0 ? x - w / z : m_eivalues.coeff(n-1).real(), 0.0); + + x = matH.coeff(n,n-1); + s = ei_abs(x) + ei_abs(z); + p = x / s; + q = z / s; + r = ei_sqrt(p * p+q * q); + p = p / r; + q = q / r; + + // Row modification + for (int j = n-1; j < nn; ++j) + { + z = matH.coeff(n-1,j); + matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j); + matH.coeffRef(n,j) = q * matH.coeff(n,j) - p * z; + } + + // Column modification + for (int i = 0; i <= n; ++i) + { + z = matH.coeff(i,n-1); + matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n); + matH.coeffRef(i,n) = q * matH.coeff(i,n) - p * z; + } + + // Accumulate transformations + for (int i = low; i <= high; ++i) + { + z = m_eivec.coeff(i,n-1); + m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n); + m_eivec.coeffRef(i,n) = q * m_eivec.coeff(i,n) - p * z; + } + } + else // Complex pair + { + m_eivalues.coeffRef(n-1) = Complex(x + p, z); + m_eivalues.coeffRef(n) = Complex(x + p, -z); + } + n = n - 2; + iter = 0; + } + else // No convergence yet + { + // Form shift + x = matH.coeff(n,n); + y = 0.0; + w = 0.0; + if (l < n) + { + y = matH.coeff(n-1,n-1); + w = matH.coeff(n,n-1) * matH.coeff(n-1,n); + } + + // Wilkinson's original ad hoc shift + if (iter == 10) + { + exshift += x; + for (int i = low; i <= n; ++i) + matH.coeffRef(i,i) -= x; + s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2)); + x = y = Scalar(0.75) * s; + w = Scalar(-0.4375) * s * s; + } + + // MATLAB's new ad hoc shift + if (iter == 30) + { + s = Scalar((y - x) / 2.0); + s = s * s + w; + if (s > 0) + { + s = ei_sqrt(s); + if (y < x) + s = -s; + s = Scalar(x - w / ((y - x) / 2.0 + s)); + for (int i = low; i <= n; ++i) + matH.coeffRef(i,i) -= s; + exshift += s; + x = y = w = Scalar(0.964); + } + } + + iter = iter + 1; // (Could check iteration count here.) + + // Look for two consecutive small sub-diagonal elements + int m = n-2; + while (m >= l) + { + z = matH.coeff(m,m); + r = x - z; + s = y - z; + p = (r * s - w) / matH.coeff(m+1,m) + matH.coeff(m,m+1); + q = matH.coeff(m+1,m+1) - z - r - s; + r = matH.coeff(m+2,m+1); + s = ei_abs(p) + ei_abs(q) + ei_abs(r); + p = p / s; + q = q / s; + r = r / s; + if (m == l) { + break; + } + if (ei_abs(matH.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) < + eps * (ei_abs(p) * (ei_abs(matH.coeff(m-1,m-1)) + ei_abs(z) + + ei_abs(matH.coeff(m+1,m+1))))) + { + break; + } + m--; + } + + for (int i = m+2; i <= n; ++i) + { + matH.coeffRef(i,i-2) = 0.0; + if (i > m+2) + matH.coeffRef(i,i-3) = 0.0; + } + + // Double QR step involving rows l:n and columns m:n + for (int k = m; k <= n-1; ++k) + { + int notlast = (k != n-1); + if (k != m) { + p = matH.coeff(k,k-1); + q = matH.coeff(k+1,k-1); + r = notlast ? matH.coeff(k+2,k-1) : Scalar(0); + x = ei_abs(p) + ei_abs(q) + ei_abs(r); + if (x != 0.0) + { + p = p / x; + q = q / x; + r = r / x; + } + } + + if (x == 0.0) + break; + + s = ei_sqrt(p * p + q * q + r * r); + + if (p < 0) + s = -s; + + if (s != 0) + { + if (k != m) + matH.coeffRef(k,k-1) = -s * x; + else if (l != m) + matH.coeffRef(k,k-1) = -matH.coeff(k,k-1); + + p = p + s; + x = p / s; + y = q / s; + z = r / s; + q = q / p; + r = r / p; + + // Row modification + for (int j = k; j < nn; ++j) + { + p = matH.coeff(k,j) + q * matH.coeff(k+1,j); + if (notlast) + { + p = p + r * matH.coeff(k+2,j); + matH.coeffRef(k+2,j) = matH.coeff(k+2,j) - p * z; + } + matH.coeffRef(k,j) = matH.coeff(k,j) - p * x; + matH.coeffRef(k+1,j) = matH.coeff(k+1,j) - p * y; + } + + // Column modification + for (int i = 0; i <= std::min(n,k+3); ++i) + { + p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1); + if (notlast) + { + p = p + z * matH.coeff(i,k+2); + matH.coeffRef(i,k+2) = matH.coeff(i,k+2) - p * r; + } + matH.coeffRef(i,k) = matH.coeff(i,k) - p; + matH.coeffRef(i,k+1) = matH.coeff(i,k+1) - p * q; + } + + // Accumulate transformations + for (int i = low; i <= high; ++i) + { + p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1); + if (notlast) + { + p = p + z * m_eivec.coeff(i,k+2); + m_eivec.coeffRef(i,k+2) = m_eivec.coeff(i,k+2) - p * r; + } + m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) - p; + m_eivec.coeffRef(i,k+1) = m_eivec.coeff(i,k+1) - p * q; + } + } // (s != 0) + } // k loop + } // check convergence + } // while (n >= low) +} + +#endif // EIGEN_REAL_SCHUR_H From 9d6afdeb22d1ccc17a2d97966163c6d8f7651047 Mon Sep 17 00:00:00 2001 From: Hauke Heibel Date: Thu, 1 Apr 2010 15:10:52 +0200 Subject: [PATCH 17/19] ei_psqrt fix for zero input --- Eigen/src/Core/arch/SSE/MathFunctions.h | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/Eigen/src/Core/arch/SSE/MathFunctions.h b/Eigen/src/Core/arch/SSE/MathFunctions.h index 3c0020248..99662eb6d 100644 --- a/Eigen/src/Core/arch/SSE/MathFunctions.h +++ b/Eigen/src/Core/arch/SSE/MathFunctions.h @@ -369,10 +369,14 @@ static EIGEN_DONT_INLINE EIGEN_UNUSED Packet4f ei_pcos(Packet4f x) // For detail see here: http://www.beyond3d.com/content/articles/8/ static EIGEN_UNUSED Packet4f ei_psqrt(Packet4f _x) { - Packet4f half = ei_pmul(_x, ei_pset1(.5f)); - Packet4f x = _mm_rsqrt_ps(_x); - x = ei_pmul(x, ei_psub(ei_pset1(1.5f), ei_pmul(half, ei_pmul(x,x)))); - return ei_pmul(_x,x); + Packet4f half = ei_pmul(_x, ei_pset1(.5f)); + + /* select only the inverse sqrt of non-zero inputs */ + Packet4f non_zero_mask = _mm_cmpgt_ps(_x, ei_pset1(std::numeric_limits::epsilon())); + Packet4f x = _mm_and_ps(non_zero_mask, _mm_rsqrt_ps(_x)); + + x = ei_pmul(x, ei_psub(ei_pset1(1.5f), ei_pmul(half, ei_pmul(x,x)))); + return ei_pmul(_x,x); } #endif // EIGEN_MATH_FUNCTIONS_SSE_H From a16a36ecf2596fadffd6474d059c334baa905c06 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Fri, 2 Apr 2010 14:32:20 +0100 Subject: [PATCH 18/19] Add tests for real and complex Schur; extend test for Hessenberg. Make a minor correction to the ComplexSchur class. --- Eigen/src/Eigenvalues/ComplexSchur.h | 6 +-- test/CMakeLists.txt | 4 +- test/hessenberg.cpp | 38 ++++++++++---- test/schur_complex.cpp | 67 +++++++++++++++++++++++++ test/schur_real.cpp | 75 ++++++++++++++++++++++++++++ 5 files changed, 176 insertions(+), 14 deletions(-) create mode 100644 test/schur_complex.cpp create mode 100644 test/schur_real.cpp diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h index 1ab2a0287..04b9d6b80 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/Eigen/src/Eigenvalues/ComplexSchur.h @@ -86,7 +86,7 @@ template class ComplexSchur /** \brief Default constructor. * - * \param [in] size The size of the matrix whose Schur decomposition will be computed. + * \param [in] size Positive integer, size of the matrix whose Schur decomposition will be computed. * * The default constructor is useful in cases in which the user * intends to perform decompositions via compute(). The \p size @@ -95,7 +95,7 @@ template class ComplexSchur * * \sa compute() for an example. */ - ComplexSchur(int size = RowsAtCompileTime==Dynamic ? 0 : RowsAtCompileTime) + ComplexSchur(int size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : m_matT(size,size), m_matU(size,size), m_isInitialized(false), m_matUisUptodate(false) {} @@ -157,7 +157,7 @@ template class ComplexSchur */ const ComplexMatrixType& matrixT() const { - ei_assert(m_isInitialized && "ComplexShur is not initialized."); + ei_assert(m_isInitialized && "ComplexSchur is not initialized."); return m_matT; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 072f63e1d..f7d5ed8e9 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -138,7 +138,9 @@ ei_add_test(qr) ei_add_test(qr_colpivoting) ei_add_test(qr_fullpivoting) ei_add_test(upperbidiagonalization) -ei_add_test(hessenberg " " "${GSL_LIBRARIES}") +ei_add_test(hessenberg) +ei_add_test(schur_real) +ei_add_test(schur_complex) ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_complex) diff --git a/test/hessenberg.cpp b/test/hessenberg.cpp index d917be357..cba9c8fda 100644 --- a/test/hessenberg.cpp +++ b/test/hessenberg.cpp @@ -2,6 +2,7 @@ // for linear algebra. // // Copyright (C) 2009 Gael Guennebaud +// Copyright (C) 2010 Jitse Niesen // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -28,19 +29,36 @@ template void hessenberg(int size = Size) { typedef Matrix MatrixType; - MatrixType m = MatrixType::Random(size,size); - HessenbergDecomposition hess(m); - VERIFY_IS_APPROX(m, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint()); + // Test basic functionality: A = U H U* and H is Hessenberg + for(int counter = 0; counter < g_repeat; ++counter) { + MatrixType m = MatrixType::Random(size,size); + HessenbergDecomposition hess(m); + VERIFY_IS_APPROX(m, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint()); + MatrixType H = hess.matrixH(); + for(int row = 2; row < size; ++row) { + for(int col = 0; col < row-1; ++col) { + VERIFY(H(row,col) == (typename MatrixType::Scalar)0); + } + } + } + + // Test whether compute() and constructor returns same result + MatrixType A = MatrixType::Random(size, size); + HessenbergDecomposition cs1; + cs1.compute(A); + HessenbergDecomposition cs2(A); + VERIFY_IS_EQUAL(cs1.matrixQ(), cs2.matrixQ()); + VERIFY_IS_EQUAL(cs1.matrixH(), cs2.matrixH()); + + // TODO: Add tests for packedMatrix() and householderCoefficients() } void test_hessenberg() { - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_1(( hessenberg,1>() )); - CALL_SUBTEST_2(( hessenberg,2>() )); - CALL_SUBTEST_3(( hessenberg,4>() )); - CALL_SUBTEST_4(( hessenberg(ei_random(1,320)) )); - CALL_SUBTEST_5(( hessenberg,Dynamic>(ei_random(1,320)) )); - } + CALL_SUBTEST_1(( hessenberg,1>() )); + CALL_SUBTEST_2(( hessenberg,2>() )); + CALL_SUBTEST_3(( hessenberg,4>() )); + CALL_SUBTEST_4(( hessenberg(ei_random(1,320)) )); + CALL_SUBTEST_5(( hessenberg,Dynamic>(ei_random(1,320)) )); } diff --git a/test/schur_complex.cpp b/test/schur_complex.cpp new file mode 100644 index 000000000..3659a074c --- /dev/null +++ b/test/schur_complex.cpp @@ -0,0 +1,67 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2010 Jitse Niesen +// +// 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 . + +#include "main.h" +#include + +template void schur(int size = MatrixType::ColsAtCompileTime) +{ + typedef typename ComplexSchur::ComplexScalar ComplexScalar; + typedef typename ComplexSchur::ComplexMatrixType ComplexMatrixType; + + // Test basic functionality: T is triangular and A = U T U* + for(int counter = 0; counter < g_repeat; ++counter) { + MatrixType A = MatrixType::Random(size, size); + ComplexSchur schurOfA(A); + ComplexMatrixType U = schurOfA.matrixU(); + ComplexMatrixType T = schurOfA.matrixT(); + for(int row = 1; row < size; ++row) { + for(int col = 0; col < row; ++col) { + VERIFY(T(row,col) == (typename MatrixType::Scalar)0); + } + } + VERIFY_IS_APPROX(A.template cast(), U * T * U.adjoint()); + } + + // Test asserts when not initialized + ComplexSchur csUninitialized; + VERIFY_RAISES_ASSERT(csUninitialized.matrixT()); + VERIFY_RAISES_ASSERT(csUninitialized.matrixU()); + + // Test whether compute() and constructor returns same result + MatrixType A = MatrixType::Random(size, size); + ComplexSchur cs1; + cs1.compute(A); + ComplexSchur cs2(A); + VERIFY_IS_EQUAL(cs1.matrixT(), cs2.matrixT()); + VERIFY_IS_EQUAL(cs1.matrixU(), cs2.matrixU()); +} + +void test_schur_complex() +{ + CALL_SUBTEST_1(( schur() )); + CALL_SUBTEST_2(( schur(ei_random(1,50)) )); + CALL_SUBTEST_3(( schur, 1, 1> >() )); + CALL_SUBTEST_4(( schur >() )); +} diff --git a/test/schur_real.cpp b/test/schur_real.cpp new file mode 100644 index 000000000..77ef5e2dc --- /dev/null +++ b/test/schur_real.cpp @@ -0,0 +1,75 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2010 Jitse Niesen +// +// 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 . + +#include "main.h" +#include + +template void verifyIsQuasiTriangular(const MatrixType& T) +{ + const int size = T.cols(); + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + + // The "zeros" in the real Schur decomposition are only approximately zero + RealScalar norm = T.norm(); + + // Check T is lower Hessenberg + for(int row = 2; row < size; ++row) { + for(int col = 0; col < row - 1; ++col) { + VERIFY_IS_MUCH_SMALLER_THAN(T(row,col), norm); + } + } + + // Check that any non-zero on the subdiagonal is followed by a zero and is + // part of a 2x2 diagonal block with imaginary eigenvalues. + for(int row = 1; row < size; ++row) { + if (!test_ei_isMuchSmallerThan(T(row,row-1), norm)) { + VERIFY(row == size-1 || test_ei_isMuchSmallerThan(T(row+1,row), norm)); + Scalar tr = T(row-1,row-1) + T(row,row); + Scalar det = T(row-1,row-1) * T(row,row) - T(row-1,row) * T(row,row-1); + VERIFY(4 * det > tr * tr); + } + } +} + +template void schur(int size = MatrixType::ColsAtCompileTime) +{ + // Test basic functionality: T is quasi-triangular and A = U T U* + for(int counter = 0; counter < g_repeat; ++counter) { + MatrixType A = MatrixType::Random(size, size); + RealSchur schurOfA(A); + MatrixType U = schurOfA.matrixU(); + MatrixType T = schurOfA.matrixT(); + verifyIsQuasiTriangular(T); + VERIFY_IS_APPROX(A, U * T * U.transpose()); + } +} + +void test_schur_real() +{ + CALL_SUBTEST_1(( schur() )); + CALL_SUBTEST_2(( schur(ei_random(1,50)) )); + CALL_SUBTEST_3(( schur >() )); + CALL_SUBTEST_4(( schur >() )); +} From 79e1ce609319e749d8d9a9aa2ffbcf0600ab1b93 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Fri, 2 Apr 2010 21:05:32 +0100 Subject: [PATCH 19/19] RealSchur and EigenSolver: some straightforward renames. --- Eigen/src/Eigenvalues/EigenSolver.h | 16 +-- Eigen/src/Eigenvalues/RealSchur.h | 155 ++++++++++++++-------------- 2 files changed, 85 insertions(+), 86 deletions(-) diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index efd31f18c..44a0fd485 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -95,21 +95,21 @@ template class EigenSolver * \c float or \c double) and just \c Scalar if #Scalar is * complex. */ - typedef std::complex Complex; + typedef std::complex ComplexScalar; /** \brief Type for vector of eigenvalues as returned by eigenvalues(). * - * This is a column vector with entries of type #Complex. + * This is a column vector with entries of type #ComplexScalar. * The length of the vector is the size of \p _MatrixType. */ - typedef Matrix EigenvalueType; + typedef Matrix EigenvalueType; /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). * - * This is a square matrix with entries of type #Complex. + * This is a square matrix with entries of type #ComplexScalar. * The size is the same as the size of \p _MatrixType. */ - typedef Matrix EigenvectorType; + typedef Matrix EigenvectorType; /** \brief Default constructor. * @@ -286,15 +286,15 @@ typename EigenSolver::EigenvectorType EigenSolver::eigen if (ei_isMuchSmallerThan(ei_abs(ei_imag(m_eivalues.coeff(j))), ei_abs(ei_real(m_eivalues.coeff(j))))) { // we have a real eigen value - matV.col(j) = m_eivec.col(j).template cast(); + matV.col(j) = m_eivec.col(j).template cast(); } else { // we have a pair of complex eigen values for (int i=0; i class RealSchur MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; typedef typename MatrixType::Scalar Scalar; - typedef std::complex::Real> Complex; - typedef Matrix EigenvalueType; + typedef std::complex::Real> ComplexScalar; + typedef Matrix EigenvalueType; /** \brief Constructor; computes Schur decomposition of given matrix. */ RealSchur(const MatrixType& matrix) - : matH(matrix.rows(),matrix.cols()), - m_eivec(matrix.rows(),matrix.cols()), + : m_matT(matrix.rows(),matrix.cols()), + m_matU(matrix.rows(),matrix.cols()), m_eivalues(matrix.rows()), m_isInitialized(false) { @@ -64,14 +64,14 @@ template class RealSchur const MatrixType& matrixU() const { ei_assert(m_isInitialized && "RealSchur is not initialized."); - return m_eivec; + return m_matU; } /** \brief Returns the quasi-triangular matrix in the Schur decomposition. */ const MatrixType& matrixT() const { ei_assert(m_isInitialized && "RealSchur is not initialized."); - return matH; + return m_matT; } /** \brief Returns vector of eigenvalues. @@ -88,8 +88,8 @@ template class RealSchur private: - MatrixType matH; - MatrixType m_eivec; + MatrixType m_matT; + MatrixType m_matU; EigenvalueType m_eivalues; bool m_isInitialized; @@ -105,8 +105,8 @@ void RealSchur::compute(const MatrixType& matrix) // Reduce to Hessenberg form // TODO skip Q if skipU = true HessenbergDecomposition hess(matrix); - matH = hess.matrixH(); - m_eivec = hess.matrixQ(); + m_matT = hess.matrixH(); + m_matU = hess.matrixQ(); // Reduce to Real Schur form hqr2(); @@ -124,26 +124,25 @@ void RealSchur::hqr2() // Fortran subroutine in EISPACK. // Initialize - int nn = m_eivec.cols(); - int n = nn-1; + const int size = m_matU.cols(); + int n = size-1; const int low = 0; - const int high = nn-1; - const Scalar eps = ei_pow(Scalar(2),ei_is_same_type::ret ? Scalar(-23) : Scalar(-52)); + const int high = size-1; Scalar exshift = 0.0; Scalar p=0,q=0,r=0,s=0,z=0,w,x,y; // Store roots isolated by balanc and compute matrix norm // FIXME to be efficient the following would requires a triangular reduxion code - // Scalar norm = matH.upper().cwiseAbs().sum() + matH.corner(BottomLeft,n,n).diagonal().cwiseAbs().sum(); + // Scalar norm = m_matT.upper().cwiseAbs().sum() + m_matT.corner(BottomLeft,n,n).diagonal().cwiseAbs().sum(); Scalar norm = 0.0; - for (int j = 0; j < nn; ++j) + for (int j = 0; j < size; ++j) { // FIXME what's the purpose of the following since the condition is always false if ((j < low) || (j > high)) { - m_eivalues.coeffRef(j) = Complex(matH.coeff(j,j), 0.0); + m_eivalues.coeffRef(j) = ComplexScalar(m_matT.coeff(j,j), 0.0); } - norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwiseAbs().sum(); + norm += m_matT.row(j).segment(std::max(j-1,0), size-std::max(j-1,0)).cwiseAbs().sum(); } // Outer loop over eigenvalue index @@ -154,10 +153,10 @@ void RealSchur::hqr2() int l = n; while (l > low) { - s = ei_abs(matH.coeff(l-1,l-1)) + ei_abs(matH.coeff(l,l)); + s = ei_abs(m_matT.coeff(l-1,l-1)) + ei_abs(m_matT.coeff(l,l)); if (s == 0.0) s = norm; - if (ei_abs(matH.coeff(l,l-1)) < eps * s) + if (ei_abs(m_matT.coeff(l,l-1)) < NumTraits::epsilon() * s) break; l--; } @@ -166,20 +165,20 @@ void RealSchur::hqr2() // One root found if (l == n) { - matH.coeffRef(n,n) = matH.coeff(n,n) + exshift; - m_eivalues.coeffRef(n) = Complex(matH.coeff(n,n), 0.0); + m_matT.coeffRef(n,n) = m_matT.coeff(n,n) + exshift; + m_eivalues.coeffRef(n) = ComplexScalar(m_matT.coeff(n,n), 0.0); n--; iter = 0; } else if (l == n-1) // Two roots found { - w = matH.coeff(n,n-1) * matH.coeff(n-1,n); - p = (matH.coeff(n-1,n-1) - matH.coeff(n,n)) * Scalar(0.5); + w = m_matT.coeff(n,n-1) * m_matT.coeff(n-1,n); + p = (m_matT.coeff(n-1,n-1) - m_matT.coeff(n,n)) * Scalar(0.5); q = p * p + w; z = ei_sqrt(ei_abs(q)); - matH.coeffRef(n,n) = matH.coeff(n,n) + exshift; - matH.coeffRef(n-1,n-1) = matH.coeff(n-1,n-1) + exshift; - x = matH.coeff(n,n); + m_matT.coeffRef(n,n) = m_matT.coeff(n,n) + exshift; + m_matT.coeffRef(n-1,n-1) = m_matT.coeff(n-1,n-1) + exshift; + x = m_matT.coeff(n,n); // Scalar pair if (q >= 0) @@ -189,10 +188,10 @@ void RealSchur::hqr2() else z = p - z; - m_eivalues.coeffRef(n-1) = Complex(x + z, 0.0); - m_eivalues.coeffRef(n) = Complex(z!=0.0 ? x - w / z : m_eivalues.coeff(n-1).real(), 0.0); + m_eivalues.coeffRef(n-1) = ComplexScalar(x + z, 0.0); + m_eivalues.coeffRef(n) = ComplexScalar(z!=0.0 ? x - w / z : m_eivalues.coeff(n-1).real(), 0.0); - x = matH.coeff(n,n-1); + x = m_matT.coeff(n,n-1); s = ei_abs(x) + ei_abs(z); p = x / s; q = z / s; @@ -201,33 +200,33 @@ void RealSchur::hqr2() q = q / r; // Row modification - for (int j = n-1; j < nn; ++j) + for (int j = n-1; j < size; ++j) { - z = matH.coeff(n-1,j); - matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j); - matH.coeffRef(n,j) = q * matH.coeff(n,j) - p * z; + z = m_matT.coeff(n-1,j); + m_matT.coeffRef(n-1,j) = q * z + p * m_matT.coeff(n,j); + m_matT.coeffRef(n,j) = q * m_matT.coeff(n,j) - p * z; } // Column modification for (int i = 0; i <= n; ++i) { - z = matH.coeff(i,n-1); - matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n); - matH.coeffRef(i,n) = q * matH.coeff(i,n) - p * z; + z = m_matT.coeff(i,n-1); + m_matT.coeffRef(i,n-1) = q * z + p * m_matT.coeff(i,n); + m_matT.coeffRef(i,n) = q * m_matT.coeff(i,n) - p * z; } // Accumulate transformations for (int i = low; i <= high; ++i) { - z = m_eivec.coeff(i,n-1); - m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n); - m_eivec.coeffRef(i,n) = q * m_eivec.coeff(i,n) - p * z; + z = m_matU.coeff(i,n-1); + m_matU.coeffRef(i,n-1) = q * z + p * m_matU.coeff(i,n); + m_matU.coeffRef(i,n) = q * m_matU.coeff(i,n) - p * z; } } else // Complex pair { - m_eivalues.coeffRef(n-1) = Complex(x + p, z); - m_eivalues.coeffRef(n) = Complex(x + p, -z); + m_eivalues.coeffRef(n-1) = ComplexScalar(x + p, z); + m_eivalues.coeffRef(n) = ComplexScalar(x + p, -z); } n = n - 2; iter = 0; @@ -235,13 +234,13 @@ void RealSchur::hqr2() else // No convergence yet { // Form shift - x = matH.coeff(n,n); + x = m_matT.coeff(n,n); y = 0.0; w = 0.0; if (l < n) { - y = matH.coeff(n-1,n-1); - w = matH.coeff(n,n-1) * matH.coeff(n-1,n); + y = m_matT.coeff(n-1,n-1); + w = m_matT.coeff(n,n-1) * m_matT.coeff(n-1,n); } // Wilkinson's original ad hoc shift @@ -249,8 +248,8 @@ void RealSchur::hqr2() { exshift += x; for (int i = low; i <= n; ++i) - matH.coeffRef(i,i) -= x; - s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2)); + m_matT.coeffRef(i,i) -= x; + s = ei_abs(m_matT.coeff(n,n-1)) + ei_abs(m_matT.coeff(n-1,n-2)); x = y = Scalar(0.75) * s; w = Scalar(-0.4375) * s * s; } @@ -267,7 +266,7 @@ void RealSchur::hqr2() s = -s; s = Scalar(x - w / ((y - x) / 2.0 + s)); for (int i = low; i <= n; ++i) - matH.coeffRef(i,i) -= s; + m_matT.coeffRef(i,i) -= s; exshift += s; x = y = w = Scalar(0.964); } @@ -279,12 +278,12 @@ void RealSchur::hqr2() int m = n-2; while (m >= l) { - z = matH.coeff(m,m); + z = m_matT.coeff(m,m); r = x - z; s = y - z; - p = (r * s - w) / matH.coeff(m+1,m) + matH.coeff(m,m+1); - q = matH.coeff(m+1,m+1) - z - r - s; - r = matH.coeff(m+2,m+1); + p = (r * s - w) / m_matT.coeff(m+1,m) + m_matT.coeff(m,m+1); + q = m_matT.coeff(m+1,m+1) - z - r - s; + r = m_matT.coeff(m+2,m+1); s = ei_abs(p) + ei_abs(q) + ei_abs(r); p = p / s; q = q / s; @@ -292,9 +291,9 @@ void RealSchur::hqr2() if (m == l) { break; } - if (ei_abs(matH.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) < - eps * (ei_abs(p) * (ei_abs(matH.coeff(m-1,m-1)) + ei_abs(z) + - ei_abs(matH.coeff(m+1,m+1))))) + if (ei_abs(m_matT.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) < + NumTraits::epsilon() * (ei_abs(p) * (ei_abs(m_matT.coeff(m-1,m-1)) + ei_abs(z) + + ei_abs(m_matT.coeff(m+1,m+1))))) { break; } @@ -303,9 +302,9 @@ void RealSchur::hqr2() for (int i = m+2; i <= n; ++i) { - matH.coeffRef(i,i-2) = 0.0; + m_matT.coeffRef(i,i-2) = 0.0; if (i > m+2) - matH.coeffRef(i,i-3) = 0.0; + m_matT.coeffRef(i,i-3) = 0.0; } // Double QR step involving rows l:n and columns m:n @@ -313,9 +312,9 @@ void RealSchur::hqr2() { int notlast = (k != n-1); if (k != m) { - p = matH.coeff(k,k-1); - q = matH.coeff(k+1,k-1); - r = notlast ? matH.coeff(k+2,k-1) : Scalar(0); + p = m_matT.coeff(k,k-1); + q = m_matT.coeff(k+1,k-1); + r = notlast ? m_matT.coeff(k+2,k-1) : Scalar(0); x = ei_abs(p) + ei_abs(q) + ei_abs(r); if (x != 0.0) { @@ -336,9 +335,9 @@ void RealSchur::hqr2() if (s != 0) { if (k != m) - matH.coeffRef(k,k-1) = -s * x; + m_matT.coeffRef(k,k-1) = -s * x; else if (l != m) - matH.coeffRef(k,k-1) = -matH.coeff(k,k-1); + m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1); p = p + s; x = p / s; @@ -348,42 +347,42 @@ void RealSchur::hqr2() r = r / p; // Row modification - for (int j = k; j < nn; ++j) + for (int j = k; j < size; ++j) { - p = matH.coeff(k,j) + q * matH.coeff(k+1,j); + p = m_matT.coeff(k,j) + q * m_matT.coeff(k+1,j); if (notlast) { - p = p + r * matH.coeff(k+2,j); - matH.coeffRef(k+2,j) = matH.coeff(k+2,j) - p * z; + p = p + r * m_matT.coeff(k+2,j); + m_matT.coeffRef(k+2,j) = m_matT.coeff(k+2,j) - p * z; } - matH.coeffRef(k,j) = matH.coeff(k,j) - p * x; - matH.coeffRef(k+1,j) = matH.coeff(k+1,j) - p * y; + m_matT.coeffRef(k,j) = m_matT.coeff(k,j) - p * x; + m_matT.coeffRef(k+1,j) = m_matT.coeff(k+1,j) - p * y; } // Column modification for (int i = 0; i <= std::min(n,k+3); ++i) { - p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1); + p = x * m_matT.coeff(i,k) + y * m_matT.coeff(i,k+1); if (notlast) { - p = p + z * matH.coeff(i,k+2); - matH.coeffRef(i,k+2) = matH.coeff(i,k+2) - p * r; + p = p + z * m_matT.coeff(i,k+2); + m_matT.coeffRef(i,k+2) = m_matT.coeff(i,k+2) - p * r; } - matH.coeffRef(i,k) = matH.coeff(i,k) - p; - matH.coeffRef(i,k+1) = matH.coeff(i,k+1) - p * q; + m_matT.coeffRef(i,k) = m_matT.coeff(i,k) - p; + m_matT.coeffRef(i,k+1) = m_matT.coeff(i,k+1) - p * q; } // Accumulate transformations for (int i = low; i <= high; ++i) { - p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1); + p = x * m_matU.coeff(i,k) + y * m_matU.coeff(i,k+1); if (notlast) { - p = p + z * m_eivec.coeff(i,k+2); - m_eivec.coeffRef(i,k+2) = m_eivec.coeff(i,k+2) - p * r; + p = p + z * m_matU.coeff(i,k+2); + m_matU.coeffRef(i,k+2) = m_matU.coeff(i,k+2) - p * r; } - m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) - p; - m_eivec.coeffRef(i,k+1) = m_eivec.coeff(i,k+1) - p * q; + m_matU.coeffRef(i,k) = m_matU.coeff(i,k) - p; + m_matU.coeffRef(i,k+1) = m_matU.coeff(i,k+1) - p * q; } } // (s != 0) } // k loop