Make MatrixPowerTriangularAtomic::computePade static because it should be.

2 files changed, 13 insertions, 19 deletions

diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
index a5e0fcaa3..f2b1a5993 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -102,7 +102,7 @@ void MatrixPowerTriangular<MatrixType>::compute(MatrixType& res, RealScalar p)
       break;
     default:
       RealScalar intpart, x = modfAndInit(p, &intpart);
-      res = m_Id;
+      res = MatrixType::Identity(m_A.rows(), m_A.cols());
       computeIntPower(res, intpart);
       computeFracPower(res, x);
   }
@@ -162,7 +162,7 @@ void MatrixPowerTriangular<MatrixType>::computeIntPower(ResultType& res, RealSca
 {
   RealScalar pp = std::abs(p);
 
-  if (p<0)  m_tmp1 = m_T.solve(m_Id);
+  if (p<0)  m_tmp1 = m_T.solve(MatrixType::Identity(m_A.rows(), m_A.cols()));
   else      m_tmp1 = m_T;
 
   while (pp >= 1) {
@@ -178,7 +178,7 @@ 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;
+    m_tmp2 = MatrixType::Identity(m_A.rows(), m_A.cols());
     computeIntPower(m_tmp2, p);
     res.noalias() = m_tmp2.template triangularView<Upper>() * b;
   }
@@ -201,7 +201,7 @@ void MatrixPowerTriangular<MatrixType>::computeIntPower(const Derived& b, Result
       return;
     }
     else {
-      m_tmp1 = m_T.solve(m_Id);
+      m_tmp1 = m_T.solve(MatrixType::Identity(m_A.rows(), m_A.cols()));
     }
 
     while (b.cols()*(pp-applyings) > m_A.cols()*squarings) {
@@ -330,7 +330,7 @@ void MatrixPower<MatrixType>::compute(MatrixType& res, RealScalar p)
       break;
     default:
       RealScalar intpart, x = modfAndInit(p, &intpart);
-      res = m_Id;
+      res = MatrixType::Identity(m_A.rows(), m_A.cols());
       computeIntPower(res, intpart);
       computeFracPower(res, x);
   }
@@ -409,7 +409,7 @@ template<typename Derived, typename ResultType>
 void MatrixPower<MatrixType>::computeIntPower(const Derived& b, ResultType& res, RealScalar p)
 {
   if (b.cols() >= m_A.cols()) {
-    m_tmp2 = m_Id;
+    m_tmp2 = MatrixType::Identity(m_A.rows(), m_A.cols());
     computeIntPower(m_tmp2, p);
     res.noalias() = m_tmp2 * b;
   }
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
index 2f54a5c5b..e26497f65 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
@@ -25,7 +25,6 @@ namespace Eigen {
 
 #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;
@@ -77,7 +76,6 @@ class MatrixPowerBase
 
     explicit MatrixPowerBase(const MatrixType& A) :
       m_A(A),
-      m_Id(MatrixType::Identity(A.rows(), A.cols())),
       m_conditionNumber(0)
     { eigen_assert(A.rows() == A.cols()); }
 
@@ -100,7 +98,6 @@ class MatrixPowerBase
     typedef Array<RealScalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> RealArray;
 
     const MatrixType& m_A;
-    const MatrixType m_Id;
     MatrixType m_tmp1, m_tmp2;
     RealScalar m_conditionNumber;
 };
@@ -286,9 +283,8 @@ class MatrixPowerTriangularAtomic
     typedef Array<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType;
 
     const MatrixType& m_A;
-    const MatrixType m_Id;
 
-    void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) const;
+    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;
 
@@ -299,8 +295,7 @@ class MatrixPowerTriangularAtomic
 
 template<typename MatrixType>
 MatrixPowerTriangularAtomic<MatrixType>::MatrixPowerTriangularAtomic(const MatrixType& T) :
-  m_A(T),
-  m_Id(MatrixType::Identity(T.rows(), T.cols()))
+  m_A(T)
 { eigen_assert(T.rows() == T.cols()); }
 
 template<typename MatrixType>
@@ -321,16 +316,15 @@ void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScala
 }
 
 template<typename MatrixType>
-void MatrixPowerTriangularAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res,
-  RealScalar p) const
+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 = (m_Id + 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()) + res).template triangularView<Upper>()
+	.solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : (p-(i>>1))/((i-1)<<1)) * IminusT).eval();
   }
-  res += m_Id;
+  res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
 }
 
 template<typename MatrixType>
@@ -374,7 +368,7 @@ void MatrixPowerTriangularAtomic<MatrixType>::computeBig(MatrixType& res, RealSc
   bool hasExtraSquareRoot = false;
 
   while (true) {
-    IminusT = m_Id - T;
+    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);