Implement complex MatrixPowerTriangular. There are still problems with real one.

2 files changed, 408 insertions, 206 deletions

diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
index 3ac11979e..a618a536f 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -12,7 +12,222 @@
 
 namespace Eigen {
 
-template<typename MatrixType> class MatrixPowerEvaluator;
+/**
+ * \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.
+ *
+ * This class is capable of computing complex upper triangular matrices raised
+ * to an arbitrary real power.
+ */
+template<typename MatrixType>
+class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<MatrixType>,MatrixType>
+{
+  public:
+    EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(MatrixPowerTriangular)
+
+    /**
+     * \brief Constructor.
+     *
+     * \param[in] A  the base of the matrix power.
+     *
+     * The class stores a reference to A, so it should not be changed
+     * (or destroyed) before evaluation.
+     */
+    explicit MatrixPowerTriangular(const MatrixType& A) : Base(A,0), m_T(Base::m_A)
+    { }
+
+  #ifdef EIGEN_PARSED_BY_DOXYGEN
+    /**
+     * \brief Returns the matrix power.
+     *
+     * \param[in] p  exponent, a real scalar.
+     * \return The expression \f$ A^p \f$, where A is specified in the
+     * constructor.
+     */
+    const MatrixPowerBaseReturnValue<MatrixPowerTriangular<MatrixType>,MatrixType> operator()(RealScalar p);
+  #endif
+
+    /**
+     * \brief Compute the matrix power.
+     *
+     * \param[in]  p    exponent, a real scalar.
+     * \param[out] res  \f$ A^p \f$ where A is specified in the
+     * 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);
+
+  private:
+    EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(MatrixPowerTriangular)
+
+    const TriangularView<MatrixType,Upper> m_T;
+
+    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<typename MatrixType>
+void MatrixPowerTriangular<MatrixType>::compute(MatrixType& res, RealScalar p)
+{
+  switch (m_A.cols()) {
+    case 0:
+      break;
+    case 1:
+      res(0,0) = std::pow(m_T.coeff(0,0), p);
+      break;
+    default:
+      RealScalar intpart, x = modfAndInit(p, &intpart);
+      res = m_Id;
+      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>::Base::RealScalar
+MatrixPowerTriangular<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
+{
+  *intpart = std::floor(x);
+  RealScalar res = x - *intpart;
+
+  if (!m_conditionNumber && res) {
+    const RealArray absTdiag = m_A.diagonal().array().abs();
+    m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff();
+  }
+
+  if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber,res)) {
+    --res;
+    ++*intpart;
+  }
+  return res;
+}
+
+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(m_Id);
+  else      m_tmp1 = m_T;
+
+  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;
+    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 = m_Id;
+    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 (m_A.cols() > 2 && 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(m_Id);
+    }
+
+    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;
+  }
+}
 
 /**
  * \ingroup MatrixFunctions_Module
@@ -44,18 +259,22 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
      *
      * \param[in] A  the base of the matrix power.
      *
-     * \warning Construct with a matrix, not a matrix expression!
+     * The class stores a reference to A, so it should not be changed
+     * (or destroyed) before evaluation.
      */
     explicit MatrixPower(const MatrixType& A) : Base(A,0)
     { }
 
+  #ifdef EIGEN_PARSED_BY_DOXYGEN
     /**
-     * \brief Return the expression \f$ A^p \f$.
+     * \brief Returns the matrix power.
      *
      * \param[in] p  exponent, a real scalar.
+     * \return The expression \f$ A^p \f$, where A is specified in the
+     * constructor.
      */
-    const MatrixPowerEvaluator<MatrixType> operator()(RealScalar p)
-    { return MatrixPowerEvaluator<MatrixType>(*this, p); }
+    const MatrixPowerBaseReturnValue<MatrixPower<MatrixType>,MatrixType> operator()(RealScalar p);
+  #endif
 
     /**
      * \brief Compute the matrix power.
@@ -242,31 +461,13 @@ void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
   }
 }
 
-template<typename Lhs, typename Rhs>
-class MatrixPowerMatrixProduct : public MatrixPowerProductBase<MatrixPowerMatrixProduct<Lhs,Rhs>,Lhs,Rhs>
-{
-  public:
-    EIGEN_MATRIX_POWER_PRODUCT_PUBLIC_INTERFACE(MatrixPowerMatrixProduct)
-
-    MatrixPowerMatrixProduct(MatrixPower<Lhs>& pow, const Rhs& 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); }
+namespace internal {
 
-    Index rows() const { return m_pow.rows(); }
-    Index cols() const { return m_b.cols(); }
+template<typename Derived>
+struct traits<MatrixPowerReturnValue<Derived> >
+{ typedef typename Derived::PlainObject ReturnType; };
 
-  private:
-    MatrixPower<Lhs>& m_pow;
-    const Rhs& m_b;
-    const RealScalar m_p;
-    MatrixPowerMatrixProduct& operator=(const MatrixPowerMatrixProduct&);
-};
+} // namespace internal
 
 /**
  * \ingroup MatrixFunctions_Module
@@ -316,10 +517,11 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
      * \param[in] b  the matrix (expression) to be applied.
      */
     template<typename OtherDerived>
-    const MatrixPowerMatrixProduct<PlainObject,OtherDerived> operator*(const MatrixBase<OtherDerived>& b) const
+    const MatrixPowerProduct<MatrixPower<PlainObject>,PlainObject,OtherDerived>
+    operator*(const MatrixBase<OtherDerived>& b) const
     {
       MatrixPower<PlainObject> Apow(m_A.eval());
-      return MatrixPowerMatrixProduct<PlainObject,OtherDerived>(Apow, b.derived(), m_p);
+      return MatrixPowerProduct<MatrixPower<PlainObject>,PlainObject,OtherDerived>(Apow, b.derived(), m_p);
     }
 
     Index rows() const { return m_A.rows(); }
@@ -331,52 +533,6 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
     MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&);
 };
 
-template<typename MatrixType>
-class MatrixPowerEvaluator : public ReturnByValue<MatrixPowerEvaluator<MatrixType> >
-{
-  public:
-    typedef typename MatrixType::RealScalar RealScalar;
-    typedef typename MatrixType::Index Index;
-
-    MatrixPowerEvaluator(MatrixPower<MatrixType>& 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 Derived>
-    const MatrixPowerMatrixProduct<MatrixType,Derived> operator*(const MatrixBase<Derived>& b) const
-    { return MatrixPowerMatrixProduct<MatrixType,Derived>(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;
-    MatrixPowerEvaluator& operator=(const MatrixPowerEvaluator&);
-};
-
-namespace internal {
-template<typename Lhs, typename Rhs>
-struct nested<MatrixPowerMatrixProduct<Lhs,Rhs> >
-{ typedef typename MatrixPowerMatrixProduct<Lhs,Rhs>::PlainObject const& type; };
-
-template<typename Derived>
-struct traits<MatrixPowerReturnValue<Derived> >
-{ typedef typename Derived::PlainObject ReturnType; };
-
-template<typename MatrixType>
-struct traits<MatrixPowerEvaluator<MatrixType> >
-{ typedef MatrixType ReturnType; };
-
-template<typename Lhs, typename Rhs>
-struct traits<MatrixPowerMatrixProduct<Lhs,Rhs> >
-: traits<MatrixPowerProductBase<MatrixPowerMatrixProduct<Lhs,Rhs>,Lhs,Rhs> >
-{ };
-} // namespace internal
-
 template<typename Derived>
 const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const
 { return MatrixPowerReturnValue<Derived>(derived(), p); }
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
index 9616659ca..b67939039 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
@@ -12,9 +12,171 @@
 
 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_Id; \
+  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, RealScalar cond) :
+      m_A(A),
+      m_Id(MatrixType::Identity(A.rows(),A.cols())),
+      m_conditionNumber(cond)
+    { 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;
+
+    const MatrixType& m_A;
+    const MatrixType m_Id;
+    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<Derived,Lhs,Rhs> > 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;
+    const Rhs& 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<int IsComplex>
-struct recompose_complex_schur
+
+template<typename Derived, typename MatrixType>
+struct traits<MatrixPowerBaseReturnValue<Derived, MatrixType> >
+{ typedef MatrixType ReturnType; };
+
+template<typename Derived, typename Lhs, typename Rhs>
+struct nested<MatrixPowerProduct<Derived,Lhs,Rhs> >
+{ typedef typename MatrixPowerProduct<Derived,Lhs,Rhs>::PlainObject const& type; };
+
+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 remove_all<MatrixPowerProduct<Derived,_Lhs,_Rhs> >::type PlainObject;
+  typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar;
+  typedef typename promote_storage_type<typename traits<Lhs>::StorageKind,
+					typename traits<Rhs>::StorageKind>::ret StorageKind;
+  typedef typename promote_index_type<typename traits<Lhs>::Index,
+				      typename traits<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<bool IsComplex> struct recompose_complex_schur;
+
+template<>
+struct recompose_complex_schur<true>
 {
   template<typename ResultType, typename MatrixType>
   static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U)
@@ -22,7 +184,7 @@ struct recompose_complex_schur
 };
 
 template<>
-struct recompose_complex_schur<0>
+struct recompose_complex_schur<false>
 {
   template<typename ResultType, typename MatrixType>
   static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U)
@@ -109,10 +271,14 @@ inline int matrix_power_get_pade_degree(long double normIminusT)
       break;
   return degree;
 }
+
 } // namespace internal
 
+template<typename MatrixType, bool IsComplex=NumTraits<typename MatrixType::RealScalar>::IsComplex>
+class MatrixPowerTriangularAtomic;
+
 template<typename MatrixType>
-class MatrixPowerTriangularAtomic
+class MatrixPowerTriangularAtomic<MatrixType,true>
 {
   private:
     enum {
@@ -136,13 +302,13 @@ class MatrixPowerTriangularAtomic
 };
 
 template<typename MatrixType>
-MatrixPowerTriangularAtomic<MatrixType>::MatrixPowerTriangularAtomic(const MatrixType& T) :
+MatrixPowerTriangularAtomic<MatrixType,true>::MatrixPowerTriangularAtomic(const MatrixType& T) :
   m_T(T),
   m_Id(MatrixType::Identity(T.rows(), T.cols()))
 { eigen_assert(T.rows() == T.cols()); }
 
 template<typename MatrixType>
-void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const
+void MatrixPowerTriangularAtomic<MatrixType,true>::compute(MatrixType& res, RealScalar p) const
 {
   switch (m_T.rows()) {
     case 0:
@@ -159,7 +325,7 @@ void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScala
 }
 
 template<typename MatrixType>
-void MatrixPowerTriangularAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res,
+void MatrixPowerTriangularAtomic<MatrixType,true>::computePade(int degree, const MatrixType& IminusT, MatrixType& res,
   RealScalar p) const
 {
   int i = degree<<1;
@@ -172,7 +338,7 @@ void MatrixPowerTriangularAtomic<MatrixType>::computePade(int degree, const Matr
 }
 
 template<typename MatrixType>
-void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
+void MatrixPowerTriangularAtomic<MatrixType,true>::compute2x2(MatrixType& res, RealScalar p) const
 {
   using std::abs;
   using std::pow;
@@ -198,7 +364,7 @@ void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealSc
 }
 
 template<typename MatrixType>
-void MatrixPowerTriangularAtomic<MatrixType>::computeBig(MatrixType& res, RealScalar p) const
+void MatrixPowerTriangularAtomic<MatrixType,true>::computeBig(MatrixType& res, RealScalar p) const
 {
   const int digits = std::numeric_limits<RealScalar>::digits;
   const RealScalar maxNormForPade = digits <=  24? 4.3386528e-1f:                           // sigle precision
@@ -212,7 +378,7 @@ void MatrixPowerTriangularAtomic<MatrixType>::computeBig(MatrixType& res, RealSc
   bool hasExtraSquareRoot=false;
 
   while (true) {
-    IminusT = MatrixType::Identity(m_T.rows(), m_T.cols()) - T;
+    IminusT = m_Id - T;
     normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
     if (normIminusT < maxNormForPade) {
       degree = internal::matrix_power_get_pade_degree(normIminusT);
@@ -234,126 +400,6 @@ void MatrixPowerTriangularAtomic<MatrixType>::computeBig(MatrixType& res, RealSc
   compute2x2(res, p);
 }
 
-#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_Id; \
-  using Base::m_tmp1; \
-  using Base::m_tmp2; \
-  using Base::m_conditionNumber;
-
-#define	EIGEN_MATRIX_POWER_PRODUCT_PUBLIC_INTERFACE(Derived) \
-  typedef MatrixPowerProductBase<Derived, Lhs, Rhs> Base; \
-  EIGEN_DENSE_PUBLIC_INTERFACE(Derived)
-
-namespace internal {
-template<typename Derived, typename _Lhs, typename _Rhs>
-struct traits<MatrixPowerProductBase<Derived,_Lhs,_Rhs> >
-{
-  typedef MatrixXpr XprKind;
-  typedef typename remove_all<_Lhs>::type Lhs;
-  typedef typename remove_all<_Rhs>::type Rhs;
-  typedef typename remove_all<Derived>::type PlainObject;
-  typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar;
-  typedef typename promote_storage_type<typename traits<Lhs>::StorageKind,
-					typename traits<Rhs>::StorageKind>::ret StorageKind;
-  typedef typename promote_index_type<typename traits<Lhs>::Index,
-				      typename traits<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
-  };
-};
-} // namespace internal
-
-template<typename Derived, typename MatrixType>
-class MatrixPowerBase
-{
-  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, RealScalar cond);
-
-    void compute(MatrixType& res, RealScalar p);
-
-    template<typename OtherDerived, typename ResultType>
-    void compute(const OtherDerived& b, ResultType& res, RealScalar p);
-    
-    Index rows() const { return m_A.rows(); }
-    Index cols() const { return m_A.cols(); }
-
-  protected:
-    typedef Array<RealScalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> RealArray;
-
-    const MatrixType& m_A;
-    const MatrixType m_Id;
-    MatrixType m_tmp1, m_tmp2;
-    RealScalar m_conditionNumber;
-};
-
-template<typename Derived, typename MatrixType>
-MatrixPowerBase<Derived,MatrixType>::MatrixPowerBase(const MatrixType& A, RealScalar cond) :
-  m_A(A),
-  m_Id(MatrixType::Identity(A.rows(),A.cols())),
-  m_conditionNumber(cond)
-{ eigen_assert(A.rows() == A.cols()); }
-
-template<typename Derived, typename MatrixType>
-void MatrixPowerBase<Derived,MatrixType>::compute(MatrixType& res, RealScalar p)
-{ static_cast<Derived*>(this)->compute(res,p); }
-
-template<typename Derived, typename MatrixType>
-template<typename OtherDerived, typename ResultType>
-void MatrixPowerBase<Derived,MatrixType>::compute(const OtherDerived& b, ResultType& res, RealScalar p)
-{ static_cast<Derived*>(this)->compute(b,res,p); }
-
-template<typename Derived, typename Lhs, typename Rhs>
-class MatrixPowerProductBase : public MatrixBase<Derived>
-{
-  public:
-    typedef MatrixBase<Derived> Base;
-    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); }
-};
-
-template<typename Derived>
-template<typename ProductDerived, typename Lhs, typename Rhs>
-Derived& MatrixBase<Derived>::lazyAssign(const MatrixPowerProductBase<ProductDerived,Lhs,Rhs>& other)
-{
-  other.derived().evalTo(derived());
-  return derived();
-}
-
 } // namespace Eigen
 
 #endif // EIGEN_MATRIX_POWER