aboutsummaryrefslogtreecommitdiffhomepage
path: root/unsupported/Eigen
diff options
context:
space:
mode:
Diffstat (limited to 'unsupported/Eigen')
-rw-r--r--unsupported/Eigen/src/MatrixFunctions/MatrixPower.h144
-rw-r--r--unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h107
2 files changed, 137 insertions, 114 deletions
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
index 15ede1c2a..996c24240 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -31,57 +31,19 @@ namespace Eigen {
* \include MatrixPower_optimal.cpp
* Output: \verbinclude MatrixPower_optimal.out
*/
-template<typename MatrixType> class MatrixPower
+template<typename MatrixType>
+class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType>
{
- private:
- static const int Rows = MatrixType::RowsAtCompileTime;
- static const int Cols = MatrixType::ColsAtCompileTime;
- static const int Options = MatrixType::Options;
- static const int MaxRows = MatrixType::MaxRowsAtCompileTime;
- static const int MaxCols = MatrixType::MaxColsAtCompileTime;
-
- typedef typename MatrixType::Scalar Scalar;
- typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
- typedef Matrix<std::complex<RealScalar>,Rows,Cols,Options,MaxRows,MaxCols> ComplexMatrix;
-
- const MatrixType* m_A;
- MatrixType m_tmp1, m_tmp2;
- ComplexMatrix m_T, m_U, m_fT;
- char m_flag;
-
- RealScalar modfAndInit(RealScalar, RealScalar*);
-
- template<typename Derived, typename ResultType>
- void apply(const Derived&, ResultType&, bool&);
-
- template<typename ResultType>
- void computeIntPower(ResultType&, RealScalar);
-
- template<typename Derived, typename ResultType>
- void computeIntPower(const Derived&, ResultType&, RealScalar);
-
- template<typename ResultType>
- void computeFracPower(ResultType&, RealScalar);
-
public:
- /**
- * \brief Constructor.
- *
- * \param[in] A the base of the matrix power.
- */
- explicit MatrixPower(const MatrixType& A);
+ EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(MatrixPower)
/**
* \brief Constructor.
*
* \param[in] A the base of the matrix power.
*/
- template<typename Derived>
- explicit MatrixPower(const MatrixBase<Derived>& A);
-
- /** \brief Destructor. */
- ~MatrixPower();
+ template<typename MatrixExpression>
+ explicit MatrixPower(const MatrixExpression& A);
/**
* \brief Return the expression \f$ A^p \f$.
@@ -111,40 +73,47 @@ template<typename MatrixType> class MatrixPower
*/
template<typename Derived, typename ResultType>
void compute(const Derived& b, ResultType& res, RealScalar p);
-
- Index rows() const { return m_A->rows(); }
- Index cols() const { return m_A->cols(); }
-};
-template<typename MatrixType>
-MatrixPower<MatrixType>::MatrixPower(const MatrixType& A) :
- m_A(&A),
- m_flag(0)
-{ /* empty body */ }
+ private:
+ using Base::m_A;
+ MatrixType m_tmp1, m_tmp2;
+ ComplexMatrix m_T, m_U, m_fT;
+ bool m_init;
-template<typename MatrixType>
-template<typename Derived>
-MatrixPower<MatrixType>::MatrixPower(const MatrixBase<Derived>& A) :
- m_A(new MatrixType(A)),
- m_flag(2)
-{ /* empty body */ }
+ RealScalar modfAndInit(RealScalar, RealScalar*);
+
+ template<typename Derived, typename ResultType>
+ void apply(const Derived&, ResultType&, bool&);
+
+ template<typename ResultType>
+ void computeIntPower(ResultType&, RealScalar);
+
+ template<typename Derived, typename ResultType>
+ void computeIntPower(const Derived&, ResultType&, RealScalar);
+
+ template<typename ResultType>
+ void computeFracPower(ResultType&, RealScalar);
+};
template<typename MatrixType>
-MatrixPower<MatrixType>::~MatrixPower()
-{ if (m_flag & 2) delete m_A; }
+template<typename MatrixExpression>
+MatrixPower<MatrixType>::MatrixPower(const MatrixExpression& A) :
+ Base(A),
+ m_init(false)
+{ /* empty body */ }
template<typename MatrixType>
void MatrixPower<MatrixType>::compute(MatrixType& res, RealScalar p)
{
- switch (m_A->cols()) {
+ switch (m_A.cols()) {
case 0:
break;
case 1:
- res(0,0) = std::pow(m_A->coeff(0,0), p);
+ res(0,0) = std::pow(m_A.coeff(0,0), p);
break;
default:
RealScalar intpart, 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);
computeFracPower(res, x);
}
@@ -154,11 +123,11 @@ template<typename MatrixType>
template<typename Derived, typename ResultType>
void MatrixPower<MatrixType>::compute(const Derived& b, ResultType& res, RealScalar p)
{
- switch (m_A->cols()) {
+ switch (m_A.cols()) {
case 0:
break;
case 1:
- res = std::pow(m_A->coeff(0,0), p) * b;
+ res = std::pow(m_A.coeff(0,0), p) * b;
break;
default:
RealScalar intpart, x = modfAndInit(p, &intpart);
@@ -168,20 +137,19 @@ void MatrixPower<MatrixType>::compute(const Derived& b, ResultType& res, RealSca
}
template<typename MatrixType>
-typename MatrixType::RealScalar MatrixPower<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
+typename MatrixPower<MatrixType>::Base::RealScalar MatrixPower<MatrixType>::modfAndInit(RealScalar x, RealScalar* intpart)
{
static RealScalar maxAbsEival, minAbsEival;
*intpart = std::floor(x);
RealScalar res = x - *intpart;
- if (!(m_flag & 1) && res) {
- const ComplexSchur<MatrixType> schurOfA(*m_A);
+ if (!m_init && res) {
+ const ComplexSchur<MatrixType> schurOfA(m_A);
m_T = schurOfA.matrixT();
m_U = schurOfA.matrixU();
- m_flag |= 1;
+ m_init = true;
- 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();
+ const RealArray absTdiag = m_T.diagonal().array().abs();
maxAbsEival = absTdiag.maxCoeff();
minAbsEival = absTdiag.minCoeff();
}
@@ -211,8 +179,8 @@ void MatrixPower<MatrixType>::computeIntPower(ResultType& res, RealScalar p)
{
RealScalar pp = std::abs(p);
- if (p<0) m_tmp1 = m_A->inverse();
- else m_tmp1 = *m_A;
+ if (p<0) m_tmp1 = m_A.inverse();
+ else m_tmp1 = m_A;
while (pp >= 1) {
if (std::fmod(pp, 2) >= 1)
@@ -226,8 +194,8 @@ 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());
+ 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;
}
@@ -241,20 +209,20 @@ void MatrixPower<MatrixType>::computeIntPower(const Derived& b, ResultType& res,
return;
}
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) {
- PartialPivLU<MatrixType> A(*m_A);
+ else if (m_A.cols() > 2 && 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();
+ 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) {
apply(b, res, init);
--applyings;
@@ -330,13 +298,13 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
* \param[in] p scalar, the exponent of the matrix power.
*/
MatrixPowerReturnValue(const Derived& A, RealScalar p)
- : m_pow(new MatrixPower<PlainObject>(A)), m_p(p), m_del(true) { }
+ : m_pow(*new MatrixPower<PlainObject>(A)), m_p(p), m_del(true) { }
MatrixPowerReturnValue(MatrixPower<PlainObject>& pow, RealScalar p)
- : m_pow(&pow), m_p(p), m_del(false) { }
+ : m_pow(pow), m_p(p), m_del(false) { }
~MatrixPowerReturnValue()
- { if (m_del) delete m_pow; }
+ { if (m_del) delete &m_pow; }
/**
* \brief Compute the matrix power.
@@ -346,17 +314,17 @@ class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Deriv
*/
template<typename ResultType>
inline void evalTo(ResultType& res) const
- { m_pow->compute(res, m_p); }
+ { m_pow.compute(res, m_p); }
template<typename OtherDerived>
const MatrixPowerMatrixProduct<PlainObject,OtherDerived> operator*(const MatrixBase<OtherDerived>& b) const
- { return MatrixPowerMatrixProduct<PlainObject,OtherDerived>(*m_pow, b.derived(), m_p); }
+ { return MatrixPowerMatrixProduct<PlainObject,OtherDerived>(m_pow, b.derived(), m_p); }
- Index rows() const { return m_pow->rows(); }
- Index cols() const { return m_pow->cols(); }
+ Index rows() const { return m_pow.rows(); }
+ Index cols() const { return m_pow.cols(); }
private:
- MatrixPower<PlainObject>* m_pow;
+ MatrixPower<PlainObject>& m_pow;
const RealScalar m_p;
const bool m_del; // whether to delete the pointer at destruction
MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&);
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
index a809609d5..ca5a604fc 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
@@ -29,30 +29,6 @@ struct recompose_complex_schur<0>
{ res = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); }
};
-template<typename Derived, typename _Lhs, typename _Rhs>
-struct traits<MatrixPowerProductBase<Derived,_Lhs,_Rhs> >
-{
- typedef MatrixXpr XprKind;
- typedef typename remove_all<_Lhs>::type Lhs;
- typedef typename remove_all<_Rhs>::type Rhs;
- typedef typename remove_all<Derived>::type PlainObject;
- typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar;
- typedef typename promote_storage_type<typename traits<Lhs>::StorageKind,
- typename traits<Rhs>::StorageKind>::ret StorageKind;
- typedef typename promote_index_type<typename traits<Lhs>::Index,
- typename traits<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<typename T>
inline int binary_powering_cost(T p, int* squarings)
{
@@ -121,7 +97,8 @@ inline int matrix_power_get_pade_degree(long double normIminusT)
}
} // namespace internal
-template<typename MatrixType, int UpLo=Upper> class MatrixPowerTriangularAtomic
+template<typename MatrixType, int UpLo=Upper>
+class MatrixPowerTriangularAtomic
{
private:
typedef typename MatrixType::Scalar Scalar;
@@ -239,10 +216,88 @@ void MatrixPowerTriangularAtomic<MatrixType,UpLo>::computeBig(MatrixType& res, R
compute2x2(res, p);
}
+#define EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(Derived) \
+ typedef MatrixPowerBase<Derived<MatrixType>,MatrixType> Base; \
+ using typename Base::Scalar; \
+ using typename Base::RealScalar; \
+ using typename Base::ComplexMatrix; \
+ using typename Base::RealArray;
+
#define EIGEN_MATRIX_POWER_PRODUCT_PUBLIC_INTERFACE(Derived) \
- typedef MatrixPowerProductBase<Derived, Lhs, Rhs > Base; \
+ typedef MatrixPowerProductBase<Derived, Lhs, Rhs> Base; \
EIGEN_DENSE_PUBLIC_INTERFACE(Derived)
+namespace internal {
+template<typename Derived, typename _Lhs, typename _Rhs>
+struct traits<MatrixPowerProductBase<Derived,_Lhs,_Rhs> >
+{
+ typedef MatrixXpr XprKind;
+ typedef typename remove_all<_Lhs>::type Lhs;
+ typedef typename remove_all<_Rhs>::type Rhs;
+ typedef typename remove_all<Derived>::type PlainObject;
+ typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar;
+ typedef typename promote_storage_type<typename traits<Lhs>::StorageKind,
+ typename traits<Rhs>::StorageKind>::ret StorageKind;
+ typedef typename promote_index_type<typename traits<Lhs>::Index,
+ typename traits<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
+ };
+};
+} // namespace internal
+
+template<typename Derived, typename MatrixType>
+class MatrixPowerBase
+{
+ protected:
+ static const int Rows = MatrixType::RowsAtCompileTime;
+ static const int Cols = MatrixType::ColsAtCompileTime;
+ static const int Options = MatrixType::Options;
+ static const int MaxRows = MatrixType::MaxRowsAtCompileTime;
+ static const int MaxCols = MatrixType::MaxColsAtCompileTime;
+
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::RealScalar RealScalar;
+ typedef typename MatrixType::Index Index;
+ typedef Matrix<std::complex<RealScalar>,Rows,Cols,Options,MaxRows,MaxCols> ComplexMatrix;
+ typedef Array<RealScalar,Rows,1,ColMajor,MaxRows> RealArray;
+
+ const MatrixType& m_A;
+ const bool m_del; // whether to delete the pointer at destruction
+
+ public:
+ explicit MatrixPowerBase(const MatrixType& A) :
+ m_A(A),
+ m_del(false)
+ { /* empty body */ }
+
+ template<typename OtherDerived>
+ explicit MatrixPowerBase(const MatrixBase<OtherDerived>& A) :
+ m_A(*new MatrixType(A)),
+ m_del(true)
+ { /* empty body */ }
+
+ ~MatrixPowerBase()
+ { if (m_del) delete &m_A; }
+
+ void compute(MatrixType& res, RealScalar p)
+ { static_cast<Derived*>(this)->compute(res,p); }
+
+ template<typename OtherDerived, typename ResultType>
+ void compute(const OtherDerived& b, ResultType& res, RealScalar p)
+ { static_cast<Derived*>(this)->compute(b,res,p); }
+
+ Index rows() const { return m_A.rows(); }
+ Index cols() const { return m_A.cols(); }
+};
+
template<typename Derived, typename Lhs, typename Rhs>
class MatrixPowerProductBase : public MatrixBase<Derived>
{
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-16 14:17:36 +0000