diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 520f23bff..f66600dbf 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -194,7 +194,7 @@ void MatrixPowerTriangular::computeIntPower(const Derived& b, Result else if (p>0) { m_tmp1 = m_T; } - else if (m_OKforLU && b.cols()*(pp-applyings) <= m_A.cols()*squarings) { + else if (b.cols()*(pp-applyings) <= m_A.cols()*squarings) { res = m_T.solve(b); for (--pp; pp >= 1; --pp) res = m_T.solve(res); @@ -301,6 +301,7 @@ class MatrixPower : public MatrixPowerBase,MatrixType> typedef Matrix, RowsAtCompileTime, ColsAtCompileTime, Options,MaxRowsAtCompileTime,MaxColsAtCompileTime> ComplexMatrix; + static const bool m_OKforLU = RowsAtCompileTime == Dynamic || RowsAtCompileTime > 4; ComplexMatrix m_T, m_U, m_fT; RealScalar modfAndInit(RealScalar, RealScalar*); diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h index 96846cb93..2f54a5c5b 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h @@ -24,7 +24,6 @@ namespace Eigen { typedef typename Base::RealArray RealArray; #define EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(Derived) \ - using Base::m_OKforLU; \ using Base::m_A; \ using Base::m_Id; \ using Base::m_tmp1; \ @@ -99,7 +98,6 @@ class MatrixPowerBase protected: typedef Array RealArray; - static const bool m_OKforLU = RowsAtCompileTime == Dynamic || RowsAtCompileTime > 4; const MatrixType& m_A; const MatrixType m_Id; @@ -284,6 +282,7 @@ class MatrixPowerTriangularAtomic }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; typedef Array ArrayType; const MatrixType& m_A; @@ -343,7 +342,7 @@ void MatrixPowerTriangularAtomic::compute2x2(MatrixType& res, RealSc ArrayType logTdiag = m_A.diagonal().array().log(); res.coeffRef(0,0) = pow(m_A.coeff(0,0), p); - for (int i=1; i < m_A.cols(); ++i) { + for (Index i=1; i < m_A.cols(); ++i) { res.coeffRef(i,i) = pow(m_A.coeff(i,i), p); if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) { res.coeffRef(i-1,i) = p * pow(m_A.coeff(i-1,i), p-1); @@ -392,7 +391,7 @@ void MatrixPowerTriangularAtomic::computeBig(MatrixType& res, RealSc for (; numberOfSquareRoots; --numberOfSquareRoots) { compute2x2(res, std::ldexp(p,-numberOfSquareRoots)); - res *= res; + res = res.template triangularView() * res; } compute2x2(res, p); }