Remove mat.pow * vec specialization, which causes segfault for mat.pow * mat.pow

4 files changed, 367 insertions, 715 deletions

diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions
index 124175841..0991817d5 100644
--- a/unsupported/Eigen/MatrixFunctions
+++ b/unsupported/Eigen/MatrixFunctions
@@ -59,7 +59,6 @@
 #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 e75fc25b4..12628df1c 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
+// Copyright (C) 2012, 2013 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
@@ -12,6 +12,248 @@
 
 namespace Eigen {
 
+namespace MatrixPowerHelper {
+
+template<typename MatrixPowerType>
+class ReturnValue : public ReturnByValue< ReturnValue<MatrixPowerType> >
+{
+  public:
+    typedef typename MatrixPowerType::PlainObject::RealScalar RealScalar;
+    typedef typename MatrixPowerType::PlainObject::Index Index;
+
+    ReturnValue(MatrixPowerType& pow, RealScalar p) : m_pow(pow), m_p(p)
+    { }
+
+    template<typename ResultType>
+    inline void evalTo(ResultType& res) const
+    { m_pow.compute(res, m_p); }
+
+    Index rows() const { return m_pow.rows(); }
+    Index cols() const { return m_pow.cols(); }
+
+  private:
+    MatrixPowerType& m_pow;
+    const RealScalar m_p;
+    ReturnValue& operator=(const ReturnValue&);
+};
+
+} // namespace MatrixPowerHelper
+
+template<typename MatrixType>
+class MatrixPowerAtomic
+{
+  private:
+    enum {
+      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
+    };
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    typedef std::complex<RealScalar> ComplexScalar;
+    typedef typename MatrixType::Index Index;
+    typedef Array< Scalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > ArrayType;
+
+    const MatrixType& m_A;
+    RealScalar m_p;
+
+    void computePade(int degree, const MatrixType& IminusT, MatrixType& res) const;
+    void compute2x2(MatrixType& res, RealScalar p) const;
+    void computeBig(MatrixType& res) const;
+    static int getPadeDegree(float normIminusT);
+    static int getPadeDegree(double normIminusT);
+    static int getPadeDegree(long double normIminusT);
+    static ComplexScalar computeSuperDiag(const ComplexScalar&, const ComplexScalar&, RealScalar p);
+    static RealScalar computeSuperDiag(RealScalar, RealScalar, RealScalar p);
+
+  public:
+    MatrixPowerAtomic(const MatrixType& T, RealScalar p);
+    void compute(MatrixType& res) const;
+};
+
+template<typename MatrixType>
+MatrixPowerAtomic<MatrixType>::MatrixPowerAtomic(const MatrixType& T, RealScalar p) :
+  m_A(T), m_p(p)
+{ eigen_assert(T.rows() == T.cols()); }
+
+template<typename MatrixType>
+void MatrixPowerAtomic<MatrixType>::compute(MatrixType& res) const
+{
+  res.resizeLike(m_A);
+  switch (m_A.rows()) {
+    case 0:
+      break;
+    case 1:
+      res(0,0) = std::pow(m_A(0,0), m_p);
+      break;
+    case 2:
+      compute2x2(res, m_p);
+      break;
+    default:
+      computeBig(res);
+  }
+}
+
+template<typename MatrixType>
+void MatrixPowerAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res) const
+{
+  int i = degree<<1;
+  res = (m_p-degree) / ((i-1)<<1) * IminusT;
+  for (--i; i; --i) {
+    res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView<Upper>()
+	.solve((i==1 ? -m_p : i&1 ? (-m_p-(i>>1))/(i<<1) : (m_p-(i>>1))/((i-1)<<1)) * IminusT).eval();
+  }
+  res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
+}
+
+// This function assumes that res has the correct size (see bug 614)
+template<typename MatrixType>
+void MatrixPowerAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
+{
+  using std::abs;
+  using std::pow;
+  
+  ArrayType logTdiag = m_A.diagonal().array().log();
+  res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
+
+  for (Index i=1; i < m_A.cols(); ++i) {
+    res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
+    if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i))
+      res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
+    else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1)))
+      res.coeffRef(i-1,i) = (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));
+    else
+      res.coeffRef(i-1,i) = computeSuperDiag(m_A.coeff(i,i), m_A.coeff(i-1,i-1), p);
+    res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
+  }
+}
+
+template<typename MatrixType>
+void MatrixPowerAtomic<MatrixType>::computeBig(MatrixType& res) 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-1L:   // double-double
+						   9.134603732914548552537150753385375e-2L; // quadruple precision
+  MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
+  RealScalar normIminusT;
+  int degree, degree2, numberOfSquareRoots = 0;
+  bool hasExtraSquareRoot = false;
+
+  /* FIXME
+   * For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite
+   * loop.  We should move 0 eigenvalues to bottom right corner.  We need not
+   * worry about tiny values (e.g. 1e-300) because they will reach 1 if
+   * repetitively sqrt'ed.
+   *
+   * If the 0 eigenvalues are semisimple, they can form a 0 matrix at the
+   * bottom right corner.
+   *
+   * [ T  A ]^p   [ T^p  (T^-1 T^p A) ]
+   * [      ]   = [                   ]
+   * [ 0  0 ]     [  0         0      ]
+   */
+  for (Index i=0; i < m_A.cols(); ++i)
+    eigen_assert(m_A(i,i) != RealScalar(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/2);
+      if (degree - degree2 <= 1 || hasExtraSquareRoot)
+	break;
+      hasExtraSquareRoot = true;
+    }
+    MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
+    T = sqrtT.template triangularView<Upper>();
+    ++numberOfSquareRoots;
+  }
+  computePade(degree, IminusT, res);
+
+  for (; numberOfSquareRoots; --numberOfSquareRoots) {
+    compute2x2(res, std::ldexp(m_p, -numberOfSquareRoots));
+    res = res.template triangularView<Upper>() * res;
+  }
+  compute2x2(res, m_p);
+}
+  
+template<typename MatrixType>
+inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(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;
+}
+
+template<typename MatrixType>
+inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(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;
+}
+
+template<typename MatrixType>
+inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(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;
+}
+
+template<typename MatrixType>
+inline typename MatrixPowerAtomic<MatrixType>::ComplexScalar
+MatrixPowerAtomic<MatrixType>::computeSuperDiag(const ComplexScalar& curr, const ComplexScalar& prev, RealScalar p)
+{
+  ComplexScalar logCurr = std::log(curr);
+  ComplexScalar logPrev = std::log(prev);
+  int unwindingNumber = std::ceil((numext::imag(logCurr - logPrev) - M_PI) / (2*M_PI));
+  ComplexScalar w = numext::atanh2(curr - prev, curr + prev) + ComplexScalar(0, M_PI*unwindingNumber);
+  return RealScalar(2) * std::exp(RealScalar(0.5) * p * (logCurr + logPrev)) * std::sinh(p * w) / (curr - prev);
+}
+
+template<typename MatrixType>
+inline typename MatrixPowerAtomic<MatrixType>::RealScalar
+MatrixPowerAtomic<MatrixType>::computeSuperDiag(RealScalar curr, RealScalar prev, RealScalar p)
+{
+  RealScalar w = numext::atanh2(curr - prev, curr + prev);
+  return 2 * std::exp(p * (std::log(curr) + std::log(prev)) / 2) * std::sinh(p * w) / (curr - prev);
+}
+
 /**
  * \ingroup MatrixFunctions_Module
  *
@@ -24,10 +266,22 @@ namespace Eigen {
  * to an arbitrary real power.
  */
 template<typename MatrixType>
-class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<MatrixType>,MatrixType>
+class MatrixPowerTriangular
 {
+  private:
+    enum {
+      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+      Options = MatrixType::Options,
+      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
+    };
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    typedef typename MatrixType::Index Index;
+
   public:
-    EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(MatrixPowerTriangular)
+    typedef MatrixType PlainObject;
 
     /**
      * \brief Constructor.
@@ -37,10 +291,9 @@ class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<Matri
      * The class stores a reference to A, so it should not be changed
      * (or destroyed) before evaluation.
      */
-    explicit MatrixPowerTriangular(const MatrixType& A) : Base(A), m_T(Base::m_A)
-    { }
+    explicit MatrixPowerTriangular(const MatrixType& A) : m_A(A), m_conditionNumber(0)
+    { eigen_assert(A.rows() == A.cols()); }
 
-  #ifdef EIGEN_PARSED_BY_DOXYGEN
     /**
      * \brief Returns the matrix power.
      *
@@ -48,8 +301,8 @@ class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<Matri
      * \return The expression \f$ A^p \f$, where A is specified in the
      * constructor.
      */
-    const MatrixPowerBaseReturnValue<MatrixPowerTriangular<MatrixType>,MatrixType> operator()(RealScalar p);
-  #endif
+    const MatrixPowerHelper::ReturnValue<MatrixPowerTriangular> operator()(RealScalar p)
+    { return MatrixPowerHelper::ReturnValue<MatrixPowerTriangular>(*this, p); }
 
     /**
      * \brief Compute the matrix power.
@@ -59,34 +312,19 @@ class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<Matri
      * constructor.
      */
     void compute(MatrixType& res, RealScalar p);
-
-    /**
-     * \brief Compute the matrix power multiplied by another matrix.
-     *
-     * \param[in]  b    a matrix with the same rows as A.
-     * \param[in]  p    exponent, a real scalar.
-     * \param[out] res  \f$ A^p b \f$, where A is specified in the
-     * constructor.
-     */
-    template<typename Derived, typename ResultType>
-    void compute(const Derived& b, ResultType& res, RealScalar p);
+    
+    Index rows() const { return m_A.rows(); }
+    Index cols() const { return m_A.cols(); }
 
   private:
-    EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(MatrixPowerTriangular)
-
-    const TriangularView<MatrixType,Upper> m_T;
-
+    typename MatrixType::Nested m_A;
+    MatrixType m_tmp;
+    RealScalar m_conditionNumber;
     RealScalar modfAndInit(RealScalar, RealScalar*);
 
-    template<typename Derived, typename ResultType>
-    void apply(const Derived&, ResultType&, bool&);
-
     template<typename ResultType>
     void computeIntPower(ResultType&, RealScalar);
 
-    template<typename Derived, typename ResultType>
-    void computeIntPower(const Derived&, ResultType&, RealScalar);
-
     template<typename ResultType>
     void computeFracPower(ResultType&, RealScalar);
 };
@@ -94,41 +332,25 @@ class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<Matri
 template<typename MatrixType>
 void MatrixPowerTriangular<MatrixType>::compute(MatrixType& res, RealScalar p)
 {
-  switch (m_A.cols()) {
+  switch (cols()) {
     case 0:
       break;
     case 1:
-      res(0,0) = std::pow(m_T.coeff(0,0), p);
+      res(0,0) = std::pow(m_A.coeff(0,0), p);
       break;
     default:
       RealScalar intpart, x = modfAndInit(p, &intpart);
-      res = MatrixType::Identity(m_A.rows(), m_A.cols());
       computeIntPower(res, intpart);
       computeFracPower(res, x);
   }
 }
 
 template<typename MatrixType>
-template<typename Derived, typename ResultType>
-void MatrixPowerTriangular<MatrixType>::compute(const Derived& b, ResultType& res, RealScalar p)
-{
-  switch (m_A.cols()) {
-    case 0:
-      break;
-    case 1:
-      res = std::pow(m_T.coeff(0,0), p) * b;
-      break;
-    default:
-      RealScalar intpart, x = modfAndInit(p, &intpart);
-      computeIntPower(b, res, intpart);
-      computeFracPower(res, x);
-  }
-}
-
-template<typename MatrixType>
 typename MatrixPowerTriangular<MatrixType>::RealScalar
 MatrixPowerTriangular<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
 {
+  typedef Array< RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > RealArray;
+
   *intpart = std::floor(x);
   RealScalar res = x - *intpart;
 
@@ -137,7 +359,7 @@ MatrixPowerTriangular<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart
     m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff();
   }
 
-  if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber,res)) {
+  if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber, res)) {
     --res;
     ++*intpart;
   }
@@ -145,87 +367,31 @@ MatrixPowerTriangular<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart
 }
 
 template<typename MatrixType>
-template<typename Derived, typename ResultType>
-void MatrixPowerTriangular<MatrixType>::apply(const Derived& b, ResultType& res, bool& init)
-{
-  if (init)
-    res = m_tmp1.template triangularView<Upper>() * res;
-  else {
-    init = true;
-    res.noalias() = m_tmp1.template triangularView<Upper>() * b;
-  }
-}
-
-template<typename MatrixType>
 template<typename ResultType>
 void MatrixPowerTriangular<MatrixType>::computeIntPower(ResultType& res, RealScalar p)
 {
   RealScalar pp = std::abs(p);
 
-  if (p<0)  m_tmp1 = m_T.solve(MatrixType::Identity(m_A.rows(), m_A.cols()));
-  else      m_tmp1 = m_T;
+  if (p<0)  m_tmp = m_A.template triangularView<Upper>().solve(MatrixType::Identity(rows(), cols()));
+  else      m_tmp = m_A.template triangularView<Upper>();
 
+  res = MatrixType::Identity(rows(), cols());
   while (pp >= 1) {
     if (std::fmod(pp, 2) >= 1)
-      res = m_tmp1.template triangularView<Upper>() * res;
-    m_tmp1 = m_tmp1.template triangularView<Upper>() * m_tmp1;
+      res.template triangularView<Upper>() = m_tmp.template triangularView<Upper>() * res;
+    m_tmp.template triangularView<Upper>() = m_tmp.template triangularView<Upper>() * m_tmp;
     pp /= 2;
   }
 }
 
 template<typename MatrixType>
-template<typename Derived, typename ResultType>
-void MatrixPowerTriangular<MatrixType>::computeIntPower(const Derived& b, ResultType& res, RealScalar p)
-{
-  if (b.cols() >= m_A.cols()) {
-    m_tmp2 = MatrixType::Identity(m_A.rows(), m_A.cols());
-    computeIntPower(m_tmp2, p);
-    res.noalias() = m_tmp2.template triangularView<Upper>() * b;
-  }
-  else {
-    RealScalar pp = std::abs(p);
-    int squarings, applyings = internal::binary_powering_cost(pp, &squarings);
-    bool init = false;
-
-    if (p==0) {
-      res = b;
-      return;
-    }
-    else if (p>0) {
-      m_tmp1 = m_T;
-    }
-    else if (b.cols()*(pp-applyings) <= m_A.cols()*squarings) {
-      res = m_T.solve(b);
-      for (--pp; pp >= 1; --pp)
-	res = m_T.solve(res);
-      return;
-    }
-    else {
-      m_tmp1 = m_T.solve(MatrixType::Identity(m_A.rows(), m_A.cols()));
-    }
-
-    while (b.cols()*(pp-applyings) > m_A.cols()*squarings) {
-      if (std::fmod(pp, 2) >= 1) {
-	apply(b, res, init);
-	--applyings;
-      }
-      m_tmp1 = m_tmp1.template triangularView<Upper>() * m_tmp1;
-      --squarings;
-      pp /= 2;
-    }
-    for (; pp >= 1; --pp)
-      apply(b, res, init);
-  }
-}
-
-template<typename MatrixType>
 template<typename ResultType>
 void MatrixPowerTriangular<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
 {
   if (p) {
     eigen_assert(m_conditionNumber);
-    MatrixPowerTriangularAtomic<MatrixType>(m_A).compute(m_tmp1, p);
-    res = m_tmp1.template triangularView<Upper>() * res;
+    MatrixPowerAtomic<MatrixType>(m_A, p).compute(m_tmp);
+    res = m_tmp * res;
   }
 }
 
@@ -249,10 +415,22 @@ void MatrixPowerTriangular<MatrixType>::computeFracPower(ResultType& res, RealSc
  * Output: \verbinclude MatrixPower_optimal.out
  */
 template<typename MatrixType>
-class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
+class MatrixPower
 {
+  private:
+    enum {
+      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+      Options = MatrixType::Options,
+      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
+    };
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    typedef typename MatrixType::Index Index;
+
   public:
-    EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(MatrixPower)
+    typedef MatrixType PlainObject;
 
     /**
      * \brief Constructor.
@@ -262,10 +440,9 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
      * The class stores a reference to A, so it should not be changed
      * (or destroyed) before evaluation.
      */
-    explicit MatrixPower(const MatrixType& A) : Base(A)
-    { }
+    explicit MatrixPower(const MatrixType& A) : m_A(A), m_conditionNumber(0)
+    { eigen_assert(A.rows() == A.cols()); }
 
-  #ifdef EIGEN_PARSED_BY_DOXYGEN
     /**
      * \brief Returns the matrix power.
      *
@@ -273,8 +450,8 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
      * \return The expression \f$ A^p \f$, where A is specified in the
      * constructor.
      */
-    const MatrixPowerBaseReturnValue<MatrixPower<MatrixType>,MatrixType> operator()(RealScalar p);
-  #endif
+    const MatrixPowerHelper::ReturnValue<MatrixPower> operator()(RealScalar p)
+    { return MatrixPowerHelper::ReturnValue<MatrixPower>(*this, p); }
 
     /**
      * \brief Compute the matrix power.
@@ -284,45 +461,45 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
      * constructor.
      */
     void compute(MatrixType& res, RealScalar p);
-
-    /**
-     * \brief Compute the matrix power multiplied by another matrix.
-     *
-     * \param[in]  b    a matrix with the same rows as A.
-     * \param[in]  p    exponent, a real scalar.
-     * \param[out] res  \f$ A^p b \f$, where A is specified in the
-     * constructor.
-     */
-    template<typename Derived, typename ResultType>
-    void compute(const Derived& b, ResultType& res, RealScalar p);
+    
+    Index rows() const { return m_A.rows(); }
+    Index cols() const { return m_A.cols(); }
 
   private:
-    EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(MatrixPower)
+    typedef std::complex<RealScalar> ComplexScalar;
+    typedef Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime,
+              MaxColsAtCompileTime > ComplexMatrix;
 
-    typedef Matrix<std::complex<RealScalar>,   RowsAtCompileTime,   ColsAtCompileTime,
-				    Options,MaxRowsAtCompileTime,MaxColsAtCompileTime> ComplexMatrix;
-    static const bool m_OKforLU = RowsAtCompileTime == Dynamic || RowsAtCompileTime > 4;    
+    typename MatrixType::Nested m_A;
+    MatrixType m_tmp;
     ComplexMatrix m_T, m_U, m_fT;
+    RealScalar m_conditionNumber;
 
     RealScalar modfAndInit(RealScalar, RealScalar*);
 
-    template<typename Derived, typename ResultType>
-    void apply(const Derived&, ResultType&, bool&);
-
     template<typename ResultType>
     void computeIntPower(ResultType&, RealScalar);
 
-    template<typename Derived, typename ResultType>
-    void computeIntPower(const Derived&, ResultType&, RealScalar);
-
     template<typename ResultType>
     void computeFracPower(ResultType&, RealScalar);
+
+    template<int Rows, int Cols, int Opt, int MaxRows, int MaxCols>
+    static void revertSchur(
+        Matrix< ComplexScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res,
+        const ComplexMatrix& T,
+        const ComplexMatrix& U);
+
+    template<int Rows, int Cols, int Opt, int MaxRows, int MaxCols>
+    static void revertSchur(
+        Matrix< RealScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res,
+        const ComplexMatrix& T,
+        const ComplexMatrix& U);
 };
 
 template<typename MatrixType>
 void MatrixPower<MatrixType>::compute(MatrixType& res, RealScalar p)
 {
-  switch (m_A.cols()) {
+  switch (cols()) {
     case 0:
       break;
     case 1:
@@ -330,32 +507,17 @@ void MatrixPower<MatrixType>::compute(MatrixType& res, RealScalar p)
       break;
     default:
       RealScalar intpart, x = modfAndInit(p, &intpart);
-      res = MatrixType::Identity(m_A.rows(), m_A.cols());
       computeIntPower(res, intpart);
       computeFracPower(res, x);
   }
 }
 
 template<typename MatrixType>
-template<typename Derived, typename ResultType>
-void MatrixPower<MatrixType>::compute(const Derived& b, ResultType& res, RealScalar p)
+typename MatrixPower<MatrixType>::RealScalar
+MatrixPower<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
 {
-  switch (m_A.cols()) {
-    case 0:
-      break;
-    case 1:
-      res = std::pow(m_A.coeff(0,0), p) * b;
-      break;
-    default:
-      RealScalar intpart, x = modfAndInit(p, &intpart);
-      computeIntPower(b, res, intpart);
-      computeFracPower(res, x);
-  }
-}
+  typedef Array< RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > RealArray;
 
-template<typename MatrixType>
-typename MatrixPower<MatrixType>::RealScalar MatrixPower<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
-{
   *intpart = std::floor(x);
   RealScalar res = x - *intpart;
 
@@ -376,99 +538,50 @@ typename MatrixPower<MatrixType>::RealScalar MatrixPower<MatrixType>::modfAndIni
 }
 
 template<typename MatrixType>
-template<typename Derived, typename ResultType>
-void MatrixPower<MatrixType>::apply(const Derived& b, ResultType& res, bool& init)
-{
-  if (init)
-    res = m_tmp1 * res;
-  else {
-    init = true;
-    res.noalias() = m_tmp1 * b;
-  }
-}
-
-template<typename MatrixType>
 template<typename ResultType>
 void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p)
 {
   RealScalar pp = std::abs(p);
 
-  if (p<0)  m_tmp1 = m_A.inverse();
-  else      m_tmp1 = m_A;
+  if (p<0)  m_tmp = m_A.inverse();
+  else      m_tmp = m_A;
 
+  res = MatrixType::Identity(rows(), cols());
   while (pp >= 1) {
     if (std::fmod(pp, 2) >= 1)
-      res = m_tmp1 * res;
-    m_tmp1 *= m_tmp1;
+      res = m_tmp * res;
+    m_tmp *= m_tmp;
     pp /= 2;
   }
 }
 
 template<typename MatrixType>
-template<typename Derived, typename ResultType>
-void MatrixPower<MatrixType>::computeIntPower(const Derived& b, ResultType& res, RealScalar p)
-{
-  if (b.cols() >= m_A.cols()) {
-    m_tmp2 = MatrixType::Identity(m_A.rows(), m_A.cols());
-    computeIntPower(m_tmp2, p);
-    res.noalias() = m_tmp2 * b;
-  }
-  else {
-    RealScalar pp = std::abs(p);
-    int squarings, applyings = internal::binary_powering_cost(pp, &squarings);
-    bool init = false;
-
-    if (p==0) {
-      res = b;
-      return;
-    }
-    else if (p>0) {
-      m_tmp1 = m_A;
-    }
-    else if (m_OKforLU && b.cols()*(pp-applyings) <= m_A.cols()*squarings) {
-      PartialPivLU<MatrixType> A(m_A);
-      res = A.solve(b);
-      for (--pp; pp >= 1; --pp)
-	res = A.solve(res);
-      return;
-    }
-    else {
-      m_tmp1 = m_A.inverse();
-    }
-
-    while (b.cols()*(pp-applyings) > m_A.cols()*squarings) {
-      if (std::fmod(pp, 2) >= 1) {
-	apply(b, res, init);
-	--applyings;
-      }
-      m_tmp1 *= m_tmp1;
-      --squarings;
-      pp /= 2;
-    }
-    for (; pp >= 1; --pp)
-      apply(b, res, init);
-  }
-}
-
-template<typename MatrixType>
 template<typename ResultType>
 void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
 {
   if (p) {
     eigen_assert(m_conditionNumber);
-    MatrixPowerTriangularAtomic<ComplexMatrix>(m_T).compute(m_fT, p);
-    internal::recompose_complex_schur<NumTraits<Scalar>::IsComplex>::run(m_tmp1, m_fT, m_U);
-    res = m_tmp1 * res;
+    MatrixPowerAtomic<ComplexMatrix>(m_T, p).compute(m_fT);
+    revertSchur(m_tmp, m_fT, m_U);
+    res = m_tmp * res;
   }
 }
 
-namespace internal {
-
-template<typename Derived>
-struct traits<MatrixPowerReturnValue<Derived> >
-{ typedef typename Derived::PlainObject ReturnType; };
+template<typename MatrixType>
+template<int Rows, int Cols, int Opt, int MaxRows, int MaxCols>
+inline void MatrixPower<MatrixType>::revertSchur(
+    Matrix< ComplexScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res,
+    const ComplexMatrix& T,
+    const ComplexMatrix& U)
+{ res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint()); }
 
-} // namespace internal
+template<typename MatrixType>
+template<int Rows, int Cols, int Opt, int MaxRows, int MaxCols>
+inline void MatrixPower<MatrixType>::revertSchur(
+    Matrix< RealScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res,
+    const ComplexMatrix& T,
+    const ComplexMatrix& U)
+{ res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }
 
 /**
  * \ingroup MatrixFunctions_Module
@@ -484,7 +597,7 @@ struct traits<MatrixPowerReturnValue<Derived> >
  * time this is the only way it is used.
  */
 template<typename Derived>
-class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Derived> >
+class MatrixPowerReturnValue : public ReturnByValue< MatrixPowerReturnValue<Derived> >
 {
   public:
     typedef typename Derived::PlainObject PlainObject;
@@ -510,21 +623,6 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
     inline void evalTo(ResultType& res) const
     { MatrixPower<PlainObject>(m_A.eval()).compute(res, m_p); }
 
-    /**
-     * \brief Return the expression \f$ A^p b \f$.
-     *
-     * \p A and \p p are specified in the constructor.
-     *
-     * \param[in] b  the matrix (expression) to be applied.
-     */
-    template<typename OtherDerived>
-    const MatrixPowerProduct<MatrixPower<PlainObject>,PlainObject,OtherDerived>
-    operator*(const MatrixBase<OtherDerived>& b) const
-    {
-      MatrixPower<PlainObject> Apow(m_A.eval());
-      return MatrixPowerProduct<MatrixPower<PlainObject>,PlainObject,OtherDerived>(Apow, b.derived(), m_p);
-    }
-
     Index rows() const { return m_A.rows(); }
     Index cols() const { return m_A.cols(); }
 
@@ -534,6 +632,18 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
     MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&);
 };
 
+namespace internal {
+
+template<typename MatrixPowerType>
+struct traits< MatrixPowerHelper::ReturnValue<MatrixPowerType> >
+{ typedef typename MatrixPowerType::PlainObject ReturnType; };
+
+template<typename Derived>
+struct traits< MatrixPowerReturnValue<Derived> >
+{ typedef typename Derived::PlainObject ReturnType; };
+
+}
+
 template<typename Derived>
 const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(const RealScalar& p) const
 { return MatrixPowerReturnValue<Derived>(derived(), p); }
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
deleted file mode 100644
index 25b3f4496..000000000
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
+++ /dev/null
@@ -1,404 +0,0 @@
-// 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 {
-
-#define EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(Derived) \
-  typedef MatrixPowerBase<Derived, MatrixType> Base; \
-  using Base::RowsAtCompileTime; \
-  using Base::ColsAtCompileTime; \
-  using Base::Options; \
-  using Base::MaxRowsAtCompileTime; \
-  using Base::MaxColsAtCompileTime; \
-  typedef typename Base::Scalar     Scalar; \
-  typedef typename Base::RealScalar RealScalar; \
-  typedef typename Base::RealArray  RealArray;
-
-#define EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(Derived) \
-  using Base::m_A; \
-  using Base::m_tmp1; \
-  using Base::m_tmp2; \
-  using Base::m_conditionNumber;
-
-template<typename Derived, typename MatrixType>
-class MatrixPowerBaseReturnValue : public ReturnByValue<MatrixPowerBaseReturnValue<Derived,MatrixType> >
-{
-  public:
-    typedef typename MatrixType::RealScalar RealScalar;
-    typedef typename MatrixType::Index Index;
-
-    MatrixPowerBaseReturnValue(Derived& pow, RealScalar p) : m_pow(pow), m_p(p)
-    { }
-
-    template<typename ResultType>
-    inline void evalTo(ResultType& res) const
-    { m_pow.compute(res, m_p); }
-
-    template<typename OtherDerived>
-    const MatrixPowerProduct<Derived,MatrixType,OtherDerived> operator*(const MatrixBase<OtherDerived>& b) const
-    { return MatrixPowerProduct<Derived,MatrixType,OtherDerived>(m_pow, b.derived(), m_p); }
-
-    Index rows() const { return m_pow.rows(); }
-    Index cols() const { return m_pow.cols(); }
-
-  private:
-    Derived& m_pow;
-    const RealScalar m_p;
-    MatrixPowerBaseReturnValue& operator=(const MatrixPowerBaseReturnValue&);
-};
-
-template<typename Derived, typename MatrixType>
-class MatrixPowerBase
-{
-  private:
-    Derived& derived() { return *static_cast<Derived*>(this); }
-
-  public:
-    enum {
-      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
-      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
-      Options = MatrixType::Options,
-      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
-      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
-    };
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
-    typedef typename MatrixType::Index Index;
-
-    explicit MatrixPowerBase(const MatrixType& A) :
-      m_A(A),
-      m_conditionNumber(0)
-    { eigen_assert(A.rows() == A.cols()); }
-
-  #ifndef EIGEN_PARSED_BY_DOXYGEN
-    const MatrixPowerBaseReturnValue<Derived,MatrixType> operator()(RealScalar p)
-    { return MatrixPowerBaseReturnValue<Derived,MatrixType>(derived(), p); }
-  #endif
-
-    void compute(MatrixType& res, RealScalar p)
-    { derived().compute(res,p); }
-
-    template<typename OtherDerived, typename ResultType>
-    void compute(const OtherDerived& b, ResultType& res, RealScalar p)
-    { derived().compute(b,res,p); }
-    
-    Index rows() const { return m_A.rows(); }
-    Index cols() const { return m_A.cols(); }
-
-  protected:
-    typedef Array<RealScalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> RealArray;
-
-    typename MatrixType::Nested m_A;
-    MatrixType m_tmp1, m_tmp2;
-    RealScalar m_conditionNumber;
-};
-
-template<typename Derived, typename Lhs, typename Rhs>
-class MatrixPowerProduct : public MatrixBase<MatrixPowerProduct<Derived,Lhs,Rhs> >
-{
-  public:
-    typedef MatrixBase<MatrixPowerProduct> Base;
-    EIGEN_DENSE_PUBLIC_INTERFACE(MatrixPowerProduct)
-
-    MatrixPowerProduct(Derived& pow, const Rhs& b, RealScalar p) :
-      m_pow(pow),
-      m_b(b),
-      m_p(p)
-    { eigen_assert(pow.cols() == b.rows()); }
-
-    template<typename ResultType>
-    inline void evalTo(ResultType& res) const
-    { m_pow.compute(m_b, res, m_p); }
-
-    inline Index rows() const { return m_pow.rows(); }
-    inline Index cols() const { return m_b.cols(); }
-
-  private:
-    Derived& m_pow;
-    typename Rhs::Nested m_b;
-    const RealScalar m_p;
-};
-
-template<typename Derived>
-template<typename MatrixPower, typename Lhs, typename Rhs>
-Derived& MatrixBase<Derived>::lazyAssign(const MatrixPowerProduct<MatrixPower,Lhs,Rhs>& other)
-{
-  other.evalTo(derived());
-  return derived();
-}
-
-namespace internal {
-
-template<typename Derived, typename MatrixType>
-struct traits<MatrixPowerBaseReturnValue<Derived, MatrixType> >
-{ typedef MatrixType ReturnType; };
-
-template<typename Derived, typename _Lhs, typename _Rhs>
-struct traits<MatrixPowerProduct<Derived,_Lhs,_Rhs> >
-{
-  typedef MatrixXpr XprKind;
-  typedef typename remove_all<_Lhs>::type Lhs;
-  typedef typename remove_all<_Rhs>::type Rhs;
-  typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar;
-  typedef Dense StorageKind;
-  typedef typename promote_index_type<typename Lhs::Index, typename Rhs::Index>::type Index;
-
-  enum {
-    RowsAtCompileTime = traits<Lhs>::RowsAtCompileTime,
-    ColsAtCompileTime = traits<Rhs>::ColsAtCompileTime,
-    MaxRowsAtCompileTime = traits<Lhs>::MaxRowsAtCompileTime,
-    MaxColsAtCompileTime = traits<Rhs>::MaxColsAtCompileTime,
-    Flags = (MaxRowsAtCompileTime==1 ? RowMajorBit : 0)
-	  | EvalBeforeNestingBit | EvalBeforeAssigningBit | NestByRefBit,
-    CoeffReadCost = 0
-  };
-};
-
-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.noalias() = 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.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }
-};
-
-template<typename Scalar, int IsComplex = NumTraits<Scalar>::IsComplex>
-struct matrix_power_unwinder
-{
-  static inline Scalar run(const Scalar& eival, const Scalar& eival0, int unwindingNumber)
-  { return numext::atanh2(eival-eival0, eival+eival0) + Scalar(0, M_PI*unwindingNumber); }
-};
-
-template<typename Scalar>
-struct matrix_power_unwinder<Scalar,0>
-{
-  static inline Scalar run(Scalar eival, Scalar eival0, int)
-  { return numext::atanh2(eival-eival0, eival+eival0); }
-};
-
-template<typename T>
-inline int binary_powering_cost(T p, int* squarings)
-{
-  int applyings=0, tmp;
-
-  frexp(p, squarings);
-  --*squarings;
-
-  while (std::frexp(p, &tmp), tmp > 0) {
-    p -= std::ldexp(static_cast<T>(0.5), tmp);
-    ++applyings;
-  }
-  return applyings;
-}
-
-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>
-class MatrixPowerTriangularAtomic
-{
-  private:
-    enum {
-      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
-      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
-    };
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
-    typedef typename MatrixType::Index Index;
-    typedef Array<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType;
-
-    const MatrixType& m_A;
-
-    static void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p);
-    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>
-MatrixPowerTriangularAtomic<MatrixType>::MatrixPowerTriangularAtomic(const MatrixType& T) :
-  m_A(T)
-{ eigen_assert(T.rows() == T.cols()); }
-
-template<typename MatrixType>
-void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const
-{
-  res.resizeLike(m_A);
-  switch (m_A.rows()) {
-    case 0:
-      break;
-    case 1:
-      res(0,0) = std::pow(m_A(0,0), p);
-      break;
-    case 2:
-      compute2x2(res, p);
-      break;
-    default:
-      computeBig(res, p);
-  }
-}
-
-template<typename MatrixType>
-void MatrixPowerTriangularAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p)
-{
-  int i = degree<<1;
-  res = (p-degree) / ((i-1)<<1) * IminusT;
-  for (--i; i; --i) {
-    res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView<Upper>()
-	.solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : (p-(i>>1))/((i-1)<<1)) * IminusT).eval();
-  }
-  res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
-}
-
-// this function assumes that res has the correct size (see bug 614)
-template<typename MatrixType>
-void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
-{
-  using std::abs;
-  using std::pow;
-  
-  ArrayType logTdiag = m_A.diagonal().array().log();
-  res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
-
-  for (Index i=1; i < m_A.cols(); ++i) {
-    res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
-    if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) {
-      res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
-    }
-    else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) {
-      res.coeffRef(i-1,i) = (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));
-    }
-    else {
-      int unwindingNumber = std::ceil((numext::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI));
-      Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber);
-      res.coeffRef(i-1,i) = RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) * std::sinh(p * w)
-          / (m_A.coeff(i,i) - m_A.coeff(i-1,i-1));
-    }
-    res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
-  }
-}
-
-template<typename MatrixType>
-void MatrixPowerTriangularAtomic<MatrixType>::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-1L:   // double-double
-						   9.134603732914548552537150753385375e-2L; // quadruple precision
-  MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
-  RealScalar normIminusT;
-  int degree, degree2, numberOfSquareRoots = 0;
-  bool hasExtraSquareRoot = false;
-
-  /* FIXME
-   * For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite
-   * loop.  We should move 0 eigenvalues to bottom right corner.  We need not
-   * worry about tiny values (e.g. 1e-300) because they will reach 1 if
-   * repetitively sqrt'ed.
-   *
-   * [ T  A ]^p   [ T^p  (T^-1 T^p A) ]
-   * [      ]   = [                   ]
-   * [ 0  0 ]     [  0         0      ]
-   */
-  for (Index i=0; i < m_A.cols(); ++i)
-    eigen_assert(m_A(i,i) != RealScalar(0));
-
-  while (true) {
-    IminusT = MatrixType::Identity(m_A.rows(), m_A.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 || hasExtraSquareRoot)
-	break;
-      hasExtraSquareRoot = true;
-    }
-    MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
-    T = sqrtT.template triangularView<Upper>();
-    ++numberOfSquareRoots;
-  }
-  computePade(degree, IminusT, res, p);
-
-  for (; numberOfSquareRoots; --numberOfSquareRoots) {
-    compute2x2(res, std::ldexp(p,-numberOfSquareRoots));
-    res = res.template triangularView<Upper>() * res;
-  }
-  compute2x2(res, p);
-}
-
-} // namespace Eigen
-
-#endif // EIGEN_MATRIX_POWER
diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp
index 4bb2b4d03..b9d513b45 100644
--- a/unsupported/test/matrix_power.cpp
+++ b/unsupported/test/matrix_power.cpp
@@ -109,62 +109,9 @@ void testExponentLaws(const MatrixType& m, double tol)
   }
 }
 
-template<typename MatrixType, typename VectorType>
-void testProduct(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 << "testProduct: error powerm = " << relerr(v2, v3) << '\n';
-    VERIFY(v2.isApprox(v3, static_cast<RealScalar>(tol)));
-  }
-}
-
-template<typename MatrixType, typename VectorType>
-void testTriangularProduct(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) {
-    generateTriangularMatrix<MatrixType>::run(m1, m.rows());
-    MatrixPowerTriangular<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 << "testTriangularProduct: error powerm = " << relerr(v2, v3) << '\n';
-    VERIFY(v2.isApprox(v3, static_cast<RealScalar>(tol)));
-  }
-}
-
-template<typename MatrixType, typename VectorType>
-void testMatrixVector(const MatrixType& m, const VectorType& v, double tol)
-{
-  testExponentLaws(m,tol);
-  testProduct(m,v,tol);
-  testTriangularProduct(m,v,tol);
-}
-
 typedef Matrix<double,3,3,RowMajor>         Matrix3dRowMajor;
 typedef Matrix<long double,Dynamic,Dynamic> MatrixXe;
-typedef Matrix<long double,Dynamic,1>       VectorXe;
-  
+ 
 void test_matrix_power()
 {
   CALL_SUBTEST_2(test2dRotation<double>(1e-13));
@@ -174,13 +121,13 @@ void test_matrix_power()
   CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
   CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14));
 
-  CALL_SUBTEST_2(testMatrixVector(Matrix2d(),         Vector2d(),    1e-13));
-  CALL_SUBTEST_7(testMatrixVector(Matrix3dRowMajor(), MatrixXd(3,5), 1e-13));
-  CALL_SUBTEST_3(testMatrixVector(Matrix4cd(),        Vector4cd(),   1e-13));
-  CALL_SUBTEST_4(testMatrixVector(MatrixXd(8,8),      VectorXd(8),   2e-12));
-  CALL_SUBTEST_1(testMatrixVector(Matrix2f(),         Vector2f(),    1e-4));
-  CALL_SUBTEST_5(testMatrixVector(Matrix3cf(),        Vector3cf(),   1e-4));
-  CALL_SUBTEST_8(testMatrixVector(Matrix4f(),         Vector4f(),    1e-4));
-  CALL_SUBTEST_6(testMatrixVector(MatrixXf(2,2),      VectorXf(2),   1e-3)); // see bug 614
-  CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7),      VectorXe(7),   1e-13));
+  CALL_SUBTEST_2(testExponentLaws(Matrix2d(),         1e-13));
+  CALL_SUBTEST_7(testExponentLaws(Matrix3dRowMajor(), 1e-13));
+  CALL_SUBTEST_3(testExponentLaws(Matrix4cd(),        1e-13));
+  CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8),      2e-12));
+  CALL_SUBTEST_1(testExponentLaws(Matrix2f(),         1e-4));
+  CALL_SUBTEST_5(testExponentLaws(Matrix3cf(),        1e-4));
+  CALL_SUBTEST_8(testExponentLaws(Matrix4f(),         1e-4));
+  CALL_SUBTEST_6(testExponentLaws(MatrixXf(2,2),      1e-3)); // see bug 614
+  CALL_SUBTEST_9(testExponentLaws(MatrixXe(7,7),      1e-13));
 }