GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-14 12:46:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

Tidy up and write dox.

Browse Source
This commit is contained in:
Chen-Pang He 2012-08-28 01:55:13 +08:00
parent 5252d823c9
commit ba4e886376
5 changed files with 178 additions and 291 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Eigen/src/Core/MatrixBase.h
View File

@ -454,8 +454,7 @@ template<typename Derived> class MatrixBase
const MatrixFunctionReturnValue<Derived> sin() const; const MatrixFunctionReturnValue<Derived> sin() const;
const MatrixSquareRootReturnValue<Derived> sqrt() const; const MatrixSquareRootReturnValue<Derived> sqrt() const;
const MatrixLogarithmReturnValue<Derived> log() const; const MatrixLogarithmReturnValue<Derived> log() const;
template <typename ExponentType> const MatrixPowerReturnValue<Derived> pow(RealScalar p) const;
const MatrixPowerReturnValue<Derived, ExponentType> pow(const ExponentType& p) const;
#ifdef EIGEN2_SUPPORT #ifdef EIGEN2_SUPPORT
template<typename ProductDerived, typename Lhs, typename Rhs> template<typename ProductDerived, typename Lhs, typename Rhs>

2
Eigen/src/Core/util/ForwardDeclarations.h
View File

@ -271,7 +271,7 @@ template<typename Derived> struct MatrixExponentialReturnValue;
template<typename Derived> class MatrixFunctionReturnValue; template<typename Derived> class MatrixFunctionReturnValue;
template<typename Derived> class MatrixSquareRootReturnValue; template<typename Derived> class MatrixSquareRootReturnValue;
template<typename Derived> class MatrixLogarithmReturnValue; template<typename Derived> class MatrixLogarithmReturnValue;
template<typename Derived, typename ExponentType> class MatrixPowerReturnValue; template<typename Derived> class MatrixPowerReturnValue;
namespace internal { namespace internal {
template <typename Scalar> template <typename Scalar>

12
unsupported/Eigen/MatrixFunctions
View File

@ -223,8 +223,7 @@ Output: \verbinclude MatrixLogarithm.out
Compute the matrix raised to arbitrary real power. Compute the matrix raised to arbitrary real power.
\code \code
template <typename ExponentType> const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const
const MatrixPowerReturnValue<Derived, ExponentType> MatrixBase<Derived>::pow(const ExponentType& p) const
\endcode \endcode
\param[in] M base of the matrix power, should be a square matrix. \param[in] M base of the matrix power, should be a square matrix.
@ -247,6 +246,15 @@ diagonal and the first super-diagonal is directly computed.
The actual work is done by the MatrixPower class, which can compute The actual work is done by the MatrixPower class, which can compute
\f$ M^p v \f$, where \p v is another matrix with the same rows as \f$ M^p v \f$, where \p v is another matrix with the same rows as
\p M. The matrix \p v is set to be the identity matrix by default. \p M. The matrix \p v is set to be the identity matrix by default.
Therefore, the expression <tt>M.pow(p) * v</tt> is specialized for
this. No temporary storage is created for the result. The code below
directly evaluates R-values into L-values without aliasing issue. Do
\b NOT try to \a optimize with noalias(). It won't compile.
\code
v = m.pow(p) * v;
m = m.pow(q);
// v2.noalias() = m.pow(p) * v1; Won't compile!
\endcode
Details of the algorithm can be found in: Nicholas J. Higham and Details of the algorithm can be found in: Nicholas J. Higham and
Lijing Lin, "A Schur-Pad&eacute; algorithm for fractional powers of a Lijing Lin, "A Schur-Pad&eacute; algorithm for fractional powers of a

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

@ -23,10 +23,9 @@ namespace Eigen {
* *
* \tparam MatrixType type of the base, expected to be an instantiation * \tparam MatrixType type of the base, expected to be an instantiation
* of the Matrix class template. * of the Matrix class template.
* \tparam IsInteger used internally to select correct specialization.
* \tparam PlainObject type of the multiplier. * \tparam PlainObject type of the multiplier.
*/ */
template <typename MatrixType, int IsInteger, typename PlainObject = MatrixType> template<typename MatrixType, typename PlainObject = MatrixType>
class MatrixPower class MatrixPower
{ {
private: private:
@ -65,11 +64,11 @@ class MatrixPower
* *
* \param[out] result \f$ A^p b \f$, as specified in the constructor. * \param[out] result \f$ A^p b \f$, as specified in the constructor.
*/ */
template <typename ResultType> void compute(ResultType& result); template<typename ResultType> void compute(ResultType& result);
private: private:
/** /**
* \brief Compute the matrix power. * \brief Compute the matrix to integral power.
* *
* If \p b is \em fatter than \p A, it computes \f$ A^{p_{\textrm int}} * If \p b is \em fatter than \p A, it computes \f$ A^{p_{\textrm int}}
* \f$ first, and then multiplies it with \p b. Otherwise, * \f$ first, and then multiplies it with \p b. Otherwise,
@ -77,7 +76,7 @@ class MatrixPower
* *
* \sa computeChainProduct(ResultType&); * \sa computeChainProduct(ResultType&);
*/ */
template <typename ResultType> template<typename ResultType>
void computeIntPower(ResultType& result); void computeIntPower(ResultType& result);
/** /**
@ -87,14 +86,14 @@ class MatrixPower
* powering or matrix chain multiplication or solving the linear system * powering or matrix chain multiplication or solving the linear system
* repetitively, according to which algorithm costs less. * repetitively, according to which algorithm costs less.
*/ */
template <typename ResultType> template<typename ResultType>
void computeChainProduct(ResultType&); void computeChainProduct(ResultType&);
/** \brief Compute the cost of binary powering. */ /** \brief Compute the cost of binary powering. */
static int computeCost(RealScalar); static int computeCost(RealScalar);
/** \brief Solve the linear system repetitively. */ /** \brief Solve the linear system repetitively. */
template <typename ResultType> template<typename ResultType>
void partialPivLuSolve(ResultType&, RealScalar); void partialPivLuSolve(ResultType&, RealScalar);
/** \brief Compute Schur decomposition of #m_A. */ /** \brief Compute Schur decomposition of #m_A. */
@ -170,76 +169,9 @@ class MatrixPower
ComplexArray m_logTdiag; ///< \brief Logarithm of the main diagonal of #m_T. ComplexArray m_logTdiag; ///< \brief Logarithm of the main diagonal of #m_T.
}; };
/** template<typename MatrixType, typename PlainObject>
* \internal \ingroup MatrixFunctions_Module template<typename ResultType>
* \brief Partial specialization for integral exponents. void MatrixPower<MatrixType,PlainObject>::compute(ResultType& result)
*/
template <typename MatrixType, typename PlainObject>
class MatrixPower<MatrixType, 1, PlainObject>
{
public:
/**
* \brief Constructor.
*
* \param[in] A the base of the matrix power.
* \param[in] p the exponent of the matrix power.
* \param[in] b the multiplier.
*/
MatrixPower(const MatrixType& A, int p, const PlainObject& b) :
m_A(A),
m_p(p),
m_b(b),
m_dimA(A.cols()),
m_dimb(b.cols())
{ /* empty body */ }
/**
* \brief Compute the matrix power.
*
* If \p b is \em fatter than \p A, it computes \f$ A^p \f$ first, and
* then multiplies it with \p b. Otherwise, #computeChainProduct
* optimizes the expression.
*
* \param[out] result \f$ A^p b \f$, as specified in the constructor.
*
* \sa computeChainProduct(ResultType&);
*/
template <typename ResultType>
void compute(ResultType& result);
private:
typedef typename MatrixType::Index Index;
const MatrixType& m_A; ///< \brief Reference to the matrix base.
const int m_p; ///< \brief The integral exponent.
const PlainObject& m_b; ///< \brief Reference to the multiplier.
const Index m_dimA; ///< \brief The dimension of #m_A, equivalent to %m_A.cols().
const Index m_dimb; ///< \brief The dimension of #m_b, equivalent to %m_b.cols().
MatrixType m_tmp; ///< \brief Used for temporary storage.
/**
* \brief Convert matrix power into chain product.
*
* This optimizes the matrix expression. It automatically chooses binary
* powering or matrix chain multiplication or solving the linear system
* repetitively, according to which algorithm costs less.
*/
template <typename ResultType>
void computeChainProduct(ResultType& result);
/** \brief Compute the cost of binary powering. */
static int computeCost(int);
/** \brief Solve the linear system repetitively. */
template <typename ResultType>
void partialPivLuSolve(ResultType&, int);
};
/******* Specialized for real exponents *******/
template <typename MatrixType, int IsInteger, typename PlainObject>
template <typename ResultType>
void MatrixPower<MatrixType,IsInteger,PlainObject>::compute(ResultType& result)
{ {
using std::floor; using std::floor;
using std::pow; using std::pow;
@ -255,19 +187,18 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::compute(ResultType& result)
computeSchurDecomposition(); computeSchurDecomposition();
getFractionalExponent(); getFractionalExponent();
computeIntPower(result); computeIntPower(result);
if (m_dimA == 2) { if (m_dimA == 2)
compute2x2(m_pFrac); compute2x2(m_pFrac);
} else { else
computeBig(); computeBig();
}
computeTmp(Scalar()); computeTmp(Scalar());
result = m_tmp * result; result = m_tmp * result;
} }
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
template <typename ResultType> template<typename ResultType>
void MatrixPower<MatrixType,IsInteger,PlainObject>::computeIntPower(ResultType& result) void MatrixPower<MatrixType,PlainObject>::computeIntPower(ResultType& result)
{ {
MatrixType tmp; MatrixType tmp;
if (m_dimb > m_dimA) { if (m_dimb > m_dimA) {
@ -280,9 +211,9 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::computeIntPower(ResultType&
} }
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
template <typename ResultType> template<typename ResultType>
void MatrixPower<MatrixType,IsInteger,PlainObject>::computeChainProduct(ResultType& result) void MatrixPower<MatrixType,PlainObject>::computeChainProduct(ResultType& result)
{ {
using std::abs; using std::abs;
using std::fmod; using std::fmod;
@ -314,8 +245,8 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::computeChainProduct(ResultTy
result = m_tmp * result; result = m_tmp * result;
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
int MatrixPower<MatrixType,IsInteger,PlainObject>::computeCost(RealScalar p) int MatrixPower<MatrixType,PlainObject>::computeCost(RealScalar p)
{ {
using std::frexp; using std::frexp;
using std::ldexp; using std::ldexp;
@ -329,25 +260,25 @@ int MatrixPower<MatrixType,IsInteger,PlainObject>::computeCost(RealScalar p)
return cost; return cost;
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
template <typename ResultType> template<typename ResultType>
void MatrixPower<MatrixType,IsInteger,PlainObject>::partialPivLuSolve(ResultType& result, RealScalar p) void MatrixPower<MatrixType,PlainObject>::partialPivLuSolve(ResultType& result, RealScalar p)
{ {
const PartialPivLU<MatrixType> Asolver(m_A); const PartialPivLU<MatrixType> Asolver(m_A);
for (; p >= RealScalar(1); p--) for (; p >= RealScalar(1); p--)
result = Asolver.solve(result); result = Asolver.solve(result);
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
void MatrixPower<MatrixType,IsInteger,PlainObject>::computeSchurDecomposition() void MatrixPower<MatrixType,PlainObject>::computeSchurDecomposition()
{ {
const ComplexSchur<MatrixType> schurOfA(m_A); const ComplexSchur<MatrixType> schurOfA(m_A);
m_T = schurOfA.matrixT(); m_T = schurOfA.matrixT();
m_U = schurOfA.matrixU(); m_U = schurOfA.matrixU();
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
void MatrixPower<MatrixType,IsInteger,PlainObject>::getFractionalExponent() void MatrixPower<MatrixType,PlainObject>::getFractionalExponent()
{ {
using std::pow; using std::pow;
typedef Array<RealScalar, Rows, 1, ColMajor, MaxRows> RealArray; typedef Array<RealScalar, Rows, 1, ColMajor, MaxRows> RealArray;
@ -365,9 +296,9 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::getFractionalExponent()
} }
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
std::complex<typename MatrixType::RealScalar> std::complex<typename MatrixType::RealScalar>
MatrixPower<MatrixType,IsInteger,PlainObject>::atanh2(const ComplexScalar& y, const ComplexScalar& x) MatrixPower<MatrixType,PlainObject>::atanh2(const ComplexScalar& y, const ComplexScalar& x)
{ {
using std::abs; using std::abs;
using std::log; using std::log;
@ -380,8 +311,8 @@ MatrixPower<MatrixType,IsInteger,PlainObject>::atanh2(const ComplexScalar& y, co
return z + z*z*z / RealScalar(3); return z + z*z*z / RealScalar(3);
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
void MatrixPower<MatrixType,IsInteger,PlainObject>::compute2x2(RealScalar p) void MatrixPower<MatrixType,PlainObject>::compute2x2(RealScalar p)
{ {
using std::abs; using std::abs;
using std::ceil; using std::ceil;
@ -413,8 +344,8 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::compute2x2(RealScalar p)
} }
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
void MatrixPower<MatrixType,IsInteger,PlainObject>::computeBig() void MatrixPower<MatrixType,PlainObject>::computeBig()
{ {
using std::ldexp; using std::ldexp;
const int digits = std::numeric_limits<RealScalar>::digits; const int digits = std::numeric_limits<RealScalar>::digits;
@ -450,8 +381,8 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::computeBig()
compute2x2(m_pFrac); compute2x2(m_pFrac);
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(float normIminusT) inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(float normIminusT)
{ {
const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f }; const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f };
int degree = 3; int degree = 3;
@ -461,8 +392,8 @@ inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(float no
return degree; return degree;
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(double normIminusT) inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(double normIminusT)
{ {
const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2,
1.239917516308172e-1, 1.999045567181744e-1, 2.789358995219730e-1 }; 1.239917516308172e-1, 1.999045567181744e-1, 2.789358995219730e-1 };
@ -473,8 +404,8 @@ inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(double n
return degree; return degree;
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(long double normIminusT) inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(long double normIminusT)
{ {
#if LDBL_MANT_DIG == 53 #if LDBL_MANT_DIG == 53
const int maxPadeDegree = 7; const int maxPadeDegree = 7;
@ -506,8 +437,8 @@ inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(long dou
break; break;
return degree; return degree;
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
void MatrixPower<MatrixType,IsInteger,PlainObject>::computePade(const int& degree, const ComplexMatrix& IminusT) void MatrixPower<MatrixType,PlainObject>::computePade(const int& degree, const ComplexMatrix& IminusT)
{ {
int i = degree << 1; int i = degree << 1;
m_fT = coeff(i) * IminusT; m_fT = coeff(i) * IminusT;
@ -518,8 +449,8 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::computePade(const int& degre
m_fT += ComplexMatrix::Identity(m_dimA, m_dimA); m_fT += ComplexMatrix::Identity(m_dimA, m_dimA);
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
inline typename MatrixType::RealScalar MatrixPower<MatrixType,IsInteger,PlainObject>::coeff(const int& i) inline typename MatrixType::RealScalar MatrixPower<MatrixType,PlainObject>::coeff(const int& i)
{ {
if (i == 1) if (i == 1)
return -m_pFrac; return -m_pFrac;
@ -529,82 +460,72 @@ inline typename MatrixType::RealScalar MatrixPower<MatrixType,IsInteger,PlainObj
return (m_pFrac - RealScalar(i >> 1)) / RealScalar((i - 1) << 1); return (m_pFrac - RealScalar(i >> 1)) / RealScalar((i - 1) << 1);
} }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
void MatrixPower<MatrixType,IsInteger,PlainObject>::computeTmp(RealScalar) void MatrixPower<MatrixType,PlainObject>::computeTmp(RealScalar)
{ m_tmp = (m_U * m_fT * m_U.adjoint()).real(); } { m_tmp = (m_U * m_fT * m_U.adjoint()).real(); }
template <typename MatrixType, int IsInteger, typename PlainObject> template<typename MatrixType, typename PlainObject>
void MatrixPower<MatrixType,IsInteger,PlainObject>::computeTmp(ComplexScalar) void MatrixPower<MatrixType,PlainObject>::computeTmp(ComplexScalar)
{ m_tmp = m_U * m_fT * m_U.adjoint(); } { m_tmp = m_U * m_fT * m_U.adjoint(); }
/******* Specialized for integral exponents *******/ /**
* \ingroup MatrixFunctions_Module
template <typename MatrixType, typename PlainObject> *
template <typename ResultType> * \brief Proxy for the matrix power multiplied by other matrix.
void MatrixPower<MatrixType,1,PlainObject>::compute(ResultType& result) *
* \tparam Derived type of the base, a matrix (expression).
* \tparam RhsDerived type of the multiplier.
*
* This class holds the arguments to the matrix power until it is
* assigned or evaluated for some other reason (so the argument
* should not be changed in the meantime). It is the return type of
* MatrixPowerReturnValue::operator*() and related functions and most
* of the time this is the only way it is used.
*/
template<typename Derived, typename RhsDerived>
class MatrixPowerProductReturnValue : public ReturnByValue<MatrixPowerProductReturnValue<Derived,RhsDerived> >
{ {
MatrixType tmp; private:
if (m_dimb > m_dimA) { typedef typename Derived::PlainObject MatrixType;
tmp = MatrixType::Identity(m_dimA, m_dimA); typedef typename RhsDerived::PlainObject PlainObject;
computeChainProduct(tmp); typedef typename RhsDerived::RealScalar RealScalar;
result.noalias() = tmp * m_b; typedef typename RhsDerived::Index Index;
} else {
result = m_b;
computeChainProduct(result);
}
}
template <typename MatrixType, typename PlainObject> public:
int MatrixPower<MatrixType,1,PlainObject>::computeCost(int p) /**
{ * \brief Constructor.
int cost = 0, tmp; *
for (tmp = p; tmp; tmp >>= 1) * \param[in] A %Matrix (expression), the base of the matrix power.
cost++; * \param[in] p scalar, the exponent of the matrix power.
for (tmp = 1; tmp <= p; tmp <<= 1) * \prarm[in] b %Matrix (expression), the multiplier.
if (tmp & p) cost++; */
return cost; MatrixPowerProductReturnValue(const Derived& A, RealScalar p, const RhsDerived& b)
} : m_A(A), m_p(p), m_b(b) { }
template <typename MatrixType, typename PlainObject> /**
template <typename ResultType> * \brief Compute the expression.
void MatrixPower<MatrixType,1,PlainObject>::partialPivLuSolve(ResultType& result, int p) *
{ * \param[out] result \f$ A^p b \f$ where \p A, \p p and \p bare as
const PartialPivLU<MatrixType> Asolver(m_A); * in the constructor.
for(; p; p--) */
result = Asolver.solve(result); template<typename ResultType>
} inline void evalTo(ResultType& result) const
{
template <typename MatrixType, typename PlainObject> const MatrixType A = m_A;
template <typename ResultType> const PlainObject b = m_b;
void MatrixPower<MatrixType,1,PlainObject>::computeChainProduct(ResultType& result) MatrixPower<MatrixType, PlainObject> mp(A, m_p, b);
{ mp.compute(result);
using std::abs;
int p = abs(m_p), cost = computeCost(p);
if (m_p < 0) {
if (p * m_dimb <= cost * m_dimA) {
partialPivLuSolve(result, p);
return;
} else {
m_tmp = m_A.inverse();
}
} else {
m_tmp = m_A;
} }
while (p * m_dimb > cost * m_dimA) { Index rows() const { return m_b.rows(); }
if (p & 1) { Index cols() const { return m_b.cols(); }
result = m_tmp * result;
cost--;
}
m_tmp *= m_tmp;
cost--;
p >>= 1;
}
for (; p; p--) private:
result = m_tmp * result; const Derived& m_A;
} const RealScalar m_p;
const RhsDerived& m_b;
MatrixPowerProductReturnValue& operator=(const MatrixPowerProductReturnValue&);
};
/** /**
* \ingroup MatrixFunctions_Module * \ingroup MatrixFunctions_Module
@ -612,7 +533,6 @@ void MatrixPower<MatrixType,1,PlainObject>::computeChainProduct(ResultType& resu
* \brief Proxy for the matrix power of some matrix (expression). * \brief Proxy for the matrix power of some matrix (expression).
* *
* \tparam Derived type of the base, a matrix (expression). * \tparam Derived type of the base, a matrix (expression).
* \tparam ExponentType type of the exponent, a scalar.
* *
* This class holds the arguments to the matrix power until it is * This class holds the arguments to the matrix power until it is
* assigned or evaluated for some other reason (so the argument * assigned or evaluated for some other reason (so the argument
@ -620,19 +540,21 @@ void MatrixPower<MatrixType,1,PlainObject>::computeChainProduct(ResultType& resu
* MatrixBase::pow() and related functions and most of the * MatrixBase::pow() and related functions and most of the
* time this is the only way it is used. * time this is the only way it is used.
*/ */
template<typename Derived, typename ExponentType> class MatrixPowerReturnValue template<typename Derived>
: public ReturnByValue<MatrixPowerReturnValue<Derived, ExponentType> > class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Derived> >
{ {
public: private:
typedef typename Derived::RealScalar RealScalar;
typedef typename Derived::Index Index; typedef typename Derived::Index Index;
public:
/** /**
* \brief Constructor. * \brief Constructor.
* *
* \param[in] A %Matrix (expression), the base of the matrix power. * \param[in] A %Matrix (expression), the base of the matrix power.
* \param[in] p scalar, the exponent of the matrix power. * \param[in] p scalar, the exponent of the matrix power.
*/ */
MatrixPowerReturnValue(const Derived& A, const ExponentType& p) MatrixPowerReturnValue(const Derived& A, RealScalar p)
: m_A(A), m_p(p) { } : m_A(A), m_p(p) { }
/** /**
@ -641,22 +563,19 @@ template<typename Derived, typename ExponentType> class MatrixPowerReturnValue
* The %MatrixPower class can optimize \f$ A^p b \f$ computing, and * The %MatrixPower class can optimize \f$ A^p b \f$ computing, and
* this method provides an elegant way to call it. * this method provides an elegant way to call it.
* *
* Unlike general matrix-matrix / matrix-vector product, this does
* \b NOT produce a temporary storage for the result. Therefore,
* the code below is \a already optimal:
* \code
* v = A.pow(p) * b;
* // v.noalias() = A.pow(p) * b; Won't compile!
* \endcode
*
* \param[in] b %Matrix (expression), the multiplier. * \param[in] b %Matrix (expression), the multiplier.
* \param[out] result \f$ A^p b \f$ where \p A and \p p are as in
* the constructor.
*/ */
template <typename OtherDerived> template<typename RhsDerived>
const typename OtherDerived::PlainObject operator*(const MatrixBase<OtherDerived>& b) const const MatrixPowerProductReturnValue<Derived,RhsDerived> operator*(const MatrixBase<RhsDerived>& b) const
{ { return MatrixPowerProductReturnValue<Derived,RhsDerived>(m_A, m_p, b.derived()); }
typedef typename OtherDerived::PlainObject PlainObject;
const int IsInteger = NumTraits<ExponentType>::IsInteger;
const typename Derived::PlainObject Aevaluated = m_A.eval();
const PlainObject bevaluated = b.eval();
PlainObject result;
MatrixPower<Derived, IsInteger, PlainObject> mp(Aevaluated, m_p, bevaluated);
mp.compute(result);
return result;
}
/** /**
* \brief Compute the matrix power. * \brief Compute the matrix power.
@ -664,14 +583,13 @@ template<typename Derived, typename ExponentType> class MatrixPowerReturnValue
* \param[out] result \f$ A^p \f$ where \p A and \p p are as in the * \param[out] result \f$ A^p \f$ where \p A and \p p are as in the
* constructor. * constructor.
*/ */
template <typename ResultType> template<typename ResultType>
inline void evalTo(ResultType& result) const inline void evalTo(ResultType& result) const
{ {
typedef typename Derived::PlainObject PlainObject; typedef typename Derived::PlainObject PlainObject;
const int IsInteger = NumTraits<ExponentType>::IsInteger; const PlainObject A = m_A;
const PlainObject Aevaluated = m_A.eval();
const PlainObject Identity = PlainObject::Identity(m_A.rows(), m_A.cols()); const PlainObject Identity = PlainObject::Identity(m_A.rows(), m_A.cols());
MatrixPower<PlainObject, IsInteger> mp(Aevaluated, m_p, Identity); MatrixPower<PlainObject> mp(A, m_p, Identity);
mp.compute(result); mp.compute(result);
} }
@ -680,25 +598,29 @@ template<typename Derived, typename ExponentType> class MatrixPowerReturnValue
private: private:
const Derived& m_A; const Derived& m_A;
const ExponentType& m_p; const RealScalar m_p;
MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&); MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&);
}; };
namespace internal { namespace internal {
template<typename Derived, typename ExponentType> template<typename Derived>
struct traits<MatrixPowerReturnValue<Derived, ExponentType> > struct traits<MatrixPowerReturnValue<Derived> >
{ {
typedef typename Derived::PlainObject ReturnType; typedef typename Derived::PlainObject ReturnType;
}; };
template<typename Derived, typename RhsDerived>
struct traits<MatrixPowerProductReturnValue<Derived,RhsDerived> >
{
typedef typename RhsDerived::PlainObject ReturnType;
};
} }
template <typename Derived> template<typename Derived>
template <typename ExponentType> const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const
const MatrixPowerReturnValue<Derived, ExponentType> MatrixBase<Derived>::pow(const ExponentType& p) const
{ {
eigen_assert(rows() == cols()); eigen_assert(rows() == cols());
return MatrixPowerReturnValue<Derived, ExponentType>(derived(), p); return MatrixPowerReturnValue<Derived>(derived(), p);
} }
} // end namespace Eigen } // end namespace Eigen

108
unsupported/test/matrix_power.cpp
View File

@ -9,7 +9,7 @@
#include "matrix_functions.h" #include "matrix_functions.h"
template <typename T> template<typename T>
void test2dRotation(double tol) void test2dRotation(double tol)
{ {
Matrix<T,2,2> A, B, C; Matrix<T,2,2> A, B, C;
@ -28,7 +28,7 @@ void test2dRotation(double tol)
} }
} }
template <typename T> template<typename T>
void test2dHyperbolicRotation(double tol) void test2dHyperbolicRotation(double tol)
{ {
Matrix<std::complex<T>,2,2> A, B, C; Matrix<std::complex<T>,2,2> A, B, C;
@ -48,55 +48,7 @@ void test2dHyperbolicRotation(double tol)
} }
} }
template <typename MatrixType> template<typename MatrixType>
void testIntPowers(const MatrixType& m, double tol)
{
typedef typename MatrixType::RealScalar RealScalar;
const MatrixType m1 = MatrixType::Random(m.rows(), m.cols());
const MatrixType identity = MatrixType::Identity(m.rows(), m.cols());
const PartialPivLU<MatrixType> solver(m1);
MatrixType m2, m3, m4;
m3 = m1.pow(0);
m4 = m1.pow(0.);
std::cout << "testIntPower: i = 0 error powerm = " << relerr(identity, m3) << " " << relerr(identity, m4) << '\n';
VERIFY(identity == m3 && identity == m4);
m3 = m1.pow(1);
m4 = m1.pow(1.);
std::cout << "testIntPower: i = 1 error powerm = " << relerr(m1, m3) << " " << relerr(m1, m4) << '\n';
VERIFY(m1 == m3 && m1 == m4);
m2.noalias() = m1 * m1;
m3 = m1.pow(2);
m4 = m1.pow(2.);
std::cout << "testIntPower: i = 2 error powerm = " << relerr(m2, m3) << " " << relerr(m2, m4) << '\n';
VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol)));
for (int i = 3; i <= 20; i++) {
m2 *= m1;
m3 = m1.pow(i);
m4 = m1.pow(RealScalar(i));
std::cout << "testIntPower: i = " << i << " error powerm = " << relerr(m2, m3) << " " << relerr (m2, m4) << '\n';
VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol)));
}
m2 = solver.inverse();
m3 = m1.pow(-1);
m4 = m1.pow(-1.);
std::cout << "testIntPower: i = -1 error powerm = " << relerr(m2, m3) << " " << relerr (m2, m4) << '\n';
VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol)));
for (int i = -2; i >= -20; i--) {
m2 = solver.solve(m2);
m3 = m1.pow(i);
m4 = m1.pow(RealScalar(i));
std::cout << "testIntPower: i = " << i << " error powerm = " << relerr(m2, m3) << " " << relerr (m2, m4) << '\n';
VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol)));
}
}
template <typename MatrixType>
void testExponentLaws(const MatrixType& m, double tol) void testExponentLaws(const MatrixType& m, double tol)
{ {
typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::RealScalar RealScalar;
@ -129,31 +81,40 @@ void testExponentLaws(const MatrixType& m, double tol)
} }
} }
template <typename MatrixType, typename VectorType> template<typename MatrixType, typename VectorType>
void testMatrixVectorProduct(const MatrixType& m, const VectorType& v, double tol) void testMatrixVectorProduct(const MatrixType& m, const VectorType& v, double tol)
{ {
typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::RealScalar RealScalar;
MatrixType m1; MatrixType m1;
VectorType v1, v2, v3; VectorType v1, v2, v3;
RealScalar pReal; RealScalar p;
signed char pInt;
for (int i = 0; i < g_repeat; i++) { for (int i = 0; i < g_repeat; i++) {
generateTestMatrix<MatrixType>::run(m1, m.rows()); generateTestMatrix<MatrixType>::run(m1, m.rows());
v1 = VectorType::Random(v.rows(), v.cols()); v1 = VectorType::Random(v.rows(), v.cols());
pReal = internal::random<RealScalar>(); p = internal::random<RealScalar>();
pInt = rand();
pInt >>= 2;
v2.noalias() = m1.pow(pReal).eval() * v1; v2.noalias() = m1.pow(p).eval() * v1;
v3.noalias() = m1.pow(pReal) * v1; v1 = m1.pow(p) * v1;
std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v3); std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v1) << '\n';
VERIFY(v2.isApprox(v3, RealScalar(tol))); VERIFY(v2.isApprox(v1, RealScalar(tol)));
}
}
v2.noalias() = m1.pow(pInt).eval() * v1; template<typename MatrixType>
v3.noalias() = m1.pow(pInt) * v1; void testAliasing(const MatrixType& m)
std::cout << " " << relerr(v2, v3) << '\n'; {
VERIFY(v2.isApprox(v3, RealScalar(tol)) || v2 == v3); typedef typename MatrixType::RealScalar RealScalar;
MatrixType m1, m2;
RealScalar p;
for (int i = 0; i < g_repeat; i++) {
generateTestMatrix<MatrixType>::run(m1, m.rows());
p = internal::random<RealScalar>();
m2 = m1.pow(p);
m1 = m1.pow(p);
VERIFY(m1 == m2);
} }
} }
@ -166,15 +127,6 @@ 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(testIntPowers(Matrix2d(), 1e-13));
CALL_SUBTEST_7(testIntPowers(Matrix<double,3,3,RowMajor>(), 1e-13));
CALL_SUBTEST_3(testIntPowers(Matrix4cd(), 1e-13));
CALL_SUBTEST_4(testIntPowers(MatrixXd(8,8), 1e-13));
CALL_SUBTEST_1(testIntPowers(Matrix2f(), 1e-4));
CALL_SUBTEST_5(testIntPowers(Matrix3cf(), 1e-4));
CALL_SUBTEST_8(testIntPowers(Matrix4f(), 1e-4));
CALL_SUBTEST_6(testIntPowers(MatrixXf(8,8), 1e-4));
CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13)); CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13));
CALL_SUBTEST_7(testExponentLaws(Matrix<double,3,3,RowMajor>(), 1e-13)); CALL_SUBTEST_7(testExponentLaws(Matrix<double,3,3,RowMajor>(), 1e-13));
CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13)); CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13));
@ -192,5 +144,11 @@ void test_matrix_power()
CALL_SUBTEST_5(testMatrixVectorProduct(Matrix3cf(), Vector3cf(), 1e-4)); CALL_SUBTEST_5(testMatrixVectorProduct(Matrix3cf(), Vector3cf(), 1e-4));
CALL_SUBTEST_8(testMatrixVectorProduct(Matrix4f(), Vector4f(), 1e-4)); CALL_SUBTEST_8(testMatrixVectorProduct(Matrix4f(), Vector4f(), 1e-4));
CALL_SUBTEST_6(testMatrixVectorProduct(MatrixXf(8,8), VectorXf(8), 1e-4)); CALL_SUBTEST_6(testMatrixVectorProduct(MatrixXf(8,8), VectorXf(8), 1e-4));
CALL_SUBTEST_10(testMatrixVectorProduct(Matrix<long double,Dynamic,Dynamic>(7,7), Matrix<long double,7,9>(), 1e-13)); CALL_SUBTEST_9(testMatrixVectorProduct(Matrix<long double,Dynamic,Dynamic>(7,7), Matrix<long double,7,9>(), 1e-13));
CALL_SUBTEST_7(testAliasing(Matrix<double,3,3,RowMajor>()));
CALL_SUBTEST_3(testAliasing(Matrix4cd()));
CALL_SUBTEST_4(testAliasing(MatrixXd(8,8)));
CALL_SUBTEST_5(testAliasing(Matrix3cf()));
CALL_SUBTEST_6(testAliasing(MatrixXf(8,8)));
} }