From eaf92ef48c0037a239b1fecbd0e0bbe5c3d48fd9 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Thu, 4 Jul 2013 04:42:02 +0800 Subject: [PATCH 1/2] Remove unreachable MatrixPowerTriangular, paving the way to future cleanups. --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 141 ------------------ 1 file changed, 141 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index cdf2272ed..ffd31f6d9 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -250,147 +250,6 @@ MatrixPowerAtomic::computeSuperDiag(RealScalar curr, RealScalar prev return 2 * std::exp(p * (std::log(curr) + std::log(prev)) / 2) * std::sinh(p * w) / (curr - prev); } -/** - * \ingroup MatrixFunctions_Module - * - * \brief Class for computing matrix powers. - * - * \tparam MatrixType type of the base, expected to be an instantiation - * of the Matrix class template. - * - * This class is capable of computing upper triangular matrices raised - * to an arbitrary real power. - */ -template -class MatrixPowerTriangular -{ - private: - 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; - - public: - typedef MatrixType PlainObject; - - /** - * \brief Constructor. - * - * \param[in] A the base of the matrix power. - * - * The class stores a reference to A, so it should not be changed - * (or destroyed) before evaluation. - */ - explicit MatrixPowerTriangular(const MatrixType& A) : m_A(A), m_conditionNumber(0) - { eigen_assert(A.rows() == A.cols()); } - - /** - * \brief Returns the matrix power. - * - * \param[in] p exponent, a real scalar. - * \return The expression \f$ A^p \f$, where A is specified in the - * constructor. - */ - const MatrixPowerRetval operator()(RealScalar p) - { return MatrixPowerRetval(*this, p); } - - /** - * \brief Compute the matrix power. - * - * \param[in] p exponent, a real scalar. - * \param[out] res \f$ A^p \f$ where A is specified in the - * constructor. - */ - void compute(MatrixType& res, RealScalar p); - - Index rows() const { return m_A.rows(); } - Index cols() const { return m_A.cols(); } - - private: - typename MatrixType::Nested m_A; - MatrixType m_tmp; - RealScalar m_conditionNumber; - RealScalar modfAndInit(RealScalar, RealScalar*); - - template - void computeIntPower(ResultType&, RealScalar); - - template - void computeFracPower(ResultType&, RealScalar); -}; - -template -void MatrixPowerTriangular::compute(MatrixType& res, RealScalar p) -{ - switch (cols()) { - case 0: - break; - case 1: - res(0,0) = std::pow(m_A.coeff(0,0), p); - break; - default: - RealScalar intpart, x = modfAndInit(p, &intpart); - computeIntPower(res, intpart); - computeFracPower(res, x); - } -} - -template -typename MatrixPowerTriangular::RealScalar -MatrixPowerTriangular::modfAndInit(RealScalar x, RealScalar* intpart) -{ - typedef Array< RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > RealArray; - - *intpart = std::floor(x); - RealScalar res = x - *intpart; - - if (!m_conditionNumber && res) { - const RealArray absTdiag = m_A.diagonal().array().abs(); - m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff(); - } - - if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber, res)) { - --res; - ++*intpart; - } - return res; -} - -template -template -void MatrixPowerTriangular::computeIntPower(ResultType& res, RealScalar p) -{ - RealScalar pp = std::abs(p); - - if (p<0) m_tmp = m_A.template triangularView().solve(MatrixType::Identity(rows(), cols())); - else m_tmp = m_A.template triangularView(); - - res = MatrixType::Identity(rows(), cols()); - while (pp >= 1) { - if (std::fmod(pp, 2) >= 1) - res.template triangularView() = m_tmp.template triangularView() * res; - m_tmp.template triangularView() = m_tmp.template triangularView() * m_tmp; - pp /= 2; - } -} - -template -template -void MatrixPowerTriangular::computeFracPower(ResultType& res, RealScalar p) -{ - if (p) { - eigen_assert(m_conditionNumber); - MatrixPowerAtomic(m_A, p).compute(m_tmp); - res = m_tmp * res; - } -} - /** * \ingroup MatrixFunctions_Module * From 3cda1deb52d887218e0220f90c90d119a9b807b1 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Thu, 4 Jul 2013 05:10:43 +0800 Subject: [PATCH 2/2] Simplify class hierarchy. --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 51 ++++++++++--------- 1 file changed, 26 insertions(+), 25 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index ffd31f6d9..c32437281 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -12,14 +12,16 @@ namespace Eigen { -template -class MatrixPowerRetval : public ReturnByValue< MatrixPowerRetval > +template class MatrixPower; + +template +class MatrixPowerRetval : public ReturnByValue< MatrixPowerRetval > { public: - typedef typename MatrixPowerType::PlainObject::RealScalar RealScalar; - typedef typename MatrixPowerType::PlainObject::Index Index; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; - MatrixPowerRetval(MatrixPowerType& pow, RealScalar p) : m_pow(pow), m_p(p) + MatrixPowerRetval(MatrixPower& pow, RealScalar p) : m_pow(pow), m_p(p) { } template @@ -30,7 +32,7 @@ class MatrixPowerRetval : public ReturnByValue< MatrixPowerRetval& m_pow; const RealScalar m_p; MatrixPowerRetval& operator=(const MatrixPowerRetval&); }; @@ -47,7 +49,7 @@ class MatrixPowerAtomic typedef typename MatrixType::RealScalar RealScalar; typedef std::complex ComplexScalar; typedef typename MatrixType::Index Index; - typedef Array< Scalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > ArrayType; + typedef Array ArrayType; const MatrixType& m_A; RealScalar m_p; @@ -276,7 +278,6 @@ class MatrixPower enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, - Options = MatrixType::Options, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; @@ -285,8 +286,6 @@ class MatrixPower typedef typename MatrixType::Index Index; public: - typedef MatrixType PlainObject; - /** * \brief Constructor. * @@ -305,8 +304,8 @@ class MatrixPower * \return The expression \f$ A^p \f$, where A is specified in the * constructor. */ - const MatrixPowerRetval operator()(RealScalar p) - { return MatrixPowerRetval(*this, p); } + const MatrixPowerRetval operator()(RealScalar p) + { return MatrixPowerRetval(*this, p); } /** * \brief Compute the matrix power. @@ -315,15 +314,16 @@ class MatrixPower * \param[out] res \f$ A^p \f$ where A is specified in the * constructor. */ - void compute(MatrixType& res, RealScalar p); + template + void compute(ResultType& res, RealScalar p); Index rows() const { return m_A.rows(); } Index cols() const { return m_A.cols(); } private: typedef std::complex ComplexScalar; - typedef Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, - MaxColsAtCompileTime > ComplexMatrix; + typedef Matrix ComplexMatrix; typename MatrixType::Nested m_A; MatrixType m_tmp; @@ -338,21 +338,22 @@ class MatrixPower template void computeFracPower(ResultType&, RealScalar); - template + template static void revertSchur( - Matrix< ComplexScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res, + Matrix& res, const ComplexMatrix& T, const ComplexMatrix& U); - template + template static void revertSchur( - Matrix< RealScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res, + Matrix& res, const ComplexMatrix& T, const ComplexMatrix& U); }; template -void MatrixPower::compute(MatrixType& res, RealScalar p) +template +void MatrixPower::compute(ResultType& res, RealScalar p) { switch (cols()) { case 0: @@ -371,7 +372,7 @@ template typename MatrixPower::RealScalar MatrixPower::modfAndInit(RealScalar x, RealScalar* intpart) { - typedef Array< RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > RealArray; + typedef Array RealArray; *intpart = std::floor(x); RealScalar res = x - *intpart; @@ -423,17 +424,17 @@ void MatrixPower::computeFracPower(ResultType& res, RealScalar p) } template -template +template inline void MatrixPower::revertSchur( - Matrix< ComplexScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res, + Matrix& res, const ComplexMatrix& T, const ComplexMatrix& U) { res.noalias() = U * (T.template triangularView() * U.adjoint()); } template -template +template inline void MatrixPower::revertSchur( - Matrix< RealScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res, + Matrix& res, const ComplexMatrix& T, const ComplexMatrix& U) { res.noalias() = (U * (T.template triangularView() * U.adjoint())).real(); }