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,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,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,Rows,Cols,Options,MaxRows,MaxCols> ComplexMatrix; - + typedef Matrix, 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::Base::RealScalar MatrixPower::modf m_U = schurOfA.matrixU(); m_init = true; - const Array 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::computeFracPower(ResultType& res, RealScalar p) { if (p) { MatrixPowerTriangularAtomic(m_T).compute(m_fT, p); - internal::recompose_complex_schur::IsComplex>::run(m_tmp1, m_fT, m_U); - res = m_tmp1 * res; + internal::recompose_complex_schur::IsComplex>::run(m_tmp2, m_fT, m_U); + res = m_tmp2 * res; } } @@ -347,9 +338,9 @@ class MatrixPowerReturnValue : public ReturnByValue -struct nested > -{ typedef typename MatrixPowerMatrixProduct::PlainObject const& type; }; +template +struct nested > +{ typedef typename MatrixPowerMatrixProduct::PlainObject const& type; }; template struct traits > 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 +template class MatrixPowerTriangularAtomic { private: + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime + }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - typedef Array ArrayType; + typedef Array 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 -MatrixPowerTriangularAtomic::MatrixPowerTriangularAtomic(const MatrixType& T) : - m_T(T) -{ eigen_assert(T.rows() == T.cols()); } +template +MatrixPowerTriangularAtomic::MatrixPowerTriangularAtomic(const MatrixType& T) : + m_T(T), + m_Id(MatrixType::Identity(T.rows(), T.cols())) +{ /* empty body */ } -template -void MatrixPowerTriangularAtomic::compute(MatrixType& res, RealScalar p) const +template +void MatrixPowerTriangularAtomic::compute(MatrixType& res, RealScalar p) const { switch (m_T.rows()) { case 0: @@ -140,21 +144,21 @@ void MatrixPowerTriangularAtomic::compute(MatrixType& res, Real } } -template -void MatrixPowerTriangularAtomic::computePade(int degree, const MatrixType& IminusT, MatrixType& res, +template +void MatrixPowerTriangularAtomic::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() - .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().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 -void MatrixPowerTriangularAtomic::compute2x2(MatrixType& res, RealScalar p) const +template +void MatrixPowerTriangularAtomic::compute2x2(MatrixType& res, RealScalar p) const { using std::abs; using std::pow; @@ -180,8 +184,8 @@ void MatrixPowerTriangularAtomic::compute2x2(MatrixType& res, R } } -template -void MatrixPowerTriangularAtomic::computeBig(MatrixType& res, RealScalar p) const +template +void MatrixPowerTriangularAtomic::computeBig(MatrixType& res, RealScalar p) const { const int digits = std::numeric_limits::digits; const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision @@ -217,9 +221,16 @@ void MatrixPowerTriangularAtomic::computeBig(MatrixType& res, R } #define EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(Derived) \ - typedef MatrixPowerBase,MatrixType> Base; \ - using typename Base::Scalar; \ - using typename Base::RealScalar; + typedef MatrixPowerBase, 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 Base; \ @@ -254,25 +265,24 @@ struct traits > template 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 - explicit MatrixPowerBase(const MatrixBase& A) : - m_A(*new MatrixType(A)), - m_del(true) - { /* empty body */ } + explicit MatrixPowerBase(const MatrixBase& 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 RealArray; + const MatrixType& m_A; + + private: + const bool m_del; // whether to delete the pointer at destruction }; template 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 -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 +void testMatrixVector(const MatrixType& m, const VectorType& v, double tol) +{ + testExponentLaws(m,tol); + testProduct(m,v,tol); +} + void test_matrix_power() { + typedef Matrix Matrix3dRowMajor; typedef Matrix MatrixXe; + typedef Matrix VectorXe; CALL_SUBTEST_2(test2dRotation(1e-13)); CALL_SUBTEST_1(test2dRotation(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(1e-5)); CALL_SUBTEST_9(test2dHyperbolicRotation(1e-14)); - CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13)); - CALL_SUBTEST_7(testExponentLaws(Matrix(), 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(), 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)); }