Tidy up and write dox.
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 576c5d2..c00c148 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h
@@ -454,8 +454,7 @@ const MatrixFunctionReturnValue<Derived> sin() const; const MatrixSquareRootReturnValue<Derived> sqrt() const; const MatrixLogarithmReturnValue<Derived> log() const; - template <typename ExponentType> - const MatrixPowerReturnValue<Derived, ExponentType> pow(const ExponentType& p) const; + const MatrixPowerReturnValue<Derived> pow(RealScalar p) const; #ifdef EIGEN2_SUPPORT template<typename ProductDerived, typename Lhs, typename Rhs>
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index da79c9e..30d32e2 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -271,7 +271,7 @@ template<typename Derived> class MatrixFunctionReturnValue; template<typename Derived> class MatrixSquareRootReturnValue; template<typename Derived> class MatrixLogarithmReturnValue; -template<typename Derived, typename ExponentType> class MatrixPowerReturnValue; +template<typename Derived> class MatrixPowerReturnValue; namespace internal { template <typename Scalar>
diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index 041d3b7..002c1f7 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions
@@ -223,8 +223,7 @@ Compute the matrix raised to arbitrary real power. \code -template <typename ExponentType> -const MatrixPowerReturnValue<Derived, ExponentType> MatrixBase<Derived>::pow(const ExponentType& p) const +const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const \endcode \param[in] M base of the matrix power, should be a square matrix. @@ -247,6 +246,15 @@ The actual work is done by the MatrixPower class, which can compute \f$ M^p v \f$, where \p v is another matrix with the same rows as \p M. The matrix \p v is set to be the identity matrix by default. +Therefore, the expression <tt>M.pow(p) * v</tt> is specialized for +this. No temporary storage is created for the result. The code below +directly evaluates R-values into L-values without aliasing issue. Do +\b NOT try to \a optimize with noalias(). It won't compile. +\code +v = m.pow(p) * v; +m = m.pow(q); +// v2.noalias() = m.pow(p) * v1; Won't compile! +\endcode Details of the algorithm can be found in: Nicholas J. Higham and Lijing Lin, "A Schur-Padé algorithm for fractional powers of a
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 4c9039c..2a46d2c 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -23,10 +23,9 @@ * * \tparam MatrixType type of the base, expected to be an instantiation * of the Matrix class template. - * \tparam IsInteger used internally to select correct specialization. * \tparam PlainObject type of the multiplier. */ -template <typename MatrixType, int IsInteger, typename PlainObject = MatrixType> +template<typename MatrixType, typename PlainObject = MatrixType> class MatrixPower { private: @@ -65,11 +64,11 @@ * * \param[out] result \f$ A^p b \f$, as specified in the constructor. */ - template <typename ResultType> void compute(ResultType& result); + template<typename ResultType> void compute(ResultType& result); private: /** - * \brief Compute the matrix power. + * \brief Compute the matrix to integral power. * * If \p b is \em fatter than \p A, it computes \f$ A^{p_{\textrm int}} * \f$ first, and then multiplies it with \p b. Otherwise, @@ -77,7 +76,7 @@ * * \sa computeChainProduct(ResultType&); */ - template <typename ResultType> + template<typename ResultType> void computeIntPower(ResultType& result); /** @@ -87,14 +86,14 @@ * powering or matrix chain multiplication or solving the linear system * repetitively, according to which algorithm costs less. */ - template <typename ResultType> + template<typename ResultType> void computeChainProduct(ResultType&); /** \brief Compute the cost of binary powering. */ static int computeCost(RealScalar); /** \brief Solve the linear system repetitively. */ - template <typename ResultType> + template<typename ResultType> void partialPivLuSolve(ResultType&, RealScalar); /** \brief Compute Schur decomposition of #m_A. */ @@ -170,76 +169,9 @@ ComplexArray m_logTdiag; ///< \brief Logarithm of the main diagonal of #m_T. }; -/** - * \internal \ingroup MatrixFunctions_Module - * \brief Partial specialization for integral exponents. - */ -template <typename MatrixType, typename PlainObject> -class MatrixPower<MatrixType, 1, PlainObject> -{ - public: - /** - * \brief Constructor. - * - * \param[in] A the base of the matrix power. - * \param[in] p the exponent of the matrix power. - * \param[in] b the multiplier. - */ - MatrixPower(const MatrixType& A, int p, const PlainObject& b) : - m_A(A), - m_p(p), - m_b(b), - m_dimA(A.cols()), - m_dimb(b.cols()) - { /* empty body */ } - - /** - * \brief Compute the matrix power. - * - * If \p b is \em fatter than \p A, it computes \f$ A^p \f$ first, and - * then multiplies it with \p b. Otherwise, #computeChainProduct - * optimizes the expression. - * - * \param[out] result \f$ A^p b \f$, as specified in the constructor. - * - * \sa computeChainProduct(ResultType&); - */ - template <typename ResultType> - void compute(ResultType& result); - - private: - typedef typename MatrixType::Index Index; - - const MatrixType& m_A; ///< \brief Reference to the matrix base. - const int m_p; ///< \brief The integral exponent. - const PlainObject& m_b; ///< \brief Reference to the multiplier. - const Index m_dimA; ///< \brief The dimension of #m_A, equivalent to %m_A.cols(). - const Index m_dimb; ///< \brief The dimension of #m_b, equivalent to %m_b.cols(). - MatrixType m_tmp; ///< \brief Used for temporary storage. - - /** - * \brief Convert matrix power into chain product. - * - * This optimizes the matrix expression. It automatically chooses binary - * powering or matrix chain multiplication or solving the linear system - * repetitively, according to which algorithm costs less. - */ - template <typename ResultType> - void computeChainProduct(ResultType& result); - - /** \brief Compute the cost of binary powering. */ - static int computeCost(int); - - /** \brief Solve the linear system repetitively. */ - template <typename ResultType> - void partialPivLuSolve(ResultType&, int); -}; - -/******* Specialized for real exponents *******/ - -template <typename MatrixType, int IsInteger, typename PlainObject> -template <typename ResultType> -void MatrixPower<MatrixType,IsInteger,PlainObject>::compute(ResultType& result) +template<typename MatrixType, typename PlainObject> +template<typename ResultType> +void MatrixPower<MatrixType,PlainObject>::compute(ResultType& result) { using std::floor; using std::pow; @@ -255,19 +187,18 @@ computeSchurDecomposition(); getFractionalExponent(); computeIntPower(result); - if (m_dimA == 2) { + if (m_dimA == 2) compute2x2(m_pFrac); - } else { + else computeBig(); - } computeTmp(Scalar()); result = m_tmp * result; } } -template <typename MatrixType, int IsInteger, typename PlainObject> -template <typename ResultType> -void MatrixPower<MatrixType,IsInteger,PlainObject>::computeIntPower(ResultType& result) +template<typename MatrixType, typename PlainObject> +template<typename ResultType> +void MatrixPower<MatrixType,PlainObject>::computeIntPower(ResultType& result) { MatrixType tmp; if (m_dimb > m_dimA) { @@ -280,9 +211,9 @@ } } -template <typename MatrixType, int IsInteger, typename PlainObject> -template <typename ResultType> -void MatrixPower<MatrixType,IsInteger,PlainObject>::computeChainProduct(ResultType& result) +template<typename MatrixType, typename PlainObject> +template<typename ResultType> +void MatrixPower<MatrixType,PlainObject>::computeChainProduct(ResultType& result) { using std::abs; using std::fmod; @@ -314,8 +245,8 @@ result = m_tmp * result; } -template <typename MatrixType, int IsInteger, typename PlainObject> -int MatrixPower<MatrixType,IsInteger,PlainObject>::computeCost(RealScalar p) +template<typename MatrixType, typename PlainObject> +int MatrixPower<MatrixType,PlainObject>::computeCost(RealScalar p) { using std::frexp; using std::ldexp; @@ -329,25 +260,25 @@ return cost; } -template <typename MatrixType, int IsInteger, typename PlainObject> -template <typename ResultType> -void MatrixPower<MatrixType,IsInteger,PlainObject>::partialPivLuSolve(ResultType& result, RealScalar p) +template<typename MatrixType, typename PlainObject> +template<typename ResultType> +void MatrixPower<MatrixType,PlainObject>::partialPivLuSolve(ResultType& result, RealScalar p) { const PartialPivLU<MatrixType> Asolver(m_A); for (; p >= RealScalar(1); p--) result = Asolver.solve(result); } -template <typename MatrixType, int IsInteger, typename PlainObject> -void MatrixPower<MatrixType,IsInteger,PlainObject>::computeSchurDecomposition() +template<typename MatrixType, typename PlainObject> +void MatrixPower<MatrixType,PlainObject>::computeSchurDecomposition() { const ComplexSchur<MatrixType> schurOfA(m_A); m_T = schurOfA.matrixT(); m_U = schurOfA.matrixU(); } -template <typename MatrixType, int IsInteger, typename PlainObject> -void MatrixPower<MatrixType,IsInteger,PlainObject>::getFractionalExponent() +template<typename MatrixType, typename PlainObject> +void MatrixPower<MatrixType,PlainObject>::getFractionalExponent() { using std::pow; typedef Array<RealScalar, Rows, 1, ColMajor, MaxRows> RealArray; @@ -365,9 +296,9 @@ } } -template <typename MatrixType, int IsInteger, typename PlainObject> +template<typename MatrixType, typename PlainObject> std::complex<typename MatrixType::RealScalar> -MatrixPower<MatrixType,IsInteger,PlainObject>::atanh2(const ComplexScalar& y, const ComplexScalar& x) +MatrixPower<MatrixType,PlainObject>::atanh2(const ComplexScalar& y, const ComplexScalar& x) { using std::abs; using std::log; @@ -380,8 +311,8 @@ return z + z*z*z / RealScalar(3); } -template <typename MatrixType, int IsInteger, typename PlainObject> -void MatrixPower<MatrixType,IsInteger,PlainObject>::compute2x2(RealScalar p) +template<typename MatrixType, typename PlainObject> +void MatrixPower<MatrixType,PlainObject>::compute2x2(RealScalar p) { using std::abs; using std::ceil; @@ -413,8 +344,8 @@ } } -template <typename MatrixType, int IsInteger, typename PlainObject> -void MatrixPower<MatrixType,IsInteger,PlainObject>::computeBig() +template<typename MatrixType, typename PlainObject> +void MatrixPower<MatrixType,PlainObject>::computeBig() { using std::ldexp; const int digits = std::numeric_limits<RealScalar>::digits; @@ -450,8 +381,8 @@ compute2x2(m_pFrac); } -template <typename MatrixType, int IsInteger, typename PlainObject> -inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(float normIminusT) +template<typename MatrixType, typename PlainObject> +inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(float normIminusT) { const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f }; int degree = 3; @@ -461,8 +392,8 @@ return degree; } -template <typename MatrixType, int IsInteger, typename PlainObject> -inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(double normIminusT) +template<typename MatrixType, typename PlainObject> +inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(double normIminusT) { const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1, 1.999045567181744e-1, 2.789358995219730e-1 }; @@ -473,8 +404,8 @@ return degree; } -template <typename MatrixType, int IsInteger, typename PlainObject> -inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(long double normIminusT) +template<typename MatrixType, typename PlainObject> +inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(long double normIminusT) { #if LDBL_MANT_DIG == 53 const int maxPadeDegree = 7; @@ -506,8 +437,8 @@ break; return degree; } -template <typename MatrixType, int IsInteger, typename PlainObject> -void MatrixPower<MatrixType,IsInteger,PlainObject>::computePade(const int& degree, const ComplexMatrix& IminusT) +template<typename MatrixType, typename PlainObject> +void MatrixPower<MatrixType,PlainObject>::computePade(const int& degree, const ComplexMatrix& IminusT) { int i = degree << 1; m_fT = coeff(i) * IminusT; @@ -518,8 +449,8 @@ m_fT += ComplexMatrix::Identity(m_dimA, m_dimA); } -template <typename MatrixType, int IsInteger, typename PlainObject> -inline typename MatrixType::RealScalar MatrixPower<MatrixType,IsInteger,PlainObject>::coeff(const int& i) +template<typename MatrixType, typename PlainObject> +inline typename MatrixType::RealScalar MatrixPower<MatrixType,PlainObject>::coeff(const int& i) { if (i == 1) return -m_pFrac; @@ -529,90 +460,79 @@ return (m_pFrac - RealScalar(i >> 1)) / RealScalar((i - 1) << 1); } -template <typename MatrixType, int IsInteger, typename PlainObject> -void MatrixPower<MatrixType,IsInteger,PlainObject>::computeTmp(RealScalar) +template<typename MatrixType, typename PlainObject> +void MatrixPower<MatrixType,PlainObject>::computeTmp(RealScalar) { m_tmp = (m_U * m_fT * m_U.adjoint()).real(); } -template <typename MatrixType, int IsInteger, typename PlainObject> -void MatrixPower<MatrixType,IsInteger,PlainObject>::computeTmp(ComplexScalar) +template<typename MatrixType, typename PlainObject> +void MatrixPower<MatrixType,PlainObject>::computeTmp(ComplexScalar) { m_tmp = m_U * m_fT * m_U.adjoint(); } -/******* Specialized for integral exponents *******/ - -template <typename MatrixType, typename PlainObject> -template <typename ResultType> -void MatrixPower<MatrixType,1,PlainObject>::compute(ResultType& result) +/** + * \ingroup MatrixFunctions_Module + * + * \brief Proxy for the matrix power multiplied by other matrix. + * + * \tparam Derived type of the base, a matrix (expression). + * \tparam RhsDerived type of the multiplier. + * + * This class holds the arguments to the matrix power until it is + * assigned or evaluated for some other reason (so the argument + * should not be changed in the meantime). It is the return type of + * MatrixPowerReturnValue::operator*() and related functions and most + * of the time this is the only way it is used. + */ +template<typename Derived, typename RhsDerived> +class MatrixPowerProductReturnValue : public ReturnByValue<MatrixPowerProductReturnValue<Derived,RhsDerived> > { - MatrixType tmp; - if (m_dimb > m_dimA) { - tmp = MatrixType::Identity(m_dimA, m_dimA); - computeChainProduct(tmp); - result.noalias() = tmp * m_b; - } else { - result = m_b; - computeChainProduct(result); - } -} + private: + typedef typename Derived::PlainObject MatrixType; + typedef typename RhsDerived::PlainObject PlainObject; + typedef typename RhsDerived::RealScalar RealScalar; + typedef typename RhsDerived::Index Index; -template <typename MatrixType, typename PlainObject> -int MatrixPower<MatrixType,1,PlainObject>::computeCost(int p) -{ - int cost = 0, tmp; - for (tmp = p; tmp; tmp >>= 1) - cost++; - for (tmp = 1; tmp <= p; tmp <<= 1) - if (tmp & p) cost++; - return cost; -} + public: + /** + * \brief Constructor. + * + * \param[in] A %Matrix (expression), the base of the matrix power. + * \param[in] p scalar, the exponent of the matrix power. + * \prarm[in] b %Matrix (expression), the multiplier. + */ + MatrixPowerProductReturnValue(const Derived& A, RealScalar p, const RhsDerived& b) + : m_A(A), m_p(p), m_b(b) { } -template <typename MatrixType, typename PlainObject> -template <typename ResultType> -void MatrixPower<MatrixType,1,PlainObject>::partialPivLuSolve(ResultType& result, int p) -{ - const PartialPivLU<MatrixType> Asolver(m_A); - for(; p; p--) - result = Asolver.solve(result); -} - -template <typename MatrixType, typename PlainObject> -template <typename ResultType> -void MatrixPower<MatrixType,1,PlainObject>::computeChainProduct(ResultType& result) -{ - using std::abs; - int p = abs(m_p), cost = computeCost(p); - - if (m_p < 0) { - if (p * m_dimb <= cost * m_dimA) { - partialPivLuSolve(result, p); - return; - } else { - m_tmp = m_A.inverse(); + /** + * \brief Compute the expression. + * + * \param[out] result \f$ A^p b \f$ where \p A, \p p and \p bare as + * in the constructor. + */ + template<typename ResultType> + inline void evalTo(ResultType& result) const + { + const MatrixType A = m_A; + const PlainObject b = m_b; + MatrixPower<MatrixType, PlainObject> mp(A, m_p, b); + mp.compute(result); } - } else { - m_tmp = m_A; - } - - while (p * m_dimb > cost * m_dimA) { - if (p & 1) { - result = m_tmp * result; - cost--; - } - m_tmp *= m_tmp; - cost--; - p >>= 1; - } - for (; p; p--) - result = m_tmp * result; -} + Index rows() const { return m_b.rows(); } + Index cols() const { return m_b.cols(); } + + private: + const Derived& m_A; + const RealScalar m_p; + const RhsDerived& m_b; + MatrixPowerProductReturnValue& operator=(const MatrixPowerProductReturnValue&); +}; /** * \ingroup MatrixFunctions_Module * * \brief Proxy for the matrix power of some matrix (expression). * - * \tparam Derived type of the base, a matrix (expression). - * \tparam ExponentType type of the exponent, a scalar. + * \tparam Derived type of the base, a matrix (expression). * * This class holds the arguments to the matrix power until it is * assigned or evaluated for some other reason (so the argument @@ -620,19 +540,21 @@ * MatrixBase::pow() and related functions and most of the * time this is the only way it is used. */ -template<typename Derived, typename ExponentType> class MatrixPowerReturnValue -: public ReturnByValue<MatrixPowerReturnValue<Derived, ExponentType> > +template<typename Derived> +class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Derived> > { - public: + private: + typedef typename Derived::RealScalar RealScalar; typedef typename Derived::Index Index; + public: /** * \brief Constructor. * * \param[in] A %Matrix (expression), the base of the matrix power. * \param[in] p scalar, the exponent of the matrix power. */ - MatrixPowerReturnValue(const Derived& A, const ExponentType& p) + MatrixPowerReturnValue(const Derived& A, RealScalar p) : m_A(A), m_p(p) { } /** @@ -641,22 +563,19 @@ * The %MatrixPower class can optimize \f$ A^p b \f$ computing, and * this method provides an elegant way to call it. * - * \param[in] b %Matrix (expression), the multiplier. - * \param[out] result \f$ A^p b \f$ where \p A and \p p are as in - * the constructor. + * Unlike general matrix-matrix / matrix-vector product, this does + * \b NOT produce a temporary storage for the result. Therefore, + * the code below is \a already optimal: + * \code + * v = A.pow(p) * b; + * // v.noalias() = A.pow(p) * b; Won't compile! + * \endcode + * + * \param[in] b %Matrix (expression), the multiplier. */ - template <typename OtherDerived> - const typename OtherDerived::PlainObject operator*(const MatrixBase<OtherDerived>& b) const - { - typedef typename OtherDerived::PlainObject PlainObject; - const int IsInteger = NumTraits<ExponentType>::IsInteger; - const typename Derived::PlainObject Aevaluated = m_A.eval(); - const PlainObject bevaluated = b.eval(); - PlainObject result; - MatrixPower<Derived, IsInteger, PlainObject> mp(Aevaluated, m_p, bevaluated); - mp.compute(result); - return result; - } + template<typename RhsDerived> + const MatrixPowerProductReturnValue<Derived,RhsDerived> operator*(const MatrixBase<RhsDerived>& b) const + { return MatrixPowerProductReturnValue<Derived,RhsDerived>(m_A, m_p, b.derived()); } /** * \brief Compute the matrix power. @@ -664,14 +583,13 @@ * \param[out] result \f$ A^p \f$ where \p A and \p p are as in the * constructor. */ - template <typename ResultType> + template<typename ResultType> inline void evalTo(ResultType& result) const { typedef typename Derived::PlainObject PlainObject; - const int IsInteger = NumTraits<ExponentType>::IsInteger; - const PlainObject Aevaluated = m_A.eval(); + const PlainObject A = m_A; const PlainObject Identity = PlainObject::Identity(m_A.rows(), m_A.cols()); - MatrixPower<PlainObject, IsInteger> mp(Aevaluated, m_p, Identity); + MatrixPower<PlainObject> mp(A, m_p, Identity); mp.compute(result); } @@ -680,25 +598,29 @@ private: const Derived& m_A; - const ExponentType& m_p; - + const RealScalar m_p; MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&); }; namespace internal { - template<typename Derived, typename ExponentType> - struct traits<MatrixPowerReturnValue<Derived, ExponentType> > + template<typename Derived> + struct traits<MatrixPowerReturnValue<Derived> > { typedef typename Derived::PlainObject ReturnType; }; + + template<typename Derived, typename RhsDerived> + struct traits<MatrixPowerProductReturnValue<Derived,RhsDerived> > + { + typedef typename RhsDerived::PlainObject ReturnType; + }; } -template <typename Derived> -template <typename ExponentType> -const MatrixPowerReturnValue<Derived, ExponentType> MatrixBase<Derived>::pow(const ExponentType& p) const +template<typename Derived> +const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const { eigen_assert(rows() == cols()); - return MatrixPowerReturnValue<Derived, ExponentType>(derived(), p); + return MatrixPowerReturnValue<Derived>(derived(), p); } } // end namespace Eigen
diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index 3c0e4f3..d891641 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp
@@ -9,7 +9,7 @@ #include "matrix_functions.h" -template <typename T> +template<typename T> void test2dRotation(double tol) { Matrix<T,2,2> A, B, C; @@ -28,7 +28,7 @@ } } -template <typename T> +template<typename T> void test2dHyperbolicRotation(double tol) { Matrix<std::complex<T>,2,2> A, B, C; @@ -48,55 +48,7 @@ } } -template <typename MatrixType> -void testIntPowers(const MatrixType& m, double tol) -{ - typedef typename MatrixType::RealScalar RealScalar; - const MatrixType m1 = MatrixType::Random(m.rows(), m.cols()); - const MatrixType identity = MatrixType::Identity(m.rows(), m.cols()); - const PartialPivLU<MatrixType> solver(m1); - MatrixType m2, m3, m4; - - m3 = m1.pow(0); - m4 = m1.pow(0.); - std::cout << "testIntPower: i = 0 error powerm = " << relerr(identity, m3) << " " << relerr(identity, m4) << '\n'; - VERIFY(identity == m3 && identity == m4); - - m3 = m1.pow(1); - m4 = m1.pow(1.); - std::cout << "testIntPower: i = 1 error powerm = " << relerr(m1, m3) << " " << relerr(m1, m4) << '\n'; - VERIFY(m1 == m3 && m1 == m4); - - m2.noalias() = m1 * m1; - m3 = m1.pow(2); - m4 = m1.pow(2.); - std::cout << "testIntPower: i = 2 error powerm = " << relerr(m2, m3) << " " << relerr(m2, m4) << '\n'; - VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol))); - - for (int i = 3; i <= 20; i++) { - m2 *= m1; - m3 = m1.pow(i); - m4 = m1.pow(RealScalar(i)); - std::cout << "testIntPower: i = " << i << " error powerm = " << relerr(m2, m3) << " " << relerr (m2, m4) << '\n'; - VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol))); - } - - m2 = solver.inverse(); - m3 = m1.pow(-1); - m4 = m1.pow(-1.); - std::cout << "testIntPower: i = -1 error powerm = " << relerr(m2, m3) << " " << relerr (m2, m4) << '\n'; - VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol))); - - for (int i = -2; i >= -20; i--) { - m2 = solver.solve(m2); - m3 = m1.pow(i); - m4 = m1.pow(RealScalar(i)); - std::cout << "testIntPower: i = " << i << " error powerm = " << relerr(m2, m3) << " " << relerr (m2, m4) << '\n'; - VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol))); - } -} - -template <typename MatrixType> +template<typename MatrixType> void testExponentLaws(const MatrixType& m, double tol) { typedef typename MatrixType::RealScalar RealScalar; @@ -129,31 +81,40 @@ } } -template <typename MatrixType, typename VectorType> +template<typename MatrixType, typename VectorType> void testMatrixVectorProduct(const MatrixType& m, const VectorType& v, double tol) { typedef typename MatrixType::RealScalar RealScalar; MatrixType m1; VectorType v1, v2, v3; - RealScalar pReal; - signed char pInt; + RealScalar p; for (int i = 0; i < g_repeat; i++) { generateTestMatrix<MatrixType>::run(m1, m.rows()); v1 = VectorType::Random(v.rows(), v.cols()); - pReal = internal::random<RealScalar>(); - pInt = rand(); - pInt >>= 2; + p = internal::random<RealScalar>(); - v2.noalias() = m1.pow(pReal).eval() * v1; - v3.noalias() = m1.pow(pReal) * v1; - std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v3); - VERIFY(v2.isApprox(v3, RealScalar(tol))); + v2.noalias() = m1.pow(p).eval() * v1; + v1 = m1.pow(p) * v1; + std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v1) << '\n'; + VERIFY(v2.isApprox(v1, RealScalar(tol))); + } +} - v2.noalias() = m1.pow(pInt).eval() * v1; - v3.noalias() = m1.pow(pInt) * v1; - std::cout << " " << relerr(v2, v3) << '\n'; - VERIFY(v2.isApprox(v3, RealScalar(tol)) || v2 == v3); +template<typename MatrixType> +void testAliasing(const MatrixType& m) +{ + typedef typename MatrixType::RealScalar RealScalar; + MatrixType m1, m2; + RealScalar p; + + for (int i = 0; i < g_repeat; i++) { + generateTestMatrix<MatrixType>::run(m1, m.rows()); + p = internal::random<RealScalar>(); + + m2 = m1.pow(p); + m1 = m1.pow(p); + VERIFY(m1 == m2); } } @@ -166,15 +127,6 @@ CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5)); CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14)); - CALL_SUBTEST_2(testIntPowers(Matrix2d(), 1e-13)); - CALL_SUBTEST_7(testIntPowers(Matrix<double,3,3,RowMajor>(), 1e-13)); - CALL_SUBTEST_3(testIntPowers(Matrix4cd(), 1e-13)); - CALL_SUBTEST_4(testIntPowers(MatrixXd(8,8), 1e-13)); - CALL_SUBTEST_1(testIntPowers(Matrix2f(), 1e-4)); - CALL_SUBTEST_5(testIntPowers(Matrix3cf(), 1e-4)); - CALL_SUBTEST_8(testIntPowers(Matrix4f(), 1e-4)); - CALL_SUBTEST_6(testIntPowers(MatrixXf(8,8), 1e-4)); - CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13)); CALL_SUBTEST_7(testExponentLaws(Matrix<double,3,3,RowMajor>(), 1e-13)); CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13)); @@ -192,5 +144,11 @@ CALL_SUBTEST_5(testMatrixVectorProduct(Matrix3cf(), Vector3cf(), 1e-4)); CALL_SUBTEST_8(testMatrixVectorProduct(Matrix4f(), Vector4f(), 1e-4)); CALL_SUBTEST_6(testMatrixVectorProduct(MatrixXf(8,8), VectorXf(8), 1e-4)); - CALL_SUBTEST_10(testMatrixVectorProduct(Matrix<long double,Dynamic,Dynamic>(7,7), Matrix<long double,7,9>(), 1e-13)); + CALL_SUBTEST_9(testMatrixVectorProduct(Matrix<long double,Dynamic,Dynamic>(7,7), Matrix<long double,7,9>(), 1e-13)); + + CALL_SUBTEST_7(testAliasing(Matrix<double,3,3,RowMajor>())); + CALL_SUBTEST_3(testAliasing(Matrix4cd())); + CALL_SUBTEST_4(testAliasing(MatrixXd(8,8))); + CALL_SUBTEST_5(testAliasing(Matrix3cf())); + CALL_SUBTEST_6(testAliasing(MatrixXf(8,8))); }