Fix a lot in MatrixPower.h
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 2dff280..f4f5b88 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -21,14 +21,12 @@ * * \brief Class for computing matrix powers. * - * \tparam MatrixType type of the base, expected to be an instantiation + * \tparam MatrixType type of the base, expected to be an instantiation * of the Matrix class template. - * \tparam ExponentType type of the exponent, a real scalar. - * \tparam PlainObject type of the multiplier. - * \tparam IsInteger used internally to select correct specialization. + * \tparam IsInteger used internally to select correct specialization. + * \tparam PlainObject type of the multiplier. */ -template <typename MatrixType, typename ExponentType, typename PlainObject = MatrixType, - int IsInteger = NumTraits<ExponentType>::IsInteger> +template <typename MatrixType, int IsInteger, typename PlainObject = MatrixType> class MatrixPower { private: @@ -93,7 +91,7 @@ void computeChainProduct(ResultType&); /** \brief Compute the cost of binary powering. */ - int computeCost(RealScalar); + static int computeCost(RealScalar); /** \brief Solve the linear system repetitively. */ template <typename ResultType> @@ -106,8 +104,8 @@ * \brief Split #m_p into integral part and fractional part. * * This method stores the integral part \f$ p_{\textrm int} \f$ into - * #m_pint and the fractional part \f$ p_{\textrm frac} \f$ into - * #m_pfrac, where #m_pfrac is in the interval \f$ (-1,1) \f$. To + * #m_pInt and the fractional part \f$ p_{\textrm frac} \f$ into + * #m_pFrac, where #m_pFrac is in the interval \f$ (-1,1) \f$. To * choose between the possibilities below, it considers the computation * of \f$ A^{p_1} \f$ and \f$ A^{p_2} \f$ and determines which of these * computations is the better conditioned. @@ -115,10 +113,10 @@ void getFractionalExponent(); /** \brief Compute atanh (inverse hyperbolic tangent) for \f$ y / x \f$. */ - ComplexScalar atanh2(const ComplexScalar& y, const ComplexScalar& x); + static ComplexScalar atanh2(const ComplexScalar& y, const ComplexScalar& x); /** \brief Compute power of 2x2 triangular matrix. */ - void compute2x2(const RealScalar& p); + void compute2x2(RealScalar p); /** * \brief Compute power of triangular matrices with size > 2. @@ -159,16 +157,16 @@ void computeTmp(RealScalar); const MatrixType& m_A; ///< \brief Reference to the matrix base. - const RealScalar& m_p; ///< \brief Reference to the real exponent. + const RealScalar m_p; ///< \brief The real 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. - RealScalar m_pint; ///< \brief Integer part of #m_p. - RealScalar m_pfrac; ///< \brief Fractional part of #m_p. + RealScalar m_pInt; ///< \brief Integral part of #m_p. + RealScalar m_pFrac; ///< \brief Fractional part of #m_p. ComplexMatrix m_T; ///< \brief Triangular part of Schur decomposition. ComplexMatrix m_U; ///< \brief Unitary part of Schur decomposition. - ComplexMatrix m_fT; ///< \brief #m_T to the power of #m_pfrac. + ComplexMatrix m_fT; ///< \brief #m_T to the power of #m_pFrac. ComplexArray m_logTdiag; ///< \brief Logarithm of the main diagonal of #m_T. }; @@ -176,8 +174,8 @@ * \internal \ingroup MatrixFunctions_Module * \brief Partial specialization for integral exponents. */ -template <typename MatrixType, typename IntExponent, typename PlainObject> -class MatrixPower<MatrixType, IntExponent, PlainObject, 1> +template <typename MatrixType, typename PlainObject> +class MatrixPower<MatrixType, 1, PlainObject> { public: /** @@ -187,7 +185,7 @@ * \param[in] p the exponent of the matrix power. * \param[in] b the multiplier. */ - MatrixPower(const MatrixType& A, const IntExponent& p, const PlainObject& b) : + MatrixPower(const MatrixType& A, int p, const PlainObject& b) : m_A(A), m_p(p), m_b(b), @@ -213,7 +211,7 @@ typedef typename MatrixType::Index Index; const MatrixType& m_A; ///< \brief Reference to the matrix base. - const IntExponent& m_p; ///< \brief Reference to the real exponent. + 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(). @@ -230,48 +228,51 @@ void computeChainProduct(ResultType& result); /** \brief Compute the cost of binary powering. */ - int computeCost(const IntExponent& p); + static int computeCost(int); /** \brief Solve the linear system repetitively. */ template <typename ResultType> - void partialPivLuSolve(ResultType&, IntExponent); + void partialPivLuSolve(ResultType&, int); }; /******* Specialized for real exponents *******/ -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> +template <typename MatrixType, int IsInteger, typename PlainObject> template <typename ResultType> -void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::compute(ResultType& result) +void MatrixPower<MatrixType,IsInteger,PlainObject>::compute(ResultType& result) { + using std::abs; using std::floor; using std::pow; - m_pint = floor(m_p); - m_pfrac = m_p - m_pint; + m_pInt = floor(m_p + RealScalar(0.5)); + m_pFrac = m_p - m_pInt; - if (m_pfrac == RealScalar(0)) + if (!m_pFrac) { computeIntPower(result); - else if (m_dimA == 1) + } else if (m_dimA == 1) result = pow(m_A(0,0), m_p) * m_b; else { computeSchurDecomposition(); getFractionalExponent(); computeIntPower(result); - if (m_dimA == 2) - compute2x2(m_pfrac); - else + if (m_dimA == 2) { + compute2x2(m_pFrac); + } else { computeBig(); + } computeTmp(Scalar()); - result *= m_tmp; + result = m_tmp * result; } } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> +template <typename MatrixType, int IsInteger, typename PlainObject> template <typename ResultType> -void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeIntPower(ResultType& result) +void MatrixPower<MatrixType,IsInteger,PlainObject>::computeIntPower(ResultType& result) { + MatrixType tmp; if (m_dimb > m_dimA) { - MatrixType tmp = MatrixType::Identity(m_A.rows(), m_A.cols()); + tmp = MatrixType::Identity(m_dimA, m_dimA); computeChainProduct(tmp); result = tmp * m_b; } else { @@ -280,18 +281,19 @@ } } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> +template <typename MatrixType, int IsInteger, typename PlainObject> template <typename ResultType> -void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeChainProduct(ResultType& result) +void MatrixPower<MatrixType,IsInteger,PlainObject>::computeChainProduct(ResultType& result) { + using std::abs; + using std::fmod; using std::frexp; using std::ldexp; - const bool pIsNegative = m_pint < RealScalar(0); - RealScalar p = pIsNegative? -m_pint: m_pint; + RealScalar p = abs(m_pInt); int cost = computeCost(p); - if (pIsNegative) { + if (m_pInt < RealScalar(0)) { if (p * m_dimb <= cost * m_dimA) { partialPivLuSolve(result, p); return; @@ -314,12 +316,13 @@ result = m_tmp * result; } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -int MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeCost(RealScalar p) +template <typename MatrixType, int IsInteger, typename PlainObject> +int MatrixPower<MatrixType,IsInteger,PlainObject>::computeCost(RealScalar p) { using std::frexp; using std::ldexp; int cost, tmp; + frexp(p, &cost); while (frexp(p, &tmp), tmp > 0) { p -= ldexp(RealScalar(0.5), tmp); @@ -328,61 +331,49 @@ return cost; } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> +template <typename MatrixType, int IsInteger, typename PlainObject> template <typename ResultType> -void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::partialPivLuSolve(ResultType& result, RealScalar p) +void MatrixPower<MatrixType,IsInteger,PlainObject>::partialPivLuSolve(ResultType& result, RealScalar p) { const PartialPivLU<MatrixType> Asolver(m_A); for (; p >= RealScalar(1); p--) result = Asolver.solve(result); } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeSchurDecomposition() +template <typename MatrixType, int IsInteger, typename PlainObject> +void MatrixPower<MatrixType,IsInteger,PlainObject>::computeSchurDecomposition() { const ComplexSchur<MatrixType> schurOfA(m_A); m_T = schurOfA.matrixT(); m_U = schurOfA.matrixU(); } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::getFractionalExponent() +template <typename MatrixType, int IsInteger, typename PlainObject> +void MatrixPower<MatrixType,IsInteger,PlainObject>::getFractionalExponent() { using std::pow; - typedef Array<RealScalar, Rows, 1, ColMajor, MaxRows> RealArray; + const ComplexArray Tdiag = m_T.diagonal(); - RealScalar maxAbsEival, minAbsEival, *begin, *end; - RealArray absTdiag; + const RealArray absTdiag = Tdiag.abs(); + const RealScalar maxAbsEival = absTdiag.maxCoeff(); + const RealScalar minAbsEival = absTdiag.minCoeff(); m_logTdiag = Tdiag.log(); - absTdiag = Tdiag.abs(); - maxAbsEival = minAbsEival = absTdiag[0]; - begin = absTdiag.data(); - end = begin + m_dimA; - - // This avoids traversing the array twice. - for (RealScalar *ptr = begin + 1; ptr < end; ptr++) { - if (*ptr > maxAbsEival) - maxAbsEival = *ptr; - else if (*ptr < minAbsEival) - minAbsEival = *ptr; - } - if (m_pfrac > RealScalar(0.5) && // This is just a shortcut. - m_pfrac > (RealScalar(1) - m_pfrac) * pow(maxAbsEival/minAbsEival, m_pfrac)) { - m_pfrac--; - m_pint++; + if (m_pFrac > RealScalar(0.5) && // This is just a shortcut. + m_pFrac > (RealScalar(1) - m_pFrac) * pow(maxAbsEival/minAbsEival, m_pFrac)) { + m_pFrac--; + m_pInt++; } } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> +template <typename MatrixType, int IsInteger, typename PlainObject> std::complex<typename MatrixType::RealScalar> -MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::atanh2(const ComplexScalar& y, const ComplexScalar& x) +MatrixPower<MatrixType,IsInteger,PlainObject>::atanh2(const ComplexScalar& y, const ComplexScalar& x) { using std::abs; using std::log; using std::sqrt; - const ComplexScalar z = y / x; if (abs(z) > sqrt(NumTraits<RealScalar>::epsilon())) @@ -391,8 +382,8 @@ return z + z*z*z / RealScalar(3); } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::compute2x2(const RealScalar& p) +template <typename MatrixType, int IsInteger, typename PlainObject> +void MatrixPower<MatrixType,IsInteger,PlainObject>::compute2x2(RealScalar p) { using std::abs; using std::ceil; @@ -407,7 +398,6 @@ ComplexScalar w; m_fT(0,0) = pow(m_T(0,0), p); - for (j = 1; j < m_dimA; j++) { i = j - 1; m_fT(j,j) = pow(m_T(j,j), p); @@ -426,8 +416,8 @@ } } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeBig() +template <typename MatrixType, int IsInteger, typename PlainObject> +void MatrixPower<MatrixType,IsInteger,PlainObject>::computeBig() { using std::ldexp; const int digits = std::numeric_limits<RealScalar>::digits; @@ -441,7 +431,7 @@ RealScalar normIminusT; while (true) { - IminusT = ComplexMatrix::Identity(m_A.rows(), m_A.cols()) - T; + IminusT = ComplexMatrix::Identity(m_dimA, m_dimA) - T; normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff(); if (normIminusT < maxNormForPade) { degree = getPadeDegree(normIminusT); @@ -457,14 +447,14 @@ computePade(degree, IminusT); for (; numberOfSquareRoots; numberOfSquareRoots--) { - compute2x2(ldexp(m_pfrac, -numberOfSquareRoots)); + compute2x2(ldexp(m_pFrac, -numberOfSquareRoots)); m_fT *= m_fT; } - compute2x2(m_pfrac); + compute2x2(m_pFrac); } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -inline int MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::getPadeDegree(float normIminusT) +template <typename MatrixType, int IsInteger, typename PlainObject> +inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(float normIminusT) { const float maxNormForPade[] = { 2.7996156e-1f /* degree = 3 */ , 4.3268868e-1f }; int degree = 3; @@ -474,8 +464,8 @@ return degree; } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -inline int MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::getPadeDegree(double normIminusT) +template <typename MatrixType, int IsInteger, typename PlainObject> +inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(double normIminusT) { const double maxNormForPade[] = { 1.882832775783710e-2 /* degree = 3 */ , 6.036100693089536e-2, 1.239372725584857e-1, 1.998030690604104e-1, 2.787629930861592e-1 }; @@ -486,8 +476,8 @@ return degree; } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -inline int MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::getPadeDegree(long double normIminusT) +template <typename MatrixType, int IsInteger, typename PlainObject> +inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(long double normIminusT) { #if LDBL_MANT_DIG == 53 const int maxPadeDegree = 7; @@ -519,45 +509,46 @@ break; return degree; } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computePade(const int& degree, const ComplexMatrix& IminusT) +template <typename MatrixType, int IsInteger, typename PlainObject> +void MatrixPower<MatrixType,IsInteger,PlainObject>::computePade(const int& degree, const ComplexMatrix& IminusT) { int i = degree << 1; m_fT = coeff(i) * IminusT; for (i--; i; i--) { - m_fT = (ComplexMatrix::Identity(m_A.rows(), m_A.cols()) + m_fT).template triangularView<Upper>() + m_fT = (ComplexMatrix::Identity(m_dimA, m_dimA) + m_fT).template triangularView<Upper>() .solve(coeff(i) * IminusT).eval(); } - m_fT += ComplexMatrix::Identity(m_A.rows(), m_A.cols()); + m_fT += ComplexMatrix::Identity(m_dimA, m_dimA); } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -inline typename MatrixType::RealScalar MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::coeff(const int& i) +template <typename MatrixType, int IsInteger, typename PlainObject> +inline typename MatrixType::RealScalar MatrixPower<MatrixType,IsInteger,PlainObject>::coeff(const int& i) { if (i == 1) - return -m_pfrac; + return -m_pFrac; else if (i & 1) - return (-m_pfrac - RealScalar(i >> 1)) / RealScalar(i << 1); + return (-m_pFrac - RealScalar(i >> 1)) / RealScalar(i << 1); else - return (m_pfrac - RealScalar(i >> 1)) / RealScalar(i-1 << 1); + return (m_pFrac - RealScalar(i >> 1)) / RealScalar((i - 1) << 1); } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeTmp(RealScalar) +template <typename MatrixType, int IsInteger, typename PlainObject> +void MatrixPower<MatrixType,IsInteger,PlainObject>::computeTmp(RealScalar) { m_tmp = (m_U * m_fT * m_U.adjoint()).real(); } -template <typename MatrixType, typename ExponentType, typename PlainObject, int IsInteger> -void MatrixPower<MatrixType,ExponentType,PlainObject,IsInteger>::computeTmp(ComplexScalar) +template <typename MatrixType, int IsInteger, typename PlainObject> +void MatrixPower<MatrixType,IsInteger,PlainObject>::computeTmp(ComplexScalar) { m_tmp = m_U * m_fT * m_U.adjoint(); } /******* Specialized for integral exponents *******/ -template <typename MatrixType, typename IntExponent, typename PlainObject> +template <typename MatrixType, typename PlainObject> template <typename ResultType> -void MatrixPower<MatrixType,IntExponent,PlainObject,1>::compute(ResultType& result) +void MatrixPower<MatrixType,1,PlainObject>::compute(ResultType& result) { + MatrixType tmp; if (m_dimb > m_dimA) { - MatrixType tmp = MatrixType::Identity(m_dimA, m_dimA); + tmp = MatrixType::Identity(m_dimA, m_dimA); computeChainProduct(tmp); result = tmp * m_b; } else { @@ -566,41 +557,43 @@ } } -template <typename MatrixType, typename IntExponent, typename PlainObject> -int MatrixPower<MatrixType,IntExponent,PlainObject,1>::computeCost(const IntExponent& p) +template <typename MatrixType, typename PlainObject> +int MatrixPower<MatrixType,1,PlainObject>::computeCost(int p) { - int cost = 0; - IntExponent tmp = p; - for (tmp = p >> 1; tmp; tmp >>= 1) + int cost = 0, tmp; + for (tmp = p; tmp; tmp >>= 1) cost++; - for (tmp = IntExponent(1); tmp <= p; tmp <<= 1) + for (tmp = 1; tmp <= p; tmp <<= 1) if (tmp & p) cost++; return cost; } -template <typename MatrixType, typename IntExponent, typename PlainObject> +template <typename MatrixType, typename PlainObject> template <typename ResultType> -void MatrixPower<MatrixType,IntExponent,PlainObject,1>::partialPivLuSolve(ResultType& result, IntExponent p) +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 IntExponent, typename PlainObject> +template <typename MatrixType, typename PlainObject> template <typename ResultType> -void MatrixPower<MatrixType,IntExponent,PlainObject,1>::computeChainProduct(ResultType& result) +void MatrixPower<MatrixType,1,PlainObject>::computeChainProduct(ResultType& result) { - const bool pIsNegative = m_p < IntExponent(0); - IntExponent p = pIsNegative? -m_p: m_p; - int cost = computeCost(p); + using std::abs; + int p = abs(m_p), cost = computeCost(p); - if (pIsNegative) { + if (m_p < 0) { if (p * m_dimb <= cost * m_dimA) { partialPivLuSolve(result, p); return; - } else { m_tmp = m_A.inverse(); } - } else { m_tmp = m_A; } + } else { + m_tmp = m_A.inverse(); + } + } else { + m_tmp = m_A; + } while (p * m_dimb > cost * m_dimA) { if (p & 1) { @@ -658,9 +651,10 @@ inline void evalTo(ResultType& result) const { typedef typename Derived::PlainObject PlainObject; + const int IsInteger = NumTraits<ExponentType>::IsInteger; const typename MatrixType::PlainObject Aevaluated = m_A.eval(); const PlainObject bevaluated = m_b.eval(); - MatrixPower<MatrixType, ExponentType, PlainObject> mp(Aevaluated, m_p, bevaluated); + MatrixPower<MatrixType, IsInteger, PlainObject> mp(Aevaluated, m_p, bevaluated); mp.compute(result); } @@ -726,9 +720,10 @@ 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 Identity = PlainObject::Identity(m_A.rows(), m_A.cols()); - MatrixPower<PlainObject, ExponentType> mp(Aevaluated, m_p, Identity); + MatrixPower<PlainObject, IsInteger> mp(Aevaluated, m_p, Identity); mp.compute(result); }
diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index f3ef571..80f65eb 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp
@@ -23,7 +23,7 @@ B << c, s, -s, c; C = A.pow(std::ldexp(angle, 1) / M_PI); - std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C, B) << "\n"; + std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C, B) << '\n'; VERIFY(C.isApprox(B, T(tol))); } } @@ -43,44 +43,117 @@ B << ch, ish, -ish, ch; C = A.pow(angle); - std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C, B) << "\n"; + std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C, B) << '\n'; VERIFY(C.isApprox(B, T(tol))); } } 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 = 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> void testExponentLaws(const MatrixType& m, double tol) { - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<Scalar>::Real RealScalar; + typedef typename MatrixType::RealScalar RealScalar; + MatrixType m1, m2, m3, m4, m5; + RealScalar x, y; - typename MatrixType::Index rows = m.rows(); - typename MatrixType::Index cols = m.cols(); - MatrixType m1, m1x, m1y, m2, m3; - RealScalar x = internal::random<RealScalar>(), y = internal::random<RealScalar>(); - double err[3]; - - for(int i = 0; i < g_repeat; i++) { + for (int i = 0; i < g_repeat; i++) { generateTestMatrix<MatrixType>::run(m1, m.rows()); - m1x = m1.pow(x); - m1y = m1.pow(y); + x = internal::random<RealScalar>(); + y = internal::random<RealScalar>(); + m2 = m1.pow(x); + m3 = m1.pow(y); - m2 = m1.pow(x + y); - m3 = m1x * m1y; - err[0] = relerr(m2, m3); - VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol))); + m4 = m1.pow(x + y); + m5 = m2 * m3; + std::cout << "testExponentLaws: error powerm = " << relerr(m4, m5); + VERIFY(m4.isApprox(m5, RealScalar(tol))); - m2 = m1.pow(x * y); - m3 = m1x.pow(y); - err[1] = relerr(m2, m3); - VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol))); + if (!NumTraits<typename MatrixType::Scalar>::IsComplex) { + m4 = m1.pow(x * y); + m5 = m2.pow(y); + std::cout << " " << relerr(m4, m5); + VERIFY(m4.isApprox(m5, RealScalar(tol))); + } - m2 = (std::abs(x) * m1).pow(y); - m3 = std::pow(std::abs(x), y) * m1y; - err[2] = relerr(m2, m3); - VERIFY(m2.isApprox(m3, static_cast<RealScalar>(tol))); + m4 = (std::abs(x) * m1).pow(y); + m5 = std::pow(std::abs(x), y) * m3; + std::cout << " " << relerr(m4, m5) << '\n'; + VERIFY(m4.isApprox(m5, RealScalar(tol))); + } +} - std::cout << "testExponentLaws: error powerm = " << err[0] << " " << err[1] << " " << err[2] << "\n"; +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; + + 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; + + v2 = m1.pow(pReal).eval() * v1; + v3 = m1.pow(pReal) * v1; + std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v3); + VERIFY(v2.isApprox(v3, RealScalar(tol))); + + v2 = m1.pow(pInt).eval() * v1; + v3 = m1.pow(pInt) * v1; + std::cout << " " << relerr(v2, v3) << '\n'; + VERIFY(v2.isApprox(v3, RealScalar(tol)) || v2 == v3); } } @@ -88,17 +161,36 @@ { CALL_SUBTEST_2(test2dRotation<double>(1e-13)); CALL_SUBTEST_1(test2dRotation<float>(2e-5)); // was 1e-5, relaxed for clang 2.8 / linux / x86-64 - CALL_SUBTEST_8(test2dRotation<long double>(1e-13)); + CALL_SUBTEST_9(test2dRotation<long double>(1e-13)); CALL_SUBTEST_2(test2dHyperbolicRotation<double>(1e-14)); CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5)); - CALL_SUBTEST_8(test2dHyperbolicRotation<long double>(1e-14)); + 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)); CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 1e-13)); CALL_SUBTEST_1(testExponentLaws(Matrix2f(), 1e-4)); CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4)); - CALL_SUBTEST_1(testExponentLaws(Matrix4f(), 1e-4)); + CALL_SUBTEST_8(testExponentLaws(Matrix4f(), 1e-4)); CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4)); - CALL_SUBTEST_9(testExponentLaws(Matrix<long double,Dynamic,Dynamic>(7,7), 1e-13)); + + CALL_SUBTEST_2(testMatrixVectorProduct(Matrix2d(), Vector2d(), 1e-13)); + CALL_SUBTEST_7(testMatrixVectorProduct(Matrix<double,3,3,RowMajor>(), Vector3d(), 1e-13)); + CALL_SUBTEST_3(testMatrixVectorProduct(Matrix4cd(), Vector4cd(), 1e-13)); + CALL_SUBTEST_4(testMatrixVectorProduct(MatrixXd(8,8), MatrixXd(8,2), 1e-13)); + CALL_SUBTEST_1(testMatrixVectorProduct(Matrix2f(), Vector2f(), 1e-4)); + 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)); }