mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-01-30 17:40:05 +08:00
Remove unreachable MatrixPowerTriangular, paving the way to future cleanups.
This commit is contained in:
parent
155fa0ca83
commit
eaf92ef48c
@ -250,147 +250,6 @@ MatrixPowerAtomic<MatrixType>::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<typename MatrixType>
|
||||
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<MatrixPowerTriangular> operator()(RealScalar p)
|
||||
{ return MatrixPowerRetval<MatrixPowerTriangular>(*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<typename ResultType>
|
||||
void computeIntPower(ResultType&, RealScalar);
|
||||
|
||||
template<typename ResultType>
|
||||
void computeFracPower(ResultType&, RealScalar);
|
||||
};
|
||||
|
||||
template<typename MatrixType>
|
||||
void MatrixPowerTriangular<MatrixType>::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 MatrixType>
|
||||
typename MatrixPowerTriangular<MatrixType>::RealScalar
|
||||
MatrixPowerTriangular<MatrixType>::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<typename MatrixType>
|
||||
template<typename ResultType>
|
||||
void MatrixPowerTriangular<MatrixType>::computeIntPower(ResultType& res, RealScalar p)
|
||||
{
|
||||
RealScalar pp = std::abs(p);
|
||||
|
||||
if (p<0) m_tmp = m_A.template triangularView<Upper>().solve(MatrixType::Identity(rows(), cols()));
|
||||
else m_tmp = m_A.template triangularView<Upper>();
|
||||
|
||||
res = MatrixType::Identity(rows(), cols());
|
||||
while (pp >= 1) {
|
||||
if (std::fmod(pp, 2) >= 1)
|
||||
res.template triangularView<Upper>() = m_tmp.template triangularView<Upper>() * res;
|
||||
m_tmp.template triangularView<Upper>() = m_tmp.template triangularView<Upper>() * m_tmp;
|
||||
pp /= 2;
|
||||
}
|
||||
}
|
||||
|
||||
template<typename MatrixType>
|
||||
template<typename ResultType>
|
||||
void MatrixPowerTriangular<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
|
||||
{
|
||||
if (p) {
|
||||
eigen_assert(m_conditionNumber);
|
||||
MatrixPowerAtomic<MatrixType>(m_A, p).compute(m_tmp);
|
||||
res = m_tmp * res;
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup MatrixFunctions_Module
|
||||
*
|
||||
|
Loading…
Reference in New Issue
Block a user