From 4be172d84f7033b24ac9fed3ffd3e2091af25dfd Mon Sep 17 00:00:00 2001 From: jdh8 Date: Wed, 15 Aug 2012 00:34:20 +0800 Subject: [PATCH] matrix power: MatrixBase::pow(RealScalar) and MatrixBase::pow(T) where T is integral type --- Eigen/src/Core/MatrixBase.h | 2 + Eigen/src/Core/util/ForwardDeclarations.h | 1 + unsupported/Eigen/MatrixFunctions | 2 +- .../Eigen/src/MatrixFunctions/MatrixPower.h | 953 ++++++++++++++++++ 4 files changed, 957 insertions(+), 1 deletion(-) create mode 100644 unsupported/Eigen/src/MatrixFunctions/MatrixPower.h diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index bb569238f..576c5d25b 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -454,6 +454,8 @@ template class MatrixBase const MatrixFunctionReturnValue sin() const; const MatrixSquareRootReturnValue sqrt() const; const MatrixLogarithmReturnValue log() const; + template + const MatrixPowerReturnValue pow(const ExponentType& p) const; #ifdef EIGEN2_SUPPORT template diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 3afcc5ddb..da79c9e0f 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -271,6 +271,7 @@ template struct MatrixExponentialReturnValue; template class MatrixFunctionReturnValue; template class MatrixSquareRootReturnValue; template class MatrixLogarithmReturnValue; +template class MatrixPowerReturnValue; namespace internal { template diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index 56ab71cd3..9917efeed 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -57,7 +57,7 @@ #include "src/MatrixFunctions/MatrixFunction.h" #include "src/MatrixFunctions/MatrixSquareRoot.h" #include "src/MatrixFunctions/MatrixLogarithm.h" - +#include "src/MatrixFunctions/MatrixPower.h" /** diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h new file mode 100644 index 000000000..69c4000cf --- /dev/null +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -0,0 +1,953 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Chen-Pang He +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_MATRIX_POWER +#define EIGEN_MATRIX_POWER + +#ifndef M_PI +#define M_PI 3.141592653589793238462643383279503L +#endif + +namespace Eigen { + +/** + * \ingroup MatrixFunctions_Module + * + * \brief Class for computing matrix powers. + * + * \tparam MatrixType type of the base, expected to be an instantiation + * of the Matrix class template. + * \tparam RealScalar type of the exponent, a real scalar. + * \tparam PlainObject type of the multiplier. + * \tparam IsInteger used internally to select correct specialization. + */ +template ::IsInteger> +class MatrixPower +{ + 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, const RealScalar& 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. + * + * \param[out] result \f$ A^p b \f$, as specified in the constructor. + */ + template void compute(ResultType& result); + + private: + typedef internal::traits Traits; + static const int Rows = Traits::RowsAtCompileTime; + static const int Cols = Traits::ColsAtCompileTime; + static const int Options = Traits::Options; + static const int MaxRows = Traits::MaxRowsAtCompileTime; + static const int MaxCols = Traits::MaxColsAtCompileTime; + + typedef typename MatrixType::Scalar Scalar; + typedef std::complex ComplexScalar; + typedef typename MatrixType::Index Index; + typedef Matrix ComplexMatrix; + typedef Array ComplexArray; + + /** + * \brief Compute the matrix power. + * + * If \c b is \em fatter than \c A, it computes \f$ A^{p_{\textrm int}} + * \f$ first, and then multiplies it with \c b. Otherwise, + * #computeChainProduct optimizes the expression. + * + * \sa computeChainProduct(ResultType&); + */ + template + void computeIntPower(ResultType& result); + + /** + * \brief Convert integral power of the matrix 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 + void computeChainProduct(ResultType&); + + /** \brief Compute the cost of binary powering. */ + int computeCost(RealScalar); + + /** \brief Solve the linear system repetitively. */ + template + void partialPivLuSolve(RealScalar, ResultType&); + + /** \brief Compute Schur decomposition of #m_A. */ + void computeSchurDecomposition(); + + /** + * \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 + * 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. + */ + void getFractionalExponent(); + + /** \brief Compute atanh (inverse hyperbolic tangent). */ + ComplexScalar atanh(const ComplexScalar& x); + + /** \brief Compute power of 2x2 triangular matrix. */ + void compute2x2(const RealScalar& p); + + /** + * \brief Compute power of triangular matrices with size > 2. + * \details This uses a Schur-Padé algorithm. + */ + void computeBig(); + + /** \brief Get suitable degree for Pade approximation. (specialized for \c RealScalar = \c double) */ + inline int getPadeDegree(double); +/* TODO + * inline int getPadeDegree(float); + * + * inline int getPadeDegree(long double); + */ + /** \brief Compute Padé approximation to matrix fractional power. */ + void computePade(int degree, const ComplexMatrix& IminusT); + + /** \brief Get a certain coefficient of the Padé approximation. */ + inline RealScalar coeff(int degree); + + /** + * \brief Store the fractional power into #m_tmp. + * + * This intended for complex matrices. + */ + void computeTmp(ComplexScalar); + + /** + * \brief Store the fractional power into #m_tmp. + * + * This is intended for real matrices. It takes the real part of + * \f$ U T^{p_{\textrm frac}} U^H \f$. + * + * \sa computeTmp(ComplexScalar); + */ + void computeTmp(RealScalar); + + const MatrixType& m_A; ///< \brief Reference to the matrix base. + const RealScalar& m_p; ///< \brief Reference to 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. + 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. + ComplexArray m_logTdiag; ///< \brief Logarithm of the main diagonal of #m_T. +}; + +/** + * \internal \ingroup MatrixFunctions_Module + * \brief Partial specialization for integral exponents. + */ +template +class MatrixPower +{ + 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, const IntExponent& 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 \c b is \em fatter than \c A, it computes \f$ A^p \f$ first, and + * then multiplies it with \c b. Otherwise, #computeChainProduct + * optimizes the expression. + * + * \param[out] result \f$ A^p b \f$, as specified in the constructor. + * + * \sa computeChainProduct(ResultType&); + */ + template void compute(ResultType& result); + + private: + 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 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 void computeChainProduct(ResultType& result); + + /** \brief Compute the cost of binary powering. */ + int computeCost(const IntExponent& p); + + /** \brief Solve the linear system repetitively. */ + template + void partialPivLuSolve(IntExponent p, ResultType& result); +}; + +/** + * \internal \ingroup MatrixFunctions_Module + * \brief Partial specialization for complex matrices raised to complex exponents. + */ +template +class MatrixPower, PlainObject, IsInteger> +{ + private: + typedef internal::traits Traits; + static const int Rows = Traits::RowsAtCompileTime; + static const int Cols = Traits::ColsAtCompileTime; + static const int Options = Traits::Options; + static const int MaxRows = Traits::MaxRowsAtCompileTime; + static const int MaxCols = Traits::MaxColsAtCompileTime; + + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::Index Index; + typedef Array ArrayType; + + 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, const Scalar& p, const PlainObject& b) : + m_A(A), + m_p(p), + m_b(b), + m_dimA(A.cols()), + m_dimb(b.cols()) + { EIGEN_STATIC_ASSERT(false, COMPLEX_POWER_OF_A_MATRIX_IS_UNDER_CONSTRUCTION) } + + /** + * \brief Compute the matrix power. + * + * \param[out] result \f$ A^p b \f$, as specified in the constructor. + */ + template void compute(ResultType& result); + + private: + /** \brief Compute Schur decomposition of #m_A. */ + void computeSchurDecomposition(); + + /** \brief Compute atanh (inverse hyperbolic tangent). */ + Scalar atanh(const Scalar& x); + + /** \brief Compute power of 2x2 triangular matrix. */ + void compute2x2(const Scalar& p); + + /** + * \brief Compute power of triangular matrices with size > 2. + * \details This uses a Schur-Padé algorithm. + */ + void computeBig(); + + /** \brief Get suitable degree for Pade approximation. (specialized for \c RealScalar = \c double) */ + inline int getPadeDegree(double); +/* TODO + * inline int getPadeDegree(float); + * + * inline int getPadeDegree(long double); + */ + /** \brief Compute Padé approximation to matrix fractional power. */ + void computePade(int degree, const MatrixType& IminusT); + + /** \brief Get a certain coefficient of the Padé approximation. */ + inline Scalar coeff(int degree); + + const MatrixType& m_A; ///< \brief Reference to the matrix base. + const Scalar& m_p; ///< \brief Reference to 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. + MatrixType m_T; ///< \brief Triangular part of Schur decomposition. + MatrixType m_U; ///< \brief Unitary part of Schur decomposition. + MatrixType m_fT; ///< \brief #m_T to the power of #m_pfrac. + ArrayType m_logTdiag; ///< \brief Logarithm of the main diagonal of #m_T. +}; + +/******* Specialized for real exponents *******/ + +template +template +void MatrixPower::compute(ResultType& result) +{ + using std::floor; + using std::pow; + + m_pint = floor(m_p); + m_pfrac = m_p - m_pint; + + if (m_pfrac == RealScalar(0)) + computeIntPower(result); + 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 + computeBig(); + computeTmp(Scalar()); + result *= m_tmp; + } +} + +template +template +void MatrixPower::computeIntPower(ResultType& result) +{ + if (m_dimb > m_dimA) { + MatrixType tmp = MatrixType::Identity(m_A.rows(), m_A.cols()); + computeChainProduct(tmp); + result = tmp * m_b; + } else { + result = m_b; + computeChainProduct(result); + } +} + +template +template +void MatrixPower::computeChainProduct(ResultType& result) +{ + using std::frexp; + using std::ldexp; + + const bool pIsNegative = m_pint < RealScalar(0); + RealScalar p = pIsNegative? -m_pint: m_pint; + int cost = computeCost(p); + + if (pIsNegative) { + if (p * m_dimb <= cost * m_dimA) { + partialPivLuSolve(p, result); + return; + } else { + m_tmp = m_A.inverse(); + } + } else { + m_tmp = m_A; + } + while (p * m_dimb > cost * m_dimA) { + if (fmod(p, RealScalar(2)) >= RealScalar(1)) { + result = m_tmp * result; + cost--; + } + m_tmp *= m_tmp; + cost--; + p = ldexp(p, -1); + } + for (; p >= RealScalar(1); p--) + result = m_tmp * result; +} + +template +int MatrixPower::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); + cost++; + } + return cost; +} + +template +template +void MatrixPower::partialPivLuSolve(RealScalar p, ResultType& result) +{ + const PartialPivLU Asolver(m_A); + for (; p >= RealScalar(1); p--) + result = Asolver.solve(result); +} + +template +void MatrixPower::computeSchurDecomposition() +{ + const ComplexSchur schurOfA(m_A); + m_T = schurOfA.matrixT(); + m_U = schurOfA.matrixU(); +} + +template +void MatrixPower::getFractionalExponent() +{ + using std::pow; + + typedef Array RealArray; + const ComplexArray Tdiag = m_T.diagonal(); + RealScalar maxAbsEival, minAbsEival, *begin, *end; + RealArray absTdiag; + + 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++; + } +} + +template +std::complex MatrixPower::atanh(const ComplexScalar& x) +{ + using std::abs; + using std::log; + using std::sqrt; + + if (abs(x) > sqrt(NumTraits::epsilon())) + return RealScalar(0.5) * log((RealScalar(1) + x) / (RealScalar(1) - x)); + else + return x + x*x*x / RealScalar(3); +} + +template +void MatrixPower::compute2x2(const RealScalar& p) +{ + using std::abs; + using std::ceil; + using std::exp; + using std::imag; + using std::ldexp; + using std::log; + using std::pow; + using std::sinh; + + int i, j, unwindingNumber; + 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); + + if (m_T(i,i) == m_T(j,j)) + m_fT(i,j) = p * pow(m_T(i,j), p - RealScalar(1)); + else if (abs(m_T(i,i)) < ldexp(abs(m_T(j,j)), -1) || abs(m_T(j,j)) < ldexp(abs(m_T(i,i)), -1)) + m_fT(i,j) = m_T(i,j) * (m_fT(j,j) - m_fT(i,i)) / (m_T(j,j) - m_T(i,i)); + else { + // computation in previous branch is inaccurate if abs(m_T(j,j)) \approx abs(m_T(i,i)) + unwindingNumber = static_cast(ceil((imag(m_logTdiag[j] - m_logTdiag[i]) - M_PI) / (2 * M_PI))); + w = atanh((m_T(j,j) - m_T(i,i)) / (m_T(j,j) + m_T(i,i))) + ComplexScalar(0, M_PI * unwindingNumber); + m_fT(i,j) = m_T(i,j) * RealScalar(2) * exp(RealScalar(0.5) * p * (m_logTdiag[j] + m_logTdiag[i])) * + sinh(p * w) / (m_T(j,j) - m_T(i,i)); + } + } +} + +template +void MatrixPower::computeBig() +{ + using std::ldexp; + + const RealScalar maxNormForPade = 2.787629930862099e-1; + int degree, degree2, numberOfSquareRoots = 0, numberOfExtraSquareRoots = 0; + ComplexMatrix IminusT, sqrtT, T = m_T; + RealScalar normIminusT; + + while (true) { + IminusT = ComplexMatrix::Identity(m_A.rows(), m_A.cols()) - T; + normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff(); + if (normIminusT < maxNormForPade) { + degree = getPadeDegree(normIminusT); + degree2 = getPadeDegree(normIminusT * RealScalar(0.5)); + if (degree - degree2 <= 1 || numberOfExtraSquareRoots) + break; + numberOfExtraSquareRoots++; + } + MatrixSquareRootTriangular(T).compute(sqrtT); + T = sqrtT; + numberOfSquareRoots++; + } + computePade(degree, IminusT); + + for (; numberOfSquareRoots; numberOfSquareRoots--) { + compute2x2(ldexp(m_pfrac, -numberOfSquareRoots)); + m_fT *= m_fT; + } + compute2x2(m_pfrac); +} + +template +inline int MatrixPower::getPadeDegree(double normIminusT) +{ + const double maxNormForPade[] = { 1.882832775783885e-2 /* degree = 3 */ , 6.036100693089764e-2, + 1.239372725584911e-1, 1.998030690604271e-1, 2.787629930862099e-1 }; + for (int degree = 3; degree <= 7; degree++) + if (normIminusT <= maxNormForPade[degree - 3]) + return degree; + assert(false); // this line should never be reached +} + +template +void MatrixPower::computePade(int degree, const ComplexMatrix& IminusT) +{ + degree <<= 1; + m_fT = coeff(degree) * IminusT; + + for (int i = degree - 1; i; i--) { + m_fT = (ComplexMatrix::Identity(m_A.rows(), m_A.cols()) + m_fT).template triangularView() + .solve(coeff(i) * IminusT).eval(); + } + m_fT += ComplexMatrix::Identity(m_A.rows(), m_A.cols()); +} + +template +inline RealScalar MatrixPower::coeff(int i) +{ + if (i == 1) + return -m_pfrac; + else if (i & 1) + return (-m_pfrac - RealScalar(i)) / RealScalar((i<<2) + 2); + else + return (m_pfrac - RealScalar(i)) / RealScalar((i<<2) - 2); +} + +template +void MatrixPower::computeTmp(RealScalar) +{ m_tmp = (m_U * m_fT * m_U.adjoint()).real(); } + +template +void MatrixPower::computeTmp(ComplexScalar) +{ m_tmp = (m_U * m_fT * m_U.adjoint()).eval(); } + +/******* Specialized for integral exponents *******/ + +template +template +void MatrixPower::compute(ResultType& result) +{ + if (m_dimb > m_dimA) { + MatrixType tmp = MatrixType::Identity(m_dimA, m_dimA); + computeChainProduct(tmp); + result = tmp * m_b; + } else { + result = m_b; + computeChainProduct(result); + } +} + +template +int MatrixPower::computeCost(const IntExponent& p) +{ + int cost = 0; + IntExponent tmp = p; + for (tmp = p >> 1; tmp; tmp >>= 1) + cost++; + for (tmp = IntExponent(1); tmp <= p; tmp <<= 1) + if (tmp & p) cost++; + return cost; +} + +template +template +void MatrixPower::partialPivLuSolve(IntExponent p, ResultType& result) +{ + const PartialPivLU Asolver(m_A); + for(; p; p--) + result = Asolver.solve(result); +} + +template +template +void MatrixPower::computeChainProduct(ResultType& result) +{ + const bool pIsNegative = m_p < IntExponent(0); + IntExponent p = pIsNegative? -m_p: m_p; + int cost = computeCost(p); + + if (pIsNegative) { + if (p * m_dimb <= cost * m_dimA) { + partialPivLuSolve(p, result); + return; + } else { m_tmp = m_A.inverse(); } + } 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; +} + +/******* Specialized for complex exponents *******/ + +template +template +void MatrixPower,PlainObject,IsInteger>::compute(ResultType& result) +{ + using std::floor; + using std::pow; + + if (m_dimA == 1) + result = pow(m_A(0,0), m_p) * m_b; + else { + computeSchurDecomposition(); + if (m_dimA == 2) + compute2x2(m_p); + else + computeBig(); + result = m_U * m_fT * m_U.adjoint(); + } +} + +template +void MatrixPower,PlainObject,IsInteger>::computeSchurDecomposition() +{ + const ComplexSchur schurOfA(m_A); + m_T = schurOfA.matrixT(); + m_U = schurOfA.matrixU(); + m_logTdiag = m_T.diagonal().array().log(); +} + +template +typename MatrixType::Scalar MatrixPower,PlainObject,IsInteger>::atanh(const Scalar& x) +{ + using std::abs; + using std::log; + using std::sqrt; + + if (abs(x) > sqrt(NumTraits::epsilon())) + return RealScalar(0.5) * log((RealScalar(1) + x) / (RealScalar(1) - x)); + else + return x + x*x*x / RealScalar(3); +} + +template +void MatrixPower,PlainObject,IsInteger>::compute2x2(const Scalar& p) +{ + using std::abs; + using std::ceil; + using std::exp; + using std::imag; + using std::ldexp; + using std::log; + using std::pow; + using std::sinh; + + int i, j, unwindingNumber; + Scalar 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); + + if (m_T(i,i) == m_T(j,j)) + m_fT(i,j) = p * pow(m_T(i,j), p - RealScalar(1)); + else if (abs(m_T(i,i)) < ldexp(abs(m_T(j,j)), -1) || abs(m_T(j,j)) < ldexp(abs(m_T(i,i)), -1)) + m_fT(i,j) = m_T(i,j) * (m_fT(j,j) - m_fT(i,i)) / (m_T(j,j) - m_T(i,i)); + else { + // computation in previous branch is inaccurate if abs(m_T(j,j)) \approx abs(m_T(i,i)) + unwindingNumber = static_cast(ceil((imag(m_logTdiag[j] - m_logTdiag[i]) - M_PI) / (2 * M_PI))); + w = atanh((m_T(j,j) - m_T(i,i)) / (m_T(j,j) + m_T(i,i))) + Scalar(0, M_PI * unwindingNumber); + m_fT(i,j) = m_T(i,j) * RealScalar(2) * exp(RealScalar(0.5) * p * (m_logTdiag[j] + m_logTdiag[i])) * + sinh(p * w) / (m_T(j,j) - m_T(i,i)); + } + } +} + +template +void MatrixPower,PlainObject,IsInteger>::computeBig() +{ + using std::abs; + using std::ceil; + using std::frexp; + using std::ldexp; + + const RealScalar maxNormForPade = 2.787629930862099e-1; + int degree, degree2, numberOfSquareRoots, numberOfExtraSquareRoots = 0; + MatrixType IminusT, sqrtT, T = m_T; + RealScalar normIminusT; + Scalar p; +/* + frexp(abs(m_p), &numberOfSquareRoots); + if (numberOfSquareRoots > 0) + p = m_p * ldexp(RealScalar(1), -numberOfSquareRoots); + else { + p = m_p; + numberOfSquareRoots = 0; + } +*/ + while (true) { + IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T; + normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff(); + if (normIminusT < maxNormForPade) { + degree = getPadeDegree(normIminusT); + degree2 = getPadeDegree(normIminusT * RealScalar(0.5)); + if (degree - degree2 <= 1 || numberOfExtraSquareRoots) + break; + numberOfExtraSquareRoots++; + } + MatrixSquareRootTriangular(T).compute(sqrtT); + T = sqrtT; + numberOfSquareRoots++; + } + computePade(degree, IminusT); + + for (; numberOfSquareRoots; numberOfSquareRoots--) { + compute2x2(p * ldexp(RealScalar(1), -numberOfSquareRoots)); + m_fT *= m_fT; + } + compute2x2(p); +} + +template +inline int MatrixPower,PlainObject,IsInteger>::getPadeDegree(double normIminusT) +{ + const double maxNormForPade[] = { 1.882832775783885e-2 /* degree = 3 */ , 6.036100693089764e-2, + 1.239372725584911e-1, 1.998030690604271e-1, 2.787629930862099e-1 }; + for (int degree = 3; degree <= 7; degree++) + if (normIminusT <= maxNormForPade[degree - 3]) + return degree; + assert(false); // this line should never be reached +} + +template +void MatrixPower,PlainObject,IsInteger>::computePade(int degree, const MatrixType& IminusT) +{ + degree <<= 1; + m_fT = coeff(degree) * IminusT; + + for (int i = degree - 1; i; i--) { + m_fT = (MatrixType::Identity(m_A.rows(), m_A.cols()) + m_fT).template triangularView() + .solve(coeff(i) * IminusT).eval(); + } + m_fT += MatrixType::Identity(m_A.rows(), m_A.cols()); +} + +template +inline typename MatrixType::Scalar MatrixPower,PlainObject,IsInteger>::coeff(int i) +{ + if (i == 1) + return -m_p; + else if (i & 1) + return (-m_p - RealScalar(i)) / RealScalar((i<<2) + 2); + else + return (m_p - RealScalar(i)) / RealScalar((i<<2) - 2); +} + +/** + * \ingroup MatrixFunctions_Module + * + * \brief Proxy for the matrix power multiplied by another matrix + * (expression). + * + * \tparam MatrixType type of the base, a matrix (expression). + * \tparam ExponentType type of the exponent, a scalar. + * \tparam Derived type of the multiplier, a matrix (expression). + * + * This class holds the arguments to the matrix expression 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 most of the time this is the + * only way it is used. + */ +template class MatrixPowerMultiplied +: public ReturnByValue > +{ + public: + typedef typename Derived::Index Index; + + /** + * \brief Constructor. + * + * \param[in] A %Matrix (expression), the base of the matrix power. + * \param[in] p scalar, the exponent of the matrix power. + * \param[in] b %Matrix (expression), the multiplier. + */ + MatrixPowerMultiplied(const MatrixType& A, const ExponentType& p, const Derived& b) + : m_A(A), m_p(p), m_b(b) { } + + /** + * \brief Compute the matrix exponential. + * + * \param[out] result \f$ A^p b \f$ where \c A ,\c p and \c b are as in + * the constructor. + */ + template + inline void evalTo(ResultType& result) const + { + typedef typename Derived::PlainObject PlainObject; + const typename MatrixType::PlainObject Aevaluated = m_A.eval(); + const PlainObject bevaluated = m_b.eval(); + MatrixPower mp(Aevaluated, m_p, bevaluated); + mp.compute(result); + } + + Index rows() const { return m_b.rows(); } + Index cols() const { return m_b.cols(); } + + private: + const MatrixType& m_A; + const ExponentType& m_p; + const Derived& m_b; + + MatrixPowerMultiplied& operator=(const MatrixPowerMultiplied&); +}; + +/** + * \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. + * + * 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 + * MatrixBase::pow() and related functions and most of the + * time this is the only way it is used. + */ +template class MatrixPowerReturnValue +: public ReturnByValue > +{ + public: + typedef typename Derived::Index Index; + + /** + * \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) + : m_A(A), m_p(p) { } + + /** + * \brief Return the matrix power multiplied by %Matrix \c b. + * + * 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 (exporession), the multiplier. + */ + template + const MatrixPowerMultiplied operator*(const MatrixBase& b) const + { return MatrixPowerMultiplied(m_A, m_p, b.derived()); } + + /** + * \brief Compute the matrix power. + * + * \param[out] result \f$ A^p \f$ where \c A and \c p are as in the + * constructor. + */ + template + inline void evalTo(ResultType& result) const + { + typedef typename Derived::PlainObject PlainObject; + const PlainObject Aevaluated = m_A.eval(); + const PlainObject Identity = PlainObject::Identity(m_A.rows(), m_A.cols()); + MatrixPower mp(Aevaluated, m_p, Identity); + mp.compute(result); + } + + Index rows() const { return m_A.rows(); } + Index cols() const { return m_A.cols(); } + + private: + const Derived& m_A; + const ExponentType& m_p; + + MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&); +}; + +namespace internal { + template + struct traits > + { + typedef typename Derived::PlainObject ReturnType; + }; + + template + struct traits > + { + typedef typename Derived::PlainObject ReturnType; + }; +} + +template +template +const MatrixPowerReturnValue MatrixBase::pow(const ExponentType& p) const +{ + eigen_assert(rows() == cols()); + return MatrixPowerReturnValue(derived(), p); +} + +} // end namespace Eigen + +#endif // EIGEN_MATRIX_POWER