Optimize MatrixPower::computeIntPower

This commit is contained in:
Chen-Pang He 2013-07-20 18:47:54 +08:00
parent 2320073e41
commit dda869051d

View File

@ -447,6 +447,8 @@ void MatrixPower<MatrixType>::compute(ResultType& res, RealScalar p)
default: default:
RealScalar intpart; RealScalar intpart;
split(p, intpart); split(p, intpart);
res = MatrixType::Identity(rows(), cols());
computeIntPower(res, intpart); computeIntPower(res, intpart);
if (p) computeFracPower(res, p); if (p) computeFracPower(res, p);
} }
@ -514,15 +516,18 @@ void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p)
{ {
RealScalar pp = std::abs(p); RealScalar pp = std::abs(p);
if (p<0) m_tmp = m_A.inverse(); if (p<0)
else m_tmp = m_A; m_tmp = m_A.inverse();
else
m_tmp = m_A;
res = MatrixType::Identity(rows(), cols()); while (true) {
while (pp >= 1) {
if (std::fmod(pp, 2) >= 1) if (std::fmod(pp, 2) >= 1)
res = m_tmp * res; res = m_tmp * res;
m_tmp *= m_tmp;
pp /= 2; pp /= 2;
if (pp < 1)
break;
m_tmp *= m_tmp;
} }
} }