From eaf92ef48c0037a239b1fecbd0e0bbe5c3d48fd9 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Thu, 4 Jul 2013 04:42:02 +0800 Subject: Remove unreachable MatrixPowerTriangular, paving the way to future cleanups. --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 141 --------------------- 1 file changed, 141 deletions(-) (limited to 'unsupported') 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 * -- cgit v1.2.3