Fix doc and tidy up

3 files changed, 80 insertions, 72 deletions

diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
index 54f275fde..64afdaab8 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -67,7 +67,6 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
      *
      * \param[in]  b    a matrix with the same rows as A.
      * \param[in]  p    exponent, a real scalar.
-     * \param[in]  noalias
      * \param[out] res  \f$ A^p b \f$, where A is specified in the
      * constructor.
      */
@@ -75,16 +74,8 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
     void compute(const Derived& b, ResultType& res, RealScalar p);
 
   private:
-    using Base::m_A;
-    
-    static const int Rows    = MatrixType::RowsAtCompileTime;
-    static const int Cols    = MatrixType::ColsAtCompileTime;
-    static const int Options = MatrixType::Options;
-    static const int MaxRows = MatrixType::MaxRowsAtCompileTime;
-    static const int MaxCols = MatrixType::MaxColsAtCompileTime;
-    
-    typedef Matrix<std::complex<RealScalar>,Rows,Cols,Options,MaxRows,MaxCols> ComplexMatrix;
-
+    typedef Matrix<std::complex<RealScalar>,   RowsAtCompileTime,   ColsAtCompileTime,
+				    Options,MaxRowsAtCompileTime,MaxColsAtCompileTime> ComplexMatrix;
     MatrixType m_tmp1, m_tmp2;
     ComplexMatrix m_T, m_U, m_fT;
     bool m_init;
@@ -158,7 +149,7 @@ typename MatrixPower<MatrixType>::Base::RealScalar MatrixPower<MatrixType>::modf
     m_U = schurOfA.matrixU();
     m_init = true;
 
-    const Array<RealScalar,Rows,1,ColMajor,MaxRows> absTdiag = m_T.diagonal().array().abs();
+    const RealArray absTdiag = m_T.diagonal().array().abs();
     maxAbsEival = absTdiag.maxCoeff();
     minAbsEival = absTdiag.minCoeff();
   }
@@ -251,8 +242,8 @@ void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
 {
   if (p) {
     MatrixPowerTriangularAtomic<ComplexMatrix>(m_T).compute(m_fT, p);
-    internal::recompose_complex_schur<NumTraits<Scalar>::IsComplex>::run(m_tmp1, m_fT, m_U);
-    res = m_tmp1 * res;
+    internal::recompose_complex_schur<NumTraits<Scalar>::IsComplex>::run(m_tmp2, m_fT, m_U);
+    res = m_tmp2 * res;
   }
 }
 
@@ -347,9 +338,9 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
 };
 
 namespace internal {
-template<typename MatrixType, typename Derived>
-struct nested<MatrixPowerMatrixProduct<MatrixType,Derived> >
-{ typedef typename MatrixPowerMatrixProduct<MatrixType,Derived>::PlainObject const& type; };
+template<typename Lhs, typename Rhs>
+struct nested<MatrixPowerMatrixProduct<Lhs,Rhs> >
+{ typedef typename MatrixPowerMatrixProduct<Lhs,Rhs>::PlainObject const& type; };
 
 template<typename Derived>
 struct traits<MatrixPowerReturnValue<Derived> >
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
index 0cb50694a..b5e8ed7ed 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
@@ -97,17 +97,20 @@ inline int matrix_power_get_pade_degree(long double normIminusT)
 }
 } // namespace internal
 
-template<typename MatrixType, int UpLo=Upper>
+template<typename MatrixType, unsigned int Mode=Upper>
 class MatrixPowerTriangularAtomic
 {
   private:
+    enum {
+      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
+    };
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
-    typedef Array<Scalar,
-		  EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime),
-		  1,ColMajor,
-		  EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::MaxRowsAtCompileTime,MatrixType::MaxColsAtCompileTime)> ArrayType;
+    typedef Array<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType;
+
     const MatrixType& m_T;
+    const MatrixType m_Id;
 
     void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) const;
     void compute2x2(MatrixType& res, RealScalar p) const;
@@ -118,13 +121,14 @@ class MatrixPowerTriangularAtomic
     void compute(MatrixType& res, RealScalar p) const;
 };
 
-template<typename MatrixType, int UpLo>
-MatrixPowerTriangularAtomic<MatrixType,UpLo>::MatrixPowerTriangularAtomic(const MatrixType& T) :
-  m_T(T)
-{ eigen_assert(T.rows() == T.cols()); }
+template<typename MatrixType, unsigned int Mode>
+MatrixPowerTriangularAtomic<MatrixType,Mode>::MatrixPowerTriangularAtomic(const MatrixType& T) :
+  m_T(T),
+  m_Id(MatrixType::Identity(T.rows(), T.cols()))
+{ /* empty body */ }
 
-template<typename MatrixType, int UpLo>
-void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute(MatrixType& res, RealScalar p) const
+template<typename MatrixType, unsigned int Mode>
+void MatrixPowerTriangularAtomic<MatrixType,Mode>::compute(MatrixType& res, RealScalar p) const
 {
   switch (m_T.rows()) {
     case 0:
@@ -140,21 +144,21 @@ void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute(MatrixType& res, Real
   }
 }
 
-template<typename MatrixType, int UpLo>
-void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computePade(int degree, const MatrixType& IminusT, MatrixType& res,
+template<typename MatrixType, unsigned int Mode>
+void MatrixPowerTriangularAtomic<MatrixType,Mode>::computePade(int degree, const MatrixType& IminusT, MatrixType& res,
   RealScalar p) const
 {
   int i = degree<<1;
   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();
+    res = (m_Id + res).template triangularView<Mode>().solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) :
+	(p-(i>>1))/((i-1)<<1)) * IminusT).eval();
   }
-  res += MatrixType::Identity(m_T.rows(), m_T.cols());
+  res += m_Id;
 }
 
-template<typename MatrixType, int UpLo>
-void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute2x2(MatrixType& res, RealScalar p) const
+template<typename MatrixType, unsigned int Mode>
+void MatrixPowerTriangularAtomic<MatrixType,Mode>::compute2x2(MatrixType& res, RealScalar p) const
 {
   using std::abs;
   using std::pow;
@@ -180,8 +184,8 @@ void MatrixPowerTriangularAtomic<MatrixType,UpLo>::compute2x2(MatrixType& res, R
   }
 }
 
-template<typename MatrixType, int UpLo>
-void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computeBig(MatrixType& res, RealScalar p) const
+template<typename MatrixType, unsigned int Mode>
+void MatrixPowerTriangularAtomic<MatrixType,Mode>::computeBig(MatrixType& res, RealScalar p) const
 {
   const int digits = std::numeric_limits<RealScalar>::digits;
   const RealScalar maxNormForPade = digits <=  24? 4.3386528e-1f:                           // sigle precision
@@ -217,9 +221,16 @@ void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computeBig(MatrixType& res, R
 }
 
 #define EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(Derived) \
-  typedef MatrixPowerBase<Derived<MatrixType>,MatrixType> Base; \
-  using typename Base::Scalar; \
-  using typename Base::RealScalar;
+  typedef MatrixPowerBase<Derived<MatrixType>, MatrixType> Base; \
+  using Base::RowsAtCompileTime; \
+  using Base::ColsAtCompileTime; \
+  using Base::Options; \
+  using Base::MaxRowsAtCompileTime; \
+  using Base::MaxColsAtCompileTime; \
+  typedef typename Base::Scalar     Scalar; \
+  typedef typename Base::RealScalar RealScalar; \
+  typedef typename Base::RealArray  RealArray; \
+  using Base::m_A;
 
 #define	EIGEN_MATRIX_POWER_PRODUCT_PUBLIC_INTERFACE(Derived) \
   typedef MatrixPowerProductBase<Derived, Lhs, Rhs> Base; \
@@ -254,25 +265,24 @@ struct traits<MatrixPowerProductBase<Derived,_Lhs,_Rhs> >
 template<typename Derived, typename MatrixType>
 class MatrixPowerBase
 {
-  protected:
+  public:
+    enum {
+      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+      Options = MatrixType::Options,
+      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
+      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
+    };
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
     typedef typename MatrixType::Index Index;
 
-    const MatrixType& m_A;
-    const bool m_del;  // whether to delete the pointer at destruction
-
-  public:
-    explicit MatrixPowerBase(const MatrixType& A) :
-      m_A(A),
-      m_del(false)
-    { /* empty body */ }
+    explicit MatrixPowerBase(const MatrixType& A)
+    : m_A(A), m_del(false) { }
 
     template<typename OtherDerived>
-    explicit MatrixPowerBase(const MatrixBase<OtherDerived>& A) :
-      m_A(*new MatrixType(A)),
-      m_del(true)
-    { /* empty body */ }
+    explicit MatrixPowerBase(const MatrixBase<OtherDerived>& A)
+    : m_A(*new MatrixType(A)), m_del(true) { }
 
     ~MatrixPowerBase()
     { if (m_del)  delete &m_A; }
@@ -286,6 +296,13 @@ class MatrixPowerBase
     
     Index rows() const { return m_A.rows(); }
     Index cols() const { return m_A.cols(); }
+
+  protected:
+    typedef Array<RealScalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> RealArray;
+    const MatrixType& m_A;
+
+  private:
+    const bool m_del;  // whether to delete the pointer at destruction
 };
 
 template<typename Derived, typename Lhs, typename Rhs>
diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp
index 7d90066c8..7945327c0 100644
--- a/unsupported/test/matrix_power.cpp
+++ b/unsupported/test/matrix_power.cpp
@@ -87,7 +87,7 @@ void testExponentLaws(const MatrixType& m, double tol)
 }
 
 template<typename MatrixType, typename VectorType>
-void testMatrixVectorProduct(const MatrixType& m, const VectorType& v, double tol)
+void testProduct(const MatrixType& m, const VectorType& v, double tol)
 {
   typedef typename MatrixType::RealScalar RealScalar;
   MatrixType m1;
@@ -108,9 +108,18 @@ void testMatrixVectorProduct(const MatrixType& m, const VectorType& v, double to
   }
 }
 
+template<typename MatrixType, typename VectorType>
+void testMatrixVector(const MatrixType& m, const VectorType& v, double tol)
+{
+  testExponentLaws(m,tol);
+  testProduct(m,v,tol);
+}
+
 void test_matrix_power()
 {
+  typedef Matrix<double,3,3,RowMajor>         Matrix3dRowMajor;
   typedef Matrix<long double,Dynamic,Dynamic> MatrixXe;
+  typedef Matrix<long double,Dynamic,1>       VectorXe;
 
   CALL_SUBTEST_2(test2dRotation<double>(1e-13));
   CALL_SUBTEST_1(test2dRotation<float>(2e-5));  // was 1e-5, relaxed for clang 2.8 / linux / x86-64
@@ -119,22 +128,13 @@ void test_matrix_power()
   CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
   CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14));
 
-  CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13));
-  CALL_SUBTEST_7(testExponentLaws(Matrix<double,3,3,RowMajor>(), 1e-13));
-  CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13));
-  CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 1e-13));
-  CALL_SUBTEST_1(testExponentLaws(Matrix2f(), 1e-4));
-  CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4));
-  CALL_SUBTEST_8(testExponentLaws(Matrix4f(), 1e-4));
-  CALL_SUBTEST_6(testExponentLaws(MatrixXf(8,8), 1e-4));
-
-  CALL_SUBTEST_2(testMatrixVectorProduct(Matrix2d(), Vector2d(), 1e-13));
-  CALL_SUBTEST_7(testMatrixVectorProduct(Matrix<double,3,3,RowMajor>(), Vector3d(), 1e-13));
-  CALL_SUBTEST_3(testMatrixVectorProduct(Matrix4cd(), Vector4cd(), 1e-13));
-  CALL_SUBTEST_4(testMatrixVectorProduct(MatrixXd(8,8), MatrixXd(8,2), 1e-13));
-  CALL_SUBTEST_1(testMatrixVectorProduct(Matrix2f(), Vector2f(), 1e-4));
-  CALL_SUBTEST_5(testMatrixVectorProduct(Matrix3cf(), Vector3cf(), 1e-4));
-  CALL_SUBTEST_8(testMatrixVectorProduct(Matrix4f(), Vector4f(), 1e-4));
-  CALL_SUBTEST_6(testMatrixVectorProduct(MatrixXf(8,8), VectorXf(8), 1e-4));
-  CALL_SUBTEST_9(testMatrixVectorProduct(MatrixXe(7,7), MatrixXe(7,9), 1e-13));
+  CALL_SUBTEST_2(testMatrixVector(Matrix2d(),         Vector2d(),    1e-13));
+  CALL_SUBTEST_7(testMatrixVector(Matrix3dRowMajor(), MatrixXd(3,5), 1e-13));
+  CALL_SUBTEST_3(testMatrixVector(Matrix4cd(),        Vector4cd(),   1e-13));
+  CALL_SUBTEST_4(testMatrixVector(MatrixXd(8,8),      VectorXd(8),   1e-13));
+  CALL_SUBTEST_1(testMatrixVector(Matrix2f(),         Vector2f(),    1e-4));
+  CALL_SUBTEST_5(testMatrixVector(Matrix3cf(),        Vector3cf(),   1e-4));
+  CALL_SUBTEST_8(testMatrixVector(Matrix4f(),         Vector4f(),    1e-4));
+  CALL_SUBTEST_6(testMatrixVector(MatrixXf(8,8),      VectorXf(8),   1e-4));
+  CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7),      VectorXe(7),   1e-13));
 }