Implement matrix power-matrix product again

4 files changed, 344 insertions, 214 deletions

diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions
index ffebe0324..1a4d42de0 100644
--- a/unsupported/Eigen/MatrixFunctions
+++ b/unsupported/Eigen/MatrixFunctions
@@ -59,6 +59,7 @@
 #include "src/MatrixFunctions/MatrixFunction.h"
 #include "src/MatrixFunctions/MatrixSquareRoot.h"
 #include "src/MatrixFunctions/MatrixLogarithm.h"
+#include "src/MatrixFunctions/MatrixPowerBase.h"
 #include "src/MatrixFunctions/MatrixPower.h"
 
 
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
index 9b61502ed..ff2e31d83 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -12,209 +12,8 @@
 
 namespace Eigen {
 
-namespace internal {
-
-template<int IsComplex>
-struct recompose_complex_schur
-{
-  template<typename ResultType, typename MatrixType>
-  static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U)
-  { res = U * (T.template triangularView<Upper>() * U.adjoint()); }
-};
-
-template<>
-struct recompose_complex_schur<0>
-{
-  template<typename ResultType, typename MatrixType>
-  static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U)
-  { res = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }
-};
-
-template<typename T>
-inline int binary_powering_cost(T p)
-{
-  int cost, tmp;
-  frexp(p, &cost);
-  while (std::frexp(p, &tmp), tmp > 0) {
-    p -= std::ldexp(static_cast<T>(0.5), tmp);
-    ++cost;
-  }
-  return cost;
-}
-
-inline int matrix_power_get_pade_degree(float normIminusT)
-{
-  const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f };
-  int degree = 3;
-  for (; degree <= 4; ++degree)
-    if (normIminusT <= maxNormForPade[degree - 3])
-      break;
-  return degree;
-}
-
-inline int matrix_power_get_pade_degree(double normIminusT)
-{
-  const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1,
-      1.999045567181744e-1, 2.789358995219730e-1 };
-  int degree = 3;
-  for (; degree <= 7; ++degree)
-    if (normIminusT <= maxNormForPade[degree - 3])
-      break;
-  return degree;
-}
-
-inline int matrix_power_get_pade_degree(long double normIminusT)
-{
-#if   LDBL_MANT_DIG == 53
-  const int maxPadeDegree = 7;
-  const double maxNormForPade[] = { 1.884160592658218e-2L /* degree = 3 */ , 6.038881904059573e-2L, 1.239917516308172e-1L,
-      1.999045567181744e-1L, 2.789358995219730e-1L };
-#elif LDBL_MANT_DIG <= 64
-  const int maxPadeDegree = 8;
-  const double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L,
-      6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L };
-#elif LDBL_MANT_DIG <= 106
-  const int maxPadeDegree = 10;
-  const double maxNormForPade[] = { 1.0007161601787493236741409687186e-4L /* degree = 3 */ ,
-      1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L,
-      2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L,
-      1.1016843812851143391275867258512e-1L };
-#else
-  const int maxPadeDegree = 10;
-  const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ ,
-      6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L,
-      9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L,
-      3.908166513900489428442993794761185e-2L, 6.266780814639442865832535460550138e-2L,
-      9.134603732914548552537150753385375e-2L };
-#endif
-  int degree = 3;
-  for (; degree <= maxPadeDegree; ++degree)
-    if (normIminusT <= maxNormForPade[degree - 3])
-      break;
-  return degree;
-}
-} // namespace internal
-
-/* (non-doc)
- * \brief Class for computing triangular matrices to fractional power.
- *
- * \tparam MatrixType  type of the base.
- */
-template<typename MatrixType, int UpLo = Upper> class MatrixPowerTriangularAtomic
-{
-  private:
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
-    typedef Array<Scalar,
-		  EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime),
-		  1,ColMajor,
-		  EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::MaxRowsAtCompileTime,MatrixType::MaxColsAtCompileTime)> ArrayType;
-    const MatrixType& m_T;
-
-    void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) const;
-    void compute2x2(MatrixType& res, RealScalar p) const;
-    void computeBig(MatrixType& res, RealScalar p) const;
-
-  public:
-    explicit MatrixPowerTriangularAtomic(const MatrixType& T);
-    void compute(MatrixType& res, RealScalar p) const;
-};
-
-template<typename MatrixType, int UpLo>
-MatrixPowerTriangularAtomic<MatrixType,UpLo>::MatrixPowerTriangularAtomic(const MatrixType& T) :
-  m_T(T)
-{ eigen_assert(T.rows() == T.cols()); }
-
-template<typename MatrixType, int UpLo>
-void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute(MatrixType& res, RealScalar p) const
-{
-  switch (m_T.rows()) {
-    case 0:
-      break;
-    case 1:
-      res(0,0) = std::pow(m_T(0,0), p);
-      break;
-    case 2:
-      compute2x2(res, p);
-      break;
-    default:
-      computeBig(res, p);
-  }
-}
-
-template<typename MatrixType, int UpLo>
-void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computePade(int degree, const MatrixType& IminusT, MatrixType& res,
-  RealScalar p) const
-{
-  int i = degree<<1;
-  res = (p-(i>>1)) / ((i-1)<<1) * IminusT;
-  for (--i; i; --i) {
-    res = (MatrixType::Identity(m_T.rows(), m_T.cols()) + res).template triangularView<UpLo>()
-	.solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : (p-(i>>1))/((i-1)<<1)) * IminusT).eval();
-  }
-  res += MatrixType::Identity(m_T.rows(), m_T.cols());
-}
-
-template<typename MatrixType, int UpLo>
-void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute2x2(MatrixType& res, RealScalar p) const
-{
-  using std::abs;
-  using std::pow;
-  
-  ArrayType logTdiag = m_T.diagonal().array().log();
-  res(0,0) = pow(m_T(0,0), p);
-
-  for (int i=1; i < m_T.cols(); ++i) {
-    res(i,i) = pow(m_T(i,i), p);
-    if (m_T(i-1,i-1) == m_T(i,i)) {
-      res(i-1,i) = p * pow(m_T(i-1,i), p-1);
-    } else if (2*abs(m_T(i-1,i-1)) < abs(m_T(i,i)) || 2*abs(m_T(i,i)) < abs(m_T(i-1,i-1))) {
-      res(i-1,i) = m_T(i-1,i) * (res(i,i)-res(i-1,i-1)) / (m_T(i,i)-m_T(i-1,i-1));
-    } else {
-      // computation in previous branch is inaccurate if abs(m_T(i,i)) \approx abs(m_T(i-1,i-1))
-      int unwindingNumber = std::ceil(((logTdiag[i]-logTdiag[i-1]).imag() - M_PI) / (2*M_PI));
-      Scalar w = internal::atanh2(m_T(i,i)-m_T(i-1,i-1), m_T(i,i)+m_T(i-1,i-1)) + Scalar(0, M_PI*unwindingNumber);
-      res(i-1,i) = m_T(i-1,i) * RealScalar(2) * std::exp(RealScalar(0.5) * p * (logTdiag[i]+logTdiag[i-1])) *
-	  std::sinh(p * w) / (m_T(i,i) - m_T(i-1,i-1));
-    }
-  }
-}
-
-template<typename MatrixType, int UpLo>
-void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computeBig(MatrixType& res, RealScalar p) const
-{
-  const int digits = std::numeric_limits<RealScalar>::digits;
-  const RealScalar maxNormForPade = digits <=  24? 4.3386528e-1f:                           // sigle precision
-				    digits <=  53? 2.789358995219730e-1:                    // double precision
-				    digits <=  64? 2.4471944416607995472e-1L:               // extended precision
-				    digits <= 106? 1.1016843812851143391275867258512e-01:   // double-double
-						   9.134603732914548552537150753385375e-02; // quadruple precision
-  int degree, degree2, numberOfSquareRoots=0, numberOfExtraSquareRoots=0;
-  MatrixType IminusT, sqrtT, T=m_T;
-  RealScalar normIminusT;
-
-  while (true) {
-    IminusT = MatrixType::Identity(m_T.rows(), m_T.cols()) - T;
-    normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
-    if (normIminusT < maxNormForPade) {
-      degree = internal::matrix_power_get_pade_degree(normIminusT);
-      degree2 = internal::matrix_power_get_pade_degree(normIminusT/2);
-      if (degree - degree2 <= 1 || numberOfExtraSquareRoots)
-	break;
-      ++numberOfExtraSquareRoots;
-    }
-    MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
-    T = sqrtT;
-    ++numberOfSquareRoots;
-  }
-  computePade(degree, IminusT, res, p);
-
-  for (; numberOfSquareRoots; --numberOfSquareRoots) {
-    compute2x2(res, std::ldexp(p,-numberOfSquareRoots));
-    res *= res;
-  }
-  compute2x2(res, p);
-}
+template<typename MatrixType>
+class MatrixPowerEvaluator;
 
 /**
  * \ingroup MatrixFunctions_Module
@@ -281,8 +80,8 @@ template<typename MatrixType> class MatrixPower
      *
      * \param[in] p  exponent, a real scalar.
      */
-    const MatrixPowerReturnValue<MatrixPower<MatrixType> > operator()(RealScalar p)
-    { return MatrixPowerReturnValue<MatrixPower<MatrixType> >(*this, p); }
+    const MatrixPowerEvaluator<MatrixType> operator()(RealScalar p)
+    { return MatrixPowerEvaluator<MatrixType>(*this, p); }
 
     /**
      * \brief Compute the matrix power.
@@ -451,6 +250,30 @@ void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
   }
 }
 
+template<typename MatrixType, typename PlainObject>
+class MatrixPowerMatrixProduct : public MatrixPowerProductBase<MatrixPowerMatrixProduct<MatrixType,PlainObject> >
+{
+  public:
+    typedef MatrixPowerProductBase<MatrixPowerMatrixProduct<MatrixType,PlainObject> > Base;
+    EIGEN_DENSE_PUBLIC_INTERFACE(MatrixPowerMatrixProduct)
+
+    MatrixPowerMatrixProduct(MatrixPower<MatrixType>& pow, const PlainObject& b, RealScalar p)
+    : m_pow(pow), m_b(b), m_p(p) { }
+
+    template<typename ResultType>
+    inline void evalTo(ResultType& res) const
+    { m_pow.compute(m_b, res, m_p); }
+
+    Index rows() const { return m_b.rows(); }
+    Index cols() const { return m_b.cols(); }
+
+  private:
+    MatrixPower<MatrixType>& m_pow;
+    const PlainObject& m_b;
+    const RealScalar m_p;
+    MatrixPowerMatrixProduct& operator=(const MatrixPowerMatrixProduct&);
+};
+
 /**
  * \ingroup MatrixFunctions_Module
  *
@@ -500,41 +323,68 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
 };
 
 template<typename MatrixType>
-class MatrixPowerReturnValue<MatrixPower<MatrixType> >
-: public ReturnByValue<MatrixPowerReturnValue<MatrixPower<MatrixType> > >
+class MatrixPowerEvaluator
+: public ReturnByValue<MatrixPowerEvaluator<MatrixType> >
 {
   public:
     typedef typename MatrixType::RealScalar RealScalar;
     typedef typename MatrixType::Index Index;
 
-    MatrixPowerReturnValue(MatrixPower<MatrixType>& ref, RealScalar p)
+    MatrixPowerEvaluator(MatrixPower<MatrixType>& ref, RealScalar p)
     : m_pow(ref), m_p(p) { }
 
     template<typename ResultType>
     inline void evalTo(ResultType& res) const
     { m_pow.compute(res, m_p); }
 
+    template<typename Derived>
+    const MatrixPowerMatrixProduct<MatrixType, typename Derived::PlainObject> operator*(const MatrixBase<Derived>& b) const
+    { return MatrixPowerMatrixProduct<MatrixType, typename Derived::PlainObject>(m_pow, b.derived(), m_p); }
+
     Index rows() const { return m_pow.rows(); }
     Index cols() const { return m_pow.cols(); }
 
   private:
     MatrixPower<MatrixType>& m_pow;
     const RealScalar m_p;
-    MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&);
+    MatrixPowerEvaluator& operator=(const MatrixPowerEvaluator&);
 };
 
 namespace internal {
+template<typename MatrixType, typename PlainObject>
+struct nested<MatrixPowerMatrixProduct<MatrixType,PlainObject> >
+{ typedef PlainObject const& type; };
+
 template<typename Derived>
 struct traits<MatrixPowerReturnValue<Derived> >
 { typedef typename Derived::PlainObject ReturnType; };
 
 template<typename MatrixType>
-struct traits<MatrixPowerReturnValue<MatrixPower<MatrixType> > >
+struct traits<MatrixPowerEvaluator<MatrixType> >
 { typedef MatrixType ReturnType; };
 
-template<typename Derived>
-struct traits<MatrixPowerProductBase<Derived> >
-{ typedef typename traits<Derived>::ReturnType ReturnType; };
+template<typename MatrixType, typename PlainObject>
+struct traits<MatrixPowerMatrixProduct<MatrixType,PlainObject> >
+{
+  typedef MatrixXpr XprKind;
+  typedef typename scalar_product_traits<typename MatrixType::Scalar, typename PlainObject::Scalar>::ReturnType Scalar;
+  typedef typename promote_storage_type<typename traits<MatrixType>::StorageKind,
+					typename traits<PlainObject>::StorageKind>::ret StorageKind;
+  typedef typename promote_index_type<typename traits<MatrixType>::Index,
+				      typename traits<PlainObject>::Index>::type Index;
+
+  enum {
+    RowsAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(traits<MatrixType>::RowsAtCompileTime,
+						    traits<PlainObject>::RowsAtCompileTime),
+    ColsAtCompileTime = traits<PlainObject>::ColsAtCompileTime,
+    MaxRowsAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(traits<MatrixType>::MaxRowsAtCompileTime,
+						       traits<PlainObject>::MaxRowsAtCompileTime),
+    MaxColsAtCompileTime = traits<PlainObject>::MaxColsAtCompileTime,
+    Flags = (MaxRowsAtCompileTime==1 ? RowMajorBit : 0)
+	  | EvalBeforeNestingBit | EvalBeforeAssigningBit | NestByRefBit,
+    CoeffReadCost = 0
+  };
+};
 }
 
 template<typename Derived>
@@ -544,6 +394,6 @@ const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) con
   return MatrixPowerReturnValue<Derived>(derived(), p);
 }
 
-} // end namespace Eigen
+} // namespace Eigen
 
 #endif // EIGEN_MATRIX_POWER
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
new file mode 100644
index 000000000..e5a1fec31
--- /dev/null
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
@@ -0,0 +1,247 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
+//
+// 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_BASE
+#define EIGEN_MATRIX_POWER_BASE
+
+namespace Eigen {
+
+namespace internal {
+template<int IsComplex>
+struct recompose_complex_schur
+{
+  template<typename ResultType, typename MatrixType>
+  static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U)
+  { res = U * (T.template triangularView<Upper>() * U.adjoint()); }
+};
+
+template<>
+struct recompose_complex_schur<0>
+{
+  template<typename ResultType, typename MatrixType>
+  static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U)
+  { res = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }
+};
+
+template<typename Derived>
+struct traits<MatrixPowerProductBase<Derived> > : traits<Derived>
+{ };
+
+template<typename T>
+inline int binary_powering_cost(T p)
+{
+  int cost, tmp;
+  frexp(p, &cost);
+  while (std::frexp(p, &tmp), tmp > 0) {
+    p -= std::ldexp(static_cast<T>(0.5), tmp);
+    ++cost;
+  }
+  return cost;
+}
+
+inline int matrix_power_get_pade_degree(float normIminusT)
+{
+  const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f };
+  int degree = 3;
+  for (; degree <= 4; ++degree)
+    if (normIminusT <= maxNormForPade[degree - 3])
+      break;
+  return degree;
+}
+
+inline int matrix_power_get_pade_degree(double normIminusT)
+{
+  const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1,
+      1.999045567181744e-1, 2.789358995219730e-1 };
+  int degree = 3;
+  for (; degree <= 7; ++degree)
+    if (normIminusT <= maxNormForPade[degree - 3])
+      break;
+  return degree;
+}
+
+inline int matrix_power_get_pade_degree(long double normIminusT)
+{
+#if   LDBL_MANT_DIG == 53
+  const int maxPadeDegree = 7;
+  const double maxNormForPade[] = { 1.884160592658218e-2L /* degree = 3 */ , 6.038881904059573e-2L, 1.239917516308172e-1L,
+      1.999045567181744e-1L, 2.789358995219730e-1L };
+#elif LDBL_MANT_DIG <= 64
+  const int maxPadeDegree = 8;
+  const double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L,
+      6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L };
+#elif LDBL_MANT_DIG <= 106
+  const int maxPadeDegree = 10;
+  const double maxNormForPade[] = { 1.0007161601787493236741409687186e-4L /* degree = 3 */ ,
+      1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L,
+      2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L,
+      1.1016843812851143391275867258512e-1L };
+#else
+  const int maxPadeDegree = 10;
+  const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ ,
+      6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L,
+      9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L,
+      3.908166513900489428442993794761185e-2L, 6.266780814639442865832535460550138e-2L,
+      9.134603732914548552537150753385375e-2L };
+#endif
+  int degree = 3;
+  for (; degree <= maxPadeDegree; ++degree)
+    if (normIminusT <= maxNormForPade[degree - 3])
+      break;
+  return degree;
+}
+} // namespace internal
+
+template<typename MatrixType, int UpLo = Upper> class MatrixPowerTriangularAtomic
+{
+  private:
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    typedef Array<Scalar,
+		  EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime),
+		  1,ColMajor,
+		  EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::MaxRowsAtCompileTime,MatrixType::MaxColsAtCompileTime)> ArrayType;
+    const MatrixType& m_T;
+
+    void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) const;
+    void compute2x2(MatrixType& res, RealScalar p) const;
+    void computeBig(MatrixType& res, RealScalar p) const;
+
+  public:
+    explicit MatrixPowerTriangularAtomic(const MatrixType& T);
+    void compute(MatrixType& res, RealScalar p) const;
+};
+
+template<typename MatrixType, int UpLo>
+MatrixPowerTriangularAtomic<MatrixType,UpLo>::MatrixPowerTriangularAtomic(const MatrixType& T) :
+  m_T(T)
+{ eigen_assert(T.rows() == T.cols()); }
+
+template<typename MatrixType, int UpLo>
+void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute(MatrixType& res, RealScalar p) const
+{
+  switch (m_T.rows()) {
+    case 0:
+      break;
+    case 1:
+      res(0,0) = std::pow(m_T(0,0), p);
+      break;
+    case 2:
+      compute2x2(res, p);
+      break;
+    default:
+      computeBig(res, p);
+  }
+}
+
+template<typename MatrixType, int UpLo>
+void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computePade(int degree, const MatrixType& IminusT, MatrixType& res,
+  RealScalar p) const
+{
+  int i = degree<<1;
+  res = (p-(i>>1)) / ((i-1)<<1) * IminusT;
+  for (--i; i; --i) {
+    res = (MatrixType::Identity(m_T.rows(), m_T.cols()) + res).template triangularView<UpLo>()
+	.solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : (p-(i>>1))/((i-1)<<1)) * IminusT).eval();
+  }
+  res += MatrixType::Identity(m_T.rows(), m_T.cols());
+}
+
+template<typename MatrixType, int UpLo>
+void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute2x2(MatrixType& res, RealScalar p) const
+{
+  using std::abs;
+  using std::pow;
+  
+  ArrayType logTdiag = m_T.diagonal().array().log();
+  res(0,0) = pow(m_T(0,0), p);
+
+  for (int i=1; i < m_T.cols(); ++i) {
+    res(i,i) = pow(m_T(i,i), p);
+    if (m_T(i-1,i-1) == m_T(i,i)) {
+      res(i-1,i) = p * pow(m_T(i-1,i), p-1);
+    } else if (2*abs(m_T(i-1,i-1)) < abs(m_T(i,i)) || 2*abs(m_T(i,i)) < abs(m_T(i-1,i-1))) {
+      res(i-1,i) = m_T(i-1,i) * (res(i,i)-res(i-1,i-1)) / (m_T(i,i)-m_T(i-1,i-1));
+    } else {
+      // computation in previous branch is inaccurate if abs(m_T(i,i)) \approx abs(m_T(i-1,i-1))
+      int unwindingNumber = std::ceil(((logTdiag[i]-logTdiag[i-1]).imag() - M_PI) / (2*M_PI));
+      Scalar w = internal::atanh2(m_T(i,i)-m_T(i-1,i-1), m_T(i,i)+m_T(i-1,i-1)) + Scalar(0, M_PI*unwindingNumber);
+      res(i-1,i) = m_T(i-1,i) * RealScalar(2) * std::exp(RealScalar(0.5) * p * (logTdiag[i]+logTdiag[i-1])) *
+	  std::sinh(p * w) / (m_T(i,i) - m_T(i-1,i-1));
+    }
+  }
+}
+
+template<typename MatrixType, int UpLo>
+void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computeBig(MatrixType& res, RealScalar p) const
+{
+  const int digits = std::numeric_limits<RealScalar>::digits;
+  const RealScalar maxNormForPade = digits <=  24? 4.3386528e-1f:                           // sigle precision
+				    digits <=  53? 2.789358995219730e-1:                    // double precision
+				    digits <=  64? 2.4471944416607995472e-1L:               // extended precision
+				    digits <= 106? 1.1016843812851143391275867258512e-01:   // double-double
+						   9.134603732914548552537150753385375e-02; // quadruple precision
+  int degree, degree2, numberOfSquareRoots=0, numberOfExtraSquareRoots=0;
+  MatrixType IminusT, sqrtT, T=m_T;
+  RealScalar normIminusT;
+
+  while (true) {
+    IminusT = MatrixType::Identity(m_T.rows(), m_T.cols()) - T;
+    normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
+    if (normIminusT < maxNormForPade) {
+      degree = internal::matrix_power_get_pade_degree(normIminusT);
+      degree2 = internal::matrix_power_get_pade_degree(normIminusT/2);
+      if (degree - degree2 <= 1 || numberOfExtraSquareRoots)
+	break;
+      ++numberOfExtraSquareRoots;
+    }
+    MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
+    T = sqrtT;
+    ++numberOfSquareRoots;
+  }
+  computePade(degree, IminusT, res, p);
+
+  for (; numberOfSquareRoots; --numberOfSquareRoots) {
+    compute2x2(res, std::ldexp(p,-numberOfSquareRoots));
+    res *= res;
+  }
+  compute2x2(res, p);
+}
+
+template<typename Derived>
+class MatrixPowerProductBase : public MatrixBase<Derived>
+{
+  public:
+    typedef MatrixBase<Derived> Base;
+    typedef typename Base::PlainObject PlainObject;
+    EIGEN_DENSE_PUBLIC_INTERFACE(MatrixPowerProductBase)
+    
+    inline Index rows() const { return derived().rows(); }
+    inline Index cols() const { return derived().cols(); }
+
+    template<typename ResultType>
+    inline void evalTo(ResultType& res) const { derived().evalTo(res); }
+
+    const PlainObject& eval() const
+    {
+      m_result.resize(rows(), cols());
+      derived().evalTo(m_result);
+      return m_result;
+    }
+
+    operator const PlainObject&() const
+    { return eval(); }
+
+  protected:
+    mutable PlainObject m_result;
+};
+
+} // namespace Eigen
+
+#endif // EIGEN_MATRIX_POWER
diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp
index c5967d2eb..7d90066c8 100644
--- a/unsupported/test/matrix_power.cpp
+++ b/unsupported/test/matrix_power.cpp
@@ -86,6 +86,28 @@ void testExponentLaws(const MatrixType& m, double tol)
   }
 }
 
+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 p;
+
+  for (int i=0; i<g_repeat; ++i) {
+    generateTestMatrix<MatrixType>::run(m1, m.rows());
+    MatrixPower<MatrixType> mpow(m1);
+
+    v1 = VectorType::Random(v.rows(), v.cols());
+    p = internal::random<RealScalar>();
+
+    v2.noalias() = mpow(p) * v1;
+    v3.noalias() = mpow(p).eval() * v1;
+    std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v3) << '\n';
+    VERIFY(v2.isApprox(v3, static_cast<RealScalar>(tol)));
+  }
+}
+
 void test_matrix_power()
 {
   typedef Matrix<long double,Dynamic,Dynamic> MatrixXe;
@@ -105,4 +127,14 @@ void test_matrix_power()
   CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4));
   CALL_SUBTEST_8(testExponentLaws(Matrix4f(), 1e-4));
   CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4));
+
+  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_9(testMatrixVectorProduct(MatrixXe(7,7), MatrixXe(7,9), 1e-13));
 }