Make use of TRMM (speed up), and remove useless condition (the triangular don't need LU)

2 files changed, 5 insertions, 5 deletions

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<MatrixType>::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<MatrixPower<MatrixType>,MatrixType>
 
     typedef Matrix<std::complex<RealScalar>,   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<RealScalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> 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<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType;
 
     const MatrixType& m_A;
@@ -343,7 +342,7 @@ void MatrixPowerTriangularAtomic<MatrixType>::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<MatrixType>::computeBig(MatrixType& res, RealSc
 
   for (; numberOfSquareRoots; --numberOfSquareRoots) {
     compute2x2(res, std::ldexp(p,-numberOfSquareRoots));
-    res *= res;
+    res = res.template triangularView<Upper>() * res;
   }
   compute2x2(res, p);
 }