From dda869051db084d22dcfc46b7732b04d77ff9b4b Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sat, 20 Jul 2013 18:47:54 +0800 Subject: [PATCH] Optimize MatrixPower::computeIntPower --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index be13095e5..7124874f7 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -447,6 +447,8 @@ void MatrixPower::compute(ResultType& res, RealScalar p) default: RealScalar intpart; split(p, intpart); + + res = MatrixType::Identity(rows(), cols()); computeIntPower(res, intpart); if (p) computeFracPower(res, p); } @@ -514,15 +516,18 @@ void MatrixPower::computeIntPower(ResultType& res, RealScalar p) { RealScalar pp = std::abs(p); - if (p<0) m_tmp = m_A.inverse(); - else m_tmp = m_A; + if (p<0) + m_tmp = m_A.inverse(); + else + m_tmp = m_A; - res = MatrixType::Identity(rows(), cols()); - while (pp >= 1) { + while (true) { if (std::fmod(pp, 2) >= 1) res = m_tmp * res; - m_tmp *= m_tmp; pp /= 2; + if (pp < 1) + break; + m_tmp *= m_tmp; } }