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

Fix bug in MatrixPower(expression) due to destruction of temporary objects. Sorry for ugly pointer manipulation but it prevents copying a PlainObject.

Browse Source
This commit is contained in:
Chen-Pang He 2012-09-23 18:49:44 +08:00
parent 963794b04a
commit 1d402dac03
1 changed files with 56 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -48,10 +48,10 @@ template<typename MatrixType> class MatrixPower
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
typedef Matrix<std::complex<RealScalar>,Rows,Cols,Options,MaxRows,MaxCols> ComplexMatrix; typedef Matrix<std::complex<RealScalar>,Rows,Cols,Options,MaxRows,MaxCols> ComplexMatrix;
const MatrixType& m_A; const MatrixType* m_A;
MatrixType m_tmp1, m_tmp2; MatrixType m_tmp1, m_tmp2;
ComplexMatrix m_T, m_U, m_fT; ComplexMatrix m_T, m_U, m_fT;
bool m_init; char m_flag;
RealScalar modfAndInit(RealScalar, RealScalar*); RealScalar modfAndInit(RealScalar, RealScalar*);
@ -75,6 +75,17 @@ template<typename MatrixType> class MatrixPower
*/ */
explicit MatrixPower(const MatrixType& A); explicit MatrixPower(const MatrixType& A);
/**
* \brief Constructor.
*
* \param[in] A the base of the matrix power.
*/
template<typename Derived>
explicit MatrixPower(const MatrixBase<Derived>& A);
/** \brief Destructor. */
~MatrixPower();
/** /**
* \brief Return the expression \f$ A^p \f$. * \brief Return the expression \f$ A^p \f$.
* *
@ -104,29 +115,39 @@ template<typename MatrixType> class MatrixPower
template<typename Derived, typename ResultType> template<typename Derived, typename ResultType>
void compute(const Derived& b, ResultType& res, RealScalar p); void compute(const Derived& b, ResultType& res, RealScalar 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(); }
}; };
template<typename MatrixType> template<typename MatrixType>
MatrixPower<MatrixType>::MatrixPower(const MatrixType& A) : MatrixPower<MatrixType>::MatrixPower(const MatrixType& A) :
m_A(A), m_A(&A),
m_init(false) m_flag(0)
{ /* empty body */ } { /* empty body */ }
template<typename MatrixType>
template<typename Derived>
MatrixPower<MatrixType>::MatrixPower(const MatrixBase<Derived>& A) :
m_A(new MatrixType(A)),
m_flag(2)
{ /* empty body */ }
template<typename MatrixType>
MatrixPower<MatrixType>::~MatrixPower()
{ if (m_flag & 2) delete m_A; }
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 (m_A->cols()) {
case 0: case 0:
break; break;
case 1: case 1:
res(0,0) = std::pow(m_A(0,0), p); res(0,0) = std::pow(m_A->coeff(0,0), p);
break; break;
default: default:
RealScalar intpart; RealScalar intpart, x = modfAndInit(p, &intpart);
RealScalar x = modfAndInit(p, &intpart); res = MatrixType::Identity(m_A->rows(), m_A->cols());
res = MatrixType::Identity(m_A.rows(),m_A.cols());
computeIntPower(res, intpart); computeIntPower(res, intpart);
computeFracPower(res, x); computeFracPower(res, x);
} }
@ -136,15 +157,14 @@ template<typename MatrixType>
template<typename Derived, typename ResultType> template<typename Derived, typename ResultType>
void MatrixPower<MatrixType>::compute(const Derived& b, ResultType& res, RealScalar p) void MatrixPower<MatrixType>::compute(const Derived& b, ResultType& res, RealScalar p)
{ {
switch (m_A.cols()) { switch (m_A->cols()) {
case 0: case 0:
break; break;
case 1: case 1:
res = std::pow(m_A(0,0), p) * b; res = std::pow(m_A->coeff(0,0), p) * b;
break; break;
default: default:
RealScalar intpart; RealScalar intpart, x = modfAndInit(p, &intpart);
RealScalar x = modfAndInit(p, &intpart);
computeIntPower(b, res, intpart); computeIntPower(b, res, intpart);
computeFracPower(res, x); computeFracPower(res, x);
} }
@ -157,11 +177,11 @@ typename MatrixType::RealScalar MatrixPower<MatrixType>::modfAndInit(RealScalar
*intpart = std::floor(x); *intpart = std::floor(x);
RealScalar res = x - *intpart; RealScalar res = x - *intpart;
if (!m_init && res) { if (!(m_flag & 1) && res) {
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();
m_init = true; m_flag |= 1;
const Array<RealScalar,EIGEN_SIZE_MIN_PREFER_FIXED(Rows,Cols),1,ColMajor,EIGEN_SIZE_MIN_PREFER_FIXED(MaxRows,MaxCols)> const Array<RealScalar,EIGEN_SIZE_MIN_PREFER_FIXED(Rows,Cols),1,ColMajor,EIGEN_SIZE_MIN_PREFER_FIXED(MaxRows,MaxCols)>
absTdiag = m_T.diagonal().array().abs(); absTdiag = m_T.diagonal().array().abs();
@ -194,8 +214,8 @@ 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_tmp1 = m_A->inverse();
else m_tmp1 = m_A; else m_tmp1 = *m_A;
while (pp >= 1) { while (pp >= 1) {
if (std::fmod(pp, 2) >= 1) if (std::fmod(pp, 2) >= 1)
@ -209,8 +229,8 @@ template<typename MatrixType>
template<typename Derived, typename ResultType> template<typename Derived, typename ResultType>
void MatrixPower<MatrixType>::computeIntPower(const Derived& b, ResultType& res, RealScalar p) void MatrixPower<MatrixType>::computeIntPower(const Derived& b, ResultType& res, RealScalar p)
{ {
if (b.cols() >= m_A.cols()) { if (b.cols() >= m_A->cols()) {
m_tmp2 = MatrixType::Identity(m_A.rows(),m_A.cols()); m_tmp2 = MatrixType::Identity(m_A->rows(), m_A->cols());
computeIntPower(m_tmp2, p); computeIntPower(m_tmp2, p);
res.noalias() = m_tmp2 * b; res.noalias() = m_tmp2 * b;
} }
@ -224,20 +244,20 @@ void MatrixPower<MatrixType>::computeIntPower(const Derived& b, ResultType& res,
return; return;
} }
else if (p>0) { else if (p>0) {
m_tmp1 = m_A; m_tmp1 = *m_A;
} }
else if (m_A.cols() > 2 && b.cols()*(pp-applyings) <= m_A.cols()*squarings) { else if (m_A->cols() > 2 && b.cols()*(pp-applyings) <= m_A->cols()*squarings) {
PartialPivLU<MatrixType> A(m_A); PartialPivLU<MatrixType> A(*m_A);
res = A.solve(b); res = A.solve(b);
for (--pp; pp >= 1; --pp) for (--pp; pp >= 1; --pp)
res = A.solve(res); res = A.solve(res);
return; return;
} }
else { else {
m_tmp1 = m_A.inverse(); m_tmp1 = m_A->inverse();
} }
while (b.cols()*(pp-applyings) > m_A.cols()*squarings) { while (b.cols()*(pp-applyings) > m_A->cols()*squarings) {
if (std::fmod(pp, 2) >= 1) { if (std::fmod(pp, 2) >= 1) {
apply(b, res, init); apply(b, res, init);
--applyings; --applyings;
@ -302,6 +322,7 @@ 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::RealScalar RealScalar; typedef typename Derived::RealScalar RealScalar;
typedef typename Derived::Index Index; typedef typename Derived::Index Index;
@ -322,7 +343,14 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
*/ */
template<typename ResultType> template<typename ResultType>
inline void evalTo(ResultType& res) const inline void evalTo(ResultType& res) const
{ MatrixPower<typename Derived::PlainObject>(m_A.eval()).compute(res, m_p); } { MatrixPower<PlainObject>(m_A).compute(res, m_p); }
template<typename OtherDerived>
const MatrixPowerMatrixProduct<PlainObject,OtherDerived> operator*(const MatrixBase<OtherDerived>& b) const
{
MatrixPower<PlainObject> Apow(m_A);
return MatrixPowerMatrixProduct<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(); }