GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-22 01:29:35 +08:00
Code Issues Packages Projects Releases Wiki Activity

Remove mat.pow * vec specialization, which causes segfault for mat.pow * mat.pow

Browse Source
This commit is contained in:
Chen-Pang He 2013-06-24 23:56:17 +08:00
parent ee8a28fb85
commit b9fc9d8f32
4 changed files with 367 additions and 715 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
unsupported/Eigen/MatrixFunctions
View File

@ -59,7 +59,6 @@
#include "src/MatrixFunctions/MatrixFunction.h" #include "src/MatrixFunctions/MatrixFunction.h"
#include "src/MatrixFunctions/MatrixSquareRoot.h" #include "src/MatrixFunctions/MatrixSquareRoot.h"
#include "src/MatrixFunctions/MatrixLogarithm.h" #include "src/MatrixFunctions/MatrixLogarithm.h"
#include "src/MatrixFunctions/MatrixPowerBase.h"
#include "src/MatrixFunctions/MatrixPower.h" #include "src/MatrixFunctions/MatrixPower.h"

600
unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library // This file is part of Eigen, a lightweight C++ template library
// for linear algebra. // for linear algebra.
// //
// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net> // Copyright (C) 2012, 2013 Chen-Pang He <jdh8@ms63.hinet.net>
// //
// This Source Code Form is subject to the terms of the Mozilla // This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed // Public License v. 2.0. If a copy of the MPL was not distributed
@ -12,6 +12,248 @@
namespace Eigen { namespace Eigen {
namespace MatrixPowerHelper {
template<typename MatrixPowerType>
class ReturnValue : public ReturnByValue< ReturnValue<MatrixPowerType> >
{
public:
typedef typename MatrixPowerType::PlainObject::RealScalar RealScalar;
typedef typename MatrixPowerType::PlainObject::Index Index;
ReturnValue(MatrixPowerType& pow, RealScalar p) : m_pow(pow), m_p(p)
{ }
template<typename ResultType>
inline void evalTo(ResultType& res) const
{ m_pow.compute(res, m_p); }
Index rows() const { return m_pow.rows(); }
Index cols() const { return m_pow.cols(); }
private:
MatrixPowerType& m_pow;
const RealScalar m_p;
ReturnValue& operator=(const ReturnValue&);
};
} // namespace MatrixPowerHelper
template<typename MatrixType>
class MatrixPowerAtomic
{
private:
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef std::complex<RealScalar> ComplexScalar;
typedef typename MatrixType::Index Index;
typedef Array< Scalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > ArrayType;
const MatrixType& m_A;
RealScalar m_p;
void computePade(int degree, const MatrixType& IminusT, MatrixType& res) const;
void compute2x2(MatrixType& res, RealScalar p) const;
void computeBig(MatrixType& res) const;
static int getPadeDegree(float normIminusT);
static int getPadeDegree(double normIminusT);
static int getPadeDegree(long double normIminusT);
static ComplexScalar computeSuperDiag(const ComplexScalar&, const ComplexScalar&, RealScalar p);
static RealScalar computeSuperDiag(RealScalar, RealScalar, RealScalar p);
public:
MatrixPowerAtomic(const MatrixType& T, RealScalar p);
void compute(MatrixType& res) const;
};
template<typename MatrixType>
MatrixPowerAtomic<MatrixType>::MatrixPowerAtomic(const MatrixType& T, RealScalar p) :
m_A(T), m_p(p)
{ eigen_assert(T.rows() == T.cols()); }
template<typename MatrixType>
void MatrixPowerAtomic<MatrixType>::compute(MatrixType& res) const
{
res.resizeLike(m_A);
switch (m_A.rows()) {
case 0:
break;
case 1:
res(0,0) = std::pow(m_A(0,0), m_p);
break;
case 2:
compute2x2(res, m_p);
break;
default:
computeBig(res);
}
}
template<typename MatrixType>
void MatrixPowerAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res) const
{
int i = degree<<1;
res = (m_p-degree) / ((i-1)<<1) * IminusT;
for (--i; i; --i) {
res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView<Upper>()
.solve((i==1 ? -m_p : i&1 ? (-m_p-(i>>1))/(i<<1) : (m_p-(i>>1))/((i-1)<<1)) * IminusT).eval();
}
res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
}
// This function assumes that res has the correct size (see bug 614)
template<typename MatrixType>
void MatrixPowerAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
{
using std::abs;
using std::pow;
ArrayType logTdiag = m_A.diagonal().array().log();
res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
for (Index i=1; i < m_A.cols(); ++i) {
res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i))
res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1)))
res.coeffRef(i-1,i) = (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));
else
res.coeffRef(i-1,i) = computeSuperDiag(m_A.coeff(i,i), m_A.coeff(i-1,i-1), p);
res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
}
}
template<typename MatrixType>
void MatrixPowerAtomic<MatrixType>::computeBig(MatrixType& res) const
{
const int digits = std::numeric_limits<RealScalar>::digits;
const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision
digits <= 53? 2.789358995219730e-1: // double precision
digits <= 64? 2.4471944416607995472e-1L: // extended precision
digits <= 106? 1.1016843812851143391275867258512e-1L: // double-double
9.134603732914548552537150753385375e-2L; // quadruple precision
MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
RealScalar normIminusT;
int degree, degree2, numberOfSquareRoots = 0;
bool hasExtraSquareRoot = false;
/* FIXME
* For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite
* loop. We should move 0 eigenvalues to bottom right corner. We need not
* worry about tiny values (e.g. 1e-300) because they will reach 1 if
* repetitively sqrt'ed.
*
* If the 0 eigenvalues are semisimple, they can form a 0 matrix at the
* bottom right corner.
*
* [ T A ]^p [ T^p (T^-1 T^p A) ]
* [ ] = [ ]
* [ 0 0 ] [ 0 0 ]
*/
for (Index i=0; i < m_A.cols(); ++i)
eigen_assert(m_A(i,i) != RealScalar(0));
while (true) {
IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
if (normIminusT < maxNormForPade) {
degree = getPadeDegree(normIminusT);
degree2 = getPadeDegree(normIminusT/2);
if (degree - degree2 <= 1 || hasExtraSquareRoot)
break;
hasExtraSquareRoot = true;
}
MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
T = sqrtT.template triangularView<Upper>();
++numberOfSquareRoots;
}
computePade(degree, IminusT, res);
for (; numberOfSquareRoots; --numberOfSquareRoots) {
compute2x2(res, std::ldexp(m_p, -numberOfSquareRoots));
res = res.template triangularView<Upper>() * res;
}
compute2x2(res, m_p);
}
template<typename MatrixType>
inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(float normIminusT)
{
const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f };
int degree = 3;
for (; degree <= 4; ++degree)
if (normIminusT <= maxNormForPade[degree - 3])
break;
return degree;
}
template<typename MatrixType>
inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(double normIminusT)
{
const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1,
1.999045567181744e-1, 2.789358995219730e-1 };
int degree = 3;
for (; degree <= 7; ++degree)
if (normIminusT <= maxNormForPade[degree - 3])
break;
return degree;
}
template<typename MatrixType>
inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(long double normIminusT)
{
#if LDBL_MANT_DIG == 53
const int maxPadeDegree = 7;
const double maxNormForPade[] = { 1.884160592658218e-2L /* degree = 3 */ , 6.038881904059573e-2L, 1.239917516308172e-1L,
1.999045567181744e-1L, 2.789358995219730e-1L };
#elif LDBL_MANT_DIG <= 64
const int maxPadeDegree = 8;
const double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L,
6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L };
#elif LDBL_MANT_DIG <= 106
const int maxPadeDegree = 10;
const double maxNormForPade[] = { 1.0007161601787493236741409687186e-4L /* degree = 3 */ ,
1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L,
2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L,
1.1016843812851143391275867258512e-1L };
#else
const int maxPadeDegree = 10;
const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ ,
6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L,
9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L,
3.908166513900489428442993794761185e-2L, 6.266780814639442865832535460550138e-2L,
9.134603732914548552537150753385375e-2L };
#endif
int degree = 3;
for (; degree <= maxPadeDegree; ++degree)
if (normIminusT <= maxNormForPade[degree - 3])
break;
return degree;
}
template<typename MatrixType>
inline typename MatrixPowerAtomic<MatrixType>::ComplexScalar
MatrixPowerAtomic<MatrixType>::computeSuperDiag(const ComplexScalar& curr, const ComplexScalar& prev, RealScalar p)
{
ComplexScalar logCurr = std::log(curr);
ComplexScalar logPrev = std::log(prev);
int unwindingNumber = std::ceil((numext::imag(logCurr - logPrev) - M_PI) / (2*M_PI));
ComplexScalar w = numext::atanh2(curr - prev, curr + prev) + ComplexScalar(0, M_PI*unwindingNumber);
return RealScalar(2) * std::exp(RealScalar(0.5) * p * (logCurr + logPrev)) * std::sinh(p * w) / (curr - prev);
}
template<typename MatrixType>
inline typename MatrixPowerAtomic<MatrixType>::RealScalar
MatrixPowerAtomic<MatrixType>::computeSuperDiag(RealScalar curr, RealScalar prev, RealScalar p)
{
RealScalar w = numext::atanh2(curr - prev, curr + prev);
return 2 * std::exp(p * (std::log(curr) + std::log(prev)) / 2) * std::sinh(p * w) / (curr - prev);
}
/** /**
* \ingroup MatrixFunctions_Module * \ingroup MatrixFunctions_Module
* *
@ -24,10 +266,22 @@ namespace Eigen {
* to an arbitrary real power. * to an arbitrary real power.
*/ */
template<typename MatrixType> template<typename MatrixType>
class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<MatrixType>,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: public:
EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(MatrixPowerTriangular) typedef MatrixType PlainObject;
/** /**
* \brief Constructor. * \brief Constructor.
@ -37,10 +291,9 @@ class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<Matri
* The class stores a reference to A, so it should not be changed * The class stores a reference to A, so it should not be changed
* (or destroyed) before evaluation. * (or destroyed) before evaluation.
*/ */
explicit MatrixPowerTriangular(const MatrixType& A) : Base(A), m_T(Base::m_A) explicit MatrixPowerTriangular(const MatrixType& A) : m_A(A), m_conditionNumber(0)
{ } { eigen_assert(A.rows() == A.cols()); }
#ifdef EIGEN_PARSED_BY_DOXYGEN
/** /**
* \brief Returns the matrix power. * \brief Returns the matrix power.
* *
@ -48,8 +301,8 @@ class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<Matri
* \return The expression \f$ A^p \f$, where A is specified in the * \return The expression \f$ A^p \f$, where A is specified in the
* constructor. * constructor.
*/ */
const MatrixPowerBaseReturnValue<MatrixPowerTriangular<MatrixType>,MatrixType> operator()(RealScalar p); const MatrixPowerHelper::ReturnValue<MatrixPowerTriangular> operator()(RealScalar p)
#endif { return MatrixPowerHelper::ReturnValue<MatrixPowerTriangular>(*this, p); }
/** /**
* \brief Compute the matrix power. * \brief Compute the matrix power.
@ -60,33 +313,18 @@ class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<Matri
*/ */
void compute(MatrixType& res, RealScalar p); void compute(MatrixType& res, RealScalar p);
/** Index rows() const { return m_A.rows(); }
* \brief Compute the matrix power multiplied by another matrix. Index cols() const { return m_A.cols(); }
*
* \param[in] b a matrix with the same rows as A.
* \param[in] p exponent, a real scalar.
* \param[out] res \f$ A^p b \f$, where A is specified in the
* constructor.
*/
template<typename Derived, typename ResultType>
void compute(const Derived& b, ResultType& res, RealScalar p);
private: private:
EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(MatrixPowerTriangular) typename MatrixType::Nested m_A;
MatrixType m_tmp;
const TriangularView<MatrixType,Upper> m_T; RealScalar m_conditionNumber;
RealScalar modfAndInit(RealScalar, RealScalar*); RealScalar modfAndInit(RealScalar, RealScalar*);
template<typename Derived, typename ResultType>
void apply(const Derived&, ResultType&, bool&);
template<typename ResultType> template<typename ResultType>
void computeIntPower(ResultType&, RealScalar); void computeIntPower(ResultType&, RealScalar);
template<typename Derived, typename ResultType>
void computeIntPower(const Derived&, ResultType&, RealScalar);
template<typename ResultType> template<typename ResultType>
void computeFracPower(ResultType&, RealScalar); void computeFracPower(ResultType&, RealScalar);
}; };
@ -94,41 +332,25 @@ class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<Matri
template<typename MatrixType> template<typename MatrixType>
void MatrixPowerTriangular<MatrixType>::compute(MatrixType& res, RealScalar p) void MatrixPowerTriangular<MatrixType>::compute(MatrixType& res, RealScalar p)
{ {
switch (m_A.cols()) { switch (cols()) {
case 0: case 0:
break; break;
case 1: case 1:
res(0,0) = std::pow(m_T.coeff(0,0), p); res(0,0) = std::pow(m_A.coeff(0,0), p);
break; break;
default: default:
RealScalar intpart, x = modfAndInit(p, &intpart); RealScalar intpart, x = modfAndInit(p, &intpart);
res = MatrixType::Identity(m_A.rows(), m_A.cols());
computeIntPower(res, intpart); computeIntPower(res, intpart);
computeFracPower(res, x); computeFracPower(res, x);
} }
} }
template<typename MatrixType>
template<typename Derived, typename ResultType>
void MatrixPowerTriangular<MatrixType>::compute(const Derived& b, ResultType& res, RealScalar p)
{
switch (m_A.cols()) {
case 0:
break;
case 1:
res = std::pow(m_T.coeff(0,0), p) * b;
break;
default:
RealScalar intpart, x = modfAndInit(p, &intpart);
computeIntPower(b, res, intpart);
computeFracPower(res, x);
}
}
template<typename MatrixType> template<typename MatrixType>
typename MatrixPowerTriangular<MatrixType>::RealScalar typename MatrixPowerTriangular<MatrixType>::RealScalar
MatrixPowerTriangular<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart) MatrixPowerTriangular<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
{ {
typedef Array< RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > RealArray;
*intpart = std::floor(x); *intpart = std::floor(x);
RealScalar res = x - *intpart; RealScalar res = x - *intpart;
@ -137,95 +359,39 @@ MatrixPowerTriangular<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart
m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff(); m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff();
} }
if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber,res)) { if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber, res)) {
--res; --res;
++*intpart; ++*intpart;
} }
return res; return res;
} }
template<typename MatrixType>
template<typename Derived, typename ResultType>
void MatrixPowerTriangular<MatrixType>::apply(const Derived& b, ResultType& res, bool& init)
{
if (init)
res = m_tmp1.template triangularView<Upper>() * res;
else {
init = true;
res.noalias() = m_tmp1.template triangularView<Upper>() * b;
}
}
template<typename MatrixType> template<typename MatrixType>
template<typename ResultType> template<typename ResultType>
void MatrixPowerTriangular<MatrixType>::computeIntPower(ResultType& res, RealScalar p) void MatrixPowerTriangular<MatrixType>::computeIntPower(ResultType& res, RealScalar p)
{ {
RealScalar pp = std::abs(p); RealScalar pp = std::abs(p);
if (p<0) m_tmp1 = m_T.solve(MatrixType::Identity(m_A.rows(), m_A.cols())); if (p<0) m_tmp = m_A.template triangularView<Upper>().solve(MatrixType::Identity(rows(), cols()));
else m_tmp1 = m_T; else m_tmp = m_A.template triangularView<Upper>();
res = MatrixType::Identity(rows(), cols());
while (pp >= 1) { while (pp >= 1) {
if (std::fmod(pp, 2) >= 1) if (std::fmod(pp, 2) >= 1)
res = m_tmp1.template triangularView<Upper>() * res; res.template triangularView<Upper>() = m_tmp.template triangularView<Upper>() * res;
m_tmp1 = m_tmp1.template triangularView<Upper>() * m_tmp1; m_tmp.template triangularView<Upper>() = m_tmp.template triangularView<Upper>() * m_tmp;
pp /= 2; pp /= 2;
} }
} }
template<typename MatrixType>
template<typename Derived, typename ResultType>
void MatrixPowerTriangular<MatrixType>::computeIntPower(const Derived& b, ResultType& res, RealScalar p)
{
if (b.cols() >= m_A.cols()) {
m_tmp2 = MatrixType::Identity(m_A.rows(), m_A.cols());
computeIntPower(m_tmp2, p);
res.noalias() = m_tmp2.template triangularView<Upper>() * b;
}
else {
RealScalar pp = std::abs(p);
int squarings, applyings = internal::binary_powering_cost(pp, &squarings);
bool init = false;
if (p==0) {
res = b;
return;
}
else if (p>0) {
m_tmp1 = m_T;
}
else if (b.cols()*(pp-applyings) <= m_A.cols()*squarings) {
res = m_T.solve(b);
for (--pp; pp >= 1; --pp)
res = m_T.solve(res);
return;
}
else {
m_tmp1 = m_T.solve(MatrixType::Identity(m_A.rows(), m_A.cols()));
}
while (b.cols()*(pp-applyings) > m_A.cols()*squarings) {
if (std::fmod(pp, 2) >= 1) {
apply(b, res, init);
--applyings;
}
m_tmp1 = m_tmp1.template triangularView<Upper>() * m_tmp1;
--squarings;
pp /= 2;
}
for (; pp >= 1; --pp)
apply(b, res, init);
}
}
template<typename MatrixType> template<typename MatrixType>
template<typename ResultType> template<typename ResultType>
void MatrixPowerTriangular<MatrixType>::computeFracPower(ResultType& res, RealScalar p) void MatrixPowerTriangular<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
{ {
if (p) { if (p) {
eigen_assert(m_conditionNumber); eigen_assert(m_conditionNumber);
MatrixPowerTriangularAtomic<MatrixType>(m_A).compute(m_tmp1, p); MatrixPowerAtomic<MatrixType>(m_A, p).compute(m_tmp);
res = m_tmp1.template triangularView<Upper>() * res; res = m_tmp * res;
} }
} }
@ -249,10 +415,22 @@ void MatrixPowerTriangular<MatrixType>::computeFracPower(ResultType& res, RealSc
* Output: \verbinclude MatrixPower_optimal.out * Output: \verbinclude MatrixPower_optimal.out
*/ */
template<typename MatrixType> template<typename MatrixType>
class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType> class MatrixPower
{ {
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: public:
EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(MatrixPower) typedef MatrixType PlainObject;
/** /**
* \brief Constructor. * \brief Constructor.
@ -262,10 +440,9 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
* The class stores a reference to A, so it should not be changed * The class stores a reference to A, so it should not be changed
* (or destroyed) before evaluation. * (or destroyed) before evaluation.
*/ */
explicit MatrixPower(const MatrixType& A) : Base(A) explicit MatrixPower(const MatrixType& A) : m_A(A), m_conditionNumber(0)
{ } { eigen_assert(A.rows() == A.cols()); }
#ifdef EIGEN_PARSED_BY_DOXYGEN
/** /**
* \brief Returns the matrix power. * \brief Returns the matrix power.
* *
@ -273,8 +450,8 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
* \return The expression \f$ A^p \f$, where A is specified in the * \return The expression \f$ A^p \f$, where A is specified in the
* constructor. * constructor.
*/ */
const MatrixPowerBaseReturnValue<MatrixPower<MatrixType>,MatrixType> operator()(RealScalar p); const MatrixPowerHelper::ReturnValue<MatrixPower> operator()(RealScalar p)
#endif { return MatrixPowerHelper::ReturnValue<MatrixPower>(*this, p); }
/** /**
* \brief Compute the matrix power. * \brief Compute the matrix power.
@ -285,44 +462,44 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
*/ */
void compute(MatrixType& res, RealScalar p); void compute(MatrixType& res, RealScalar p);
/** Index rows() const { return m_A.rows(); }
* \brief Compute the matrix power multiplied by another matrix. Index cols() const { return m_A.cols(); }
*
* \param[in] b a matrix with the same rows as A.
* \param[in] p exponent, a real scalar.
* \param[out] res \f$ A^p b \f$, where A is specified in the
* constructor.
*/
template<typename Derived, typename ResultType>
void compute(const Derived& b, ResultType& res, RealScalar p);
private: private:
EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(MatrixPower) typedef std::complex<RealScalar> ComplexScalar;
typedef Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime,
MaxColsAtCompileTime > ComplexMatrix;
typedef Matrix<std::complex<RealScalar>, RowsAtCompileTime, ColsAtCompileTime, typename MatrixType::Nested m_A;
Options,MaxRowsAtCompileTime,MaxColsAtCompileTime> ComplexMatrix; MatrixType m_tmp;
static const bool m_OKforLU = RowsAtCompileTime == Dynamic || RowsAtCompileTime > 4;
ComplexMatrix m_T, m_U, m_fT; ComplexMatrix m_T, m_U, m_fT;
RealScalar m_conditionNumber;
RealScalar modfAndInit(RealScalar, RealScalar*); RealScalar modfAndInit(RealScalar, RealScalar*);
template<typename Derived, typename ResultType>
void apply(const Derived&, ResultType&, bool&);
template<typename ResultType> template<typename ResultType>
void computeIntPower(ResultType&, RealScalar); void computeIntPower(ResultType&, RealScalar);
template<typename Derived, typename ResultType>
void computeIntPower(const Derived&, ResultType&, RealScalar);
template<typename ResultType> template<typename ResultType>
void computeFracPower(ResultType&, RealScalar); void computeFracPower(ResultType&, RealScalar);
template<int Rows, int Cols, int Opt, int MaxRows, int MaxCols>
static void revertSchur(
Matrix< ComplexScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res,
const ComplexMatrix& T,
const ComplexMatrix& U);
template<int Rows, int Cols, int Opt, int MaxRows, int MaxCols>
static void revertSchur(
Matrix< RealScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res,
const ComplexMatrix& T,
const ComplexMatrix& U);
}; };
template<typename MatrixType> template<typename MatrixType>
void MatrixPower<MatrixType>::compute(MatrixType& res, RealScalar p) void MatrixPower<MatrixType>::compute(MatrixType& res, RealScalar p)
{ {
switch (m_A.cols()) { switch (cols()) {
case 0: case 0:
break; break;
case 1: case 1:
@ -330,32 +507,17 @@ void MatrixPower<MatrixType>::compute(MatrixType& res, RealScalar p)
break; break;
default: default:
RealScalar intpart, x = modfAndInit(p, &intpart); RealScalar intpart, x = modfAndInit(p, &intpart);
res = MatrixType::Identity(m_A.rows(), m_A.cols());
computeIntPower(res, intpart); computeIntPower(res, intpart);
computeFracPower(res, x); computeFracPower(res, x);
} }
} }
template<typename MatrixType> template<typename MatrixType>
template<typename Derived, typename ResultType> typename MatrixPower<MatrixType>::RealScalar
void MatrixPower<MatrixType>::compute(const Derived& b, ResultType& res, RealScalar p) MatrixPower<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
{ {
switch (m_A.cols()) { typedef Array< RealScalar, RowsAtCompileTime, 1, ColMajor, MaxRowsAtCompileTime > RealArray;
case 0:
break;
case 1:
res = std::pow(m_A.coeff(0,0), p) * b;
break;
default:
RealScalar intpart, x = modfAndInit(p, &intpart);
computeIntPower(b, res, intpart);
computeFracPower(res, x);
}
}
template<typename MatrixType>
typename MatrixPower<MatrixType>::RealScalar MatrixPower<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
{
*intpart = std::floor(x); *intpart = std::floor(x);
RealScalar res = x - *intpart; RealScalar res = x - *intpart;
@ -375,100 +537,51 @@ typename MatrixPower<MatrixType>::RealScalar MatrixPower<MatrixType>::modfAndIni
return res; return res;
} }
template<typename MatrixType>
template<typename Derived, typename ResultType>
void MatrixPower<MatrixType>::apply(const Derived& b, ResultType& res, bool& init)
{
if (init)
res = m_tmp1 * res;
else {
init = true;
res.noalias() = m_tmp1 * b;
}
}
template<typename MatrixType> template<typename MatrixType>
template<typename ResultType> template<typename ResultType>
void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p) void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p)
{ {
RealScalar pp = std::abs(p); RealScalar pp = std::abs(p);
if (p<0) m_tmp1 = m_A.inverse(); if (p<0) m_tmp = m_A.inverse();
else m_tmp1 = m_A; else m_tmp = m_A;
res = MatrixType::Identity(rows(), cols());
while (pp >= 1) { while (pp >= 1) {
if (std::fmod(pp, 2) >= 1) if (std::fmod(pp, 2) >= 1)
res = m_tmp1 * res; res = m_tmp * res;
m_tmp1 *= m_tmp1; m_tmp *= m_tmp;
pp /= 2; pp /= 2;
} }
} }
template<typename MatrixType>
template<typename Derived, typename ResultType>
void MatrixPower<MatrixType>::computeIntPower(const Derived& b, ResultType& res, RealScalar p)
{
if (b.cols() >= m_A.cols()) {
m_tmp2 = MatrixType::Identity(m_A.rows(), m_A.cols());
computeIntPower(m_tmp2, p);
res.noalias() = m_tmp2 * b;
}
else {
RealScalar pp = std::abs(p);
int squarings, applyings = internal::binary_powering_cost(pp, &squarings);
bool init = false;
if (p==0) {
res = b;
return;
}
else if (p>0) {
m_tmp1 = m_A;
}
else if (m_OKforLU && b.cols()*(pp-applyings) <= m_A.cols()*squarings) {
PartialPivLU<MatrixType> A(m_A);
res = A.solve(b);
for (--pp; pp >= 1; --pp)
res = A.solve(res);
return;
}
else {
m_tmp1 = m_A.inverse();
}
while (b.cols()*(pp-applyings) > m_A.cols()*squarings) {
if (std::fmod(pp, 2) >= 1) {
apply(b, res, init);
--applyings;
}
m_tmp1 *= m_tmp1;
--squarings;
pp /= 2;
}
for (; pp >= 1; --pp)
apply(b, res, init);
}
}
template<typename MatrixType> template<typename MatrixType>
template<typename ResultType> template<typename ResultType>
void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p) void MatrixPower<MatrixType>::computeFracPower(ResultType& res, RealScalar p)
{ {
if (p) { if (p) {
eigen_assert(m_conditionNumber); eigen_assert(m_conditionNumber);
MatrixPowerTriangularAtomic<ComplexMatrix>(m_T).compute(m_fT, p); MatrixPowerAtomic<ComplexMatrix>(m_T, p).compute(m_fT);
internal::recompose_complex_schur<NumTraits<Scalar>::IsComplex>::run(m_tmp1, m_fT, m_U); revertSchur(m_tmp, m_fT, m_U);
res = m_tmp1 * res; res = m_tmp * res;
} }
} }
namespace internal { template<typename MatrixType>
template<int Rows, int Cols, int Opt, int MaxRows, int MaxCols>
inline void MatrixPower<MatrixType>::revertSchur(
Matrix< ComplexScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res,
const ComplexMatrix& T,
const ComplexMatrix& U)
{ res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint()); }
template<typename Derived> template<typename MatrixType>
struct traits<MatrixPowerReturnValue<Derived> > template<int Rows, int Cols, int Opt, int MaxRows, int MaxCols>
{ typedef typename Derived::PlainObject ReturnType; }; inline void MatrixPower<MatrixType>::revertSchur(
Matrix< RealScalar, Rows, Cols, Opt, MaxRows, MaxCols >& res,
} // namespace internal const ComplexMatrix& T,
const ComplexMatrix& U)
{ res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }
/** /**
* \ingroup MatrixFunctions_Module * \ingroup MatrixFunctions_Module
@ -484,7 +597,7 @@ struct traits<MatrixPowerReturnValue<Derived> >
* time this is the only way it is used. * time this is the only way it is used.
*/ */
template<typename Derived> template<typename Derived>
class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Derived> > class MatrixPowerReturnValue : public ReturnByValue< MatrixPowerReturnValue<Derived> >
{ {
public: public:
typedef typename Derived::PlainObject PlainObject; typedef typename Derived::PlainObject PlainObject;
@ -510,21 +623,6 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
inline void evalTo(ResultType& res) const inline void evalTo(ResultType& res) const
{ MatrixPower<PlainObject>(m_A.eval()).compute(res, m_p); } { MatrixPower<PlainObject>(m_A.eval()).compute(res, m_p); }
/**
* \brief Return the expression \f$ A^p b \f$.
*
* \p A and \p p are specified in the constructor.
*
* \param[in] b the matrix (expression) to be applied.
*/
template<typename OtherDerived>
const MatrixPowerProduct<MatrixPower<PlainObject>,PlainObject,OtherDerived>
operator*(const MatrixBase<OtherDerived>& b) const
{
MatrixPower<PlainObject> Apow(m_A.eval());
return MatrixPowerProduct<MatrixPower<PlainObject>,PlainObject,OtherDerived>(Apow, b.derived(), m_p);
}
Index rows() const { return m_A.rows(); } Index rows() const { return m_A.rows(); }
Index cols() const { return m_A.cols(); } Index cols() const { return m_A.cols(); }
@ -534,6 +632,18 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&); MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&);
}; };
namespace internal {
template<typename MatrixPowerType>
struct traits< MatrixPowerHelper::ReturnValue<MatrixPowerType> >
{ typedef typename MatrixPowerType::PlainObject ReturnType; };
template<typename Derived>
struct traits< MatrixPowerReturnValue<Derived> >
{ typedef typename Derived::PlainObject ReturnType; };
}
template<typename Derived> template<typename Derived>
const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(const RealScalar& p) const const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(const RealScalar& p) const
{ return MatrixPowerReturnValue<Derived>(derived(), p); } { return MatrixPowerReturnValue<Derived>(derived(), p); }

404
unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
View File

@ -1,404 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef EIGEN_MATRIX_POWER_BASE
#define EIGEN_MATRIX_POWER_BASE
namespace Eigen {
#define EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(Derived) \
typedef MatrixPowerBase<Derived, 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;
#define EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(Derived) \
using Base::m_A; \
using Base::m_tmp1; \
using Base::m_tmp2; \
using Base::m_conditionNumber;
template<typename Derived, typename MatrixType>
class MatrixPowerBaseReturnValue : public ReturnByValue<MatrixPowerBaseReturnValue<Derived,MatrixType> >
{
public:
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
MatrixPowerBaseReturnValue(Derived& pow, RealScalar p) : m_pow(pow), m_p(p)
{ }
template<typename ResultType>
inline void evalTo(ResultType& res) const
{ m_pow.compute(res, m_p); }
template<typename OtherDerived>
const MatrixPowerProduct<Derived,MatrixType,OtherDerived> operator*(const MatrixBase<OtherDerived>& b) const
{ return MatrixPowerProduct<Derived,MatrixType,OtherDerived>(m_pow, b.derived(), m_p); }
Index rows() const { return m_pow.rows(); }
Index cols() const { return m_pow.cols(); }
private:
Derived& m_pow;
const RealScalar m_p;
MatrixPowerBaseReturnValue& operator=(const MatrixPowerBaseReturnValue&);
};
template<typename Derived, typename MatrixType>
class MatrixPowerBase
{
private:
Derived& derived() { return *static_cast<Derived*>(this); }
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;
explicit MatrixPowerBase(const MatrixType& A) :
m_A(A),
m_conditionNumber(0)
{ eigen_assert(A.rows() == A.cols()); }
#ifndef EIGEN_PARSED_BY_DOXYGEN
const MatrixPowerBaseReturnValue<Derived,MatrixType> operator()(RealScalar p)
{ return MatrixPowerBaseReturnValue<Derived,MatrixType>(derived(), p); }
#endif
void compute(MatrixType& res, RealScalar p)
{ derived().compute(res,p); }
template<typename OtherDerived, typename ResultType>
void compute(const OtherDerived& b, ResultType& res, RealScalar p)
{ derived().compute(b,res,p); }
Index rows() const { return m_A.rows(); }
Index cols() const { return m_A.cols(); }
protected:
typedef Array<RealScalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> RealArray;
typename MatrixType::Nested m_A;
MatrixType m_tmp1, m_tmp2;
RealScalar m_conditionNumber;
};
template<typename Derived, typename Lhs, typename Rhs>
class MatrixPowerProduct : public MatrixBase<MatrixPowerProduct<Derived,Lhs,Rhs> >
{
public:
typedef MatrixBase<MatrixPowerProduct> Base;
EIGEN_DENSE_PUBLIC_INTERFACE(MatrixPowerProduct)
MatrixPowerProduct(Derived& pow, const Rhs& b, RealScalar p) :
m_pow(pow),
m_b(b),
m_p(p)
{ eigen_assert(pow.cols() == b.rows()); }
template<typename ResultType>
inline void evalTo(ResultType& res) const
{ m_pow.compute(m_b, res, m_p); }
inline Index rows() const { return m_pow.rows(); }
inline Index cols() const { return m_b.cols(); }
private:
Derived& m_pow;
typename Rhs::Nested m_b;
const RealScalar m_p;
};
template<typename Derived>
template<typename MatrixPower, typename Lhs, typename Rhs>
Derived& MatrixBase<Derived>::lazyAssign(const MatrixPowerProduct<MatrixPower,Lhs,Rhs>& other)
{
other.evalTo(derived());
return derived();
}
namespace internal {
template<typename Derived, typename MatrixType>
struct traits<MatrixPowerBaseReturnValue<Derived, MatrixType> >
{ typedef MatrixType ReturnType; };
template<typename Derived, typename _Lhs, typename _Rhs>
struct traits<MatrixPowerProduct<Derived,_Lhs,_Rhs> >
{
typedef MatrixXpr XprKind;
typedef typename remove_all<_Lhs>::type Lhs;
typedef typename remove_all<_Rhs>::type Rhs;
typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar;
typedef Dense StorageKind;
typedef typename promote_index_type<typename Lhs::Index, typename Rhs::Index>::type Index;
enum {
RowsAtCompileTime = traits<Lhs>::RowsAtCompileTime,
ColsAtCompileTime = traits<Rhs>::ColsAtCompileTime,
MaxRowsAtCompileTime = traits<Lhs>::MaxRowsAtCompileTime,
MaxColsAtCompileTime = traits<Rhs>::MaxColsAtCompileTime,
Flags = (MaxRowsAtCompileTime==1 ? RowMajorBit : 0)
| EvalBeforeNestingBit | EvalBeforeAssigningBit | NestByRefBit,
CoeffReadCost = 0
};
};
template<int IsComplex>
struct recompose_complex_schur
{
template<typename ResultType, typename MatrixType>
static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U)
{ res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint()); }
};
template<>
struct recompose_complex_schur<0>
{
template<typename ResultType, typename MatrixType>
static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U)
{ res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }
};
template<typename Scalar, int IsComplex = NumTraits<Scalar>::IsComplex>
struct matrix_power_unwinder
{
static inline Scalar run(const Scalar& eival, const Scalar& eival0, int unwindingNumber)
{ return numext::atanh2(eival-eival0, eival+eival0) + Scalar(0, M_PI*unwindingNumber); }
};
template<typename Scalar>
struct matrix_power_unwinder<Scalar,0>
{
static inline Scalar run(Scalar eival, Scalar eival0, int)
{ return numext::atanh2(eival-eival0, eival+eival0); }
};
template<typename T>
inline int binary_powering_cost(T p, int* squarings)
{
int applyings=0, tmp;
frexp(p, squarings);
--*squarings;
while (std::frexp(p, &tmp), tmp > 0) {
p -= std::ldexp(static_cast<T>(0.5), tmp);
++applyings;
}
return applyings;
}
inline int matrix_power_get_pade_degree(float normIminusT)
{
const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f };
int degree = 3;
for (; degree <= 4; ++degree)
if (normIminusT <= maxNormForPade[degree - 3])
break;
return degree;
}
inline int matrix_power_get_pade_degree(double normIminusT)
{
const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1,
1.999045567181744e-1, 2.789358995219730e-1 };
int degree = 3;
for (; degree <= 7; ++degree)
if (normIminusT <= maxNormForPade[degree - 3])
break;
return degree;
}
inline int matrix_power_get_pade_degree(long double normIminusT)
{
#if LDBL_MANT_DIG == 53
const int maxPadeDegree = 7;
const double maxNormForPade[] = { 1.884160592658218e-2L /* degree = 3 */ , 6.038881904059573e-2L, 1.239917516308172e-1L,
1.999045567181744e-1L, 2.789358995219730e-1L };
#elif LDBL_MANT_DIG <= 64
const int maxPadeDegree = 8;
const double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L,
6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L };
#elif LDBL_MANT_DIG <= 106
const int maxPadeDegree = 10;
const double maxNormForPade[] = { 1.0007161601787493236741409687186e-4L /* degree = 3 */ ,
1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L,
2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L,
1.1016843812851143391275867258512e-1L };
#else
const int maxPadeDegree = 10;
const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ ,
6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L,
9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L,
3.908166513900489428442993794761185e-2L, 6.266780814639442865832535460550138e-2L,
9.134603732914548552537150753385375e-2L };
#endif
int degree = 3;
for (; degree <= maxPadeDegree; ++degree)
if (normIminusT <= maxNormForPade[degree - 3])
break;
return degree;
}
} // namespace internal
template<typename MatrixType>
class MatrixPowerTriangularAtomic
{
private:
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
typedef Array<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType;
const MatrixType& m_A;
static void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p);
void compute2x2(MatrixType& res, RealScalar p) const;
void computeBig(MatrixType& res, RealScalar p) const;
public:
explicit MatrixPowerTriangularAtomic(const MatrixType& T);
void compute(MatrixType& res, RealScalar p) const;
};
template<typename MatrixType>
MatrixPowerTriangularAtomic<MatrixType>::MatrixPowerTriangularAtomic(const MatrixType& T) :
m_A(T)
{ eigen_assert(T.rows() == T.cols()); }
template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const
{
res.resizeLike(m_A);
switch (m_A.rows()) {
case 0:
break;
case 1:
res(0,0) = std::pow(m_A(0,0), p);
break;
case 2:
compute2x2(res, p);
break;
default:
computeBig(res, p);
}
}
template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p)
{
int i = degree<<1;
res = (p-degree) / ((i-1)<<1) * IminusT;
for (--i; i; --i) {
res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView<Upper>()
.solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : (p-(i>>1))/((i-1)<<1)) * IminusT).eval();
}
res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
}
// this function assumes that res has the correct size (see bug 614)
template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
{
using std::abs;
using std::pow;
ArrayType logTdiag = m_A.diagonal().array().log();
res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
for (Index i=1; i < m_A.cols(); ++i) {
res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) {
res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
}
else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) {
res.coeffRef(i-1,i) = (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));
}
else {
int unwindingNumber = std::ceil((numext::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI));
Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber);
res.coeffRef(i-1,i) = RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) * std::sinh(p * w)
/ (m_A.coeff(i,i) - m_A.coeff(i-1,i-1));
}
res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
}
}
template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType>::computeBig(MatrixType& res, RealScalar p) const
{
const int digits = std::numeric_limits<RealScalar>::digits;
const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision
digits <= 53? 2.789358995219730e-1: // double precision
digits <= 64? 2.4471944416607995472e-1L: // extended precision
digits <= 106? 1.1016843812851143391275867258512e-1L: // double-double
9.134603732914548552537150753385375e-2L; // quadruple precision
MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
RealScalar normIminusT;
int degree, degree2, numberOfSquareRoots = 0;
bool hasExtraSquareRoot = false;
/* FIXME
* For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite
* loop. We should move 0 eigenvalues to bottom right corner. We need not
* worry about tiny values (e.g. 1e-300) because they will reach 1 if
* repetitively sqrt'ed.
*
* [ T A ]^p [ T^p (T^-1 T^p A) ]
* [ ] = [ ]
* [ 0 0 ] [ 0 0 ]
*/
for (Index i=0; i < m_A.cols(); ++i)
eigen_assert(m_A(i,i) != RealScalar(0));
while (true) {
IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
if (normIminusT < maxNormForPade) {
degree = internal::matrix_power_get_pade_degree(normIminusT);
degree2 = internal::matrix_power_get_pade_degree(normIminusT/2);
if (degree - degree2 <= 1 || hasExtraSquareRoot)
break;
hasExtraSquareRoot = true;
}
MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
T = sqrtT.template triangularView<Upper>();
++numberOfSquareRoots;
}
computePade(degree, IminusT, res, p);
for (; numberOfSquareRoots; --numberOfSquareRoots) {
compute2x2(res, std::ldexp(p,-numberOfSquareRoots));
res = res.template triangularView<Upper>() * res;
}
compute2x2(res, p);
}
} // namespace Eigen
#endif // EIGEN_MATRIX_POWER

71
unsupported/test/matrix_power.cpp
View File

@ -109,61 +109,8 @@ void testExponentLaws(const MatrixType& m, double tol)
} }
} }
template<typename MatrixType, typename VectorType>
void testProduct(const MatrixType& m, const VectorType& v, double tol)
{
typedef typename MatrixType::RealScalar RealScalar;
MatrixType m1;
VectorType v1, v2, v3;
RealScalar p;
for (int i=0; i < g_repeat; ++i) {
generateTestMatrix<MatrixType>::run(m1, m.rows());
MatrixPower<MatrixType> mpow(m1);
v1 = VectorType::Random(v.rows(), v.cols());
p = internal::random<RealScalar>();
v2.noalias() = mpow(p) * v1;
v3.noalias() = mpow(p).eval() * v1;
std::cout << "testProduct: error powerm = " << relerr(v2, v3) << '\n';
VERIFY(v2.isApprox(v3, static_cast<RealScalar>(tol)));
}
}
template<typename MatrixType, typename VectorType>
void testTriangularProduct(const MatrixType& m, const VectorType& v, double tol)
{
typedef typename MatrixType::RealScalar RealScalar;
MatrixType m1;
VectorType v1, v2, v3;
RealScalar p;
for (int i=0; i < g_repeat; ++i) {
generateTriangularMatrix<MatrixType>::run(m1, m.rows());
MatrixPowerTriangular<MatrixType> mpow(m1);
v1 = VectorType::Random(v.rows(), v.cols());
p = internal::random<RealScalar>();
v2.noalias() = mpow(p) * v1;
v3.noalias() = mpow(p).eval() * v1;
std::cout << "testTriangularProduct: error powerm = " << relerr(v2, v3) << '\n';
VERIFY(v2.isApprox(v3, static_cast<RealScalar>(tol)));
}
}
template<typename MatrixType, typename VectorType>
void testMatrixVector(const MatrixType& m, const VectorType& v, double tol)
{
testExponentLaws(m,tol);
testProduct(m,v,tol);
testTriangularProduct(m,v,tol);
}
typedef Matrix<double,3,3,RowMajor> Matrix3dRowMajor; typedef Matrix<double,3,3,RowMajor> Matrix3dRowMajor;
typedef Matrix<long double,Dynamic,Dynamic> MatrixXe; typedef Matrix<long double,Dynamic,Dynamic> MatrixXe;
typedef Matrix<long double,Dynamic,1> VectorXe;
void test_matrix_power() void test_matrix_power()
{ {
@ -174,13 +121,13 @@ void test_matrix_power()
CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5)); CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14)); CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14));
CALL_SUBTEST_2(testMatrixVector(Matrix2d(), Vector2d(), 1e-13)); CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13));
CALL_SUBTEST_7(testMatrixVector(Matrix3dRowMajor(), MatrixXd(3,5), 1e-13)); CALL_SUBTEST_7(testExponentLaws(Matrix3dRowMajor(), 1e-13));
CALL_SUBTEST_3(testMatrixVector(Matrix4cd(), Vector4cd(), 1e-13)); CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13));
CALL_SUBTEST_4(testMatrixVector(MatrixXd(8,8), VectorXd(8), 2e-12)); CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 2e-12));
CALL_SUBTEST_1(testMatrixVector(Matrix2f(), Vector2f(), 1e-4)); CALL_SUBTEST_1(testExponentLaws(Matrix2f(), 1e-4));
CALL_SUBTEST_5(testMatrixVector(Matrix3cf(), Vector3cf(), 1e-4)); CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4));
CALL_SUBTEST_8(testMatrixVector(Matrix4f(), Vector4f(), 1e-4)); CALL_SUBTEST_8(testExponentLaws(Matrix4f(), 1e-4));
CALL_SUBTEST_6(testMatrixVector(MatrixXf(2,2), VectorXf(2), 1e-3)); // see bug 614 CALL_SUBTEST_6(testExponentLaws(MatrixXf(2,2), 1e-3)); // see bug 614
CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7), VectorXe(7), 1e-13)); CALL_SUBTEST_9(testExponentLaws(MatrixXe(7,7), 1e-13));
} }