From dd8034bd1c414d7b2058d020203f5c82dc5b54b6 Mon Sep 17 00:00:00 2001
From: Chen-Pang He <jdh8@ms63.hinet.net>
Date: Sat, 22 Sep 2012 17:37:14 +0800
Subject: Fix cost evaluation. (chain product for integral power)

---
 .../Eigen/src/MatrixFunctions/MatrixPower.h        | 30 +++++++++++++++-------
 .../Eigen/src/MatrixFunctions/MatrixPowerBase.h    | 23 ++++++++++-------
 2 files changed, 35 insertions(+), 18 deletions(-)

(limited to 'unsupported')

diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
index ff2e31d83..2ed319641 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -157,7 +157,7 @@ typename MatrixType::RealScalar MatrixPower<MatrixType>::modfAndInit(RealScalar
   *intpart = std::floor(x);
   RealScalar res = x - *intpart;
 
-  if (!m_init && res) { // !init && res
+  if (!m_init && res) {
     const ComplexSchur<MatrixType> schurOfA(m_A);
     m_T = schurOfA.matrixT();
     m_U = schurOfA.matrixU();
@@ -209,29 +209,41 @@ template<typename MatrixType>
 template<typename PlainObject, typename ResultType>
 void MatrixPower<MatrixType>::computeIntPower(const PlainObject& b, ResultType& res, RealScalar p)
 {
-  if (b.cols() > m_A.cols()) {
+  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 {
+  }
+  else {
     RealScalar pp = std::abs(p);
-    int cost = internal::binary_powering_cost(pp);
+    int squarings, applyings = internal::binary_powering_cost(pp, &squarings);
     bool init = false;
 
     if (p==0) {
       res = b;
       return;
     }
-    if (p<0)  m_tmp1 = m_A.inverse();
-    else      m_tmp1 = m_A;
+    else if (p>0) {
+      m_tmp1 = m_A;
+    }
+    else if (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 > m_A.cols()*cost) {
+    while (b.cols() * (pp - applyings) > m_A.cols() * squarings) {
       if (std::fmod(pp, 2) >= 1) {
 	apply(b, res, init);
-	--cost;
+	--applyings;
       }
       m_tmp1 *= m_tmp1;
-      --cost;
+      --squarings;
       pp /= 2;
     }
     for (; pp >= 1; --pp)
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
index e5a1fec31..eac43fa52 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
@@ -34,15 +34,18 @@ struct traits<MatrixPowerProductBase<Derived> > : traits<Derived>
 { };
 
 template<typename T>
-inline int binary_powering_cost(T p)
+inline int binary_powering_cost(T p, int* squarings)
 {
-  int cost, tmp;
-  frexp(p, &cost);
+  int applyings=0, tmp;
+
+  if (frexp(p, squarings) != 0.5);
+    --*squarings;
+
   while (std::frexp(p, &tmp), tmp > 0) {
     p -= std::ldexp(static_cast<T>(0.5), tmp);
-    ++cost;
+    ++applyings;
   }
-  return cost;
+  return applyings;
 }
 
 inline int matrix_power_get_pade_degree(float normIminusT)
@@ -145,7 +148,7 @@ void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computePade(int degree, const
   RealScalar p) const
 {
   int i = degree<<1;
-  res = (p-(i>>1)) / ((i-1)<<1) * IminusT;
+  res = (p-degree) / ((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();
@@ -166,9 +169,11 @@ void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute2x2(MatrixType& res, R
     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))) {
+    }
+    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 {
+    }
+    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);
@@ -187,9 +192,9 @@ void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computeBig(MatrixType& res, R
 				    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;
+  int degree, degree2, numberOfSquareRoots=0, numberOfExtraSquareRoots=0;
 
   while (true) {
     IminusT = MatrixType::Identity(m_T.rows(), m_T.cols()) - T;
-- 
cgit v1.2.3